phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
siir.cc
Go to the documentation of this file.
1 // SIIR: Two-strain SIR model (C++)
2 #include "master.h"
3 #include "popul_proc.h"
4 #include "generics.h"
5 #include "internal.h"
6 
7 static int strain1 = 0;
8 static int strain2 = 1;
9 
11 typedef struct {
12  int S;
13  int I1;
14  int I2;
15  int R;
16  double N;
17 } siir_state_t;
18 
20 typedef struct {
21  double Beta1;
22  double Beta2;
23  double gamma;
24  double psi1;
25  double psi2;
26  double sigma12;
27  double sigma21;
28  double omega;
29  int S0;
30  int I1_0;
31  int I2_0;
32  int R0;
34 
37 
38 template<>
39 std::string siir_proc_t::yaml (std::string tab) const {
40  std::string t = tab + " ";
41  std::string p = tab + "parameter:\n"
42  + YAML_PARAM(Beta1)
43  + YAML_PARAM(Beta2)
44  + YAML_PARAM(gamma)
45  + YAML_PARAM(psi1)
46  + YAML_PARAM(psi2)
47  + YAML_PARAM(sigma12)
48  + YAML_PARAM(sigma21)
49  + YAML_PARAM(omega)
50  + YAML_PARAM(S0)
51  + YAML_PARAM(I1_0)
52  + YAML_PARAM(I2_0)
53  + YAML_PARAM(R0);
54  std::string s = tab + "state:\n"
55  + YAML_STATE(S)
56  + YAML_STATE(I1)
57  + YAML_STATE(I2)
58  + YAML_STATE(R)
59  + YAML_STATE(N);
60  return p+s;
61 }
62 
63 template<>
64 void siir_proc_t::update_params (double *p, int n) {
65  int m = 0;
66  PARAM_SET(Beta1);
67  PARAM_SET(Beta2);
69  PARAM_SET(psi1);
70  PARAM_SET(psi2);
71  PARAM_SET(sigma12);
72  PARAM_SET(sigma21);
74  if (m != n) err("wrong number of parameters!");
75 }
76 
77 template<>
78 void siir_proc_t::update_IVPs (double *p, int n) {
79  int m = 0;
80  PARAM_SET(S0);
81  PARAM_SET(I1_0);
82  PARAM_SET(I2_0);
83  PARAM_SET(R0);
84  if (m != n) err("wrong number of initial-value parameters!");
85 }
86 
87 template<>
88 double siir_proc_t::event_rates (double *rate, int n) const {
89  int m = 0;
90  double total = 0;
91  RATE_CALC(params.Beta1 * state.S * state.I1 / state.N);
92  RATE_CALC(params.Beta2 * state.S * state.I2 / state.N);
93  RATE_CALC(params.gamma * state.I1);
94  RATE_CALC(params.gamma * state.I2);
95  RATE_CALC(params.psi1 * state.I1);
96  RATE_CALC(params.psi2 * state.I2);
97  RATE_CALC(params.sigma12 * state.I1);
98  RATE_CALC(params.sigma21 * state.I2);
99  RATE_CALC(params.omega * state.R);
100  if (m != n) err("wrong number of events!");
101  return total;
102 }
103 
104 template<>
105 void siir_genealogy_t::rinit (void) {
106  state.S = params.S0;
107  state.I1 = params.I1_0;
108  state.I2 = params.I2_0;
109  state.R = params.R0;
110  state.N = double(params.S0+params.I1_0+params.I2_0+params.R0);
111  graft(strain1,params.I1_0);
112  graft(strain2,params.I2_0);
113 }
114 
115 template<>
116 void siir_genealogy_t::jump (int event) {
117  switch (event) {
118  case 0:
119  state.S -= 1; state.I1 += 1; birth(strain1,strain1);
120  break;
121  case 1:
122  state.S -= 1; state.I2 += 1; birth(strain2,strain2);
123  break;
124  case 2:
125  state.I1 -= 1; state.R += 1; death(strain1);
126  break;
127  case 3:
128  state.I2 -= 1; state.R += 1; death(strain2);
129  break;
130  case 4:
131  sample(strain1);
132  break;
133  case 5:
134  sample(strain2);
135  break;
136  case 6:
137  state.I1 -= 1; state.I2 += 1; migrate(strain1,strain2);
138  break;
139  case 7:
140  state.I1 += 1; state.I2 -= 1; migrate(strain2,strain1);
141  break;
142  case 8:
143  state.S += 1; state.R -= 1;
144  break;
145  default: // #nocov
146  assert(0); // #nocov
147  break; // #nocov
148  }
149 }
150 
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 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 omega
Definition: seirs_pomp.c:27
#define S
Definition: seirs_pomp.c:33
static int strain1
Definition: siir.cc:7
static int strain2
Definition: siir.cc:8
SIIR process parameters.
Definition: siir.cc:20
double psi2
Definition: siir.cc:25
double sigma21
Definition: siir.cc:27
double omega
Definition: siir.cc:28
double Beta1
Definition: siir.cc:21
double sigma12
Definition: siir.cc:26
double psi1
Definition: siir.cc:24
double Beta2
Definition: siir.cc:22
double gamma
Definition: siir.cc:23
SIIR process state.
Definition: siir.cc:11
int R
Definition: siir.cc:15
int S
Definition: siir.cc:12
int I2
Definition: siir.cc:14
double N
Definition: siir.cc:16
int I1
Definition: siir.cc:13
#define I2_0
#define psi1
#define I2
#define I1_0
#define psi2
#define I1