phylopomp
Phylodynamics for POMPs
Loading...
Searching...
No Matches
Functions
decls.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

void R_init_phylopomp (DllInfo *)
 
void lbdp_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 Latent-state initializer (rinit).
 
void lbdp_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void lbdp_dmeas (double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double)
 Measurement model likelihood (dmeasure).
 
void seirs_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 
void seirs_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void seirs_dmeas (double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double)
 Measurement model likelihood (dmeasure).
 
void sirs_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 Latent-state initializer (rinit).
 
void sirs_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void sirs_dmeas (double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double)
 Measurement model likelihood (dmeasure).
 
void strains_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 Latent-state initializer (rinit).
 
void strains_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void strains_dmeas (double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double)
 Measurement model likelihood (dmeasure).
 
void twospecies_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 
void twospecies_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void twospecies_dmeas (double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double)
 Measurement model likelihood (dmeasure).
 

Function Documentation

◆ lbdp_dmeas()

void lbdp_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 )
extern

Measurement model likelihood (dmeasure).

Definition at line 167 of file lbdp_pomp.c.

180 {
181 assert(!ISNAN(ll));
182 lik = (give_log) ? ll : exp(ll);
183}
lik
#define lik
Definition lbdp_pomp.c:164
ll
#define ll
Definition lbdp_pomp.c:10

◆ lbdp_gill()

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

Latent-state process simulator (rprocess).

This integrates the filter equation.

Definition at line 72 of file lbdp_pomp.c.

82 {
83 double tstep = 0.0, tmax = t + dt;
84 const int *nodetype = get_userdata_int("nodetype");
85 const int *sat = get_userdata_int("saturation");
86 int parent = (int) nearbyint(node);
87
88#ifndef NDEBUG
89 int nnode = *get_userdata_int("nnode");
90 assert(parent>=0);
91 assert(parent<=nnode);
92#endif
93
94 ll = 0;
95
96 // singular portion of filter equation
97 switch (nodetype[parent]) {
98 default: // non-genealogical event #nocov
99 break; // #nocov
100 case 0: // root
101 ell += 1;
102 break;
103 case 1: // sample
104 assert(n >= ell);
105 assert(ell >= 0);
106 if (sat[parent] == 1) { // s=1
107 ll += log(psi);
108 } else if (sat[parent] == 0) { // s=0
109 ell -= 1;
110 double drate = chi*n;
111 double trate = drate+psi*(n-ell);
112 ll += (trate > 0) ? log(trate) : R_NegInf;
113 if (trate > 0 && unif_rand() < drate/trate) n -= 1;
114 } else {
115 assert(0); // #nocov
116 ll += R_NegInf; // #nocov
117 }
118 break;
119 case 2: // branch point s=2
120 ll = 0;
121 assert(n >= 0);
122 assert(ell > 0);
123 assert(sat[parent]==2);
124 n += 1; ell += 1;
125 ll += log(2*lambda/n);
126 break;
127 }
128
129 if (tmax > t) {
130
131 // Gillespie steps:
132 int event;
133 double penalty = 0;
134 double rate[2];
135
136 double event_rate = EVENT_RATES;
137 tstep = exp_rand()/event_rate;
138
139 while (t + tstep < tmax) {
140 event = rcateg(event_rate,rate,2);
141 assert(event>=0 && event<2);
142 ll -= penalty*tstep;
143 switch (event) {
144 case 0: // birth
145 n += 1;
146 break;
147 case 1: // death
148 n -= 1;
149 break;
150 default: // #nocov
151 assert(0); // #nocov
152 break; // #nocov
153 }
154 t += tstep;
155 event_rate = EVENT_RATES;
156 tstep = exp_rand()/event_rate;
157 }
158 tstep = tmax - t;
159 ll -= penalty*tstep;
160 }
161 node += 1;
162}
get_userdata_int
get_userdata_int_t * get_userdata_int
Definition init.c:7
rcateg
static int rcateg(double erate, double *rate, int nrate)
Definition internal.h:85
chi
#define chi
Definition lbdp_pomp.c:7
ell
#define ell
Definition lbdp_pomp.c:11
n
#define n
Definition lbdp_pomp.c:9
lambda
#define lambda
Definition lbdp_pomp.c:4
psi
#define psi
Definition lbdp_pomp.c:6
EVENT_RATES
#define EVENT_RATES
Definition lbdp_pomp.c:14
node
#define node
Definition lbdp_pomp.c:12
Here is the call graph for this function:

◆ lbdp_rinit()

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

Latent-state initializer (rinit).

Definition at line 53 of file lbdp_pomp.c.

62 {
63 n = nearbyint(n0);
64 ll = 0;
65 ell = 0;
66 node = 0;
67}
n0
#define n0
Definition lbdp_pomp.c:8

◆ R_init_phylopomp()

void R_init_phylopomp ( DllInfo * info)
extern

Definition at line 58 of file init.c.

58 {
59 // Register routines
60 R_registerRoutines(info,NULL,callMethods,NULL,extMethods);
61 R_useDynamicSymbols(info,TRUE);
62 // R_useDynamicSymbols(info,FALSE);
63 // R_forceSymbols(info,TRUE);
64 get_userdata = (get_userdata_t*) R_GetCCallable("pomp","get_userdata");
65 get_userdata_double = (get_userdata_double_t*) R_GetCCallable("pomp","get_userdata_double");
66 get_userdata_int = (get_userdata_int_t*) R_GetCCallable("pomp","get_userdata_int");
67}
extMethods
static const R_CallMethodDef extMethods[]
Definition init.c:52
callMethods
static const R_CallMethodDef callMethods[]
Definition init.c:32
get_userdata
get_userdata_t * get_userdata
Definition init.c:5
get_userdata_double
get_userdata_double_t * get_userdata_double
Definition init.c:6

◆ seirs_dmeas()

void seirs_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 )
extern

Measurement model likelihood (dmeasure).

Definition at line 351 of file seirs_pomp.c.

364 {
365 assert(!ISNAN(ll));
366 lik = (give_log) ? ll : exp(ll);
367}

◆ seirs_gill()

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

Simulator for the latent-state process (rprocess).

This is the Gillespie algorithm applied to the solution of the filter equation for the SEIRS process.

Definition at line 143 of file seirs_pomp.c.

153 {
154 double tstep = 0.0, tmax = t + dt;
155 double *color = &COLOR;
156 const int nsample = *get_userdata_int("nsample");
157 const int *nodetype = get_userdata_int("nodetype");
158 const int *nodedeme = get_userdata_int("deme");
159 const int *lineage = get_userdata_int("lineage");
160 const int *sat = get_userdata_int("saturation");
161 const int *index = get_userdata_int("index");
162 const int *child = get_userdata_int("child");
163
164 int parent = (int) nearbyint(node);
165
166#ifndef NDEBUG
167 int nnode = *get_userdata_int("nnode");
168 assert(parent>=0);
169 assert(parent<=nnode);
170#endif
171
172 int parlin = lineage[parent];
173 int parcol = color[parlin];
174 int deme = nodedeme[parent];
175 assert(parlin >= 0 && parlin < nsample);
176
177 ll = 0;
178
179 // singular portion of filter equation
180 switch (nodetype[parent]) {
181 default: // non-genealogical event #nocov
182 break; // #nocov
183 case 0: // root
184 // color lineages by sampling without replacement
185 assert(sat[parent]==1);
186 int c = child[index[parent]];
187 assert(lineage[parent]==lineage[c]);
188 if (E-ellE + I-ellI > 0) {
189 double x = (E-ellE)/(E-ellE + I-ellI);
190 if (unif_rand() < x) { // lineage is put into E deme
191 color[lineage[c]] = Exposed;
192 ellE += 1;
193 ll -= log(x);
194 } else { // lineage is put into I deme
195 color[lineage[c]] = Infected;
196 ellI += 1;
197 ll -= log(1-x);
198 }
199 } else { // more roots than infectives
200 ll += R_NegInf; // this is incompatible with the genealogy
201 // the following keeps the state valid
202 if (unif_rand() < 0.5) { // lineage is put into E deme
203 color[lineage[c]] = Exposed;
204 ellE += 1; E += 1;
205 // ll -= log(0.5);
206 } else { // lineage is put into I deme
207 color[lineage[c]] = Infected;
208 ellI += 1; I += 1;
209 // ll -= log(0.5);
210 }
211 }
212 break;
213 case 1: // sample
214 // If parent is not in deme I, likelihood = 0.
215 assert(deme==Infected);
216 if (parcol != deme) {
217 ll += R_NegInf;
218 color[parlin] = Infected;
219 // the following keeps the state valid
220 ellE -= 1; ellI += 1;
221 E -= 1; I += 1;
222 }
223 if (sat[parent] == 1) { // s=(0,1)
224 int c = child[index[parent]];
225 color[lineage[c]] = Infected;
226 ll += log(psi);
227 } else if (sat[parent] == 0) { // s=(0,0)
228 ellI -= 1;
229 ll += log(psi*(I-ellI));
230 } else {
231 assert(0); // #nocov
232 ll += R_NegInf; // #nocov
233 }
234 color[parlin] = R_NaReal;
235 break;
236 case 2: // branch point s=(1,1)
237 // If parent is not in deme I, likelihood = 0.
238 if (parcol != Infected) {
239 ll += R_NegInf;
240 color[parlin] = Infected;
241 // the following keeps the state valid
242 ellE -= 1; ellI += 1;
243 E -= 1; I += 1;
244 }
245 assert(sat[parent]==2);
246 ll += (S > 0 && I > 0) ? log(Beta*S/N/(E+1)) : R_NegInf;
247 S -= 1; E += 1;
248 ellE += 1;
249 S = (S > 0) ? S : 0;
250 int c1 = child[index[parent]];
251 int c2 = child[index[parent]+1];
252 assert(c1 != c2);
253 assert(lineage[c1] != lineage[c2]);
254 assert(lineage[c1] != parlin || lineage[c2] != parlin);
255 assert(lineage[c1] == parlin || lineage[c2] == parlin);
256 if (unif_rand() < 0.5) {
257 color[lineage[c1]] = Exposed;
258 color[lineage[c2]] = Infected;
259 } else {
260 color[lineage[c1]] = Infected;
261 color[lineage[c2]] = Exposed;
262 }
263 ll -= log(0.5);
264 break;
265 }
266
267 // continuous portion of filter equation:
268 // take Gillespie steps to the end of the interval
269 if (tmax > t) {
270
271 double rate[nrate], logpi[nrate];
272 int event;
273 double event_rate = 0;
274 double penalty = 0;
275
276 event_rate = EVENT_RATES;
277 tstep = exp_rand()/event_rate;
278
279 while (t + tstep < tmax) {
280 event = rcateg(event_rate,rate,nrate);
281 assert(event>=0 && event<nrate);
282 ll -= penalty*tstep + logpi[event];
283 switch (event) {
284 case 0: // transmission, s=(0,0)
285 assert(S>=1);
286 S -= 1; E += 1;
287 ll += log(1-ellI/I)+log(1-ellE/E);
288 assert(!ISNAN(ll));
289 break;
290 case 1: // transmission, s=(0,1)
291 assert(S>=1);
292 S -= 1; E += 1;
293 ll += log(1-ellE/E)-log(I);
294 assert(!ISNAN(ll));
295 break;
296 case 2: // transmission, s=(1,0)
297 assert(S>=1);
298 change_color(color,nsample,random_choice(ellI),Infected,Exposed);
299 ellE += 1; ellI -= 1;
300 S -= 1; E += 1;
301 ll += log(1-ellI/I)-log(E);
302 assert(!ISNAN(ll));
303 break;
304 case 3: // progression, s=(0,0)
305 assert(E>=1);
306 E -= 1; I += 1;
307 ll += log(1-ellI/I);
308 assert(!ISNAN(ll));
309 break;
310 case 4: // progression, s=(0,1)
311 assert(E>=1);
312 change_color(color,nsample,random_choice(ellE),Exposed,Infected);
313 ellE -= 1; ellI += 1;
314 E -= 1; I += 1;
315 ll -= log(I);
316 assert(!ISNAN(ll));
317 break;
318 case 5: // recovery
319 assert(I>=1);
320 I -= 1; R += 1;
321 assert(!ISNAN(ll));
322 break;
323 case 6: // waning
324 assert(R>=1);
325 R -= 1; S += 1;
326 assert(!ISNAN(ll));
327 break;
328 default: // #nocov
329 assert(0); // #nocov
330 ll += R_NegInf; // #nocov
331 break; // #nocov
332 }
333
334 ellE = nearbyint(ellE);
335 ellI = nearbyint(ellI);
336
337 t += tstep;
338 event_rate = EVENT_RATES;
339 tstep = exp_rand()/event_rate;
340
341 }
342 tstep = tmax - t;
343 ll -= penalty*tstep;
344 }
345 node += 1;
346}
nsample
SEXP nsample(TYPE &X)
Definition generics.h:12
deme
static const int deme
Definition lbdp.cc:7
N
#define N
Definition seirs_pomp.c:36
ellE
#define ellE
Definition seirs_pomp.c:43
E
#define E
Definition seirs_pomp.c:38
change_color
static void change_color(double *color, int nsample, int n, int from, int to)
Definition seirs_pomp.c:13
COLOR
#define COLOR
Definition seirs_pomp.c:45
random_choice
static int random_choice(double n)
Definition seirs_pomp.c:9
R
#define R
Definition seirs_pomp.c:40
I
#define I
Definition seirs_pomp.c:39
Beta
#define Beta
Definition seirs_pomp.c:27
ellI
#define ellI
Definition seirs_pomp.c:44
nrate
static const int nrate
Definition seirs_pomp.c:7
Infected
#define Infected
Definition seirs_pomp.c:5
Exposed
#define Exposed
Definition seirs_pomp.c:4
S
#define S
Definition seirs_pomp.c:37
Here is the call graph for this function:

◆ seirs_rinit()

void seirs_rinit ( double * __x,
const double * __p,
double t0,
const int * __stateindex,
const int * __parindex,
const int * __covindex,
const double * __covars )
extern

Latent-state initializer (rinit component).

The state variables include S, E, I, R plus 'ellE' and 'ellI' (numbers of E- and I-deme lineages), the accumulated weight ('ll'), the current node number ('node'), and the coloring of each lineage ('COLOR').

Definition at line 118 of file seirs_pomp.c.

127 {
128 double adj = N/(S0+E0+I0+R0);
129 S = nearbyint(S0*adj);
130 E = nearbyint(E0*adj);
131 I = nearbyint(I0*adj);
132 R = nearbyint(R0*adj);
133 ellE = 0;
134 ellI = 0;
135 ll = 0;
136 node = 0;
137}
R0
#define R0
Definition seirs_pomp.c:35
S0
#define S0
Definition seirs_pomp.c:32
I0
#define I0
Definition seirs_pomp.c:34
E0
#define E0
Definition seirs_pomp.c:33

◆ 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 )
extern

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}

◆ 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 )
extern

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 ll = 0;
114
115 // singular portion of filter equation
116 switch (nodetype[parent]) {
117 default: // non-genealogical event #nocov
118 break; // #nocov
119 case 0: // root
120 ellI += 1;
121 break;
122 case 1: // sample
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 assert(S >= 0);
137 assert(I >= 0);
138 assert(ellI > 0);
139 assert(sat[parent]==2);
140 ll += (I > 0 && I >= ellI) ? log(Beta*S*I/N) : R_NegInf;
141 S -= 1; I += 1;
142 ellI += 1;
143 ll -= log(I*(I-1)/2);
144 S = (S > 0) ? S : 0;
145 break;
146 }
147
148 if (tmax > t) {
149
150 // take Gillespie steps to the end of the interval:
151 int event;
152 double penalty = 0;
153 double rate[nrate];
154
155 double event_rate = EVENT_RATES;
156 tstep = exp_rand()/event_rate;
157
158 while (t + tstep < tmax) {
159 event = rcateg(event_rate,rate,nrate);
160 assert(event>=0 && event<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}
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 )
extern

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}

◆ strains_dmeas()

void strains_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 )
extern

Measurement model likelihood (dmeasure).

Definition at line 294 of file strains_pomp.c.

307 {
308 assert(!ISNAN(ll));
309 lik = (give_log) ? ll : exp(ll);
310}

◆ strains_gill()

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

Latent-state process simulator (rprocess).

This integrates the filter equation.

Definition at line 140 of file strains_pomp.c.

150 {
151 double tstep = 0.0, tmax = t + dt;
152 const int *nodetype = get_userdata_int("nodetype");
153 const int *sat = get_userdata_int("saturation");
154 const int *deme = get_userdata_int("deme");
155 const int *index = get_userdata_int("index");
156 const int *child = get_userdata_int("child");
157
158 int parent = (int) nearbyint(node);
159 int c = child[index[parent]];
160
161#ifndef NDEBUG
162 const int *lineage = get_userdata_int("lineage");
163 int nnode = *get_userdata_int("nnode");
164 assert(parent>=0);
165 assert(parent<=nnode);
166#endif
167
168 ll = 0;
169
170 // singular portion of filter equation
171 switch (nodetype[parent]) {
172 default: // non-genealogical event #nocov
173 break; // #nocov
174 case 0: // root
175 assert(sat[parent]==1);
176 assert(lineage[parent]==lineage[c]);
177 switch (deme[c]) {
178 case STRAIN1:
179 ellI1 += 1; break;
180 case STRAIN2:
181 ellI2 += 1; break;
182 case STRAIN3:
183 ellI3 += 1; break;
184 default: // #nocov
185 assert(0); break; // #nocov
186 }
187 break;
188 case 1: // sample
189 assert(sat[parent]==0);
190 switch (deme[parent]) {
191 case STRAIN1:
192 assert(I1 >= ellI1);
193 assert(ellI1 >= 0);
194 ellI1 -= 1; I1 -= 1;
195 ll += log(psi1);
196 break;
197 case STRAIN2:
198 assert(I2 >= ellI2);
199 assert(ellI2 >= 0);
200 ellI2 -= 1; I2 -= 1;
201 ll += log(psi2);
202 break;
203 case STRAIN3:
204 assert(I3 >= ellI3);
205 assert(ellI3 >= 0);
206 ellI3 -= 1; I3 -= 1;
207 ll += log(psi3);
208 break;
209 default: // #nocov
210 assert(0); break; // #nocov
211 }
212 break;
213 case 2: // branch point s=(1,1)
214 assert(S >= 0);
215 if (sat[parent]!=2) break;
216 switch (deme[parent]) {
217 case STRAIN1:
218 assert(I1 >= 0);
219 assert(ellI1 > 0);
220 ll += (I1 > 0 && I1 >= ellI1) ? log(Beta1*S*I1/N) : R_NegInf;
221 S -= 1; I1 += 1; ellI1 += 1;
222 ll -= log(I1*(I1-1)/2);
223 break;
224 case STRAIN2:
225 assert(I2 >= 0);
226 assert(ellI2 > 0);
227 ll += (I2 > 0 && I2 >= ellI2) ? log(Beta2*S*I2/N) : R_NegInf;
228 S -= 1; I2 += 1; ellI2 += 1;
229 ll -= log(I2*(I2-1)/2);
230 break;
231 case STRAIN3:
232 assert(I3 >= 0);
233 assert(ellI3 > 0);
234 ll += (I3 > 0 && I3 >= ellI3) ? log(Beta3*S*I3/N) : R_NegInf;
235 S -= 1; I3 += 1; ellI3 += 1;
236 ll -= log(I3*(I3-1)/2);
237 break;
238 default: // #nocov
239 assert(0); break; // #nocov
240 }
241 S = (S > 0) ? S : 0;
242 break;
243 }
244
245 if (tmax > t) {
246
247 // take Gillespie steps to the end of the interval:
248 int event;
249 double penalty = 0;
250 double rate[nrate];
251
252 double event_rate = EVENT_RATES;
253 tstep = exp_rand()/event_rate;
254
255 while (t + tstep < tmax) {
256 event = rcateg(event_rate,rate,nrate);
257 assert(event>=0 && event<nrate);
258 ll -= penalty*tstep;
259 switch (event) {
260 case 0: // strain-1 transmission
261 S -= 1; I1 += 1;
262 break;
263 case 1: // strain-1 recovery
264 I1 -= 1; R += 1;
265 break;
266 case 2: // strain-2 transmission
267 S -= 1; I2 += 1;
268 break;
269 case 3: // strain-2 recovery
270 I2 -= 1; R += 1;
271 break;
272 case 4: // strain-3 transmission
273 S -= 1; I3 += 1;
274 break;
275 case 5: // strain-3 recovery
276 I3 -= 1; R += 1;
277 break;
278 default: // #nocov
279 assert(0); break; // #nocov
280 }
281 t += tstep;
282 event_rate = EVENT_RATES;
283 tstep = exp_rand()/event_rate;
284 }
285 tstep = tmax - t;
286 ll -= penalty*tstep;
287 }
288 node += 1;
289}
psi1
#define psi1
Definition strains_pomp.c:14
Beta3
#define Beta3
Definition strains_pomp.c:12
psi3
#define psi3
Definition strains_pomp.c:16
STRAIN3
#define STRAIN3
Definition strains_pomp.c:8
ellI1
#define ellI1
Definition strains_pomp.c:29
ellI2
#define ellI2
Definition strains_pomp.c:30
ellI3
#define ellI3
Definition strains_pomp.c:31
STRAIN1
#define STRAIN1
Definition strains_pomp.c:6
STRAIN2
#define STRAIN2
Definition strains_pomp.c:7
I2
#define I2
Definition strains_pomp.c:25
Beta2
#define Beta2
Definition strains_pomp.c:11
Beta1
#define Beta1
Definition strains_pomp.c:10
psi2
#define psi2
Definition strains_pomp.c:15
I3
#define I3
Definition strains_pomp.c:26
I1
#define I1
Definition strains_pomp.c:24
Here is the call graph for this function:

◆ strains_rinit()

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

Latent-state initializer (rinit).

Definition at line 114 of file strains_pomp.c.

123 {
124 double m = N/(S_0+I1_0+I2_0+I3_0+R_0);
125 S = nearbyint(S_0*m);
126 I1 = nearbyint(I1_0*m);
127 I2 = nearbyint(I2_0*m);
128 I3 = nearbyint(I3_0*m);
129 R = nearbyint(R_0*m);
130 ll = 0;
131 ellI1 = 0;
132 ellI2 = 0;
133 ellI3 = 0;
134 node = 0;
135}
I2_0
#define I2_0
Definition strains_pomp.c:19
S_0
#define S_0
Definition strains_pomp.c:17
I3_0
#define I3_0
Definition strains_pomp.c:20
I1_0
#define I1_0
Definition strains_pomp.c:18
R_0
#define R_0
Definition strains_pomp.c:21

◆ twospecies_dmeas()

void twospecies_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 )
extern

Measurement model likelihood (dmeasure).

Definition at line 637 of file twospecies_pomp.c.

650 {
651 assert(!ISNAN(ll));
652 lik = (give_log) ? ll : exp(ll);
653}

◆ twospecies_gill()

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

Simulator for the latent-state process (rprocess).

This is the Gillespie algorithm applied to the solution of the filter equation for the TwoSpecies model. It advances the state from time t to time t+dt.

A tricky aspect of this function is that it must return a "valid" state even when the state is incompatible with the genealogy. In such a case, we set the log likelihood (ll) to R_NegInf, but the state must remain valid. Hence insertion of extra infectives, etc.

FIXME: At the moment, the following codes exclude the possibility of importation of infection.

Definition at line 290 of file twospecies_pomp.c.

300 {
301 double tstep = 0.0, tmax = t + dt;
302 double *color = &COLOR;
303 const int nsample = *get_userdata_int("nsample");
304 const int *nodetype = get_userdata_int("nodetype");
305 const int *nodedeme = get_userdata_int("deme");
306 const int *lineage = get_userdata_int("lineage");
307 const int *sat = get_userdata_int("saturation");
308 const int *index = get_userdata_int("index");
309 const int *child = get_userdata_int("child");
310
311 int parent = (int) nearbyint(node);
312
313#ifndef NDEBUG
314 int nnode = *get_userdata_int("nnode");
315 assert(parent>=0);
316 assert(parent<=nnode);
317#endif
318
319 int parlin = lineage[parent];
320 int parcol = color[parlin];
321 int deme = nodedeme[parent];
322 assert(parlin >= 0 && parlin < nsample);
323 assert(nearbyint(N1)==nearbyint(S1+I1+R1));
324 assert(nearbyint(N2)==nearbyint(S2+I2+R2));
325 assert(check_color(color,nsample,ell1,ell2));
326
327 ll = 0;
328
329 // singular portion of filter equation
330 switch (nodetype[parent]) {
331 default: // non-genealogical event #nocov
332 break; // #nocov
333 case 0: // root
334 // color lineages by sampling without replacement
335 assert(sat[parent]==1);
336 int c = child[index[parent]];
337 assert(parlin==lineage[c]);
338 if (I1-ell1+I2-ell2 > 0) {
339 double x = (I1-ell1)/(I1-ell1 + I2-ell2);
340 if (unif_rand() < x) { // lineage is put into I1 deme
341 color[lineage[c]] = host1;
342 ell1 += 1;
343 ll -= log(x);
344 } else { // lineage is put into I2 deme
345 color[lineage[c]] = host2;
346 ell2 += 1;
347 ll -= log(1-x);
348 }
349 } else { // more roots than infectives
350 ll += R_NegInf; // this is incompatible with the genealogy
351 // the following keeps the state valid
352 if (unif_rand() < 0.5) { // lineage is put into I1 deme
353 color[lineage[c]] = host1;
354 ell1 += 1; I1 += 1; N1 += 1;
355 } else { // lineage is put into I2 deme
356 color[lineage[c]] = host2;
357 ell2 += 1; I2 += 1; N2 += 1;
358 }
359 }
360 assert(nearbyint(N1)==nearbyint(S1+I1+R1));
361 assert(nearbyint(N2)==nearbyint(S2+I2+R2));
362 assert(check_color(color,nsample,ell1,ell2));
363 break;
364 case 1: // sample
365 if (parcol != deme) { // parent color does not match the observed deme
366 ll += R_NegInf;
367 }
368 if (sat[parent] == 0) { // s=(0,0)
369 if (parcol == host1) {
370 ell1 -= 1;
371 if (C1 < 1 && unif_rand() > C1) {
372 ll += log(psi1*(I1-ell1));
373 } else {
374 ll += log(psi1*I1);
375 I1 -= 1; N1 -= 1;
376 }
377 } else if (parcol == host2) {
378 ell2 -= 1;
379 if (C2 < 1 && unif_rand() > C2) {
380 ll += log(psi2*(I2-ell2));
381 } else {
382 ll += log(psi2*I2);
383 I2 -= 1; N2 -= 1;
384 }
385 } else {
386 assert(0); // #nocov
387 ll += R_NegInf; // #nocov
388 }
389 } else if (sat[parent] == 1) {
390 int c = child[index[parent]];
391 color[lineage[c]] = parcol;
392 if (parcol==host1) {
393 ll += log(psi1*(1-C1)); // s=(1,0)
394 } else if (parcol==host2) {
395 ll += log(psi2*(1-C2)); // s=(0,1)
396 } else {
397 assert(0); // #nocov
398 ll += R_NegInf; // #nocov
399 }
400 } else {
401 assert(0); // #nocov
402 ll += R_NegInf; // #nocov
403 }
404 color[parlin] = R_NaReal;
405 assert(nearbyint(N1)==nearbyint(S1+I1+R1));
406 assert(nearbyint(N2)==nearbyint(S2+I2+R2));
407 assert(check_color(color,nsample,ell1,ell2));
408 break;
409 case 2: // branch point
410 assert(sat[parent]==2);
411 if (parcol == host1) { // parent is in I1
412 assert(S1>=0 && S2 >=0 && I1>=ell1 && ell1>=0);
413 double lambda11 = Beta11*S1*I1/N1;
414 double lambda21 = Beta21*S2*I1/N1;
415 double lambda = lambda11+lambda21;
416 double x = (lambda > 0) ? lambda11/lambda : 0;
417 int c1 = child[index[parent]];
418 int c2 = child[index[parent]+1];
419 assert(c1 != c2);
420 assert(lineage[c1] != lineage[c2]);
421 assert(lineage[c1] != parlin || lineage[c2] != parlin);
422 assert(lineage[c1] == parlin || lineage[c2] == parlin);
423 if (unif_rand() < x) { // s = (2,0)
424 color[lineage[c1]] = host1;
425 color[lineage[c2]] = host1;
426 if (S1 > 0) {
427 S1 -= 1; I1 += 1; ell1 += 1;
428 ll += log(lambda)-log(I1*(I1-1)/2);
429 } else {
430 // the genealogy is incompatible with the state.
431 // nevertheless, the state remains valid.
432 I1 += 1; N1 += 1; ell1 += 1; // #nocov
433 ll += R_NegInf; // #nocov
434 }
435 } else { // s = (1,1)
436 if (unif_rand() < 0.5) {
437 color[lineage[c1]] = host1;
438 color[lineage[c2]] = host2;
439 } else {
440 color[lineage[c1]] = host2;
441 color[lineage[c2]] = host1;
442 }
443 ll -= log(0.5);
444 if (S2 > 0) {
445 S2 -= 1; I2 += 1; ell2 += 1;
446 ll += log(lambda)-log(I1*I2);
447 } else {
448 // the genealogy is incompatible with the state.
449 // nevertheless, the state remains valid.
450 I2 += 1; N2 += 1; ell2 += 1;
451 ll += R_NegInf;
452 }
453 }
454 } else if (parcol == host2) { // parent is in I2
455 assert(S1>=0 && S2 >=0 && I2>=ell2);
456 double lambda12 = Beta12*S1*I2/N2;
457 double lambda22 = Beta22*S2*I2/N2;
458 double lambda = lambda12+lambda22;
459 double x = (lambda > 0) ? lambda22/lambda : 0;
460 int c1 = child[index[parent]];
461 int c2 = child[index[parent]+1];
462 assert(c1 != c2);
463 assert(lineage[c1] != lineage[c2]);
464 assert(lineage[c1] != parlin || lineage[c2] != parlin);
465 assert(lineage[c1] == parlin || lineage[c2] == parlin);
466 if (unif_rand() < x) { // s = (0,2)
467 color[lineage[c1]] = host2;
468 color[lineage[c2]] = host2;
469 if (S2 > 0) {
470 S2 -= 1; I2 += 1; ell2 += 1;
471 ll += log(lambda)-log(I2*(I2-1)/2);
472 } else {
473 // the genealogy is incompatible with the state.
474 // nevertheless, the state remains valid.
475 I2 += 1; N2 += 1; ell2 += 1;
476 ll += R_NegInf;
477 }
478 } else { // s = (1,1)
479 if (unif_rand() < 0.5) {
480 color[lineage[c1]] = host1;
481 color[lineage[c2]] = host2;
482 } else {
483 color[lineage[c1]] = host2;
484 color[lineage[c2]] = host1;
485 }
486 ll -= log(0.5);
487 if (S1 > 0) {
488 S1 -= 1; I1 += 1; ell1 += 1;
489 ll += log(lambda)-log(I1*I2);
490 } else {
491 // the genealogy is incompatible with the state.
492 // nevertheless, the state remains valid.
493 I1 += 1; N1 += 1; ell1 += 1;
494 ll += R_NegInf;
495 }
496 }
497 } else {
498 assert(0); // #nocov
499 ll += R_NegInf; // #nocov
500 }
501 assert(nearbyint(N1)==nearbyint(S1+I1+R1));
502 assert(nearbyint(N2)==nearbyint(S2+I2+R2));
503 assert(check_color(color,nsample,ell1,ell2));
504 break;
505 }
506
507 // continuous portion of filter equation:
508 // take Gillespie steps to the end of the interval.
509 if (tmax > t) {
510
511 double rate[nrate], logpi[nrate];
512 int event;
513 double event_rate = 0;
514 double penalty = 0;
515
516 event_rate = EVENT_RATES;
517 tstep = exp_rand()/event_rate;
518
519 while (t + tstep < tmax) {
520 event = rcateg(event_rate,rate,nrate);
521 assert(event>=0 && event<nrate);
522 ll -= penalty*tstep + logpi[event];
523 switch (event) {
524 case 0: // 0: Trans_11, s = (0,0),(1,0)
525 assert(S1>=1 && I1>=0);
526 S1 -= 1; I1 += 1;
527 break;
528 case 1: // 1: Trans_22, s = (0,0),(0,1)
529 assert(S2>=1 && I2>=0);
530 S2 -= 1; I2 += 1;
531 break;
532 case 2: // 2: Trans_21, s = (0,0),(1,0)
533 assert(S2>=1 && I1>=0);
534 S2 -= 1; I2 += 1;
535 ll += log(1-ell2/I2);
536 assert(!ISNAN(ll));
537 break;
538 case 3: // 3: Trans_21, s = (0,1)
539 assert(S2>=1 && I1>=0);
540 S2 -= 1; I2 += 1;
541 change_color(color,nsample,random_choice(ell1),host1,host2);
542 ell1 -= 1; ell2 += 1;
543 ll += log(1-ell1/I1)-log(I2);
544 assert(check_color(color,nsample,ell1,ell2));
545 assert(!ISNAN(ll));
546 break;
547 case 4: // 4: Trans_12, s = (0,0),(0,1)
548 assert(S1>=1 && I2>=0);
549 S1 -= 1; I1 += 1;
550 ll += log(1-ell1/I1);
551 assert(!ISNAN(ll));
552 break;
553 case 5: // 5: Trans_12, s = (1,0)
554 assert(S1>=1 && I2>=0);
555 S1 -= 1; I1 += 1;
556 change_color(color,nsample,random_choice(ell2),host2,host1);
557 ell2 -= 1; ell1 += 1;
558 ll += log(1-ell2/I2)-log(I1);
559 assert(check_color(color,nsample,ell1,ell2));
560 assert(!ISNAN(ll));
561 break;
562 case 6: // 6: Recov_1
563 assert(I1>=1);
564 I1 -= 1; R1 += 1;
565 break;
566 case 7: // 7: Recov_2
567 assert(I2>=1);
568 I2 -= 1; R2 += 1;
569 break;
570 case 8: // 8: Wane_1
571 assert(R1>=1);
572 R1 -= 1; S1 += 1;
573 break;
574 case 9: // 9: Wane_2
575 assert(R2>=1);
576 R2 -= 1; S2 += 1;
577 break;
578 case 10: // 10: Birth_1
579 assert(N1>=1);
580 S1 += 1; N1 += 1;
581 break;
582 case 11: // 11: Birth_2
583 assert(N2>=1);
584 S2 += 1; N2 += 1;
585 break;
586 case 12: // 12: Death_S1
587 assert(S1>=1 && N1>=1);
588 S1 -= 1; N1 -= 1;
589 break;
590 case 13: // 13: Death_I1
591 assert(I1>=1 && N1>=1);
592 I1 -= 1; N1 -= 1;
593 break;
594 case 14: // 14: Death_R1
595 assert(R1>=1 && N1>=1);
596 R1 -= 1; N1 -= 1;
597 break;
598 case 15: // 15: Death_S2
599 assert(S2>=1 && N2>=1);
600 S2 -= 1; N2 -= 1;
601 break;
602 case 16: // 16: Death_I2
603 assert(I2>=1 && N2>=1);
604 I2 -= 1; N2 -= 1;
605 break;
606 case 17: // 17: Death_R2
607 assert(R2>=1 && N2>=1);
608 R2 -= 1; N2 -= 1;
609 break;
610 default: // #nocov
611 assert(0); // #nocov
612 ll += R_NegInf; // #nocov
613 break; // #nocov
614 }
615
616 ell1 = nearbyint(ell1);
617 ell2 = nearbyint(ell2);
618
619 assert(nearbyint(N1)==nearbyint(S1+I1+R1));
620 assert(nearbyint(N2)==nearbyint(S2+I2+R2));
621 assert(check_color(color,nsample,ell1,ell2));
622
623 t += tstep;
624 event_rate = EVENT_RATES;
625 tstep = exp_rand()/event_rate;
626
627 }
628 tstep = tmax - t;
629 ll -= penalty*tstep;
630 }
631 node += 1;
632}
ell1
#define ell1
Definition twospecies_pomp.c:73
change_color
static void change_color(double *color, int nsample, int n, int from, int to)
Definition twospecies_pomp.c:12
Beta21
#define Beta21
Definition twospecies_pomp.c:43
check_color
static int check_color(double *color, int nsample, double size1, double size2)
Definition twospecies_pomp.c:26
Beta12
#define Beta12
Definition twospecies_pomp.c:42
C1
#define C1
Definition twospecies_pomp.c:55
Beta22
#define Beta22
Definition twospecies_pomp.c:44
random_choice
static int random_choice(double n)
Definition twospecies_pomp.c:8
S1
#define S1
Definition twospecies_pomp.c:63
R2
#define R2
Definition twospecies_pomp.c:68
R1
#define R1
Definition twospecies_pomp.c:67
host1
#define host1
Definition twospecies_pomp.c:5
Beta11
#define Beta11
Definition twospecies_pomp.c:41
N1
#define N1
Definition twospecies_pomp.c:69
host2
#define host2
Definition twospecies_pomp.c:6
N2
#define N2
Definition twospecies_pomp.c:70
S2
#define S2
Definition twospecies_pomp.c:64
C2
#define C2
Definition twospecies_pomp.c:56
ell2
#define ell2
Definition twospecies_pomp.c:74
Here is the call graph for this function:

◆ twospecies_rinit()

void twospecies_rinit ( double * __x,
const double * __p,
double t0,
const int * __stateindex,
const int * __parindex,
const int * __covindex,
const double * __covars )
extern

Latent-state initializer (rinit component).

The state variables include S, E, I, R plus 'ellE' and 'ellI' (numbers of E- and I-deme lineages), the accumulated weight ('ll'), the current node number ('node'), and the coloring of each lineage ('COLOR').

Definition at line 249 of file twospecies_pomp.c.

258 {
259 double adj;
260 N1 = S1_0+I1_0+R1_0;
261 N2 = S2_0+I2_0+R2_0;
262 adj = N1/(S1_0+I1_0+R1_0);
263 S1 = nearbyint(S1_0*adj);
264 I1 = nearbyint(I1_0*adj);
265 R1 = N1-S1-I1;
266 adj = N2/(S2_0+I2_0+R2_0);
267 S2 = nearbyint(S2_0*adj);
268 I2 = nearbyint(I2_0*adj);
269 R2 = N2-S2-I2;
270 ell1 = 0;
271 ell2 = 0;
272 ll = 0;
273 node = 0;
274}
S2_0
#define S2_0
Definition twospecies_pomp.c:58
S1_0
#define S1_0
Definition twospecies_pomp.c:57
R2_0
#define R2_0
Definition twospecies_pomp.c:62
R1_0
#define R1_0
Definition twospecies_pomp.c:61