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 163 of file lbdp_pomp.c.

176 {
177 assert(!ISNAN(ll));
178 lik = (give_log) ? ll : exp(ll);
179}
lik
#define lik
Definition lbdp_pomp.c:160
ll
#define ll
Definition lbdp_pomp.c:9

◆ 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 71 of file lbdp_pomp.c.

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

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

◆ R_init_phylopomp()

void R_init_phylopomp ( DllInfo * info)
extern

Definition at line 52 of file init.c.

52 {
53 // Register routines
54 R_registerRoutines(info,NULL,callMethods,NULL,extMethods);
55 R_useDynamicSymbols(info,TRUE);
56 // R_useDynamicSymbols(info,FALSE);
57 // R_forceSymbols(info,TRUE);
58 get_userdata = (get_userdata_t*) R_GetCCallable("pomp","get_userdata");
59 get_userdata_double = (get_userdata_double_t*) R_GetCCallable("pomp","get_userdata_double");
60 get_userdata_int = (get_userdata_int_t*) R_GetCCallable("pomp","get_userdata_int");
61}
extMethods
static const R_CallMethodDef extMethods[]
Definition init.c:47
callMethods
static const R_CallMethodDef callMethods[]
Definition init.c:29
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 347 of file seirs_pomp.c.

360 {
361 assert(!ISNAN(ll));
362 lik = (give_log) ? ll : exp(ll);
363}

◆ 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 139 of file seirs_pomp.c.

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

123 {
124 double adj = N/(S0+E0+I0+R0);
125 S = nearbyint(S0*adj);
126 E = nearbyint(E0*adj);
127 I = nearbyint(I0*adj);
128 R = nearbyint(R0*adj);
129 ellE = 0;
130 ellI = 0;
131 ll = 0;
132 node = 0;
133}
R0
#define R0
Definition seirs_pomp.c:31
S0
#define S0
Definition seirs_pomp.c:28
I0
#define I0
Definition seirs_pomp.c:30
E0
#define E0
Definition seirs_pomp.c:29

◆ 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
118 break;
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 290 of file strains_pomp.c.

303 {
304 assert(!ISNAN(ll));
305 lik = (give_log) ? ll : exp(ll);
306}

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