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 
4 static 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 
26 static 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
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
#define N
Definition: sirs_pomp.c:13
#define R0
Definition: sirs_pomp.c:12
#define gamma
Definition: sirs_pomp.c:7
#define R
Definition: sirs_pomp.c:16
#define S0
Definition: sirs_pomp.c:10
#define I
Definition: sirs_pomp.c:15
#define Beta
Definition: sirs_pomp.c:6
#define I0
Definition: sirs_pomp.c:11
#define ellI
Definition: sirs_pomp.c:18
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
static const int nrate
Definition: sirs_pomp.c:4
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
#define lik
Definition: sirs_pomp.c:186
#define ll
Definition: sirs_pomp.c:17
#define omega
Definition: sirs_pomp.c:9
#define psi
Definition: sirs_pomp.c:8
#define EVENT_RATES
Definition: sirs_pomp.c:21
#define node
Definition: sirs_pomp.c:19
#define S
Definition: sirs_pomp.c:14