phylopomp
Phylodynamics for POMPs
Loading...
Searching...
No Matches
Macros | Functions | Variables
strains_pomp.c File Reference
#include "pomplink.h"
#include "internal.h"
Include dependency graph for strains_pomp.c:

Go to the source code of this file.

Macros

#define STRAIN1   1
 
#define STRAIN2   2
 
#define STRAIN3   3
 
#define Beta1   (__p[__parindex[0]])
 
#define Beta2   (__p[__parindex[1]])
 
#define Beta3   (__p[__parindex[2]])
 
#define gamma   (__p[__parindex[3]])
 
#define chi   (__p[__parindex[4]])
 
#define S_0   (__p[__parindex[5]])
 
#define I1_0   (__p[__parindex[6]])
 
#define I2_0   (__p[__parindex[7]])
 
#define I3_0   (__p[__parindex[8]])
 
#define R_0   (__p[__parindex[9]])
 
#define POP   (__p[__parindex[10]])
 
#define S   (__x[__stateindex[0]])
 
#define I1   (__x[__stateindex[1]])
 
#define I2   (__x[__stateindex[2]])
 
#define I3   (__x[__stateindex[3]])
 
#define R   (__x[__stateindex[4]])
 
#define ll   (__x[__stateindex[5]])
 
#define ellI1   (__x[__stateindex[6]])
 
#define ellI2   (__x[__stateindex[7]])
 
#define ellI3   (__x[__stateindex[8]])
 
#define node   (__x[__stateindex[9]])
 
#define EVENT_RATES
 
#define lik   (__lik[0])
 

Functions

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)
 
void strains_rinit (double *__x, const double *__p, double t, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars)
 Latent-state initializer (rinit).
 
void strains_gill (double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)
 
void strains_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).
 

Variables

static const int nrate = 6
 

Macro Definition Documentation

◆ Beta1

#define Beta1   (__p[__parindex[0]])

Definition at line 10 of file strains_pomp.c.

◆ Beta2

#define Beta2   (__p[__parindex[1]])

Definition at line 11 of file strains_pomp.c.

◆ Beta3

#define Beta3   (__p[__parindex[2]])

Definition at line 12 of file strains_pomp.c.

◆ chi

#define chi   (__p[__parindex[4]])

Definition at line 14 of file strains_pomp.c.

◆ ellI1

#define ellI1   (__x[__stateindex[6]])

Definition at line 27 of file strains_pomp.c.

◆ ellI2

#define ellI2   (__x[__stateindex[7]])

Definition at line 28 of file strains_pomp.c.

◆ ellI3

#define ellI3   (__x[__stateindex[8]])

Definition at line 29 of file strains_pomp.c.

◆ EVENT_RATES

#define EVENT_RATES
Value:
event_rates(__x,__p,t, \
__stateindex,__parindex,__covindex, \
__covars,rate,logpi,&penalty) \
event_rates
static double event_rates(double *__x, const double *__p, double t, const int *__stateindex, const int *__parindex, double *rate, double *penalty)
Definition lbdp_pomp.c:19

Definition at line 32 of file strains_pomp.c.

32#define EVENT_RATES \
33 event_rates(__x,__p,t, \
34 __stateindex,__parindex,__covindex, \
35 __covars,rate,logpi,&penalty) \
36

◆ gamma

#define gamma   (__p[__parindex[3]])

Definition at line 13 of file strains_pomp.c.

◆ I1

#define I1   (__x[__stateindex[1]])

Definition at line 22 of file strains_pomp.c.

◆ I1_0

#define I1_0   (__p[__parindex[6]])

Definition at line 16 of file strains_pomp.c.

◆ I2

#define I2   (__x[__stateindex[2]])

Definition at line 23 of file strains_pomp.c.

◆ I2_0

#define I2_0   (__p[__parindex[7]])

Definition at line 17 of file strains_pomp.c.

◆ I3

#define I3   (__x[__stateindex[3]])

Definition at line 24 of file strains_pomp.c.

◆ I3_0

#define I3_0   (__p[__parindex[8]])

Definition at line 18 of file strains_pomp.c.

◆ lik

#define lik   (__lik[0])

Definition at line 287 of file strains_pomp.c.

◆ ll

#define ll   (__x[__stateindex[5]])

Definition at line 26 of file strains_pomp.c.

◆ node

#define node   (__x[__stateindex[9]])

Definition at line 30 of file strains_pomp.c.

◆ POP

#define POP   (__p[__parindex[10]])

Definition at line 20 of file strains_pomp.c.

◆ R

#define R   (__x[__stateindex[4]])

Definition at line 25 of file strains_pomp.c.

◆ R_0

#define R_0   (__p[__parindex[9]])

Definition at line 19 of file strains_pomp.c.

◆ S

#define S   (__x[__stateindex[0]])

Definition at line 21 of file strains_pomp.c.

◆ S_0

#define S_0   (__p[__parindex[5]])

Definition at line 15 of file strains_pomp.c.

◆ STRAIN1

#define STRAIN1   1

Definition at line 6 of file strains_pomp.c.

◆ STRAIN2

#define STRAIN2   2

Definition at line 7 of file strains_pomp.c.

◆ STRAIN3

#define STRAIN3   3

Definition at line 8 of file strains_pomp.c.

Function Documentation

◆ event_rates()

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

Definition at line 37 of file strains_pomp.c.

49 {
50 double event_rate = 0;
51 double alpha, pi;
52 *penalty = 0;
53 assert(S >= 0);
54 assert(I1 >= ellI1);
55 assert(I2 >= ellI2);
56 assert(I3 >= ellI3);
57 assert(ellI1 >= 0);
58 assert(ellI2 >= 0);
59 assert(ellI3 >= 0);
60 // 0: strain-1 transmission with saturation 0 or 1
61 alpha = Beta1*S*I1/POP;
62 pi = (I1 > 0) ? 1-ellI1*(ellI1-1)/I1/(I1+1) : 0;
63 event_rate += (*rate = alpha*pi); rate++;
64 *logpi = 0; logpi++; // this cancels the boost exactly
65 *penalty += alpha*(1-pi);
66 // 1: strain 1 recovery
67 alpha = gamma*I1;
68 pi = (I1 > ellI1) ? 1 : 0;
69 event_rate += (*rate = alpha*pi); rate++;
70 *logpi = log(pi); logpi++;
71 *penalty += alpha*(1-pi);
72 // strain 1 sampling
73 alpha = chi*I1;
74 *penalty += alpha;
75 // 2: strain-2 transmission with saturation 0 or 1
76 alpha = Beta2*S*I2/POP;
77 pi = (I2 > 0) ? 1-ellI2*(ellI2-1)/I2/(I2+1) : 0;
78 event_rate += (*rate = alpha*pi); rate++;
79 *logpi = 0; logpi++; // this cancels the boost exactly
80 *penalty += alpha*(1-pi);
81 // 3: strain 2 recovery
82 alpha = gamma*I2;
83 pi = (I2 > ellI2) ? 1 : 0;
84 event_rate += (*rate = alpha*pi); rate++;
85 *logpi = log(pi); logpi++;
86 *penalty += alpha*(1-pi);
87 // strain 2 sampling
88 alpha = chi*I2;
89 *penalty += alpha;
90 // 4: strain-3 transmission with saturation 0 or 1
91 alpha = Beta3*S*I3/POP;
92 pi = (I3 > 0) ? 1-ellI3*(ellI3-1)/I3/(I3+1) : 0;
93 event_rate += (*rate = alpha*pi); rate++;
94 *logpi = 0; logpi++; // this cancels the boost exactly
95 *penalty += alpha*(1-pi);
96 // 5: strain 3 recovery
97 alpha = gamma*I3;
98 pi = (I3 > ellI3) ? 1 : 0;
99 event_rate += (*rate = alpha*pi); rate++;
100 *logpi = log(pi); logpi++;
101 *penalty += alpha*(1-pi);
102 // strain 3 sampling
103 alpha = chi*I3;
104 *penalty += alpha;
105 assert(R_FINITE(event_rate));
106 return event_rate;
107}
chi
#define chi
Definition lbdp_pomp.c:7
POP
#define POP
Definition seirs_pomp.c:37
gamma
#define gamma
Definition seirs_pomp.c:29
S
#define S
Definition seirs_pomp.c:38
Beta3
#define Beta3
Definition strains_pomp.c:12
ellI1
#define ellI1
Definition strains_pomp.c:27
ellI2
#define ellI2
Definition strains_pomp.c:28
ellI3
#define ellI3
Definition strains_pomp.c:29
I2
#define I2
Definition strains_pomp.c:23
Beta2
#define Beta2
Definition strains_pomp.c:11
Beta1
#define Beta1
Definition strains_pomp.c:10
I3
#define I3
Definition strains_pomp.c:24
I1
#define I1
Definition strains_pomp.c:22

◆ strains_dmeas()

void strains_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 at line 290 of file strains_pomp.c.

303 {
304 assert(!ISNAN(ll));
305 lik = (give_log) ? ll : exp(ll);
306}
lik
#define lik
Definition lbdp_pomp.c:165
ll
#define ll
Definition lbdp_pomp.c:10

◆ strains_gill()

void strains_gill ( double * __x,
const double * __p,
const int * __stateindex,
const int * __parindex,
const int * __covindex,
const double * __covars,
double t,
double dt )

Latent-state process simulator (rprocess).

This integrates the filter equation.

Definition at line 136 of file strains_pomp.c.

146 {
147 double tstep = 0.0, tmax = t + dt;
148 const int *nodetype = get_userdata_int("nodetype");
149 const int *sat = get_userdata_int("saturation");
150 const int *deme = get_userdata_int("deme");
151 const int *index = get_userdata_int("index");
152 const int *child = get_userdata_int("child");
153
154 int parent = (int) nearbyint(node);
155 int c = child[index[parent]];
156
157#ifndef NDEBUG
158 const int *lineage = get_userdata_int("lineage");
159 int nnode = *get_userdata_int("nnode");
160 assert(parent>=0);
161 assert(parent<=nnode);
162#endif
163
164 ll = 0;
165
166 // singular portion of filter equation
167 switch (nodetype[parent]) {
168 default: // non-genealogical event #nocov
169 break; // #nocov
170 case 0: // root
171 assert(sat[parent]==1);
172 assert(lineage[parent]==lineage[c]);
173 switch (deme[c]) {
174 case STRAIN1:
175 ellI1 += 1; break;
176 case STRAIN2:
177 ellI2 += 1; break;
178 case STRAIN3:
179 ellI3 += 1; break;
180 default: // #nocov
181 assert(0); break; // #nocov
182 }
183 break;
184 case 1: // sample
185 assert(sat[parent]==0);
186 switch (deme[parent]) {
187 case STRAIN1:
188 assert(I1 >= ellI1);
189 assert(ellI1 >= 0);
190 ll += log(chi*I1);
191 ellI1 -= 1; I1 -= 1;
192 break;
193 case STRAIN2:
194 assert(I2 >= ellI2);
195 assert(ellI2 >= 0);
196 ll += log(chi*I2);
197 ellI2 -= 1; I2 -= 1;
198 break;
199 case STRAIN3:
200 assert(I3 >= ellI3);
201 assert(ellI3 >= 0);
202 ll += log(chi*I3);
203 ellI3 -= 1; I3 -= 1;
204 break;
205 default: // #nocov
206 assert(0); break; // #nocov
207 }
208 break;
209 case 2: // branch point s=(1,1)
210 assert(S >= 0);
211 if (sat[parent]!=2) break;
212 switch (deme[parent]) {
213 case STRAIN1:
214 assert(I1 >= 0);
215 assert(ellI1 > 0);
216 ll += (I1 > 0 && I1 >= ellI1) ? log(Beta1*S*I1/POP) : R_NegInf;
217 S -= 1; I1 += 1; ellI1 += 1;
218 ll -= log(I1*(I1-1)/2);
219 break;
220 case STRAIN2:
221 assert(I2 >= 0);
222 assert(ellI2 > 0);
223 ll += (I2 > 0 && I2 >= ellI2) ? log(Beta2*S*I2/POP) : R_NegInf;
224 S -= 1; I2 += 1; ellI2 += 1;
225 ll -= log(I2*(I2-1)/2);
226 break;
227 case STRAIN3:
228 assert(I3 >= 0);
229 assert(ellI3 > 0);
230 ll += (I3 > 0 && I3 >= ellI3) ? log(Beta3*S*I3/POP) : R_NegInf;
231 S -= 1; I3 += 1; ellI3 += 1;
232 ll -= log(I3*(I3-1)/2);
233 break;
234 default: // #nocov
235 assert(0); break; // #nocov
236 }
237 S = (S > 0) ? S : 0;
238 break;
239 }
240
241 if (tmax > t && R_FINITE(ll)) {
242
243 // take Gillespie steps to the end of the interval:
244 int event;
245 double penalty = 0;
246 double rate[nrate], logpi[nrate];
247
248 double event_rate = EVENT_RATES;
249 tstep = exp_rand()/event_rate;
250
251 while (t + tstep < tmax) {
252 event = rcateg(event_rate,rate,nrate);
253 assert(event>=0 && event<nrate);
254 ll -= penalty*tstep + logpi[event];
255 switch (event) {
256 case 0: // strain-1 transmission
257 S -= 1; I1 += 1;
258 break;
259 case 1: // strain-1 recovery
260 I1 -= 1; R += 1;
261 break;
262 case 2: // strain-2 transmission
263 S -= 1; I2 += 1;
264 break;
265 case 3: // strain-2 recovery
266 I2 -= 1; R += 1;
267 break;
268 case 4: // strain-3 transmission
269 S -= 1; I3 += 1;
270 break;
271 case 5: // strain-3 recovery
272 I3 -= 1; R += 1;
273 break;
274 default: // #nocov
275 assert(0); break; // #nocov
276 }
277 t += tstep;
278 event_rate = EVENT_RATES;
279 tstep = exp_rand()/event_rate;
280 }
281 tstep = tmax - t;
282 ll -= penalty*tstep;
283 }
284 node += 1;
285}
get_userdata_int
get_userdata_int_t * get_userdata_int
Definition init.c:7
rcateg
static int rcateg(double erate, double *rate, int nrate)
Definition internal.h:85
deme
static const int deme
Definition lbdp.cc:7
EVENT_RATES
#define EVENT_RATES
Definition lbdp_pomp.c:14
node
#define node
Definition lbdp_pomp.c:12
R
#define R
Definition seirs_pomp.c:41
nrate
static const int nrate
Definition seirs_pomp.c:7
STRAIN3
#define STRAIN3
Definition strains_pomp.c:8
STRAIN1
#define STRAIN1
Definition strains_pomp.c:6
STRAIN2
#define STRAIN2
Definition strains_pomp.c:7
Here is the call graph for this function:

◆ strains_rinit()

void strains_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 at line 110 of file strains_pomp.c.

119 {
120 double m = POP/(S_0+I1_0+I2_0+I3_0+R_0);
121 S = nearbyint(S_0*m);
122 I1 = nearbyint(I1_0*m);
123 I2 = nearbyint(I2_0*m);
124 I3 = nearbyint(I3_0*m);
125 R = nearbyint(R_0*m);
126 ll = 0;
127 ellI1 = 0;
128 ellI2 = 0;
129 ellI3 = 0;
130 node = 0;
131}
I2_0
#define I2_0
Definition strains_pomp.c:17
S_0
#define S_0
Definition strains_pomp.c:15
I3_0
#define I3_0
Definition strains_pomp.c:18
I1_0
#define I1_0
Definition strains_pomp.c:16
R_0
#define R_0
Definition strains_pomp.c:19

Variable Documentation

◆ nrate

const int nrate = 6
static

Definition at line 4 of file strains_pomp.c.