Commit 85620fed authored by Thomas  OBADIA's avatar Thomas OBADIA 💬

Merge branch 'devel-reconcile-v4' into 'master'

Merge devel-reconcile-v4 into default branch

See merge request !3
parents 52e6322c 9983f393
......@@ -4,6 +4,7 @@ Output/
Release/
Debug/
Pv_mod/Debug/
Pv_mod/Pv_model.o
*.vcxproj*
*.sln
IB_mod_batch.bat
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -39,6 +39,7 @@ public:
////////////////////////////////////////////////////
// 4.1.1.1. Class constructor
Individual(double a, double zeta)
{
age = a;
......@@ -52,6 +53,7 @@ public:
// Delete unwanted copy constructors
Individual(Individual&) = delete;
void operator= (Individual&) = delete;
// Allow default move constructors
Individual(Individual&&) = default;
Individual& operator= (Individual&&) = default;
......@@ -79,9 +81,10 @@ public:
double age; // Person's age
double zeta_het; // Heterogeneity in exposure to mosquitoes
bool gender; // 0 = male; 1 = female
bool G6PD_def; // Is the person G6PD deficient? 0 = normal; 1 = deficient (homozygous); 2 = deficient (heterozygous - female only)
double G6PD_activity; // G6PD activity - quantitative measurement
int gender; // 0 = male (domestic); 1 = male (occupational); 2 = female
int G6PD_def; // Is the person G6PD deficient? 0 = normal; 1 = deficient (homozygous); 2 = deficient (heterozygous - female only)
double G6PD_activity; // G6PD activity - quantitative measurement (true value)
double G6PD_read; // G6PD activity - quantitative measurement (from dx)
bool CYP2D6; // Does the person have CYP2D6 phenotype? 0 = normal; 1 = low metabolizer
......@@ -120,7 +123,6 @@ public:
int Hyp;
////////////////////////////////////////////////////
// Indicator for competing hazards move
......@@ -133,10 +135,14 @@ public:
bool I_PCR_new; // New PCR-detectable infection (I_LM, I_D & T included here)
bool I_LM_new; // New LM-detectable infection (I_D & T included here)
bool I_D_new; // New clinical episode (treated or untreated)
bool BS_treat; // Indicator for chloroquine (or other blood-stage drug) treatment
bool CQ_treat; // Indicator for chloroquine (or other blood-stage drug) treatment
bool PQ_treat; // Indicator for primaquine treatment
bool TQ_treat; // Indicator for tafenoquine treatment
bool G6PD_test; // New G6PD test administered
bool CQ_effective; // Effective blood-stage treatment (CQ)
bool PQ_effective; // Effective 8-aminoquinoline treatment (PQ)
bool PQ_overtreat; // Over-treatment with 8-aminoquinolines (PQ) - defined as someone without hypnozoites being treated
bool PQ_overtreat_9m; // Over-treatment with 8-aminoquinolines (PQ) - defined as someone without BS infection in last 9 mths being treated
......@@ -186,9 +192,9 @@ public:
bool LLIN; // Does the person have an LLIN?
double LLIN_age; // How old is the LLIN
double r_LLIN[N_spec]; // repellency
double d_LLIN[N_spec]; // killing effect
double s_LLIN[N_spec]; // survival
double r_LLIN[N_mosq]; // repellency
double d_LLIN[N_mosq]; // killing effect
double s_LLIN[N_mosq]; // survival
///////////////////////////////
......@@ -197,9 +203,9 @@ public:
bool IRS; // Is the person protected by IRS
double IRS_age; // How long ago was their house sprayed?
double r_IRS[N_spec]; // repellency
double d_IRS[N_spec]; // killing effect
double s_IRS[N_spec]; // survival
double r_IRS[N_mosq]; // repellency
double d_IRS[N_mosq]; // killing effect
double s_IRS[N_mosq]; // survival
///////////////////////////////
......@@ -220,9 +226,9 @@ public:
///////////////////////////////////////////////////
// Individual-level effect of vector control
double z_VC[N_spec]; // probability of mosquito being repelled from this individual during a single feeding attempt
double y_VC[N_spec]; // probability of mosquito feeding on this individual during a single attempt
double w_VC[N_spec]; // probability of mosquito feeding and surviving on this individual during a single feeding attempt
double z_VC[N_mosq]; // probability of mosquito being repelled from this individual during a single feeding attempt
double y_VC[N_mosq]; // probability of mosquito feeding on this individual during a single attempt
double w_VC[N_mosq]; // probability of mosquito feeding and surviving on this individual during a single feeding attempt
};
#endif
This diff is collapsed.
This diff is collapsed.
......@@ -23,22 +23,22 @@
// //
//////////////////////////////////////////////////////////////////////////////
void mosq_derivs(const double t, double(&yM)[N_spec][N_M_comp], double(&dyMdt)[N_spec][N_M_comp], Params& theta, Population& POP)
void mosq_derivs(const double t, double(&yM)[N_mosq][N_M_comp], double(&dyMdt)[N_mosq][N_M_comp], Params& theta, Population& POP)
{
double Karry_seas_inv[N_spec];
double Karry_seas_inv[N_mosq];
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
Karry_seas_inv[g] = 1.0 / (theta.Karry[g] * (theta.dry_seas[g] + (1 - theta.dry_seas[g])*pow(0.5 + 0.5*cos(0.01721421*(t - theta.t_peak_seas[g])), theta.kappa_seas[g]) / theta.denom_seas[g]));
Karry_seas_inv[v] = 1.0 / (theta.Karry[v] * (theta.dry_seas[v] + (1 - theta.dry_seas[v])*pow(0.5 + 0.5*cos(0.01721421*(t - theta.t_peak_seas[v])), theta.kappa_seas[v]) / theta.denom_seas[v]));
//Karry_seas_inv[g] = 1.0/theta.Karry[g];
//Karry_seas_inv[v] = 1.0/theta.Karry[v];
dyMdt[g][0] = POP.beta_VC[g] * (yM[g][3] + yM[g][4] + yM[g][5]) - yM[g][0] / theta.d_E_larvae - yM[g][0] * theta.mu_E0*(1.0 + (yM[g][0] + yM[g][1])*Karry_seas_inv[g]);
dyMdt[g][1] = yM[g][0] / theta.d_E_larvae - yM[g][1] / theta.d_L_larvae - yM[g][1] * theta.mu_L0*(1.0 + theta.gamma_larvae*(yM[g][0] + yM[g][1])*Karry_seas_inv[g]);
dyMdt[g][2] = yM[g][1] / theta.d_L_larvae - yM[g][2] / theta.d_pupae - yM[g][2] * theta.mu_P;
dyMdt[g][3] = 0.5*yM[g][2] / theta.d_pupae - theta.lam_M[g] * yM[g][3] - POP.mu_M_VC[g] * yM[g][3];
dyMdt[g][4] = +theta.lam_M[g] * yM[g][3] - theta.lam_S_M_track[g][0] * POP.exp_muM_tauM_VC[g] - POP.mu_M_VC[g] * yM[g][4];
dyMdt[g][5] = +theta.lam_S_M_track[g][0] * POP.exp_muM_tauM_VC[g] - POP.mu_M_VC[g] * yM[g][5];
dyMdt[v][0] = POP.beta_VC[v] * (yM[v][3] + yM[v][4] + yM[v][5]) - yM[v][0] / theta.d_E_larvae - yM[v][0] * theta.mu_E0*(1.0 + (yM[v][0] + yM[v][1])*Karry_seas_inv[v]);
dyMdt[v][1] = yM[v][0] / theta.d_E_larvae - yM[v][1] / theta.d_L_larvae - yM[v][1] * theta.mu_L0*(1.0 + theta.gamma_larvae*(yM[v][0] + yM[v][1])*Karry_seas_inv[v]);
dyMdt[v][2] = yM[v][1] / theta.d_L_larvae - yM[v][2] / theta.d_pupae - yM[v][2] * theta.mu_P;
dyMdt[v][3] = 0.5*yM[v][2] / theta.d_pupae - theta.lam_M[v] * yM[v][3] - POP.mu_M_VC[v] * yM[v][3];
dyMdt[v][4] = +theta.lam_M[v] * yM[v][3] - theta.lam_S_M_track[v][0] * POP.exp_muM_tauM_VC[v] - POP.mu_M_VC[v] * yM[v][4];
dyMdt[v][5] = +theta.lam_S_M_track[v][0] * POP.exp_muM_tauM_VC[v] - POP.mu_M_VC[v] * yM[v][5];
}
}
......@@ -48,9 +48,9 @@ void mosq_derivs(const double t, double(&yM)[N_spec][N_M_comp], double(&dyMdt)[N
// //
//////////////////////////////////////////////////////////////////////////////
void mosq_rk4(const double t, const double t_step_mosq, double(&yM)[N_spec][N_M_comp], Params& theta, Population& POP)
void mosq_rk4(const double t, const double t_step_mosq, double(&yM)[N_mosq][N_M_comp], Params& theta, Population& POP)
{
double k1_yM[N_spec][N_M_comp], k2_yM[N_spec][N_M_comp], k3_yM[N_spec][N_M_comp], k4_yM[N_spec][N_M_comp], yM_temp[N_spec][N_M_comp];
double k1_yM[N_mosq][N_M_comp], k2_yM[N_mosq][N_M_comp], k3_yM[N_mosq][N_M_comp], k4_yM[N_mosq][N_M_comp], yM_temp[N_mosq][N_M_comp];
//////////////////////////
......@@ -62,11 +62,11 @@ void mosq_rk4(const double t, const double t_step_mosq, double(&yM)[N_spec][N_M_
//////////////////////////
// step 2
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
for (int k = 0; k < N_M_comp; k++)
{
yM_temp[g][k] = yM[g][k] + 0.5*t_step_mosq*k1_yM[g][k];
yM_temp[v][k] = yM[v][k] + 0.5*t_step_mosq*k1_yM[v][k];
}
}
......@@ -76,11 +76,11 @@ void mosq_rk4(const double t, const double t_step_mosq, double(&yM)[N_spec][N_M_
//////////////////////////
// step 3
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
for (int k = 0; k < N_M_comp; k++)
{
yM_temp[g][k] = yM[g][k] + 0.5*t_step_mosq*k2_yM[g][k];
yM_temp[v][k] = yM[v][k] + 0.5*t_step_mosq*k2_yM[v][k];
}
}
......@@ -90,11 +90,11 @@ void mosq_rk4(const double t, const double t_step_mosq, double(&yM)[N_spec][N_M_
//////////////////////////
// step 4
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
for (int k = 0; k < N_M_comp; k++)
{
yM_temp[g][k] = yM[g][k] + t_step_mosq * k3_yM[g][k];
yM_temp[v][k] = yM[v][k] + t_step_mosq * k3_yM[v][k];
}
}
......@@ -104,11 +104,11 @@ void mosq_rk4(const double t, const double t_step_mosq, double(&yM)[N_spec][N_M_
//////////////////////////
// output
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
for (int k = 0; k < N_M_comp; k++)
{
yM[g][k] = yM[g][k] + t_step_mosq * k1_yM[g][k] / 6.0 + t_step_mosq * k2_yM[g][k] / 3.0 + t_step_mosq * k3_yM[g][k] / 3.0 + t_step_mosq * k4_yM[g][k] / 6.0;
yM[v][k] = yM[v][k] + t_step_mosq * k1_yM[v][k] / 6.0 + t_step_mosq * k2_yM[v][k] / 3.0 + t_step_mosq * k3_yM[v][k] / 3.0 + t_step_mosq * k4_yM[v][k] / 6.0;
}
}
......@@ -130,13 +130,13 @@ void mosquito_step(double t, Params& theta, Population& POP)
//////////////////////////////////
// Set up mosquito state vector
double yM[N_spec][N_M_comp];
double yM[N_mosq][N_M_comp];
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
for (int k = 0; k < N_M_comp; k++)
{
yM[g][k] = POP.yM[g][k];
yM[v][k] = POP.yM[v][k];
}
}
......@@ -147,16 +147,16 @@ void mosquito_step(double t, Params& theta, Population& POP)
//////////////////////////////////
// Force of infection on mosquitoes
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
theta.lam_M[g] = 0.0;
theta.lam_M[v] = 0.0;
}
for (int n = 0; n < POP.N_pop; n++)
{
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
theta.lam_M[g] = theta.lam_M[g] + POP.lam_n[n][g] * (theta.c_PCR*POP.people[n].I_PCR + theta.c_LM*POP.people[n].I_LM +
theta.lam_M[v] = theta.lam_M[v] + POP.lam_n[n][v] * (theta.c_PCR*POP.people[n].I_PCR + theta.c_LM*POP.people[n].I_LM +
theta.c_D*POP.people[n].I_D + theta.c_T*POP.people[n].T);
}
}
......@@ -169,18 +169,18 @@ void mosquito_step(double t, Params& theta, Population& POP)
{
mosq_rk4(t, t_step_mosq, yM, theta, POP);
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
theta.lam_S_M_track[g].push_back(theta.lam_M[g] * yM[g][3]);
theta.lam_S_M_track[g].erase(theta.lam_S_M_track[g].begin());
theta.lam_S_M_track[v].push_back(theta.lam_M[v] * yM[v][3]);
theta.lam_S_M_track[v].erase(theta.lam_S_M_track[v].begin());
}
}
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
for (int k = 0; k < N_M_comp; k++)
{
POP.yM[g][k] = yM[g][k];
POP.yM[v][k] = yM[v][k];
}
}
}
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
......@@ -40,7 +40,12 @@ public:
////////////////////////////////////////////////////
// Initialise a population of individuals at equilibrium
void equi_pop_setup(Params& theta);
void pop_setup(Params& theta);
void pop_at_equil(Params& theta);
void ind_at_equil(Params& theta);
////////////////////////////////////////////////////
// Update human individuals
......@@ -77,23 +82,23 @@ public:
//
// Depends dynamically on vector control
double yM[N_spec][N_M_comp]; // mosquito state
double yM[N_mosq][N_M_comp]; // mosquito state
double SUM_pi_w[N_spec];
double SUM_pi_z[N_spec];
double SUM_pi_w[N_mosq];
double SUM_pi_z[N_mosq];
double delta_1_VC[N_spec]; // duration spent foraging
double delta_VC[N_spec]; // duration of gonotrophic cycle = duration between blood meals
double delta_1_VC[N_mosq]; // duration spent foraging
double delta_VC[N_mosq]; // duration of gonotrophic cycle = duration between blood meals
double Z_VC[N_spec]; // average probability of mosquito repeating during a single attempt
double W_VC[N_spec]; // average probability of successfully feeding on a human during a single attempt
double Q_VC[N_spec]; // human blood index
double p_1_VC[N_spec]; // probability of surviving foraging
double mu_M_VC[N_spec]; // mosquito death rate
double aa_VC[N_spec]; // mosquito biting rate
double exp_muM_tauM_VC[N_spec]; // probability of surviving sporogony
double beta_VC[N_spec]; // egg oviposition rate
double Z_VC[N_mosq]; // average probability of mosquito repeating during a single attempt
double W_VC[N_mosq]; // average probability of successfully feeding on a human during a single attempt
double Q_VC[N_mosq]; // human blood index
double p_1_VC[N_mosq]; // probability of surviving foraging
double mu_M_VC[N_mosq]; // mosquito death rate
double aa_VC[N_mosq]; // mosquito biting rate
double exp_muM_tauM_VC[N_mosq]; // probability of surviving sporogony
double beta_VC[N_mosq]; // egg oviposition rate
////////////////////////////////////////////////////////////////
......@@ -101,15 +106,16 @@ public:
int yH[N_H_comp]; // Human compartmental states
int prev_all[12]; // Contains {N_pop, PvPR_PCR, PvPR_LM, Pv_clin, PvHR, PvHR_batches, new_PCR, new_LM, new_D, new_ACT, new_PQ, new_TQ}
int prev_U5[12]; // Contains {N_pop, PvPR_PCR, PvPR_LM, Pv_clin, PvHR, PvHR_batches, new_PCR, new_LM, new_D, new_ACT, new_PQ, new_TQ}
int prev_U10[12]; // Contains {N_pop, PvPR_PCR, PvPR_LM, Pv_clin, PvHR, PvHR_batches, new_PCR, new_LM, new_D, new_ACT, new_PQ, new_TQ}
int prev_all[15]; // Contains {N_pop, PvPR_PCR, PvPR_LM, Pv_clin, PvHR, PvHR_batches, new_PCR, new_LM, new_D, new_CQ, new_PQ, new_TQ, G6PD_test, PQ_effective, TQ_effective}
int prev_U5[15]; // Contains {N_pop, PvPR_PCR, PvPR_LM, Pv_clin, PvHR, PvHR_batches, new_PCR, new_LM, new_D, new_CQ, new_PQ, new_TQ, G6PD_test, PQ_effective, TQ_effective}
int prev_U10[15]; // Contains {N_pop, PvPR_PCR, PvPR_LM, Pv_clin, PvHR, PvHR_batches, new_PCR, new_LM, new_D, new_CQ, new_PQ, new_TQ, G6PD_test, PQ_effective, TQ_effective}
double EIR_t; // EIR
double EIR_dom_t; // EIR (domestic)
double EIR_occ_t; // EIR (occupational)
int LLIN_cov_t; // LLIN coverage
int IRS_cov_t; // IRS coverage
int IVM_cov_t; // IVM coverage
int BS_treat_t; // Indicator for chloroquine (or other blood-stage drug) treatment
int CQ_treat_t; // Indicator for chloroquine (or other blood-stage drug) treatment
int PQ_treat_t; // Indicator for primaquine treatment
int TQ_treat_t; // Indicator for tafenoquine treatment
int pregnant_t; // Coverage with front-line treatment (primaquine or tafenoquine)
......@@ -122,6 +128,12 @@ public:
int TQ_overtreat_t; // Over-treatment with 8-aminoquinolines (TQ) - defined as someone without hypnozoites being treated
int TQ_overtreat_9m_t; // Over-treatment with 8-aminoquinolines (TQ) - defined as someone without BS infection in last 9 mths being treated
int cases_M_O16_t; // Detected cases in males over 16
int cases_M_U16_t; // Detected cases in males under 16
int cases_F_O16_t; // Detected cases in females over 16
int cases_F_U16_t; // Detected cases in females under 16
int cases_preg_t; // Detected cases in pregnant women
double A_par_mean_t; // Average anti-parasite immunity
double A_clin_mean_t; // Average anti-clinical immunity
......@@ -147,6 +159,11 @@ public:
double x_age_het[N_age][N_het];
double w_age_het[N_age][N_het];
double w_gen_het_at_birth[N_gen][N_het];
double x_gen_age_het[N_mosq][N_gen][N_age][N_het];
double w_gen_age_het[N_mosq][N_gen][N_age][N_het];
////////////////////////////////////////
// Ageing rates
......@@ -156,10 +173,48 @@ public:
double P_age_bite; // Proportional change for age-dependent biting
////////////////////////////////////////
// Equilibrium vectors
vector<vector<vector<vector<vector<double>>>>> yH_eq;
vector<vector<vector<double>>> lam_eq;
vector<vector<vector<vector<double>>>> A_par_eq;
vector<vector<vector<double>>> A_par_eq_mean;
vector<vector<vector<vector<double>>>> A_clin_eq;
vector<vector<vector<double>>> A_clin_eq_mean;
vector<vector<vector<vector<double>>>> phi_LM_eq;
vector<vector<vector<vector<double>>>> phi_D_eq;
vector<vector<vector<vector<double>>>> r_PCR_eq;
vector<vector<vector<vector<double>>>> yH_eq_cumsum;
double denom_w_occup;
////////////////////////////////////////
// Equilibrium moves to zero HPZ due to PQ
int equil_setup_count;
vector<vector<vector<vector<vector<double>>>>> theta_HPZzero;
vector<vector<vector<vector<double>>>> theta_HPZzero_weight;
////////////////////////////////////////
// Treatment effectiveness
vector<vector<double>> treat_PQeff;
////////////////////////////////////////
// Indicator for age group of 20 year old woman
int index_age_20;
int index_preg_age_low;
int index_preg_age_high;
int index_occup_age_low;
int index_occup_age_high;
};
#endif
......@@ -55,10 +55,10 @@ Simulation::Simulation(SimTimes times) :
yM_t.resize(N_time);
for (int i = 0; i < N_time; i++)
{
yM_t[i].resize(N_spec);
for (int g = 0; g < N_spec; g++)
yM_t[i].resize(N_mosq);
for (int v = 0; v < N_mosq; v++)
{
yM_t[i][g].resize(N_M_comp);
yM_t[i][v].resize(N_M_comp);
}
}
......@@ -66,27 +66,28 @@ Simulation::Simulation(SimTimes times) :
prev_all.resize(N_time);
for (int i = 0; i < N_time; i++)
{
prev_all[i].resize(12);
prev_all[i].resize(15);
}
prev_U5.resize(N_time);
for (int i = 0; i < N_time; i++)
{
prev_U5[i].resize(12);
prev_U5[i].resize(15);
}
prev_U10.resize(N_time);
for (int i = 0; i < N_time; i++)
{
prev_U10[i].resize(12);
prev_U10[i].resize(15);
}
EIR_t.resize(N_time);
EIR_dom_t.resize(N_time);
EIR_occ_t.resize(N_time);
LLIN_cov_t.resize(N_time);
IRS_cov_t.resize(N_time);
IVM_cov_t.resize(N_time);
BS_treat_t.resize(N_time);
CQ_treat_t.resize(N_time);
PQ_treat_t.resize(N_time);
TQ_treat_t.resize(N_time);
pregnant_t.resize(N_time);
......@@ -97,6 +98,12 @@ Simulation::Simulation(SimTimes times) :
TQ_overtreat_t.resize(N_time);
TQ_overtreat_9m_t.resize(N_time);
cases_M_O16_t.resize(N_time);
cases_M_U16_t.resize(N_time);
cases_F_O16_t.resize(N_time);
cases_F_U16_t.resize(N_time);
cases_preg_t.resize(N_time);
A_par_mean_t.resize(N_time);
A_clin_mean_t.resize(N_time);
}
......@@ -136,13 +143,13 @@ void Simulation::run(Params& theta, Population& POP, Intervention& INTVEN)
for (int k = 0; k < N_M_comp; k++)
{
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
yM_t[i][g][k] = POP.yM[g][k];
yM_t[i][v][k] = POP.yM[v][k];
}
}
for (int k = 0; k < 12; k++)
for (int k = 0; k < 15; k++)
{
prev_all[i][k] = POP.prev_all[k];
prev_U5[i][k] = POP.prev_U5[k];
......@@ -153,7 +160,7 @@ void Simulation::run(Params& theta, Population& POP, Intervention& INTVEN)
LLIN_cov_t[i] = POP.LLIN_cov_t;
IRS_cov_t[i] = POP.IRS_cov_t;
IVM_cov_t[i] = POP.IVM_cov_t;
BS_treat_t[i] = POP.BS_treat_t;
CQ_treat_t[i] = POP.CQ_treat_t;
PQ_treat_t[i] = POP.PQ_treat_t;
TQ_treat_t[i] = POP.TQ_treat_t;
pregnant_t[i] = POP.pregnant_t;
......@@ -164,11 +171,14 @@ void Simulation::run(Params& theta, Population& POP, Intervention& INTVEN)
TQ_overtreat_t[i] = POP.TQ_overtreat_t;
TQ_overtreat_9m_t[i] = POP.TQ_overtreat_9m_t;
EIR_t[i] = 0.0;
for (int g = 0; g < N_spec; g++)
{
EIR_t[i] = EIR_t[i] + POP.aa_VC[g] * POP.yM[g][5];
}
EIR_dom_t[i] = EIR_dom_t[i] + POP.aa_VC[0] * POP.yM[0][5];
EIR_occ_t[i] = EIR_occ_t[i] + POP.aa_VC[1] * POP.yM[1][5];
cases_M_O16_t[i] = POP.cases_M_O16_t;
cases_M_U16_t[i] = POP.cases_M_U16_t;
cases_F_O16_t[i] = POP.cases_F_O16_t;
cases_F_U16_t[i] = POP.cases_F_U16_t;
cases_preg_t[i] = POP.cases_preg_t;
A_par_mean_t[i] = POP.A_par_mean_t;
A_clin_mean_t[i] = POP.A_clin_mean_t;
......@@ -198,37 +208,38 @@ void Simulation::write_output(const char *output_File)
output_Stream << yH_t[i][k] << "\t";
}
for (int g = 0; g < N_spec; g++)
for (int v = 0; v < N_mosq; v++)
{
// Write only compartments S, E and I in mosquitoes
for (int k = 3; k < N_M_comp; k++)
// Write all compartments in mosquitoes
// for (int k = 0; k < N_M_comp; k++)
// for (int k = 3; k < N_M_comp; k++)
// Write all compartments in mosquitoes
for (int k = 0; k < N_M_comp; k++)
{
output_Stream << yM_t[i][g][k] << "\t";
output_Stream << yM_t[i][v][k] << "\t";
}
}
for (int k = 0; k < 12; k++)
for (int k = 0; k < 15; k++)
{
output_Stream << prev_all[i][k] << "\t";
}
// Write output for age categories U5 and U10
/*for (int k = 0; k<12; k++)
/*for (int k = 0; k<15; k++)
{
output_Stream << prev_U5[i][k] << "\t";
}
for (int k = 0; k<12; k++)
for (int k = 0; k<15; k++)
{
output_Stream << prev_U10[i][k] << "\t";
}*/
output_Stream << EIR_t[i] << "\t";
output_Stream << EIR_dom_t[i] << "\t";
output_Stream << EIR_occ_t[i] << "\t";
output_Stream << LLIN_cov_t[i] << "\t";
output_Stream << IRS_cov_t[i] << "\t";
output_Stream << IVM_cov_t[i] << "\t";
output_Stream << BS_treat_t[i] << "\t";
output_Stream << CQ_treat_t[i] << "\t";
output_Stream << PQ_treat_t[i] << "\t";
output_Stream << TQ_treat_t[i] << "\t";
......@@ -239,11 +250,17 @@ void Simulation::write_output(const char *output_File)
output_Stream << TQ_overtreat_9m_t[i] << "\t";
// Write number of pregnant women
// output_Stream << pregnant_t[i] << "\t";
output_Stream << pregnant_t[i] << "\t";
output_Stream << cases_M_O16_t[i] << "\t";
output_Stream << cases_M_U16_t[i] << "\t";
output_Stream << cases_F_O16_t[i] << "\t";
output_Stream << cases_F_U16_t[i] << "\t";
output_Stream << cases_preg_t[i] << "\t";
// Write A_par_mean_t and A_clin_mean_t
/*output_Stream << A_par_mean_t[i] << "\t";
output_Stream << A_clin_mean_t[i] << "\t";*/
output_Stream << A_par_mean_t[i] << "\t";
output_Stream << A_clin_mean_t[i] << "\t";
output_Stream << endl;
}
......
......@@ -76,20 +76,22 @@ public:
vector<vector<vector<double>>> yM_t;
vector<double> EIR_t;
vector<double> EIR_dom_t;
vector<double> EIR_occ_t;
vector<vector<int>> prev_all; // Contains {N_pop, PvPR_PCR, PvPR_LM, Pv_clin, PvHR, PvHR_batches, new_PCR, new_LM, new_D, new_T}
vector<vector<int>> prev_U5; // Contains {N_pop, PvPR_PCR, PvPR_LM, Pv_clin, PvHR, PvHR_batches, new_PCR, new_LM, new_D, new_T}
vector<vector<int>> prev_U10; // Contains {N_pop, PvPR_PCR, PvPR_LM, Pv_clin, PvHR, PvHR_batches, new_PCR, new_LM, new_D, new_T}
////////////////////////////////////////
//////////////////////////////////////////
// Tracking coverage over time
vector<int> LLIN_cov_t;
vector<int> IRS_cov_t;
vector<int> IVM_cov_t;
vector<int> BS_treat_t;
vector<int> CQ_treat_t;
vector<int> PQ_treat_t;
vector<int> TQ_treat_t;
vector<int> pregnant_t;
......@@ -100,6 +102,12 @@ public:
vector<int> TQ_overtreat_t;
vector<int> TQ_overtreat_9m_t;
vector<int> cases_M_O16_t;
vector<int> cases_M_U16_t;
vector<int> cases_F_O16_t;
vector<int> cases_F_U16_t;
vector<int> cases_preg_t;
//////////////////////////////////////////
// Tracking immunity over time
......
......@@ -7,7 +7,8 @@
/// Institut Pasteur ///
/// michael.white@pasteur.fr ///
/// ///
/// With contributions from Dr Thomas Obadia and Diggory Hardy ///
/// With contributions from Dr Thomas Obadia, Dr Narimane Nekkab ///
/// and Diggory Hardy ///
/// ///
/// Please feel free to use and modify if you wish. However, ///
/// please provide appropriate acknowledgement and get in touch ///
......@@ -88,7 +89,7 @@ int main(int argc, char** argv)
////////////////////////////////////////////
// do we have the correct command line?
if (argc != 4 + N_spec_max)
if (argc != 4 + N_mosq)
{
std::cout << "Incorrect command line.\n";
return 0;
......@@ -96,74 +97,109 @@ int main(int argc, char** argv)
const char* parameter_File = argv[1];
const char* mosquito_File[N_spec_max];
for (int g = 0; g < N_spec_max; g++)
const char* mosquito_File[N_mosq];
for (int v = 0; v < N_mosq; v++)
{
mosquito_File[g] = argv[2 + g];
mosquito_File[v] = argv[2 + v];
}
const char* coverage_File = argv[5];
const char* output_File = argv[6];
const char* coverage_File = argv[2 + N_mosq];
const char* output_File = argv[3 + N_mosq];
////////////////////////////////////////////
// //
// 1.1.2. Initialise objects //
// 1.1.3. Initialise objects //
// //
////////////////////////////////////////////
Population PNG_pop;