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