phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sirs_pomp.c File Reference
#include "pomplink.h"
#include "internal.h"
Include dependency graph for sirs_pomp.c:

Go to the source code of this file.

Macros

#define Beta   (__p[__parindex[0]])
 
#define gamma   (__p[__parindex[1]])
 
#define psi   (__p[__parindex[2]])
 
#define omega   (__p[__parindex[3]])
 
#define S0   (__p[__parindex[4]])
 
#define I0   (__p[__parindex[5]])
 
#define R0   (__p[__parindex[6]])
 
#define N   (__p[__parindex[7]])
 
#define S   (__x[__stateindex[0]])
 
#define I   (__x[__stateindex[1]])
 
#define R   (__x[__stateindex[2]])
 
#define ll   (__x[__stateindex[3]])
 
#define ellI   (__x[__stateindex[4]])
 
#define node   (__x[__stateindex[5]])
 
#define EVENT_RATES
 
#define lik   (__lik[0])
 

Functions

static double event_rates (double *__x, const double *__p, double t, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double *rate, double *penalty)
 
void sirs_rinit (double *__x, const double *__p, double t, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars)
 Latent-state initializer (rinit). More...
 
void sirs_gill (double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)
 
void sirs_dmeas (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)
 Measurement model likelihood (dmeasure). More...
 

Variables

static const int nrate = 3
 

Macro Definition Documentation

◆ Beta

#define Beta   (__p[__parindex[0]])

Definition at line 6 of file sirs_pomp.c.

◆ ellI

#define ellI   (__x[__stateindex[4]])

Definition at line 18 of file sirs_pomp.c.

◆ EVENT_RATES

#define EVENT_RATES
Value:
event_rates(__x,__p,t, \
__stateindex,__parindex,__covindex, \
__covars,rate,&penalty) \
static double event_rates(double *__x, const double *__p, double t, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double *rate, double *penalty)
Definition: sirs_pomp.c:27

Definition at line 21 of file sirs_pomp.c.

◆ gamma

#define gamma   (__p[__parindex[1]])

Definition at line 7 of file sirs_pomp.c.

◆ I

#define I   (__x[__stateindex[1]])

Definition at line 15 of file sirs_pomp.c.

◆ I0

#define I0   (__p[__parindex[5]])

Definition at line 11 of file sirs_pomp.c.

◆ lik

#define lik   (__lik[0])

Definition at line 186 of file sirs_pomp.c.

◆ ll

#define ll   (__x[__stateindex[3]])

Definition at line 17 of file sirs_pomp.c.

◆ N

#define N   (__p[__parindex[7]])

Definition at line 13 of file sirs_pomp.c.

◆ node

#define node   (__x[__stateindex[5]])

Definition at line 19 of file sirs_pomp.c.

◆ omega

#define omega   (__p[__parindex[3]])

Definition at line 9 of file sirs_pomp.c.

◆ psi

#define psi   (__p[__parindex[2]])

Definition at line 8 of file sirs_pomp.c.

◆ R

#define R   (__x[__stateindex[2]])

Definition at line 16 of file sirs_pomp.c.

◆ R0

#define R0   (__p[__parindex[6]])

Definition at line 12 of file sirs_pomp.c.

◆ S

#define S   (__x[__stateindex[0]])

Definition at line 14 of file sirs_pomp.c.

◆ S0

#define S0   (__p[__parindex[4]])

Definition at line 10 of file sirs_pomp.c.

Function Documentation

◆ event_rates()

static double event_rates ( double *  __x,
const double *  __p,
double  t,
const int *  __stateindex,
const int *  __parindex,
const int *  __covindex,
const double *  __covars,
double *  rate,
double *  penalty 
)
static

Definition at line 26 of file sirs_pomp.c.

37  {
38  double event_rate = 0;
39  double alpha, disc;
40  *penalty = 0;
41  assert(I >= ellI);
42  assert(ellI >= 0);
43  assert(S >= 0);
44  // 0: transmission with saturation 0 or 1
45  alpha = Beta*S*I/N;
46  disc = (I > 0) ? ellI*(ellI-1)/I/(I+1) : 1;
47  event_rate += (*rate = alpha*(1-disc)); rate++;
48  *penalty += alpha*disc;
49  // 1: recovery
50  alpha = gamma*I;
51  if (I > ellI) {
52  event_rate += (*rate = alpha); rate++;
53  } else {
54  *rate = 0; rate++;
55  *penalty += alpha;
56  }
57  // 2: loss of immunity
58  alpha = omega*R;
59  event_rate += (*rate = alpha); rate++;
60  // sampling
61  alpha = psi*I;
62  *penalty += alpha;
63  assert(R_FINITE(event_rate));
64  return event_rate;
65 }
#define N
Definition: sirs_pomp.c:13
#define gamma
Definition: sirs_pomp.c:7
#define R
Definition: sirs_pomp.c:16
#define I
Definition: sirs_pomp.c:15
#define Beta
Definition: sirs_pomp.c:6
#define ellI
Definition: sirs_pomp.c:18
#define omega
Definition: sirs_pomp.c:9
#define psi
Definition: sirs_pomp.c:8
#define S
Definition: sirs_pomp.c:14

◆ sirs_dmeas()

void sirs_dmeas ( 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 
)

Measurement model likelihood (dmeasure).

Definition at line 189 of file sirs_pomp.c.

202  {
203  assert(!ISNAN(ll));
204  lik = (give_log) ? ll : exp(ll);
205 }
#define lik
Definition: sirs_pomp.c:186
#define ll
Definition: sirs_pomp.c:17

◆ sirs_gill()

void sirs_gill ( double *  __x,
const double *  __p,
const int *  __stateindex,
const int *  __parindex,
const int *  __covindex,
const double *  __covars,
double  t,
double  dt 
)

Latent-state process simulator (rprocess).

This integrates the filter equation.

Definition at line 90 of file sirs_pomp.c.

100  {
101  double tstep = 0.0, tmax = t + dt;
102  const int *nodetype = get_userdata_int("nodetype");
103  const int *sat = get_userdata_int("saturation");
104 
105  int parent = (int) nearbyint(node);
106 
107 #ifndef NDEBUG
108  int nnode = *get_userdata_int("nnode");
109  assert(parent>=0);
110  assert(parent<=nnode);
111 #endif
112 
113  // singular portion of filter equation
114  switch (nodetype[parent]) {
115  default: // non-genealogical event
116  break;
117  case 0: // root
118  ll = 0;
119  ellI += 1;
120  break;
121  case 1: // sample
122  ll = 0;
123  assert(I >= ellI);
124  assert(ellI >= 0);
125  if (sat[parent] == 1) {
126  ll += log(psi);
127  } else if (sat[parent] == 0) {
128  ellI -= 1;
129  ll += log(psi*(I-ellI));
130  } else {
131  assert(0); // #nocov
132  ll += R_NegInf; // #nocov
133  }
134  break;
135  case 2: // branch point s=(1,1)
136  ll = 0;
137  assert(S >= 0);
138  assert(I >= 0);
139  assert(ellI > 0);
140  assert(sat[parent]==2);
141  ll += (I > 0 && I >= ellI) ? log(Beta*S*I/N) : R_NegInf;
142  S -= 1; I += 1;
143  ellI += 1;
144  ll -= log(I*(I-1)/2);
145  S = (S > 0) ? S : 0;
146  break;
147  }
148 
149  if (tmax > t) {
150 
151  // take Gillespie steps to the end of the interval:
152  int event;
153  double penalty = 0;
154  double rate[nrate];
155 
156  double event_rate = EVENT_RATES;
157  tstep = exp_rand()/event_rate;
158 
159  while (t + tstep < tmax) {
160  event = rcateg(event_rate,rate,nrate);
161  ll -= penalty*tstep;
162  switch (event) {
163  case 0: // transmission
164  S -= 1; I += 1;
165  break;
166  case 1: // recovery
167  I -= 1; R += 1;
168  break;
169  case 2: // loss of immunity
170  R -= 1; S += 1;
171  break;
172  default: // #nocov
173  assert(0); // #nocov
174  break; // #nocov
175  }
176  t += tstep;
177  event_rate = EVENT_RATES;
178  tstep = exp_rand()/event_rate;
179  }
180  tstep = tmax - t;
181  ll -= penalty*tstep;
182  }
183  node += 1;
184 }
get_userdata_int_t * get_userdata_int
Definition: init.c:7
static int rcateg(double erate, double *rate, int nrate)
Definition: internal.h:76
static const int nrate
Definition: sirs_pomp.c:4
#define EVENT_RATES
Definition: sirs_pomp.c:21
#define node
Definition: sirs_pomp.c:19
Here is the call graph for this function:

◆ sirs_rinit()

void sirs_rinit ( double *  __x,
const double *  __p,
double  t,
const int *  __stateindex,
const int *  __parindex,
const int *  __covindex,
const double *  __covars 
)

Latent-state initializer (rinit).

Definition at line 68 of file sirs_pomp.c.

77  {
78  double m = N/(S0+I0+R0);
79  S = nearbyint(S0*m);
80  I = nearbyint(I0*m);
81  R = nearbyint(R0*m);
82  ll = 0;
83  ellI = 0;
84  node = 0;
85 }
#define R0
Definition: sirs_pomp.c:12
#define S0
Definition: sirs_pomp.c:10
#define I0
Definition: sirs_pomp.c:11

Variable Documentation

◆ nrate

const int nrate = 3
static

Definition at line 4 of file sirs_pomp.c.