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
4static const int nrate = 7;
5
6static inline int random_choice (double n) {
7 return floor(R_unif_index(n));
8}
9
10static 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
48static 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 lik
Definition lbdp_pomp.c:159
#define n
Definition lbdp_pomp.c:8
#define ll
Definition lbdp_pomp.c:9
#define psi
Definition lbdp_pomp.c:6
#define EVENT_RATES
Definition lbdp_pomp.c:13
#define node
Definition lbdp_pomp.c:11
#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
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 omega
Definition seirs_pomp.c:27
#define S
Definition seirs_pomp.c:33