phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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). More...
 
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). More...
 
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). More...
 
void sirs_rinit (double *, const double *, double, const int *, const int *, const int *, const double *)
 Latent-state initializer (rinit). More...
 
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). More...
 
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). More...
 

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 
)

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

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_t * get_userdata_int
Definition: init.c:7
static int rcateg(double erate, double *rate, int nrate)
Definition: internal.h:76
#define ell
Definition: lbdp_pomp.c:10
#define n
Definition: lbdp_pomp.c:8
#define lambda
Definition: lbdp_pomp.c:4
#define psi
Definition: lbdp_pomp.c:6
#define EVENT_RATES
Definition: lbdp_pomp.c:13
#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 
)

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 }
#define n0
Definition: lbdp_pomp.c:7

◆ R_init_phylopomp()

void R_init_phylopomp ( DllInfo *  info)

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 }
static const R_CallMethodDef extMethods[]
Definition: init.c:43
static const R_CallMethodDef callMethods[]
Definition: init.c:27
get_userdata_t * get_userdata
Definition: init.c:5
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 
)

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 }
#define lik
Definition: seirs_pomp.c:341
#define ll
Definition: seirs_pomp.c:37

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

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);
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);
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 }
SEXP nsample(TYPE &X)
Definition: generics.h:12
#define N
Definition: seirs_pomp.c:32
#define ellE
Definition: seirs_pomp.c:39
#define E
Definition: seirs_pomp.c:34
static void change_color(double *color, int nsample, int n, int from, int to)
Definition: seirs_pomp.c:10
#define COLOR
Definition: seirs_pomp.c:41
static int random_choice(double n)
Definition: seirs_pomp.c:6
#define R
Definition: seirs_pomp.c:36
#define I
Definition: seirs_pomp.c:35
#define Beta
Definition: seirs_pomp.c:23
#define ellI
Definition: seirs_pomp.c:40
static const int nrate
Definition: seirs_pomp.c:4
#define psi
Definition: seirs_pomp.c:26
#define EVENT_RATES
Definition: seirs_pomp.c:43
#define node
Definition: seirs_pomp.c:38
#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 
)

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 }
#define R0
Definition: seirs_pomp.c:31
#define S0
Definition: seirs_pomp.c:28
#define I0
Definition: seirs_pomp.c:30
#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 
)

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 }
#define lik
Definition: sirs_pomp.c:186
#define ll
Definition: sirs_pomp.c:17

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

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 }
#define N
Definition: sirs_pomp.c:13
#define R
Definition: sirs_pomp.c:16
#define I
Definition: sirs_pomp.c:15
#define Beta
Definition: sirs_pomp.c:6
#define ellI
Definition: sirs_pomp.c:18
static const int nrate
Definition: sirs_pomp.c:4
#define psi
Definition: sirs_pomp.c:8
#define EVENT_RATES
Definition: sirs_pomp.c:21
#define node
Definition: sirs_pomp.c:19
#define S
Definition: sirs_pomp.c:14
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 
)

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 }
#define R0
Definition: sirs_pomp.c:12
#define S0
Definition: sirs_pomp.c:10
#define I0
Definition: sirs_pomp.c:11

◆ 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 }
#define lik
#define ll

◆ 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;
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;
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 }
#define psi1
#define ell1
static void change_color(double *color, int nsample, int n, int from, int to)
#define Beta21
static int check_color(double *color, int nsample, double size1, double size2)
#define Beta12
#define C1
#define COLOR
#define Beta22
static int random_choice(double n)
#define S1
#define R2
#define R1
#define I2
#define host1
static const int nrate
#define Beta11
#define N1
#define host2
#define N2
#define S2
#define C2
#define psi2
#define I1
#define ell2
#define EVENT_RATES
#define node
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 }
#define I2_0
#define S2_0
#define S1_0
#define R2_0
#define R1_0
#define I1_0