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
7static int strain1 = 0;
8static int strain2 = 1;
9
11typedef struct {
12 int S;
13 int I1;
14 int I2;
15 int R;
16 double N;
18
20typedef 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
38template<>
39std::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)
47 + YAML_PARAM(sigma12)
48 + YAML_PARAM(sigma21)
50 + YAML_PARAM(S0)
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
63template<>
64void siir_proc_t::update_params (double *p, int n) {
65 int m = 0;
66 PARAM_SET(Beta1);
67 PARAM_SET(Beta2);
71 PARAM_SET(sigma12);
72 PARAM_SET(sigma21);
74 if (m != n) err("wrong number of parameters!");
75}
76
77template<>
78void siir_proc_t::update_IVPs (double *p, int n) {
79 int m = 0;
84 if (m != n) err("wrong number of initial-value parameters!");
85}
86
87template<>
88double 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
104template<>
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
115template<>
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:
132 break;
133 case 5:
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 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
double event_rates(double *rate, int n) const
Definition siir.cc:88
std::string yaml(std::string tab) const
Definition siir.cc:39
#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
master_t< siir_proc_t, 2 > siir_genealogy_t
Definition siir.cc:36
popul_proc_t< siir_state_t, siir_parameters_t, 9 > siir_proc_t
Definition siir.cc:35
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
double N
Definition siir.cc:16
#define I2_0
#define psi1
#define I2
#define I1_0
#define psi2
#define I1