#include "pomplink.h"
#include "internal.h"
Go to the source code of this file.
|
#define | Beta (__p[__parindex[0]]) |
|
#define | sigma (__p[__parindex[1]]) |
|
#define | gamma (__p[__parindex[2]]) |
|
#define | psi (__p[__parindex[3]]) |
|
#define | omega (__p[__parindex[4]]) |
|
#define | S0 (__p[__parindex[5]]) |
|
#define | E0 (__p[__parindex[6]]) |
|
#define | I0 (__p[__parindex[7]]) |
|
#define | R0 (__p[__parindex[8]]) |
|
#define | N (__p[__parindex[9]]) |
|
#define | S (__x[__stateindex[0]]) |
|
#define | E (__x[__stateindex[1]]) |
|
#define | I (__x[__stateindex[2]]) |
|
#define | R (__x[__stateindex[3]]) |
|
#define | ll (__x[__stateindex[4]]) |
|
#define | node (__x[__stateindex[5]]) |
|
#define | ellE (__x[__stateindex[6]]) |
|
#define | ellI (__x[__stateindex[7]]) |
|
#define | COLOR (__x[__stateindex[8]]) |
|
#define | EVENT_RATES |
|
#define | lik (__lik[0]) |
|
|
static int | random_choice (double n) |
|
static void | change_color (double *color, int nsample, int n, int from, int to) |
|
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 *logpi, double *penalty) |
|
void | seirs_rinit (double *__x, const double *__p, double t0, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars) |
|
void | seirs_gill (double *__x, const double *__p, const int *__stateindex, const int *__parindex, const int *__covindex, const double *__covars, double t, double dt) |
|
void | seirs_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). More...
|
|
|
static const int | nrate = 7 |
|
◆ Beta
#define Beta (__p[__parindex[0]]) |
◆ COLOR
#define COLOR (__x[__stateindex[8]]) |
#define E (__x[__stateindex[1]]) |
◆ E0
#define E0 (__p[__parindex[6]]) |
◆ ellE
#define ellE (__x[__stateindex[6]]) |
◆ ellI
#define ellI (__x[__stateindex[7]]) |
◆ EVENT_RATES
Value:
__stateindex,__parindex,__covindex, \
__covars,rate,logpi,&penalty) \
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 *logpi, double *penalty)
Definition at line 43 of file seirs_pomp.c.
◆ gamma
#define gamma (__p[__parindex[2]]) |
#define I (__x[__stateindex[2]]) |
◆ I0
#define I0 (__p[__parindex[7]]) |
◆ lik
◆ ll
#define ll (__x[__stateindex[4]]) |
#define N (__p[__parindex[9]]) |
◆ node
#define node (__x[__stateindex[5]]) |
◆ omega
#define omega (__p[__parindex[4]]) |
◆ psi
#define psi (__p[__parindex[3]]) |
#define R (__x[__stateindex[3]]) |
◆ R0
#define R0 (__p[__parindex[8]]) |
#define S (__x[__stateindex[0]]) |
◆ S0
#define S0 (__p[__parindex[5]]) |
◆ sigma
#define sigma (__p[__parindex[1]]) |
◆ change_color()
static void change_color |
( |
double * |
color, |
|
|
int |
nsample, |
|
|
int |
n, |
|
|
int |
from, |
|
|
int |
to |
|
) |
| |
|
static |
Definition at line 10 of file seirs_pomp.c.
15 if (!ISNA(color[i]) && nearbyint(color[i]) == from)
n--;
19 assert(nearbyint(color[i]) == from);
◆ 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 * |
logpi, |
|
|
double * |
penalty |
|
) |
| |
|
static |
Definition at line 48 of file seirs_pomp.c.
61 double event_rate = 0;
67 pi = (
I > 0) ? 1-
ellI/
I : 0;
69 event_rate += (*rate = alpha*pi); rate++;
70 *logpi = log(pi); logpi++;
73 event_rate += (*rate = alpha*pi); rate++;
74 *logpi = log(pi)-log(
ellI); logpi++;
76 event_rate += (*rate = alpha*pi); rate++;
77 *logpi = log(pi)-log(
ellI); logpi++;
83 event_rate += (*rate = alpha*(1-pi)); rate++;
84 *logpi = log(1-pi); logpi++;
86 event_rate += (*rate = alpha*pi); rate++;
87 *logpi = log(pi)-log(
ellE); logpi++;
92 event_rate += (*rate = alpha); rate++;
100 event_rate += (*rate =
omega*
R); rate++;
104 assert(R_FINITE(event_rate));
◆ random_choice()
static int random_choice |
( |
double |
n | ) |
|
|
inlinestatic |
Definition at line 6 of file seirs_pomp.c.
7 return floor(R_unif_index(
n));
◆ seirs_dmeas()
void seirs_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 344 of file seirs_pomp.c.
359 lik = (give_log) ?
ll : exp(
ll);
◆ seirs_gill()
void seirs_gill |
( |
double * |
__x, |
|
|
const double * |
__p, |
|
|
const int * |
__stateindex, |
|
|
const int * |
__parindex, |
|
|
const int * |
__covindex, |
|
|
const double * |
__covars, |
|
|
double |
t, |
|
|
double |
dt |
|
) |
| |
Simulator for the latent-state process (rprocess).
This is the Gillespie algorithm applied to the solution of the filter equation for the SEIRS process.
Definition at line 139 of file seirs_pomp.c.
150 double tstep = 0.0, tmax = t + dt;
151 double *color = &
COLOR;
159 int parent = (int) nearbyint(
node);
164 assert(parent<=nnode);
167 int parlin = lineage[parent];
168 int parcol = color[parlin];
169 assert(parlin >= 0 && parlin <
nsample);
172 switch (nodetype[parent]) {
178 assert(sat[parent]==1);
179 int c = child[index[parent]];
180 assert(lineage[parent]==lineage[c]);
183 if (unif_rand() < x) {
184 color[lineage[c]] = 0;
188 color[lineage[c]] = 1;
195 if (unif_rand() < 0.5) {
196 color[lineage[c]] = 0;
200 color[lineage[c]] = 1;
216 if (sat[parent] == 1) {
217 int c = child[index[parent]];
218 color[lineage[c]] = 1;
220 }
else if (sat[parent] == 0) {
227 color[parlin] = R_NaReal;
239 assert(sat[parent]==2);
240 ll += (
S > 0 &&
I > 0) ? log(
Beta*
S/
N/(
E+1)) : R_NegInf;
244 int c1 = child[index[parent]];
245 int c2 = child[index[parent]+1];
247 assert(lineage[c1] != lineage[c2]);
248 assert(lineage[c1] != parlin || lineage[c2] != parlin);
249 assert(lineage[c1] == parlin || lineage[c2] == parlin);
250 if (unif_rand() < 0.5) {
251 color[lineage[c1]] = 0;
252 color[lineage[c2]] = 1;
254 color[lineage[c1]] = 1;
255 color[lineage[c2]] = 0;
267 double event_rate = 0;
271 tstep = exp_rand()/event_rate;
273 while (t + tstep < tmax) {
275 ll -= penalty*tstep + logpi[event];
332 tstep = exp_rand()/event_rate;
get_userdata_int_t * get_userdata_int
static int rcateg(double erate, double *rate, int nrate)
static void change_color(double *color, int nsample, int n, int from, int to)
static int random_choice(double n)
◆ seirs_rinit()
void seirs_rinit |
( |
double * |
__x, |
|
|
const double * |
__p, |
|
|
double |
t0, |
|
|
const int * |
__stateindex, |
|
|
const int * |
__parindex, |
|
|
const int * |
__covindex, |
|
|
const double * |
__covars |
|
) |
| |
Latent-state initializer (rinit component).
The state variables include S, E, I, R plus 'ellE' and 'ellI' (numbers of E- and I-deme lineages), the accumulated weight ('ll'), the current node number ('node'), and the coloring of each lineage ('COLOR').
Definition at line 114 of file seirs_pomp.c.
125 S = nearbyint(
S0*adj);
126 E = nearbyint(
E0*adj);
127 I = nearbyint(
I0*adj);
128 R = nearbyint(
R0*adj);
◆ nrate