Dlaydon/vary betas

ci/install_macos_dependencies.sh

@@ -1,3 +1,6 @@
 #!/bin/sh -e
 
 brew install llvm cmake libomp python3 coreutils
+export CMAKE_CXX_FLAGS="$CMAKE_CXX_FLAGS -Xpreprocessor -fopenmp"
+export OpenMP_C_LIB_NAMES="libomp"
+export OpenMP_C_FLAGS="-fopenmp"

covid-sim.vcxproj

@@ -21,13 +21,13 @@
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Debug|x64'" Label="Configuration">
     <ConfigurationType>Application</ConfigurationType>
     <UseDebugLibraries>true</UseDebugLibraries>
-    <PlatformToolset>v142</PlatformToolset>
+    <PlatformToolset>v143</PlatformToolset>
     <CharacterSet>MultiByte</CharacterSet>
   </PropertyGroup>
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|x64'" Label="Configuration">
     <ConfigurationType>Application</ConfigurationType>
     <UseDebugLibraries>false</UseDebugLibraries>
-    <PlatformToolset>v142</PlatformToolset>
+    <PlatformToolset>v143</PlatformToolset>
     <WholeProgramOptimization>true</WholeProgramOptimization>
     <CharacterSet>MultiByte</CharacterSet>
   </PropertyGroup>

src/CalcInfSusc.cpp

@@ -28,7 +28,7 @@ double CalcPlaceInf(int person, int PlaceType, unsigned short int TimeStepNow)
 		*	((Hosts[person].digitalContactTraced == 1) ? P.Efficacies[DigContactTracing][PlaceType] : 1.0)
 		*	((HOST_QUARANTINED(person) && (!Hosts[person].care_home_resident) && (Hosts[person].digitalContactTraced != 1) && (!(HOST_ISOLATED(person)))) ? P.Efficacies[HomeQuarantine][PlaceType] : 1.0)
 		*	((Hosts[person].is_case() && (!Hosts[person].care_home_resident)) ? P.SymptPlaceTypeContactRate[PlaceType] : 1.0)
-		*	P.PlaceTypeTrans[PlaceType] / P.PlaceTypeGroupSizeParam1[PlaceType] * CalcPersonInf(person, TimeStepNow);
+		*	CalcPersonInf(person, TimeStepNow);
 }
 
 double CalcSpatialInf(int person, unsigned short int TimeStepNow)

src/Constants.h

@@ -110,10 +110,10 @@ const int MAX_NUM_CFR_CHANGE_TIMES = 100; /**< To allow IFR to scale over time.
 #define _I64(x) static_cast<int64_t>(x)
 
 // Settings (numbers not arbitrary - don't change without checking)
-// const int PrimarySchool		= 0; 
-// const int SecondarySchool	= 1; 
-// const int University		= 2; 
-// const int Workplace			= 3; 
+const int PrimarySchool		= 0; 
+const int SecondarySchool	= 1; 
+const int University		= 2; 
+const int Workplace			= 3; 
 const int House				= MAX_NUM_PLACE_TYPES;      // Max number of potential place types.
 const int Spatial			= MAX_NUM_PLACE_TYPES + 1;  // community
 

src/Models/Cell.h

@@ -12,7 +12,9 @@ struct Cell
 {
 	int n; /**< number of people in cell */
 	int S, L, I, R, D; /**< S, L, I, R, D are numbers of Susceptible, Latently infected, Infectious, Recovered and Dead people in cell */
-	int cumTC, S0, tot_treat, tot_vacc;
+	int cumTC, S0;
+	int tot_treat; //// total number of people in this cell who have been treated. Initialized in CovidSim.cpp::InitModel and incremented in DoTreat, DoProph and DoProphNoDelay
+	int tot_vacc; //// total number of people in this cell who have been vaccinated. Initialized in CovidSim.cpp::InitModel and incremented in DoVacc or DoVaccNoDelay
 	int* members, *susceptible, *latent, *infected; /**< pointers to people in cell. e.g. *susceptible identifies where the final susceptible member of cell is. */
 	int* InvCDF;
 	float tot_prob, *cum_trans, *max_trans;

src/Param.h

@@ -36,6 +36,7 @@ struct Param
 	/**< Time-step defintions. Differentiates between length of time between model updates (ModelTimeStep) and length of time between calculating model outputs (OutputTimeStep). */
 	double SimulationDuration;				/**< The number of days to run for */
 	double ModelTimeStep;					/**< The length of a time step, in days */
+	double TimeStepsPerDay; 
 	double OutputTimeStep;					/**< The length of time in days between calculating model outputs, in days. Note ModelTimeStep <= OutputTimeStep. */
 	int NumModelTimeStepsPerOutputTimeStep;	/**< Number of time steps between samples. NumModelTimeStepsPerOutputTimeStep = OutputTimeStep / ModelTimeStep */
 	int NumOutputTimeSteps;					/**< Total number of time output steps that will be made. NumOutputTimeSteps = SimulationDuration / OutputTimeStep */ 
@@ -108,13 +109,23 @@ struct Param
 
 	CovidSim::Geometry::BoundingBox2d SpatialBoundingBox;
 	double** LocationInitialInfection;
-	double InitialInfectionsAdminUnitWeight[MAX_NUM_SEED_LOCATIONS], InitialInfectionCalTime, TimeStepsPerDay;
+	double InitialInfectionsAdminUnitWeight[MAX_NUM_SEED_LOCATIONS], InitialInfectionCalTime;
 	double FalsePositiveRate, FalsePositivePerCapitaIncidence, FalsePositiveAgeRate[NUM_AGE_GROUPS];
 	double SeroConvMaxSens, SeroConvP1, SeroConvP2, SeroConvSpec, InfPrevSurveyScale;
 
 	double AirportTrafficScale;
 
-	double R0, R0scale, LocalBeta;
+	/**< Betas array (average number of infectious contacts per day at time t in given region for particular setting)
+	// by default, betas do not vary over time or between regions, unless P.VaryBetasOverTimeByRegion == true
+	// setting refers to either place, household and spatial/community
+	// Flow is to read in mobility data (proportional change in mobility type by region over time).
+	// If VaryBetasOverTimeByRegion == false, rather than read in mobility data, first set all betas to 1, then multipy each day and regions value by P.LocalBeta, P.HouseholdTrans etc. 
+	*/
+	double*** Betas; //// indexed by i) time; ii) region / admin unit; iii) setting.
+	bool VaryBetasOverTimeByRegion = false;
+
+	double R0, R0scale;
+	double LocalBeta; //// Spatial/community beta: (roughly the average number of infectious spatial/community contacts per person per unit time)
 	double infectious_prof[INFPROF_RES + 1], infectiousness[MAX_INFECTIOUS_STEPS];
 	double InfectiousPeriod; // In days. Mean of icdf (inverse cumulative distribution function).
 	double R0household, R0places, R0spatial;
@@ -140,7 +151,8 @@ struct Param
 	int AbsenteeismPlaceClosure; // Default 0 (off), 1 (on) track place closures in more detail
 	int MaxAbsentTime; // In days.  Max number of days absent, range [0, MAX_ABSENT_TIME].  Default 0 if !P.AbsenteeismPlaceClosure, otherwise MAX_ABSENT_TIME
 	int InvJourneyDurationDistrib[1025], InvLocalJourneyDurationDistrib[1025];
-	double HouseholdTrans, HouseholdTransPow;
+	double HouseholdTrans; //// household beta (roughly the average number of infectious household contacts per person per unit time)
+	double HouseholdTransPow;
 	double** HouseholdSizeDistrib; // [MAX_ADUNITS] [MAX_HOUSEHOLD_SIZE]
 	double HouseholdDenomLookup[MAX_HOUSEHOLD_SIZE];
 	int PlaceTypeAgeMin[MAX_NUM_PLACE_TYPES], PlaceTypeAgeMax[MAX_NUM_PLACE_TYPES], PlaceTypeMaxAgeRead[MAX_NUM_PLACE_TYPES];
@@ -149,7 +161,8 @@ struct Param
 	int PlaceTypeNearestNeighb[MAX_NUM_PLACE_TYPES], PlaceTypeKernelType[MAX_NUM_PLACE_TYPES];
 	double PlaceTypePropAgeGroup[MAX_NUM_PLACE_TYPES], PlaceTypePropAgeGroup2[MAX_NUM_PLACE_TYPES];
 	double PlaceTypePropAgeGroup3[MAX_NUM_PLACE_TYPES], PlaceTypeKernelShape[MAX_NUM_PLACE_TYPES], PlaceTypeKernelScale[MAX_NUM_PLACE_TYPES];
-	double PlaceTypeKernelP3[MAX_NUM_PLACE_TYPES], PlaceTypeKernelP4[MAX_NUM_PLACE_TYPES], PlaceTypeTrans[MAX_NUM_PLACE_TYPES];
+	double PlaceTypeKernelP3[MAX_NUM_PLACE_TYPES], PlaceTypeKernelP4[MAX_NUM_PLACE_TYPES];
+	double PlaceTypeTrans[MAX_NUM_PLACE_TYPES]; //// indexed by place time place type betas (roughly average number of infectious place contacts per person per unit time).
 	double PlaceTypeMeanSize[MAX_NUM_PLACE_TYPES], PlaceTypePropBetweenGroupLinks[MAX_NUM_PLACE_TYPES], PlaceTypeSizeSD[MAX_NUM_PLACE_TYPES]; //added PlaceTypeSizeSD for lognormal distribution - ggilani 09/02/17
 	double PlaceTypeSizePower[MAX_NUM_PLACE_TYPES], PlaceTypeSizeOffset[MAX_NUM_PLACE_TYPES], PlaceTypeSizeMax[MAX_NUM_PLACE_TYPES], PlaceTypeSizeMin[MAX_NUM_PLACE_TYPES];
 	double PlaceTypeGroupSizeParam1[MAX_NUM_PLACE_TYPES], PlaceExclusivityMatrix[MAX_NUM_PLACE_TYPES * MAX_NUM_PLACE_TYPES]; //changed PlaceExclusivityMatrix from [MAX_NUM_PLACE_TYPES][MAX_NUM_PLACE_TYPES]
@@ -263,6 +276,7 @@ struct Param
 
 	int NumInterventionClasses = 0, NumInfectionSettings = 0; // initialized in CovidSim.cpp::main.
 
+
 	///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****
 	///// **** VARIABLE EFFICACIES OVER TIME
 	///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****

src/ReadParams.cpp

@@ -2272,30 +2272,30 @@ void Params::ReadParams(std::string const& ParamFile, std::string const& PrePara
 		//Add origin-destination matrix parameter
 		P->DoOriginDestinationMatrix = Params::get_int(params, pre_params, "Output origin destination matrix", 0, P);
 
-		P->MeanChildAgeGap = Params::get_int(params, pre_params, adm_params, "Mean child age gap", 2, P);
-		P->MinAdultAge = Params::get_int(params, pre_params, adm_params, "Min adult age", 19, P);
-		P->MaxMFPartnerAgeGap = Params::get_int(params, pre_params, adm_params, "Max MF partner age gap", 5, P);
-		P->MaxFMPartnerAgeGap = Params::get_int(params, pre_params, adm_params, "Max FM partner age gap", 5, P);
-		P->MinParentAgeGap = Params::get_int(params, pre_params, adm_params, "Min parent age gap", 19, P);
-		P->MaxParentAgeGap = Params::get_int(params, pre_params, adm_params, "Max parent age gap", 44, P);
-		P->MaxChildAge = Params::get_int(params, pre_params, adm_params, "Max child age", 20, P);
-		P->OneChildTwoPersProb = Params::get_double(params, pre_params, adm_params, "One Child Two Pers Prob", 0.08, P);
-		P->TwoChildThreePersProb = Params::get_double(params, pre_params, adm_params, "Two Child Three Pers Prob", 0.11, P);
-		P->OnePersHouseProbOld = Params::get_double(params, pre_params, adm_params, "One Pers House Prob Old", 0.5, P);
-		P->TwoPersHouseProbOld = Params::get_double(params, pre_params, adm_params, "Two Pers House Prob Old", 0.5, P);
-		P->OnePersHouseProbYoung = Params::get_double(params, pre_params, adm_params, "One Pers House Prob Young", 0.23, P);
-		P->TwoPersHouseProbYoung = Params::get_double(params, pre_params, adm_params, "Two Pers House Prob Young", 0.23, P);
-		P->OneChildProbYoungestChildUnderFive = Params::get_double(params, pre_params, adm_params, "One Child Prob Youngest Child Under Five", 0.5, P);
-		P->TwoChildrenProbYoungestUnderFive = Params::get_double(params, pre_params, adm_params, "Two Children Prob Youngest Under Five", 0.0, P);
-		P->TwoChildrenProbYoungestUnderFive = Params::get_double(params, pre_params, adm_params, "Prob Youngest Child Under Five", 0, P);
-		P->ZeroChildThreePersProb = Params::get_double(params, pre_params, adm_params, "Zero Child Three Pers Prob", 0.25, P);
-		P->OneChildFourPersProb = Params::get_double(params, pre_params, adm_params, "One Child Four Pers Prob", 0.2, P);
-		P->YoungAndSingleSlope = Params::get_double(params, pre_params, adm_params, "Young And Single Slope", 0.7, P);
-		P->YoungAndSingle = Params::get_int(params, pre_params, adm_params, "Young And Single", 36, P);
-		P->NoChildPersAge = Params::get_int(params, pre_params, adm_params, "No Child Pers Age", 44, P);
-		P->OldPersAge = Params::get_int(params, pre_params, adm_params, "Old Pers Age", 60, P);
-		P->ThreeChildFivePersProb = Params::get_double(params, pre_params, adm_params, "Three Child Five Pers Prob", 0.5, P);
-		P->OlderGenGap = Params::get_int(params, pre_params, adm_params, "Older Gen Gap", 19, P);
+		P->MeanChildAgeGap						= Params::get_int(params, pre_params, adm_params, "Mean child age gap", 2, P);
+		P->MinAdultAge							= Params::get_int(params, pre_params, adm_params, "Min adult age", 19, P);
+		P->MaxMFPartnerAgeGap					= Params::get_int(params, pre_params, adm_params, "Max MF partner age gap", 5, P);
+		P->MaxFMPartnerAgeGap					= Params::get_int(params, pre_params, adm_params, "Max FM partner age gap", 5, P);
+		P->MinParentAgeGap						= Params::get_int(params, pre_params, adm_params, "Min parent age gap", 19, P);
+		P->MaxParentAgeGap						= Params::get_int(params, pre_params, adm_params, "Max parent age gap", 44, P);
+		P->MaxChildAge							= Params::get_int(params, pre_params, adm_params, "Max child age", 20, P);
+		P->OneChildTwoPersProb					= Params::get_double(params, pre_params, adm_params, "One Child Two Pers Prob", 0.08, P);
+		P->TwoChildThreePersProb				= Params::get_double(params, pre_params, adm_params, "Two Child Three Pers Prob", 0.11, P);
+		P->OnePersHouseProbOld					= Params::get_double(params, pre_params, adm_params, "One Pers House Prob Old", 0.5, P);
+		P->TwoPersHouseProbOld					= Params::get_double(params, pre_params, adm_params, "Two Pers House Prob Old", 0.5, P);
+		P->OnePersHouseProbYoung				= Params::get_double(params, pre_params, adm_params, "One Pers House Prob Young", 0.23, P);
+		P->TwoPersHouseProbYoung				= Params::get_double(params, pre_params, adm_params, "Two Pers House Prob Young", 0.23, P);
+		P->OneChildProbYoungestChildUnderFive	= Params::get_double(params, pre_params, adm_params, "One Child Prob Youngest Child Under Five", 0.5, P);
+		P->TwoChildrenProbYoungestUnderFive		= Params::get_double(params, pre_params, adm_params, "Two Children Prob Youngest Under Five", 0.0, P);
+		P->TwoChildrenProbYoungestUnderFive		= Params::get_double(params, pre_params, adm_params, "Prob Youngest Child Under Five", 0, P);
+		P->ZeroChildThreePersProb				= Params::get_double(params, pre_params, adm_params, "Zero Child Three Pers Prob", 0.25, P);
+		P->OneChildFourPersProb					= Params::get_double(params, pre_params, adm_params, "One Child Four Pers Prob", 0.2, P);
+		P->YoungAndSingleSlope					= Params::get_double(params, pre_params, adm_params, "Young And Single Slope", 0.7, P);
+		P->YoungAndSingle						= Params::get_int(params, pre_params, adm_params, "Young And Single", 36, P);
+		P->NoChildPersAge						= Params::get_int(params, pre_params, adm_params, "No Child Pers Age", 44, P);
+		P->OldPersAge							= Params::get_int(params, pre_params, adm_params, "Old Pers Age", 60, P);
+		P->ThreeChildFivePersProb				= Params::get_double(params, pre_params, adm_params, "Three Child Five Pers Prob", 0.5, P);
+		P->OlderGenGap							= Params::get_int(params, pre_params, adm_params, "Older Gen Gap", 19, P);
 	}
 
 	if (P->DoOneGen != 0) P->DoOneGen = 1;
@@ -2322,7 +2322,7 @@ void Params::ReadParams(std::string const& ParamFile, std::string const& PrePara
 		if (P->KeyWorkerProphCellIncThresh < 1) P->KeyWorkerProphCellIncThresh = 1;
 	*/
 
-	//// Make unsigned short versions of various intervention variables. And scaled them by number of timesteps per day
+	//// Make unsigned short versions of various intervention variables. And scale them by number of timesteps per day
 	P->usHQuarantineHouseDuration = ((unsigned short int) (P->HQuarantineHouseDuration * P->TimeStepsPerDay));
 	P->usVaccTimeToEfficacy = ((unsigned short int) (P->VaccTimeToEfficacy * P->TimeStepsPerDay));
 	P->usVaccTimeEfficacySwitch = ((unsigned short int) (P->VaccTimeEfficacySwitch * P->TimeStepsPerDay));
@@ -2344,14 +2344,12 @@ void Params::ReadParams(std::string const& ParamFile, std::string const& PrePara
 			P->cosx[i] = cos(t);
 		}
 	}
-	if (P->R0scale != 1.0)
+	if (P->R0scale != 1.0) // if scaling R0, scale household, place and spatial betas
 	{
-		P->HouseholdTrans *= P->R0scale;
-		if (P->FixLocalBeta) P->LocalBeta *= P->R0scale;
-		P->R0 *= P->R0scale;
-		for (int place_type = 0; place_type < P->NumPlaceTypes; place_type++) {
-			P->PlaceTypeTrans[place_type] *= P->R0scale;
-		}
+		P->HouseholdTrans *= P->R0scale; // household
+		if (P->FixLocalBeta) P->LocalBeta *= P->R0scale; // spatial
+		P->R0 *= P->R0scale; // R0
+		for (int place_type = 0; place_type < P->NumPlaceTypes; place_type++)	P->PlaceTypeTrans[place_type] *= P->R0scale; // place
 		Files::xfprintf_stderr("Rescaled transmission coefficients by factor of %lg\n", P->R0scale);
 	}
 	Files::xfprintf_stderr("Parameters read\n");