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 psi1   (__p[__parindex[4]])
 
#define psi2   (__p[__parindex[5]])
 
#define psi3   (__p[__parindex[6]])
 
#define S_0   (__p[__parindex[7]])
 
#define I1_0   (__p[__parindex[8]])
 
#define I2_0   (__p[__parindex[9]])
 
#define I3_0   (__p[__parindex[10]])
 
#define R_0   (__p[__parindex[11]])
 
#define N   (__p[__parindex[12]])
 
#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 *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.

◆ ellI1

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

Definition at line 29 of file strains_pomp.c.

◆ ellI2

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

Definition at line 30 of file strains_pomp.c.

◆ ellI3

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

Definition at line 31 of file strains_pomp.c.

◆ EVENT_RATES

#define EVENT_RATES
Value:
event_rates(__x,__p,t, \
__stateindex,__parindex,__covindex, \
__covars,rate,&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 34 of file strains_pomp.c.

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

◆ gamma

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

Definition at line 13 of file strains_pomp.c.

◆ I1

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

Definition at line 24 of file strains_pomp.c.

◆ I1_0

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

Definition at line 18 of file strains_pomp.c.

◆ I2

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

Definition at line 25 of file strains_pomp.c.

◆ I2_0

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

Definition at line 19 of file strains_pomp.c.

◆ I3

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

Definition at line 26 of file strains_pomp.c.

◆ I3_0

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

Definition at line 20 of file strains_pomp.c.

◆ lik

#define lik   (__lik[0])

Definition at line 291 of file strains_pomp.c.

◆ ll

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

Definition at line 28 of file strains_pomp.c.

◆ N

#define N   (__p[__parindex[12]])

Definition at line 22 of file strains_pomp.c.

◆ node

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

Definition at line 32 of file strains_pomp.c.

◆ psi1

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

Definition at line 14 of file strains_pomp.c.

◆ psi2

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

Definition at line 15 of file strains_pomp.c.

◆ psi3

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

Definition at line 16 of file strains_pomp.c.

◆ R

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

Definition at line 27 of file strains_pomp.c.

◆ R_0

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

Definition at line 21 of file strains_pomp.c.

◆ S

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

Definition at line 23 of file strains_pomp.c.

◆ S_0

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

Definition at line 17 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 * penalty )
static

Definition at line 39 of file strains_pomp.c.

50 {
51 double event_rate = 0;
52 double alpha, disc;
53 *penalty = 0;
54 assert(S >= 0);
55 assert(I1 >= ellI1);
56 assert(I2 >= ellI2);
57 assert(I3 >= ellI3);
58 assert(ellI1 >= 0);
59 assert(ellI2 >= 0);
60 assert(ellI3 >= 0);
61 // 0: strain-1 transmission with saturation 0 or 1
62 alpha = Beta1*S*I1/N;
63 disc = (I1 > 0) ? ellI1*(ellI1-1)/I1/(I1+1) : 1;
64 event_rate += (*rate = alpha*(1-disc)); rate++;
65 *penalty += alpha*disc;
66 // 1: strain 1 recovery
67 alpha = gamma*I1;
68 if (I1 > ellI1) {
69 event_rate += (*rate = alpha); rate++;
70 } else {
71 *rate = 0; rate++;
72 *penalty += alpha;
73 }
74 // strain 1 sampling
75 alpha = psi1*I1;
76 *penalty += alpha;
77 // 2: strain-2 transmission with saturation 0 or 1
78 alpha = Beta2*S*I2/N;
79 disc = (I2 > 0) ? ellI2*(ellI2-1)/I2/(I2+1) : 1;
80 event_rate += (*rate = alpha*(1-disc)); rate++;
81 *penalty += alpha*disc;
82 // 3: strain 2 recovery
83 alpha = gamma*I2;
84 if (I2 > ellI2) {
85 event_rate += (*rate = alpha); rate++;
86 } else {
87 *rate = 0; rate++;
88 *penalty += alpha;
89 }
90 // strain 2 sampling
91 alpha = psi2*I2;
92 *penalty += alpha;
93 // 4: strain-3 transmission with saturation 0 or 1
94 alpha = Beta3*S*I3/N;
95 disc = (I3 > 0) ? ellI3*(ellI3-1)/I3/(I3+1) : 1;
96 event_rate += (*rate = alpha*(1-disc)); rate++;
97 *penalty += alpha*disc;
98 // 5: strain 3 recovery
99 alpha = gamma*I3;
100 if (I3 > ellI3) {
101 event_rate += (*rate = alpha); rate++;
102 } else {
103 *rate = 0; rate++;
104 *penalty += alpha;
105 }
106 // strain 3 sampling
107 alpha = psi3*I3;
108 *penalty += alpha;
109 assert(R_FINITE(event_rate));
110 return event_rate;
111}
N
#define N
Definition seirs_pomp.c:36
gamma
#define gamma
Definition seirs_pomp.c:29
S
#define S
Definition seirs_pomp.c:37
psi1
#define psi1
Definition strains_pomp.c:14
Beta3
#define Beta3
Definition strains_pomp.c:12
psi3
#define psi3
Definition strains_pomp.c:16
ellI1
#define ellI1
Definition strains_pomp.c:29
ellI2
#define ellI2
Definition strains_pomp.c:30
ellI3
#define ellI3
Definition strains_pomp.c:31
I2
#define I2
Definition strains_pomp.c:25
Beta2
#define Beta2
Definition strains_pomp.c:11
Beta1
#define Beta1
Definition strains_pomp.c:10
psi2
#define psi2
Definition strains_pomp.c:15
I3
#define I3
Definition strains_pomp.c:26
I1
#define I1
Definition strains_pomp.c:24

◆ 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 294 of file strains_pomp.c.

307 {
308 assert(!ISNAN(ll));
309 lik = (give_log) ? ll : exp(ll);
310}
lik
#define lik
Definition lbdp_pomp.c:164
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 140 of file strains_pomp.c.

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

123 {
124 double m = N/(S_0+I1_0+I2_0+I3_0+R_0);
125 S = nearbyint(S_0*m);
126 I1 = nearbyint(I1_0*m);
127 I2 = nearbyint(I2_0*m);
128 I3 = nearbyint(I3_0*m);
129 R = nearbyint(R_0*m);
130 ll = 0;
131 ellI1 = 0;
132 ellI2 = 0;
133 ellI3 = 0;
134 node = 0;
135}
I2_0
#define I2_0
Definition strains_pomp.c:19
S_0
#define S_0
Definition strains_pomp.c:17
I3_0
#define I3_0
Definition strains_pomp.c:20
I1_0
#define I1_0
Definition strains_pomp.c:18
R_0
#define R_0
Definition strains_pomp.c:21

Variable Documentation

◆ nrate

const int nrate = 6
static

Definition at line 4 of file strains_pomp.c.