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
17static 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