phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
seirs_pomp.c
Go to the documentation of this file.
1 #include "pomplink.h"
2 #include "internal.h"
3 
4 static const int nrate = 7;
5 
6 static inline int random_choice (double n) {
7  return floor(R_unif_index(n));
8 }
9 
10 static void change_color (double *color, int nsample,
11  int n, int from, int to) {
12  int i = -1;
13  while (n >= 0 && i < nsample) {
14  i++;
15  if (!ISNA(color[i]) && nearbyint(color[i]) == from) n--;
16  }
17  assert(i < nsample);
18  assert(n == -1);
19  assert(nearbyint(color[i]) == from);
20  color[i] = to;
21 }
22 
23 #define Beta (__p[__parindex[0]])
24 #define sigma (__p[__parindex[1]])
25 #define gamma (__p[__parindex[2]])
26 #define psi (__p[__parindex[3]])
27 #define omega (__p[__parindex[4]])
28 #define S0 (__p[__parindex[5]])
29 #define E0 (__p[__parindex[6]])
30 #define I0 (__p[__parindex[7]])
31 #define R0 (__p[__parindex[8]])
32 #define N (__p[__parindex[9]])
33 #define S (__x[__stateindex[0]])
34 #define E (__x[__stateindex[1]])
35 #define I (__x[__stateindex[2]])
36 #define R (__x[__stateindex[3]])
37 #define ll (__x[__stateindex[4]])
38 #define node (__x[__stateindex[5]])
39 #define ellE (__x[__stateindex[6]])
40 #define ellI (__x[__stateindex[7]])
41 #define COLOR (__x[__stateindex[8]])
42 
43 #define EVENT_RATES \
44  event_rates(__x,__p,t, \
45  __stateindex,__parindex,__covindex, \
46  __covars,rate,logpi,&penalty) \
47 
48 static double event_rates
49 (
50  double *__x,
51  const double *__p,
52  double t,
53  const int *__stateindex,
54  const int *__parindex,
55  const int *__covindex,
56  const double *__covars,
57  double *rate,
58  double *logpi,
59  double *penalty
60  ) {
61  double event_rate = 0;
62  double alpha, pi;
63  *penalty = 0;
64  // 0: transmission, s=(0,0)
65  assert(S>=0 && I>=0);
66  alpha = (N > 0) ? Beta*S*I/N : 0;
67  pi = (I > 0) ? 1-ellI/I : 0;
68  assert(I >= ellI);
69  event_rate += (*rate = alpha*pi); rate++;
70  *logpi = log(pi); logpi++;
71  // 1: transmission, s=(0,1)
72  pi = (1-pi)/2;
73  event_rate += (*rate = alpha*pi); rate++;
74  *logpi = log(pi)-log(ellI); logpi++;
75  // 2: transmission, s=(1,0)
76  event_rate += (*rate = alpha*pi); rate++;
77  *logpi = log(pi)-log(ellI); logpi++;
78  // 3: progression, s=(0,0)
79  assert(E>=0);
80  alpha = sigma*E;
81  pi = (E > 0) ? ellE/E : 1;
82  assert(E >= ellE);
83  event_rate += (*rate = alpha*(1-pi)); rate++;
84  *logpi = log(1-pi); logpi++;
85  // 4: progression, s=(0,1)
86  event_rate += (*rate = alpha*pi); rate++;
87  *logpi = log(pi)-log(ellE); logpi++;
88  // 5: recovery
89  assert(I>=0);
90  alpha = gamma*I;
91  if (I > ellI) {
92  event_rate += (*rate = alpha); rate++;
93  *logpi = 0; logpi++;
94  } else {
95  *rate = 0; rate++;
96  *logpi = 0; logpi++;
97  *penalty += alpha;
98  }
99  // 6: waning
100  event_rate += (*rate = omega*R); rate++;
101  *logpi = 0; logpi++;
102  // 7: sampling (Q = 0)
103  *penalty += psi*I;
104  assert(R_FINITE(event_rate));
105  return event_rate;
106 }
107 
115 (
116  double *__x,
117  const double *__p,
118  double t0,
119  const int *__stateindex,
120  const int *__parindex,
121  const int *__covindex,
122  const double *__covars
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 }
134 
140 (
141  double *__x,
142  const double *__p,
143  const int *__stateindex,
144  const int *__parindex,
145  const int *__covindex,
146  const double *__covars,
147  double t,
148  double dt
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 }
340 
341 # define lik (__lik[0])
342 
345 (
346  double *__lik,
347  const double *__y,
348  const double *__x,
349  const double *__p,
350  int give_log,
351  const int *__obsindex,
352  const int *__stateindex,
353  const int *__parindex,
354  const int *__covindex,
355  const double *__covars,
356  double t
357  ) {
358  assert(!ISNAN(ll));
359  lik = (give_log) ? ll : exp(ll);
360 }
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 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 R0
Definition: seirs_pomp.c:31
#define gamma
Definition: seirs_pomp.c:25
#define COLOR
Definition: seirs_pomp.c:41
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)
Definition: seirs_pomp.c:49
static int random_choice(double n)
Definition: seirs_pomp.c:6
#define R
Definition: seirs_pomp.c:36
#define S0
Definition: seirs_pomp.c:28
#define I
Definition: seirs_pomp.c:35
#define sigma
Definition: seirs_pomp.c:24
#define Beta
Definition: seirs_pomp.c:23
#define I0
Definition: seirs_pomp.c:30
#define E0
Definition: seirs_pomp.c:29
#define ellI
Definition: seirs_pomp.c:40
void seirs_gill(double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)
Definition: seirs_pomp.c:140
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: seirs_pomp.c:345
static const int nrate
Definition: seirs_pomp.c:4
#define lik
Definition: seirs_pomp.c:341
void seirs_rinit(double *__x, const double *__p, double t0, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars)
Definition: seirs_pomp.c:115
#define ll
Definition: seirs_pomp.c:37
#define omega
Definition: seirs_pomp.c:27
#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