phylopomp
Phylodynamics for POMPs
All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
s2i2r2.cc
Go to the documentation of this file.
1 // S2I2R2: Two-host infection model with waning, immigration, and demography. (C++)
2 #include "master.h"
3 #include "popul_proc.h"
4 #include "generics.h"
5 #include "internal.h"
6 
7 static int host1 = 0;
8 static int host2 = 1;
9 static int outside = 2;
10 
12 typedef struct {
13  int S1;
14  int I1;
15  int R1;
16  int S2;
17  int I2;
18  int R2;
19  double N1;
20  double N2;
22 
24 typedef struct {
25  double Beta11;
26  double Beta12;
27  double Beta22;
28  double gamma1;
29  double gamma2;
30  double psi1;
31  double psi2;
32  double omega1;
33  double omega2;
34  double b1;
35  double b2;
36  double d1;
37  double d2;
38  double iota1;
39  double iota2;
40  int S1_0;
41  int S2_0;
42  int I1_0;
43  int I2_0;
44  int R1_0;
45  int R2_0;
47 
50 
51 template<>
52 std::string s2i2r2_proc_t::yaml (std::string tab) const {
53  std::string t = tab + " ";
54  std::string p = tab + "parameter:\n"
60  + YAML_PARAM(psi1)
61  + YAML_PARAM(psi2)
64  + YAML_PARAM(b1)
65  + YAML_PARAM(b2)
66  + YAML_PARAM(d1)
67  + YAML_PARAM(d2)
68  + YAML_PARAM(iota1)
69  + YAML_PARAM(iota2)
70  + YAML_PARAM(S1_0)
71  + YAML_PARAM(S2_0)
72  + YAML_PARAM(I1_0)
73  + YAML_PARAM(I2_0)
74  + YAML_PARAM(R1_0)
75  + YAML_PARAM(R2_0);
76  std::string s = tab + "state:\n"
77  + YAML_STATE(S1)
78  + YAML_STATE(I1)
79  + YAML_STATE(R1)
80  + YAML_STATE(S2)
81  + YAML_STATE(I2)
82  + YAML_STATE(R2)
83  + YAML_STATE(N1)
84  + YAML_STATE(N2);
85  return p+s;
86 }
87 
88 template<>
89 void s2i2r2_proc_t::update_params (double *p, int n) {
90  int m = 0;
96  PARAM_SET(psi1);
97  PARAM_SET(psi2);
100  PARAM_SET(b1);
101  PARAM_SET(b2);
102  PARAM_SET(d1);
103  PARAM_SET(d2);
104  PARAM_SET(iota1);
105  PARAM_SET(iota2);
106  if (m != n) err("wrong number of parameters!");
107 }
108 
109 template<>
110 void s2i2r2_proc_t::update_IVPs (double *p, int n) {
111  int m = 0;
112  PARAM_SET(S1_0);
113  PARAM_SET(S2_0);
114  PARAM_SET(I1_0);
115  PARAM_SET(I2_0);
116  PARAM_SET(R1_0);
117  PARAM_SET(R2_0);
118  if (m != n) err("wrong number of initial-value parameters!");
119 }
120 
121 template<>
122 double s2i2r2_proc_t::event_rates (double *rate, int n) const {
123  int m = 0;
124  double total = 0;
125  RATE_CALC(params.Beta11 * state.I1 / state.N1 * state.S1);
126  RATE_CALC(params.Beta22 * state.I2 / state.N2 * state.S2);
127  RATE_CALC(params.Beta12 * state.I2 / state.N2 * state.S1);
128  RATE_CALC(params.gamma1 * state.I1);
129  RATE_CALC(params.gamma2 * state.I2);
130  RATE_CALC(params.omega1 * state.R1);
131  RATE_CALC(params.omega2 * state.R2);
132  RATE_CALC(params.psi1 * state.I1);
133  RATE_CALC(params.psi2 * state.I2);
134  RATE_CALC(params.iota1 * state.S1);
135  RATE_CALC(params.iota2 * state.S2);
136  RATE_CALC(params.d1 * state.S1);
137  RATE_CALC(params.d2 * state.S2);
138  RATE_CALC(params.d1 * state.I1);
139  RATE_CALC(params.d2 * state.I2);
140  RATE_CALC(params.d1 * state.R1);
141  RATE_CALC(params.d2 * state.R2);
142  RATE_CALC(params.b1 * state.N1);
143  RATE_CALC(params.b2 * state.N2);
144  if (m != n) err("wrong number of events!");
145  return total;
146 }
147 
148 template<>
149 void s2i2r2_genealogy_t::rinit (void) {
150  state.S1 = params.S1_0;
151  state.I1 = params.I1_0;
152  state.R1 = params.R1_0;
153  state.S2 = params.S2_0;
154  state.I2 = params.I2_0;
155  state.R2 = params.R2_0;
156  state.N1 = double(params.S1_0+params.I1_0+params.R1_0);
157  state.N2 = double(params.S2_0+params.I2_0+params.R2_0);
158  graft(host1,params.I1_0);
159  graft(host2,params.I2_0);
160 }
161 
162 template<>
163 void s2i2r2_genealogy_t::jump (int event) {
164  switch (event) {
165  case 0:
166  state.S1 -= 1; state.I1 += 1; birth(host1,host1);
167  break;
168  case 1:
169  state.S2 -= 1; state.I2 += 1; birth(host2,host2);
170  break;
171  case 2:
172  state.S1 -= 1; state.I1 += 1; birth(host2,host1);
173  break;
174  case 3:
175  state.I1 -= 1; state.R1 += 1; death(host1);
176  break;
177  case 4:
178  state.I2 -= 1; state.R2 += 1; death(host2);
179  break;
180  case 5:
181  state.R1 -= 1; state.S1 += 1;
182  break;
183  case 6:
184  state.R2 -= 1; state.S2 += 1;
185  break;
186  case 7:
187  sample(host1);
188  break;
189  case 8:
190  sample(host2);
191  break;
192  case 9:
193  state.S1 -= 1; state.I1 += 1; graft(outside); migrate(outside,host1);
194  break;
195  case 10:
196  state.S2 -= 1; state.I2 += 1; graft(outside); migrate(outside,host2);
197  break;
198  case 11:
199  state.S1 -= 1; state.N1 -= 1;
200  break;
201  case 12:
202  state.S2 -= 1; state.N2 -= 1;
203  break;
204  case 13:
205  state.I1 -= 1; state.N1 -= 1; death(host1);
206  break;
207  case 14:
208  state.I2 -= 1; state.N2 -= 1; death(host2);
209  break;
210  case 15:
211  state.R1 -= 1; state.N1 -= 1;
212  break;
213  case 16:
214  state.R2 -= 1; state.N2 -= 1;
215  break;
216  case 17:
217  state.S1 += 1; state.N1 += 1;
218  break;
219  case 18:
220  state.S2 += 1; state.N2 += 1;
221  break;
222  default: // #nocov
223  assert(0); // #nocov
224  break; // #nocov
225  }
226 }
227 
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
static int host2
Definition: s2i2r2.cc:8
static int host1
Definition: s2i2r2.cc:7
static int outside
Definition: s2i2r2.cc:9
S2I2R2 process parameters.
Definition: s2i2r2.cc:24
S2I2R2 process state.
Definition: s2i2r2.cc:12
double N1
Definition: s2i2r2.cc:19
double N2
Definition: s2i2r2.cc:20
#define I2_0
#define psi1
#define S2_0
#define d1
#define Beta12
#define S1_0
#define b1
#define Beta22
#define S1
#define R2_0
#define R2
#define R1
#define I2
#define R1_0
#define Beta11
#define I1_0
#define N1
#define omega2
#define gamma1
#define N2
#define b2
#define S2
#define psi2
#define d2
#define I1
#define gamma2
#define omega1