phylopomp
Phylodynamics for POMPs
•All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Functions
decls.h File Reference
This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

void R_init_phylopomp (DllInfo *)
 
void lbdp_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 Latent-state initializer (rinit).
 
void lbdp_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void lbdp_dmeas (double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double)
 Measurement model likelihood (dmeasure).
 
void seirs_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 
void seirs_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void seirs_dmeas (double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double)
 Measurement model likelihood (dmeasure).
 
void sirs_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 Latent-state initializer (rinit).
 
void sirs_gill (double *, const double *, const int *, const int *, const int *, const double *, double, double)
 
void sirs_dmeas (double *, const double *, const double *, const double *, int, const int *, const int *, const int *, const int *, const double *, double)
 Measurement model likelihood (dmeasure).
 
void 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 162 of file lbdp_pomp.c.

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

◆ lbdp_gill()

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

Latent-state process simulator (rprocess).

This integrates the filter equation.

Definition at line 71 of file lbdp_pomp.c.

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

◆ lbdp_rinit()

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

Latent-state initializer (rinit).

Definition at line 52 of file lbdp_pomp.c.

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

◆ R_init_phylopomp()

void R_init_phylopomp ( DllInfo * info)
extern

Definition at line 48 of file init.c.

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

357 {
358 assert(!ISNAN(ll));
359 lik = (give_log) ? ll : exp(ll);
360}

◆ seirs_gill()

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

Simulator for the latent-state process (rprocess).

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

Definition at line 139 of file seirs_pomp.c.

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

◆ seirs_rinit()

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

Latent-state initializer (rinit component).

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

Definition at line 114 of file seirs_pomp.c.

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

◆ sirs_dmeas()

void sirs_dmeas ( double * __lik,
const double * __y,
const double * __x,
const double * __p,
int give_log,
const int * __obsindex,
const int * __stateindex,
const int * __parindex,
const int * __covindex,
const double * __covars,
double t )
extern

Measurement model likelihood (dmeasure).

Definition at line 189 of file sirs_pomp.c.

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

◆ sirs_gill()

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

Latent-state process simulator (rprocess).

This integrates the filter equation.

Definition at line 90 of file sirs_pomp.c.

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

◆ sirs_rinit()

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

Latent-state initializer (rinit).

Definition at line 68 of file sirs_pomp.c.

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

◆ 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 634 of file twospecies_pomp.c.

647 {
648 assert(!ISNAN(ll));
649 lik = (give_log) ? ll : exp(ll);
650}

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