phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
sirs_pomp.c
Go to the documentation of this file.
1#include "pomplink.h"
2#include "internal.h"
3
4static const int nrate = 3;
5
6#define Beta (__p[__parindex[0]])
7#define gamma (__p[__parindex[1]])
8#define psi (__p[__parindex[2]])
9#define omega (__p[__parindex[3]])
10#define S0 (__p[__parindex[4]])
11#define I0 (__p[__parindex[5]])
12#define R0 (__p[__parindex[6]])
13#define N (__p[__parindex[7]])
14#define S (__x[__stateindex[0]])
15#define I (__x[__stateindex[1]])
16#define R (__x[__stateindex[2]])
17#define ll (__x[__stateindex[3]])
18#define ellI (__x[__stateindex[4]])
19#define node (__x[__stateindex[5]])
20
21#define EVENT_RATES \
22 event_rates(__x,__p,t, \
23 __stateindex,__parindex,__covindex, \
24 __covars,rate,&penalty) \
25
26static double event_rates
27(
28 double *__x,
29 const double *__p,
30 double t,
31 const int *__stateindex,
32 const int *__parindex,
33 const int *__covindex,
34 const double *__covars,
35 double *rate,
36 double *penalty
37 ) {
38 double event_rate = 0;
39 double alpha, disc;
40 *penalty = 0;
41 assert(I >= ellI);
42 assert(ellI >= 0);
43 assert(S >= 0);
44 // 0: transmission with saturation 0 or 1
45 alpha = Beta*S*I/N;
46 disc = (I > 0) ? ellI*(ellI-1)/I/(I+1) : 1;
47 event_rate += (*rate = alpha*(1-disc)); rate++;
48 *penalty += alpha*disc;
49 // 1: recovery
50 alpha = gamma*I;
51 if (I > ellI) {
52 event_rate += (*rate = alpha); rate++;
53 } else {
54 *rate = 0; rate++;
55 *penalty += alpha;
56 }
57 // 2: loss of immunity
58 alpha = omega*R;
59 event_rate += (*rate = alpha); rate++;
60 // sampling
61 alpha = psi*I;
62 *penalty += alpha;
63 assert(R_FINITE(event_rate));
64 return event_rate;
65}
66
69(
70 double *__x,
71 const double *__p,
72 double t,
73 const int *__stateindex,
74 const int *__parindex,
75 const int *__covindex,
76 const double *__covars
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}
86
91(
92 double *__x,
93 const double *__p,
94 const int *__stateindex,
95 const int *__parindex,
96 const int *__covindex,
97 const double *__covars,
98 double t,
99 double dt
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}
185
186# define lik (__lik[0])
187
190(
191 double *__lik,
192 const double *__y,
193 const double *__x,
194 const double *__p,
195 int give_log,
196 const int *__obsindex,
197 const int *__stateindex,
198 const int *__parindex,
199 const int *__covindex,
200 const double *__covars,
201 double t
202 ){
203 assert(!ISNAN(ll));
204 lik = (give_log) ? ll : exp(ll);
205}
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 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 R0
Definition seirs_pomp.c:31
#define gamma
Definition seirs_pomp.c:25
#define R
Definition seirs_pomp.c:36
#define S0
Definition seirs_pomp.c:28
#define I
Definition seirs_pomp.c:35
#define Beta
Definition seirs_pomp.c:23
#define I0
Definition seirs_pomp.c:30
#define ellI
Definition seirs_pomp.c:40
static const int nrate
Definition seirs_pomp.c:4
#define omega
Definition seirs_pomp.c:27
#define S
Definition seirs_pomp.c:33
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 sirs_pomp.c:190
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 *penalty)
Definition sirs_pomp.c:27
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 sirs_pomp.c:69
void sirs_gill(double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)
Definition sirs_pomp.c:91