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   0
 
#define STRAIN2   1
 
#define STRAIN3   2
 
#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:18

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 287 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   0

Definition at line 6 of file strains_pomp.c.

◆ STRAIN2

#define STRAIN2   1

Definition at line 7 of file strains_pomp.c.

◆ STRAIN3

#define STRAIN3   2

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:32
gamma
#define gamma
Definition seirs_pomp.c:25
S
#define S
Definition seirs_pomp.c:33
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 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:160
ll
#define ll
Definition lbdp_pomp.c:9

◆ 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
156 int parent = (int) nearbyint(node);
157
158#ifndef NDEBUG
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
169 break;
170 case 0: // root
171 switch (deme[parent]) {
172 case STRAIN1:
173 ellI1 += 1; break;
174 case STRAIN2:
175 ellI2 += 1; break;
176 case STRAIN3:
177 ellI3 += 1; break;
178 default:
179 assert(0); break;
180 }
181 break;
182 case 1: // sample
183 assert(sat[parent]==0);
184 switch (deme[parent]) {
185 case STRAIN1:
186 assert(I1 >= ellI1);
187 assert(ellI1 >= 0);
188 ellI1 -= 1; I1 -= 1;
189 ll += log(psi1);
190 break;
191 case STRAIN2:
192 assert(I2 >= ellI2);
193 assert(ellI2 >= 0);
194 ellI2 -= 1; I2 -= 1;
195 ll += log(psi2);
196 break;
197 case STRAIN3:
198 assert(I3 >= ellI3);
199 assert(ellI3 >= 0);
200 ellI3 -= 1; I3 -= 1;
201 ll += log(psi3);
202 break;
203 default:
204 assert(0); break;
205 }
206 break;
207 case 2: // branch point s=(1,1)
208 assert(S >= 0);
209 if (sat[parent]!=2) break;
210 switch (deme[parent]) {
211 case STRAIN1:
212 assert(I1 >= 0);
213 assert(ellI1 > 0);
214 ll += (I1 > 0 && I1 >= ellI1) ? log(Beta1*S*I1/N) : R_NegInf;
215 S -= 1; I1 += 1; ellI1 += 1;
216 ll -= log(I1*(I1-1)/2);
217 break;
218 case STRAIN2:
219 assert(I2 >= 0);
220 assert(ellI2 > 0);
221 ll += (I2 > 0 && I2 >= ellI2) ? log(Beta2*S*I2/N) : R_NegInf;
222 S -= 1; I2 += 1; ellI2 += 1;
223 ll -= log(I2*(I2-1)/2);
224 break;
225 case STRAIN3:
226 assert(I3 >= 0);
227 assert(ellI3 > 0);
228 ll += (I3 > 0 && I3 >= ellI3) ? log(Beta3*S*I3/N) : R_NegInf;
229 S -= 1; I3 += 1; ellI3 += 1;
230 ll -= log(I3*(I3-1)/2);
231 break;
232 default:
233 assert(0);
234 break;
235 }
236 S = (S > 0) ? S : 0;
237 break;
238 }
239
240 if (tmax > t) {
241
242 // take Gillespie steps to the end of the interval:
243 int event;
244 double penalty = 0;
245 double rate[nrate];
246
247 double event_rate = EVENT_RATES;
248 tstep = exp_rand()/event_rate;
249
250 while (t + tstep < tmax) {
251 event = rcateg(event_rate,rate,nrate);
252 assert(event>=0 && event<nrate);
253 ll -= penalty*tstep;
254 switch (event) {
255 case 0: // strain-1 transmission
256 S -= 1; I1 += 1;
257 break;
258 case 1: // strain-1 recovery
259 I1 -= 1; R += 1;
260 break;
261 case 2: // strain-2 transmission
262 S -= 1; I2 += 1;
263 break;
264 case 3: // strain-2 recovery
265 I2 -= 1; R += 1;
266 break;
267 case 4: // strain-3 transmission
268 S -= 1; I3 += 1;
269 break;
270 case 5: // strain-3 recovery
271 I3 -= 1; R += 1;
272 break;
273 default: // #nocov
274 assert(0); // #nocov
275 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:76
deme
static const int deme
Definition lbdp.cc:7
EVENT_RATES
#define EVENT_RATES
Definition lbdp_pomp.c:13
node
#define node
Definition lbdp_pomp.c:11
R
#define R
Definition seirs_pomp.c:36
nrate
static const int nrate
Definition seirs_pomp.c:4
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.