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 bdei_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 
void bdei_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void bdei_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 bdss_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 
void bdss_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void bdss_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 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 si2rs_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 
void si2rs_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void si2rs_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

◆ bdei_dmeas()

void bdei_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 320 of file bdei_pomp.c.

333 {
334 assert(!ISNAN(ll));
335 lik = (give_log) ? ll : exp(ll);
336}
lik
#define lik
Definition bdei_pomp.c:317
ll
#define ll
Definition bdei_pomp.c:35

◆ bdei_gill()

void bdei_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 BDEI process.

Definition at line 138 of file bdei_pomp.c.

148 {
149 double tstep = 0.0, tmax = t + dt;
150 double *color = &COLOR;
151 const int nsample = *get_userdata_int("nsample");
152 const int *nodetype = get_userdata_int("nodetype");
153 const int *lineage = get_userdata_int("lineage");
154 const int *sat = get_userdata_int("saturation");
155 const int *index = get_userdata_int("index");
156 const int *child = get_userdata_int("child");
157
158 int parent = (int) nearbyint(node);
159
160#ifndef NDEBUG
161 int nnode = *get_userdata_int("nnode");
162 assert(parent>=0);
163 assert(parent<=nnode);
164#endif
165
166 int parlin = lineage[parent];
167 int parcol = color[parlin];
168 assert(parlin >= 0 && parlin < nsample);
169
170 ll = 0;
171
172 // singular portion of filter equation
173 switch (nodetype[parent]) {
174 default:
175 break;
176 case 0: // root
177 assert(sat[parent]==1);
178 {
179 int c = child[index[parent]];
180 assert(lineage[parent]==lineage[c]);
181 if (E-ellE + I-ellI > 0) {
182 double x = (E-ellE)/(E-ellE + I-ellI);
183 if (unif_rand() < x) {
184 color[lineage[c]] = Exposed;
185 ellE += 1;
186 ll -= log(x);
187 } else {
188 color[lineage[c]] = Infected;
189 ellI += 1;
190 ll -= log(1-x);
191 }
192 } else {
193 ll += R_NegInf;
194 if (unif_rand() < 0.5) {
195 color[lineage[c]] = Exposed;
196 ellE += 1; E += 1;
197 } else {
198 color[lineage[c]] = Infected;
199 ellI += 1; I += 1;
200 }
201 }
202 }
203 break;
204 case 1: // sample
205 if (parcol != Infected) {
206 ll += R_NegInf;
207 color[parlin] = Infected;
208 ellE -= 1; ellI += 1;
209 E -= 1; I += 1;
210 }
211 if (sat[parent] == 0) { // destructive sample
212 ellI -= 1;
213 ll += (I > 0) ? log(chi*I) : R_NegInf;
214 I -= 1;
215 } else {
216 assert(0); // #nocov
217 ll += R_NegInf; // #nocov
218 }
219 color[parlin] = R_NaReal;
220 break;
221 case 2: // branch point: I births E
222 if (parcol != Infected) {
223 ll += R_NegInf;
224 color[parlin] = Infected;
225 ellE -= 1; ellI += 1;
226 E -= 1; I += 1;
227 }
228 assert(sat[parent]==2);
229 ll += (I > 0) ? log(lambda/(E+1)) : R_NegInf;
230 E += 1;
231 ellE += 1;
232 {
233 int c1 = child[index[parent]];
234 int c2 = child[index[parent]+1];
235 assert(c1 != c2);
236 assert(lineage[c1] != lineage[c2]);
237 assert(lineage[c1] != parlin || lineage[c2] != parlin);
238 assert(lineage[c1] == parlin || lineage[c2] == parlin);
239 if (unif_rand() < 0.5) {
240 color[lineage[c1]] = Exposed;
241 color[lineage[c2]] = Infected;
242 } else {
243 color[lineage[c1]] = Infected;
244 color[lineage[c2]] = Exposed;
245 }
246 ll -= log(0.5);
247 }
248 break;
249 }
250
251 // regular portion of filter equation
252 if (tmax > t && R_FINITE(ll)) {
253
254 double rate[nrate], logpi[nrate];
255 int event;
256 double event_rate = 0;
257 double penalty = 0;
258
259 event_rate = EVENT_RATES;
260 tstep = exp_rand()/event_rate;
261
262 while (t + tstep < tmax) {
263 event = rcateg(event_rate,rate,nrate);
264 assert(event>=0 && event<nrate);
265 ll -= penalty*tstep + logpi[event];
266 switch (event) {
267 case 0: // birth, s=(0,0) or s=(0,1)
268 E += 1;
269 ll += log(1-ellE/E);
270 assert(!ISNAN(ll));
271 break;
272 case 1: // birth, s=(1,0)
273 change_color(color,nsample,random_choice(ellI),Infected,Exposed);
274 ellE += 1; ellI -= 1;
275 E += 1;
276 ll += log(1-ellI/I)-log(E);
277 assert(!ISNAN(ll));
278 break;
279 case 2: // progression, s=(0,0)
280 assert(E>=1);
281 E -= 1; I += 1;
282 ll += log(1-ellI/I);
283 assert(!ISNAN(ll));
284 break;
285 case 3: // progression, s=(0,1)
286 assert(E>=1);
287 change_color(color,nsample,random_choice(ellE),Exposed,Infected);
288 ellE -= 1; ellI += 1;
289 E -= 1; I += 1;
290 ll -= log(I);
291 assert(!ISNAN(ll));
292 break;
293 case 4: // death
294 assert(I>=1);
295 I -= 1;
296 assert(!ISNAN(ll));
297 break;
298 default: // #nocov
299 assert(0); // #nocov
300 break; // #nocov
301 }
302
303 ellE = nearbyint(ellE);
304 ellI = nearbyint(ellI);
305
306 t += tstep;
307 event_rate = EVENT_RATES;
308 tstep = exp_rand()/event_rate;
309
310 }
311 tstep = tmax - t;
312 ll -= penalty*tstep;
313 }
314 node += 1;
315}
ellE
#define ellE
Definition bdei_pomp.c:37
E
#define E
Definition bdei_pomp.c:33
change_color
static void change_color(double *color, int nsample, int n, int from, int to)
Definition bdei_pomp.c:13
COLOR
#define COLOR
Definition bdei_pomp.c:39
random_choice
static int random_choice(double n)
Definition bdei_pomp.c:9
I
#define I
Definition bdei_pomp.c:34
ellI
#define ellI
Definition bdei_pomp.c:38
chi
#define chi
Definition bdei_pomp.c:29
nrate
static const int nrate
Definition bdei_pomp.c:7
Infected
#define Infected
Definition bdei_pomp.c:5
lambda
#define lambda
Definition bdei_pomp.c:27
Exposed
#define Exposed
Definition bdei_pomp.c:4
EVENT_RATES
#define EVENT_RATES
Definition bdei_pomp.c:41
node
#define node
Definition bdei_pomp.c:36
nsample
SEXP nsample(TYPE &X)
Definition generics.h:12
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
Here is the call graph for this function:

◆ bdei_rinit()

void bdei_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 E, I 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 106 of file bdei_pomp.c.

115 {
116 if (!R_FINITE(sigma) || sigma < 0 ||
117 !R_FINITE(lambda) || lambda < 0 ||
118 !R_FINITE(mu) || mu < 0 || !R_FINITE(chi) || chi < 0)
119 err("BDEI rates must be finite and nonnegative.");
120 if (!R_FINITE(POP) || POP < 0 ||
121 !R_FINITE(E0) || E0 < 0 || !R_FINITE(I0) || I0 < 0 ||
122 !R_FINITE(E0+I0) || E0+I0 <= 0)
123 err("BDEI initial values must be finite and nonnegative, with E0+I0 > 0.");
124 double adj = POP/(E0+I0);
125 if (!R_FINITE(adj)) err("BDEI: pop/(E0+I0) must be finite.");
126 E = nearbyint(E0*adj);
127 I = nearbyint(I0*adj);
128 ellE = 0;
129 ellI = 0;
130 ll = 0;
131 node = 0;
132}
POP
#define POP
Definition bdei_pomp.c:30
mu
#define mu
Definition bdei_pomp.c:28
sigma
#define sigma
Definition bdei_pomp.c:26
I0
#define I0
Definition bdei_pomp.c:32
E0
#define E0
Definition bdei_pomp.c:31
err
#define err(...)
Definition internal.h:18

◆ bdss_dmeas()

void bdss_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 392 of file bdss_pomp.c.

405 {
406 assert(!ISNAN(ll));
407 lik = (give_log) ? ll : exp(ll);
408}

◆ bdss_gill()

void bdss_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 BDSS process.

Definition at line 162 of file bdss_pomp.c.

172 {
173 double tstep = 0.0, tmax = t + dt;
174 double *color = &COLOR;
175 const int nsample = *get_userdata_int("nsample");
176 const int *nodetype = get_userdata_int("nodetype");
177 const int *lineage = get_userdata_int("lineage");
178 const int *sat = get_userdata_int("saturation");
179 const int *index = get_userdata_int("index");
180 const int *child = get_userdata_int("child");
181
182 int parent = (int) nearbyint(node);
183
184#ifndef NDEBUG
185 int nnode = *get_userdata_int("nnode");
186 assert(parent>=0);
187 assert(parent<=nnode);
188#endif
189
190 int parlin = lineage[parent];
191 int parcol = color[parlin];
192 assert(parlin >= 0 && parlin < nsample);
193
194 ll = 0;
195
196 // singular portion of filter equation
197 switch (nodetype[parent]) {
198 default:
199 break;
200 case 0: // root
201 assert(sat[parent]==1);
202 {
203 int c = child[index[parent]];
204 assert(lineage[parent]==lineage[c]);
205 if (N-ell1 + S-ell2 > 0) {
206 double x = (N-ell1)/(N-ell1 + S-ell2);
207 if (unif_rand() < x) {
208 color[lineage[c]] = Normal;
209 ell1 += 1;
210 ll -= log(x);
211 } else {
212 color[lineage[c]] = Superspreader;
213 ell2 += 1;
214 ll -= log(1-x);
215 }
216 } else {
217 ll += R_NegInf;
218 if (unif_rand() < 0.5) {
219 color[lineage[c]] = Normal;
220 ell1 += 1; N += 1;
221 } else {
222 color[lineage[c]] = Superspreader;
223 ell2 += 1; S += 1;
224 }
225 }
226 }
227 break;
228 case 1: // sample (destructive only)
229 if (sat[parent] == 0) {
230 if (parcol == Normal) {
231 ell1 -= 1;
232 ll += (N > 0) ? log(chi*N) : R_NegInf;
233 N -= 1;
234 } else if (parcol == Superspreader) {
235 ell2 -= 1;
236 ll += (S > 0) ? log(chi*S) : R_NegInf;
237 S -= 1;
238 } else {
239 ll += R_NegInf;
240 }
241 } else {
242 assert(0); // #nocov
243 ll += R_NegInf; // #nocov
244 }
245 color[parlin] = R_NaReal;
246 break;
247 case 2: // branch point
248 assert(sat[parent]==2);
249 {
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
257 if (parcol == Normal) {
258 double lam = lambda_nn + lambda_ns;
259 if (lam <= 0) { ll += R_NegInf; break; }
260 double p_nn = lambda_nn/lam;
261 if (unif_rand() < p_nn) {
262 // N->N birth: both children Normal
263 N += 1; ell1 += 1;
264 ll += log(2*lam/N);
265 if (N*(N-1) <= 0) ll += R_NegInf;
266 color[lineage[c1]] = Normal;
267 color[lineage[c2]] = Normal;
268 } else {
269 // N->S birth: one N, one S
270 S += 1; ell2 += 1;
271 ll += log(lam/S);
272 if (N*S <= 0) ll += R_NegInf;
273 if (unif_rand() < 0.5) {
274 color[lineage[c1]] = Normal;
275 color[lineage[c2]] = Superspreader;
276 } else {
277 color[lineage[c1]] = Superspreader;
278 color[lineage[c2]] = Normal;
279 }
280 ll -= log(0.5);
281 }
282 } else if (parcol == Superspreader) {
283 double lam = lambda_sn + lambda_ss;
284 if (lam <= 0) { ll += R_NegInf; break; }
285 double p_ss = lambda_ss/lam;
286 if (unif_rand() < p_ss) {
287 // S->S birth: both children Superspreader
288 S += 1; ell2 += 1;
289 ll += log(2*lam/S);
290 if (S*(S-1) <= 0) ll += R_NegInf;
291 color[lineage[c1]] = Superspreader;
292 color[lineage[c2]] = Superspreader;
293 } else {
294 // S->N birth: one S, one N
295 N += 1; ell1 += 1;
296 ll += log(lam/N);
297 if (N*S <= 0) ll += R_NegInf;
298 if (unif_rand() < 0.5) {
299 color[lineage[c1]] = Normal;
300 color[lineage[c2]] = Superspreader;
301 } else {
302 color[lineage[c1]] = Superspreader;
303 color[lineage[c2]] = Normal;
304 }
305 ll -= log(0.5);
306 }
307 } else {
308 ll += R_NegInf;
309 }
310 }
311 break;
312 }
313
314 // regular portion of filter equation
315 if (tmax > t && R_FINITE(ll)) {
316
317 double rate[nrate], logpi[nrate];
318 int event;
319 double event_rate = 0;
320 double penalty = 0;
321
322 event_rate = EVENT_RATES;
323 tstep = exp_rand()/event_rate;
324
325 while (t + tstep < tmax) {
326 event = rcateg(event_rate,rate,nrate);
327 assert(event>=0 && event<nrate);
328 ll -= penalty*tstep + logpi[event];
329 switch (event) {
330 case 0: // N->N birth
331 N += 1;
332 assert(!ISNAN(ll));
333 break;
334 case 1: // N->S birth, s=(0,0) or s=(0,1)
335 S += 1;
336 ll += log(1-ell2/S);
337 assert(!ISNAN(ll));
338 break;
339 case 2: // N->S birth, s=(1,0)
340 change_color(color,nsample,random_choice(ell1),Normal,Superspreader);
341 ell1 -= 1; ell2 += 1;
342 S += 1;
343 ll += log(1-ell1/N)-log(S);
344 assert(!ISNAN(ll));
345 break;
346 case 3: // S->S birth
347 S += 1;
348 assert(!ISNAN(ll));
349 break;
350 case 4: // S->N birth, s=(0,0) or s=(1,0)
351 N += 1;
352 ll += log(1-ell1/N);
353 assert(!ISNAN(ll));
354 break;
355 case 5: // S->N birth, s=(0,1)
356 change_color(color,nsample,random_choice(ell2),Superspreader,Normal);
357 ell2 -= 1; ell1 += 1;
358 N += 1;
359 ll += log(1-ell2/S)-log(N);
360 assert(!ISNAN(ll));
361 break;
362 case 6: // death N
363 assert(N>=1);
364 N -= 1;
365 break;
366 case 7: // death S
367 assert(S>=1);
368 S -= 1;
369 break;
370 default: // #nocov
371 assert(0); // #nocov
372 break; // #nocov
373 }
374
375 ell1 = nearbyint(ell1);
376 ell2 = nearbyint(ell2);
377
378 t += tstep;
379 event_rate = EVENT_RATES;
380 tstep = exp_rand()/event_rate;
381
382 }
383 tstep = tmax - t;
384 ll -= penalty*tstep;
385 }
386 node += 1;
387}
N
#define N
Definition bdss_pomp.c:35
ell1
#define ell1
Definition bdss_pomp.c:39
change_color
static void change_color(double *color, int nsample, int n, int from, int to)
Definition bdss_pomp.c:13
lambda_ns
#define lambda_ns
Definition bdss_pomp.c:27
random_choice
static int random_choice(double n)
Definition bdss_pomp.c:9
Normal
#define Normal
Definition bdss_pomp.c:4
lambda_nn
#define lambda_nn
Definition bdss_pomp.c:26
Superspreader
#define Superspreader
Definition bdss_pomp.c:5
lambda_sn
#define lambda_sn
Definition bdss_pomp.c:28
ell2
#define ell2
Definition bdss_pomp.c:40
S
#define S
Definition bdss_pomp.c:36
lambda_ss
#define lambda_ss
Definition bdss_pomp.c:29
Here is the call graph for this function:

◆ bdss_rinit()

void bdss_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 N, S plus 'ell1' and 'ell2' (numbers of N- and S-deme lineages), the accumulated weight ('ll'), the current node number ('node'), and the coloring of each lineage ('COLOR').

Definition at line 128 of file bdss_pomp.c.

137 {
138 if (!R_FINITE(lambda_nn) || lambda_nn < 0 ||
139 !R_FINITE(lambda_ns) || lambda_ns < 0 ||
140 !R_FINITE(lambda_sn) || lambda_sn < 0 ||
141 !R_FINITE(lambda_ss) || lambda_ss < 0 ||
142 !R_FINITE(mu) || mu < 0 || !R_FINITE(chi) || chi < 0)
143 err("BDSS rates must be finite and nonnegative.");
144 if (!R_FINITE(POP) || POP < 0 ||
145 !R_FINITE(N0) || N0 < 0 || !R_FINITE(S0) || S0 < 0 ||
146 !R_FINITE(N0+S0) || N0+S0 <= 0)
147 err("BDSS initial values must be finite and nonnegative, with N0+S0 > 0.");
148 double adj = POP/(N0+S0);
149 if (!R_FINITE(adj)) err("BDSS: pop/(N0+S0) must be finite.");
150 N = nearbyint(N0*adj);
151 S = nearbyint(S0*adj);
152 ell1 = 0;
153 ell2 = 0;
154 ll = 0;
155 node = 0;
156}
N0
#define N0
Definition bdss_pomp.c:33
S0
#define S0
Definition bdss_pomp.c:34

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

181 {
182 assert(!ISNAN(ll));
183 lik = (give_log) ? ll : exp(ll);
184}

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

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

◆ R_init_phylopomp()

void R_init_phylopomp ( DllInfo * info)
extern

Definition at line 70 of file init.c.

70 {
71 // Register routines
72 R_registerRoutines(info,NULL,callMethods,NULL,extMethods);
73 R_useDynamicSymbols(info,TRUE);
74 // R_useDynamicSymbols(info,FALSE);
75 // R_forceSymbols(info,TRUE);
76 get_userdata = (get_userdata_t*) R_GetCCallable("pomp","get_userdata");
77 get_userdata_double = (get_userdata_double_t*) R_GetCCallable("pomp","get_userdata_double");
78 get_userdata_int = (get_userdata_int_t*) R_GetCCallable("pomp","get_userdata_int");
79}
extMethods
static const R_CallMethodDef extMethods[]
Definition init.c:64
callMethods
static const R_CallMethodDef callMethods[]
Definition init.c:38
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 349 of file seirs_pomp.c.

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

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

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

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

◆ si2rs_dmeas()

void si2rs_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 403 of file si2r_pomp.c.

416 {
417 assert(!ISNAN(ll));
418 lik = (give_log) ? ll : exp(ll);
419}

◆ si2rs_gill()

void si2rs_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 SI2RS process.

Definition at line 183 of file si2r_pomp.c.

193 {
194 double tstep = 0.0, tmax = t + dt;
195 double *color = &COLOR;
196 const int nsample = *get_userdata_int("nsample");
197 const int *nodetype = get_userdata_int("nodetype");
198 const int *lineage = get_userdata_int("lineage");
199 const int *index = get_userdata_int("index");
200 const int *child = get_userdata_int("child");
201
202 int parent = (int) nearbyint(node);
203
204#ifndef NDEBUG
205 const int *sat = get_userdata_int("saturation");
206 int nnode = *get_userdata_int("nnode");
207 assert(parent>=0);
208 assert(parent<=nnode);
209#endif
210
211 int parlin = lineage[parent];
212 int parcol = color[parlin];
213 assert(parlin >= 0 && parlin < nsample);
214
215 ll = 0;
216
217 // singular portion of filter equation
218 switch (nodetype[parent]) {
219 default: // non-genealogical event #nocov
220 break; // #nocov
221 case 0: // root
222 // color lineages by sampling without replacement
223 assert(sat[parent]==1);
224 int c = child[index[parent]];
225 assert(lineage[parent]==lineage[c]);
226 if (IL-ellL + IH-ellH > 0) {
227 double x = (IL-ellL)/(IL-ellL + IH-ellH);
228 if (unif_rand() < x) { // lineage is put into Low deme
229 color[lineage[c]] = Low;
230 ellL += 1;
231 ll -= log(x);
232 } else { // lineage is put into High deme
233 color[lineage[c]] = High;
234 ellH += 1;
235 ll -= log(1-x);
236 }
237 assert(!ISNAN(ll));
238 } else { // more roots than infectives
239 ll += R_NegInf; // this is incompatible with the genealogy
240 // the following keeps the state valid
241 color[lineage[c]] = Low;
242 ellL += 1; IL += 1;
243 }
244 break;
245 case 1: // sample
246 assert(sat[parent]==0);
247 if (parcol == Low) {
248 assert(ellL>=1 && IL >= ellL);
249 ll += log(chi*IL);
250 ellL -= 1; IL -= 1; R += 1;
251 } else if (parcol == High) {
252 assert(ellH>=1 && IH >= ellH);
253 ll += log(chi*IH);
254 ellH -= 1; IH -= 1; R += 1;
255 } else {
256 assert(0); // #nocov
257 }
258 color[parlin] = R_NaReal;
259 break;
260 case 2:
261 assert(sat[parent]==2);
262 int c1 = child[index[parent]];
263 int c2 = child[index[parent]+1];
264 assert(c1 != c2);
265 assert(lineage[c1] != lineage[c2]);
266 assert(lineage[c1] != parlin || lineage[c2] != parlin);
267 assert(lineage[c1] == parlin || lineage[c2] == parlin);
268 if (parcol == Low) {
269 assert(ellL >= 1 && IL >= ellL);
270 if (S >= 1 && POP > 0) {
271 ll += log(Beta*S*IL/POP);
272 S -= 1; IL += 1;
273 ellL += 1;
274 color[lineage[c1]] = Low;
275 color[lineage[c2]] = Low;
276 ll -= log(IL*(IL-1)/2);
277 } else {
278 ll += R_NegInf;
279 IL += 1; ellL += 1;
280 color[lineage[c1]] = Low;
281 color[lineage[c2]] = Low;
282 }
283 } else if (parcol == High) {
284 assert(ellH >= 1 && IH >= ellH);
285 if (S>=1 && POP > 0) {
286 ll += log(kappa*Beta*S*IH/POP);
287 S -= 1; IL += 1;
288 ellL += 1;
289 if (unif_rand() < 0.5) {
290 color[lineage[c1]] = Low;
291 color[lineage[c2]] = High;
292 } else {
293 color[lineage[c1]] = High;
294 color[lineage[c2]] = Low;
295 }
296 ll -= log(IL*IH/2);
297 assert(!ISNAN(ll));
298 } else {
299 ll += R_NegInf;
300 IL += 1; ellL += 1;
301 color[lineage[c1]] = Low;
302 color[lineage[c2]] = High;
303 }
304 } else {
305 assert(0); // #nocov
306 }
307 break;
308 }
309
310 // continuous portion of filter equation:
311 // take Gillespie steps to the end of the interval
312 if (tmax > t && R_FINITE(ll)) {
313
314 double rate[nrate], logpi[nrate];
315 int event;
316 double event_rate = 0;
317 double decay = 0;
318
319 event_rate = EVENT_RATES;
320 tstep = exp_rand()/event_rate;
321
322 while (t + tstep < tmax) {
323 event = rcateg(event_rate,rate,nrate);
324 assert(event>=0 && event<nrate);
325 ll -= decay*tstep + logpi[event];
326 switch (event) {
327 case 0: // TL, s=(0,0) or s=(1,0)
328 assert(S>=1 && IL>=1);
329 S -= 1; IL += 1;
330 ll += log(1-ellL*(ellL-1)/IL/(IL-1));
331 break;
332 case 1: // TH, s=(0,0) or s=(0,1)
333 assert(S>=1 && IH >= 1);
334 S -= 1; IL += 1;
335 ll += log(1-ellL/IL);
336 break;
337 case 2: // TH, s=(1,0)
338 assert(S>=1 && IH >= 1);
339 S -= 1; IL += 1;
340 change_color(color,nsample,random_choice(ellH),High,Low);
341 ellH -= 1; ellL += 1;
342 ll += log((1-ellH/IH)/IL);
343 break;
344 case 3: // L, s=(0,0)
345 assert(IL>=1);
346 IL -= 1; IH += 1;
347 ll += log(1-ellH/IH);
348 break;
349 case 4: // L, s=(0,1)
350 assert(IL>=1);
351 change_color(color,nsample,random_choice(ellL),Low,High);
352 ellL -= 1; ellH += 1;
353 IL -= 1; IH += 1;
354 ll -= log(IH);
355 break;
356 case 5: // H, s=(0,0)
357 assert(IH>=1);
358 IL += 1; IH -= 1;
359 ll += log(1-ellL/IL);
360 break;
361 case 6: // H, s=(1,0)
362 assert(IH>=1);
363 change_color(color,nsample,random_choice(ellH),High,Low);
364 ellL += 1; ellH -= 1;
365 IL += 1; IH -= 1;
366 ll -= log(IL);
367 break;
368 case 7: // RL
369 assert(IL>=1);
370 IL -= 1; R += 1;
371 break;
372 case 8: // RH
373 assert(IH>=1);
374 IH -= 1; R += 1;
375 break;
376 case 9: // W
377 assert(R>=1);
378 R -= 1; S += 1;
379 break;
380 default: // #nocov
381 assert(0); // #nocov
382 ll += R_NegInf; // #nocov
383 break; // #nocov
384 }
385
386 ellL = nearbyint(ellL);
387 ellH = nearbyint(ellH);
388
389 t += tstep;
390 event_rate = EVENT_RATES;
391 tstep = exp_rand()/event_rate;
392
393 }
394 tstep = tmax - t;
395 ll -= decay*tstep;
396 }
397 node += 1;
398}
kappa
#define kappa
Definition si2r_pomp.c:27
ellH
#define ellH
Definition si2r_pomp.c:45
change_color
static void change_color(double *color, int nsample, int n, int from, int to)
Definition si2r_pomp.c:13
High
#define High
Definition si2r_pomp.c:5
IH
#define IH
Definition si2r_pomp.c:40
ellL
#define ellL
Definition si2r_pomp.c:44
random_choice
static int random_choice(double n)
Definition si2r_pomp.c:9
Low
#define Low
Definition si2r_pomp.c:4
IL
#define IL
Definition si2r_pomp.c:39
Here is the call graph for this function:

◆ si2rs_rinit()

void si2rs_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, IL, IH, R plus 'ellL' and 'ellH' (numbers of L- and H-deme lineages), the accumulated weight ('ll'), the current node number ('node'), and the coloring of each lineage ('COLOR').

Definition at line 158 of file si2r_pomp.c.

167 {
168 double adj = POP/(S0+IL0+IH0+R0);
169 S = nearbyint(S0*adj);
170 IL = nearbyint(IL0*adj);
171 IH = nearbyint(IH0*adj);
172 R = nearbyint(R0*adj);
173 ellL = 0;
174 ellH = 0;
175 ll = 0;
176 node = 0;
177}
IL0
#define IL0
Definition si2r_pomp.c:35
IH0
#define IH0
Definition si2r_pomp.c:36

◆ 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 206 of file sirs_pomp.c.

219 {
220 assert(!ISNAN(ll));
221 lik = (give_log) ? ll : exp(ll);
222}

◆ 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 91 of file sirs_pomp.c.

101 {
102 double tstep = 0.0, tmax = t + dt;
103 const int *nodetype = get_userdata_int("nodetype");
104 const int *sat = get_userdata_int("saturation");
105
106 int parent = (int) nearbyint(node);
107
108#ifndef NDEBUG
109 int nnode = *get_userdata_int("nnode");
110 assert(parent>=0);
111 assert(parent<=nnode);
112#endif
113
114 ll = 0;
115
116 // singular portion of filter equation
117 switch (nodetype[parent]) {
118 default: // non-genealogical event #nocov
119 break; // #nocov
120 case 0: // root
121 if (I - ellI > 0)
122 ellI += 1;
123 else {
124 ellI += 1;
125 I += 1;
126 ll += R_NegInf;
127 }
128 break;
129 case 1: // sample
130 assert(I >= ellI);
131 assert(ellI >= 0);
132 if (sat[parent] == 1) {
133 ll += log(psi);
134 } else if (sat[parent] == 0) {
135 if (psi+chi > 0) {
136 ll += log((psi+chi)*I);
137 if (unif_rand() < psi/(psi+chi)) { // non-destructive sample
138 ellI -= 1;
139 ll += log(1-ellI/I);
140 } else { // destructive sample
141 ellI -= 1; I -= 1; R += 1;
142 }
143 } else {
144 ll += R_NegInf;
145 ellI -= 1; I -= 1; R += 1;
146 }
147 } else {
148 assert(0); // #nocov
149 ll += R_NegInf; // #nocov
150 }
151 break;
152 case 2: // branch point s=(1,1)
153 assert(S >= 0);
154 assert(I >= 0);
155 assert(ellI > 0);
156 assert(sat[parent]==2);
157 ll += (I > 0 && I >= ellI) ? log(Beta*S*I/POP) : R_NegInf;
158 S -= 1; I += 1;
159 ellI += 1;
160 ll -= log(I*(I-1)/2);
161 S = (S > 0) ? S : 0;
162 break;
163 }
164
165 if (tmax > t) {
166
167 // take Gillespie steps to the end of the interval:
168 int event;
169 double decay = 0;
170 double rate[nrate];
171
172 double event_rate = EVENT_RATES;
173 tstep = exp_rand()/event_rate;
174
175 while (t + tstep < tmax) {
176 event = rcateg(event_rate,rate,nrate);
177 assert(event>=0 && event<nrate);
178 ll -= decay*tstep;
179 switch (event) {
180 case 0: // transmission
181 S -= 1; I += 1;
182 break;
183 case 1: // recovery
184 I -= 1; R += 1;
185 break;
186 case 2: // loss of immunity
187 R -= 1; S += 1;
188 break;
189 default: // #nocov
190 assert(0); // #nocov
191 break; // #nocov
192 }
193 t += tstep;
194 event_rate = EVENT_RATES;
195 tstep = exp_rand()/event_rate;
196 }
197 tstep = tmax - t;
198 ll -= decay*tstep;
199 }
200 node += 1;
201}
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 69 of file sirs_pomp.c.

78 {
79 double m = POP/(S0+I0+R0);
80 S = nearbyint(S0*m);
81 I = nearbyint(I0*m);
82 R = nearbyint(R0*m);
83 ll = 0;
84 ellI = 0;
85 node = 0;
86}

◆ 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 136 of file strains_pomp.c.

146 {
147 double tstep = 0.0, tmax = t + dt;
148 const int *nodetype = get_userdata_int("nodetype");
149 const int *sat = get_userdata_int("saturation");
150 const int *deme = get_userdata_int("deme");
151 const int *index = get_userdata_int("index");
152 const int *child = get_userdata_int("child");
153
154 int parent = (int) nearbyint(node);
155 int c = child[index[parent]];
156
157#ifndef NDEBUG
158 const int *lineage = get_userdata_int("lineage");
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 #nocov
169 break; // #nocov
170 case 0: // root
171 assert(sat[parent]==1);
172 assert(lineage[parent]==lineage[c]);
173 switch (deme[c]) {
174 case STRAIN1:
175 ellI1 += 1; break;
176 case STRAIN2:
177 ellI2 += 1; break;
178 case STRAIN3:
179 ellI3 += 1; break;
180 default: // #nocov
181 assert(0); break; // #nocov
182 }
183 break;
184 case 1: // sample
185 assert(sat[parent]==0);
186 switch (deme[parent]) {
187 case STRAIN1:
188 assert(I1 >= ellI1);
189 assert(ellI1 >= 0);
190 ll += log(chi*I1);
191 ellI1 -= 1; I1 -= 1;
192 break;
193 case STRAIN2:
194 assert(I2 >= ellI2);
195 assert(ellI2 >= 0);
196 ll += log(chi*I2);
197 ellI2 -= 1; I2 -= 1;
198 break;
199 case STRAIN3:
200 assert(I3 >= ellI3);
201 assert(ellI3 >= 0);
202 ll += log(chi*I3);
203 ellI3 -= 1; I3 -= 1;
204 break;
205 default: // #nocov
206 assert(0); break; // #nocov
207 }
208 break;
209 case 2: // branch point s=(1,1)
210 assert(S >= 0);
211 if (sat[parent]!=2) break;
212 switch (deme[parent]) {
213 case STRAIN1:
214 assert(I1 >= 0);
215 assert(ellI1 > 0);
216 ll += (I1 > 0 && I1 >= ellI1) ? log(Beta1*S*I1/POP) : R_NegInf;
217 S -= 1; I1 += 1; ellI1 += 1;
218 ll -= log(I1*(I1-1)/2);
219 break;
220 case STRAIN2:
221 assert(I2 >= 0);
222 assert(ellI2 > 0);
223 ll += (I2 > 0 && I2 >= ellI2) ? log(Beta2*S*I2/POP) : R_NegInf;
224 S -= 1; I2 += 1; ellI2 += 1;
225 ll -= log(I2*(I2-1)/2);
226 break;
227 case STRAIN3:
228 assert(I3 >= 0);
229 assert(ellI3 > 0);
230 ll += (I3 > 0 && I3 >= ellI3) ? log(Beta3*S*I3/POP) : R_NegInf;
231 S -= 1; I3 += 1; ellI3 += 1;
232 ll -= log(I3*(I3-1)/2);
233 break;
234 default: // #nocov
235 assert(0); break; // #nocov
236 }
237 S = (S > 0) ? S : 0;
238 break;
239 }
240
241 if (tmax > t && R_FINITE(ll)) {
242
243 // take Gillespie steps to the end of the interval:
244 int event;
245 double penalty = 0;
246 double rate[nrate], logpi[nrate];
247
248 double event_rate = EVENT_RATES;
249 tstep = exp_rand()/event_rate;
250
251 while (t + tstep < tmax) {
252 event = rcateg(event_rate,rate,nrate);
253 assert(event>=0 && event<nrate);
254 ll -= penalty*tstep + logpi[event];
255 switch (event) {
256 case 0: // strain-1 transmission
257 S -= 1; I1 += 1;
258 break;
259 case 1: // strain-1 recovery
260 I1 -= 1; R += 1;
261 break;
262 case 2: // strain-2 transmission
263 S -= 1; I2 += 1;
264 break;
265 case 3: // strain-2 recovery
266 I2 -= 1; R += 1;
267 break;
268 case 4: // strain-3 transmission
269 S -= 1; I3 += 1;
270 break;
271 case 5: // strain-3 recovery
272 I3 -= 1; R += 1;
273 break;
274 default: // #nocov
275 assert(0); 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}
deme
static const int deme
Definition lbdp.cc:7
Beta3
#define Beta3
Definition strains_pomp.c:12
STRAIN3
#define STRAIN3
Definition strains_pomp.c:8
ellI1
#define ellI1
Definition strains_pomp.c:27
ellI2
#define ellI2
Definition strains_pomp.c:28
ellI3
#define ellI3
Definition strains_pomp.c:29
STRAIN1
#define STRAIN1
Definition strains_pomp.c:6
STRAIN2
#define STRAIN2
Definition strains_pomp.c:7
I2
#define I2
Definition strains_pomp.c:23
Beta2
#define Beta2
Definition strains_pomp.c:11
Beta1
#define Beta1
Definition strains_pomp.c:10
I3
#define I3
Definition strains_pomp.c:24
I1
#define I1
Definition strains_pomp.c:22
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 110 of file strains_pomp.c.

119 {
120 double m = POP/(S_0+I1_0+I2_0+I3_0+R_0);
121 S = nearbyint(S_0*m);
122 I1 = nearbyint(I1_0*m);
123 I2 = nearbyint(I2_0*m);
124 I3 = nearbyint(I3_0*m);
125 R = nearbyint(R_0*m);
126 ll = 0;
127 ellI1 = 0;
128 ellI2 = 0;
129 ellI3 = 0;
130 node = 0;
131}
I2_0
#define I2_0
Definition strains_pomp.c:17
S_0
#define S_0
Definition strains_pomp.c:15
I3_0
#define I3_0
Definition strains_pomp.c:18
I1_0
#define I1_0
Definition strains_pomp.c:16
R_0
#define R_0
Definition strains_pomp.c:19

◆ 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}
psi1
#define psi1
Definition twospecies_pomp.c:47
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
psi2
#define psi2
Definition twospecies_pomp.c:48
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