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

Go to the source code of this file.

Macros

#define Beta   (__p[__parindex[0]])
 
#define gamma   (__p[__parindex[1]])
 
#define psi   (__p[__parindex[2]])
 
#define omega   (__p[__parindex[3]])
 
#define S0   (__p[__parindex[4]])
 
#define I0   (__p[__parindex[5]])
 
#define R0   (__p[__parindex[6]])
 
#define N   (__p[__parindex[7]])
 
#define S   (__x[__stateindex[0]])
 
#define I   (__x[__stateindex[1]])
 
#define R   (__x[__stateindex[2]])
 
#define ll   (__x[__stateindex[3]])
 
#define ellI   (__x[__stateindex[4]])
 
#define node   (__x[__stateindex[5]])
 
#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 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).
 
void sirs_gill (double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)
 
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).
 

Variables

static const int nrate = 3
 

Macro Definition Documentation

◆ Beta

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

Definition at line 6 of file sirs_pomp.c.

◆ ellI

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

Definition at line 18 of file sirs_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 21 of file sirs_pomp.c.

21#define EVENT_RATES \
22 event_rates(__x,__p,t, \
23 __stateindex,__parindex,__covindex, \
24 __covars,rate,&penalty) \
25

◆ gamma

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

Definition at line 7 of file sirs_pomp.c.

◆ I

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

Definition at line 15 of file sirs_pomp.c.

◆ I0

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

Definition at line 11 of file sirs_pomp.c.

◆ lik

#define lik   (__lik[0])

Definition at line 186 of file sirs_pomp.c.

◆ ll

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

Definition at line 17 of file sirs_pomp.c.

◆ N

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

Definition at line 13 of file sirs_pomp.c.

◆ node

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

Definition at line 19 of file sirs_pomp.c.

◆ omega

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

Definition at line 9 of file sirs_pomp.c.

◆ psi

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

Definition at line 8 of file sirs_pomp.c.

◆ R

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

Definition at line 16 of file sirs_pomp.c.

◆ R0

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

Definition at line 12 of file sirs_pomp.c.

◆ S

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

Definition at line 14 of file sirs_pomp.c.

◆ S0

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

Definition at line 10 of file sirs_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 26 of file sirs_pomp.c.

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}
psi
#define psi
Definition lbdp_pomp.c:6
N
#define N
Definition seirs_pomp.c:32
gamma
#define gamma
Definition seirs_pomp.c:25
R
#define R
Definition seirs_pomp.c:36
I
#define I
Definition seirs_pomp.c:35
Beta
#define Beta
Definition seirs_pomp.c:23
ellI
#define ellI
Definition seirs_pomp.c:40
omega
#define omega
Definition seirs_pomp.c:27
S
#define S
Definition seirs_pomp.c:33

◆ sirs_dmeas()

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 at line 189 of file sirs_pomp.c.

202 {
203 assert(!ISNAN(ll));
204 lik = (give_log) ? ll : exp(ll);
205}
lik
#define lik
Definition lbdp_pomp.c:159
ll
#define ll
Definition lbdp_pomp.c:9

◆ sirs_gill()

void sirs_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 90 of file sirs_pomp.c.

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}
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
EVENT_RATES
#define EVENT_RATES
Definition lbdp_pomp.c:13
node
#define node
Definition lbdp_pomp.c:11
nrate
static const int nrate
Definition seirs_pomp.c:4
Here is the call graph for this function:

◆ sirs_rinit()

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 at line 68 of file sirs_pomp.c.

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}
R0
#define R0
Definition seirs_pomp.c:31
S0
#define S0
Definition seirs_pomp.c:28
I0
#define I0
Definition seirs_pomp.c:30

Variable Documentation

◆ nrate

const int nrate = 3
static

Definition at line 4 of file sirs_pomp.c.