phylopomp
Phylodynamics for POMPs
Loading...
Searching...
No Matches
Macros | Functions | Variables
twospecies_pomp.c File Reference
#include "pomplink.h"
#include "internal.h"
Include dependency graph for twospecies_pomp.c:

Go to the source code of this file.

Macros

#define host1   0
 
#define host2   1
 
#define Beta11   (__p[__parindex[0]])
 
#define Beta12   (__p[__parindex[1]])
 
#define Beta21   (__p[__parindex[2]])
 
#define Beta22   (__p[__parindex[3]])
 
#define gamma1   (__p[__parindex[4]])
 
#define gamma2   (__p[__parindex[5]])
 
#define psi1   (__p[__parindex[6]])
 
#define psi2   (__p[__parindex[7]])
 
#define omega1   (__p[__parindex[8]])
 
#define omega2   (__p[__parindex[9]])
 
#define b1   (__p[__parindex[10]])
 
#define b2   (__p[__parindex[11]])
 
#define d1   (__p[__parindex[12]])
 
#define d2   (__p[__parindex[13]])
 
#define C1   (__p[__parindex[14]])
 
#define C2   (__p[__parindex[15]])
 
#define S1_0   (__p[__parindex[16]])
 
#define S2_0   (__p[__parindex[17]])
 
#define I1_0   (__p[__parindex[18]])
 
#define I2_0   (__p[__parindex[19]])
 
#define R1_0   (__p[__parindex[20]])
 
#define R2_0   (__p[__parindex[21]])
 
#define S1   (__x[__stateindex[0]])
 
#define S2   (__x[__stateindex[1]])
 
#define I1   (__x[__stateindex[2]])
 
#define I2   (__x[__stateindex[3]])
 
#define R1   (__x[__stateindex[4]])
 
#define R2   (__x[__stateindex[5]])
 
#define N1   (__x[__stateindex[6]])
 
#define N2   (__x[__stateindex[7]])
 
#define ll   (__x[__stateindex[8]])
 
#define node   (__x[__stateindex[9]])
 
#define ell1   (__x[__stateindex[10]])
 
#define ell2   (__x[__stateindex[11]])
 
#define COLOR   (__x[__stateindex[12]])
 
#define EVENT_RATES
 
#define lik   (__lik[0])
 

Functions

static int random_choice (double n)
 
static void change_color (double *color, int nsample, int n, int from, int to)
 
static int check_color (double *color, int nsample, double size1, double size2)
 
static double event_rates (double *__x, const double *__p, double t, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double *rate, double *logpi, double *penalty)
 
void twospecies_rinit (double *__x, const double *__p, double t0, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars)
 
void twospecies_gill (double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)
 
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)
 Measurement model likelihood (dmeasure).
 

Variables

static const int nrate = 18
 

Macro Definition Documentation

◆ b1

#define b1   (__p[__parindex[10]])

Definition at line 51 of file twospecies_pomp.c.

◆ b2

#define b2   (__p[__parindex[11]])

Definition at line 52 of file twospecies_pomp.c.

◆ Beta11

#define Beta11   (__p[__parindex[0]])

Definition at line 41 of file twospecies_pomp.c.

◆ Beta12

#define Beta12   (__p[__parindex[1]])

Definition at line 42 of file twospecies_pomp.c.

◆ Beta21

#define Beta21   (__p[__parindex[2]])

Definition at line 43 of file twospecies_pomp.c.

◆ Beta22

#define Beta22   (__p[__parindex[3]])

Definition at line 44 of file twospecies_pomp.c.

◆ C1

#define C1   (__p[__parindex[14]])

Definition at line 55 of file twospecies_pomp.c.

◆ C2

#define C2   (__p[__parindex[15]])

Definition at line 56 of file twospecies_pomp.c.

◆ COLOR

#define COLOR   (__x[__stateindex[12]])

Definition at line 75 of file twospecies_pomp.c.

◆ d1

#define d1   (__p[__parindex[12]])

Definition at line 53 of file twospecies_pomp.c.

◆ d2

#define d2   (__p[__parindex[13]])

Definition at line 54 of file twospecies_pomp.c.

◆ ell1

#define ell1   (__x[__stateindex[10]])

Definition at line 73 of file twospecies_pomp.c.

◆ ell2

#define ell2   (__x[__stateindex[11]])

Definition at line 74 of file twospecies_pomp.c.

◆ EVENT_RATES

#define EVENT_RATES
Value:
event_rates(__x,__p,t, \
__stateindex,__parindex,__covindex, \
__covars,rate,logpi,&penalty) \
event_rates
static double event_rates(double *__x, const double *__p, double t, const int *__stateindex, const int *__parindex, double *rate, double *penalty)
Definition lbdp_pomp.c:18

Definition at line 77 of file twospecies_pomp.c.

77#define EVENT_RATES \
78 event_rates(__x,__p,t, \
79 __stateindex,__parindex,__covindex, \
80 __covars,rate,logpi,&penalty) \
81

◆ gamma1

#define gamma1   (__p[__parindex[4]])

Definition at line 45 of file twospecies_pomp.c.

◆ gamma2

#define gamma2   (__p[__parindex[5]])

Definition at line 46 of file twospecies_pomp.c.

◆ host1

#define host1   0

Definition at line 5 of file twospecies_pomp.c.

◆ host2

#define host2   1

Definition at line 6 of file twospecies_pomp.c.

◆ I1

#define I1   (__x[__stateindex[2]])

Definition at line 65 of file twospecies_pomp.c.

◆ I1_0

#define I1_0   (__p[__parindex[18]])

Definition at line 59 of file twospecies_pomp.c.

◆ I2

#define I2   (__x[__stateindex[3]])

Definition at line 66 of file twospecies_pomp.c.

◆ I2_0

#define I2_0   (__p[__parindex[19]])

Definition at line 60 of file twospecies_pomp.c.

◆ lik

#define lik   (__lik[0])

Definition at line 631 of file twospecies_pomp.c.

◆ ll

#define ll   (__x[__stateindex[8]])

Definition at line 71 of file twospecies_pomp.c.

◆ N1

#define N1   (__x[__stateindex[6]])

Definition at line 69 of file twospecies_pomp.c.

◆ N2

#define N2   (__x[__stateindex[7]])

Definition at line 70 of file twospecies_pomp.c.

◆ node

#define node   (__x[__stateindex[9]])

Definition at line 72 of file twospecies_pomp.c.

◆ omega1

#define omega1   (__p[__parindex[8]])

Definition at line 49 of file twospecies_pomp.c.

◆ omega2

#define omega2   (__p[__parindex[9]])

Definition at line 50 of file twospecies_pomp.c.

◆ psi1

#define psi1   (__p[__parindex[6]])

Definition at line 47 of file twospecies_pomp.c.

◆ psi2

#define psi2   (__p[__parindex[7]])

Definition at line 48 of file twospecies_pomp.c.

◆ R1

#define R1   (__x[__stateindex[4]])

Definition at line 67 of file twospecies_pomp.c.

◆ R1_0

#define R1_0   (__p[__parindex[20]])

Definition at line 61 of file twospecies_pomp.c.

◆ R2

#define R2   (__x[__stateindex[5]])

Definition at line 68 of file twospecies_pomp.c.

◆ R2_0

#define R2_0   (__p[__parindex[21]])

Definition at line 62 of file twospecies_pomp.c.

◆ S1

#define S1   (__x[__stateindex[0]])

Definition at line 63 of file twospecies_pomp.c.

◆ S1_0

#define S1_0   (__p[__parindex[16]])

Definition at line 57 of file twospecies_pomp.c.

◆ S2

#define S2   (__x[__stateindex[1]])

Definition at line 64 of file twospecies_pomp.c.

◆ S2_0

#define S2_0   (__p[__parindex[17]])

Definition at line 58 of file twospecies_pomp.c.

Function Documentation

◆ change_color()

static void change_color ( double * color,
int nsample,
int n,
int from,
int to )
static

Definition at line 12 of file twospecies_pomp.c.

13 {
14 int i = -1;
15 while (n >= 0 && i < nsample) {
16 i++;
17 if (!ISNA(color[i]) && nearbyint(color[i]) == from) n--;
18 }
19 assert(i < nsample);
20 assert(n == -1);
21 assert(nearbyint(color[i]) == from);
22 color[i] = to;
23}
nsample
SEXP nsample(TYPE &X)
Definition generics.h:12
n
#define n
Definition lbdp_pomp.c:8
Here is the call graph for this function:
Here is the caller graph for this function:

◆ check_color()

static int check_color ( double * color,
int nsample,
double size1,
double size2 )
static

Definition at line 26 of file twospecies_pomp.c.

27 {
28 int n1 = 0, n2 = 0;
29 int s1 = (int) nearbyint(size1);
30 int s2 = (int) nearbyint(size2);
31 for (int i = 0; i < nsample; i++) {
32 if (!ISNA(color[i])) {
33 if (nearbyint(color[i]) == host1) n1++;
34 if (nearbyint(color[i]) == host2) n2++;
35 }
36 }
37 return (n1==s1) && (n2==s2);
38}
host1
#define host1
Definition twospecies_pomp.c:5
host2
#define host2
Definition twospecies_pomp.c:6
Here is the call graph for this function:
Here is the caller graph for this function:

◆ event_rates()

static double event_rates ( double * __x,
const double * __p,
double t,
const int * __stateindex,
const int * __parindex,
const int * __covindex,
const double * __covars,
double * rate,
double * logpi,
double * penalty )
static

Definition at line 84 of file twospecies_pomp.c.

96 {
97 double event_rate = 0;
98 double alpha, pi, disc;
99 *penalty = 0;
100 // 0: Trans_11, s = (0,0),(1,0)
101 assert(S1>=0 && I1>=ell1 && ell1>=0);
102 alpha = (N1 > 0) ? Beta11*S1*I1/N1 : 0;
103 disc = (I1 > 0) ? ell1*(ell1-1)/I1/(I1+1) : 1;
104 *penalty += alpha*disc;
105 event_rate += (*rate = alpha*(1-disc)); rate++;
106 *logpi = 0; logpi++;
107 assert(R_FINITE(event_rate));
108 // 1: Trans_22, s = (0,0),(0,1)
109 assert(S2>=0 && I2>=ell2 && ell2>=0);
110 alpha = (N2 > 0) ? Beta22*S2*I2/N2 : 0;
111 disc = (I2 > 0) ? ell2*(ell2-1)/I2/(I2+1) : 1;
112 *penalty += alpha*disc;
113 event_rate += (*rate = alpha*(1-disc)); rate++;
114 *logpi = 0; logpi++;
115 assert(R_FINITE(event_rate));
116 // 2: Trans_21, s = (0,0),(1,0)
117 // s = (0,0): pi = (I1-ell1)/I1, m = 1,
118 // s = (1,0): pi = 1/(2*I1), m = ell1
119 // total pi = (I1-0.5*ell1)/I1
120 assert(S2>=0 && I1>=ell1 && ell1>=0);
121 alpha = (N1 > 0) ? Beta21*S2*I1/N1 : 0;
122 pi = (I1 > 0) ? (I1-0.5*ell1)/I1 : 0;
123 event_rate += (*rate = alpha*pi); rate++;
124 *logpi = log(pi); logpi++;
125 assert(R_FINITE(event_rate));
126 // 3: Trans_21, s = (0,1)
127 // s = (0,1): pi = 1/(2*I1), m = ell1
128 pi = 1-pi;
129 event_rate += (*rate = alpha*pi); rate++;
130 *logpi = log(pi)-log(ell1); logpi++;
131 assert(R_FINITE(event_rate));
132 // 4: Trans_12, s = (0,0),(0,1)
133 assert(S1>=0 && I2>=ell2 && ell2>=0);
134 alpha = (N2 > 0) ? Beta12*S1*I2/N2 : 0;
135 pi = (I2 > 0) ? (I2-ell2*0.5)/I2 : 0;
136 event_rate += (*rate = alpha*pi); rate++;
137 *logpi = log(pi); logpi++;
138 assert(R_FINITE(event_rate));
139 // 5: Trans_12, s = (1,0)
140 pi = 1-pi;
141 event_rate += (*rate = alpha*pi); rate++;
142 *logpi = log(pi)-log(ell2); logpi++;
143 assert(R_FINITE(event_rate));
144 // 6: Recov_1
145 assert(I1>=0 && ell1>=0);
146 alpha = gamma1*I1;
147 if (I1 > ell1) {
148 event_rate += (*rate = alpha); rate++;
149 } else {
150 *rate = 0; rate++;
151 *penalty += alpha;
152 }
153 *logpi = 0; logpi++;
154 assert(R_FINITE(event_rate));
155 // 7: Recov_2
156 assert(I2>=0 && ell2>=0);
157 alpha = gamma2*I2;
158 if (I2 > ell2) {
159 event_rate += (*rate = alpha); rate++;
160 } else {
161 *rate = 0; rate++;
162 *penalty += alpha;
163 }
164 *logpi = 0; logpi++;
165 assert(R_FINITE(event_rate));
166 // 8: Wane_1
167 assert(R1>=0);
168 alpha = omega1*R1;
169 event_rate += (*rate = alpha); rate++;
170 *logpi = 0; logpi++;
171 assert(R_FINITE(event_rate));
172 // 9: Wane_2
173 assert(R2>=0);
174 alpha = omega2*R2;
175 event_rate += (*rate = alpha); rate++;
176 *logpi = 0; logpi++;
177 assert(R_FINITE(event_rate));
178 // 10: Birth_1
179 assert(N1>=0);
180 alpha = b1*N1;
181 event_rate += (*rate = alpha); rate++;
182 *logpi = 0; logpi++;
183 assert(R_FINITE(event_rate));
184 // 11: Birth_2
185 assert(N2>=0);
186 alpha = b2*N2;
187 event_rate += (*rate = alpha); rate++;
188 *logpi = 0; logpi++;
189 assert(R_FINITE(event_rate));
190 // 12: Death_S1
191 assert(S1>=0);
192 alpha = d1*S1;
193 event_rate += (*rate = alpha); rate++;
194 *logpi = 0; logpi++;
195 assert(R_FINITE(event_rate));
196 // 13: Death_I1
197 assert(I1>=0);
198 alpha = d1*I1;
199 if (I1 > ell1) {
200 event_rate += (*rate = alpha); rate++;
201 } else {
202 *rate = 0; rate++;
203 *penalty += alpha;
204 }
205 *logpi = 0; logpi++;
206 assert(R_FINITE(event_rate));
207 // 14: Death_R1
208 assert(R1>=0);
209 alpha = d1*R1;
210 event_rate += (*rate = alpha); rate++;
211 *logpi = 0; logpi++;
212 assert(R_FINITE(event_rate));
213 // 15: Death_S2
214 assert(S2>=0);
215 alpha = d2*S2;
216 event_rate += (*rate = alpha); rate++;
217 *logpi = 0; logpi++;
218 assert(R_FINITE(event_rate));
219 // 16: Death_I2
220 assert(I2>=0);
221 alpha = d2*I2;
222 if (I2 > ell2) {
223 event_rate += (*rate = alpha); rate++;
224 } else {
225 *rate = 0; rate++;
226 *penalty += alpha;
227 }
228 *logpi = 0; logpi++;
229 assert(R_FINITE(event_rate));
230 // 17: Death_R2
231 assert(R2>=0);
232 alpha = d2*R2;
233 event_rate += (*rate = alpha); rate++;
234 *logpi = 0; logpi++;
235 assert(R_FINITE(event_rate));
236 // Sample_1 (Q = 0) // NB: (1-C1)*psi1+C1*psi1 = psi1
237 *penalty += psi1*I1;
238 // Sample_2 (Q = 0) // NB: (1-C2)*psi2+C2*psi2 = psi2
239 *penalty += psi2*I2;
240 return event_rate;
241}
psi1
#define psi1
Definition twospecies_pomp.c:47
ell1
#define ell1
Definition twospecies_pomp.c:73
d1
#define d1
Definition twospecies_pomp.c:53
Beta21
#define Beta21
Definition twospecies_pomp.c:43
Beta12
#define Beta12
Definition twospecies_pomp.c:42
b1
#define b1
Definition twospecies_pomp.c:51
Beta22
#define Beta22
Definition twospecies_pomp.c:44
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
Beta11
#define Beta11
Definition twospecies_pomp.c:41
N1
#define N1
Definition twospecies_pomp.c:69
omega2
#define omega2
Definition twospecies_pomp.c:50
gamma1
#define gamma1
Definition twospecies_pomp.c:45
N2
#define N2
Definition twospecies_pomp.c:70
b2
#define b2
Definition twospecies_pomp.c:52
S2
#define S2
Definition twospecies_pomp.c:64
psi2
#define psi2
Definition twospecies_pomp.c:48
d2
#define d2
Definition twospecies_pomp.c:54
I1
#define I1
Definition twospecies_pomp.c:65
ell2
#define ell2
Definition twospecies_pomp.c:74
gamma2
#define gamma2
Definition twospecies_pomp.c:46
omega1
#define omega1
Definition twospecies_pomp.c:49

◆ random_choice()

static int random_choice ( double n)
inlinestatic

Definition at line 8 of file twospecies_pomp.c.

8 {
9 return floor(R_unif_index(n));
10}
Here is the caller graph for this function:

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

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}
lik
#define lik
Definition lbdp_pomp.c:159
ll
#define ll
Definition lbdp_pomp.c:9

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

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}
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
lambda
#define lambda
Definition lbdp_pomp.c:4
EVENT_RATES
#define EVENT_RATES
Definition lbdp_pomp.c:13
node
#define node
Definition lbdp_pomp.c:11
COLOR
#define COLOR
Definition seirs_pomp.c:41
nrate
static const int nrate
Definition seirs_pomp.c:4
change_color
static void change_color(double *color, int nsample, int n, int from, int to)
Definition twospecies_pomp.c:12
check_color
static int check_color(double *color, int nsample, double size1, double size2)
Definition twospecies_pomp.c:26
C1
#define C1
Definition twospecies_pomp.c:55
random_choice
static int random_choice(double n)
Definition twospecies_pomp.c:8
C2
#define C2
Definition twospecies_pomp.c:56
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 )

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

Variable Documentation

◆ nrate

const int nrate = 18
static

Definition at line 4 of file twospecies_pomp.c.