phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
lbdp_pomp.c
Go to the documentation of this file.
1 #include "pomplink.h"
2 #include "internal.h"
3 
4 #define lambda (__p[__parindex[0]])
5 #define mu (__p[__parindex[1]])
6 #define psi (__p[__parindex[2]])
7 #define n0 (__p[__parindex[3]])
8 #define n (__x[__stateindex[0]])
9 #define ll (__x[__stateindex[1]])
10 #define ell (__x[__stateindex[2]])
11 #define node (__x[__stateindex[3]])
12 
13 #define EVENT_RATES \
14  event_rates(__x,__p,t, \
15  __stateindex,__parindex,rate,&penalty) \
16 
17 static double event_rates
18 (
19  double *__x,
20  const double *__p,
21  double t,
22  const int *__stateindex,
23  const int *__parindex,
24  double *rate,
25  double *penalty
26  ) {
27  double event_rate = 0;
28  double alpha, disc;
29  *penalty = 0;
30  assert(n >= ell);
31  assert(ell >= 0);
32  // birth with saturation 0 or 1
33  alpha = lambda*n;
34  disc = (n > 0) ? ell*(ell-1)/n/(n+1) : 1;
35  event_rate += (*rate = alpha*(1-disc)); rate++;
36  *penalty += alpha*disc;
37  // death
38  alpha = mu*n;
39  if (n > ell) {
40  event_rate += (*rate = alpha); rate++;
41  } else {
42  *rate = 0; rate++;
43  *penalty += alpha;
44  }
45  // sampling
46  *penalty += psi*n;
47  assert(R_FINITE(event_rate));
48  return event_rate;
49 }
50 
53 (
54  double *__x,
55  const double *__p,
56  double t,
57  const int *__stateindex,
58  const int *__parindex,
59  const int *__covindex,
60  const double *__covars
61  ){
62  n = nearbyint(n0);
63  ll = 0;
64  ell = 0;
65  node = 0;
66 }
67 
72 (
73  double *__x,
74  const double *__p,
75  const int *__stateindex,
76  const int *__parindex,
77  const int *__covindex,
78  const double *__covars,
79  double t,
80  double dt
81  ){
82  double tstep = 0.0, tmax = t + dt;
83  const int *nodetype = get_userdata_int("nodetype");
84  const int *sat = get_userdata_int("saturation");
85  int parent = (int) nearbyint(node);
86 
87 #ifndef NDEBUG
88  int nnode = *get_userdata_int("nnode");
89  assert(parent>=0);
90  assert(parent<=nnode);
91 #endif
92 
93  // singular portion of filter equation
94  switch (nodetype[parent]) {
95  default: // non-genealogical event
96  break;
97  case 0: // root
98  ll = 0;
99  ell += 1;
100  break;
101  case 1: // sample
102  ll = 0;
103  assert(n >= ell);
104  assert(ell >= 0);
105  if (sat[parent] == 1) { // s=1
106  ll += log(psi);
107  } else if (sat[parent] == 0) { // s=0
108  ell -= 1;
109  ll += log(psi*(n-ell));
110  } else {
111  assert(0); // #nocov
112  ll += R_NegInf; // #nocov
113  }
114  break;
115  case 2: // branch point s=2
116  ll = 0;
117  assert(n >= 0);
118  assert(ell > 0);
119  assert(sat[parent]==2);
120  n += 1; ell += 1;
121  ll += log(2*lambda/n);
122  break;
123  }
124 
125  if (tmax > t) {
126 
127  // Gillespie steps:
128  int event;
129  double penalty = 0;
130  double rate[2];
131 
132  double event_rate = EVENT_RATES;
133  tstep = exp_rand()/event_rate;
134 
135  while (t + tstep < tmax) {
136  event = rcateg(event_rate,rate,2);
137  ll -= penalty*tstep;
138  switch (event) {
139  case 0: // birth
140  n += 1;
141  break;
142  case 1: // death
143  n -= 1;
144  break;
145  default: // #nocov
146  assert(0); // #nocov
147  break; // #nocov
148  }
149  t += tstep;
150  event_rate = EVENT_RATES;
151  tstep = exp_rand()/event_rate;
152  }
153  tstep = tmax - t;
154  ll -= penalty*tstep;
155  }
156  node += 1;
157 }
158 
159 # define lik (__lik[0])
160 
163 (
164  double *__lik,
165  const double *__y,
166  const double *__x,
167  const double *__p,
168  int give_log,
169  const int *__obsindex,
170  const int *__stateindex,
171  const int *__parindex,
172  const int *__covindex,
173  const double *__covars,
174  double t
175  ) {
176  assert(!ISNAN(ll));
177  lik = (give_log) ? ll : exp(ll);
178 }
get_userdata_int_t * get_userdata_int
Definition: init.c:7
static int rcateg(double erate, double *rate, int nrate)
Definition: internal.h:76
#define mu
Definition: lbdp_pomp.c:5
void lbdp_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: lbdp_pomp.c:53
#define n0
Definition: lbdp_pomp.c:7
void lbdp_gill(double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt)
Definition: lbdp_pomp.c:72
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
#define ell
Definition: lbdp_pomp.c:10
#define lik
Definition: lbdp_pomp.c:159
#define n
Definition: lbdp_pomp.c:8
#define ll
Definition: lbdp_pomp.c:9
void lbdp_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: lbdp_pomp.c:163
#define lambda
Definition: lbdp_pomp.c:4
#define psi
Definition: lbdp_pomp.c:6
#define EVENT_RATES
Definition: lbdp_pomp.c:13
#define node
Definition: lbdp_pomp.c:11