phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
si2r.cc
Go to the documentation of this file.
1 // SI2R: Two-deme model of superspreading (C++)
2 #include "master.h"
3 #include "popul_proc.h"
4 #include "generics.h"
5 #include "internal.h"
6 
7 static int ordinary = 0;
8 static int superspreader = 1;
9 
11 typedef struct {
12  int S;
13  int I1;
14  int I2;
15  int R;
16  double N;
17 } si2r_state_t;
18 
20 typedef struct {
21  double Beta;
22  double mu;
23  double gamma;
24  double omega;
25  double psi1;
26  double psi2;
27  double sigma12;
28  double sigma21;
29  int S0;
30  int I0;
31  int R0;
33 
36 
37 template<>
38 std::string si2r_proc_t::yaml (std::string tab) const {
39  std::string t = tab + " ";
40  std::string p = tab + "parameter:\n"
41  + YAML_PARAM(Beta)
42  + YAML_PARAM(mu)
43  + YAML_PARAM(gamma)
44  + YAML_PARAM(omega)
45  + YAML_PARAM(psi1)
46  + YAML_PARAM(psi2)
47  + YAML_PARAM(sigma12)
48  + YAML_PARAM(sigma21)
49  + YAML_PARAM(S0)
50  + YAML_PARAM(I0)
51  + YAML_PARAM(R0);
52  std::string s = tab + "state:\n"
53  + YAML_STATE(S)
54  + YAML_STATE(I1)
55  + YAML_STATE(I2)
56  + YAML_STATE(R)
57  + YAML_STATE(N);
58  return p+s;
59 }
60 
61 template<>
62 void si2r_proc_t::update_params (double *p, int n) {
63  int m = 0;
64  PARAM_SET(Beta);
65  PARAM_SET(mu);
68  PARAM_SET(psi1);
69  PARAM_SET(psi2);
70  PARAM_SET(sigma12);
71  PARAM_SET(sigma21);
72  if (m != n) err("wrong number of parameters!");
73 }
74 
75 template<>
76 void si2r_proc_t::update_IVPs (double *p, int n) {
77  int m = 0;
78  PARAM_SET(S0);
79  PARAM_SET(I0);
80  PARAM_SET(R0);
81  if (m != n) err("wrong number of initial-value parameters!");
82 }
83 
84 template<>
85 double si2r_proc_t::event_rates (double *rate, int n) const {
86  int m = 0;
87  double total = 0;
88  RATE_CALC(params.Beta * state.S * state.I1 / state.N);
89  RATE_CALC(params.Beta * state.S * state.I2 / state.N);
90  RATE_CALC(params.gamma * state.I1);
91  RATE_CALC(params.gamma * state.I2);
92  RATE_CALC(params.psi1 * state.I1);
93  RATE_CALC(params.psi2 * state.I2);
94  RATE_CALC(params.sigma12 * state.I1);
95  RATE_CALC(params.sigma21 * state.I2);
96  RATE_CALC(params.omega * state.R);
97  if (m != n) err("wrong number of events!");
98  return total;
99 }
100 
101 template<>
102 void si2r_genealogy_t::rinit (void) {
103  state.S = params.S0;
104  state.I1 = params.I0;
105  state.I2 = 0;
106  state.R = params.R0;
107  state.N = double(params.S0+params.I0+params.R0);
108  graft(ordinary,params.I0);
109 }
110 
111 template<>
112 void si2r_genealogy_t::jump (int event) {
113  switch (event) {
114  case 0:
115  state.S -= 1; state.I1 += 1; birth(ordinary,ordinary);
116  break;
117  case 1:
118  {
119  int n = 1+int(rgeom(1.0/params.mu));
120  if (state.S >= n) {
121  state.S -= n; state.I1 += n;
123  } else {
124  birth(superspreader,ordinary,state.S);
125  state.I1 += state.S; state.S = 0;
126  }
127  }
128  break;
129  case 2:
130  state.I1 -= 1; state.R += 1; death(ordinary);
131  break;
132  case 3:
133  state.I2 -= 1; state.R += 1; death(superspreader);
134  break;
135  case 4:
136  sample(ordinary);
137  break;
138  case 5:
140  break;
141  case 6:
142  state.I1 -= 1; state.I2 += 1; migrate(ordinary,superspreader);
143  break;
144  case 7:
145  state.I1 += 1; state.I2 -= 1; migrate(superspreader,ordinary);
146  break;
147  case 8:
148  state.R -= 1; state.S += 1;
149  break;
150  default: // #nocov
151  assert(0); // #nocov
152  break; // #nocov
153  }
154 }
155 
Encodes the master process.
Definition: master.h:22
void birth(name_t i=0, name_t j=0, int n=1)
n births into deme j with parent in deme i
Definition: master.h:142
void sample(name_t i=0, int n=1)
sample in deme i
Definition: master.h:166
void death(name_t i=0)
death in deme i
Definition: master.h:153
void jump(int e)
makes a jump
Definition: lbdp.cc:72
void rinit(void)
initialize the state
Definition: lbdp.cc:66
void graft(name_t i=0, int m=1)
new root in deme i
Definition: master.h:159
void migrate(name_t i=0, name_t j=0)
migration from deme i to deme j
Definition: master.h:185
Population process class.
Definition: popul_proc.h:20
void update_IVPs(double *, int)
set initial-value parameters
Definition: lbdp.cc:48
void update_params(double *, int)
set parameters
Definition: lbdp.cc:39
std::string yaml(std::string tab) const
machine/human readable info
Definition: lbdp.cc:26
parameters_t params
Definition: popul_proc.h:35
double event_rates(double *rate, int n) const
compute event rates
Definition: lbdp.cc:55
state_t state
Definition: popul_proc.h:34
#define GENERICS(X, TYPE)
Definition: generics.h:149
#define err(...)
Definition: internal.h:18
#define mu
Definition: lbdp_pomp.c:5
#define n
Definition: lbdp_pomp.c:8
#define YAML_PARAM(X)
Definition: popul_proc.h:180
#define RATE_CALC(X)
Definition: popul_proc.h:179
#define YAML_STATE(X)
Definition: popul_proc.h:181
#define PARAM_SET(X)
Definition: popul_proc.h:177
#define N
Definition: seirs_pomp.c:32
#define R0
Definition: seirs_pomp.c:31
#define gamma
Definition: seirs_pomp.c:25
#define R
Definition: seirs_pomp.c:36
#define S0
Definition: seirs_pomp.c:28
#define Beta
Definition: seirs_pomp.c:23
#define I0
Definition: seirs_pomp.c:30
#define omega
Definition: seirs_pomp.c:27
#define S
Definition: seirs_pomp.c:33
static int ordinary
Definition: si2r.cc:7
static int superspreader
Definition: si2r.cc:8
SI2R process parameters.
Definition: si2r.cc:20
double gamma
Definition: si2r.cc:23
double psi1
Definition: si2r.cc:25
double omega
Definition: si2r.cc:24
double psi2
Definition: si2r.cc:26
double Beta
Definition: si2r.cc:21
double mu
Definition: si2r.cc:22
double sigma21
Definition: si2r.cc:28
double sigma12
Definition: si2r.cc:27
SI2R process state.
Definition: si2r.cc:11
int I1
Definition: si2r.cc:13
int I2
Definition: si2r.cc:14
int R
Definition: si2r.cc:15
double N
Definition: si2r.cc:16
int S
Definition: si2r.cc:12
#define psi1
#define I2
#define psi2
#define I1