From b293e21e8ee87f7efdfe50d47a9fb245c1fc727b Mon Sep 17 00:00:00 2001 From: chris-j-tang Date: Sat, 9 May 2020 15:49:34 -0400 Subject: [PATCH 1/7] Removing Cells/CellLookup/Mcells/McellLookup calloc --- src/CovidSim.cpp | 156 ++++++++++++++++---------------- src/Dist.cpp | 15 ++-- src/Model.h | 10 ++- src/SetupModel.cpp | 215 ++++++++++++++++++++++----------------------- 4 files changed, 199 insertions(+), 197 deletions(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index 6324a2390..1f1698826 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -71,10 +71,11 @@ param P; person* Hosts; household* Households; popvar State, StateT[MAX_NUM_THREADS]; -cell* Cells; // Cells[i] is the i'th cell -cell ** CellLookup; // CellLookup[i] is a pointer to the i'th populated cell -microcell* Mcells, ** McellLookup; -place** Places; +std::vector Cells(P.NC); // Cells.at(i) is the i'th cell +std::vector CellLookup(P.NCP); // CellLookup[i] is a pointer to the i'th populated cell +std::vector Mcells(P.NMC); +std::vector McellLookup(P.NMCP); +std::vector> Places; adminunit AdUnits[MAX_ADUNITS]; //// Time Series defs: //// TimeSeries is an array of type results, used to store (unsurprisingly) a time series of every quantity in results. Mostly used in RecordSample. @@ -2567,25 +2568,25 @@ void InitModel(int run) // passing run number so we can save run number in the i { for (i = tn; i < P.NC; i+=P.NumThreads) { - if ((Cells[i].tot_treat != 0) || (Cells[i].tot_vacc != 0) || (Cells[i].S != Cells[i].n) || (Cells[i].D > 0) || (Cells[i].R > 0)) + if ((Cells.at(i).tot_treat != 0) || (Cells.at(i).tot_vacc != 0) || (Cells.at(i).S != Cells.at(i).n) || (Cells.at(i).D > 0) || (Cells.at(i).R > 0)) { - for (j = 0; j < Cells[i].n; j++) + for (j = 0; j < Cells.at(i).n; j++) { - k = Cells[i].members[j]; - Cells[i].susceptible[j] = k; //added this in here instead + k = Cells.at(i).members[j]; + Cells.at(i).susceptible[j] = k; //added this in here instead Hosts[k].listpos = j; } - Cells[i].S = Cells[i].n; - Cells[i].L = Cells[i].I = Cells[i].R = Cells[i].cumTC = Cells[i].D = 0; - Cells[i].infected = Cells[i].latent = Cells[i].susceptible + Cells[i].S; - Cells[i].tot_treat = Cells[i].tot_vacc = 0; - for (l = 0; l < MAX_INTERVENTION_TYPES; l++) Cells[i].CurInterv[l] = -1; + Cells.at(i).S = Cells.at(i).n; + Cells.at(i).L = Cells.at(i).I = Cells.at(i).R = Cells.at(i).cumTC = Cells.at(i).D = 0; + Cells.at(i).infected = Cells.at(i).latent = Cells.at(i).susceptible + Cells.at(i).S; + Cells.at(i).tot_treat = Cells.at(i).tot_vacc = 0; + for (l = 0; l < MAX_INTERVENTION_TYPES; l++) Cells.at(i).CurInterv[l] = -1; // Next loop needs to count down for DoImmune host list reordering to work if(!P.DoPartialImmunity) - for (j = Cells[i].n - 1; j >= 0; j--) + for (j = Cells.at(i).n - 1; j >= 0; j--) { - k = Cells[i].members[j]; + k = Cells.at(i).members[j]; if (P.DoWholeHouseholdImmunity) { // note that this breaks determinism of runs if executed due to reordering of Cell members list each realisation @@ -2618,15 +2619,15 @@ void InitModel(int run) // passing run number so we can save run number in the i #pragma omp parallel for private(i,j,k,l) schedule(static,500) for (l = 0; l < P.NMCP; l++) { - i = (int)(McellLookup[l] - Mcells); - Mcells[i].vacc_start_time = Mcells[i].treat_start_time = USHRT_MAX - 1; - Mcells[i].treat_end_time = 0; - Mcells[i].treat_trig = Mcells[i].vacc_trig = Mcells[i].vacc = Mcells[i].treat = 0; - Mcells[i].place_trig = Mcells[i].move_trig = Mcells[i].socdist_trig = Mcells[i].keyworkerproph_trig = - Mcells[i].placeclose = Mcells[i].moverest = Mcells[i].socdist = Mcells[i].keyworkerproph = 0; - Mcells[i].move_start_time = USHRT_MAX - 1; - Mcells[i].place_end_time = Mcells[i].move_end_time = - Mcells[i].socdist_end_time = Mcells[i].keyworkerproph_end_time = 0; + i = (int)(McellLookup.at(l) - &Mcells[0]); + Mcells.at(i).vacc_start_time = Mcells.at(i).treat_start_time = USHRT_MAX - 1; + Mcells.at(i).treat_end_time = 0; + Mcells.at(i).treat_trig = Mcells.at(i).vacc_trig = Mcells.at(i).vacc = Mcells.at(i).treat = 0; + Mcells.at(i).place_trig = Mcells.at(i).move_trig = Mcells.at(i).socdist_trig = Mcells.at(i).keyworkerproph_trig = + Mcells.at(i).placeclose = Mcells.at(i).moverest = Mcells.at(i).socdist = Mcells.at(i).keyworkerproph = 0; + Mcells.at(i).move_start_time = USHRT_MAX - 1; + Mcells.at(i).place_end_time = Mcells.at(i).move_end_time = + Mcells.at(i).socdist_end_time = Mcells.at(i).keyworkerproph_end_time = 0; } if (P.DoPlaces) #pragma omp parallel for private(m,l) schedule(static,1) @@ -2771,7 +2772,7 @@ void SeedInfection(double t, int* nsi, int rf, int run) //adding run number to p m = 0; for (k = 0; (k < nsi[i]) && (m < 10000); k++) { - l = Mcells[j].members[(int)(ranf() * ((double)Mcells[j].n))]; //// randomly choose member of microcell j. Name this member l + l = Mcells.at(j).members[(int)(ranf() * ((double)Mcells.at(j).n))]; //// randomly choose member of microcell j. Name this member l if (Hosts[l].inf == InfStat_Susceptible) //// If Host l is uninfected. { if (CalcPersonSusc(l, 0, 0, 0) > 0) @@ -2801,13 +2802,13 @@ void SeedInfection(double t, int* nsi, int rf, int run) //adding run number to p { l = (int)(ranf() * ((double)P.PopSize)); j = Hosts[l].mcell; - //fprintf(stderr,"%i ",AdUnits[Mcells[j].adunit].id); - } while ((Mcells[j].n < nsi[i]) || (Mcells[j].n > P.MaxPopDensForInitialInfection) - || (Mcells[j].n < P.MinPopDensForInitialInfection) - || ((P.InitialInfectionsAdminUnit[i] > 0) && ((AdUnits[Mcells[j].adunit].id % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor != (P.InitialInfectionsAdminUnit[i] % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor))); + //fprintf(stderr,"%i ",AdUnits[Mcells.at(j).adunit].id); + } while ((Mcells.at(j).n < nsi[i]) || (Mcells.at(j).n > P.MaxPopDensForInitialInfection) + || (Mcells.at(j).n < P.MinPopDensForInitialInfection) + || ((P.InitialInfectionsAdminUnit[i] > 0) && ((AdUnits[Mcells.at(j).adunit].id % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor != (P.InitialInfectionsAdminUnit[i] % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor))); for (k = 0; (k < nsi[i]) && (m < 10000); k++) { - l = Mcells[j].members[(int)(ranf() * ((double)Mcells[j].n))]; + l = Mcells.at(j).members[(int)(ranf() * ((double)Mcells.at(j).n))]; if (Hosts[l].inf == InfStat_Susceptible) { if (CalcPersonSusc(l, 0, 0, 0) > 0) @@ -2839,11 +2840,11 @@ void SeedInfection(double t, int* nsi, int rf, int run) //adding run number to p { l = (int)(ranf() * ((double)P.PopSize)); j = Hosts[l].mcell; - //fprintf(stderr,"@@ %i %i ",AdUnits[Mcells[j].adunit].id, (int)(AdUnits[Mcells[j].adunit].id / P.CountryDivisor)); - } while ((Mcells[j].n == 0) || (Mcells[j].n > P.MaxPopDensForInitialInfection) - || (Mcells[j].n < P.MinPopDensForInitialInfection) - || ((P.InitialInfectionsAdminUnit[i] > 0) && ((AdUnits[Mcells[j].adunit].id % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor != (P.InitialInfectionsAdminUnit[i] % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor))); - l = Mcells[j].members[(int)(ranf() * ((double)Mcells[j].n))]; + //fprintf(stderr,"@@ %i %i ",AdUnits[Mcells.at(j).adunit].id, (int)(AdUnits[Mcells.at(j).adunit].id / P.CountryDivisor)); + } while ((Mcells.at(j).n == 0) || (Mcells.at(j).n > P.MaxPopDensForInitialInfection) + || (Mcells.at(j).n < P.MinPopDensForInitialInfection) + || ((P.InitialInfectionsAdminUnit[i] > 0) && ((AdUnits[Mcells.at(j).adunit].id % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor != (P.InitialInfectionsAdminUnit[i] % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor))); + l = Mcells.at(j).members[(int)(ranf() * ((double)Mcells.at(j).n))]; if (Hosts[l].inf == InfStat_Susceptible) { if (CalcPersonSusc(l, 0, 0, 0) > 0) @@ -3997,10 +3998,10 @@ void LoadSnapshot(void) if (!(Array_tot_prob = (float*)malloc(P.NCP * sizeof(float)))) ERR_CRITICAL("Unable to temp allocate cell storage\n"); for (i = 0; i < P.NCP; i++) { - Array_InvCDF[i] = Cells[i].InvCDF; - Array_max_trans[i] = Cells[i].max_trans; - Array_cum_trans[i] = Cells[i].cum_trans; - Array_tot_prob[i] = Cells[i].tot_prob; + Array_InvCDF[i] = Cells.at(i).InvCDF; + Array_max_trans[i] = Cells.at(i).max_trans; + Array_cum_trans[i] = Cells.at(i).cum_trans; + Array_tot_prob[i] = Cells.at(i).tot_prob; } fread_big((void*)& i, sizeof(int), 1, dat); if (i != P.PopSize) ERR_CRITICAL_FMT("Incorrect N (%i %i) in snapshot file.\n", P.PopSize, i); @@ -4026,9 +4027,9 @@ void LoadSnapshot(void) fprintf(stderr, "."); zfread_big((void*)Households, sizeof(household), (size_t)P.NH, dat); fprintf(stderr, "."); - zfread_big((void*)Cells, sizeof(cell), (size_t)P.NC, dat); + zfread_big((void*)&Cells.front(), sizeof(cell), (size_t)P.NC, dat); fprintf(stderr, "."); - zfread_big((void*)Mcells, sizeof(microcell), (size_t)P.NMC, dat); + zfread_big((void*)&Mcells.front(), sizeof(microcell), (size_t)P.NMC, dat); fprintf(stderr, "."); zfread_big((void*)State.CellMemberArray, sizeof(int), (size_t)P.PopSize, dat); fprintf(stderr, "."); @@ -4036,25 +4037,25 @@ void LoadSnapshot(void) fprintf(stderr, "."); for (i = 0; i < P.NC; i++) { - if (Cells[i].n > 0) + if (Cells.at(i).n > 0) { - Cells[i].members += CM_offset; - Cells[i].susceptible += CSM_offset; - Cells[i].latent += CSM_offset; - Cells[i].infected += CSM_offset; + Cells.at(i).members += CM_offset; + Cells.at(i).susceptible += CSM_offset; + Cells.at(i).latent += CSM_offset; + Cells.at(i).infected += CSM_offset; } - for (j = 0; j < MAX_INTERVENTION_TYPES; j++) Cells[i].CurInterv[j] = -1; // turn interventions off in loaded image + for (j = 0; j < MAX_INTERVENTION_TYPES; j++) Cells.at(i).CurInterv[j] = -1; // turn interventions off in loaded image } for (i = 0; i < P.NMC; i++) - if (Mcells[i].n > 0) - Mcells[i].members += CM_offset; + if (Mcells.at(i).n > 0) + Mcells.at(i).members += CM_offset; for (i = 0; i < P.NCP; i++) { - Cells[i].InvCDF = Array_InvCDF[i]; - Cells[i].max_trans = Array_max_trans[i]; - Cells[i].cum_trans = Array_cum_trans[i]; - Cells[i].tot_prob = Array_tot_prob[i]; + Cells.at(i).InvCDF = Array_InvCDF[i]; + Cells.at(i).max_trans = Array_max_trans[i]; + Cells.at(i).cum_trans = Array_cum_trans[i]; + Cells.at(i).tot_prob = Array_tot_prob[i]; } free(Array_tot_prob); free(Array_cum_trans); @@ -4101,9 +4102,9 @@ void SaveSnapshot(void) fprintf(stderr, "## %i\n", i++); zfwrite_big((void*)Households, sizeof(household), (size_t)P.NH, dat); fprintf(stderr, "## %i\n", i++); - zfwrite_big((void*)Cells, sizeof(cell), (size_t)P.NC, dat); + zfwrite_big((void*)&Cells.front(), sizeof(cell), (size_t)P.NC, dat); fprintf(stderr, "## %i\n", i++); - zfwrite_big((void*)Mcells, sizeof(microcell), (size_t)P.NMC, dat); + zfwrite_big((void*)&Mcells.front(), sizeof(microcell), (size_t)P.NMC, dat); fprintf(stderr, "## %i\n", i++); zfwrite_big((void*)State.CellMemberArray, sizeof(int), (size_t)P.PopSize, dat); @@ -4123,12 +4124,12 @@ void UpdateProbs(int DoPlace) #pragma omp parallel for private(j) schedule(static,500) for (j = 0; j < P.NCP; j++) { - CellLookup[j]->tot_prob = 0; - CellLookup[j]->S0 = CellLookup[j]->S + CellLookup[j]->L + CellLookup[j]->I; + Cells.at(j).tot_prob = 0; + Cells.at(j).S0 = Cells.at(j).S + Cells.at(j).L + Cells.at(j).I; if (P.DoDeath) { - CellLookup[j]->S0 += CellLookup[j]->n / 5; - if ((CellLookup[j]->n < 100) || (CellLookup[j]->S0 > CellLookup[j]->n)) CellLookup[j]->S0 = CellLookup[j]->n; + Cells.at(j).S0 += Cells.at(j).n / 5; + if ((Cells.at(j).n < 100) || (Cells.at(j).S0 > Cells.at(j).n)) Cells.at(j).S0 = Cells.at(j).n; } } } @@ -4137,8 +4138,8 @@ void UpdateProbs(int DoPlace) #pragma omp parallel for private(j) schedule(static,500) for (j = 0; j < P.NCP; j++) { - CellLookup[j]->S0 = CellLookup[j]->S; - CellLookup[j]->tot_prob = 0; + Cells.at(j).S0 = Cells.at(j).S; + Cells.at(j).tot_prob = 0; } } #pragma omp parallel for private(j) schedule(static,500) @@ -4146,21 +4147,21 @@ void UpdateProbs(int DoPlace) { int m, k; float t; - CellLookup[j]->cum_trans[0] = ((float)(CellLookup[0]->S0)) * CellLookup[j]->max_trans[0]; - t = ((float)CellLookup[0]->n) * CellLookup[j]->max_trans[0]; + Cells.at(j).cum_trans[0] = ((float)(Cells.front().S0)) * Cells.at(j).max_trans[0]; + t = ((float)Cells.front().n) * Cells.at(j).max_trans[0]; for (m = 1; m < P.NCP; m++) { - CellLookup[j]->cum_trans[m] = CellLookup[j]->cum_trans[m - 1] + ((float)(CellLookup[m]->S0)) * CellLookup[j]->max_trans[m]; - t += ((float)CellLookup[m]->n) * CellLookup[j]->max_trans[m]; + Cells.at(j).cum_trans[m] = Cells.at(j).cum_trans[m - 1] + ((float)(Cells.at(m).S0)) * Cells.at(j).max_trans[m]; + t += ((float)Cells.at(m).n) * Cells.at(j).max_trans[m]; } - CellLookup[j]->tot_prob = CellLookup[j]->cum_trans[P.NCP - 1]; + Cells.at(j).tot_prob = Cells.at(j).cum_trans[P.NCP - 1]; for (m = 0; m < P.NCP; m++) - CellLookup[j]->cum_trans[m] /= CellLookup[j]->tot_prob; - CellLookup[j]->tot_prob /= t; + Cells.at(j).cum_trans[m] /= Cells.at(j).tot_prob; + Cells.at(j).tot_prob /= t; for (k = m = 0; k <= 1024; k++) { - while (CellLookup[j]->cum_trans[m] * 1024 < ((float)k)) m++; - CellLookup[j]->InvCDF[k] = m; + while (Cells.at(j).cum_trans[m] * 1024 < ((float)k)) m++; + Cells.at(j).InvCDF[k] = m; } } } @@ -4658,7 +4659,7 @@ void RecordSample(double t, int n) State.NumPlacesClosed[i] = numPC; TimeSeries[n].PropPlacesClosed[i] = ((double)numPC) / ((double)P.Nplace[i]); } - for (i = k = 0; i < P.NMC; i++) if (Mcells[i].socdist == 2) k++; + for (i = k = 0; i < P.NMC; i++) if (Mcells.at(i).socdist == 2) k++; TimeSeries[n].PropSocDist=((double)k)/((double)P.NMC); //update contact number distribution in State @@ -5061,13 +5062,14 @@ void CalcOriginDestMatrix_adunit() #pragma omp parallel for private(tn,i,j,k,l,m,p,total_flow,mcl_from,mcl_to,cl_from,cl_to,cl_from_mcl,cl_to_mcl,flow) schedule(static) //reduction(+:s,t2) for (tn = 0; tn < P.NumThreads; tn++) { + auto Cells_front = &Cells.front(); for (i = tn; i < P.NCP; i += P.NumThreads) { //reset pop density matrix to zero double pop_dens_from[MAX_ADUNITS] = {}; //find index of cell from which flow travels - cl_from = CellLookup[i] - Cells; + cl_from = CellLookup[i] - Cells_front; cl_from_mcl = (cl_from / P.nch) * P.NMCL * P.nmch + (cl_from % P.nch) * P.NMCL; //loop over microcells in these cells to find populations in each admin unit and so flows @@ -5091,17 +5093,17 @@ void CalcOriginDestMatrix_adunit() double pop_dens_to[MAX_ADUNITS] = {}; //find index of cell which flow travels to - cl_to = CellLookup[j] - Cells; + cl_to = CellLookup[j] - Cells_front; cl_to_mcl = (cl_to / P.nch) * P.NMCL * P.nmch + (cl_to % P.nch) * P.NMCL; //calculate distance and kernel between the cells //total_flow=Cells[cl_from].max_trans[j]*Cells[cl_from].n*Cells[cl_to].n; if (j == 0) { - total_flow = Cells[cl_from].cum_trans[j] * Cells[cl_from].n; + total_flow = Cells.at(cl_from).cum_trans[j] * Cells.at(cl_from).n; } else { - total_flow = (Cells[cl_from].cum_trans[j] - Cells[cl_from].cum_trans[j - 1]) * Cells[cl_from].n; + total_flow = (Cells.at(cl_from).cum_trans[j] - Cells.at(cl_from).cum_trans[j - 1]) * Cells.at(cl_from).n; } //loop over microcells within destination cell @@ -5111,10 +5113,10 @@ void CalcOriginDestMatrix_adunit() { //get index of microcell mcl_to = cl_to_mcl + p + m * P.nmch; - if (Mcells[mcl_to].n > 0) + if (Mcells.at(mcl_to).n > 0) { //get proportion of each population of cell that exists in each admin unit - pop_dens_to[Mcells[mcl_to].adunit] += (((double)Mcells[mcl_to].n) / ((double)Cells[cl_to].n)); + pop_dens_to[Mcells.at(mcl_to).adunit] += (((double)Mcells.at(mcl_to).n) / ((double)Cells.at(cl_to).n)); } } } diff --git a/src/Dist.cpp b/src/Dist.cpp index 5e4fe65e5..40879a83b 100644 --- a/src/Dist.cpp +++ b/src/Dist.cpp @@ -62,8 +62,9 @@ double dist2_cc(cell* a, cell* b) double x, y; int l, m; - l = (int)(a - Cells); - m = (int)(b - Cells); + auto Cells_front = &Cells.front(); + l = (int)(a - Cells_front); + m = (int)(b - Cells_front); if (P.DoUTM_coords) return dist2UTM(P.cwidth * fabs((double)(l / P.nch)), P.cheight * fabs((double)(l % P.nch)), P.cwidth * fabs((double)(m / P.nch)), P.cheight * fabs((double)(m % P.nch))); @@ -84,8 +85,9 @@ double dist2_cc_min(cell* a, cell* b) double x, y; int l, m, i, j; - l = (int)(a - Cells); - m = (int)(b - Cells); + auto Cells_front = &Cells.front(); + l = (int)(a - Cells_front); + m = (int)(b - Cells_front); i = l; j = m; if (P.DoUTM_coords) { @@ -155,8 +157,9 @@ double dist2_mm(microcell* a, microcell* b) double x, y; int l, m; - l = (int)(a - Mcells); - m = (int)(b - Mcells); + auto Mcells_front = &Mcells.front(); + l = (int)(a - Mcells_front); + m = (int)(b - Mcells_front); if (P.DoUTM_coords) return dist2UTM(P.mcwidth * fabs((double)(l / P.nmch)), P.mcheight * fabs((double)(l % P.nmch)), P.mcwidth * fabs((double)(m / P.nmch)), P.mcheight * fabs((double)(m % P.nmch))); diff --git a/src/Model.h b/src/Model.h index 801b185c3..fcc98768a 100644 --- a/src/Model.h +++ b/src/Model.h @@ -3,6 +3,8 @@ #ifndef COVIDSIM_MODEL_H_INCLUDED_ #define COVIDSIM_MODEL_H_INCLUDED_ +#include + #include "Country.h" #include "MachineDefines.h" #include "Constants.h" @@ -334,9 +336,11 @@ typedef struct ADMINUNIT { extern person* Hosts; extern household* Households; extern popvar State, StateT[MAX_NUM_THREADS]; -extern cell* Cells, ** CellLookup; -extern microcell* Mcells, ** McellLookup; -extern place** Places; +extern std::vector Cells; +extern std::vector CellLookup; +extern std::vector Mcells; +extern std::vector McellLookup; +extern std::vector> Places; extern adminunit AdUnits[MAX_ADUNITS]; //// Time Series defs: diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp index ce6f9f12d..37e2ccb43 100644 --- a/src/SetupModel.cpp +++ b/src/SetupModel.cpp @@ -297,9 +297,11 @@ void SetupModel(char* DensityFile, char* NetworkFile, char* SchoolFile, char* Re StratifyPlaces(); for (i = 0; i < P.NC; i++) { - Cells[i].S = Cells[i].n; - Cells[i].L = Cells[i].I = Cells[i].R = 0; - //Cells[i].susceptible=Cells[i].members; //added this line + cell cell_i; + cell_i.S = cell_i.n; + cell_i.L = cell_i.I = cell_i.R = 0; + Cells.push_back(cell_i); + //Cells.at(i).susceptible=Cells.at(i).members; //added this line } for (i = 0; i < P.PopSize; i++) Hosts[i].keyworker = 0; P.KeyWorkerNum = P.KeyWorkerIncHouseNum = m = l = 0; @@ -577,7 +579,7 @@ void SetupModel(char* DensityFile, char* NetworkFile, char* SchoolFile, char* Re for (j = 0; j <= MAX_HOUSEHOLD_SIZE; j++) inf_household_av[i][j] = case_household_av[i][j] = 0; DoInitUpdateProbs = 1; - for (i = 0; i < P.NC; i++) Cells[i].tot_treat = 1; //This makes sure InitModel intialises the cells. + for (auto cell_i : Cells) cell_i.tot_treat = 1; //This makes sure InitModel intialises the cells. P.NRactE = P.NRactNE = 0; for (i = 0; i < P.PopSize; i++) Hosts[i].esocdist_comply = (ranf() < P.EnhancedSocDistProportionCompliant[HOST_AGE_GROUP(i)]) ? 1 : 0; if (!P.EnhancedSocDistClusterByHousehold) @@ -641,11 +643,12 @@ void SetupModel(char* DensityFile, char* NetworkFile, char* SchoolFile, char* Re PeakHeightSum = PeakHeightSS = PeakTimeSum = PeakTimeSS = 0; i = (P.ncw / 2) * P.nch + P.nch / 2; j = (P.ncw / 2 + 2) * P.nch + P.nch / 2; - fprintf(stderr, "UTM dist horiz=%lg %lg\n", sqrt(dist2_cc(Cells + i, Cells + j)), sqrt(dist2_cc(Cells + j, Cells + i))); + auto Cells_front = &Cells.front(); + fprintf(stderr, "UTM dist horiz=%lg %lg\n", sqrt(dist2_cc(Cells_front + i, Cells_front + j)), sqrt(dist2_cc(Cells_front + j, Cells_front + i))); j = (P.ncw / 2) * P.nch + P.nch / 2 + 2; - fprintf(stderr, "UTM dist vert=%lg %lg\n", sqrt(dist2_cc(Cells + i, Cells + j)), sqrt(dist2_cc(Cells + j, Cells + i))); + fprintf(stderr, "UTM dist vert=%lg %lg\n", sqrt(dist2_cc(Cells_front + i, Cells_front + j)), sqrt(dist2_cc(Cells_front + j, Cells_front + i))); j = (P.ncw / 2 + 2) * P.nch + P.nch / 2 + 2; - fprintf(stderr, "UTM dist diag=%lg %lg\n", sqrt(dist2_cc(Cells + i, Cells + j)), sqrt(dist2_cc(Cells + j, Cells + i))); + fprintf(stderr, "UTM dist diag=%lg %lg\n", sqrt(dist2_cc(Cells_front + i, Cells_front + j)), sqrt(dist2_cc(Cells_front + j, Cells_front + i))); //if(P.OutputBitmap) //{ @@ -664,22 +667,18 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) const char delimiters[] = " \t,"; FILE* dat = NULL, *dat2; bin_file rec; - double *mcell_dens; - int *mcell_adunits, *mcell_num, *mcell_country; - - if (!(Cells = (cell*)calloc(P.NC, sizeof(cell)))) ERR_CRITICAL("Unable to allocate cell storage\n"); - if (!(Mcells = (microcell*)calloc(P.NMC, sizeof(microcell)))) ERR_CRITICAL("Unable to allocate cell storage\n"); - if (!(mcell_num = (int*)malloc(P.NMC * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); - if (!(mcell_dens = (double*)malloc(P.NMC * sizeof(double)))) ERR_CRITICAL("Unable to allocate cell storage\n"); - if (!(mcell_country = (int*)malloc(P.NMC * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); - if (!(mcell_adunits = (int*)malloc(P.NMC * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); + std::vector mcell_dens(P.NMC); // Pre-allocates vector of P.NMC double elements + std::vector mcell_adunits(P.NMC), mcell_num(P.NMC), mcell_country(P.NMC); for (j = 0; j < P.NMC; j++) { - Mcells[j].n = 0; - mcell_adunits[j] = -1; - mcell_dens[j] = 0; - mcell_num[j] = mcell_country[j] = 0; + microcell mcell_j; + mcell_j.n = 0; + Mcells.push_back(mcell_j); + mcell_adunits.push_back(-1); + mcell_dens.push_back(0.0); + mcell_num.push_back(0); + mcell_country.push_back(0); } if (P.DoAdUnits) for (i = 0; i < MAX_ADUNITS; i++) @@ -740,19 +739,19 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) if (l < P.NMC) { mr++; - mcell_dens[l] += t; - mcell_country[l] = country; + mcell_dens.at(l) += t; + mcell_country.at(l) = country; //fprintf(stderr,"mcell %i, country %i, pop %lg\n",l,country,t); - mcell_num[l]++; + mcell_num.at(l)++; if (P.DoAdUnits) { - mcell_adunits[l] = P.AdunitLevel1Lookup[m]; - if (mcell_adunits[l] < 0) fprintf(stderr, "Cell %i has adunits<0\n", l); + mcell_adunits.at(l) = P.AdunitLevel1Lookup[m]; + if (mcell_adunits.at(l) < 0) fprintf(stderr, "Cell %i has adunits<0\n", l); P.PopByAdunit[P.AdunitLevel1Lookup[m]][0] += t; } else - mcell_adunits[l] = 0; - if ((P.OutputDensFile) && (P.DoBin) && (mcell_adunits[l] >= 0)) + mcell_adunits.at(l) = 0; + if ((P.OutputDensFile) && (P.DoBin) && (mcell_adunits.at(l) >= 0)) { if (rn2 < rn) BF[rn2] = rec; rn2++; @@ -771,19 +770,19 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) P.DoBin = 1; P.BinFileLen = 0; for (l = 0; l < P.NMC; l++) - if (mcell_adunits[l] >= 0) P.BinFileLen++; + if (mcell_adunits.at(l) >= 0) P.BinFileLen++; if (!(BinFileBuf = (void*)malloc(P.BinFileLen * sizeof(bin_file)))) ERR_CRITICAL("Unable to allocate binary file buffer\n"); BF = (bin_file*)BinFileBuf; fprintf(stderr, "Binary density file should contain %i cells.\n", (int)P.BinFileLen); rn = 0; for (l = 0; l < P.NMC; l++) - if (mcell_adunits[l] >= 0) + if (mcell_adunits.at(l) >= 0) { BF[rn].x = (double)(P.mcwidth * (((double)(l / P.nmch)) + 0.5)) + P.SpatialBoundingBox[0]; //x BF[rn].y = (double)(P.mcheight * (((double)(l % P.nmch)) + 0.5)) + P.SpatialBoundingBox[1]; //y - BF[rn].ad = (P.DoAdUnits) ? (AdUnits[mcell_adunits[l]].id) : 0; - BF[rn].pop = mcell_dens[l]; - BF[rn].cnt = mcell_country[l]; + BF[rn].ad = (P.DoAdUnits) ? (AdUnits[mcell_adunits.at(l)].id) : 0; + BF[rn].pop = mcell_dens.at(l); + BF[rn].cnt = mcell_country.at(l); rn++; } } @@ -804,17 +803,17 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) maxd = 0; for (i = 0; i < P.NMC; i++) { - if (mcell_num[i] > 0) + if (mcell_num.at(i) > 0) { - mcell_dens[i] /= ((double)mcell_num[i]); - Mcells[i].country = (unsigned short)mcell_country[i]; + mcell_dens[i] /= ((double)mcell_num.at(i)); + Mcells.at(i).country = (unsigned short)mcell_country[i]; if (P.DoAdUnits) - Mcells[i].adunit = mcell_adunits[i]; + Mcells.at(i).adunit = mcell_adunits[i]; else - Mcells[i].adunit = 0; + Mcells.at(i).adunit = 0; } else - Mcells[i].adunit = -1; + Mcells.at(i).adunit = -1; maxd += mcell_dens[i]; } } @@ -823,7 +822,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) for (i = 0; i < P.NMC; i++) { mcell_dens[i] = 1.0; - Mcells[i].country = 1; + Mcells.at(i).country = 1; } maxd = ((double)P.NMC); } @@ -937,7 +936,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) maxd = 0; for (i = 0; i < P.NMC; i++) { - if (mcell_num[i] > 0) + if (mcell_num.at(i) > 0) { if (mcell_adunits[i] < 0) ERR_CRITICAL_FMT("Cell %i has adunits < 0 (indexing PopByAdunit)\n", i); mcell_dens[i] *= P.PopByAdunit[mcell_adunits[i]][1] / (1e-10 + P.PopByAdunit[mcell_adunits[i]][0]); @@ -956,12 +955,12 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) { s = mcell_dens[i] / maxd / t; if (s > 1.0) s = 1.0; - m += (Mcells[i].n = (int)ignbin_mt((long)(P.PopSize - m), s, 0)); + m += (Mcells.at(i).n = (int)ignbin_mt((long)(P.PopSize - m), s, 0)); t -= mcell_dens[i] / maxd; - if (Mcells[i].n > 0) { + if (Mcells.at(i).n > 0) { P.NMCP++; if (mcell_adunits[i] < 0) ERR_CRITICAL_FMT("Cell %i has adunits < 0 (indexing AdUnits)\n", i); - AdUnits[mcell_adunits[i]].n += Mcells[i].n; + AdUnits[mcell_adunits[i]].n += Mcells.at(i).n; } } Mcells[P.NMC - 1].n = P.PopSize - m; @@ -971,21 +970,16 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) AdUnits[mcell_adunits[P.NMC - 1]].n += Mcells[P.NMC - 1].n; } - free(mcell_dens); - free(mcell_num); - free(mcell_country); - free(mcell_adunits); t = 0.0; - if (!(McellLookup = (microcell * *)malloc(P.NMCP * sizeof(microcell*)))) ERR_CRITICAL("Unable to allocate cell storage\n"); if (!(mcl = (int*)malloc(P.PopSize * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); State.CellMemberArray = mcl; P.NCP = 0; for (i = i2 = j2 = 0; i < P.NC; i++) { - Cells[i].n = 0; + Cells.at(i).n = 0; k = (i / P.nch) * P.NMCL * P.nmch + (i % P.nch) * P.NMCL; - Cells[i].members = mcl + j2; + Cells.at(i).members = mcl + j2; for (l = 0; l < P.NMCL; l++) for (m = 0; m < P.NMCL; m++) { @@ -994,19 +988,18 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) { Mcells[j].members = mcl + j2; //if(!(Mcells[j].members=(int *) calloc(Mcells[j].n,sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); //replaced line above with this to ensure members don't get mixed across microcells - McellLookup[i2++] = Mcells + j; - Cells[i].n += Mcells[j].n; + McellLookup[i2++] = &Mcells.front() + j; + Cells.at(i).n += Mcells[j].n; j2 += Mcells[j].n; } } - if (Cells[i].n > 0) P.NCP++; + if (Cells.at(i).n > 0) P.NCP++; } fprintf(stderr, "Number of hosts assigned = %i\n", j2); if (!P.DoAdUnits) P.AdunitLevel1Lookup[0] = 0; fprintf(stderr, "Number of cells with non-zero population = %i\n", P.NCP); fprintf(stderr, "Number of microcells with non-zero population = %i\n", P.NMCP); - if (!(CellLookup = (cell * *)malloc(P.NCP * sizeof(cell*)))) ERR_CRITICAL("Unable to allocate cell storage\n"); if (!(mcl = (int*)malloc(P.PopSize * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); State.CellSuscMemberArray = mcl; i2 = k = 0; @@ -1034,8 +1027,8 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) } for (i = 0; i < P.NC; i++) { - Cells[i].cumTC = 0; - for (j = 0; j < Cells[i].n; j++) Cells[i].members[j] = -1; + Cells.at(i).cumTC = 0; + for (j = 0; j < Cells.at(i).n; j++) Cells.at(i).members[j] = -1; } fprintf(stderr, "Cells assigned\n"); for (i = 0; i <= MAX_HOUSEHOLD_SIZE; i++) denom_household[i] = 0; @@ -1232,7 +1225,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) maxd = ((double)P.PopSize); last_i = 0; for (i = 0; i < P.NMC; i++) - if (Mcells[i].n > 0) last_i = i; + if (Mcells.at(i).n > 0) last_i = i; fprintf(stderr, "Allocating place/age groups...\n"); for (k = 0; k < NUM_AGE_GROUPS * AGE_GROUP_WIDTH; k++) { @@ -1279,7 +1272,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) for (i = 0; i < m; i++) if (!(Places[j][i].AvailByAge = (unsigned short int*) malloc(P.PlaceTypeMaxAgeRead[j] * sizeof(unsigned short int)))) ERR_CRITICAL("Unable to allocate place storage\n"); P.Nplace[j] = 0; - for (i = 0; i < P.NMC; i++) Mcells[i].np[j] = 0; + for (i = 0; i < P.NMC; i++) Mcells.at(i).np[j] = 0; } mr = 0; while (!feof(dat)) @@ -1291,7 +1284,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) if ((x >= P.SpatialBoundingBox[0]) && (x < P.SpatialBoundingBox[2]) && (y >= P.SpatialBoundingBox[1]) && (y < P.SpatialBoundingBox[3])) { i = P.nch * ((int)(Places[j][P.Nplace[j]].loc_x / P.cwidth)) + ((int)(Places[j][P.Nplace[j]].loc_y / P.cheight)); - if (Cells[i].n == 0) mr++; + if (Cells.at(i).n == 0) mr++; Places[j][P.Nplace[j]].n = m; i = (int)(Places[j][P.Nplace[j]].loc_x / P.mcwidth); k = (int)(Places[j][P.Nplace[j]].loc_y / P.mcheight); @@ -1306,10 +1299,10 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) fprintf(stderr, "%i schools read (%i in empty cells) \n", P.Nplace[j], mr); for (i = 0; i < P.NMC; i++) for (j = 0; j < P.nsp; j++) - if (Mcells[i].np[j] > 0) + if (Mcells.at(i).np[j] > 0) { - if (!(Mcells[i].places[j] = (int*)malloc(Mcells[i].np[j] * sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n"); - Mcells[i].np[j] = 0; + if (!(Mcells.at(i).places[j] = (int*)malloc(Mcells.at(i).np[j] * sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n"); + Mcells.at(i).np[j] = 0; } for (j = 0; j < P.nsp; j++) { @@ -1345,19 +1338,19 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) t = 1.0; for (m = i = k = 0; i < P.NMC; i++) { - s = ((double) Mcells[i].n) / maxd / t; + s = ((double) Mcells.at(i).n) / maxd / t; if (s > 1.0) s = 1.0; if (i == last_i) m += (Mcells[last_i].np[j2] = P.Nplace[j2] - m); else - m += (Mcells[i].np[j2] = (int)ignbin_mt((long)(P.Nplace[j2] - m), s, tn)); - t -= ((double)Mcells[i].n) / maxd; - if (Mcells[i].np[j2] > 0) + m += (Mcells.at(i).np[j2] = (int)ignbin_mt((long)(P.Nplace[j2] - m), s, tn)); + t -= ((double)Mcells.at(i).n) / maxd; + if (Mcells.at(i).np[j2] > 0) { - if (!(Mcells[i].places[j2] = (int*)malloc(Mcells[i].np[j2] * sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n"); + if (!(Mcells.at(i).places[j2] = (int*)malloc(Mcells.at(i).np[j2] * sizeof(int)))) ERR_CRITICAL("Unable to allocate place storage\n"); x = (double)(i / P.nmch); y = (double)(i % P.nmch); - for (j = 0; j < Mcells[i].np[j2]; j++) + for (j = 0; j < Mcells.at(i).np[j2]; j++) { xh = P.mcwidth * (ranf_mt(tn) + x); yh = P.mcheight * (ranf_mt(tn) + y); @@ -1365,8 +1358,8 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) Places[j2][k].loc_y = (float)yh; Places[j2][k].n = 0; Places[j2][k].mcell = i; - Places[j2][k].country = Mcells[i].country; - Mcells[i].places[j2][j] = k; + Places[j2][k].country = Mcells.at(i).country; + Mcells.at(i).places[j2][j] = k; k++; } } @@ -1383,8 +1376,8 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) } /* for (j2 = 0; j2 < P.PlaceTypeNum; j2++) for (i =0; i < P.NMC; i++) - if ((Mcells[i].np[j2]>0) && (Mcells[i].n == 0)) - fprintf(stderr, "\n##~ %i %i %i \n", i, j2, Mcells[i].np[j2]); + if ((Mcells.at(i).np[j2]>0) && (Mcells.at(i).n == 0)) + fprintf(stderr, "\n##~ %i %i %i \n", i, j2, Mcells.at(i).np[j2]); */ fprintf(stderr, "Places assigned\n"); } l = 0; @@ -1439,9 +1432,9 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) for (i = 0; i < P.NC; i++) { - Cells[i].cumTC = 0; - Cells[i].S = Cells[i].n; - Cells[i].L = Cells[i].I = 0; + Cells.at(i).cumTC = 0; + Cells.at(i).S = Cells.at(i).n; + Cells.at(i).L = Cells.at(i).I = 0; } fprintf(stderr, "Allocated cell and host memory\n"); fprintf(stderr, "Assigned hosts to cells\n"); @@ -1468,14 +1461,14 @@ void SetupAirports(void) for (i = 0; i < P.Nairports; i++) Airports[i].num_mcell = 0; cur = base; for (i = 0; i < P.NMC; i++) - if (Mcells[i].n > 0) + if (Mcells.at(i).n > 0) { - Mcells[i].AirportList = cur; + Mcells.at(i).AirportList = cur; cur += NNA; } #pragma omp parallel for private(i,j,k,l,x,y,t,tmin) schedule(static,10000) for (i = 0; i < P.NMC; i++) - if (Mcells[i].n > 0) + if (Mcells.at(i).n > 0) { if (i % 10000 == 0) fprintf(stderr, "\n%i ", i); x = (((double)(i / P.nmch)) + 0.5) * P.mcwidth; @@ -1488,26 +1481,26 @@ void SetupAirports(void) t = numKernel(dist2_raw(x, y, Airports[j].loc_x, Airports[j].loc_y)) * Airports[j].total_traffic; if (k < NNA) { - Mcells[i].AirportList[k].id = j; - Mcells[i].AirportList[k].prob = (float)t; + Mcells.at(i).AirportList[k].id = j; + Mcells.at(i).AirportList[k].prob = (float)t; if (t < tmin) { tmin = t; l = k; } k++; } else if (t > tmin) { - Mcells[i].AirportList[l].id = j; - Mcells[i].AirportList[l].prob = (float)t; + Mcells.at(i).AirportList[l].id = j; + Mcells.at(i).AirportList[l].prob = (float)t; tmin = 1e20; for (k = 0; k < NNA; k++) - if (Mcells[i].AirportList[k].prob < tmin) + if (Mcells.at(i).AirportList[k].prob < tmin) { - tmin = Mcells[i].AirportList[k].prob; + tmin = Mcells.at(i).AirportList[k].prob; l = k; } } } for (j = 0; j < NNA; j++) - Airports[Mcells[i].AirportList[j].id].num_mcell++; + Airports[Mcells.at(i).AirportList[j].id].num_mcell++; } cur = Airports[0].DestMcells; fprintf(stderr, "Microcell airport lists collated.\n"); @@ -1519,24 +1512,24 @@ void SetupAirports(void) } #pragma omp parallel for private(i,j,k,l,t,tmin) schedule(static,10000) for (i = 0; i < P.NMC; i++) - if (Mcells[i].n > 0) + if (Mcells.at(i).n > 0) { if (i % 10000 == 0) fprintf(stderr, "\n%i ", i); t = 0; for (j = 0; j < NNA; j++) { - t += Mcells[i].AirportList[j].prob; - k = Mcells[i].AirportList[j].id; + t += Mcells.at(i).AirportList[j].prob; + k = Mcells.at(i).AirportList[j].id; #pragma omp critical (airport) l = (Airports[k].num_mcell++); Airports[k].DestMcells[l].id = i; - Airports[k].DestMcells[l].prob = Mcells[i].AirportList[j].prob * ((float)Mcells[i].n); + Airports[k].DestMcells[l].prob = Mcells.at(i).AirportList[j].prob * ((float)Mcells.at(i).n); } tmin = 0; for (j = 0; j < NNA; j++) { - Mcells[i].AirportList[j].prob = (float)(tmin + Mcells[i].AirportList[j].prob / t); - tmin = Mcells[i].AirportList[j].prob; + Mcells.at(i).AirportList[j].prob = (float)(tmin + Mcells.at(i).AirportList[j].prob / t); + tmin = Mcells.at(i).AirportList[j].prob; } } fprintf(stderr, "Airport microcell lists collated.\n"); @@ -1815,9 +1808,9 @@ void AssignPeopleToPlaces(void) fprintf(stderr, "Assigning people to places....\n"); for (i = 0; i < P.NC; i++) { - Cells[i].infected = Cells[i].susceptible; - if (!(Cells[i].susceptible = (int*)calloc(Cells[i].n, sizeof(int)))) ERR_CRITICAL("Unable to allocate state storage\n"); - Cells[i].cumTC = Cells[i].n; + Cells.at(i).infected = Cells.at(i).susceptible; + if (!(Cells.at(i).susceptible = (int*)calloc(Cells.at(i).n, sizeof(int)))) ERR_CRITICAL("Unable to allocate state storage\n"); + Cells.at(i).cumTC = Cells.at(i).n; } //PropPlaces initialisation is only valid for non-overlapping places. @@ -1987,7 +1980,7 @@ void AssignPeopleToPlaces(void) } m += ((int)Places[tp][i].treat_end_time); } - for (i = 0; i < P.NC; i++) Cells[i].L = Cells[i].I = Cells[i].R = 0; + for (i = 0; i < P.NC; i++) Cells.at(i).L = Cells.at(i).I = Cells.at(i).R = 0; } t = ((double)m) / ((double)P.Nplace[tp]); fprintf(stderr, "Adjusting place weights (Capacity=%i Demand=%i Av place size=%lg)\n", m, cnt, t); @@ -2013,20 +2006,20 @@ void AssignPeopleToPlaces(void) } if (P.PlaceTypeNearestNeighb[tp] == 0) { - for (i = 0; i < P.NC; i++) Cells[i].S = 0; + for (i = 0; i < P.NC; i++) Cells.at(i).S = 0; for (j = 0; j < P.Nplace[tp]; j++) { i = P.nch * ((int)(Places[tp][j].loc_x / P.cwidth)) + ((int)(Places[tp][j].loc_y / P.cheight)); - Cells[i].S += (int)Places[tp][j].treat_end_time; + Cells.at(i).S += (int)Places[tp][j].treat_end_time; } for (i = 0; i < P.NC; i++) { - if (Cells[i].S > Cells[i].cumTC) + if (Cells.at(i).S > Cells.at(i).cumTC) { - free(Cells[i].susceptible); - if (!(Cells[i].susceptible = (int*)calloc(Cells[i].S, sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); + free(Cells.at(i).susceptible); + if (!(Cells.at(i).susceptible = (int*)calloc(Cells.at(i).S, sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); } - Cells[i].S = 0; + Cells.at(i).S = 0; } for (j = 0; j < P.Nplace[tp]; j++) { @@ -2034,8 +2027,8 @@ void AssignPeopleToPlaces(void) k = (int)Places[tp][j].treat_end_time; for (j2 = 0; j2 < k; j2++) { - Cells[i].susceptible[Cells[i].S] = j; - Cells[i].S++; + Cells.at(i).susceptible[Cells.at(i).S] = j; + Cells.at(i).S++; } } } @@ -2215,8 +2208,8 @@ void AssignPeopleToPlaces(void) do { s = ranf(); - l = Cells[i].InvCDF[(int)floor(s * 1024)]; - while (Cells[i].cum_trans[l] < s) l++; + l = Cells.at(i).InvCDF[(int)floor(s * 1024)]; + while (Cells.at(i).cum_trans[l] < s) l++; ct = CellLookup[l]; m = (int)(ranf() * ((double)ct->S)); j = -1; @@ -2236,7 +2229,7 @@ void AssignPeopleToPlaces(void) ERR_CRITICAL("Out of bounds place link\n"); } t = dist2_raw(Households[Hosts[k].hh].loc_x, Households[Hosts[k].hh].loc_y, Places[tp][j].loc_x, Places[tp][j].loc_y); - s = ((double)ct->S) / ((double)ct->S0) * numKernel(t) / Cells[i].max_trans[l]; + s = ((double)ct->S) / ((double)ct->S0) * numKernel(t) / Cells.at(i).max_trans[l]; if ((P.DoAdUnits) && (P.InhibitInterAdunitPlaceAssignment[tp] > 0)) { if (Mcells[Hosts[k].mcell].adunit != Mcells[Places[tp][j].mcell].adunit) s *= (1 - P.InhibitInterAdunitPlaceAssignment[tp]); @@ -2277,11 +2270,11 @@ void AssignPeopleToPlaces(void) } for (i = 0; i < P.NC; i++) { - Cells[i].n = Cells[i].cumTC; - Cells[i].cumTC = 0; - Cells[i].S = Cells[i].I = Cells[i].L = Cells[i].R = 0; - free(Cells[i].susceptible); - Cells[i].susceptible = Cells[i].infected; + Cells.at(i).n = Cells.at(i).cumTC; + Cells.at(i).cumTC = 0; + Cells.at(i).S = Cells.at(i).I = Cells.at(i).L = Cells.at(i).R = 0; + free(Cells.at(i).susceptible); + Cells.at(i).susceptible = Cells.at(i).infected; } P.KernelScale = P.MoveKernelScale; P.KernelShape = P.MoveKernelShape; From 58a5b772c3318d6251e202705283a15cc8bc3cf4 Mon Sep 17 00:00:00 2001 From: chris-j-tang Date: Sat, 9 May 2020 15:54:08 -0400 Subject: [PATCH 2/7] Finished all Cell & Mcell -> vector --- src/SetupModel.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp index 37e2ccb43..2c9edf5e1 100644 --- a/src/SetupModel.cpp +++ b/src/SetupModel.cpp @@ -1003,10 +1003,11 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) if (!(mcl = (int*)malloc(P.PopSize * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); State.CellSuscMemberArray = mcl; i2 = k = 0; + auto Cells_front = &Cells.front(); for (j = 0; j < P.NC; j++) if (Cells[j].n > 0) { - CellLookup[i2++] = Cells + j; + CellLookup[i2++] = Cells_front + j; Cells[j].susceptible = mcl + k; k += Cells[j].n; } @@ -1031,11 +1032,12 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) for (j = 0; j < Cells.at(i).n; j++) Cells.at(i).members[j] = -1; } fprintf(stderr, "Cells assigned\n"); + auto Mcells_front = &Mcells.front(); for (i = 0; i <= MAX_HOUSEHOLD_SIZE; i++) denom_household[i] = 0; P.NH = 0; for (i = j2 = 0; j2 < P.NMCP; j2++) { - j = (int)(McellLookup[j2] - Mcells); + j = (int)(McellLookup[j2] - Mcells_front); l = ((j / P.nmch) / P.NMCL) * P.nch + ((j % P.nmch) / P.NMCL); ad = ((P.DoAdunitDemog) && (P.DoAdUnits)) ? Mcells[j].adunit : 0; for (k = 0; k < Mcells[j].n;) @@ -1068,9 +1070,10 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) if (P.DoHouseholds) fprintf(stderr, "Household sizes assigned to %i people\n", i); #pragma omp parallel for private(tn,j2,j,i,k,x,y,xh,yh,i2,m) schedule(static,1) for (tn = 0; tn < P.NumThreads; tn++) + auto Mcells_front = &Mcells.front(); for (j2 = tn; j2 < P.NMCP; j2 += P.NumThreads) { - j = (int)(McellLookup[j2] - Mcells); + j = (int)(McellLookup[j2] - Mcells_front); x = (double)(j / P.nmch); y = (double)(j % P.nmch); i = Mcells[j].members[0]; From 12c41c694acc77c620128d432e2f2330793d239c Mon Sep 17 00:00:00 2001 From: chris-j-tang Date: Sat, 9 May 2020 15:59:15 -0400 Subject: [PATCH 3/7] Cells[...] to Cells.at(...) in SetUpModel.cpp --- src/SetupModel.cpp | 44 ++++++++++++++++++++++---------------------- 1 file changed, 22 insertions(+), 22 deletions(-) diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp index 2c9edf5e1..9925da3a6 100644 --- a/src/SetupModel.cpp +++ b/src/SetupModel.cpp @@ -1005,11 +1005,11 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) i2 = k = 0; auto Cells_front = &Cells.front(); for (j = 0; j < P.NC; j++) - if (Cells[j].n > 0) + if (Cells.at(j).n > 0) { CellLookup[i2++] = Cells_front + j; - Cells[j].susceptible = mcl + k; - k += Cells[j].n; + Cells.at(j).susceptible = mcl + k; + k += Cells.at(j).n; } if (i2 > P.NCP) fprintf(stderr, "######## Over-run on CellLookup array NCP=%i i2=%i ###########\n", P.NCP, i2); i2 = 0; @@ -1054,8 +1054,8 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) // fprintf(stderr,"%i ",i+i2); Hosts[i + i2].listpos = m; //used temporarily to store household size Mcells[j].members[k + i2] = i + i2; - Cells[l].susceptible[Cells[l].cumTC] = i + i2; - Cells[l].members[Cells[l].cumTC++] = i + i2; + Cells.at(l).susceptible[Cells.at(l).cumTC] = i + i2; + Cells.at(l).members[Cells.at(l).cumTC++] = i + i2; Hosts[i + i2].pcell = l; Hosts[i + i2].mcell = j; Hosts[i + i2].hh = P.NH; @@ -1385,7 +1385,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) } l = 0; for (j = 0; j < P.NC; j++) - if (l < Cells[j].n) l = Cells[j].n; + if (l < Cells.at(j).n) l = Cells.at(j).n; if (!(SamplingQueue = (int**)malloc(P.NumThreads * sizeof(int*)))) ERR_CRITICAL("Unable to allocate state storage\n"); P.InfQueuePeakLength = P.PopSize / P.NumThreads / INF_QUEUE_SCALE; #pragma omp parallel for private(i,k) schedule(static,1) @@ -1922,53 +1922,53 @@ void AssignPeopleToPlaces(void) { j = (int)(Places[tp][i].loc_x / P.cwidth); k = j * P.nch + ((int)(Places[tp][i].loc_y / P.cheight)); - Cells[k].I += (int)Places[tp][i].treat_end_time; + Cells.at(k).I += (int)Places[tp][i].treat_end_time; } for (k = 0; k < P.NC; k++) { i = k % P.nch; j = k / P.nch; - f2 = Cells[k].I; f3 = Cells[k].S; + f2 = Cells.at(k).I; f3 = Cells.at(k).S; if ((i > 0) && (j > 0)) { - f2 += Cells[(j - 1) * P.nch + (i - 1)].I; f3 += Cells[(j - 1) * P.nch + (i - 1)].S; + f2 += Cells.at((j - 1) * P.nch + (i - 1)).I; f3 += Cells.at((j - 1) * P.nch + (i - 1)).S; } if (i > 0) { - f2 += Cells[j * P.nch + (i - 1)].I; f3 += Cells[j * P.nch + (i - 1)].S; + f2 += Cells.at(j * P.nch + (i - 1)).I; f3 += Cells.at(j * P.nch + (i - 1)).S; } if ((i > 0) && (j < P.ncw - 1)) { - f2 += Cells[(j + 1) * P.nch + (i - 1)].I; f3 += Cells[(j + 1) * P.nch + (i - 1)].S; + f2 += Cells.at((j + 1) * P.nch + (i - 1)).I; f3 += Cells.at((j + 1) * P.nch + (i - 1)).S; } if (j > 0) { - f2 += Cells[(j - 1) * P.nch + i].I; f3 += Cells[(j - 1) * P.nch + i].S; + f2 += Cells.at((j - 1) * P.nch + i).I; f3 += Cells.at((j - 1) * P.nch + i).S; } if (j < P.ncw - 1) { - f2 += Cells[(j + 1) * P.nch + i].I; f3 += Cells[(j + 1) * P.nch + i].S; + f2 += Cells.at((j + 1) * P.nch + i).I; f3 += Cells.at((j + 1) * P.nch + i).S; } if ((i < P.nch - 1) && (j > 0)) { - f2 += Cells[(j - 1) * P.nch + (i + 1)].I; f3 += Cells[(j - 1) * P.nch + (i + 1)].S; + f2 += Cells.at((j - 1) * P.nch + (i + 1)).I; f3 += Cells.at((j - 1) * P.nch + (i + 1)).S; } if (i < P.nch - 1) { - f2 += Cells[j * P.nch + (i + 1)].I; f3 += Cells[j * P.nch + (i + 1)].S; + f2 += Cells.at(j * P.nch + (i + 1)).I; f3 += Cells.at(j * P.nch + (i + 1)).S; } if ((i < P.nch - 1) && (j < P.ncw - 1)) { - f2 += Cells[(j + 1) * P.nch + (i + 1)].I; f3 += Cells[(j + 1) * P.nch + (i + 1)].S; + f2 += Cells.at((j + 1) * P.nch + (i + 1)).I; f3 += Cells.at((j + 1) * P.nch + (i + 1)).S; } - Cells[k].L = f3; Cells[k].R = f2; + Cells.at(k).L = f3; Cells.at(k).R = f2; } m = f2 = f3 = f4 = 0; for (k = 0; k < P.NC; k++) - if ((Cells[k].S > 0) && (Cells[k].I == 0)) + if ((Cells.at(k).S > 0) && (Cells.at(k).I == 0)) { - f2 += Cells[k].S; f3++; - if (Cells[k].R == 0) f4 += Cells[k].S; + f2 += Cells.at(k).S; f3++; + if (Cells.at(k).R == 0) f4 += Cells.at(k).S; } fprintf(stderr, "Demand in cells with no places=%i in %i cells\nDemand in cells with no places <=1 cell away=%i\n", f2, f3, f4); for (i = 0; i < P.Nplace[tp]; i++) @@ -1976,9 +1976,9 @@ void AssignPeopleToPlaces(void) { j = (int)(Places[tp][i].loc_x / P.cwidth); k = j * P.nch + ((int)(Places[tp][i].loc_y / P.cheight)); - if ((Cells[k].L > 0) && (Cells[k].R > 0)) + if ((Cells.at(k).L > 0) && (Cells.at(k).R > 0)) { - s = ((double)Cells[k].L) / ((double)Cells[k].R); + s = ((double)Cells.at(k).L) / ((double)Cells.at(k).R); Places[tp][i].treat_end_time = (unsigned short)ceil(Places[tp][i].treat_end_time * s); } m += ((int)Places[tp][i].treat_end_time); From c4665495126b5928261d36b3615d18624b04707c Mon Sep 17 00:00:00 2001 From: chris-j-tang Date: Sat, 9 May 2020 16:06:16 -0400 Subject: [PATCH 4/7] Keep Places as place** (for now) --- src/CovidSim.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index 1f1698826..461b9fa61 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -75,7 +75,7 @@ std::vector Cells(P.NC); // Cells.at(i) is the i'th cell std::vector CellLookup(P.NCP); // CellLookup[i] is a pointer to the i'th populated cell std::vector Mcells(P.NMC); std::vector McellLookup(P.NMCP); -std::vector> Places; +place** Places; adminunit AdUnits[MAX_ADUNITS]; //// Time Series defs: //// TimeSeries is an array of type results, used to store (unsurprisingly) a time series of every quantity in results. Mostly used in RecordSample. From b60635a84c4c3ca1365c00360c58b38b95f3cc70 Mon Sep 17 00:00:00 2001 From: chris-j-tang Date: Sat, 9 May 2020 16:06:52 -0400 Subject: [PATCH 5/7] Model.h's Places is correct type --- src/Model.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Model.h b/src/Model.h index fcc98768a..892158974 100644 --- a/src/Model.h +++ b/src/Model.h @@ -340,7 +340,7 @@ extern std::vector Cells; extern std::vector CellLookup; extern std::vector Mcells; extern std::vector McellLookup; -extern std::vector> Places; +extern place** Places; extern adminunit AdUnits[MAX_ADUNITS]; //// Time Series defs: From e61472a904fed564264381d0cafae3ec60a9da3f Mon Sep 17 00:00:00 2001 From: chris-j-tang Date: Sat, 9 May 2020 16:17:21 -0400 Subject: [PATCH 6/7] Sweep.cpp has correct Mcells and Cells --- src/Sweep.cpp | 244 +++++++++++++++++++++++++------------------------- 1 file changed, 123 insertions(+), 121 deletions(-) diff --git a/src/Sweep.cpp b/src/Sweep.cpp index d707261f2..ec27adc6c 100644 --- a/src/Sweep.cpp +++ b/src/Sweep.cpp @@ -100,8 +100,8 @@ void TravelDepartSweep(double t) l = Airports[i].Inv_DestMcells[(int)floor(s * 1024)]; while (Airports[i].DestMcells[l].prob < s) l++; l = Airports[i].DestMcells[l].id; - k = (int)(ranf_mt(tn) * ((double)Mcells[l].n)); - i2 = Mcells[l].members[k]; + k = (int)(ranf_mt(tn) * ((double)Mcells.at(l).n)); + i2 = Mcells.at(l).members[k]; if ((abs(Hosts[i2].inf) < InfStat_InfectiousAsymptomaticNotCase) && (Hosts[i2].inf != InfStat_Case)) { d2 = HOST_AGE_GROUP(i2); @@ -190,7 +190,7 @@ void TravelDepartSweep(double t) for (i = tn; i < P.Nplace[P.HotelPlaceType]; i += P.NumThreads) { c = ((int)(Places[P.HotelPlaceType][i].loc_x / P.cwidth)) * P.nch + ((int)(Places[P.HotelPlaceType][i].loc_y / P.cheight)); - n = (int)ignpoi_mt(nl * Cells[c].tot_prob, tn); + n = (int)ignpoi_mt(nl * Cells.at(c).tot_prob, tn); if (Places[P.HotelPlaceType][i].n + n > mps) { nsk += (Places[P.HotelPlaceType][i].n + n - mps); @@ -202,8 +202,8 @@ void TravelDepartSweep(double t) { f = 0; s = ranf_mt(tn); - l = Cells[c].InvCDF[(int)floor(s * 1024)]; - while (Cells[c].cum_trans[l] < s) l++; + l = Cells.at(c).InvCDF[(int)floor(s * 1024)]; + while (Cells.at(c).cum_trans[l] < s) l++; ct = CellLookup[l]; m = (int)(ranf_mt(tn) * ((double)ct->S0)); if (m < (ct->S + ct->L)) @@ -232,7 +232,7 @@ void TravelDepartSweep(double t) } if (f2) { - s = numKernel(s2) / Cells[c].max_trans[l]; + s = numKernel(s2) / Cells.at(c).max_trans[l]; if (ranf_mt(tn) >= s) { #pragma omp critical @@ -304,6 +304,7 @@ void InfectSweep(double t, int run) //added run number as argument in order to r for (tn = 0; tn < P.NumThreads; tn++) for (b = tn; b < P.NCP; b += P.NumThreads) //// loop over (in parallel) all populated cells. Loop 1) { + auto Mcells_front = &Mcells.front(); c = CellLookup[b]; s5 = 0; ///// spatial infectiousness summed over all infectious people in loop below for (j = 0; j < c->I; j++) //// Loop over all infectious people c->I in cell c. Loop 1a) @@ -314,8 +315,8 @@ void InfectSweep(double t, int run) //added run number as argument in order to r si = Hosts + ci; //calculate flag for digital contact tracing here at the beginning for each individual - fct = ((P.DoDigitalContactTracing) && (t >= AdUnits[Mcells[si->mcell].adunit].DigitalContactTracingTimeStart) - && (t < AdUnits[Mcells[si->mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ci].digitalContactTracingUser == 1)); // && (ts <= (Hosts[ci].detected_time + P.usCaseIsolationDelay))); + fct = ((P.DoDigitalContactTracing) && (t >= AdUnits[Mcells.at(si->mcell).adunit].DigitalContactTracingTimeStart) + && (t < AdUnits[Mcells.at(si->mcell).adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ci].digitalContactTracingUser == 1)); // && (ts <= (Hosts[ci].detected_time + P.usCaseIsolationDelay))); //// Household FOI component if (hbeta > 0) @@ -356,14 +357,14 @@ void InfectSweep(double t, int run) //added run number as argument in order to r { if (!HOST_ABSENT(ci)) { - mi = Mcells + si->mcell; + mi = Mcells_front + si->mcell; for (k = 0; k < P.PlaceTypeNum; k++) //// loop over all place types { l = si->PlaceLinks[k]; if (l >= 0) //// l>=0 means if place type k is relevant to person si. (Now allowing for partial attendance). { s3 = fp * seasonality * CalcPlaceInf(ci, k, ts); - mp = Mcells + Places[k][l].mcell; + mp = Mcells_front + Places[k][l].mcell; if (bm) { if ((dist2_raw(Households[si->hh].loc_x, Households[si->hh].loc_y, @@ -414,7 +415,7 @@ void InfectSweep(double t, int run) //added run number as argument in order to r if ((Hosts[ci].ncontacts < P.MaxDigitalContactsToTrace) && (ranf_mt(tn) = AdUnits[Mcells[si->mcell].adunit].DigitalContactTracingTimeStart) - && (t < AdUnits[Mcells[si->mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ci].digitalContactTracingUser == 1)); // && (ts <= (Hosts[ci].detected_time + P.usCaseIsolationDelay))); + fct = ((P.DoDigitalContactTracing) && (t >= AdUnits[Mcells.at(si->mcell).adunit].DigitalContactTracingTimeStart) + && (t < AdUnits[Mcells.at(si->mcell).adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ci].digitalContactTracingUser == 1)); // && (ts <= (Hosts[ci].detected_time + P.usCaseIsolationDelay))); //// decide on infectee @@ -631,8 +632,8 @@ void InfectSweep(double t, int run) //added run number as argument in order to r { if ((!Hosts[i3].Travelling) && ((c != ct) || (Hosts[i3].hh != si->hh))) //// if potential infectee not travelling, and either is not part of cell c or doesn't share a household with infector. { - mi = Mcells + si->mcell; - mt = Mcells + Hosts[i3].mcell; + mi = Mcells_front + si->mcell; + mt = Mcells_front + Hosts[i3].mcell; s = CalcSpatialSusc(i3, ts, ci, tn); //so this person is a contact - but might not be infected. if we are doing digital contact tracing, we want to add the person to the contacts list, if both are users @@ -644,7 +645,7 @@ void InfectSweep(double t, int run) //added run number as argument in order to r if ((Hosts[ci].ncontacts 0) && (ranf_mt(tn) < P.FalsePositiveRate)) @@ -810,7 +811,7 @@ void IncubRecoverySweep(double t, int run) } //once host recovers, will no longer make contacts for contact tracing - if we are doing contact tracing and case was infectious when contact tracing was active, increment state vector - if ((P.DoDigitalContactTracing) && (Hosts[ci].latent_time>= AdUnits[Mcells[Hosts[ci].mcell].adunit].DigitalContactTracingTimeStart) && (Hosts[ci].recovery_or_death_time < AdUnits[Mcells[Hosts[ci].mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ci].digitalContactTracingUser == 1) && (P.OutputDigitalContactDist)) + if ((P.DoDigitalContactTracing) && (Hosts[ci].latent_time>= AdUnits[Mcells.at(Hosts[ci].mcell).adunit].DigitalContactTracingTimeStart) && (Hosts[ci].recovery_or_death_time < AdUnits[Mcells.at(Hosts[ci].mcell).adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[ci].digitalContactTracingUser == 1) && (P.OutputDigitalContactDist)) { if (Hosts[ci].ncontacts > MAX_CONTACTS) Hosts[ci].ncontacts = MAX_CONTACTS; //increment bin in State corresponding to this number of contacts @@ -1274,8 +1275,9 @@ int TreatSweep(double t) for (tn = 0; tn < P.NumThreads; tn++) for (bs = tn; bs < P.NMCP; bs += P.NumThreads) //// loop over populated microcells { - b = (int)(McellLookup[bs] - Mcells); //// microcell number - adi = (P.DoAdUnits) ? Mcells[b].adunit : -1; + auto Mcells_front = &Mcells.front(); + b = (int)(McellLookup[bs] - &Mcells.front()); //// microcell number + adi = (P.DoAdUnits) ? Mcells.at(b).adunit : -1; ad = (P.DoAdUnits) ? AdUnits[adi].id : 0; //// Code block goes through various types of treatments/interventions (vaccination/movement restrictions etc.), @@ -1289,20 +1291,20 @@ int TreatSweep(double t) //// **** //// **** //// **** //// **** TREATMENT //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** - if ((Mcells[b].treat == 2) && (ts >= Mcells[b].treat_end_time)) + if ((Mcells.at(b).treat == 2) && (ts >= Mcells.at(b).treat_end_time)) { f = 1; - Mcells[b].treat = 0; + Mcells.at(b).treat = 0; } - if ((Mcells[b].treat == 1) && (ts >= Mcells[b].treat_start_time)) + if ((Mcells.at(b).treat == 1) && (ts >= Mcells.at(b).treat_start_time)) { f = 1; - Mcells[b].treat = 2; - Mcells[b].treat_trig = 0; - Mcells[b].treat_end_time = tstf; - for (i = 0; i < Mcells[b].n; i++) + Mcells.at(b).treat = 2; + Mcells.at(b).treat_trig = 0; + Mcells.at(b).treat_end_time = tstf; + for (i = 0; i < Mcells.at(b).n; i++) { - l = Mcells[b].members[i]; + l = Mcells.at(b).members[i]; if ((!HOST_TO_BE_TREATED(l)) && ((P.TreatPropRadial == 1) || (ranf_mt(tn) < P.TreatPropRadial))) DoProphNoDelay(l, ts, tn, 1); } @@ -1316,10 +1318,10 @@ int TreatSweep(double t) } else { - trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.TreatCellIncThresh)) / P.IncThreshPop)) : (int)P.TreatCellIncThresh; - f2 = (Mcells[b].treat_trig >= trig_thresh); + trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells.at(b).n * P.TreatCellIncThresh)) / P.IncThreshPop)) : (int)P.TreatCellIncThresh; + f2 = (Mcells.at(b).treat_trig >= trig_thresh); } - if ((t >= P.TreatTimeStart) && (Mcells[b].treat == 0) && (f2) && (P.TreatRadius2 > 0)) + if ((t >= P.TreatTimeStart) && (Mcells.at(b).treat == 0) && (f2) && (P.TreatRadius2 > 0)) { minx = (b / P.nmch); miny = (b % P.nmch); k = b; @@ -1334,17 +1336,17 @@ int TreatSweep(double t) if ((minx >= 0) && (minx < P.nmcw) && (miny >= 0) && (miny < P.nmch)) { if (P.TreatByAdminUnit) - f4 = (AdUnits[Mcells[k].adunit].id / P.TreatAdminUnitDivisor == ad2); + f4 = (AdUnits[Mcells.at(k).adunit].id / P.TreatAdminUnitDivisor == ad2); else - f4 = ((r = dist2_mm(Mcells + b, Mcells + k)) < P.TreatRadius2); + f4 = ((r = dist2_mm(Mcells_front + b, Mcells_front + k)) < P.TreatRadius2); if (f4) { f = f2 = 1; - if ((Mcells[k].n > 0) && (Mcells[k].treat == 0)) + if ((Mcells.at(k).n > 0) && (Mcells.at(k).treat == 0)) { - Mcells[k].treat_start_time = tstb; - Mcells[k].treat = 1; - maxx += Mcells[k].n; + Mcells.at(k).treat_start_time = tstb; + Mcells.at(k).treat = 1; + maxx += Mcells.at(k).n; } } } @@ -1376,15 +1378,15 @@ int TreatSweep(double t) //// vaccinates proportion VaccProp of people in microcell (or at least adds them to geovacc_queue). - if ((Mcells[b].vacc == 1) && (ts >= Mcells[b].vacc_start_time)) + if ((Mcells.at(b).vacc == 1) && (ts >= Mcells.at(b).vacc_start_time)) { f = 1; - Mcells[b].vacc_trig = 0; - //if(State.cumVG+P.NumThreads*Mcells[b].n= trig_thresh); + trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells.at(b).n * P.VaccCellIncThresh)) / P.IncThreshPop)) : (int)P.VaccCellIncThresh; + f2 = (Mcells.at(b).treat_trig >= trig_thresh); } - if ((!P.DoMassVacc) && (P.VaccRadius2 > 0) && (t >= P.VaccTimeStartGeo) && (Mcells[b].vacc == 0) && (f2)) //changed from VaccTimeStart to VaccTimeStarGeo + if ((!P.DoMassVacc) && (P.VaccRadius2 > 0) && (t >= P.VaccTimeStartGeo) && (Mcells.at(b).vacc == 0) && (f2)) //changed from VaccTimeStart to VaccTimeStarGeo { minx = (b / P.nmch); miny = (b % P.nmch); k = b; @@ -1422,20 +1424,20 @@ int TreatSweep(double t) { if (P.VaccByAdminUnit) { - f4 = (AdUnits[Mcells[k].adunit].id / P.VaccAdminUnitDivisor == ad2); + f4 = (AdUnits[Mcells.at(k).adunit].id / P.VaccAdminUnitDivisor == ad2); r = 1e20; } else - f4 = ((r = dist2_mm(Mcells + b, Mcells + k)) < P.VaccRadius2); + f4 = ((r = dist2_mm(Mcells_front + b, Mcells_front + k)) < P.VaccRadius2); if (f4) { f = f2 = 1; if (r < P.VaccMinRadius2) - Mcells[k].vacc = 3; - else if ((Mcells[k].n > 0) && (Mcells[k].vacc == 0)) + Mcells.at(k).vacc = 3; + else if ((Mcells.at(k).n > 0) && (Mcells.at(k).vacc == 0)) { - Mcells[k].vacc_start_time = tsvb; - Mcells[k].vacc = 1; + Mcells.at(k).vacc_start_time = tsvb; + Mcells.at(k).vacc = 1; } } } @@ -1475,26 +1477,26 @@ int TreatSweep(double t) } else { - trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.PlaceCloseCellIncStopThresh)) / P.IncThreshPop)) : P.PlaceCloseCellIncStopThresh; - f2 = (Mcells[b].treat_trig < trig_thresh); + trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells.at(b).n * P.PlaceCloseCellIncStopThresh)) / P.IncThreshPop)) : P.PlaceCloseCellIncStopThresh; + f2 = (Mcells.at(b).treat_trig < trig_thresh); } - if ((Mcells[b].placeclose == 2) && ((f2) || (ts >= Mcells[b].place_end_time))) //// if place closure has started, the places in this microcell are closed, and either stop threshold has been reached or place_end_time has passed, go through block + if ((Mcells.at(b).placeclose == 2) && ((f2) || (ts >= Mcells.at(b).place_end_time))) //// if place closure has started, the places in this microcell are closed, and either stop threshold has been reached or place_end_time has passed, go through block { f = 1; - Mcells[b].placeclose = P.DoPlaceCloseOnceOnly; - Mcells[b].place_end_time = ts; - Mcells[b].place_trig = 0; + Mcells.at(b).placeclose = P.DoPlaceCloseOnceOnly; + Mcells.at(b).place_end_time = ts; + Mcells.at(b).place_trig = 0; if (f2) { for (j2 = 0; j2 < P.PlaceTypeNum; j2++) if (j2 != P.HotelPlaceType) - for (i2 = 0; i2 < Mcells[b].np[j2]; i2++) - DoPlaceOpen(j2, Mcells[b].places[j2][i2], ts, tn); + for (i2 = 0; i2 < Mcells.at(b).np[j2]; i2++) + DoPlaceOpen(j2, Mcells.at(b).places[j2][i2], ts, tn); } } - if ((P.DoPlaces) && (t >= P.PlaceCloseTimeStart) && (Mcells[b].placeclose == 0)) + if ((P.DoPlaces) && (t >= P.PlaceCloseTimeStart) && (Mcells.at(b).placeclose == 0)) { ///// note that here f2 bool asks whether trigger has exceeded threshold in order to close places for first time.A few blocks up meaning was almost the opposite: asking whether trigger lower than stop threshold. @@ -1509,34 +1511,34 @@ int TreatSweep(double t) } else { - trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.PlaceCloseCellIncThresh)) / P.IncThreshPop)) : P.PlaceCloseCellIncThresh; - f2 = (Mcells[b].treat_trig >= trig_thresh); + trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells.at(b).n * P.PlaceCloseCellIncThresh)) / P.IncThreshPop)) : P.PlaceCloseCellIncThresh; + f2 = (Mcells.at(b).treat_trig >= trig_thresh); } - if (((P.PlaceCloseByAdminUnit) && (AdUnits[Mcells[b].adunit].place_close_trig < USHRT_MAX - 1) - && (((double)AdUnits[Mcells[b].adunit].place_close_trig) / ((double)AdUnits[Mcells[b].adunit].NP) > P.PlaceCloseAdunitPropThresh)) + if (((P.PlaceCloseByAdminUnit) && (AdUnits[Mcells.at(b).adunit].place_close_trig < USHRT_MAX - 1) + && (((double)AdUnits[Mcells.at(b).adunit].place_close_trig) / ((double)AdUnits[Mcells.at(b).adunit].NP) > P.PlaceCloseAdunitPropThresh)) || ((!P.PlaceCloseByAdminUnit) && (f2))) { - // if(P.PlaceCloseByAdminUnit) AdUnits[Mcells[b].adunit].place_close_trig=USHRT_MAX-1; // This means schools only close once + // if(P.PlaceCloseByAdminUnit) AdUnits[Mcells.at(b).adunit].place_close_trig=USHRT_MAX-1; // This means schools only close once interventionFlag = 1; - if ((P.DoInterventionDelaysByAdUnit)&&((t <= AdUnits[Mcells[b].adunit].PlaceCloseTimeStart) || (t >= (AdUnits[Mcells[b].adunit].PlaceCloseTimeStart + AdUnits[Mcells[b].adunit].PlaceCloseDuration)))) + if ((P.DoInterventionDelaysByAdUnit)&&((t <= AdUnits[Mcells.at(b).adunit].PlaceCloseTimeStart) || (t >= (AdUnits[Mcells.at(b).adunit].PlaceCloseTimeStart + AdUnits[Mcells.at(b).adunit].PlaceCloseDuration)))) interventionFlag = 0; if ((interventionFlag == 1)&&((!P.PlaceCloseByAdminUnit) || (ad > 0))) { ad2 = ad / P.PlaceCloseAdminUnitDivisor; - if ((Mcells[b].n > 0) && (Mcells[b].placeclose == 0)) + if ((Mcells.at(b).n > 0) && (Mcells.at(b).placeclose == 0)) { //if doing intervention delays and durations by admin unit based on global triggers if (P.DoInterventionDelaysByAdUnit) - Mcells[b].place_end_time = (unsigned short int) ceil(P.TimeStepsPerDay * (t + P.PlaceCloseDelayMean + AdUnits[Mcells[b].adunit].PlaceCloseDuration)); + Mcells.at(b).place_end_time = (unsigned short int) ceil(P.TimeStepsPerDay * (t + P.PlaceCloseDelayMean + AdUnits[Mcells.at(b).adunit].PlaceCloseDuration)); else - Mcells[b].place_end_time = tspf; - Mcells[b].place_trig = 0; - Mcells[b].placeclose = 2; + Mcells.at(b).place_end_time = tspf; + Mcells.at(b).place_trig = 0; + Mcells.at(b).placeclose = 2; for (j2 = 0; j2 < P.PlaceTypeNum; j2++) if (j2 != P.HotelPlaceType) - for (i2 = 0; i2 < Mcells[b].np[j2]; i2++) - DoPlaceClose(j2, Mcells[b].places[j2][i2], ts, tn, 1); + for (i2 = 0; i2 < Mcells.at(b).np[j2]; i2++) + DoPlaceClose(j2, Mcells.at(b).places[j2][i2], ts, tn, 1); } } } @@ -1547,17 +1549,17 @@ int TreatSweep(double t) //// **** //// **** //// **** //// **** MOVEMENT RESTRICTIONS //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** - if ((Mcells[b].moverest == 2) && (ts >= Mcells[b].move_end_time)) + if ((Mcells.at(b).moverest == 2) && (ts >= Mcells.at(b).move_end_time)) { f = 1; - Mcells[b].moverest = P.DoMoveRestrOnceOnly; + Mcells.at(b).moverest = P.DoMoveRestrOnceOnly; } - if ((Mcells[b].moverest == 1) && (ts >= Mcells[b].move_start_time)) + if ((Mcells.at(b).moverest == 1) && (ts >= Mcells.at(b).move_start_time)) { f = 1; - Mcells[b].moverest = 2; - Mcells[b].move_trig = 0; - Mcells[b].move_end_time = tsmf; + Mcells.at(b).moverest = 2; + Mcells.at(b).move_trig = 0; + Mcells.at(b).move_end_time = tsmf; } if (P.DoGlobalTriggers) f2 = (global_trig >= P.MoveRestrCellIncThresh); @@ -1568,11 +1570,11 @@ int TreatSweep(double t) } else { - trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.MoveRestrCellIncThresh)) / P.IncThreshPop)) : P.MoveRestrCellIncThresh; - f2 = (Mcells[b].treat_trig >= trig_thresh); + trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells.at(b).n * P.MoveRestrCellIncThresh)) / P.IncThreshPop)) : P.MoveRestrCellIncThresh; + f2 = (Mcells.at(b).treat_trig >= trig_thresh); } - if ((t >= P.MoveRestrTimeStart) && (Mcells[b].moverest == 0) && (f2)) + if ((t >= P.MoveRestrTimeStart) && (Mcells.at(b).moverest == 0) && (f2)) { minx = (b / P.nmch); miny = (b % P.nmch); k = b; @@ -1586,16 +1588,16 @@ int TreatSweep(double t) if ((minx >= 0) && (minx < P.nmcw) && (miny >= 0) && (miny < P.nmch)) { if (P.MoveRestrByAdminUnit) - f4 = (AdUnits[Mcells[k].adunit].id / P.MoveRestrAdminUnitDivisor == ad2); + f4 = (AdUnits[Mcells.at(k).adunit].id / P.MoveRestrAdminUnitDivisor == ad2); else - f4 = ((r = dist2_mm(Mcells + b, Mcells + k)) < P.MoveRestrRadius2); + f4 = ((r = dist2_mm(Mcells_front + b, Mcells_front + k)) < P.MoveRestrRadius2); if (f4) { f = f2 = 1; - if ((Mcells[k].n > 0) && (Mcells[k].moverest == 0)) + if ((Mcells.at(k).n > 0) && (Mcells.at(k).moverest == 0)) { - Mcells[k].move_start_time = tsmb; - Mcells[k].moverest = 1; + Mcells.at(k).move_start_time = tsmb; + Mcells.at(k).moverest = 1; } } } @@ -1634,17 +1636,17 @@ int TreatSweep(double t) } else { - trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.SocDistCellIncStopThresh)) / P.IncThreshPop)) : P.SocDistCellIncStopThresh; - f2 = (Mcells[b].treat_trig < trig_thresh); + trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells.at(b).n * P.SocDistCellIncStopThresh)) / P.IncThreshPop)) : P.SocDistCellIncStopThresh; + f2 = (Mcells.at(b).treat_trig < trig_thresh); } //// if: policy of social distancing has started AND this microcell cell has been labelled to as undergoing social distancing, AND either trigger not reached (note definition of f2 changes in next few lines) or end time has passed. - if ((t >= P.SocDistTimeStart) && (Mcells[b].socdist == 2) && ((f2) || (ts >= Mcells[b].socdist_end_time))) + if ((t >= P.SocDistTimeStart) && (Mcells.at(b).socdist == 2) && ((f2) || (ts >= Mcells.at(b).socdist_end_time))) { f = 1; - Mcells[b].socdist = P.DoSocDistOnceOnly; //// i.e. if P.DoSocDistOnceOnly set to false, socdist set to 0 here, hence will be picked up upon some subsequent call to TreatSweep if required trigger passes threshold. - Mcells[b].socdist_trig = 0; //// reset trigger - Mcells[b].socdist_end_time = ts; //// record end time. + Mcells.at(b).socdist = P.DoSocDistOnceOnly; //// i.e. if P.DoSocDistOnceOnly set to false, socdist set to 0 here, hence will be picked up upon some subsequent call to TreatSweep if required trigger passes threshold. + Mcells.at(b).socdist_trig = 0; //// reset trigger + Mcells.at(b).socdist_end_time = ts; //// record end time. } if (P.DoGlobalTriggers) f2 = (global_trig >= P.SocDistCellIncThresh); @@ -1655,29 +1657,29 @@ int TreatSweep(double t) } else { - trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.SocDistCellIncThresh)) / P.IncThreshPop)) : P.SocDistCellIncThresh; - f2 = (Mcells[b].treat_trig >= trig_thresh); + trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells.at(b).n * P.SocDistCellIncThresh)) / P.IncThreshPop)) : P.SocDistCellIncThresh; + f2 = (Mcells.at(b).treat_trig >= trig_thresh); } - if ((t >= P.SocDistTimeStart) && (Mcells[b].socdist == 0) && (f2)) + if ((t >= P.SocDistTimeStart) && (Mcells.at(b).socdist == 0) && (f2)) { //some code to try and deal with intervention delays and durations by admin unit based on global triggers interventionFlag = 1; if (P.DoInterventionDelaysByAdUnit) - if ((t <= AdUnits[Mcells[b].adunit].SocialDistanceTimeStart) || - (t >= (AdUnits[Mcells[b].adunit].SocialDistanceTimeStart + AdUnits[Mcells[b].adunit].SocialDistanceDuration))) //// i.e. if outside window of social distancing for this admin unit. + if ((t <= AdUnits[Mcells.at(b).adunit].SocialDistanceTimeStart) || + (t >= (AdUnits[Mcells.at(b).adunit].SocialDistanceTimeStart + AdUnits[Mcells.at(b).adunit].SocialDistanceDuration))) //// i.e. if outside window of social distancing for this admin unit. interventionFlag = 0; if (interventionFlag == 1) - if ((Mcells[b].n > 0) && (Mcells[b].socdist == 0)) //// if microcell populated and not currently undergoing social distancing + if ((Mcells.at(b).n > 0) && (Mcells.at(b).socdist == 0)) //// if microcell populated and not currently undergoing social distancing { - Mcells[b].socdist = 2; //// update flag to denote that cell is undergoing social distancing - Mcells[b].socdist_trig = 0; /// reset trigger + Mcells.at(b).socdist = 2; //// update flag to denote that cell is undergoing social distancing + Mcells.at(b).socdist_trig = 0; /// reset trigger //// set (admin-specific) social distancing end time. if (P.DoInterventionDelaysByAdUnit) - Mcells[b].socdist_end_time = (unsigned short int) ceil(P.TimeStepsPerDay * (t + AdUnits[Mcells[b].adunit].SocialDistanceDuration)); + Mcells.at(b).socdist_end_time = (unsigned short int) ceil(P.TimeStepsPerDay * (t + AdUnits[Mcells.at(b).adunit].SocialDistanceDuration)); else - Mcells[b].socdist_end_time = tssdf; + Mcells.at(b).socdist_end_time = tssdf; } } @@ -1686,10 +1688,10 @@ int TreatSweep(double t) //// **** //// **** //// **** //// **** KEY-WORKER PROPHYLAXIS //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** - if ((Mcells[b].keyworkerproph == 2) && (ts >= Mcells[b].keyworkerproph_end_time)) + if ((Mcells.at(b).keyworkerproph == 2) && (ts >= Mcells.at(b).keyworkerproph_end_time)) { f = 1; - Mcells[b].keyworkerproph = P.DoKeyWorkerProphOnceOnly; + Mcells.at(b).keyworkerproph = P.DoKeyWorkerProphOnceOnly; } if (P.DoGlobalTriggers) f2 = (global_trig >= P.KeyWorkerProphCellIncThresh); @@ -1700,10 +1702,10 @@ int TreatSweep(double t) } else { - trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells[b].n * P.KeyWorkerProphCellIncThresh)) / P.IncThreshPop)) : P.KeyWorkerProphCellIncThresh; - f2 = (Mcells[b].treat_trig >= trig_thresh); + trig_thresh = (P.DoPerCapitaTriggers) ? ((int)ceil(((double)(Mcells.at(b).n * P.KeyWorkerProphCellIncThresh)) / P.IncThreshPop)) : P.KeyWorkerProphCellIncThresh; + f2 = (Mcells.at(b).treat_trig >= trig_thresh); } - if ((P.DoPlaces) && (t >= P.KeyWorkerProphTimeStart) && (Mcells[b].keyworkerproph == 0) && (f2)) + if ((P.DoPlaces) && (t >= P.KeyWorkerProphTimeStart) && (Mcells.at(b).keyworkerproph == 0) && (f2)) { minx = (b / P.nmch); miny = (b % P.nmch); k = b; @@ -1712,17 +1714,17 @@ int TreatSweep(double t) do { if ((minx >= 0) && (minx < P.nmcw) && (miny >= 0) && (miny < P.nmch)) - if (dist2_mm(Mcells + b, Mcells + k) < P.KeyWorkerProphRadius2) + if (dist2_mm(Mcells_front + b, Mcells_front + k) < P.KeyWorkerProphRadius2) { f = f2 = 1; - if ((Mcells[k].n > 0) && (Mcells[k].keyworkerproph == 0)) + if ((Mcells.at(k).n > 0) && (Mcells.at(k).keyworkerproph == 0)) { - Mcells[k].keyworkerproph = 2; - Mcells[k].keyworkerproph_trig = 0; - Mcells[k].keyworkerproph_end_time = tskwpf; - for (i2 = 0; i2 < Mcells[k].n; i2++) + Mcells.at(k).keyworkerproph = 2; + Mcells.at(k).keyworkerproph_trig = 0; + Mcells.at(k).keyworkerproph_end_time = tskwpf; + for (i2 = 0; i2 < Mcells.at(k).n; i2++) { - j2 = Mcells[k].members[i2]; + j2 = Mcells.at(k).members[i2]; if ((Hosts[j2].keyworker) && (!HOST_TO_BE_TREATED(j2))) DoProphNoDelay(j2, ts, tn, nckwp); } From f8305039bc35c07a981cc57d768b41e7c73fec17 Mon Sep 17 00:00:00 2001 From: chris-j-tang Date: Sat, 9 May 2020 16:29:30 -0400 Subject: [PATCH 7/7] Mcells to vector and .at() in SetupModel.cpp --- src/SetupModel.cpp | 64 +++++++++++++++++++++++----------------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp index 9925da3a6..460dbc87d 100644 --- a/src/SetupModel.cpp +++ b/src/SetupModel.cpp @@ -967,7 +967,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) if (Mcells[P.NMC - 1].n > 0) { P.NMCP++; - AdUnits[mcell_adunits[P.NMC - 1]].n += Mcells[P.NMC - 1].n; + AdUnits[mcell_adunits[P.NMC - 1]].n += Mcells.back().n; } t = 0.0; @@ -984,13 +984,13 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) for (m = 0; m < P.NMCL; m++) { j = k + m + l * P.nmch; - if (Mcells[j].n > 0) + if (Mcells.at(j).n > 0) { - Mcells[j].members = mcl + j2; - //if(!(Mcells[j].members=(int *) calloc(Mcells[j].n,sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); //replaced line above with this to ensure members don't get mixed across microcells + Mcells.at(j).members = mcl + j2; + //if(!(Mcells.at(j).members=(int *) calloc(Mcells.at(j).n,sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); //replaced line above with this to ensure members don't get mixed across microcells McellLookup[i2++] = &Mcells.front() + j; - Cells.at(i).n += Mcells[j].n; - j2 += Mcells[j].n; + Cells.at(i).n += Mcells.at(j).n; + j2 += Mcells.at(j).n; } } if (Cells.at(i).n > 0) P.NCP++; @@ -1007,7 +1007,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) for (j = 0; j < P.NC; j++) if (Cells.at(j).n > 0) { - CellLookup[i2++] = Cells_front + j; + CellLookup.at(i2++) = Cells_front + j; Cells.at(j).susceptible = mcl + k; k += Cells.at(j).n; } @@ -1018,7 +1018,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) fprintf(stderr, "sizeof(person)=%i\n", (int) sizeof(person)); for (i = 0; i < P.NCP; i++) { - cell *c = CellLookup[i]; + cell *c = CellLookup.at(i); if (c->n > 0) { if (!(c->InvCDF = (int*)malloc(1025 * sizeof(int)))) ERR_CRITICAL("Unable to allocate cell storage\n"); @@ -1039,21 +1039,21 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) { j = (int)(McellLookup[j2] - Mcells_front); l = ((j / P.nmch) / P.NMCL) * P.nch + ((j % P.nmch) / P.NMCL); - ad = ((P.DoAdunitDemog) && (P.DoAdUnits)) ? Mcells[j].adunit : 0; - for (k = 0; k < Mcells[j].n;) + ad = ((P.DoAdunitDemog) && (P.DoAdUnits)) ? Mcells.at(j).adunit : 0; + for (k = 0; k < Mcells.at(j).n;) { m = 1; if (P.DoHouseholds) { s = ranf_mt(0); - while ((s > P.HouseholdSizeDistrib[ad][m - 1]) && (k + m < Mcells[j].n) && (m < MAX_HOUSEHOLD_SIZE)) m++; + while ((s > P.HouseholdSizeDistrib[ad][m - 1]) && (k + m < Mcells.at(j).n) && (m < MAX_HOUSEHOLD_SIZE)) m++; } denom_household[m]++; for (i2 = 0; i2 < m; i2++) { // fprintf(stderr,"%i ",i+i2); Hosts[i + i2].listpos = m; //used temporarily to store household size - Mcells[j].members[k + i2] = i + i2; + Mcells.at(j).members[k + i2] = i + i2; Cells.at(l).susceptible[Cells.at(l).cumTC] = i + i2; Cells.at(l).members[Cells.at(l).cumTC++] = i + i2; Hosts[i + i2].pcell = l; @@ -1076,10 +1076,10 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) j = (int)(McellLookup[j2] - Mcells_front); x = (double)(j / P.nmch); y = (double)(j % P.nmch); - i = Mcells[j].members[0]; + i = Mcells.at(j).members[0]; if (j % 100 == 0) - fprintf(stderr, "%i=%i (%i %i) \r", j, Mcells[j].n, Mcells[j].adunit, (AdUnits[Mcells[j].adunit].id % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor); - for (k = 0; k < Mcells[j].n;) + fprintf(stderr, "%i=%i (%i %i) \r", j, Mcells.at(j).n, Mcells.at(j).adunit, (AdUnits[Mcells.at(j).adunit].id % P.AdunitLevel1Mask) / P.AdunitLevel1Divisor); + for (k = 0; k < Mcells.at(j).n;) { m = Hosts[i].listpos; xh = P.mcwidth * (ranf_mt(tn) + x); @@ -1120,7 +1120,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) AgeDistAd[i][j] = 0; for (i = 0; i < P.PopSize; i++) { - k = (P.DoAdunitDemog) ? Mcells[Hosts[i].mcell].adunit : 0; + k = (P.DoAdunitDemog) ? Mcells.at(Hosts[i].mcell).adunit : 0; AgeDistAd[k][HOST_AGE_GROUP(i)]++; } // normalize AgeDistAd[i][j], so it's the proportion of people in adunit i that are in age group j @@ -1166,7 +1166,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) for (tn = 0; tn < P.NumThreads; tn++) for (i = tn; i < P.PopSize; i += P.NumThreads) { - m = (P.DoAdunitDemog) ? Mcells[Hosts[i].mcell].adunit : 0; + m = (P.DoAdunitDemog) ? Mcells.at(Hosts[i].mcell).adunit : 0; j = HOST_AGE_GROUP(i); s = ranf_mt(tn); // probabilistic age adjustment by one age category (5 years) @@ -1204,9 +1204,9 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) double rand_r = 0.0; //added these variables so that if initial infection location is empty we can search the 10km neighbourhood to find a suitable cell double rand_theta = 0.0; int counter = 0; - if (Mcells[j].n < P.NumInitialInfections[0]) + if (Mcells.at(j).n < P.NumInitialInfections[0]) { - while (Mcells[j].n < P.NumInitialInfections[0] && counter < 100) + while (Mcells.at(j).n < P.NumInitialInfections[0] && counter < 100) { rand_r = ranf(); rand_theta = ranf(); rand_r = 0.083 * sqrt(rand_r); rand_theta = 2 * PI * rand_theta; //rand_r is multiplied by 0.083 as this is roughly equal to 10km in decimal degrees @@ -1221,7 +1221,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) P.LocationInitialInfection[0][1] = P.LocationInitialInfection[0][1] + rand_r * sin(rand_theta); } } - if (Mcells[j].n < P.NumInitialInfections[0]) + if (Mcells.at(j).n < P.NumInitialInfections[0]) ERR_CRITICAL("Too few people in seed microcell to start epidemic with required number of initial infectionz.\n"); } fprintf(stderr, "Checking cells...\n"); @@ -1292,7 +1292,7 @@ void SetupPopulation(char* DensityFile, char* SchoolFile, char* RegDemogFile) i = (int)(Places[j][P.Nplace[j]].loc_x / P.mcwidth); k = (int)(Places[j][P.Nplace[j]].loc_y / P.mcheight); j2 = i * P.nmch + k; - Mcells[j2].np[j]++; + Mcells.at(j2).np[j]++; Places[j][P.Nplace[j]].mcell = j2; P.Nplace[j]++; if (P.Nplace[j] % 1000 == 0) fprintf(stderr, "%i read \r", P.Nplace[j]); @@ -1554,13 +1554,13 @@ void SetupAirports(void) } l = 0; for (j = 0; j < Airports[i].num_mcell; j++) - l += Mcells[Airports[i].DestMcells[j].id].np[P.HotelPlaceType]; + l += Mcells.at(Airports[i].DestMcells[j].id).np[P.HotelPlaceType]; if (l < 10) { fprintf(stderr, "(%i ", l); l = 0; for (j = 0; j < Airports[i].num_mcell; j++) - l += Mcells[Airports[i].DestMcells[j].id].n; + l += Mcells.at(Airports[i].DestMcells[j].id).n; fprintf(stderr, "%i %i) ", Airports[i].num_mcell, l); } } @@ -1635,7 +1635,7 @@ void AssignHouseholdAges(int n, int pers, int tn) int i, j, k, l, nc, ad; int a[MAX_HOUSEHOLD_SIZE + 2]; - ad = ((P.DoAdunitDemog) && (P.DoAdUnits)) ? Mcells[Hosts[pers].mcell].adunit : 0; + ad = ((P.DoAdunitDemog) && (P.DoAdUnits)) ? Mcells.at(Hosts[pers].mcell).adunit : 0; if (!P.DoHouseholds) { for (i = 0; i < n; i++) @@ -1833,7 +1833,7 @@ void AssignPeopleToPlaces(void) cnt = 0; for (a = 0; a < P.NCP; a++) { - cell *c = CellLookup[a]; + cell *c = CellLookup.at(a); c->n = 0; for (j = 0; j < c->cumTC; j++) { @@ -1857,7 +1857,7 @@ void AssignPeopleToPlaces(void) j2 = 0; for (a = 0; a < P.NCP; a++) { - cell *c = CellLookup[a]; + cell *c = CellLookup.at(a); for (j = 0; j < c->n; j++) { PeopleArray[j2] = c->susceptible[j]; @@ -2074,9 +2074,9 @@ void AssignPeopleToPlaces(void) if ((mx >= 0) && (my >= 0) && (mx < P.nmcw) && (my < P.nmch)) { ic = mx * P.nmch + my; - if (Mcells[ic].country == Mcells[Hosts[i].mcell].country) + if (Mcells[ic].country == Mcells.at(Hosts[i].mcell).country) { - for (cnt = 0; cnt < Mcells[ic].np[tp]; cnt++) + for (cnt = 0; cnt < Mcells.at(ic).np[tp]; cnt++) { if (Mcells[ic].places[tp][cnt] >= P.Nplace[tp]) fprintf(stderr, "#%i %i %i ", tp, ic, cnt); t = dist2_raw(Households[Hosts[i].hh].loc_x, Households[Hosts[i].hh].loc_y, @@ -2100,7 +2100,7 @@ void AssignPeopleToPlaces(void) { if (k < nn) { - NearestPlaces[tn][k] = Mcells[ic].places[tp][cnt]; + NearestPlaces[tn][k] = Mcells.at(ic).places[tp][cnt]; NearestPlacesProb[tn][k] = s; k++; } @@ -2116,7 +2116,7 @@ void AssignPeopleToPlaces(void) if (s > t) { NearestPlacesProb[tn][j2] = s; - NearestPlaces[tn][j2] = Mcells[ic].places[tp][cnt]; + NearestPlaces[tn][j2] = Mcells.at(ic).places[tp][cnt]; } } } @@ -2213,7 +2213,7 @@ void AssignPeopleToPlaces(void) s = ranf(); l = Cells.at(i).InvCDF[(int)floor(s * 1024)]; while (Cells.at(i).cum_trans[l] < s) l++; - ct = CellLookup[l]; + ct = CellLookup.at(l); m = (int)(ranf() * ((double)ct->S)); j = -1; #pragma omp critical @@ -2235,7 +2235,7 @@ void AssignPeopleToPlaces(void) s = ((double)ct->S) / ((double)ct->S0) * numKernel(t) / Cells.at(i).max_trans[l]; if ((P.DoAdUnits) && (P.InhibitInterAdunitPlaceAssignment[tp] > 0)) { - if (Mcells[Hosts[k].mcell].adunit != Mcells[Places[tp][j].mcell].adunit) s *= (1 - P.InhibitInterAdunitPlaceAssignment[tp]); + if (Mcells[Hosts[k].mcell].adunit != Mcells.at(Places[tp][j].mcell).adunit) s *= (1 - P.InhibitInterAdunitPlaceAssignment[tp]); } if (ranf() < s) {