Merge pull request #388 from ozmorph/assign_people_to_places_optimiza…

…tion

US Regression Test Optimization

covid-sim.vcxproj

@@ -104,7 +104,6 @@
     <ClCompile Include="src\Sweep.cpp" />
     <ClCompile Include="src\Update.cpp" />
     <ClCompile Include="src\Param.cpp" />
-    <ClCompile Include="src\MicroCellPosition.cpp" />
     <ClCompile Include="src\Direction.cpp" />
     <ClCompile Include="src\Geometry\Vector2.cpp" />
     <ClCompile Include="src\InverseCdf.cpp" />

src/Bitmap.cpp

@@ -61,16 +61,16 @@ void CaptureBitmap()
 			if (Mcells[i].n > 0)
 			{
 				f = 0;
-				if ((i < P.NMC - 1) && (i / P.get_number_of_micro_cells_high() == (i + 1) / P.get_number_of_micro_cells_high()) && (Mcells[i + 1].n > 0) && ((Mcells[i].country != Mcells[i + 1].country)
+				if ((i < P.NMC - 1) && (i / P.total_microcells_high_ == (i + 1) / P.total_microcells_high_) && (Mcells[i + 1].n > 0) && ((mcell_country[i] != mcell_country[i + 1])
 					|| ((P.DoAdunitBoundaryOutput) && ((AdUnits[Mcells[i].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor != (AdUnits[Mcells[i + 1].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor)))) f = 1;
-				if ((i > 0) && (i / P.get_number_of_micro_cells_high() == (i - 1) / P.get_number_of_micro_cells_high()) && (Mcells[i - 1].n > 0) && (Mcells[i].country != Mcells[i - 1].country)) f = 1;
-				if ((i < P.NMC - P.get_number_of_micro_cells_high()) && (Mcells[i + P.get_number_of_micro_cells_high()].n > 0) && ((Mcells[i].country != Mcells[i + P.get_number_of_micro_cells_high()].country)
-					|| ((P.DoAdunitBoundaryOutput) && ((AdUnits[Mcells[i].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor != (AdUnits[Mcells[i + P.get_number_of_micro_cells_high()].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor)))) f = 1;
-				if ((i >= P.get_number_of_micro_cells_high()) && (Mcells[i - P.get_number_of_micro_cells_high()].n > 0) && (Mcells[i].country != Mcells[i - P.get_number_of_micro_cells_high()].country)) f = 1;
+				if ((i > 0) && (i / P.total_microcells_high_ == (i - 1) / P.total_microcells_high_) && (Mcells[i - 1].n > 0) && (mcell_country[i] != mcell_country[i - 1])) f = 1;
+				if ((i < P.NMC - P.total_microcells_high_) && (Mcells[i + P.total_microcells_high_].n > 0) && ((mcell_country[i] != mcell_country[i + P.total_microcells_high_])
+					|| ((P.DoAdunitBoundaryOutput) && ((AdUnits[Mcells[i].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor != (AdUnits[Mcells[i + P.total_microcells_high_].adunit].id % P.AdunitLevel1Mask) / P.AdunitBitmapDivisor)))) f = 1;
+				if ((i >= P.total_microcells_high_) && (Mcells[i - P.total_microcells_high_].n > 0) && (mcell_country[i] != mcell_country[i - P.total_microcells_high_])) f = 1;
 				if (f)
 				{
-					x = (int)(P.in_microcells_.width * (((double)(i / P.get_number_of_micro_cells_high())) + 0.5) * P.scale.x) - P.bmin.x;
-					y = (int)(P.in_microcells_.height * (((double)(i % P.get_number_of_micro_cells_high())) + 0.5) * P.scale.y) - P.bmin.y;
+					x = (int)(P.in_microcells_.width * (((double)(i / P.total_microcells_high_)) + 0.5) * P.scale.x) - P.bmin.x;
+					y = (int)(P.in_microcells_.height * (((double)(i % P.total_microcells_high_)) + 0.5) * P.scale.y) - P.bmin.y;
 					if ((x >= 0) && (x < P.b.width) && (y >= 0) && (y < P.b.height))
 					{
 						j = y * bmh->width + x;

src/CMakeLists.txt

@@ -3,8 +3,7 @@
 # Set up the IDE
 set(MAIN_SRC_FILES CovidSim.cpp BinIO.cpp Rand.cpp Error.cpp Dist.cpp
   Kernels.cpp Bitmap.cpp SetupModel.cpp CalcInfSusc.cpp Sweep.cpp Update.cpp
-  Param.cpp MicroCellPosition.cpp Direction.cpp InverseCdf.cpp Memory.cpp
-  CLI.cpp)
+  Param.cpp Direction.cpp InverseCdf.cpp Memory.cpp CLI.cpp)
 set(MAIN_HDR_FILES CovidSim.h BinIO.h Rand.h Constants.h Country.h Error.h
   Dist.h Kernels.h Bitmap.h Model.h Param.h SetupModel.h ModelMacros.h
   InfStat.h CalcInfSusc.h Sweep.h Update.h MicroCellPosition.hpp Direction.hpp

src/CovidSim.cpp

@@ -4,9 +4,10 @@
 
 #include <algorithm>
 #include <cctype>
-#include <cstdlib>
 #include <cerrno>
 #include <cstddef>
+#include <cstdint>
+#include <cstdlib>
 #include <fstream>
 #include <iostream>
 #include <string>
@@ -77,11 +78,13 @@ int GetInputParameter3(FILE*, const char*, const char*, void*, int, int, int);
 
 Param P;
 Person* Hosts;
+std::vector<PersonQuarantine> HostsQuarantine;
 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;
+std::vector<uint16_t> mcell_country;
 Place** Places;
 AdminUnit AdUnits[MAX_ADUNITS];
 //// Time Series defs:
@@ -2470,7 +2473,7 @@ void InitModel(int run) // passing run number so we can save run number in the i
 					StateT[Thread].cumInf_age_adunit [AgeGroup][Adunit] = 0;
 				}
 
-
+	std::fill(HostsQuarantine.begin(), HostsQuarantine.end(), PersonQuarantine());
 #pragma omp parallel for schedule(static,1) default(none) \
 		shared(P, Hosts)
 	for (int tn = 0; tn < P.NumThreads; tn++)
@@ -2479,9 +2482,8 @@ void InitModel(int run) // passing run number so we can save run number in the i
 			Hosts[k].absent_start_time = USHRT_MAX - 1;
 			Hosts[k].absent_stop_time = 0;
 			if (P.DoAirports) Hosts[k].PlaceLinks[P.HotelPlaceType] = -1;
-			Hosts[k].vacc_start_time = Hosts[k].treat_start_time = Hosts[k].quar_start_time = Hosts[k].isolation_start_time = Hosts[k].absent_start_time = Hosts[k].dct_start_time = Hosts[k].dct_trigger_time = USHRT_MAX - 1;
+			Hosts[k].vacc_start_time = Hosts[k].treat_start_time = Hosts[k].isolation_start_time = Hosts[k].absent_start_time = Hosts[k].dct_start_time = Hosts[k].dct_trigger_time = USHRT_MAX - 1;
 			Hosts[k].treat_stop_time = Hosts[k].absent_stop_time = Hosts[k].dct_end_time = 0;
-			Hosts[k].quar_comply = 2;
 			Hosts[k].to_die = 0;
 			Hosts[k].Travelling = 0;
 			Hosts[k].detected = 0; //set detected to zero initially: ggilani - 19/02/15
@@ -2717,7 +2719,7 @@ void SeedInfection(double t, int* NumSeedingInfections_byLocation, int rf, int r
 		{
 			k = (int)(P.LocationInitialInfection[i][0] / P.in_microcells_.width);
 			l = (int)(P.LocationInitialInfection[i][1] / P.in_microcells_.height);
-			j = k * P.get_number_of_micro_cells_high() + l;
+			j = k * P.total_microcells_high_ + l;
 			m = 0;
 			for (k = 0; (k < NumSeedingInfections_byLocation[i]) && (m < 10000); k++)
 			{
@@ -5178,7 +5180,7 @@ void RecordInfTypes(void)
 			{
 				if (Hosts[i].latent_time * P.TimeStep <= P.SampleTime)
 					TimeSeries[(int)(Hosts[i].latent_time * P.TimeStep / P.SampleStep)].Rdenom++;
-				infcountry[Mcells[Hosts[i].mcell].country]++;
+				infcountry[mcell_country[Hosts[i].mcell]]++;
 				if (abs(Hosts[i].inf) < InfStat_Recovered)
 					l = -1;
 				else if (l >= 0)
@@ -5330,15 +5332,15 @@ void CalcOriginDestMatrix_adunit()
 
 			//find index of cell from which flow travels
 			ptrdiff_t cl_from = CellLookup[i] - Cells;
-			ptrdiff_t cl_from_mcl = (cl_from / P.nch) * P.NMCL * P.get_number_of_micro_cells_high() + (cl_from % P.nch) * P.NMCL;
+			ptrdiff_t cl_from_mcl = (cl_from / P.nch) * P.NMCL * P.total_microcells_high_ + (cl_from % P.nch) * P.NMCL;
 
 			//loop over microcells in these cells to find populations in each admin unit and so flows
 			for (int k = 0; k < P.NMCL; k++)
 			{
 				for (int l = 0; l < P.NMCL; l++)
 				{
 					//get index of microcell
-					ptrdiff_t mcl_from = cl_from_mcl + l + k * P.get_number_of_micro_cells_high();
+					ptrdiff_t mcl_from = cl_from_mcl + l + k * P.total_microcells_high_;
 					if (Mcells[mcl_from].n > 0)
 					{
 						//get proportion of each population of cell that exists in each admin unit
@@ -5354,7 +5356,7 @@ void CalcOriginDestMatrix_adunit()
 
 				//find index of cell which flow travels to
 				ptrdiff_t cl_to = CellLookup[j] - Cells;
-				ptrdiff_t cl_to_mcl = (cl_to / P.nch) * P.NMCL * P.get_number_of_micro_cells_high() + (cl_to % P.nch) * P.NMCL;
+				ptrdiff_t cl_to_mcl = (cl_to / P.nch) * P.NMCL * P.total_microcells_high_ + (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;
 				double total_flow;
@@ -5373,7 +5375,7 @@ void CalcOriginDestMatrix_adunit()
 					for (int p = 0; p < P.NMCL; p++)
 					{
 						//get index of microcell
-						ptrdiff_t mcl_to = cl_to_mcl + p + m * P.get_number_of_micro_cells_high();
+						ptrdiff_t mcl_to = cl_to_mcl + p + m * P.total_microcells_high_;
 						if (Mcells[mcl_to].n > 0)
 						{
 							//get proportion of each population of cell that exists in each admin unit

src/Dist.cpp

@@ -152,12 +152,12 @@ double dist2_mm(Microcell* a, Microcell* b)
 	l = (int)(a - Mcells);
 	m = (int)(b - Mcells);
 	if (P.DoUTM_coords)
-		return dist2UTM(P.in_microcells_.width * fabs((double)(l / P.get_number_of_micro_cells_high())), P.in_microcells_.height * fabs((double)(l % P.get_number_of_micro_cells_high())),
-			P.in_microcells_.width * fabs((double)(m / P.get_number_of_micro_cells_high())), P.in_microcells_.height * fabs((double)(m % P.get_number_of_micro_cells_high())));
+		return dist2UTM(P.in_microcells_.width * fabs((double)(l / P.total_microcells_high_)), P.in_microcells_.height * fabs((double)(l % P.total_microcells_high_)),
+			P.in_microcells_.width * fabs((double)(m / P.total_microcells_high_)), P.in_microcells_.height * fabs((double)(m % P.total_microcells_high_)));
 	else
 	{
-		x = P.in_microcells_.width * fabs((double)(l / P.get_number_of_micro_cells_high() - m / P.get_number_of_micro_cells_high()));
-		y = P.in_microcells_.height * fabs((double)(l % P.get_number_of_micro_cells_high() - m % P.get_number_of_micro_cells_high()));
+		x = P.in_microcells_.width * fabs((double)(l / P.total_microcells_high_ - m / P.total_microcells_high_));
+		y = P.in_microcells_.height * fabs((double)(l % P.total_microcells_high_ - m % P.total_microcells_high_));
 		return periodic_xy(x, y);
 	}
 }

src/MicroCellPosition.hpp

@@ -1,13 +1,22 @@
 #pragma once
 
 #include "Direction.hpp"
+#include "Error.h"
 
-class MicroCellPosition {
-public:
+struct MicroCellPosition
+{
 	int x;
 	int y;
 
-	MicroCellPosition operator+(Direction direction) const;
-	MicroCellPosition& operator+=(Direction direction);
+	inline MicroCellPosition& operator+=(Direction direction)
+	{
+		switch (direction) {
+			case Right: this->x += 1; break;
+			case Up: this->y -= 1; break;
+			case Left: this->x -= 1; break;
+			case Down: this->y += 1; break;
+			default: ERR_CRITICAL("Unknown direction");
+		}
+		return *this;
+	}
 };
-

src/Model.h

@@ -1,7 +1,9 @@
 #ifndef COVIDSIM_MODEL_H_INCLUDED_
 #define COVIDSIM_MODEL_H_INCLUDED_
 
+#include <cstdint>
 #include <cstddef>
+#include <vector>
 
 #include "Country.h"
 #include "Constants.h"
@@ -291,10 +293,12 @@ struct AdminUnit
 #pragma pack(pop)
 
 extern Person* Hosts;
+extern std::vector<PersonQuarantine> HostsQuarantine;
 extern Household* Households;
 extern PopVar State, StateT[MAX_NUM_THREADS];
 extern Cell* Cells, ** CellLookup;
 extern Microcell* Mcells, ** McellLookup;
+extern std::vector<uint16_t> mcell_country;
 extern Place** Places;
 extern AdminUnit AdUnits[MAX_ADUNITS];
 

src/ModelMacros.h

@@ -10,8 +10,8 @@
 #define HOST_TO_BE_VACCED(x) (Hosts[x].vacc_start_time < USHRT_MAX - 1)
 #define HOST_VACCED(x) (Hosts[x].vacc_start_time+P.usVaccTimeToEfficacy<=ts)
 #define HOST_VACCED_SWITCH(x) (Hosts[x].vacc_start_time >= P.usVaccTimeEfficacySwitch)
-#define HOST_QUARANTINED(x) ((Hosts[x].quar_comply == 1) && (Hosts[x].quar_start_time + P.usHQuarantineHouseDuration > ts) && (Hosts[x].quar_start_time <= ts))
-#define HOST_TO_BE_QUARANTINED(x) ((Hosts[x].quar_start_time + P.usHQuarantineHouseDuration > ts) && (Hosts[x].quar_comply < 2))
+#define HOST_QUARANTINED(x) ((HostsQuarantine[x].comply == 1) && (HostsQuarantine[x].start_time + P.usHQuarantineHouseDuration > ts) && (HostsQuarantine[x].start_time <= ts))
+#define HOST_TO_BE_QUARANTINED(x) ((HostsQuarantine[x].start_time + P.usHQuarantineHouseDuration > ts) && (HostsQuarantine[x].comply < 2))
 #define HOST_ISOLATED(x) ((Hosts[x].isolation_start_time + P.usCaseIsolationDelay <= ts) && (Hosts[x].isolation_start_time + P.usCaseIsolationDelay + P.usCaseIsolationDuration > ts))
 #define HOST_ABSENT(x) ((Hosts[x].absent_start_time <= ts) && (Hosts[x].absent_stop_time > ts))
 

src/Models/Microcell.h

@@ -16,7 +16,6 @@ struct Microcell
 	*/
 	int n /*Number of people in microcell*/, adunit;
 	int* members;
-	unsigned short int country;
 
 	int* places[NUM_PLACE_TYPES];
 	unsigned short int np[NUM_PLACE_TYPES];

src/Models/Person.h

@@ -1,5 +1,8 @@
 #pragma once
 
+#include <cstdint>
+#include <limits>
+
 #include "../Country.h"
 #include "../InfStat.h"
 
@@ -21,7 +24,6 @@ struct Person
 
 	unsigned char Travelling;	// Range up to MAX_TRAVEL_TIME
 	unsigned char age;
-	unsigned char quar_comply;		// can be 0, 1, or 2
 	unsigned char num_treats;		// set to 0 and tested < 2. but never modified?
 	Severity Severity_Current, Severity_Final; //// Note we allow Severity_Final to take values: Severity_Mild, Severity_ILI, Severity_SARI, Severity_Critical (not e.g. Severity_Dead or Severity_RecoveringFromCritical)
 
@@ -31,7 +33,7 @@ struct Person
 
 	unsigned short int detected_time; //added hospitalisation flag: ggilani 28/10/2014, added flag to determined whether this person's infection is detected or not
 	unsigned short int absent_start_time, absent_stop_time;
-	unsigned short int quar_start_time, isolation_start_time;
+	unsigned short int isolation_start_time;
 	unsigned short int infection_time, latent_time;		// Set in DoInfect function. infection time is time of infection; latent_time is a misnomer - it is the time at which person become infectious (i.e. infection time + latent period for this person). latent_time will also refer to time of onset with ILI or Mild symptomatic disease.
 	unsigned short int recovery_or_death_time;	// set in DoIncub function
 	unsigned short int SARI_time, Critical_time, RecoveringFromCritical_time; //// /*mild_time, ILI_time,*/ Time of infectiousness onset same for asymptomatic, Mild, and ILI infection so don't need mild_time etc.
@@ -43,3 +45,13 @@ struct Person
 	int ncontacts; //added this in to record total number of contacts each index case records: ggilani 13/04/20
 
 };
+
+struct PersonQuarantine
+{
+	uint8_t  comply;		// can be 0, 1, 2
+	uint16_t start_time;	// timestep quarantine is started
+
+	// don't remove the extra parentheses around std::numeric_limits<uint16_t>::max
+	// because it conflicts with the max() preprocessor macro in MSVC builds
+	PersonQuarantine() : comply(2), start_time((std::numeric_limits<uint16_t>::max)()-1) {} 
+};

src/Param.cpp

@@ -1,28 +1,20 @@
 #include "Param.h"
 
-int Param::get_number_of_micro_cells_wide() const {
-	return this->ncw * this->NMCL;
-}
-
-int Param::get_number_of_micro_cells_high() const {
-	return this->nch * this->NMCL;
-}
-
 MicroCellPosition Param::get_micro_cell_position_from_cell_index(int cell_index) const {
-	int x = cell_index / this->get_number_of_micro_cells_high();
-	int y = cell_index % this->get_number_of_micro_cells_high();
+	int x = cell_index / this->total_microcells_high_;
+	int y = cell_index % this->total_microcells_high_;
 	return {x, y};
 }
 
-bool Param::is_in_bounds(MicroCellPosition position) const {
+bool Param::is_in_bounds(MicroCellPosition const& position) const {
 	return position.x >= 0
 		&& position.y >= 0
-		&& position.x < this->get_number_of_micro_cells_wide()
-		&& position.y < this->get_number_of_micro_cells_high();
+		&& position.x < this->total_microcells_wide_
+		&& position.y < this->total_microcells_high_;
 }
 
-int Param::get_micro_cell_index_from_position(MicroCellPosition position) const {
-	int x = (position.x + this->get_number_of_micro_cells_wide()) % this->get_number_of_micro_cells_wide();
-	int y = (position.y + this->get_number_of_micro_cells_high()) % this->get_number_of_micro_cells_high();
-	return x * this->get_number_of_micro_cells_high() + y;
+int Param::get_micro_cell_index_from_position(MicroCellPosition const& position) const {
+	int x = (position.x + this->total_microcells_wide_) % this->total_microcells_wide_;
+	int y = (position.y + this->total_microcells_high_) % this->total_microcells_high_;
+	return x * this->total_microcells_high_ + y;
 }

src/Param.h

@@ -45,12 +45,11 @@ struct Param
 	int NMCL; // Number of microcells wide/high a cell is; i.e. NMC = NC * NMCL * NMCL
 	int NCP; /**< Number of populated cells  */
 	int NMCP, ncw, nch, DoUTM_coords, nsp, DoSeasonality, DoCorrectAgeDist, DoPartialImmunity;
+	int total_microcells_wide_, total_microcells_high_;
 
-	int get_number_of_micro_cells_wide() const;
-	int get_number_of_micro_cells_high() const;
 	MicroCellPosition get_micro_cell_position_from_cell_index(int cell_index) const;
-	int get_micro_cell_index_from_position(MicroCellPosition position) const;
-	bool is_in_bounds(MicroCellPosition position) const;
+	int get_micro_cell_index_from_position(MicroCellPosition const& position) const;
+	bool is_in_bounds(MicroCellPosition const& position) const;
 
 	int DoAdUnits, NumAdunits, DoAdunitBoundaries, AdunitLevel1Divisor, AdunitLevel1Mask, AdunitBitmapDivisor, CountryDivisor;
 	int DoAdunitOutput, DoAdunitBoundaryOutput, DoCorrectAdunitPop, DoSpecifyPop, AdunitLevel1Lookup[ADUNIT_LOOKUP_SIZE];