pomp C API

Overview

pomp provides C entry points to a number of facilities for model specification. Below, we describe these, providing the prototypical function calls and descriptions of their arguments.

The final section describes the prototypes for the basic model components. Users wishing to write libraries to hold basic model components must furnish functions of these prototypes that perform the basic model component computations.

Random variables

Beta-binomial distribution

double rbetabinom(double size, double prob, double theta);
double dbetabinom(double x, double size, double prob, 
  double theta, int give_log);

\(X\) is said to be Beta-binomially distributed with size \(n\), mean probability \(p\), and dispersion parameter \(\theta\) if \[P \sim \mathrm{Beta}\left(\theta\,p,\theta\,(1-p)\right)\] and \[X|P \sim \mathrm{Binomial}\left(n,P\right).\] If \(X\sim\mathrm{BetaBinomial}\left(n,p,\theta\right)\), then \[\mathbb{E}\left[X\right]=n\,p\qquad\text{and}\qquad\mathrm{Var}\left[X\right]=n\,p\,(1-p)\,\frac{\theta+n}{\theta+1}.\]

Multinomial distribution

The R C API provides a simulator for the multinomial distribution, rmultinom. See the “Rmath.h” header file for information on this facility. pomp provides an evaluator for the probability mass function (dmultinom) for this distribution.

double dmultinom(int m, const double *prob, double *x, int give_log);

Input:

The return value is the probability or log probability (as requested).

Euler-multinomial distribution

The Euler multinomial approximation of a continuous-time, stochastic compartmental model is as follows. Suppose a compartment has occupancy \(N_t\) at time \(t\) and that there are \(K\) ways of exiting the compartment, with per capita rates (hazards) \(\mu_1,\dots,\mu_K\), respectively.

**Diagram:** *A single compartment within a compartmental model. Here, there are $K=2$ paths out of the compartment.*

Diagram: A single compartment within a compartmental model. Here, there are \(K=2\) paths out of the compartment.

To make the Euler multinomial approximation, we approximate the total exit rate as constant over a small interval \([t,t+\Delta{t})\). Let the random variable \(\Delta{n_k}\), \(k=1,\dots,K\), be the number that exit by path \(k\) in this time interval and \(\Delta{n_0}\) be the number that remain. Under this assumption, the vector of numbers of exits, \((\Delta{n_{0}},\Delta{n_{1}},\dots,\Delta{n_{K}})\) is multinomially distributed with size \(N_t\) and probabilities \((p_k)_{k=0}^K\), where \[p_0 = \exp\left(-\sum\!\mu_i\,\Delta{t}\right),\] and \[p_k = \frac{\mu_k}{\sum\!\mu_i}\,\left(1-p_0\right),\qquad k=1,\dots,K.\] By way of shorthand, we say that \(\Delta{n}=(\Delta{n_k})_{k=1}^K\) is Euler-multinomially distributed with size \(N_t\), rates \(\mu=(\mu_k)_{k=1}^K\), and time-step \(\Delta{t}\) and we write \[\Delta{n} \sim \mathrm{Eulermultinom}\left(N_t,\mu,\Delta{t}\right).\]

The pomp C API provides three functions that relate to the Euler-multinomial distribution. Their descriptions follow.

Simulate an Euler-multinomial random variable

The reulermultinom function draws a random sample from this distribution. Using the notation above, one has to pack the \(K\) rates \(\mu_1,\dots,\mu_K\) into contiguous memory locations and retrieve the results in (a different set of) contiguous memory locations. For example, if rate is a pointer to \(K\) contiguous memory locations holding the rates and dn is a pointer to \(K\) contiguous memory locations ready to hold the results, then

reulermultinom(K,N,rate,dt,dn);

will result in a random sample from the Euler multinomial distribution (with timestep dt) being stored in dn[0], …, dn[K-1]. In the foregoing, we’ve assumed that the quantities \(N_t\) and \(K\) are stored in the integer variables N and K, respectively, and that the double precision variable dt holds the timestep.

The prototype is:

void reulermultinom(int m, double size, const double *rate,
  double dt, double *trans);

Input:

  • m, a positive integer, is number of potential transitions (“deaths”).
  • size, a positive integer, is the number of individuals at risk.
  • rate is a pointer to the vector of transition (“death”) rates.
  • dt, a positive real number, is the duration of time interval.
  • trans is a pointer to the vector that will hold the random deviate.

Output:

On return, trans[0], …, trans[m-1] will be the numbers of individuals making each of the respective transitions.

See ?reulermultinom for more on the Euler-multinomial distributions.

NB: reulermultinom does not call GetRNGstate() or PutRNGstate() internally. This must be done by the calling program. But note that when reulermultinom is called inside a pomp rprocess, there is no need to call either GetRNGState() or PutRNGState(); this is handled by pomp.

Probability distribution of an Euler-multinomial random variable

If \(\Delta{n} \sim \mathrm{Eulermultinom}\left(N_t,\mu,\Delta{t}\right)\), then the probability it takes a specific value can be computed using the C function deulermultinom. Its prototype is:

double deulermultinom(int m, double size, const double *rate,
  double dt, double *trans, int give_log);

Input:

  • m, a positive integer, is the number of potential transitions (“deaths”).
  • size, a positive integer, is the number of individuals at risk.
  • rate is a pointer to vector of transition (“death”) rates.
  • dt, a positive real number, is the duration of time interval.
  • trans is pointer to vector containing the data, which are numbers of individuals making the respective transitions.
  • give_log is an integer:
    • give_log=0 requests that the probability be returned.
    • give_log=1 requests that the log probability to be returned.

Output:

The value returned is the probability or log probability (as requested).

See also ?deulermultinom.

Expectation of an Euler-multinomial random variable

If \(\Delta{n} \sim \mathrm{Eulermultinom}\left(N_t,\mu,\Delta{t}\right)\), then the expectation of its \(i\)-th component is \[\mathbb{E}\left[\Delta{n}_i\right]=p_k N_t,\] where \(p_k\) is as defined above. The C function eeulermultinom computes this. Its prototype is:

void eeulermultinom(int m, double size, const double *rate,
  double dt, double *trans);

Input:

The parameters m, size, rate, and dt have the same meaning as above.

Output:

After a call to eeulermultinom, trans points to an array of doubles holding the expected values of the Euler-multinomial random variables.

Gamma white noise

double rgammawn(double sigma, double dt);

Corresponding to the R function rgammawn, this C function draws a single random increment of a gamma white-noise process. This will have expectation equal to dt and variance sigma^2*dt.

In particular, when dW = rgammawn(sigma,dt); is executed, mu*dW/dt is a candidate for a random rate process within an Euler-multinomial context, i.e., mu*dW will have expectation mu*dt and variance mu*sigma^2*dt.

Splines

void bspline_basis_eval(double x, double *knots, int degree,
  int nbasis, double *y);

void bspline_basis_eval_deriv(double x, double *knots, int degree,
  int nbasis, int deriv, double *y);

void periodic_bspline_basis_eval(double x, double period, int degree,
  int nbasis, double *y);

void periodic_bspline_basis_eval_deriv(double x, double period,
  int degree, int nbasis, int deriv, double *y);

These functions work with ordinary and periodic B-spline basis functions. For ordinary splines, one passes the spline knots in knots, together with the degree (degree) desired. Note that knots must point to an array of length at least nbasis + degree + 1. The first form evaluates the nbasis B-spline basis functions at x, the values being returned in the array pointed to by y. The second evaluates the order-deriv derivative of each basis function.

For periodic splines, the period and smoothness of the functions are given by period and degree, respectively. The fundamental domain is always assumed to be \([0,T]\), where \(T=\)period. The first form evaluates nbasis periodic B-spline basis functions at x, the values being returned in the array pointed to by y. The second evaluates the order-deriv derivative of each basis function.

Transformations

Logit transformation

double logit(double p);
double expit(double x);

The logit transformation is defined by \[x = \log\frac{p}{1-p}.\] Its inverse is therefore \[p = \frac{1}{1+e^{-x}}.\]

Log-barycentric transformation

void to_log_barycentric(double *xt, const double *x, int n);
void from_log_barycentric(double *xt, const double *x, int n);

The log-barycentric transformation takes a vector \(X\in\mathbb{R}^n_+\) to a vector \(Y\in\mathbb{R}^n\), where \[Y_i = \log\frac{X_i}{\sum_j\!X_j}.\] For every \(c>0\), this transformation maps the simplex \(\{X\in\mathbb{R}^n_+\;\vert\;\sum_i\!X_i = c\}\) bijectively onto \(\mathbb{R}^n\).

The pseudo-inverse transformation takes \(\mathbb{R}^n\) to the unit simplex \(S=\{X\in\mathbb{R}^n_+\;\vert\;\sum_i\!X_i=1\}\). Specifically, \[X_i = \frac{e^{Y_i}}{\sum_j\!e^{Y_j}}.\]

Note that if \(T:\mathbb{R}^n_+\to\mathbb{R}^n\) is the log-barycentric transformation so defined, \(U\) is the pseudo-inverse, and \(Id\) denotes the identity map, then \(T\circ U=Id:\mathbb{R}^n\to\mathbb{R}^n\) but \(U\circ T\ne Id\). However, if \(T\) is restricted to the unit simplex S, then \(U\circ{T\vert_{S}}=Id:S\to S\).

Input:

On return, xt[0], …, xt[n-1] will contain the transformed coordinates.

Convenience functions

Vector dot product

double dot_product(int n, const double *x, const double *y);

The return value is the dot (inner) product of the n-vectors x and y. By definition, \[x{\cdot}y = \sum_i\!x_i\,y_i.\]

Exponential/geometric rate conversion

double exp2geom_rate_correction(double R, double dt);

This function computes \(r\) such that if \[N \sim \mathrm{Geometric}\left(p=1-e^{-r\,\Delta{t}}\right)\] and \[T \sim \mathrm{Exponential}\left(\mathrm{rate}=R\right),\] then \(\mathbb{E}\left[N\,\Delta{t}\right] = \mathbb{E}\left[T\right]\). That is, \(r = \log{(1+R\,\Delta{t})}/{\Delta{t}}\) is the rate for an Euler process that gives the same expected waiting time as the exponential process it approximates. In particular \(r \to R\) as \(\Delta{t} \to 0\).

Access to the userdata

const int *get_userdata_int(const char *name);
const double *get_userdata_double(const char *name);
const SEXP get_userdata(const char *name);

The first function returns a pointer to the integer element of the userdata with the given name. The second retrieves a pointer to a double-precision element of userdata by name. The third form returns a general S expression (SEXP). See ?userdata for more information.

Prototypes for basic model components

pomp provides a facility whereby model codes can be compiled into a dynamically linked library for use in pomp objects. Specifically, basic model components are coded as C functions with the following prototypes.

NB: These functions should not be used within C snippets!

Indices

Each of the following functions is supplied one or more of the stateindex, parindex, covindex, obsindex, vmatindex arguments. Each of these is an integer vector: the integers within are indices giving the positions of specific model variables, according to the user’s specification, the latter being given by means of the statenames, paramnames, covarnames, and obsnames arguments. See ?pomp for more explanation. Thus, for example, within the body of a function of prototype pomp_rinit (see below),

 x[stateindex[0]];
 x[stateindex[3]];
 p[parindex[2]];
 covars[covindex[1]];

refer to the first state variable, the fourth state variable, the third parameter, and the second covariate, respectively.

rinit

void pomp_rinit (double *x, const double *p, double t0,
  const int *stateindex, const int *parindex, const int *covindex,
  const double *covars);

Description:

NB: There is no need to call GetRNGstate() or PutRNGstate() in the body of the user-defined function. The RNG is initialized before any call to this function, and the RNG state is written afterward. Inclusion of these calls in the user-defined function may result in significant slowdown.

dinit

void pomp_dinit (double *loglik, const double *x, const double *p, double t0,
  const int *stateindex, const int *parindex, const int *covindex,
  const double *covars);

Description:

rprocess

step.fun as used by euler and onestep

void pomp_onestep_sim (double *x, const double *p,
  const int *stateindex, const int *parindex, const int *covindex,
  const double *covars, double t, double dt);

Description:

  • p is the parameter vector.
  • stateindex, parindex, covindex: see Indices, above.
  • covars is the vector of covariates.
  • t is the time at the beginning of the step.
  • dt is the step size (duration of the interval).
  • x is the vector, that will, on return, contain a draw from the state process at time t+dt.

NB: There is no need to call GetRNGstate() or PutRNGstate() in the body of the user-defined function. The RNG is initialized before any call to this function, and the RNG state is written afterward. Inclusion of these calls in the user-defined function may result in significant slowdown.

rate.fun as used by gillespie

double pomp_ssa_rate_fn (int event, double t, const double *x, const double *p,
  const int *stateindex, const int *parindex, const int *covindex, 
  const double *covars);

Description:

  • event is an integer specifying the number of the reaction whose rate is desired (the first is event is 1, not 0).
  • t is the current time.
  • x is the vector of state variables.
  • p is the vector of parameters.
  • stateindex, parindex, covindex: see Indices, above.
  • covars is the vector of (possibly interpolated) covariates at time t.

The function returns the rate of the requested reaction.

dprocess

void pomp_dprocess (double *loglik,
  const double *x1, const double *x2, double t1, double t2, const double *p,
  const int *stateindex, const int *parindex, const int *covindex,
  const double *covars);

Description:

skeleton

void pomp_skeleton (double *f, const double *x, const double *p,
  const int *stateindex, const int *parindex, const int *covindex,
  const double *covars, double t);

Description:

rmeasure

void pomp_rmeasure (double *y, const double *x, const double *p,
  const int *obsindex, const int *stateindex, const int *parindex,
  const int *covindex, const double *covars, double t);

Description:

NB: There is no need to call GetRNGstate() or PutRNGstate() in the body of the user-defined function. The RNG is initialized before any call to this function, and the RNG state is written afterward. Inclusion of these calls in the user-defined function may result in significant slowdown.

dmeasure

void pomp_dmeasure (double *lik, const double *y, const double *x, const double *p,
  int give_log, const int *obsindex, const int *stateindex, const int *parindex,
  const int *covindex, const double *covars, double t);

Description:

emeasure

void pomp_emeasure (double *e, const double *x, const double *p,
  const int *obsindex, const int *stateindex, const int *parindex, const int *covindex,
  const double *covars, double t);

Description:

vmeasure

void pomp_vmeasure (double *v, const double *x, const double *p,
  const int *vmatindex, const int *stateindex, const int *parindex, const int *covindex,
  const double *covars, double t);

Description:

rprior

void pomp_rprior (double *p, const int *parindex);

Description:

On return, p will contain a new random draw from the prior distribution.

NB: There is no need to call GetRNGstate() or PutRNGstate() in the body of the user-defined function. The RNG is initialized before any call to this function, and the RNG state is written afterward. Inclusion of these calls in the user-defined function may result in significant slowdown.

dprior

void pomp_dprior (double *lik, const double *p, int give_log, const int *parindex);

Description:

partrans

void pomp_transform (double *pt, const double *p, const int *parindex);

Description: