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 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

◆ 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}
lik
#define lik
Definition lbdp_pomp.c:165
ll
#define ll
Definition lbdp_pomp.c:10

◆ lbdp_gill()

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

Latent-state process simulator (rprocess).

This integrates the filter equation.

Definition at line 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}
get_userdata_int
get_userdata_int_t * get_userdata_int
Definition init.c:7
rcateg
static int rcateg(double erate, double *rate, int nrate)
Definition internal.h:85
chi
#define chi
Definition lbdp_pomp.c:7
ell
#define ell
Definition lbdp_pomp.c:11
n
#define n
Definition lbdp_pomp.c:9
lambda
#define lambda
Definition lbdp_pomp.c:4
psi
#define psi
Definition lbdp_pomp.c:6
EVENT_RATES
#define EVENT_RATES
Definition lbdp_pomp.c:14
node
#define node
Definition lbdp_pomp.c:12
Here is the call graph for this function:

◆ lbdp_rinit()

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

Latent-state initializer (rinit).

Definition at line 53 of file lbdp_pomp.c.

62 {
63 n = nearbyint(n0);
64 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 64 of file init.c.

64 {
65 // Register routines
66 R_registerRoutines(info,NULL,callMethods,NULL,extMethods);
67 R_useDynamicSymbols(info,TRUE);
68 // R_useDynamicSymbols(info,FALSE);
69 // R_forceSymbols(info,TRUE);
70 get_userdata = (get_userdata_t*) R_GetCCallable("pomp","get_userdata");
71 get_userdata_double = (get_userdata_double_t*) R_GetCCallable("pomp","get_userdata_double");
72 get_userdata_int = (get_userdata_int_t*) R_GetCCallable("pomp","get_userdata_int");
73}
extMethods
static const R_CallMethodDef extMethods[]
Definition init.c:58
callMethods
static const R_CallMethodDef callMethods[]
Definition init.c:35
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 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 (unif_rand() < psi/(psi+chi)) { // non-destructive sample
227 ll += log(I-ellI);
228 } else { // destructive sample
229 ll += log(I);
230 I -= 1;
231 }
232 } else {
233 assert(0); // #nocov
234 ll += R_NegInf; // #nocov
235 }
236 color[parlin] = R_NaReal;
237 break;
238 case 2: // branch point s=(1,1)
239 // If parent is not in deme I, likelihood = 0.
240 if (parcol != Infected) {
241 ll += R_NegInf;
242 color[parlin] = Infected;
243 // the following keeps the state valid
244 ellE -= 1; ellI += 1;
245 E -= 1; I += 1;
246 }
247 assert(sat[parent]==2);
248 ll += (S > 0 && I > 0) ? log(Beta*S/POP/(E+1)) : R_NegInf;
249 S -= 1; E += 1;
250 ellE += 1;
251 S = (S > 0) ? S : 0;
252 int c1 = child[index[parent]];
253 int c2 = child[index[parent]+1];
254 assert(c1 != c2);
255 assert(lineage[c1] != lineage[c2]);
256 assert(lineage[c1] != parlin || lineage[c2] != parlin);
257 assert(lineage[c1] == parlin || lineage[c2] == parlin);
258 if (unif_rand() < 0.5) {
259 color[lineage[c1]] = Exposed;
260 color[lineage[c2]] = Infected;
261 } else {
262 color[lineage[c1]] = Infected;
263 color[lineage[c2]] = Exposed;
264 }
265 ll -= log(0.5);
266 break;
267 }
268
269 // continuous portion of filter equation:
270 // take Gillespie steps to the end of the interval
271 if (tmax > t && R_FINITE(ll)) {
272
273 double rate[nrate], logpi[nrate];
274 int event;
275 double event_rate = 0;
276 double penalty = 0;
277
278 event_rate = EVENT_RATES;
279 tstep = exp_rand()/event_rate;
280
281 while (t + tstep < tmax) {
282 event = rcateg(event_rate,rate,nrate);
283 assert(event>=0 && event<nrate);
284 ll -= penalty*tstep + logpi[event];
285 switch (event) {
286 case 0: // transmission, s=(0,0) or s=(0,1)
287 assert(S>=1);
288 S -= 1; E += 1;
289 ll += log(1-ellE/E);
290 assert(!ISNAN(ll));
291 break;
292 case 1: // transmission, s=(1,0)
293 assert(S>=1);
294 change_color(color,nsample,random_choice(ellI),Infected,Exposed);
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 2: // 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 3: // progression, s=(0,1)
307 assert(E>=1);
308 change_color(color,nsample,random_choice(ellE),Exposed,Infected);
309 ellE -= 1; ellI += 1;
310 E -= 1; I += 1;
311 ll -= log(I);
312 assert(!ISNAN(ll));
313 break;
314 case 4: // recovery
315 assert(I>=1);
316 I -= 1; R += 1;
317 assert(!ISNAN(ll));
318 break;
319 case 5: // 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
ellE
#define ellE
Definition seirs_pomp.c:44
E
#define E
Definition seirs_pomp.c:39
POP
#define POP
Definition seirs_pomp.c:37
change_color
static void change_color(double *color, int nsample, int n, int from, int to)
Definition seirs_pomp.c:13
COLOR
#define COLOR
Definition seirs_pomp.c:46
random_choice
static int random_choice(double n)
Definition seirs_pomp.c:9
R
#define R
Definition seirs_pomp.c:41
I
#define I
Definition seirs_pomp.c:40
Beta
#define Beta
Definition seirs_pomp.c:27
ellI
#define ellI
Definition seirs_pomp.c:45
nrate
static const int nrate
Definition seirs_pomp.c:7
Infected
#define Infected
Definition seirs_pomp.c:5
Exposed
#define Exposed
Definition seirs_pomp.c:4
S
#define S
Definition seirs_pomp.c:38
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
S0
#define S0
Definition seirs_pomp.c:33
I0
#define I0
Definition seirs_pomp.c:35
E0
#define E0
Definition seirs_pomp.c:34

◆ 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 400 of file si2r_pomp.c.

413 {
414 assert(!ISNAN(ll));
415 lik = (give_log) ? ll : exp(ll);
416}

◆ 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 181 of file si2r_pomp.c.

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

165 {
166 double adj = POP/(S0+IL0+IH0+R0);
167 S = nearbyint(S0*adj);
168 IL = nearbyint(IL0*adj);
169 IH = nearbyint(IH0*adj);
170 R = nearbyint(R0*adj);
171 ellL = 0;
172 ellH = 0;
173 ll = 0;
174 node = 0;
175}
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 189 of file sirs_pomp.c.

202 {
203 assert(!ISNAN(ll));
204 lik = (give_log) ? ll : exp(ll);
205}

◆ sirs_gill()

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

Latent-state process simulator (rprocess).

This integrates the filter equation.

Definition at line 90 of file sirs_pomp.c.

100 {
101 double tstep = 0.0, tmax = t + dt;
102 const int *nodetype = get_userdata_int("nodetype");
103 const int *sat = get_userdata_int("saturation");
104
105 int parent = (int) nearbyint(node);
106
107#ifndef NDEBUG
108 int nnode = *get_userdata_int("nnode");
109 assert(parent>=0);
110 assert(parent<=nnode);
111#endif
112
113 ll = 0;
114
115 // singular portion of filter equation
116 switch (nodetype[parent]) {
117 default: // non-genealogical event #nocov
118 break; // #nocov
119 case 0: // root
120 ellI += 1;
121 break;
122 case 1: // sample
123 assert(I >= ellI);
124 assert(ellI >= 0);
125 if (sat[parent] == 1) {
126 ll += log(psi);
127 } else if (sat[parent] == 0) {
128 ellI -= 1;
129 ll += log(psi*(I-ellI));
130 } else {
131 assert(0); // #nocov
132 ll += R_NegInf; // #nocov
133 }
134 break;
135 case 2: // branch point s=(1,1)
136 assert(S >= 0);
137 assert(I >= 0);
138 assert(ellI > 0);
139 assert(sat[parent]==2);
140 ll += (I > 0 && I >= ellI) ? log(Beta*S*I/POP) : 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 = POP/(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 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
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
psi2
#define psi2
Definition twospecies_pomp.c:48
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