rgsl

R gsl solver interface

---

Project maintained by a-kramer Hosted on GitHub Pages — Theme by mattgraham

rgsl

This R package solves a series of initial value problems given as an ordinary differential equation, an initial state and a numerical parameter vector. The C code calls gsl odeiv2 module functions to solve the problem or problems.

This project was funded by the EU as part of the Human Brain Project and supported by EBRAINS infrastructure.

Install

Using the remotes package:

remotes::install_github("a-kramer/rgsl")

Ensure that the GNU Scientific Library (gsl) is installed in your system and pkg-config can find it:

$ pkg-config --libs gsl

The above command should print something similar to:

-lgsl -lgslcblas -lm

Purpose

In this package we consider the initial value problem:

ẋ = f(x, t; p)x(t₀) = x₀

where p is a parameter vector.

We think that it is helpful to distinguish between constants (parameters that never change), unknown parameters, and input parameters.

Type	Symbol	Description
constants	NONE	never change, are built into the model
unknown	k	may be inferred from data
known	u	these distinguish different experiments

The overall vector p contains both the unknown k and the known u, e.g. as p=c(k,u) (aka a direct sum).

To determine the parameters k, the model f needs to be simulated in the provided experimental scenarios and the parameters are found through optimisation, or sampling (or perhaps another, similar method).

Typically, these methods need to solve the model many times. This is especially true for Bayesian methods (sampling). So, solving a batch of similar problems in one function call is desirable. In the next two Sections we discuss two common Scenarios.

Solver Calls

The high level solver interface functions assume that the user has a list of N simulation experiments and uses this call structure:

library("rgsl")
comment(NameOfModel) <- "path/to/sharedLibrary.so"
y <- r_gsl_odeiv2_outer(NameOfModel,experiments,p)

Where p refers to the parameters of the model. This may be a matrix of M parameter vectors, each column will result in a simulation of every experiment. The result is the outer product of M×M simulations.

Each list item of experiments is itself a list of named properties that define a simulation run, e.g.:

experiments[[1]][["time"]] <- seq(0,1,length.out=100)
experiments[[1]][["initial_value"]] <- c(0,0,0)
experiments[[1]][["input"]] <- c(1,2,3,4,5)

Return Values

a list y of the same size as the experiment list, each entry has a [["state"]] component and a [["func"]] (the function r_gsl_odeiv2_outer_state_only suppresses this) component (it is filled in if _func exists). Each y[[i]]$state is a 3d-array, where the second index corresponds to the output time points and the third enumerates the M different parameter vectors (possibly only 1).

The func component contains the output functions, as calculated by the ${MODEL}_func function. If that function is undefined, use

rgsl::r_gsl_odeiv2_outer_state_only # ${MODEL}_func is not used

Structure of Simulation Experiments

experiments:

• time simulation output time
• initial_value y(t=t₀)
• input partial parameter vector
• events OPTIONAL
  • time event occurence times (nt)
  • label transformation label, a 0-based offset (used as int) which selects the transformation to apply
  • dose a scalar floating point variable (used as double)

Models

This package requires the model to be provided as a shared library: ${MODEL}.so. The function names within the file need to correspond to the model name:

${MODEL}_vf();
${MODEL}_jac();

But, when calling the solvers in R, e.g.r_gsl_odeiv2_outer(modelName,...) the location of the shared library can be indicated through a comment on the modelName:

modelName <- "HarmonicOscillator"
comment(modelName) <- "../HaOs1.so"
y <- r_gsl_odeiv2_outer(modelName,...)

Otherwise, the path defaults to: paste0("./",modelName,".so")

The arguments of these functions conform to the official GSL documentation.

The model functions can be generated by icpm-kth/RPN-derivative/sh/ode.sh:

int $MODEL_vf(double t, const double y_[], double f_[], void *par)
int $MODEL_event(double t, double y_[], void *par, int EventLabel, double dose)
int $MODEL_jac(double t, const double y_[], double *jac_, double *dfdt_, void *par)
int $MODEL_jacp(double t, const double y_[], double *jacp_, double *dfdt_, void *par)
int $MODEL_func(double t, const double y_[], double *func_, void *par)
int $MODEL_funcJac(double t, const double y_[], double *funcJac_, void *par)
int $MODEL_funcJacp(double t, const double y_[], double *funcJacp_, void *par)
int $MODEL_default(double t, void *par)
int $MODEL_init(double t, double *y_, void *par)

where event can be called to perform a transformation, as described in the table of transformations, with EventLabel indicating which event to perform (a row offset, starting at 0).

We make some additional demands: when called with a NULL argument, the functions return the expected size of the output buffer:

[...]
size_t numPar=$MODEL_default(0,NULL);
double p =  malloc(sizeof(double)*numPar);
[...] // assign values to p
size_t numStateVar=$MODEL_vf(0,NULL,NULL,NULL);
double *y0 = malloc(sizeof(double)*numStateVar);
double *v = malloc(sizeof(double)*numStateVar);
[...] // assign values to y0
$MODEL_vf(0,y0,v,par);

The above code is OK. But, there is usually no need to call these functions explicitly. It is only important to know that this is allowed and the shared library model functions must not crash when used like this.

Enums

For the most part, the functions access the elements of vectors using enums:

enum stateVariable { _Cag, _Cal, ... , numStateVar }; /* state variable indexes  */

Similar enum definitions are made for parameters, events, and output functions:

enum param { ... , numParam }; /* parameter indexes  */
enum eventLabel { ... , numEvents }; /* event name indexes */
enum func { ... , numFunc }; /* parameter indexes  */

They are used in the event function (in a switch) and in some of the other functions – but, not in jacobians.

Example File

See the example file HarmonicOscillator to inspect an auto-generated right-hand-side and Jacobian function compatible with the solvers from the gsl.

The example script HarmonicOscillator.R contains a demo() function.

These C code files (for the gsl module odeiv2) can be automatically obtained using from icpm-kth/RPN-derivative, VFGEN or written manually.

SBtabVFGEN

This project integrates well with SBtabVFGEN as this package prints vfgen files from source files written using the SBtab format. SBtab is a format in systems biology and is suited for writing/drafting biological models. The workflow could be:

  User written:                           generated          Simulation
  +-----------+      +-----------+      +------------+      +----------+
  |           |      |           |      |  (CVODE)   |      |          |
  | SBtab (M) +--+-->+   VFGEN   +----->+  ODE code  +--+-->+   rgsl   |
  |           |  |   |           |      | +jacobian  |  |   |          |
  +-----------+  |   +-----------+      +------------+  |   +----------+
                 |                                      |
            +----+--------------+                 +----------------+
            |                   |                 |                |
            | sbtab_to_vfgen()  |                 | r_gsl_odeiv2() |
            |                   |                 |                |
            +-------------------+                 +----------------+