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 θ if P∼Beta(θp,θ(1−p)) and X|P∼Binomial(n,P). If X∼BetaBinomial(n,p,θ), then E[X]=npandVar[X]=np(1−p)θ+nθ+1.
Multinomial distribution
double dmultinom(int m, const double *prob, double *x, int give_log);
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.
Input:
m
is a positive integer, is the dimension of the random variable.prob
is a pointer to anm
-vector of probabilities.x
is a pointer to anm
-vector containing the data.give_log
is an integer:give_log=1
if the log probability is desired.give_log=0
if probability is desired.
The return value is the probability or log probability (as requested).
Euler-multinomial distribution
Simulate Euler-multinomial transitions
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.
On return, trans[0]
, …, trans[m-1]
will be the numbers of individuals making each of the respective transitions.
See ?reulermultinom
and FAQ 3.6 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.
Compute probabilities of Euler-multinomial transitions
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=1
if log probability is desired;give_log=0
if probability is desired.
The value returned is the probability or log probability (as requested).
See ?deulermultinom
and FAQ 3.6 for more on the Euler-multinomial distributions.
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=logp1−p. Its inverse is therefore p=11+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∈Rn+ to a vector Y∈Rn, where Yi=logXi∑jXj. For every c>0, this transformation maps the simplex {X∈Rn+|∑iXi=c} bijectively onto Rn.
The pseudo-inverse transformation takes Rn to the unit simplex S={X∈Rn+|∑iXi=1}. Specifically, Xi=eYi∑jeYj.
Note that if T:Rn+→Rn is the log-barycentric transformation so defined, U is the pseudo-inverse, and Id denotes the identity map, then T∘U=Id:Rn→Rn but U∘T≠Id. However, if T is restricted to the unit simplex S, then U∘T|S=Id:S→S.
Input:
x
is a pointer to vector of parameters to be tranformed either to or from log barycentric coordinates.n
is the length of this vector.xt
is a pointer to the vector that will hold the results.
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⋅y=∑ixiyi.
Exponential/geometric rate conversion
double exp2geom_rate_correction(double R, double dt);
This function computes r such that if N∼Geometric(p=1−e−rΔt) and T∼Exponential(rate=R), then E[NΔt]=E[T]. That is, r is the rate for an Euler process that gives the same expected waiting time as the exponential process it approximates. In particular r→R as Δt→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),
0]];
x[stateindex[3]];
x[stateindex[2]];
p[parindex[1]]; covars[covindex[
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:
p
is a pointer to parameter vector.t0
is the zero time.stateindex
,parindex
,covindex
: see Indices, above.covars
is a pointer to a vector containing the (possibly interpolated) values of the covariates at the current time.x
is a vector that will, on return, contain a draw from the initial-state 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.
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:
loglik
is a pointer to the scalar that will, on return, contain the log probability density.x
is the state vector at timet
.p
is a pointer to parameter vector.t0
is the zero time.stateindex
,parindex
,covindex
: see Indices, above.covars
is a pointer to a vector containing the (possibly interpolated) values of the covariates at the current time.
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 timet
+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 timet
.
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:
t1
,t2
are the times at the beginning and end of the interval, respectively.x1
,x2
are the state vectors at timet1
andt2
, respectively.p
is the parameter vector.stateindex
,parindex
,covindex
: see Indices, above.covars
is the vector of (possibly interpolated) covariates at timet
.loglik
is a pointer to the scalar that will, on return, contain the log probability density.
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:
t
is the time.x
is the state vector at timet
.p
is the parameter parameter vector.stateindex
,parindex
,covindex
: see Indices, above.covars
is the vector of (possibly interpolated) covariates at timet
.f
is a vector, of the same length asx
, that will, on return, contain the value of the map or vectorfield.
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:
t
is the time at the beginning of the Euler step.x
is the state vector at timet
.p
is the parameter vector.obsindex
,stateindex
,parindex
,covindex
: see Indices, above.covars
is the vector of (possibly interpolated) covariates at timet
.y
is a vector that will, on return, contain the simulated observations.
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:
y
is the vector of observables at timet
.x
is the state vector at timet
.p
is the parameter vector.give_log
is an integer:give_log=1
if log probability is desired;give_log=0
if probability is desired.
obsindex
,stateindex
,parindex
,covindex
: see Indices, above.covars
is the vector of (possibly interpolated) covariates at timet
.lik
is a pointer to a scalar that will, on return, contain the requested likelihood or log likelihood.
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:
t
is the time.x
is the state vector at timet
.p
is the parameter parameter vector.stateindex
,parindex
,covindex
: see Indices, above.covars
is the vector of (possibly interpolated) covariates at timet
.e
is a vector, of the same length asx
, that will, on return, contain the expected value of the observed variables.
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:
t
is the time.x
is the state vector at timet
.p
is the parameter parameter vector.stateindex
,parindex
,covindex
,vmatindex
: see Indices, above.covars
is the vector of (possibly interpolated) covariates at timet
.v
points to a square matrix that will, on return, contain the covariance matrix of the observed variables. In particular, if X is the latent state and thenobs
observables are Yi, thenv[vmatindex[i+nobs*j]]
contains Cov[Yi,Yj|X].- It is the user’s responsibility to ensure that the returned covariance matrix is symmetric: this is not checked.
rprior
void pomp_rprior (double *p, const int *parindex);
Description:
p
is the parameter vector.parindex
: see Indices, above.
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:
p
is the parameter vector.give_log
is an integer:give_log=1
if log probability is desired;give_log=0
if probability is desired.
parindex
: see Indices, above.lik
is a pointer to a scalar that will, on return, contain the requested probability density or log probability density.
partrans
void pomp_transform (double *pt, const double *p, const int *parindex);
Description:
p
is the parameter vector.parindex
: see Indices, above.pt
is the vector wherein the results will be returned.