API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Sep 29th, 2016
134
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 105.25 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream> // used for text output for debugging, and for general runtime text
  2. #include <fstream> // used for saving and loading files
  3. #include <sstream> // used to load files
  4. #include <string> // used for text output for debugging and for saving files
  5. #include <vector> // used for vector class to store data
  6. #include <random> // used for RNG
  7. #include <math.h> // used for rounding and other math functions
  8. #include <functional> // used to bind random distribution and generator
  9. #include <limits> // used for representing NaNs
  10. #include <chrono> // used to seed RNG and time code
  11.  
  12. // non standard libraries
  13. #include "boost/date_time/gregorian/gregorian.hpp" // used for dates
  14. #include "Eigen/Eigen" // used for matrix class
  15.  
  16. using std::vector;
  17. using namespace std::chrono;
  18. using namespace Eigen;
  19. using namespace boost::gregorian;
  20.  
  21. // used for loading csv data
  22. const int num_simulations = 1000; // total number of simulations
  23.  
  24. typedef Matrix<double, num_simulations, 1> MatrixS1d;
  25. void load_params(std::ifstream &ifs, Ref<MatrixS1d> my_params)
  26. {
  27. double d;
  28. int i(0);
  29.  
  30. for (int i = 0; i != num_simulations; ++i)
  31. {
  32. ifs >> d; // load data
  33. my_params(i) = d;
  34. ifs.ignore(1, ',');
  35. }
  36. };
  37. template<class T = std::mt19937, std::size_t N = T::state_size>
  38. auto ProperlySeededRandomEngine() -> typename std::enable_if<!!N, T>::type {
  39. typename T::result_type random_data[N];
  40. std::random_device source;
  41. std::generate(std::begin(random_data), std::end(random_data), std::ref(source));
  42. std::seed_seq seeds(std::begin(random_data), std::end(random_data));
  43. T seededEngine(seeds);
  44. seededEngine.discard(700000);
  45. return seededEngine;
  46. }
  47.  
  48. int main()
  49. {
  50. std::cout << "Model initializing: main" << std::endl;
  51.  
  52. // intialize IOFormat for saving file
  53. IOFormat CSV_Fmt(FullPrecision, DontAlignCols, ",", "\n", "", "", "", ""); // comma seperated values format
  54.  
  55. auto start_time = steady_clock::now(); // start time of program in milliseconds
  56.  
  57. // unsigned seed = static_cast<int>(std::chrono::high_resolution_clock::now().time_since_epoch().count()); // seed RNG based on current time
  58. // initialize random number generator
  59. //std::mt19937 generator(seed);
  60. auto generator = ProperlySeededRandomEngine();
  61. double average_year_days = 365.242; // average number of days in a year
  62. double average_month_days = average_year_days / 12.0; // average number of days in a month
  63.  
  64. // initialize uniform distribution
  65. std::uniform_real_distribution<double> uniform(0.0, 1.0);
  66. auto sample_uniform_distribution = std::bind(uniform, generator);
  67.  
  68. MatrixS1d prevalence_high_risk; // HIV prevalence among high risk partners
  69. MatrixS1d prevalence_low_risk_mult; // multiplier for low risk partners
  70. MatrixS1d incidence_high_risk; // HIV incidence among high risk partners
  71. MatrixS1d condom_frequency_long; // condom use frequency among long term partners
  72. MatrixS1d condom_frequency_short; // condom use frequency among short term partners
  73. MatrixS1d condom_efficacy; // condom efficacy
  74. MatrixS1d hiv_prevalence_clients; // HIV prevalence among clients of sex workers
  75. MatrixS1d hiv_risk_asymptomatic; // HIV risk per act with asymptomatic partner not on ART
  76. MatrixS1d relative_phase_risk_acute; // relative HIV acquisition risk with acute partner
  77. MatrixS1d relative_phase_risk_late; // relative HIV acquisition risk with late partner
  78. MatrixS1d relative_anal_risk; // relative infectiousness on anal intercourse
  79. MatrixS1d anal_acts_fraction_clients; // proportion of sex acts that are anal intercourse for clients
  80. MatrixS1d anal_condom_frequency_clients; // probability for using condom during anal intercourse for clients
  81. MatrixS1d vaginal_condom_frequency_clients; // probability of using condom during vaginal intercourse for clients
  82.  
  83. { // calibrate model
  84. std::cout << "Loading parameters...";
  85. std::string title = "calibrated_parameters_8_31.csv";
  86. std::ifstream file(title);
  87. if (file.is_open()){
  88. load_params(file, prevalence_high_risk);
  89. load_params(file, prevalence_low_risk_mult);
  90. load_params(file, incidence_high_risk);
  91. load_params(file, condom_frequency_long);
  92. load_params(file, condom_frequency_short);
  93. load_params(file, condom_efficacy);
  94. load_params(file, hiv_prevalence_clients);
  95. load_params(file, hiv_risk_asymptomatic);
  96. load_params(file, relative_phase_risk_acute);
  97. load_params(file, relative_phase_risk_late);
  98. load_params(file, relative_anal_risk);
  99. load_params(file, anal_acts_fraction_clients);
  100. load_params(file, anal_condom_frequency_clients);
  101. load_params(file, vaginal_condom_frequency_clients);
  102. file.close();
  103. std::cout << "done" << std::endl;
  104. }
  105. }
  106.  
  107. // epidemic parameters
  108. Vector3i phase_duration(5, 108, 24); // duration of HIV phases (0-acute, 1-asymptomatic, 2-late) in months
  109. int hiv_duration_months = phase_duration.sum();
  110. std::uniform_int_distribution<int> days_infected_distribution(0, std::lrint(average_month_days *
  111. hiv_duration_months)); // number of days before current date when partner is infected
  112. auto sample_days_infected_distribution = std::bind(days_infected_distribution, generator);
  113. double mixing_degree = 0.56; // degree of assortative mixing
  114. double fraction_high_risk_partners = 0.35; // probability to mix with high risk partners
  115. double art_efficacy = 0.9; // efficacy of ART in reducing HIV infectivity
  116. int art_progression = 3; // effect of ART on phase duration
  117. Vector4d art_drop_rate(0, 0, 0, 0); // rate of ART drop by phase (0-acute, 1-asymptomatic, 2-late, 3-extended)
  118. Vector4d art_initiation_rate(0, .1, .5, .5); // rate of ART initiation by phase (0-acute, 1-asymptomatic, 2-late, 3-extended)
  119. art_initiation_rate /= average_year_days; // convert to daily rate
  120. double fraction_high_risk_particpants = 0.25; // proportion of non FSW which are high risk participants
  121. double fraction_low_risk_participants = 1 - fraction_high_risk_particpants; // proportion of non FSW which are low risk participants
  122. Vector4d fraction_fsw_participants(0.0, 0.2, 0.5, 1); // different fractions of FSW by scenario
  123. Vector4d fraction_non_fsw(1 - fraction_fsw_participants(0),
  124. 1 - fraction_fsw_participants(1),
  125. 1 - fraction_fsw_participants(2),
  126. 1 - fraction_fsw_participants(3));
  127. // participants by scenario
  128. vector< std::discrete_distribution<int> > participant_risk_distribution{
  129. { fraction_non_fsw(0)*fraction_low_risk_participants, fraction_non_fsw(0)*fraction_high_risk_particpants, fraction_fsw_participants(0) }, // low, high, fsw
  130. { fraction_non_fsw(1)*fraction_low_risk_participants, fraction_non_fsw(1)*fraction_high_risk_particpants, fraction_fsw_participants(1) },
  131. { fraction_non_fsw(2)*fraction_low_risk_participants, fraction_non_fsw(2)*fraction_high_risk_particpants, fraction_fsw_participants(2) },
  132. { fraction_non_fsw(3)*fraction_low_risk_participants, fraction_non_fsw(3)*fraction_high_risk_particpants, fraction_fsw_participants(3) } };
  133.  
  134. vector <decltype(std::bind(participant_risk_distribution[0], generator))> sample_participant_risk_distribution;
  135.  
  136. for (int index = 0; index != 4; ++index)
  137. sample_participant_risk_distribution.push_back(std::bind(participant_risk_distribution[index], generator));
  138.  
  139.  
  140. Matrix<double, 5, 3> drop_out_rates_annual; // drop out rates for participants, 0 - low risk, 1 - high risk, 2 - female sex workers
  141. drop_out_rates_annual << 0.05, 0.05, 0.00,
  142. 0.05, 0.05, 0.05,
  143. 0.05, 0.05, 0.10,
  144. 0.05, 0.05, 0.15,
  145. 0.05, 0.05, 0.20;
  146. Matrix<double, 5, 3> drop_out_rates = drop_out_rates_annual / 12;
  147. // double art_fraction = 0.25; // fraction of infected partners on ART
  148. Vector3d art_fraction_stage(0, 0.172685, 0.65); // fraction of total ART in each stage // fraction of total ART in each stage
  149. double anal_frequency = 0.4; // fraction of sex acts which are anal for partners who practice anal
  150. double fraction_partners_practice_anal = 0.2; // fraction of partners who practice anal sex
  151.  
  152. Vector3d mixing_high_risk(mixing_degree*fraction_high_risk_partners, (1 - mixing_degree) + mixing_degree*fraction_high_risk_partners,
  153. fraction_high_risk_partners); // probability to mix with high-risk partner by (0-low risk, 1-high risk, 2-sex worker)
  154. // participants
  155. months short_to_long_partnership(9); // time to turn from a short to a long term partnership
  156. // months min_duration_partnership(1); // minimum duration of partnership
  157. months min_duration_partnership(0); // minimum duration of partnership
  158. Vector3d partner_acquisition_rate_nopartners(1.0 / (9 * average_month_days), 1.0 / (3 * average_month_days),
  159. 1.0 / (9 * average_month_days)); // daily partner acquisition rate for participants with 0 partners (0-low risk, 1-high risk,
  160. // 2-sex worker)
  161. Matrix <double, 3, 5> relative_partner_acquisition_rate; // relative partner acquisition rate, rows: low risk, high risk, sex worker,
  162. // columns: 0L0S, 0L1S, 0L2S, 1L0S, 1L1S (0L1S means 0 long 1 short term partnership)
  163. relative_partner_acquisition_rate << 1.00, 0.00, 0.00, 0.00, 0.00,
  164. 1.00, 0.54, 0.00, 0.17, 0.00,
  165. 1.00, 0.00, 0.00, 0.00, 0.00;
  166.  
  167. //vector< std::discrete_distribution<int> > initial_partner_distribution{ { 0.107, 0.294, 0.00, 0.599, 0.00 },
  168. //{ 0.080, 0.600, 0.12, 0.140, 0.06 },
  169. //{ 0.107, 0.294, 0.00, 0.599, 0.00 } };
  170. vector< std::discrete_distribution<int> > initial_partner_distribution{ { 0.00, 0.35, 0.00, 0.65, 0.00 },
  171. { 0.00, 0.620, 0.14, 0.160, 0.08 },
  172. { 0.00, 0.35, 0.00, 0.65, 0.00 } }; // initial partner
  173. // distribution of participants, rows: low risk, high risk, fsw, columns: 0l0s, 0l1s, 0l2s, 1l0s, 1l1s
  174. // (0l1s means 0 long 1 short term partnership)
  175. // construct vector of functions which sample initial_partner_distribution
  176. vector < decltype(std::bind(initial_partner_distribution[0], generator)) > sample_initial_partner_distribution;
  177.  
  178. for (int index = 0; index != 3; ++index)
  179. sample_initial_partner_distribution.push_back(std::bind(initial_partner_distribution[index], generator));
  180.  
  181. Vector2d partner_dissolution_rate_onepartner(1.0 / (10 * average_year_days), 1.0 / (1 * average_year_days)); // rate of partner dissolution
  182. // for participants with 1 partner (0- 1L0S, 1- 0L1S) (0L1S means 0 long 1 short term partnership)
  183. Matrix <int, 2, 5> relative_partner_dissolution_rate; // relative rate of partner dissolution for participants rows: long term,
  184. // short term, columns: 0L0S, 0L1S, 0L2S, 1L0S, 1L1S (0L1S means 0 long 1 short term partnership)
  185. relative_partner_dissolution_rate << 1, 1, 1, 1, 4,
  186. 2, 1, 2, 2, 2;
  187. Vector3d sex_frequency(0.2, 0.2, 0.2); // daily probability of sex act for long-term, short-term, and sex worker partners
  188. Vector3d extended_infection_mult(1 / 6.0, 1 / 6.0, 1 / 6.0); // probability to have sex with partner who is past late phase is ~ 1/30;
  189. // trial parameters
  190. const int num_participants = 2000; // number of participants in each arm
  191. date start_date(2012, Jan, 1); // start date of the trial
  192. int num_enrollment_interval = 1; // length of enrollment interval
  193. date end_enrollment_date = start_date ; // end of enrollment date
  194. int num_months_trial = 12; // length of trial in months
  195. int num_total_months = num_enrollment_interval + num_months_trial;
  196. date max_end_date = end_enrollment_date + months(num_months_trial); // end date of the trial
  197. int max_length_of_trial_days = (max_end_date - start_date).days(); // length of trial in days
  198.  
  199. int target_infection_number = 120; // target infection number for trial to end
  200.  
  201. // PrEP scenarios
  202. Matrix < double, 5, Dynamic > prep_efficacy(5, max_length_of_trial_days);
  203. prep_efficacy.row(0) = VectorXd::Constant(max_length_of_trial_days, 0.6).transpose(); // 60% fixed efficacy
  204. prep_efficacy.row(1) = VectorXd::Constant(max_length_of_trial_days, 0.8).transpose(); // 80% fixed efficacy
  205. prep_efficacy.row(2).block(0, 0, 1, 60) = VectorXd::LinSpaced(60, 1, 0.6).transpose(); // 100%-60% waning efficacy over 2 months (152 days ~ 5 average months)
  206.  
  207. double diff = prep_efficacy(2, 0) - prep_efficacy(3, 1);
  208. for (int j = 60; j != 70; ++j)
  209. {
  210. double next_val = prep_efficacy(2, j - 1) - diff;
  211. if (next_val > 0)
  212. prep_efficacy(2, j) = next_val;
  213. else
  214. prep_efficacy(2, j) = 0;
  215. }
  216.  
  217. prep_efficacy.row(3).block(0, 0, 1, 30) = VectorXd::Constant(30, 1).transpose(); // 100% over 30 days
  218. prep_efficacy.row(3).block(0, 30, 1, 30) = VectorXd::LinSpaced(30, 1, 0.2).transpose(); // 100%-20% over 30 days
  219. diff = prep_efficacy(3, 30) - prep_efficacy(3, 31);
  220. for (int j = 60; j != 70; ++j)
  221. {
  222. double next_val = prep_efficacy(3,j-1)-diff;
  223. if (next_val > 0)
  224. prep_efficacy(3, j) = next_val;
  225. else
  226. prep_efficacy(3, j) = 0;
  227. }
  228. prep_efficacy.row(4) = VectorXd::Constant(max_length_of_trial_days, 1.0).transpose(); // 100% fixed efficacy
  229. // behavioral parameters
  230. const int max_partners = 2; // maximum number of partners a participant can have
  231. double incidence_low_risk_mult = 0.75; // incidence multiplier for low risk partners
  232.  
  233.  
  234. std::uniform_int_distribution<int> short_days_distribution(30, 213); // number of days to subtract when creating a new short
  235. // partnership
  236. auto sample_short_days_distribution = std::bind(short_days_distribution, generator);
  237. // sample parameters
  238.  
  239. std::cout << "Initialization complete!" << std::endl;
  240. std::cout << "Model start!" << std::endl;
  241.  
  242. // typedef for ease of use
  243. typedef Matrix<int, num_participants, max_partners> MatrixPPi;
  244. typedef Matrix<double, num_participants, max_partners> MatrixPPd;
  245. typedef Matrix<date, num_participants, max_partners> MatrixPPdate;
  246. typedef Matrix<date, num_participants, 1> MatrixP1date;
  247. typedef Matrix <int, num_participants, 1> MatrixP1i;
  248. typedef Matrix <int, num_participants, 2> MatrixP2i;
  249. typedef Matrix <double, num_participants, 1> MatrixP1d;
  250.  
  251. for (int prep_scenario = 4; prep_scenario != 5; ++prep_scenario){
  252. for (int drop_out_scenario = 1; drop_out_scenario != 2; ++drop_out_scenario){ // drop out rate changes
  253. for (int enrollment_scenario = 0; enrollment_scenario != 4; ++enrollment_scenario){ // proportion of FSW enrolled changes
  254. // saved statistics
  255. Matrix <double, num_simulations, 10> sim_observed_incidence; // save incidence statistic
  256. Matrix <double, num_simulations, 1> sim_observed_efficacy; // save observed efficacy statistic
  257. Matrix <int, num_simulations, 1> sim_trial_duration_days; // save duration of each trial in days statistic
  258. MatrixXi sim_trial_duration_months(num_simulations, 1); // save duration of each trial in days statistic
  259. Matrix <double, num_simulations, 5> sim_exposed_fraction; // save exposed fraction statistic
  260. MatrixXi sim_num_infections(num_simulations, 8); // save number of infections statistic
  261. MatrixXi sim_num_infections_sex_work(num_simulations, 2); // save infections from sex work statistic
  262. MatrixXi sim_follow_up_days(num_simulations, 8); // save follow up time statistic
  263. MatrixXi sim_num_exposed(num_simulations, 8); // save number of exposed statistic
  264. MatrixXi sim_num_exposures(num_simulations, 8); // save number of exposures statistic
  265. MatrixXi sim_num_exposures_sex_work(num_simulations, 2); // save number of exposures statistic
  266. MatrixXi sim_cum_drop_out_gen(num_simulations, max_length_of_trial_days); // save cumulative drop out by low and high risk participants over the course of the trial
  267. MatrixXi sim_cum_drop_out_fsw(num_simulations, max_length_of_trial_days); // save cumulative drop out sex workers over the course of the trial
  268. Matrix <date, num_simulations, 2> sim_date_of_infection; // save date of infection for each arm
  269. MatrixXi sim_monthly_exposed(num_simulations, num_total_months * 8);
  270. MatrixXi sim_monthly_active(num_simulations, 2*num_total_months);
  271. MatrixXi sim_monthly_num_exposures(num_simulations, num_total_months * 8);
  272.  
  273. std::poisson_distribution<> num_clients_distribution(10.0 / 7.0);
  274. auto sample_num_clients_distribution = std::bind(num_clients_distribution, generator);
  275.  
  276. for (int sim_id = 0; sim_id != num_simulations; ++sim_id) // iterate over simulations
  277. {
  278. Vector2d num_infection_less(0,0);
  279. // sampled parameters
  280. Vector2d prevalence_partners(hiv_prevalence_clients(sim_id)*prevalence_low_risk_mult(sim_id), hiv_prevalence_clients(sim_id)); // HIV
  281. // prevalence of partners (0-low risk, 1-high risk)
  282. Vector2d incidence_partners(incidence_high_risk(sim_id)*incidence_low_risk_mult, incidence_high_risk(sim_id)); // HIV incidence of
  283. // partners (0-low risk, 1-high risk, 2-sex worker)
  284.  
  285. Vector2d condom_frequency(condom_frequency_long(sim_id), condom_frequency_short(sim_id)); // probability use condom (1-long,
  286. // 2-short, 3-casual)
  287. Vector4d relative_phase_risk(relative_phase_risk_acute(sim_id), 1.0, relative_phase_risk_late(sim_id), relative_phase_risk_late(sim_id)); // relative HIV
  288. // risk by HIV phase of partner (0-acute, 1-asymptomatic, 2-late, 3-extended late)
  289. date end_date = max_end_date; // end_date is max end_date by default
  290.  
  291. // Trial vector initialization
  292.  
  293. Matrix <int, 4, 2> num_infections = Matrix <int, 4, 2>::Zero(); // contains the number of infections, rows: 0-low risk,
  294. // 1-high risk, 2-sex worker, 3-total, columns: 0-CONTROL arm, 1-ACTIVE arm
  295. Vector2i num_infections_sex_work = Vector2i::Zero(); // number of infections from sex work 0-CONTROL arm, 1-ACTIVE arm
  296. Matrix <int, 4, 2> follow_up_days = Matrix <int, 4, 2>::Zero(); // contains the total follow up time, rows: 0-low risk,
  297. // 1-high risk, 2-sex worker, 3-total, columns: 0-CONTROL arm, 1-ACTIVE arm
  298. Matrix <int, 4, 2> num_exposed = Matrix <int, 4, 2>::Zero(); // contains the number of participants exposed, rows: 0-low risk,
  299. // 1-high risk, 2-sex worker, 3-total, columns: 0-CONTROL arm, 1-ACTIVE arm
  300. Matrix <int, 4, 2> num_exposures = Matrix <int, 4, 2>::Zero(); // contains the number of exposures, rows: 0-low risk,
  301. // 1-high risk, 2-sex worker, 3-total, columns: 0-CONTROL arm, 1-ACTIVE arm
  302. Vector2i num_exposures_sex_work = Vector2i::Zero(); // 0-control 1-active
  303. MatrixXi participant_monthly_hiv_exposure = MatrixXi::Zero(num_participants, 2);// 0 if not exposed to HIV, 1 if exposed to HIV for ACTIVE and CONTROL arms
  304. MatrixXi participant_monthly_num_exposures = MatrixXi::Zero(num_participants, 2);
  305. MatrixXi monthly_exposures = MatrixXi::Zero(8, num_total_months);
  306. Matrix <int, 8, Dynamic> monthly_exposed = Matrix <int, 8, Dynamic>::Zero(8, num_total_months);
  307. int trial_duration_months = 0;
  308. VectorXi cum_drop_out_gen = VectorXi::Zero(max_length_of_trial_days); // cumulative drop out by low and high risk participants over the course of the trial
  309. VectorXi cum_drop_out_fsw = VectorXi::Zero(max_length_of_trial_days); // cumulative drop out sex workers over the course of the trial
  310. VectorXi total_do = VectorXi::Zero(2); // total drop out 0-low and high risk, 1-sex workers
  311. Vector4i num_participants_risk_group = Vector4i::Zero(); // number of participants in each risk group,
  312. // 0-low risk, 1-high risk, 2-sex worker, 3-total
  313. Vector4i num_enrolled_risk_group = Vector4i::Zero(); // number of participants enrolled in each risk group, 0-low risk,
  314. // 1-high risk, 2-sex worker, 3-total
  315. num_participants_risk_group(3) = num_participants; // total number of participants in all risk groups
  316. // participant vector intialization
  317. MatrixP1date participant_enrollment_date; // enrollment date of participant
  318. MatrixP1date participant_max_follow_up_date; // last day of follow up for participant
  319. Matrix <date, num_participants, 2> participant_date_of_infection;
  320. MatrixP2i participant_hiv_status = MatrixP2i::Zero(); // sets participants HIV status to false for ACTIVE and CONTROL arms
  321. MatrixP2i partcipant_hiv_exposure = MatrixP2i::Zero(); // 0 if not exposed to HIV, 1 if exposed to HIV for ACTIVE and CONTROL arms
  322. MatrixP2i participant_active = MatrixP2i::Zero(); // sets participants to deactive; they become active as enrollment proceeds
  323. MatrixXi participant_monthly_active = MatrixXi::Zero(num_total_months, 2); // number of participants active each month
  324. MatrixP1i participant_prep_days_count = MatrixP1i::Zero(); // number of days since last PrEP dose.
  325. MatrixP1i participant_prep_month_count = MatrixP1i::Zero(); // number of months since last PrEP dose.
  326. int prep_interval = 2; // number of months between PrEP doses
  327. // they only become active in CONTROL xor in ACTIVE
  328. MatrixP1i participant_arm; // the arm particpant is in, 0-CONTROl, 1-ACTIVE
  329. // partner vector initialization
  330. MatrixPPi partner_active = MatrixPPi::Zero(); // active partners initialize to false
  331. MatrixP1i participant_days_enrolled = MatrixP1i::Zero(); // number of days participants have been enrolled
  332. // initiate participants
  333.  
  334. // partner vector initialization
  335. MatrixPPd partner_anal_probability = MatrixPPd::Zero(); // probability for partner to have anal sex
  336. MatrixPPi partner_risk_group; // risk group of partners for each participant
  337. MatrixPPi partner_type; // type of partnership 0-long, 1-short
  338. MatrixPPdate partner_start_date; // date of partnership start
  339. MatrixPPd partner_sex_probability; // probability for sex with partner
  340. MatrixPPd partner_condom_probability; // probability to use condom
  341. MatrixPPd partner_dissolution_probability; // probability to break up
  342. MatrixPPi partner_hiv_status; // HIV status of partner 0-negative 1-positive
  343. MatrixPPdate partner_date_of_infection; // date of infection of partner
  344. MatrixPPi partner_art_status; // ART status of partner 0-no ART 1- on ART
  345. MatrixP1i participant_risk_group; // risk group of participant
  346. MatrixP1i number_long_partners; // number of long term partners
  347. MatrixP1i number_short_partners; // number of short term partners
  348. MatrixP1i number_partners; // total number of partners
  349.  
  350. try // create enrollment blocks with equal number of risk groups and enroll participants accordingly
  351. {
  352. if ((num_participants % (2 * num_enrollment_interval)) != 0)
  353. {
  354. throw std::invalid_argument("Enrollment intervals can not be evenly split by control/active arm");
  355. }
  356.  
  357. int enrollment_size = num_participants / (2 * num_enrollment_interval);
  358.  
  359. int participant_id = 0;
  360. for (int enrollment_month = 0; enrollment_month != num_enrollment_interval; ++enrollment_month)
  361. {
  362. date enrollment_date = start_date + months(enrollment_month); // current enrollment date
  363. date max_follow_date = enrollment_date + months(num_months_trial); // maximum follow up date of participant
  364. for (int i = 0; i != enrollment_size; ++i) // add participants to CONTROL arm
  365. {
  366. participant_arm(participant_id) = 0; // set to CONTROL arm
  367. participant_risk_group(participant_id) = sample_participant_risk_distribution[enrollment_scenario]();
  368. participant_enrollment_date(participant_id) = enrollment_date; // set enrollment date
  369. participant_max_follow_up_date(participant_id) = max_follow_date; // set maximum follow up date
  370. ++num_participants_risk_group(participant_risk_group(participant_id)); // add one to risk group total
  371.  
  372. ++participant_id;
  373. }
  374. for (int i = 0; i != enrollment_size; ++i) // add participants to ACTIVE arm
  375. {
  376. participant_arm(participant_id) = 1; // set to ACTIVE arm
  377. participant_risk_group(participant_id) = sample_participant_risk_distribution[enrollment_scenario](); // set risk group
  378. participant_enrollment_date(participant_id) = enrollment_date; // set enrollment date
  379. participant_max_follow_up_date(participant_id) = max_follow_date; // set maximum follow up date
  380. ++num_participants_risk_group(participant_risk_group(participant_id)); // add one to risk group total
  381.  
  382. ++participant_id;
  383. }
  384. }
  385. }
  386. catch (const std::invalid_argument& e)
  387. {
  388. std::cerr << e.what() << std::endl; // output error message
  389. }
  390.  
  391.  
  392.  
  393.  
  394.  
  395. for (int participant_id = 0; participant_id != num_participants; ++participant_id)
  396. {
  397. // initialize partners
  398.  
  399. int initial_partners = sample_initial_partner_distribution[participant_risk_group(participant_id)](); // sample initial
  400. // partner distribution based on risk group
  401.  
  402. switch (initial_partners) // determine number of partners
  403. {
  404. case 0: number_long_partners(participant_id) = 0; // 0L0S
  405. number_short_partners(participant_id) = 0;
  406. break;
  407. case 1: number_long_partners(participant_id) = 0; // 0L1S
  408. number_short_partners(participant_id) = 1;
  409. break;
  410. case 2: number_long_partners(participant_id) = 0; // 0L2S
  411. number_short_partners(participant_id) = 2;
  412. break;
  413. case 3: number_long_partners(participant_id) = 1; // 1L0S
  414. number_short_partners(participant_id) = 0;
  415. break;
  416. case 4: number_long_partners(participant_id) = 1; // 1L1S
  417. number_short_partners(participant_id) = 1;
  418. }
  419.  
  420. // total number of partners
  421. number_partners(participant_id) = number_long_partners(participant_id) + number_short_partners(participant_id);
  422.  
  423. int current_partner = 0; // current partner
  424. for (int partner_id = 0; partner_id != number_long_partners(participant_id); ++partner_id) // add long partners
  425. {
  426. partner_risk_group(participant_id, partner_id) = (sample_uniform_distribution() <
  427. mixing_high_risk(participant_risk_group(participant_id))); // determine risk group of partner, 0-low, 1-high
  428. partner_type(participant_id, partner_id) = 0; // partner type is long
  429. partner_start_date(participant_id, partner_id) = start_date - years(1); // start one year before start date of
  430. // trial define partner sex probability as the sex frequency for long partners (0-long)
  431. partner_sex_probability(participant_id, partner_id) = sex_frequency(0) * (0.5 + sample_uniform_distribution());
  432. partner_condom_probability(participant_id, partner_id) = condom_frequency(0) *
  433. (0.5 + sample_uniform_distribution()); // define partner condom probability as the condom frequency for long
  434. // partners (0-long)
  435. partner_dissolution_probability(participant_id, partner_id) = partner_dissolution_rate_onepartner(0); // define
  436. // partner dissolution probability to be baseline probability for long partners (0-long)
  437. if (sample_uniform_distribution() < fraction_partners_practice_anal) // if partner practices anal sex
  438. partner_anal_probability(participant_id, partner_id) = anal_frequency *(0.5 + sample_uniform_distribution());
  439.  
  440. partner_hiv_status(participant_id, partner_id) = (sample_uniform_distribution() <
  441. prevalence_partners(partner_risk_group(participant_id, partner_id))); // determine if partner is HIV positive
  442. if (partner_hiv_status(participant_id, partner_id)) // if partner is HIV positive
  443. {
  444. // partner date of infection is a random number of days before start date
  445. partner_date_of_infection(participant_id, partner_id) = start_date - days(sample_days_infected_distribution());
  446.  
  447. int partner_phase;
  448. if (start_date < (partner_date_of_infection(participant_id, partner_id) +
  449. months(phase_duration(0)))) // if partner is in acute phase
  450. partner_phase = 0;
  451. else if (start_date < (partner_date_of_infection(participant_id, partner_id) + months(phase_duration(0)) +
  452. months(phase_duration(1)))) // if partner is in asymptomatic phase
  453. partner_phase = 1;
  454. else // if partner is in late phase
  455. partner_phase = 2;
  456.  
  457. partner_art_status(participant_id, partner_id) = (sample_uniform_distribution() < art_fraction_stage(partner_phase));
  458.  
  459. }
  460. else
  461. partner_date_of_infection(participant_id, partner_id) = (max_end_date + months(1)); // if partner is not HIV
  462. // positive, infection date is AFTER trial is over by default
  463.  
  464. partner_active(participant_id, partner_id) = 1; // partner becomes active
  465. ++current_partner; // move to next partner
  466. }
  467.  
  468. for (int partner_id = current_partner; partner_id != number_partners(participant_id); ++partner_id) // add short
  469. { // partners
  470. partner_risk_group(participant_id, partner_id) = (sample_uniform_distribution() <
  471. mixing_high_risk(participant_risk_group(participant_id))); // determine risk group of partner, 0-low, 1-high
  472. partner_type(participant_id, partner_id) = 1; // partner type is short
  473. partner_start_date(participant_id, partner_id) = start_date - days(sample_short_days_distribution()); // start
  474. // random days before start of trial
  475. // define partner sex probability as the sex frequency for long partners (1-short)
  476. partner_sex_probability(participant_id, partner_id) = sex_frequency(1)* (0.5 + sample_uniform_distribution());
  477. partner_condom_probability(participant_id, partner_id) = condom_frequency(1) *
  478. (0.5 + sample_uniform_distribution()); // define partner condom probability as the condom frequency for long
  479. // partners (1-short)
  480. partner_dissolution_probability(participant_id, partner_id) = partner_dissolution_rate_onepartner(1); // define
  481. // partner dissolution probability to be baseline probability for long partners (1-short)
  482. if (sample_uniform_distribution() < fraction_partners_practice_anal) // if partner practices anal sex
  483. partner_anal_probability(participant_id, partner_id) = anal_frequency *(0.5 + sample_uniform_distribution());
  484.  
  485. partner_hiv_status(participant_id, partner_id) = (sample_uniform_distribution() <
  486. prevalence_partners(partner_risk_group(participant_id, partner_id))); // determine if partner is HIV positive
  487. if (partner_hiv_status(participant_id, partner_id)) // if partner is HIV positive
  488. {
  489. // partner date of infection is a random number of days before start date
  490. partner_date_of_infection(participant_id, partner_id) = start_date - days(sample_days_infected_distribution());
  491.  
  492. int partner_phase;
  493. if (start_date < (partner_date_of_infection(participant_id, partner_id) +
  494. months(phase_duration(0)))) // if partner is in acute phase
  495. partner_phase = 0;
  496. else if (start_date < (partner_date_of_infection(participant_id, partner_id) + months(phase_duration(0)) +
  497. months(phase_duration(1)))) // if partner is in asymptomatic phase
  498. partner_phase = 1;
  499. else // if partner is in late phase
  500. partner_phase = 2;
  501.  
  502. partner_art_status(participant_id, partner_id) = (sample_uniform_distribution() < art_fraction_stage(partner_phase));
  503. }
  504. else
  505. partner_date_of_infection(participant_id, partner_id) = max_end_date + months(1); // if partner is not HIV
  506. // positive, infection date is AFTER trial is over by default
  507.  
  508. partner_active(participant_id, partner_id) = 1; // partner becomes active
  509. }
  510. }// end participant initialization
  511.  
  512. int number_days = 0; // number of days trial has gone on
  513. int number_month = 0;
  514. // start trial
  515. for (month_iterator month_itr(start_date); month_itr != max_end_date; ++month_itr) // iterate over months in the trial
  516. {
  517. month_iterator temp_month_itr = month_itr;
  518. ++temp_month_itr;
  519. date next_month = *temp_month_itr; // start of next month
  520.  
  521. // check enrollment
  522. for (int participant_id = 0; participant_id != num_participants; ++participant_id)
  523. {
  524. if (participant_enrollment_date(participant_id) == *month_itr) // if enrollment date is current date
  525. {
  526. ++num_enrolled_risk_group(participant_risk_group(participant_id));
  527. ++num_enrolled_risk_group(3);
  528. // participant becomes active
  529. if (participant_arm(participant_id) == 0)
  530. participant_active(participant_id, 0) = 1; // participant is in CONTROL arm
  531. else
  532. participant_active(participant_id, 1) = 1; // participant is in ACTIVE arm
  533. }
  534. if (participant_active(participant_id, 0))
  535. ++participant_monthly_active(number_month, 0);
  536. if (participant_active(participant_id, 1))
  537. ++participant_monthly_active(number_month, 1);
  538. }
  539. for (day_iterator day_itr(*month_itr, 1); day_itr != next_month; ++day_itr) // iterate over days in each month
  540. {
  541. date current_date = *day_itr; // current date
  542. // iterate through participants
  543. std::cout << current_date << " : " << num_infection_less(1) << " / " << num_infection_less(0) << std::endl;
  544. for (int participant_id = 0; participant_id != num_participants; ++participant_id)
  545. {
  546.  
  547. if (participant_active(participant_id, 0) || participant_active(participant_id, 1)) // if either CONTROL arm
  548. { // or ACTIVE arm is active
  549. // search for partners
  550. int partners = 3 * number_long_partners(participant_id) + number_short_partners(participant_id); // 0-0L0S,
  551. // 1-0L1S, 2-0L2S, 3-1L0S, 4-1L1S
  552.  
  553. double search_probability = partner_acquisition_rate_nopartners(participant_risk_group(participant_id)) *
  554. relative_partner_acquisition_rate(participant_risk_group(participant_id), partners); // probability to
  555. // find new partner
  556. int new_partner = -1; // -1 means NO new partner, if there is a new partner, this will store the id. New
  557. // partners always have sex with participant
  558. if ((sample_uniform_distribution() < search_probability))// randomly find new partner
  559. {
  560. if (partner_active(participant_id, 0)) // if first partner is active
  561. new_partner = 1; // new partner is second partner
  562. else
  563. new_partner = 0; // otherwise new partner is first partner
  564.  
  565. number_partners(participant_id) += 1; // add one to the total number of partners
  566. number_short_partners(participant_id) += 1; // update number of short partners (all partners start a
  567. // short term)
  568. partner_risk_group(participant_id, new_partner) = (sample_uniform_distribution() <
  569. mixing_high_risk(participant_risk_group(participant_id))); // determine risk group of partner,
  570. // 0-low, 1-high
  571. partner_type(participant_id, new_partner) = 1; // partner type is short
  572. partner_start_date(participant_id, new_partner) = current_date; // start random days before current date
  573. partner_sex_probability(participant_id, new_partner) = sex_frequency(1) * (0.5 + sample_uniform_distribution()); // define partner sex probability as the sex frequency for
  574. // short partners (1-short)
  575. partner_condom_probability(participant_id, new_partner) = condom_frequency(1) *
  576. (0.5 + sample_uniform_distribution()); // define partner condom probability as the condom frequency
  577. // for short partners (1-short)
  578. // define partner dissolution probability to be baseline probability for short partners (1-short)
  579. partner_dissolution_probability(participant_id, new_partner) = partner_dissolution_rate_onepartner(1);
  580.  
  581. if (sample_uniform_distribution() < fraction_partners_practice_anal) // if partner practices anal sex
  582. partner_anal_probability(participant_id, new_partner) = anal_frequency *
  583. (0.5 + sample_uniform_distribution());
  584.  
  585. partner_hiv_status(participant_id, new_partner) = (sample_uniform_distribution() <
  586. prevalence_partners(partner_risk_group(participant_id, new_partner))); // determine if partner is
  587. // HIV positive
  588. if (partner_hiv_status(participant_id, new_partner)) // if partner is HIV positive
  589. {
  590. // partner date of infection is a random number of days before start date
  591. partner_date_of_infection(participant_id, new_partner) = current_date - days(sample_days_infected_distribution());
  592.  
  593. int partner_phase;
  594. if (current_date < (partner_date_of_infection(participant_id, new_partner) +
  595. months(phase_duration(0)))) // if partner is in acute phase
  596. partner_phase = 0;
  597. else if (current_date < (partner_date_of_infection(participant_id, new_partner) + months(phase_duration(0)) +
  598. months(phase_duration(1)))) // if partner is in asymptomatic phase
  599. partner_phase = 1;
  600. else // if partner is in late phase
  601. partner_phase = 2;
  602.  
  603. partner_art_status(participant_id, new_partner) = (sample_uniform_distribution() < art_fraction_stage(partner_phase));
  604. }
  605. else
  606. partner_date_of_infection(participant_id, new_partner) = max_end_date + months(1); // if partner is not
  607. // HIV positive, infection date is AFTER trial is over by default
  608.  
  609. partner_active(participant_id, new_partner) = 1; // partner becomes active
  610. } // end search partner
  611.  
  612. // have partnership sex, break up with partners and update from short to long partners
  613. for (int partner_id = 0; partner_id != 2; ++partner_id)
  614. {
  615. if (partner_active(participant_id, partner_id)) // only work with active partners
  616. {
  617. ++num_infection_less(0);
  618. // have sex with partner
  619. double temp_sex_probability = partner_sex_probability(participant_id, partner_id); // temporary
  620. // probability storage
  621. //if (number_partners(participant_id) > 1)
  622. //{
  623. // temp_sex_probability = temp_sex_probability / number_partners(participant_id); // reduce
  624. // // probability for sex by number of partners
  625. //}
  626. // determine HIV phase partner is in
  627. Vector3i partner_phase_duration = phase_duration; // base phase duration
  628.  
  629. int partner_phase;
  630. if (current_date < (partner_date_of_infection(participant_id, partner_id) +
  631. months(partner_phase_duration(0)))) // if partner is in acute phase
  632. partner_phase = 0;
  633. else if (current_date < (partner_date_of_infection(participant_id, partner_id) + months(partner_phase_duration(0)) +
  634. months(partner_phase_duration(1)))) // if partner is in asymptomatic phase
  635. partner_phase = 1;
  636. else if (current_date < (partner_date_of_infection(participant_id, partner_id) + months(hiv_duration_months))) // if partner is in late phase
  637. partner_phase = 2;
  638. else
  639. partner_phase = 3; // partner past late phase
  640.  
  641. if (partner_phase == 3)
  642. {
  643. temp_sex_probability *= extended_infection_mult(partner_type(participant_id, partner_id)); //reduce sex probability if after late phase
  644. }
  645.  
  646. if ((sample_uniform_distribution() < temp_sex_probability) || (partner_id == new_partner)) // if
  647. { // sample true for sex OR if there is a new partner
  648. if (partner_hiv_status(participant_id, partner_id)) // if partner is HIV positive
  649. {
  650.  
  651. int use_condom = (sample_uniform_distribution() <
  652. partner_condom_probability(participant_id, partner_id)); // determine if condom is used
  653.  
  654. double hiv_risk = hiv_risk_asymptomatic(sim_id) *
  655. (1 - use_condom * condom_efficacy(sim_id)); // risk for hiv transmission taking condom use into
  656. // account
  657.  
  658. if ((sample_uniform_distribution() < partner_anal_probability(participant_id, partner_id)))
  659. { // if sex act is anal
  660. hiv_risk = hiv_risk * relative_anal_risk(sim_id);
  661. }
  662.  
  663.  
  664. // determine HIV phase partner is in
  665. if (partner_art_status(participant_id, partner_id)) // if partner is on ART
  666. {
  667. partner_phase_duration(1) = art_progression * partner_phase_duration(1); // ART
  668. // lengthens the phases
  669. partner_phase_duration(2) = art_progression * partner_phase_duration(2);
  670. hiv_risk = (1 - art_efficacy) * hiv_risk; // ART reduces HIV transmission risk
  671. }
  672.  
  673. // HIV transmission
  674. if (participant_active(participant_id, 0) && !participant_hiv_status(participant_id, 0)) // if control arm is active and not infected
  675. {
  676. ++participant_monthly_num_exposures(participant_id, 0);
  677. if (!partcipant_hiv_exposure(participant_id, 0)) // if participant is has not been exposed
  678. {
  679. partcipant_hiv_exposure(participant_id, 0) = 1; //participant is exposed to hiv
  680. ++num_exposed(participant_risk_group(participant_id), 0);
  681. ++num_exposed(3, 0);
  682. }
  683. if (!participant_monthly_hiv_exposure(participant_id, 0)) // if participant is has not been exposed
  684. {
  685. participant_monthly_hiv_exposure(participant_id, 0) = 1; //participant is exposed to hiv
  686. ++monthly_exposed(participant_risk_group(participant_id), number_month);
  687. ++monthly_exposed(3, number_month);
  688. }
  689. ++num_exposures(participant_risk_group(participant_id), 0);
  690. ++num_exposures(3, 0);
  691. ++num_infection_less(1);
  692.  
  693.  
  694. if (sample_uniform_distribution() < (hiv_risk * relative_phase_risk(partner_phase))) // if
  695. { // HIV is transmitted
  696. participant_hiv_status(participant_id, 0) = 1; // participant becomes infected
  697. participant_date_of_infection(participant_id, 0) = current_date;
  698. }
  699. }
  700. // HIV transmission with PrEP use
  701. // if with PrEP participant is NOT infected (ACTIVE ARM)
  702. if (participant_active(participant_id, 1) && !participant_hiv_status(participant_id, 1)) // if active arm is active and not infected
  703. {
  704. ++participant_monthly_num_exposures(participant_id, 1);
  705. if (!partcipant_hiv_exposure(participant_id, 1)) // if participant is has not been exposed
  706. {
  707. partcipant_hiv_exposure(participant_id, 1) = 1; //participant is exposed to hiv
  708. ++num_exposed(participant_risk_group(participant_id), 1);
  709. ++num_exposed(3, 1);
  710. }
  711. if (!participant_monthly_hiv_exposure(participant_id, 1)) // if participant is has not been exposed
  712. {
  713. participant_monthly_hiv_exposure(participant_id, 1) = 1; //participant is exposed to hiv
  714. ++monthly_exposed(participant_risk_group(participant_id) + 4, number_month);
  715. ++monthly_exposed(7, number_month);
  716. }
  717. ++num_exposures(participant_risk_group(participant_id), 1);
  718. ++num_exposures(3, 1);
  719. // efficacy based on days enrolled
  720. double protection_efficacy =
  721. prep_efficacy(prep_scenario, participant_prep_days_count(participant_id));
  722.  
  723. if (sample_uniform_distribution() < (hiv_risk * relative_phase_risk(partner_phase) *
  724. (1 - protection_efficacy))) // if HIV is transmitted
  725. {
  726. participant_hiv_status(participant_id, 1) = 1; // participant becomes infected
  727. participant_date_of_infection(participant_id, 1) = current_date;
  728. }
  729. }
  730. }
  731. }
  732. // break up with partners
  733. if (current_date > (partner_start_date(participant_id, partner_id) + min_duration_partnership)) // if partnership is longer than minimum partnership duration
  734. {
  735. double probability_to_break = partner_dissolution_probability(participant_id, partner_id) *
  736. relative_partner_dissolution_rate(partner_type(participant_id, partner_id), partners);
  737.  
  738. // (partners is 0-0L0S, 1-0L1S, 2-0L2S, 3-1L0S, 4-1L1S)
  739. if (sample_uniform_distribution() < probability_to_break)
  740. {
  741.  
  742. if (partner_type(participant_id, partner_id) == 0) // if partner is long remove long partner
  743. number_long_partners(participant_id) -= 1;
  744. else
  745. number_short_partners(participant_id) -= 1; // if partner is short remove short partner
  746.  
  747. partner_active(participant_id, partner_id) = 0; // partner is no longer active
  748. }
  749. }
  750. // update partnership from short to long
  751. else if ((number_long_partners(participant_id) == 0) &&
  752. (partner_type(participant_id, partner_id) == 1)) // if number of long partners is 0 and partner
  753. { // is short (can only have 1 long partner)
  754. if (current_date > (partner_start_date(participant_id, partner_id) +
  755. short_to_long_partnership)) // if partnership has lasted long enough
  756. {
  757. partner_type(participant_id, partner_id) = 0; // partner type is long
  758. partner_dissolution_probability(participant_id, partner_id) =
  759. partner_dissolution_rate_onepartner(0); // update dissolution probability to long
  760. partner_sex_probability(participant_id, partner_id) = sex_frequency(0) * (0.5 + sample_uniform_distribution()); // update sex probability to long
  761. // (with randomness)
  762. partner_condom_probability(participant_id, partner_id) = condom_frequency(0) *
  763. (0.5 + sample_uniform_distribution()); // update condom use probability to long
  764. // (with randomness)
  765.  
  766. number_long_partners(participant_id) += 1; // add one more long partner
  767. number_short_partners(participant_id) -= 1; // remove short partner
  768. }
  769. }
  770.  
  771. if (partner_art_status(participant_id, partner_id))
  772. {
  773. partner_art_status(participant_id, partner_id) = (sample_uniform_distribution() > art_drop_rate(partner_phase)); // drop art
  774. }
  775. else
  776. {
  777. partner_art_status(participant_id, partner_id) = (sample_uniform_distribution() < art_initiation_rate(partner_phase)); // initiate art
  778. }
  779.  
  780. // partner acquire HIV outside of relationship
  781. if ((!partner_hiv_status(participant_id, partner_id)) && (sample_uniform_distribution() <
  782. (incidence_partners(partner_risk_group(participant_id, partner_id)) / average_year_days)))
  783. {
  784. partner_hiv_status(participant_id, partner_id) = 1; // partner gets infected
  785. partner_date_of_infection(participant_id, partner_id) = current_date; // date of infection is today
  786.  
  787. }
  788. }
  789. } // end partnership update
  790.  
  791. number_partners(participant_id) = number_long_partners(participant_id) +
  792. number_short_partners(participant_id); // update number of partners after possible break ups
  793.  
  794. // practice sex work
  795. if (participant_risk_group(participant_id) == 2) // if participant is sex worker
  796. {
  797. int num_clients = sample_num_clients_distribution();
  798. if (num_clients > 6)
  799. num_clients = 6;
  800. for (int client_id = 0; client_id != num_clients; ++client_id)
  801. {
  802.  
  803. // HIV status of partner
  804. if (sample_uniform_distribution() < hiv_prevalence_clients(sim_id)) // if partner is HIV positive
  805. {
  806. // partner date of infection is a random number of days before current date
  807. date client_partner_date_of_infection = current_date - days(sample_days_infected_distribution());
  808. int client_partner_art_status = 0;
  809. int partner_phase;
  810. Vector3i partner_phase_duration = phase_duration; // base phase duration
  811.  
  812. if (current_date < (client_partner_date_of_infection +
  813. months(partner_phase_duration(0)))) // if partner is in acute phase
  814. partner_phase = 0;
  815. else if (current_date < (client_partner_date_of_infection + months(partner_phase_duration(0)) +
  816. months(partner_phase_duration(1)))) // if partner is in asymptomatic phase
  817. partner_phase = 1;
  818. else // if partner is in late phase
  819. partner_phase = 2;
  820.  
  821. client_partner_art_status = (sample_uniform_distribution() < art_fraction_stage(partner_phase));
  822.  
  823. double hiv_risk = hiv_risk_asymptomatic(sim_id);
  824.  
  825.  
  826. if (client_partner_art_status) // if partner is on ART
  827. {
  828. partner_phase_duration(1) = art_progression * partner_phase_duration(1); // ART lengthens the
  829. // phases
  830. partner_phase_duration(2) = art_progression * partner_phase_duration(2);
  831. hiv_risk = (1 - art_efficacy) * hiv_risk; // ART reduces HIV transmission risk
  832. }
  833. int use_condom;
  834. if (sample_uniform_distribution() < anal_acts_fraction_clients(sim_id)) // if anal sex act
  835. {
  836. use_condom = (sample_uniform_distribution() < anal_condom_frequency_clients(sim_id));
  837. hiv_risk = hiv_risk * (1 - use_condom * condom_efficacy(sim_id));
  838. hiv_risk = hiv_risk * relative_anal_risk(sim_id);
  839. }
  840. else
  841. {
  842. use_condom = (sample_uniform_distribution() < vaginal_condom_frequency_clients(sim_id));
  843. hiv_risk = hiv_risk * (1 - use_condom * condom_efficacy(sim_id));
  844. }
  845.  
  846. if (current_date < (client_partner_date_of_infection +
  847. months(partner_phase_duration(0)))) // if partner is in acute phase
  848. partner_phase = 0;
  849. else if (current_date < (client_partner_date_of_infection + months(partner_phase_duration(0)) +
  850. months(partner_phase_duration(1)))) // if partner is in asymptomatic phase
  851. partner_phase = 1;
  852. else // if partner is in late phase
  853. partner_phase = 2;
  854.  
  855. // HIV transmission without PrEP use
  856. // if without PrEP participant is NOT infected (CONTROL ARM)
  857. if (participant_active(participant_id, 0) && !participant_hiv_status(participant_id, 0)) // if control arm is active and not infected
  858. {
  859. ++participant_monthly_num_exposures(participant_id, 0);
  860. if (!partcipant_hiv_exposure(participant_id, 0)) // if participant is has not been exposed
  861. {
  862. partcipant_hiv_exposure(participant_id, 0) = 1; //participant is exposed to hiv
  863. ++num_exposed(participant_risk_group(participant_id), 0);
  864. ++num_exposed(3, 0);
  865. }
  866. if (!participant_monthly_hiv_exposure(participant_id, 0)) // if participant is has not been exposed
  867. {
  868. participant_monthly_hiv_exposure(participant_id, 0) = 1; //participant is exposed to hiv
  869. ++monthly_exposed(participant_risk_group(participant_id), number_month);
  870. ++monthly_exposed(3, number_month);
  871. }
  872. ++num_exposures(participant_risk_group(participant_id), 0);
  873. ++num_exposures(3, 0);
  874. ++num_exposures_sex_work(0);
  875. if (sample_uniform_distribution() < (hiv_risk * relative_phase_risk(partner_phase))) // if
  876. { // HIV is transmitted
  877. participant_hiv_status(participant_id, 0) = 1; // participant becomes infected
  878. participant_date_of_infection(participant_id, 0) = current_date;
  879. ++num_infections_sex_work(0);
  880. }
  881. }
  882. // HIV transmission with PrEP use
  883. // if with PrEP participant is NOT infected (ACTIVE ARM)
  884. if (participant_active(participant_id, 1) && !participant_hiv_status(participant_id, 1)) // if active arm is active and not infected
  885. {
  886. ++participant_monthly_num_exposures(participant_id, 1);
  887. if (!partcipant_hiv_exposure(participant_id, 1)) // if participant is has not been exposed
  888. {
  889. partcipant_hiv_exposure(participant_id, 1) = 1; //participant is exposed to hiv
  890. ++num_exposed(participant_risk_group(participant_id), 1);
  891. ++num_exposed(3, 1);
  892. }
  893. if (!participant_monthly_hiv_exposure(participant_id, 1)) // if participant is has not been exposed
  894. {
  895. participant_monthly_hiv_exposure(participant_id, 1) = 1; //participant is exposed to hiv
  896. ++monthly_exposed(participant_risk_group(participant_id) + 4, number_month);
  897. ++monthly_exposed(7, number_month);
  898. }
  899. ++num_exposures(participant_risk_group(participant_id), 1);
  900. ++num_exposures(3, 1);
  901. ++num_exposures_sex_work(1);
  902. // efficacy based on days enrolled
  903. double protection_efficacy =
  904. prep_efficacy(prep_scenario, participant_prep_days_count(participant_id));
  905. if (sample_uniform_distribution() < (hiv_risk * relative_phase_risk(partner_phase) *
  906. (1 - protection_efficacy))) // if HIV is transmitted
  907. {
  908. participant_hiv_status(participant_id, 1) = 1; // participant becomes infected
  909. participant_date_of_infection(participant_id, 1) = current_date;
  910. ++num_infections_sex_work(1);
  911. }
  912. }
  913. }
  914. }
  915. } // end casual sex
  916. ++participant_days_enrolled(participant_id); // participant has been enrolled for another day
  917.  
  918. } // end active participant check
  919. if (participant_active(participant_id, 1))
  920. {
  921. ++participant_prep_days_count(participant_id);
  922. }
  923. }// participant loop
  924. cum_drop_out_gen(number_days) += total_do(0); // increment cumulative drop out
  925. cum_drop_out_fsw(number_days) += total_do(1); // increment cumulative drop out
  926. ++number_days;
  927. }// end days
  928.  
  929. ++trial_duration_months;
  930. // monthly check up with participants, if HIV positive, remove from trial
  931. int num_days = (next_month - *month_itr).days(); //total follow up time for participant
  932. for (int participant_id = 0; participant_id != num_participants; ++participant_id)
  933. {
  934. if (participant_active(participant_id, 0)) // check if participant is active in CONTROL arm
  935. {
  936. // check if participant drops out of trial
  937. if (sample_uniform_distribution() <
  938. drop_out_rates(drop_out_scenario, participant_risk_group(participant_id))){
  939. // participant becomes inactive regardless of which arm they are in
  940. participant_active(participant_id, 0) = 0; // participant becomes inactive in CONTROL arm
  941. if (participant_risk_group(participant_id) == 2)
  942. ++total_do(1);
  943. else
  944. ++total_do(0);
  945. }
  946.  
  947. ++monthly_exposures(participant_risk_group(participant_id), number_month);
  948. ++monthly_exposures(3, number_month);
  949.  
  950. participant_monthly_hiv_exposure(participant_id, 0) = 0;
  951. participant_monthly_num_exposures(participant_id, 0) = 0;
  952. follow_up_days(participant_risk_group(participant_id), 0) += num_days; // total number of days since start date
  953. // for participant's risk group
  954. follow_up_days(3, 0) += num_days; // total number of days since start date for all risk groups
  955.  
  956. if (participant_hiv_status(participant_id, 0)) // check if participant is HIV positive
  957. {
  958. ++num_infections(participant_risk_group(participant_id), 0); // increase infection counter for
  959. // participant's risk group
  960. ++num_infections(3, 0); // increase infection counter for total
  961. participant_active(participant_id, 0) = 0; // remove participant from trial
  962. }
  963. if (next_month >= participant_max_follow_up_date(participant_id))
  964. {
  965. participant_active(participant_id, 0) = 0; // remove participant from trial
  966. }
  967.  
  968. }
  969.  
  970. if (participant_active(participant_id, 1)) // check if participant is active in ACTIVE arm
  971. {
  972. // check if participant drops out of trial
  973. if (sample_uniform_distribution() <
  974. drop_out_rates(drop_out_scenario, participant_risk_group(participant_id))){
  975. // participant becomes inactive regardless of which arm they are in
  976. participant_active(participant_id, 1) = 0; // participant becomes inactive in ACTIVE arm
  977. if (participant_risk_group(participant_id) == 2)
  978. ++total_do(1);
  979. else
  980. ++total_do(0);
  981. }
  982.  
  983. ++monthly_exposures(participant_risk_group(participant_id) + 4, number_month);
  984. ++monthly_exposures(7, number_month);
  985.  
  986. participant_monthly_hiv_exposure(participant_id, 1) = 0;
  987. participant_monthly_num_exposures(participant_id, 1) = 0;
  988. ++participant_prep_month_count(participant_id);
  989. if (participant_prep_month_count(participant_id) >= prep_interval)
  990. {
  991. participant_prep_days_count(participant_id) = 0;
  992. participant_prep_month_count(participant_id) = 0;
  993. }
  994.  
  995.  
  996. follow_up_days(participant_risk_group(participant_id), 1) += num_days; // total number of days since start date
  997. // for participant's risk group
  998. follow_up_days(3, 1) += num_days; // total number of days since start date for all risk groups
  999. if (participant_hiv_status(participant_id, 1)) // check if participant is HIV positive
  1000. {
  1001. ++num_infections(participant_risk_group(participant_id), 1); // increase infection counter for
  1002. // participant's risk group
  1003. ++num_infections(3, 1); // increase infection counter for total
  1004. participant_active(participant_id, 1) = 0; // remove participant from trial
  1005. }
  1006.  
  1007. if (next_month >= participant_max_follow_up_date(participant_id))
  1008. {
  1009. participant_active(participant_id, 1) = 0; // remove participant from trial
  1010. }
  1011.  
  1012. }
  1013.  
  1014. }
  1015.  
  1016. if ((num_infections(3, 0) + num_infections(3, 1)) >= target_infection_number) // if infection number reaches target
  1017. {
  1018. end_date = next_month; // end date is next month
  1019. break; // exit loop
  1020. }
  1021. ++number_month;
  1022. } // end months / end of trial
  1023.  
  1024. for (int y = 0; y != num_total_months; ++y) // number exposed monthly
  1025. {
  1026. sim_monthly_active(sim_id, y) = participant_monthly_active(y, 0);
  1027. sim_monthly_active(sim_id, num_total_months+y) = participant_monthly_active(y, 1);
  1028. for (int x = 0; x != 8; ++x)
  1029. {
  1030. sim_monthly_exposed(sim_id, x + y * 8) = monthly_exposed(x, y);
  1031. }
  1032. }
  1033.  
  1034. for (int y = 0; y != num_total_months; ++y) // number exposures monthly
  1035. {
  1036. for (int x = 0; x != 8; ++x)
  1037. {
  1038. sim_monthly_num_exposures(sim_id, x + y * 8) = monthly_exposures(x, y);
  1039. }
  1040. }
  1041.  
  1042.  
  1043. for (int y = 0; y != 2; ++y) // calculate incidence
  1044. {
  1045. for (int x = 0; x != 4; ++x)
  1046. {
  1047. if (follow_up_days(x, y) != 0)
  1048. {
  1049. sim_observed_incidence(sim_id, x + 4 * y) = num_infections(x, y) /
  1050. (follow_up_days(x, y) / average_year_days); // calculate observed yearly incidence
  1051. }
  1052. else
  1053. {
  1054. //sim_observed_incidence(sim_id, x + 4 * y) = std::numeric_limits<double>::quiet_NaN(); // incidence not observed!!
  1055. sim_observed_incidence(sim_id, x + 4 * y) = 0; // incidence not observed!!
  1056. }
  1057. }
  1058. }
  1059.  
  1060.  
  1061. for (int y = 0; y != 2; ++y) // calculate incidence
  1062. {
  1063. for (int x = 0; x != 4; ++x)
  1064. {
  1065. if (follow_up_days(x, y) != 0)
  1066. {
  1067. sim_observed_incidence(sim_id, x + 4 * y) = num_infections(x, y) /
  1068. (follow_up_days(x, y) / average_year_days); // calculate observed yearly incidence
  1069. }
  1070. else
  1071. {
  1072. //sim_observed_incidence(sim_id, x + 4 * y) = std::numeric_limits<double>::quiet_NaN(); // incidence not observed!!
  1073. sim_observed_incidence(sim_id, x + 4 * y) = 0; // incidence not observed!!
  1074. }
  1075. }
  1076. }
  1077.  
  1078. if ((follow_up_days(0, 0) + follow_up_days(1, 0)) != 0){
  1079. sim_observed_incidence(sim_id, 8) = (num_infections(0, 0) + num_infections(1, 0)) / ((follow_up_days(0, 0) +
  1080. follow_up_days(1, 0)) / average_year_days); // calculate low and high risk only for CONTROL group
  1081. }
  1082. else{
  1083. // sim_observed_incidence(sim_id, 8) = std::numeric_limits<double>::quiet_NaN();
  1084. sim_observed_incidence(sim_id, 8) = 0;
  1085. }
  1086.  
  1087. if ((follow_up_days(0, 1) + follow_up_days(1, 1)) != 0){
  1088. sim_observed_incidence(sim_id, 9) = (num_infections(0, 1) + num_infections(1, 1)) / ((follow_up_days(0, 1) +
  1089. follow_up_days(1, 1)) / average_year_days); // calculate low and high risk only for ACTIVE group
  1090. }
  1091. else{
  1092. // sim_observed_incidence(sim_id, 9) = std::numeric_limits<double>::quiet_NaN();
  1093. sim_observed_incidence(sim_id, 9) = 0;
  1094. }
  1095.  
  1096.  
  1097. if (sim_observed_incidence(sim_id, 3) != 0)
  1098. {
  1099. sim_observed_efficacy(sim_id) = 1 - (sim_observed_incidence(sim_id, 7) /
  1100. sim_observed_incidence(sim_id, 3)); // calculate observed efficacy
  1101. }
  1102. else{
  1103. //sim_observed_efficacy(sim_id) = std::numeric_limits<double>::quiet_NaN(); // efficacy not observed!!
  1104. sim_observed_efficacy(sim_id) = 0; // efficacy not observed!!
  1105. }
  1106.  
  1107. sim_trial_duration_days(sim_id) = (end_date - start_date).days(); // trial duration
  1108. sim_trial_duration_months(sim_id) = trial_duration_months;
  1109. for (int x = 0; x != 4; ++x) // calculate exposed fraction
  1110. {
  1111. if (num_participants_risk_group(x) != 0)
  1112. {
  1113. sim_exposed_fraction(sim_id, x) = (double)(num_exposed(x, 0) + num_exposed(x, 1)) / (num_enrolled_risk_group(x));
  1114. } // calculate exposed fraction
  1115. else{
  1116. // sim_exposed_fraction(sim_id, x) = std::numeric_limits<double>::quiet_NaN(); // exposed fraction not observed!!
  1117. sim_exposed_fraction(sim_id, x) = 0; // exposed fraction not observed!!
  1118. }
  1119. }
  1120. if (num_enrolled_risk_group(0) + num_enrolled_risk_group(1) != 0){
  1121. sim_exposed_fraction(sim_id, 4) = (double)(num_exposed(0, 0) + num_exposed(0, 1) + num_exposed(1, 0) + num_exposed(1, 1))
  1122. / ((num_enrolled_risk_group(0) + num_enrolled_risk_group(1))); // calculate low and high risk only
  1123. }
  1124. else{
  1125. // sim_exposed_fraction(sim_id, 4) = std::numeric_limits<double>::quiet_NaN();
  1126. sim_exposed_fraction(sim_id, 4) = 0;
  1127. }
  1128.  
  1129.  
  1130. for (int y = 0; y != 2; ++y) // save number of infections
  1131. {
  1132. for (int x = 0; x != 4; ++x)
  1133. {
  1134. sim_num_infections(sim_id, x + 4 * y) = num_infections(x, y);
  1135. }
  1136. }
  1137.  
  1138. sim_num_infections_sex_work(sim_id, 0) = num_infections_sex_work(0);
  1139. sim_num_infections_sex_work(sim_id, 1) = num_infections_sex_work(1);
  1140.  
  1141. for (int y = 0; y != 2; ++y) // save follow up time
  1142. {
  1143. for (int x = 0; x != 4; ++x)
  1144. {
  1145. sim_follow_up_days(sim_id, x + 4 * y) = follow_up_days(x, y);
  1146. }
  1147. }
  1148.  
  1149. for (int y = 0; y != 2; ++y) // save number of exposed
  1150. {
  1151. for (int x = 0; x != 4; ++x)
  1152. {
  1153. sim_num_exposed(sim_id, x + 4 * y) = num_exposed(x, y);
  1154. }
  1155. }
  1156.  
  1157. for (int y = 0; y != 2; ++y) // save number of exposures
  1158. {
  1159. for (int x = 0; x != 4; ++x)
  1160. {
  1161. sim_num_exposures(sim_id, x + 4 * y) = num_exposures(x, y);
  1162. }
  1163. }
  1164. sim_num_exposures_sex_work(sim_id, 0) = num_exposures_sex_work(0);
  1165. sim_num_exposures_sex_work(sim_id, 1) = num_exposures_sex_work(1);
  1166. // save cumulative drop out
  1167. sim_cum_drop_out_gen.block(sim_id, 0, 1, max_length_of_trial_days) = cum_drop_out_gen.transpose();
  1168. sim_cum_drop_out_fsw.block(sim_id, 0, 1, max_length_of_trial_days) = cum_drop_out_fsw.transpose();
  1169.  
  1170. } // end of simulations
  1171. // write to file
  1172.  
  1173. {
  1174. std::string prep_type; // type of prep, fixed vs waning
  1175. std::string scenarioName = "";
  1176. std::string plotsLoc = "plots" + scenarioName + "/";
  1177. switch (prep_scenario)
  1178. {
  1179. case 0:
  1180. {
  1181. prep_type = "60f" + scenarioName;
  1182. break;
  1183. }
  1184. case 1:
  1185. {
  1186. prep_type = "80f" + scenarioName;
  1187. break;
  1188. }
  1189. case 2:
  1190. {
  1191. prep_type = "100_60w" + scenarioName;
  1192. break;
  1193. }
  1194. case 3:
  1195. {
  1196. prep_type = "100_100_20w" + scenarioName;
  1197. break;
  1198. }
  1199. case 4:
  1200. {
  1201. prep_type = "100f" + scenarioName;
  1202. }
  1203. }
  1204. {
  1205. // save monthly exposed
  1206. std::string title = "monthly_exposed_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1207. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1208. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1209. + "_e"
  1210. + prep_type
  1211. + ".csv";
  1212. std::cout << "Saving file: " << title << "...";
  1213. std::ofstream file(plotsLoc + title);
  1214. if (file.is_open())
  1215. {
  1216. file << sim_monthly_exposed.format(CSV_Fmt);
  1217. file.close();
  1218. std::cout << "done" << std::endl;
  1219. }
  1220. else
  1221. std::cout << "file could not be saved" << std::endl;
  1222. }
  1223. {
  1224. // save monthly exposures
  1225. std::string title = "monthly_exposures_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1226. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1227. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1228. + "_e"
  1229. + prep_type
  1230. + ".csv";
  1231. std::cout << "Saving file: " << title << "...";
  1232. std::ofstream file(plotsLoc + title);
  1233. if (file.is_open())
  1234. {
  1235. file << sim_monthly_num_exposures.format(CSV_Fmt);
  1236. file.close();
  1237. std::cout << "done" << std::endl;
  1238. }
  1239. else
  1240. std::cout << "file could not be saved" << std::endl;
  1241. }
  1242. {
  1243. // save monthly active
  1244. std::string title = "monthly_active_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1245. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1246. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1247. + "_e"
  1248. + prep_type
  1249. + ".csv";
  1250. std::cout << "Saving file: " << title << "...";
  1251. std::ofstream file(plotsLoc + title);
  1252. if (file.is_open())
  1253. {
  1254. file << sim_monthly_active.format(CSV_Fmt);
  1255. file.close();
  1256. std::cout << "done" << std::endl;
  1257. }
  1258. else
  1259. std::cout << "file could not be saved" << std::endl;
  1260. }
  1261. {
  1262. // save incidence
  1263. std::string title = "incidence_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1264. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1265. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1266. + "_e"
  1267. + prep_type
  1268. + ".csv";
  1269. std::cout << "Saving file: " << title << "...";
  1270. std::ofstream file(plotsLoc + title);
  1271. if (file.is_open())
  1272. {
  1273. file << sim_observed_incidence.format(CSV_Fmt);
  1274. file.close();
  1275. std::cout << "done" << std::endl;
  1276. }
  1277. else
  1278. std::cout << "file could not be saved" << std::endl;
  1279. }
  1280. {
  1281. // save efficacy
  1282. std::string title = "efficacy_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1283. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1284. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1285. + "_e"
  1286. + prep_type
  1287. + ".csv";
  1288. std::cout << "Saving file: " << title << "...";
  1289. std::ofstream file(plotsLoc + title);
  1290. if (file.is_open())
  1291. {
  1292. file << sim_observed_efficacy.format(CSV_Fmt);
  1293. file.close();
  1294. std::cout << "done" << std::endl;
  1295. }
  1296. else
  1297. std::cout << "file could not be saved" << std::endl;
  1298. }
  1299. {
  1300. // save trial duration
  1301. std::string title = "trial_duration_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1302. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1303. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1304. + "_e"
  1305. + prep_type
  1306. + ".csv";
  1307. std::cout << "Saving file: " << title << "...";
  1308. std::ofstream file(plotsLoc + title);
  1309. if (file.is_open())
  1310. {
  1311. file << sim_trial_duration_days.format(CSV_Fmt);
  1312. file.close();
  1313. std::cout << "done" << std::endl;
  1314. }
  1315. else
  1316. std::cout << "file could not be saved" << std::endl;
  1317. }
  1318. {
  1319. // save trial duration
  1320. std::string title = "trial_duration_months_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1321. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1322. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1323. + "_e"
  1324. + prep_type
  1325. + ".csv";
  1326. std::cout << "Saving file: " << title << "...";
  1327. std::ofstream file(plotsLoc + title);
  1328. if (file.is_open())
  1329. {
  1330. file << sim_trial_duration_months.format(CSV_Fmt);
  1331. file.close();
  1332. std::cout << "done" << std::endl;
  1333. }
  1334. else
  1335. std::cout << "file could not be saved" << std::endl;
  1336. }
  1337. {
  1338. // save exposed fraction
  1339. std::string title = "exposed_fraction_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1340. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1341. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1342. + "_e"
  1343. + prep_type
  1344. + ".csv";
  1345. std::cout << "Saving file: " << title << "...";
  1346. std::ofstream file(plotsLoc + title);
  1347. if (file.is_open())
  1348. {
  1349. file << sim_exposed_fraction.format(CSV_Fmt);
  1350. file.close();
  1351. std::cout << "done" << std::endl;
  1352. }
  1353. else
  1354. std::cout << "file could not be saved" << std::endl;
  1355. }
  1356. {
  1357. // save number of infections
  1358. std::string title = "num_infections_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1359. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1360. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1361. + "_e"
  1362. + prep_type
  1363. + ".csv";
  1364. std::cout << "Saving file: " << title << "...";
  1365. std::ofstream file(plotsLoc + title);
  1366. if (file.is_open())
  1367. {
  1368. file << sim_num_infections.format(CSV_Fmt);
  1369. file.close();
  1370. std::cout << "done" << std::endl;
  1371. }
  1372. else
  1373. std::cout << "file could not be saved" << std::endl;
  1374. }
  1375. {
  1376. // save number of infections
  1377. std::string title = "num_infections_sex_work_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1378. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1379. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1380. + "_e"
  1381. + prep_type
  1382. + ".csv";
  1383. std::cout << "Saving file: " << title << "...";
  1384. std::ofstream file(plotsLoc + title);
  1385. if (file.is_open())
  1386. {
  1387. file << sim_num_infections_sex_work.format(CSV_Fmt);
  1388. file.close();
  1389. std::cout << "done" << std::endl;
  1390. }
  1391. else
  1392. std::cout << "file could not be saved" << std::endl;
  1393. }
  1394. {
  1395. // save total follow up
  1396. std::string title = "follow_up_days_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1397. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1398. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1399. + "_e"
  1400. + prep_type
  1401. + ".csv";
  1402. std::cout << "Saving file: " << title << "...";
  1403. std::ofstream file(plotsLoc + title);
  1404. if (file.is_open())
  1405. {
  1406. file << sim_follow_up_days.format(CSV_Fmt);
  1407. file.close();
  1408. std::cout << "done" << std::endl;
  1409. }
  1410. else
  1411. std::cout << "file could not be saved" << std::endl;
  1412. }
  1413. {
  1414. // save number of exposed
  1415. std::string title = "num_exposed_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1416. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1417. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1418. + "_e"
  1419. + prep_type
  1420. + ".csv";
  1421. std::cout << "Saving file: " << title << "...";
  1422. std::ofstream file(plotsLoc + title);
  1423. if (file.is_open())
  1424. {
  1425. file << sim_num_exposed.format(CSV_Fmt);
  1426. file.close();
  1427. std::cout << "done" << std::endl;
  1428. }
  1429. else
  1430. std::cout << "file could not be saved" << std::endl;
  1431. }
  1432. {
  1433. // save number of exposures
  1434. std::string title = "num_exposures_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1435. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1436. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1437. + "_e"
  1438. + prep_type
  1439. + ".csv";
  1440. std::cout << "Saving file: " << title << "...";
  1441. std::ofstream file(plotsLoc + title);
  1442. if (file.is_open())
  1443. {
  1444. file << sim_num_exposures.format(CSV_Fmt);
  1445. file.close();
  1446. std::cout << "done" << std::endl;
  1447. }
  1448. else
  1449. std::cout << "file could not be saved" << std::endl;
  1450. }
  1451. {
  1452. // save number of exposures
  1453. std::string title = "num_exposures_sex_work_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1454. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1455. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1456. + "_e"
  1457. + prep_type
  1458. + ".csv";
  1459. std::cout << "Saving file: " << title << "...";
  1460. std::ofstream file(plotsLoc + title);
  1461. if (file.is_open())
  1462. {
  1463. file << sim_num_exposures_sex_work.format(CSV_Fmt);
  1464. file.close();
  1465. std::cout << "done" << std::endl;
  1466. }
  1467. else
  1468. std::cout << "file could not be saved" << std::endl;
  1469. }
  1470. {
  1471. // save cumulative general drop out
  1472. std::string title = "cum_gen_dropout_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1473. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1474. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1475. + "_e"
  1476. + prep_type
  1477. + ".csv";
  1478. std::cout << "Saving file: " << title << "...";
  1479. std::ofstream file(plotsLoc + title);
  1480. if (file.is_open())
  1481. {
  1482. file << sim_cum_drop_out_gen.format(CSV_Fmt);
  1483. file.close();
  1484. std::cout << "done" << std::endl;
  1485. }
  1486. else
  1487. std::cout << "file could not be saved" << std::endl;
  1488. }
  1489. {
  1490. // save cumulative general drop out
  1491. std::string title = "cum_fsw_dropout_fsw" + std::to_string(static_cast<int>(100 * fraction_fsw_participants(enrollment_scenario)))
  1492. + "_do" + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 0))) + "_"
  1493. + std::to_string(static_cast<int>(100 * drop_out_rates_annual(drop_out_scenario, 2)))
  1494. + "_e"
  1495. + prep_type
  1496. + ".csv";
  1497. std::cout << "Saving file: " << title << "...";
  1498. std::ofstream file(plotsLoc + title);
  1499. if (file.is_open())
  1500. {
  1501. file << sim_cum_drop_out_fsw.format(CSV_Fmt);
  1502. file.close();
  1503. std::cout << "done" << std::endl;
  1504. }
  1505. else
  1506. std::cout << "file could not be saved" << std::endl;
  1507. }
  1508. }
  1509. }// end enrollment scenarios
  1510. }// end drop out scenario
  1511. }// end prep scenarios
  1512. auto end_time = steady_clock::now(); // end time of program in milliseconds
  1513. auto program_duration = duration_cast<milliseconds>(end_time - start_time); // total duration of program in milliseconds
  1514.  
  1515. std::cout << "The model has finished running!" << std::endl;
  1516. std::cout << "Total runtime: " << program_duration.count() << " ms" << std::endl;
  1517. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!