phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
twospecies_pomp.c
Go to the documentation of this file.
1 #include "pomplink.h"
2 #include "internal.h"
3 
4 static const int nrate = 18;
5 #define host1 0
6 #define host2 1
7 
8 static inline int random_choice (double n) {
9  return floor(R_unif_index(n));
10 }
11 
12 static void change_color (double *color, int nsample,
13  int n, int from, int to) {
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 }
24 
25 #ifndef NDEBUG
26 static int check_color (double *color, int nsample,
27  double size1, double size2) {
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 }
39 #endif
40 
41 #define Beta11 (__p[__parindex[0]])
42 #define Beta12 (__p[__parindex[1]])
43 #define Beta21 (__p[__parindex[2]])
44 #define Beta22 (__p[__parindex[3]])
45 #define gamma1 (__p[__parindex[4]])
46 #define gamma2 (__p[__parindex[5]])
47 #define psi1 (__p[__parindex[6]])
48 #define psi2 (__p[__parindex[7]])
49 #define omega1 (__p[__parindex[8]])
50 #define omega2 (__p[__parindex[9]])
51 #define b1 (__p[__parindex[10]])
52 #define b2 (__p[__parindex[11]])
53 #define d1 (__p[__parindex[12]])
54 #define d2 (__p[__parindex[13]])
55 #define C1 (__p[__parindex[14]])
56 #define C2 (__p[__parindex[15]])
57 #define S1_0 (__p[__parindex[16]])
58 #define S2_0 (__p[__parindex[17]])
59 #define I1_0 (__p[__parindex[18]])
60 #define I2_0 (__p[__parindex[19]])
61 #define R1_0 (__p[__parindex[20]])
62 #define R2_0 (__p[__parindex[21]])
63 #define S1 (__x[__stateindex[0]])
64 #define S2 (__x[__stateindex[1]])
65 #define I1 (__x[__stateindex[2]])
66 #define I2 (__x[__stateindex[3]])
67 #define R1 (__x[__stateindex[4]])
68 #define R2 (__x[__stateindex[5]])
69 #define N1 (__x[__stateindex[6]])
70 #define N2 (__x[__stateindex[7]])
71 #define ll (__x[__stateindex[8]])
72 #define node (__x[__stateindex[9]])
73 #define ell1 (__x[__stateindex[10]])
74 #define ell2 (__x[__stateindex[11]])
75 #define COLOR (__x[__stateindex[12]])
76 
77 #define EVENT_RATES \
78  event_rates(__x,__p,t, \
79  __stateindex,__parindex,__covindex, \
80  __covars,rate,logpi,&penalty) \
81 
82 // FIXME: At the moment, the following codes exclude the possibility of
83 // importation of infection.
84 static double event_rates
85 (
86  double *__x,
87  const double *__p,
88  double t,
89  const int *__stateindex,
90  const int *__parindex,
91  const int *__covindex,
92  const double *__covars,
93  double *rate,
94  double *logpi,
95  double *penalty
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 }
242 
250 (
251  double *__x,
252  const double *__p,
253  double t0,
254  const int *__stateindex,
255  const int *__parindex,
256  const int *__covindex,
257  const double *__covars
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 }
275 
291 (
292  double *__x,
293  const double *__p,
294  const int *__stateindex,
295  const int *__parindex,
296  const int *__covindex,
297  const double *__covars,
298  double t,
299  double dt
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 }
630 
631 # define lik (__lik[0])
632 
635 (
636  double *__lik,
637  const double *__y,
638  const double *__x,
639  const double *__p,
640  int give_log,
641  const int *__obsindex,
642  const int *__stateindex,
643  const int *__parindex,
644  const int *__covindex,
645  const double *__covars,
646  double t
647  ) {
648  assert(!ISNAN(ll));
649  lik = (give_log) ? ll : exp(ll);
650 }
SEXP nsample(TYPE &X)
Definition: generics.h:12
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 n
Definition: lbdp_pomp.c:8
#define lambda
Definition: lbdp_pomp.c:4
#define I2_0
#define psi1
#define ell1
#define S2_0
static void change_color(double *color, int nsample, int n, int from, int to)
#define d1
#define Beta21
static int check_color(double *color, int nsample, double size1, double size2)
#define Beta12
#define S1_0
#define C1
#define COLOR
#define b1
#define Beta22
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 int random_choice(double n)
#define S1
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).
void twospecies_rinit(double *__x, const double *__p, double t0, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars)
#define R2_0
#define R2
#define R1
#define I2
#define host1
static const int nrate
#define lik
#define R1_0
#define Beta11
#define I1_0
#define N1
#define omega2
#define host2
#define gamma1
#define N2
#define b2
#define S2
#define C2
#define psi2
#define ll
#define d2
#define I1
#define ell2
void twospecies_gill(double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)
#define EVENT_RATES
#define node
#define gamma2
#define omega1