4 static const int nrate = 18;
9 return floor(R_unif_index(
n));
13 int n,
int from,
int to) {
17 if (!ISNA(color[i]) && nearbyint(color[i]) == from)
n--;
21 assert(nearbyint(color[i]) == from);
27 double size1,
double size2) {
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++;
37 return (n1==s1) && (n2==s2);
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]])
78 event_rates(__x,__p,t, \
79 __stateindex,__parindex,__covindex, \
80 __covars,rate,logpi,&penalty) \
89 const int *__stateindex,
90 const int *__parindex,
91 const int *__covindex,
92 const double *__covars,
97 double event_rate = 0;
98 double alpha, pi, disc;
104 *penalty += alpha*disc;
105 event_rate += (*rate = alpha*(1-disc)); rate++;
107 assert(R_FINITE(event_rate));
112 *penalty += alpha*disc;
113 event_rate += (*rate = alpha*(1-disc)); rate++;
115 assert(R_FINITE(event_rate));
123 event_rate += (*rate = alpha*pi); rate++;
124 *logpi = log(pi); logpi++;
125 assert(R_FINITE(event_rate));
129 event_rate += (*rate = alpha*pi); rate++;
130 *logpi = log(pi)-log(
ell1); logpi++;
131 assert(R_FINITE(event_rate));
136 event_rate += (*rate = alpha*pi); rate++;
137 *logpi = log(pi); logpi++;
138 assert(R_FINITE(event_rate));
141 event_rate += (*rate = alpha*pi); rate++;
142 *logpi = log(pi)-log(
ell2); logpi++;
143 assert(R_FINITE(event_rate));
148 event_rate += (*rate = alpha); rate++;
154 assert(R_FINITE(event_rate));
159 event_rate += (*rate = alpha); rate++;
165 assert(R_FINITE(event_rate));
169 event_rate += (*rate = alpha); rate++;
171 assert(R_FINITE(event_rate));
175 event_rate += (*rate = alpha); rate++;
177 assert(R_FINITE(event_rate));
181 event_rate += (*rate = alpha); rate++;
183 assert(R_FINITE(event_rate));
187 event_rate += (*rate = alpha); rate++;
189 assert(R_FINITE(event_rate));
193 event_rate += (*rate = alpha); rate++;
195 assert(R_FINITE(event_rate));
200 event_rate += (*rate = alpha); rate++;
206 assert(R_FINITE(event_rate));
210 event_rate += (*rate = alpha); rate++;
212 assert(R_FINITE(event_rate));
216 event_rate += (*rate = alpha); rate++;
218 assert(R_FINITE(event_rate));
223 event_rate += (*rate = alpha); rate++;
229 assert(R_FINITE(event_rate));
233 event_rate += (*rate = alpha); rate++;
235 assert(R_FINITE(event_rate));
254 const int *__stateindex,
255 const int *__parindex,
256 const int *__covindex,
257 const double *__covars
294 const int *__stateindex,
295 const int *__parindex,
296 const int *__covindex,
297 const double *__covars,
301 double tstep = 0.0, tmax = t + dt;
302 double *color = &
COLOR;
310 int parent = (int) nearbyint(
node);
315 assert(parent<=nnode);
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));
326 switch (nodetype[parent]) {
332 assert(sat[parent]==1);
333 int c = child[index[parent]];
334 assert(parlin==lineage[c]);
337 if (unif_rand() < x) {
338 color[lineage[c]] =
host1;
342 color[lineage[c]] =
host2;
349 if (unif_rand() < 0.5) {
350 color[lineage[c]] =
host1;
354 color[lineage[c]] =
host2;
359 assert(nearbyint(
N1)==nearbyint(
S1+
I1+
R1));
360 assert(nearbyint(
N2)==nearbyint(
S2+
I2+
R2));
365 if (sat[parent] == 0) {
366 if (parcol ==
host1) {
368 if (
C1 < 1 && unif_rand() >
C1) {
374 }
else if (parcol ==
host2) {
376 if (
C2 < 1 && unif_rand() >
C2) {
386 }
else if (sat[parent] == 1) {
387 int c = child[index[parent]];
388 color[lineage[c]] = parcol;
391 }
else if (parcol==
host2) {
401 color[parlin] = R_NaReal;
402 assert(nearbyint(
N1)==nearbyint(
S1+
I1+
R1));
403 assert(nearbyint(
N2)==nearbyint(
S2+
I2+
R2));
408 assert(sat[parent]==2);
409 if (parcol ==
host1) {
413 double lambda = lambda11+lambda21;
415 int c1 = child[index[parent]];
416 int c2 = child[index[parent]+1];
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) {
422 color[lineage[c1]] =
host1;
423 color[lineage[c2]] =
host1;
434 if (unif_rand() < 0.5) {
435 color[lineage[c1]] =
host1;
436 color[lineage[c2]] =
host2;
438 color[lineage[c1]] =
host2;
439 color[lineage[c2]] =
host1;
452 }
else if (parcol ==
host2) {
456 double lambda = lambda12+lambda22;
458 int c1 = child[index[parent]];
459 int c2 = child[index[parent]+1];
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) {
465 color[lineage[c1]] =
host2;
466 color[lineage[c2]] =
host2;
477 if (unif_rand() < 0.5) {
478 color[lineage[c1]] =
host1;
479 color[lineage[c2]] =
host2;
481 color[lineage[c1]] =
host2;
482 color[lineage[c2]] =
host1;
499 assert(nearbyint(
N1)==nearbyint(
S1+
I1+
R1));
500 assert(nearbyint(
N2)==nearbyint(
S2+
I2+
R2));
511 double event_rate = 0;
515 tstep = exp_rand()/event_rate;
517 while (t + tstep < tmax) {
519 ll -= penalty*tstep + logpi[event];
522 assert(
S1>=1 &&
I1>=0);
526 assert(
S2>=1 &&
I2>=0);
530 assert(
S2>=1 &&
I1>=0);
536 assert(
S2>=1 &&
I1>=0);
545 assert(
S1>=1 &&
I2>=0);
551 assert(
S1>=1 &&
I2>=0);
584 assert(
S1>=1 &&
N1>=1);
588 assert(
I1>=1 &&
N1>=1);
592 assert(
R1>=1 &&
N1>=1);
596 assert(
S2>=1 &&
N2>=1);
600 assert(
I2>=1 &&
N2>=1);
604 assert(
R2>=1 &&
N2>=1);
616 assert(nearbyint(
N1)==nearbyint(
S1+
I1+
R1));
617 assert(nearbyint(
N2)==nearbyint(
S2+
I2+
R2));
622 tstep = exp_rand()/event_rate;
631 # define lik (__lik[0])
641 const int *__obsindex,
642 const int *__stateindex,
643 const int *__parindex,
644 const int *__covindex,
645 const double *__covars,
649 lik = (give_log) ?
ll : exp(
ll);
get_userdata_int_t * get_userdata_int
static int rcateg(double erate, double *rate, int nrate)
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)
static int random_choice(double n)
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)
void twospecies_gill(double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)