From de4936c3921e3068dba29d2d21d6beccec2afa0a Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Thu, 30 Jun 2022 16:17:32 +0100 Subject: [PATCH 01/20] Betas array added --- src/CovidSim.cpp | 38 +++++++++++++++++++++++++--- src/Models/Cell.h | 4 ++- src/Param.h | 18 +++++++++++--- src/ReadParams.cpp | 62 ++++++++++++++++++++++------------------------ 4 files changed, 82 insertions(+), 40 deletions(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index b480596a0..0a85843ca 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -67,6 +67,34 @@ void CalibrationThresholdCheck(double, int); void CalcLikelihood(int, std::string const&, std::string const&); void CalcOriginDestMatrix_adunit(void); //added function to calculate origin destination matrix: ggilani 28/01/15 +void AllocateMemForBetasArray() // called in main (once per fitting run) +{ + P.Betas = new double** [P.SimulationDuration](); + for (int Day = 0; Day < P.SimulationDuration; Day++) + { + P.Betas[Day] = new double* [P.NumAdunits](); + for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) P.Betas[Day][AdUnit] = new double[P.NumInfectionSettings](); + } +} +void InitBetasArray() // called in InitModel (every realistaion/parameter guess iteration). +{ + if (P.VaryBetasOverTimeByRegion == false) + { + //// if not varying by region or time, assign them to be the same for every simulation day and every admin unit. + for (int Day = 0; Day < P.SimulationDuration; Day++) + for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) + { + // place (by type) + for (int PlaceType = 0; PlaceType < P.NumPlaceTypes; PlaceType++) + P.Betas[Day][AdUnit][PlaceType] = P.PlaceTypeTrans[PlaceType]; + // Household + P.Betas[Day][AdUnit][House] = P.HouseholdTrans; + // Spatial/Community + P.Betas[Day][AdUnit][Spatial] = P.LocalBeta; + } + } +} + ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** GLOBAL VARIABLES (some structures in CovidSim.h file and some containers) - memory allocated later. ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// @@ -261,9 +289,7 @@ int main(int argc, char* argv[]) P.NumThreads = 1; #endif if (pre_param_file.empty()) - { pre_param_file = std::string(".." DIRECTORY_SEPARATOR "Pre_") + param_file; - } //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** READ IN PARAMETERS, DATA ETC. @@ -287,15 +313,17 @@ int main(int argc, char* argv[]) for (int i = 0; i < MAX_ADUNITS; i++) AdUnits[i].NI = 0; for (auto const& int_file : InterventionFiles) ReadInterventions(int_file); - Files::xfprintf_stderr("Model setup in %lf seconds\n", ((double) clock() - cl) / CLOCKS_PER_SEC); // Allocate memory for Efficacies array P.NumInfectionSettings = MAX_NUM_PLACE_TYPES + 2; // Maximum number of place types, plus household, and spatial - P.NumInterventionClasses = 6; // CI, HQ, PC, SD, Enhanced Social Distancing, DCT + P.NumInterventionClasses = 6; // CI, HQ, PC, SD, Enhanced Social Distancing, DCT P.Efficacies = new double*[P.NumInterventionClasses](); for (int intervention = 0; intervention < P.NumInterventionClasses; intervention++) P.Efficacies[intervention] = new double[P.NumInfectionSettings](); + // Allocate memory for Betas array (dynamic allocation must be done after P.NumAdunits set in SetupModel.cpp::SetupPopulation) + AllocateMemForBetasArray(); + //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** RUN MODEL @@ -1180,6 +1208,8 @@ void InitModel(int run) // passing run number so we can save run number in the i //// Add all of the above to P.Efficacies array. UpdateEfficacyArray(); + //// InitializeBetasArray (either to be same for all days and regions), or otherwise depending on how we are modelling them). Memory allocated in main using + InitBetasArray(); // Initialize CFR scalings P.CFR_Critical_Scale_Current = P.CFR_TimeScaling_Critical[0]; diff --git a/src/Models/Cell.h b/src/Models/Cell.h index 06be1026f..f1de24c37 100644 --- a/src/Models/Cell.h +++ b/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; diff --git a/src/Param.h b/src/Param.h index 4ed49471f..d6a49e903 100644 --- a/src/Param.h +++ b/src/Param.h @@ -114,7 +114,16 @@ struct Param 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 + */ + 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 +149,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 +159,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 +274,7 @@ struct Param int NumInterventionClasses = 0, NumInfectionSettings = 0; // initialized in CovidSim.cpp::main. + ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** VARIABLE EFFICACIES OVER TIME ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** diff --git a/src/ReadParams.cpp b/src/ReadParams.cpp index 580cb7f98..eacd59253 100644 --- a/src/ReadParams.cpp +++ b/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"); From f2e892cb13606ca99b3952288c64b944b75128c6 Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Fri, 1 Jul 2022 14:44:14 +0100 Subject: [PATCH 02/20] replaced HH beta with array in InfectSweep --- src/CovidSim.cpp | 12 ++-- src/Param.h | 3 +- src/Sweep.cpp | 141 ++++++++++++++++++++++++++--------------------- 3 files changed, 87 insertions(+), 69 deletions(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index 0a85843ca..9102af34f 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -69,8 +69,8 @@ void CalcOriginDestMatrix_adunit(void); //added function to calculate origin des void AllocateMemForBetasArray() // called in main (once per fitting run) { - P.Betas = new double** [P.SimulationDuration](); - for (int Day = 0; Day < P.SimulationDuration; Day++) + P.Betas = new double** [(int)P.SimulationDuration + 1](); + for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) { P.Betas[Day] = new double* [P.NumAdunits](); for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) P.Betas[Day][AdUnit] = new double[P.NumInfectionSettings](); @@ -78,10 +78,10 @@ void AllocateMemForBetasArray() // called in main (once per fitting run) } void InitBetasArray() // called in InitModel (every realistaion/parameter guess iteration). { - if (P.VaryBetasOverTimeByRegion == false) - { + //if (P.VaryBetasOverTimeByRegion == false) // while you're setting up, leave this condition out. Put back in later once mobility data included. + //{ //// if not varying by region or time, assign them to be the same for every simulation day and every admin unit. - for (int Day = 0; Day < P.SimulationDuration; Day++) + for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) { // place (by type) @@ -92,7 +92,7 @@ void InitBetasArray() // called in InitModel (every realistaion/parameter guess // Spatial/Community P.Betas[Day][AdUnit][Spatial] = P.LocalBeta; } - } + //} } ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// diff --git a/src/Param.h b/src/Param.h index d6a49e903..4b02b9d12 100644 --- a/src/Param.h +++ b/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,7 +109,7 @@ 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; diff --git a/src/Sweep.cpp b/src/Sweep.cpp index 9aef2fd73..2aa49c52c 100644 --- a/src/Sweep.cpp +++ b/src/Sweep.cpp @@ -292,10 +292,10 @@ void InfectSweep(double t, int run) // added run number as argument in order to int CellQueue; unsigned short int TimeStepNow = (unsigned short int) (P.TimeStepsPerDay * t); // TimeStepNow = the timestep number of the start of the current day + int Day = (int)floor(t); // integer value of current day (used for indexing). double fp = P.ModelTimeStep / (1 - P.FalsePositiveRate); // fp = false positive double seasonality = (P.DoSeasonality) ? (P.Seasonality[((int)t) % DAYS_PER_YEAR]) : 1.0; // if doing seasonality, pick seasonality from P.Seasonality array using day number in year. Otherwise set to 1. double SpatialSeasonal_Beta = seasonality * fp * P.LocalBeta; - double Household_Beta = (P.DoHouseholds) ? (seasonality * fp * P.HouseholdTrans) : 0; // if doing households, Household_Beta = seasonality * fp * P.HouseholdTrans, else Household_Beta = 0 // Establish if movement restrictions are in place on current day - store in BlanketMoveRestrInPlace, 0:false, 1:true int BlanketMoveRestrInPlace = ((P.DoBlanketMoveRestr) && (t >= P.MoveRestrTimeStart) && (t < P.MoveRestrTimeStart + P.MoveRestrDuration)); @@ -303,20 +303,30 @@ void InfectSweep(double t, int run) // added run number as argument in order to FILE* stderr_shared = stderr; #pragma omp parallel for private(CellQueue) schedule(static,1) default(none) \ - shared(t, P, CellLookup, Hosts, AdUnits, Households, Places, SamplingQueue, Cells, Mcells, StateT, Household_Beta, SpatialSeasonal_Beta, seasonality, TimeStepNow, fp, BlanketMoveRestrInPlace, stderr_shared) + shared(t, P, CellLookup, Hosts, AdUnits, Households, Places, SamplingQueue, Cells, Mcells, StateT, SpatialSeasonal_Beta, seasonality, TimeStepNow, fp, BlanketMoveRestrInPlace, stderr_shared) for (int ThreadNum = 0; ThreadNum < P.NumThreads; ThreadNum++) + { + Cell* ThisCell; + double Household_Beta; // if doing households, Household_Beta = seasonality * fp * P.HouseholdTrans, else Household_Beta = 0 + double SpatialInf_AllPeopleThisCell; ///// spatial infectiousness summed over all infectious people in loop below + for (int CellIndex = ThreadNum; CellIndex < P.NumPopulatedCells; CellIndex += P.NumThreads) //// loop over (in parallel) all populated cells. Loop 1) { - Cell* ThisCell = CellLookup[CellIndex]; // select Cell given by index CellIndex - double SpatialInf_AllPeopleThisCell = 0; ///// spatial infectiousness summed over all infectious people in loop below - + ThisCell = CellLookup[CellIndex]; // select Cell given by index CellIndex + SpatialInf_AllPeopleThisCell = 0; ///// spatial infectiousness summed over all infectious people in loop below + + //// quantities that will be used in loop below. + Person* InfectiousPerson; + int InfectiousPersonIndex; + bool DigiContactTrace_ThisPersonNow; + //// Loop over array of indices of infectious people ThisCell->I in cell ThisCell. Loop 1a) for (int InfectiousPersonIndex_ThisCell = 0; InfectiousPersonIndex_ThisCell < ThisCell->I; InfectiousPersonIndex_ThisCell++) { //// get InfectiousPersonIndex from InfectiousPersonIndex_ThisCell - int InfectiousPersonIndex = ThisCell->infected[InfectiousPersonIndex_ThisCell]; + InfectiousPersonIndex = ThisCell->infected[InfectiousPersonIndex_ThisCell]; //// get InfectiousPerson from Hosts (array of people) corresponding to InfectiousPersonIndex, using pointer arithmetic. - Person* InfectiousPerson = Hosts + InfectiousPersonIndex; + InfectiousPerson = Hosts + InfectiousPersonIndex; //evaluate flag for digital contact tracing (DigiContactTrace_ThisPersonNow) here at the beginning for each individual // DigiContactTrace_ThisPersonNow = 1 if: @@ -325,11 +335,20 @@ void InfectSweep(double t, int run) // added run number as argument in order to // AND Day number (t) is less than the end day for contact tracing in this administrative unit (ie. contact tracing has not ended) // AND the selected host is a digital contact tracing user // otherwise DigiContactTrace_ThisPersonNow = 0 - bool DigiContactTrace_ThisPersonNow = ((P.DoDigitalContactTracing) && (t >= AdUnits[Mcells[InfectiousPerson->mcell].adunit].DigitalContactTracingTimeStart) + DigiContactTrace_ThisPersonNow = ((P.DoDigitalContactTracing) && (t >= AdUnits[Mcells[InfectiousPerson->mcell].adunit].DigitalContactTracingTimeStart) && (t < AdUnits[Mcells[InfectiousPerson->mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[InfectiousPersonIndex].digitalContactTracingUser == 1)); // && (TimeStepNow <= (Hosts[InfectiousPersonIndex].detected_time + P.usCaseIsolationDelay))); // BEGIN HOUSEHOLD INFECTIONS - + + int AdminUnit = Mcells[Hosts[InfectiousPersonIndex].mcell].adunit; + + if (P.DoHouseholds) + { + Household_Beta = (P.VaryBetasOverTimeByRegion) ? P.Betas[Day][AdminUnit][House] : P.HouseholdTrans; + Household_Beta *= seasonality * fp; + } + else Household_Beta = 0; + if (Household_Beta > 0) { // For InfectiousPerson's household (InfectiousPerson->hh), @@ -353,14 +372,14 @@ void InfectSweep(double t, int run) // added run number as argument in order to // if individuals in the household are absent from places (ie. AtLeastOnePersonAbsent from test immediately above), scale up the infectiousness (Household_Infectiousness) of the household if (AtLeastOnePersonAbsent) { Household_Infectiousness *= P.Efficacies[PlaceClosure][House]; }/* NumPCD++;}*/ //// if people in your household are absent from places, person InfectiousPerson/InfectiousPersonIndex is more infectious to them, as they spend more time at home. - + // Loop over household members for (int HouseholdMember = FirstHouseholdMember; HouseholdMember < LastHouseholdMember; HouseholdMember++) //// loop over all people in household { if (Hosts[HouseholdMember].is_susceptible() && (!Hosts[HouseholdMember].Travelling)) //// if people in household uninfected/susceptible and not travelling { double Household_FOI = Household_Infectiousness * CalcHouseSusc(HouseholdMember, TimeStepNow, InfectiousPersonIndex); //// Household force of infection (FOI = infectiousness x susceptibility) from person InfectiousPersonIndex/InfectiousPerson on fellow household member - + // Force of Infection (Household_FOI) > random value between 0 and 1 if (ranf_mt(ThreadNum) < Household_FOI) { @@ -375,7 +394,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to } // if more than one person in household }// if Household_Beta > 0 // END HOUSHOLD INFECTIONS - + // BEGIN PLACE INFECTIONS if (P.DoPlaces) // if places functionality is enabled { @@ -404,33 +423,30 @@ void InfectSweep(double t, int run) // added run number as argument in order to } // else if movement restrictions in effect in either household microcell or place microcell else if ((Microcell_ThisPerson->moverest != Microcell_ThisPersonsPlaceLink->moverest) && ((Microcell_ThisPerson->moverest == TreatStat::Treated) || (Microcell_ThisPersonsPlaceLink->moverest == TreatStat::Treated))) - { - // multiply infectiousness of place by movement restriction effect - Place_Infectiousness *= P.MoveRestrEffect; - } - + Place_Infectiousness *= P.MoveRestrEffect; // multiply infectiousness of place by movement restriction effect + // BEGIN NON-HOTEL INFECTIONS - + // if linked place isn't a hotel and selected host isn't travelling if ((PlaceType != P.HotelPlaceType) && (!InfectiousPerson->Travelling)) { // PlaceGroupLink_index is index of group (of place type PlaceType) that selected host is linked to int PlaceGroupLink_index = (InfectiousPerson->PlaceGroupLinks[PlaceType]); - + // calculate infectiousness (PlaceInfectiousness_Scaled_DCT_copy) // which varies if contact tracing is in place // if contact tracing isn't in place PlaceInfectiousness_Scaled_DCT_copy is a copy of Place_Infectiousness // if contact tracing is in place, PlaceInfectiousness_Scaled_DCT_copy is Place_Infectiousness * P.ScalingFactorPlaceDigitalContacts // in either case PlaceInfectiousness_Scaled_DCT_copy is capped at 1 - + // if contact tracing double Place_Infectiousness_DCT_copy, PlaceInfectiousness_Scaled_DCT_copy; if (DigiContactTrace_ThisPersonNow) { // copy Place_Infectiousness - Place_Infectiousness_DCT_copy = Place_Infectiousness; + Place_Infectiousness_DCT_copy = Place_Infectiousness; // multiply Place_Infectiousness_DCT_copy by P.ScalingFactorPlaceDigitalContacts - PlaceInfectiousness_Scaled_DCT_copy = Place_Infectiousness_DCT_copy * P.ScalingFactorPlaceDigitalContacts; + PlaceInfectiousness_Scaled_DCT_copy = Place_Infectiousness_DCT_copy * P.ScalingFactorPlaceDigitalContacts; // cap Place_Infectiousness_DCT_copy at 1 if (Place_Infectiousness_DCT_copy > 1) Place_Infectiousness_DCT_copy = 1; // cap at 1 @@ -446,7 +462,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to } // if infectiousness is < 0, we have an error - end the program - int NumPotentialInfecteesPlaceGroup; + int NumPotentialInfecteesPlaceGroup; if (PlaceInfectiousness_Scaled_DCT_copy < 0) { Files::xfprintf(stderr_shared, "@@@ %lg\n", PlaceInfectiousness_Scaled_DCT_copy); @@ -462,14 +478,14 @@ void InfectSweep(double t, int run) // added run number as argument in order to { NumPotentialInfecteesPlaceGroup = (int)ignbin_mt((int32_t)Places[PlaceType][PlaceLink].group_size[PlaceGroupLink_index], PlaceInfectiousness_Scaled_DCT_copy, ThreadNum); } - + // if potential infectees > 0 - if (NumPotentialInfecteesPlaceGroup > 0) + if (NumPotentialInfecteesPlaceGroup > 0) { // pick NumPotentialInfecteesPlaceGroup members of place PlaceType,PlaceLink and add them to sampling queue for thread ThreadNum SampleWithoutReplacement(ThreadNum, NumPotentialInfecteesPlaceGroup, Places[PlaceType][PlaceLink].group_size[PlaceGroupLink_index]); //// changes thread-specific SamplingQueue. } - + // loop over sampling queue of potential infectees for (int PotentialInfectee_ThisPlaceGroupIndex = 0; PotentialInfectee_ThisPlaceGroupIndex < NumPotentialInfecteesPlaceGroup; PotentialInfectee_ThisPlaceGroupIndex++) { @@ -483,14 +499,14 @@ void InfectSweep(double t, int run) // added run number as argument in order to if ((PlaceType == P.CareHomePlaceType) && ((!Hosts[InfectiousPersonIndex].care_home_resident) || (!Hosts[PotentialInfectee_PlaceGroup].care_home_resident))) PlaceSusceptibility *= P.CareHomeWorkerGroupScaling; //these are all place group contacts to be tracked for digital contact tracing - add to StateT queue for contact tracing //if infectee is also a user, add them as a contact - + if ((DigiContactTrace_ThisPersonNow) && (Hosts[PotentialInfectee_PlaceGroup].digitalContactTracingUser) && (InfectiousPersonIndex != PotentialInfectee_PlaceGroup) && (!HOST_ABSENT(PotentialInfectee_PlaceGroup))) { // scale place susceptibility by proportion who self isolate and store as PlaceSusceptibility_DCT_scaled double PlaceSusceptibility_DCT_scaled = P.ProportionDigitalContactsIsolate * PlaceSusceptibility; // if random number < PlaceSusceptibility_DCT_scaled // AND number of contacts of InfectiousPersonIndex(!) is less than maximum digital contact to trace - if ((Hosts[InfectiousPersonIndex].ncontacts < P.MaxDigitalContactsToTrace) && (ranf_mt(ThreadNum) Travelling)) @@ -549,7 +565,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to // if contact tracing in place, multiply Place_Infectiousness_scaled = Place_Infectiousness*scalingfactor, otherwise Place_Infectiousness_scaled = Place_Infectiousness double Place_Infectiousness_scaled = (DigiContactTrace_ThisPersonNow) ? (Place_Infectiousness * P.ScalingFactorPlaceDigitalContacts) : Place_Infectiousness; // Place_Infectiousness_scaled shouldn't be less than 0 so generate error if it is - int NumPotentialInfecteesHotel; + int NumPotentialInfecteesHotel; if (Place_Infectiousness_scaled < 0) { ERR_CRITICAL_FMT("@@@ %lg\n", Place_Infectiousness); @@ -571,14 +587,14 @@ void InfectSweep(double t, int run) // added run number as argument in order to // calculate place susceptibility PlaceSusceptibility double PlaceSusceptibility = CalcPlaceSusc(PotentialInfectee_Hotel, PlaceType, TimeStepNow); // use group structure to model multiple care homes with shared staff - in which case residents of one "group" don't mix with those in another, only staff do. - if ((Hosts[InfectiousPersonIndex].care_home_resident) && (Hosts[PotentialInfectee_Hotel].care_home_resident) && (Hosts[InfectiousPersonIndex].PlaceGroupLinks[PlaceType]!= Hosts[PotentialInfectee_Hotel].PlaceGroupLinks[PlaceType])) PlaceSusceptibility *= P.CareHomeResidentPlaceScaling; + if ((Hosts[InfectiousPersonIndex].care_home_resident) && (Hosts[PotentialInfectee_Hotel].care_home_resident) && (Hosts[InfectiousPersonIndex].PlaceGroupLinks[PlaceType] != Hosts[PotentialInfectee_Hotel].PlaceGroupLinks[PlaceType])) PlaceSusceptibility *= P.CareHomeResidentPlaceScaling; // allow care home staff to have lowere contacts in care homes - to allow for PPE/environmental contamination. if ((PlaceType == P.CareHomePlaceType) && ((!Hosts[InfectiousPersonIndex].care_home_resident) || (!Hosts[PotentialInfectee_Hotel].care_home_resident))) PlaceSusceptibility *= P.CareHomeWorkerGroupScaling; - + //these are all place group contacts to be tracked for digital contact tracing - add to StateT queue for contact tracing //if infectee is also a user, add them as a contact - + // if contact tracing in place AND potential infectee PotentialInfectee_Hotel is a contact tracing user AND PotentialInfectee_Hotel isn't absent AND PotentialInfectee_Hotel isn't InfectiousPersonIndex (suspect this should be si) if ((DigiContactTrace_ThisPersonNow) && (Hosts[PotentialInfectee_Hotel].digitalContactTracingUser) && (InfectiousPersonIndex != PotentialInfectee_Hotel) && (!HOST_ABSENT(PotentialInfectee_Hotel))) @@ -604,7 +620,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to Microcell* MicroCell_PotentialInfectee_Hotel = Mcells + Hosts[PotentialInfectee_Hotel].mcell; //if doing digital contact tracing, scale down susceptibility here - PlaceSusceptibility *= CalcPersonSusc(PotentialInfectee_Hotel, TimeStepNow, InfectiousPersonIndex)*Place_Infectiousness/Place_Infectiousness_scaled; + PlaceSusceptibility *= CalcPersonSusc(PotentialInfectee_Hotel, TimeStepNow, InfectiousPersonIndex) * Place_Infectiousness / Place_Infectiousness_scaled; // if blanket movement restrictions are in place if (BlanketMoveRestrInPlace) { @@ -622,42 +638,42 @@ void InfectSweep(double t, int run) // added run number as argument in order to // multiply susceptibility by movement restriction effect PlaceSusceptibility *= P.MoveRestrEffect; } - + // ** do infections ** - + // is susceptibility is 1 (ie infect everyone) or random number is less than susceptibility if ((PlaceSusceptibility == 1) || (ranf_mt(ThreadNum) < PlaceSusceptibility)) { // explicitly cast to short to resolve level 4 warning const short int infect_type = static_cast (2 + PlaceType + MAX_NUM_PLACE_TYPES + INFECT_TYPE_MASK * (1 + InfectiousPerson->infect_type / INFECT_TYPE_MASK)); - - AddInfections(ThreadNum, Hosts[PotentialInfectee_Hotel].pcell% P.NumThreads, InfectiousPersonIndex, PotentialInfectee_Hotel, infect_type); + + AddInfections(ThreadNum, Hosts[PotentialInfectee_Hotel].pcell % P.NumThreads, InfectiousPersonIndex, PotentialInfectee_Hotel, infect_type); } // susceptibility test } // PotentialInfectee_Hotel uninfected and not absent. } // loop over sampling queue } // selected host InfectiousPerson is not travelling or selected link is to a hotel - + // ** END HOTEL INFECTIONS ** - + } // if place link relevant } // loop over place types } // if host isn't absent } // if places functionality enabled - + // END PLACE INFECTIONS - + // BEGIN SPATIAL INFECTIONS - + //// First determine spatial FOI component (SpatialInf_AllPeopleThisCell) - + // if seasonality beta > 0 // do spatial infections //// ie sum spatial infectiousness over all infected people, the infections from which are allocated after loop over infected people. - if (SpatialSeasonal_Beta > 0) + if (SpatialSeasonal_Beta > 0) { - double SpatialInf_ThisPerson; + double SpatialInf_ThisPerson; if (InfectiousPerson->Travelling) //// if host currently away from their cell, they cannot add to their cell's spatial infectiousness. - SpatialInf_ThisPerson = 0; + SpatialInf_ThisPerson = 0; else { // calculate SpatialInf_ThisPerson based on host and timestep @@ -689,8 +705,8 @@ void InfectSweep(double t, int run) // added run number as argument in order to } } } // loop over infectious people in cell - - + + //// Now allocate spatial infections using Force Of Infection (SpatialInf_AllPeopleThisCell) calculated above if (SpatialInf_AllPeopleThisCell > 0) //// if spatial infectiousness positive { @@ -698,7 +714,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to int NumPotentialCelltoCellInfections = (int)ignpoi_mt(SpatialInf_AllPeopleThisCell * SpatialSeasonal_Beta * ((double)ThisCell->tot_prob), ThreadNum); //// number people this cell's population might infect elsewhere. poisson random number based on spatial infectiousness s5, SpatialSeasonal_Beta (seasonality) and this cell's "probability" (guessing this is a function of its population and geographical size). // NumInfectiousThisCell = number of infectious people in cell ThisCell int NumInfectiousThisCell = ThisCell->I; - + if (NumPotentialCelltoCellInfections > 0) //// this block normalises cumulative infectiousness cell_inf by person. s5 is the total cumulative spatial infectiousness. Reason is so that infector can be chosen using ranf_mt, which returns random number between 0 and 1. { //// normalise by cumulative spatial infectiousness. @@ -706,17 +722,17 @@ void InfectSweep(double t, int run) // added run number as argument in order to //// does same as the above loop just a slightly faster calculation. i.e. StateT[ThreadNum].cell_inf[NumInfectiousThisCell - 1] / SpatialInf_AllPeopleThisCell would equal 1 or -1 anyway. StateT[ThreadNum].cell_inf[NumInfectiousThisCell - 1] = (StateT[ThreadNum].cell_inf[NumInfectiousThisCell - 1] < 0) ? -1.0f : 1.0f; } - + //// loop over infections to dole out. roughly speaking, this determines which infectious person in cell ThisCell infects which person elsewhere. - for (int Infection = 0; Infection < NumPotentialCelltoCellInfections; Infection++) + for (int Infection = 0; Infection < NumPotentialCelltoCellInfections; Infection++) { //// decide on infector PotentialInfector_Spatial/InfectiousPerson from cell ThisCell. int PotentialInfector_CellIndex; // PotentialInfector_CellIndex = index of infector - + if (NumInfectiousThisCell == 1) PotentialInfector_CellIndex = 0; // if only one infectious person in cell, infector index is first in cell (person 0) // if more than one infectious person in cell pick an infectious person (given by index PotentialInfector_CellIndex) //// roughly speaking, this determines which infectious person in cell ThisCell infects which person elsewhere - else + else { int StepSize; double RandomNum = ranf_mt(ThreadNum); ///// choose random number between 0 and 1 @@ -725,12 +741,12 @@ void InfectSweep(double t, int run) // added run number as argument in order to do { if (StepSize > 1) StepSize /= 2; //// amount StepSize to change PotentialInfector_CellIndex by reduced by half. Basically a binary search saying, keep amending potential infector PotentialInfector_CellIndex until either PotentialInfector_CellIndex less than zero or more than number of infected people until you find PotentialInfector_CellIndex s.t. spatial infectiousness "matches" RandomNum. - if ((PotentialInfector_CellIndex > 0) && (fabs(StateT[ThreadNum].cell_inf[PotentialInfector_CellIndex - 1]) >= RandomNum)) + if ((PotentialInfector_CellIndex > 0) && (fabs(StateT[ThreadNum].cell_inf[PotentialInfector_CellIndex - 1]) >= RandomNum)) { PotentialInfector_CellIndex -= StepSize; if (PotentialInfector_CellIndex == 0) KeepSearching = 0; } - else if ((PotentialInfector_CellIndex < NumInfectiousThisCell - 1) && (fabs(StateT[ThreadNum].cell_inf[PotentialInfector_CellIndex]) < RandomNum)) + else if ((PotentialInfector_CellIndex < NumInfectiousThisCell - 1) && (fabs(StateT[ThreadNum].cell_inf[PotentialInfector_CellIndex]) < RandomNum)) { PotentialInfector_CellIndex += StepSize; if (PotentialInfector_CellIndex == NumInfectiousThisCell - 1) KeepSearching = 0; @@ -749,7 +765,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to && (t < AdUnits[Mcells[PotentialInfector_Spatial->mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[PotentialInfector_Index].digitalContactTracingUser == 1)); // && (TimeStepNow <= (Hosts[PotentialInfector_Spatial].detected_time + P.usCaseIsolationDelay))); //// decide on infectee - + int KeepSearchingForCellToInfect = 1; // do the following while KeepSearchingForCellToInfect = 1 do { @@ -766,10 +782,10 @@ void InfectSweep(double t, int run) // added run number as argument in order to ///// pick random person SusceptiblePerson within susceptibles of cell ct (S0 initial number susceptibles within cell). int SusceptiblePerson = (int)(ranf_mt(ThreadNum) * ((double)ct->S0)); int PotentialInfectee_Spatial = ct->susceptible[SusceptiblePerson]; - + double PotentialInfectorInfecteeDistSquared = dist2(Hosts + PotentialInfectee_Spatial, Hosts + PotentialInfector_Index); /// calculate distance squared between PotentialInfectee_Spatial and PotentialInfector_Spatial double AcceptProb = P.KernelLookup.num(PotentialInfectorInfecteeDistSquared) / ThisCell->max_trans[l]; //// acceptance probability - + // initialise KeepSearchingForCellToInfect = 0 (KeepSearchingForCellToInfect = 1 is the while condition for this loop) KeepSearchingForCellToInfect = 0; // if random number greater than acceptance probablility or infectee is dead @@ -814,8 +830,8 @@ void InfectSweep(double t, int run) // added run number as argument in order to //scale down susceptibility so we don't over accept Spatial_Susc /= P.ScalingFactorSpatialDigitalContacts; } - - + + if (SusceptiblePerson < ct->S) // only consider susceptible people as possible infectees { Spatial_Susc *= CalcPersonSusc(PotentialInfectee_Spatial, TimeStepNow, PotentialInfector_Index); @@ -844,7 +860,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to { // explicitly cast to short to resolve level 4 warning const short int infect_type = static_cast(2 + 2 * MAX_NUM_PLACE_TYPES + INFECT_TYPE_MASK * (1 + PotentialInfector_Spatial->infect_type / INFECT_TYPE_MASK)); - + AddInfections(ThreadNum, CellQueue, PotentialInfector_Index, PotentialInfectee_Spatial, infect_type); } } @@ -855,6 +871,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to } // loop over infections doled out by cell } // SpatialInf_AllPeopleThisCell > 0 } + } #pragma omp parallel for schedule(static,1) default(none) \ From b2c4bc7fb06a13ff3f12d17a4b93a7449efef938 Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Fri, 1 Jul 2022 14:49:26 +0100 Subject: [PATCH 03/20] taken Place Betas out of CalcInfSusc.cpp::CalcPlaceInf and put them in Sweep.cpp::InfectSweep --- src/CalcInfSusc.cpp | 2 +- src/Sweep.cpp | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/src/CalcInfSusc.cpp b/src/CalcInfSusc.cpp index 44caf215b..4b6e1e386 100644 --- a/src/CalcInfSusc.cpp +++ b/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) diff --git a/src/Sweep.cpp b/src/Sweep.cpp index 2aa49c52c..bc8b4741f 100644 --- a/src/Sweep.cpp +++ b/src/Sweep.cpp @@ -409,7 +409,10 @@ void InfectSweep(double t, int run) // added run number as argument in order to if (PlaceLink >= 0) //// PlaceLink >= 0 means if place type PlaceType is relevant to InfectiousPerson. (Now allowing for partial attendance). { // Place_Infectiousness - double Place_Infectiousness = fp * seasonality * CalcPlaceInf(InfectiousPersonIndex, PlaceType, TimeStepNow); + double Place_Infectiousness = CalcPlaceInf(InfectiousPersonIndex, PlaceType, TimeStepNow); + Place_Infectiousness *= P.PlaceTypeTrans[PlaceType]; + Place_Infectiousness /= P.PlaceTypeGroupSizeParam1[PlaceType]; + Place_Infectiousness *= fp * seasonality; // select microcell of the place linked to InfectiousPerson with link PlaceLink Microcell* Microcell_ThisPersonsPlaceLink = Mcells + Places[PlaceType][PlaceLink].mcell; // if blanket movement restrictions are in place on current day @@ -423,7 +426,9 @@ void InfectSweep(double t, int run) // added run number as argument in order to } // else if movement restrictions in effect in either household microcell or place microcell else if ((Microcell_ThisPerson->moverest != Microcell_ThisPersonsPlaceLink->moverest) && ((Microcell_ThisPerson->moverest == TreatStat::Treated) || (Microcell_ThisPersonsPlaceLink->moverest == TreatStat::Treated))) + { Place_Infectiousness *= P.MoveRestrEffect; // multiply infectiousness of place by movement restriction effect + } // BEGIN NON-HOTEL INFECTIONS From 5ac4966a9b0fe068344e3e848e8e4237bbbb58d2 Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Fri, 1 Jul 2022 15:17:36 +0100 Subject: [PATCH 04/20] replaced place and spatial betas with Betas array in InfectSweep --- src/Sweep.cpp | 37 ++++++++++++++++--------------------- 1 file changed, 16 insertions(+), 21 deletions(-) diff --git a/src/Sweep.cpp b/src/Sweep.cpp index bc8b4741f..6ed6a4098 100644 --- a/src/Sweep.cpp +++ b/src/Sweep.cpp @@ -295,7 +295,6 @@ void InfectSweep(double t, int run) // added run number as argument in order to int Day = (int)floor(t); // integer value of current day (used for indexing). double fp = P.ModelTimeStep / (1 - P.FalsePositiveRate); // fp = false positive double seasonality = (P.DoSeasonality) ? (P.Seasonality[((int)t) % DAYS_PER_YEAR]) : 1.0; // if doing seasonality, pick seasonality from P.Seasonality array using day number in year. Otherwise set to 1. - double SpatialSeasonal_Beta = seasonality * fp * P.LocalBeta; // Establish if movement restrictions are in place on current day - store in BlanketMoveRestrInPlace, 0:false, 1:true int BlanketMoveRestrInPlace = ((P.DoBlanketMoveRestr) && (t >= P.MoveRestrTimeStart) && (t < P.MoveRestrTimeStart + P.MoveRestrDuration)); @@ -303,11 +302,12 @@ void InfectSweep(double t, int run) // added run number as argument in order to FILE* stderr_shared = stderr; #pragma omp parallel for private(CellQueue) schedule(static,1) default(none) \ - shared(t, P, CellLookup, Hosts, AdUnits, Households, Places, SamplingQueue, Cells, Mcells, StateT, SpatialSeasonal_Beta, seasonality, TimeStepNow, fp, BlanketMoveRestrInPlace, stderr_shared) + shared(t, P, CellLookup, Hosts, AdUnits, Households, Places, SamplingQueue, Cells, Mcells, StateT, seasonality, TimeStepNow, fp, BlanketMoveRestrInPlace, stderr_shared) for (int ThreadNum = 0; ThreadNum < P.NumThreads; ThreadNum++) { Cell* ThisCell; double Household_Beta; // if doing households, Household_Beta = seasonality * fp * P.HouseholdTrans, else Household_Beta = 0 + double SpatialSeasonal_Beta; double SpatialInf_AllPeopleThisCell; ///// spatial infectiousness summed over all infectious people in loop below for (int CellIndex = ThreadNum; CellIndex < P.NumPopulatedCells; CellIndex += P.NumThreads) //// loop over (in parallel) all populated cells. Loop 1) @@ -327,6 +327,8 @@ void InfectSweep(double t, int run) // added run number as argument in order to InfectiousPersonIndex = ThisCell->infected[InfectiousPersonIndex_ThisCell]; //// get InfectiousPerson from Hosts (array of people) corresponding to InfectiousPersonIndex, using pointer arithmetic. InfectiousPerson = Hosts + InfectiousPersonIndex; + int AdUnit_ThisPerson = Mcells[InfectiousPerson->mcell].adunit; + //evaluate flag for digital contact tracing (DigiContactTrace_ThisPersonNow) here at the beginning for each individual // DigiContactTrace_ThisPersonNow = 1 if: @@ -335,20 +337,15 @@ void InfectSweep(double t, int run) // added run number as argument in order to // AND Day number (t) is less than the end day for contact tracing in this administrative unit (ie. contact tracing has not ended) // AND the selected host is a digital contact tracing user // otherwise DigiContactTrace_ThisPersonNow = 0 - DigiContactTrace_ThisPersonNow = ((P.DoDigitalContactTracing) && (t >= AdUnits[Mcells[InfectiousPerson->mcell].adunit].DigitalContactTracingTimeStart) - && (t < AdUnits[Mcells[InfectiousPerson->mcell].adunit].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && (Hosts[InfectiousPersonIndex].digitalContactTracingUser == 1)); // && (TimeStepNow <= (Hosts[InfectiousPersonIndex].detected_time + P.usCaseIsolationDelay))); + DigiContactTrace_ThisPersonNow = ((P.DoDigitalContactTracing) && + (t >= AdUnits[AdUnit_ThisPerson].DigitalContactTracingTimeStart) && + (t < AdUnits[AdUnit_ThisPerson].DigitalContactTracingTimeStart + P.DigitalContactTracingPolicyDuration) && + (Hosts[InfectiousPersonIndex].digitalContactTracingUser == 1)); // && (TimeStepNow <= (Hosts[InfectiousPersonIndex].detected_time + P.usCaseIsolationDelay))); // BEGIN HOUSEHOLD INFECTIONS - int AdminUnit = Mcells[Hosts[InfectiousPersonIndex].mcell].adunit; - - if (P.DoHouseholds) - { - Household_Beta = (P.VaryBetasOverTimeByRegion) ? P.Betas[Day][AdminUnit][House] : P.HouseholdTrans; - Household_Beta *= seasonality * fp; - } - else Household_Beta = 0; + Household_Beta = (P.DoHouseholds) ? P.Betas[Day][AdUnit_ThisPerson][House] * seasonality * fp : 0; if (Household_Beta > 0) { // For InfectiousPerson's household (InfectiousPerson->hh), @@ -410,9 +407,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to { // Place_Infectiousness double Place_Infectiousness = CalcPlaceInf(InfectiousPersonIndex, PlaceType, TimeStepNow); - Place_Infectiousness *= P.PlaceTypeTrans[PlaceType]; - Place_Infectiousness /= P.PlaceTypeGroupSizeParam1[PlaceType]; - Place_Infectiousness *= fp * seasonality; + Place_Infectiousness *= P.Betas[Day][AdUnit_ThisPerson][PlaceType] * fp * seasonality / P.PlaceTypeGroupSizeParam1[PlaceType]; // select microcell of the place linked to InfectiousPerson with link PlaceLink Microcell* Microcell_ThisPersonsPlaceLink = Mcells + Places[PlaceType][PlaceLink].mcell; // if blanket movement restrictions are in place on current day @@ -514,16 +509,14 @@ void InfectSweep(double t, int run) // added run number as argument in order to if ((Hosts[InfectiousPersonIndex].ncontacts < P.MaxDigitalContactsToTrace) && (ranf_mt(ThreadNum) < PlaceSusceptibility_DCT_scaled)) { Hosts[InfectiousPersonIndex].ncontacts++; //add to number of contacts made - int AdminUnit = Mcells[Hosts[PotentialInfectee_PlaceGroup].mcell].adunit; - if ((StateT[ThreadNum].ndct_queue[AdminUnit] < AdUnits[AdminUnit].n)) + int AdUnit = Mcells[Hosts[PotentialInfectee_PlaceGroup].mcell].adunit; + if ((StateT[ThreadNum].ndct_queue[AdUnit] < AdUnits[AdUnit].n)) { //find adunit for contact and add both contact and infectious host to lists - storing both so I can set times later. - StateT[ThreadNum].dct_queue[AdminUnit][StateT[ThreadNum].ndct_queue[AdminUnit]++] = { PotentialInfectee_PlaceGroup, InfectiousPersonIndex, TimeStepNow }; + StateT[ThreadNum].dct_queue[AdUnit][StateT[ThreadNum].ndct_queue[AdUnit]++] = { PotentialInfectee_PlaceGroup, InfectiousPersonIndex, TimeStepNow }; } else - { - Files::xfprintf(stderr_shared, "No more space in queue! Thread: %i, AdUnit: %i\n", ThreadNum, AdminUnit); - } + Files::xfprintf(stderr_shared, "No more space in queue! Thread: %i, AdUnit: %i\n", ThreadNum, AdUnit); } } @@ -674,6 +667,8 @@ void InfectSweep(double t, int run) // added run number as argument in order to // if seasonality beta > 0 // do spatial infections //// ie sum spatial infectiousness over all infected people, the infections from which are allocated after loop over infected people. + + SpatialSeasonal_Beta = seasonality * fp * P.Betas[Day][AdUnit_ThisPerson][Spatial]; if (SpatialSeasonal_Beta > 0) { double SpatialInf_ThisPerson; From 82cefc655e61a4b3f96257824a881f7a867bc724 Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Fri, 1 Jul 2022 15:22:37 +0100 Subject: [PATCH 05/20] moved Beta and Efficacy array allocation a little earlier in CovidSim.cpp --- src/CovidSim.cpp | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index 9102af34f..d74c5caa8 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -309,20 +309,22 @@ int main(int argc, char* argv[]) ///// initialize model (for all realisations). SetupModel(density_file, out_density_file, load_network_file, save_network_file, school_file, reg_demog_file, output_file_base); + + // Allocate memory for Efficacies array + P.NumInfectionSettings = MAX_NUM_PLACE_TYPES + 2; // Maximum number of place types, plus household, and spatial + P.NumInterventionClasses = 6; // CI, HQ, PC, SD, Enhanced Social Distancing, DCT + P.Efficacies = new double* [P.NumInterventionClasses](); + for (int intervention = 0; intervention < P.NumInterventionClasses; intervention++) P.Efficacies[intervention] = new double[P.NumInfectionSettings](); + + // Allocate memory for Betas array (dynamic allocation must be done after P.NumAdunits set in SetupModel.cpp::SetupPopulation) + AllocateMemForBetasArray(); + InitTransmissionCoeffs(); for (int i = 0; i < MAX_ADUNITS; i++) AdUnits[i].NI = 0; for (auto const& int_file : InterventionFiles) ReadInterventions(int_file); Files::xfprintf_stderr("Model setup in %lf seconds\n", ((double) clock() - cl) / CLOCKS_PER_SEC); - // Allocate memory for Efficacies array - P.NumInfectionSettings = MAX_NUM_PLACE_TYPES + 2; // Maximum number of place types, plus household, and spatial - P.NumInterventionClasses = 6; // CI, HQ, PC, SD, Enhanced Social Distancing, DCT - P.Efficacies = new double*[P.NumInterventionClasses](); - for (int intervention = 0; intervention < P.NumInterventionClasses; intervention++) P.Efficacies[intervention] = new double[P.NumInfectionSettings](); - - // Allocate memory for Betas array (dynamic allocation must be done after P.NumAdunits set in SetupModel.cpp::SetupPopulation) - AllocateMemForBetasArray(); //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** From d436208504e34f884f7f6b49eb47852ca0b9be8d Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Fri, 1 Jul 2022 15:26:29 +0100 Subject: [PATCH 06/20] replace Person with person in InitTransmissionCoeffs to avoid ambiguity --- src/SetupModel.cpp | 80 ++++++++++++++++++++++++---------------------- 1 file changed, 41 insertions(+), 39 deletions(-) diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp index 3372a6bef..9e0e2aff8 100644 --- a/src/SetupModel.cpp +++ b/src/SetupModel.cpp @@ -666,30 +666,32 @@ void InitTransmissionCoeffs(void) #pragma omp parallel for private(Household_Infectiousness,Spatial_Infectiousness,quantile,LatentToSympDelay,ProbSurvive) schedule(static,1) reduction(+:HH_Infections,SpatialInfections,HH_SAR_Denom) default(none) shared(P, Households, Hosts, Mcells) for (int Thread = 0; Thread < P.NumThreads; Thread++) // loop over threads { - for (int Person = Thread; Person < P.PopSize; Person += P.NumThreads) // loop over people + for (int person = Thread; person < P.PopSize; person += P.NumThreads) // loop over people { + //int AdUnit = Mcells[person->mcell].adunit; + // assign susceptibility of each host. if (P.SusceptibilitySD == 0) - Hosts[Person].susc = (float)((P.DoPartialImmunity) ? (1.0 - P.InitialImmunity[HOST_AGE_GROUP(Person)]) : 1.0); + Hosts[person].susc = (float)((P.DoPartialImmunity) ? (1.0 - P.InitialImmunity[HOST_AGE_GROUP(person)]) : 1.0); else - Hosts[Person].susc = (float)(((P.DoPartialImmunity) ? (1.0 - P.InitialImmunity[HOST_AGE_GROUP(Person)]) : 1.0) * gen_gamma_mt(1 / (P.SusceptibilitySD * P.SusceptibilitySD), 1 / (P.SusceptibilitySD * P.SusceptibilitySD), Thread)); + Hosts[person].susc = (float)(((P.DoPartialImmunity) ? (1.0 - P.InitialImmunity[HOST_AGE_GROUP(person)]) : 1.0) * gen_gamma_mt(1 / (P.SusceptibilitySD * P.SusceptibilitySD), 1 / (P.SusceptibilitySD * P.SusceptibilitySD), Thread)); // assign infectiousness of each host. if (P.InfectiousnessSD == 0) - Hosts[Person].infectiousness = (float)P.AgeInfectiousness[HOST_AGE_GROUP(Person)]; + Hosts[person].infectiousness = (float)P.AgeInfectiousness[HOST_AGE_GROUP(person)]; else - Hosts[Person].infectiousness = (float)(P.AgeInfectiousness[HOST_AGE_GROUP(Person)] * gen_gamma_mt(1 / (P.InfectiousnessSD * P.InfectiousnessSD), 1 / (P.InfectiousnessSD * P.InfectiousnessSD), Thread)); + Hosts[person].infectiousness = (float)(P.AgeInfectiousness[HOST_AGE_GROUP(person)] * gen_gamma_mt(1 / (P.InfectiousnessSD * P.InfectiousnessSD), 1 / (P.InfectiousnessSD * P.InfectiousnessSD), Thread)); // scale infectiousness by symptomatic or asymptomatic multiplier - if (ranf_mt(Thread) < P.ProportionSymptomatic[HOST_AGE_GROUP(Person)]) // if symptomatic, scale by Symptomatic Infectiousness (and make negative)... - Hosts[Person].infectiousness *= (float)(-P.SymptInfectiousness); + if (ranf_mt(Thread) < P.ProportionSymptomatic[HOST_AGE_GROUP(person)]) // if symptomatic, scale by Symptomatic Infectiousness (and make negative)... + Hosts[person].infectiousness *= (float)(-P.SymptInfectiousness); else // ... or if asymptomatic - Hosts[Person].infectiousness *= (float)P.AsymptInfectiousness; + Hosts[person].infectiousness *= (float)P.AsymptInfectiousness; // choose recovery_or_death_time from infectious period quantiles (inverse cumulative distribution function). Will reset this later for each person in Update::DoIncub. int j = (int)floor((quantile = ranf_mt(Thread) * CDF_RES)); quantile -= ((double)j); - Hosts[Person].recovery_or_death_time = (unsigned short int) floor(0.5 - (P.InfectiousPeriod * log(quantile * P.infectious_icdf[j + 1] + (1.0 - quantile) * P.infectious_icdf[j]) / P.ModelTimeStep)); + Hosts[person].recovery_or_death_time = (unsigned short int) floor(0.5 - (P.InfectiousPeriod * log(quantile * P.infectious_icdf[j + 1] + (1.0 - quantile) * P.infectious_icdf[j]) / P.ModelTimeStep)); // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** Household Infections @@ -698,14 +700,14 @@ void InitTransmissionCoeffs(void) { // choose multiplier of infectiousness if (P.NoInfectiousnessSDinHH) - Household_Infectiousness = ((Hosts[Person].infectiousness < 0) ? P.SymptInfectiousness : P.AsymptInfectiousness); + Household_Infectiousness = ((Hosts[person].infectiousness < 0) ? P.SymptInfectiousness : P.AsymptInfectiousness); else - Household_Infectiousness = fabs(Hosts[Person].infectiousness); + Household_Infectiousness = fabs(Hosts[person].infectiousness); // Care home residents less likely to infect via "household" contacts. - if (Hosts[Person].care_home_resident) Household_Infectiousness *= P.CareHomeResidentHouseholdScaling; - Household_Infectiousness *= P.ModelTimeStep * P.HouseholdTrans * P.HouseholdDenomLookup[Households[Hosts[Person].hh].nhr - 1]; + if (Hosts[person].care_home_resident) Household_Infectiousness *= P.CareHomeResidentHouseholdScaling; + Household_Infectiousness *= P.ModelTimeStep * P.HouseholdTrans * P.HouseholdDenomLookup[Households[Hosts[person].hh].nhr - 1]; ProbSurvive = 1.0; - for (int InfectiousDay = 0; InfectiousDay < (int)Hosts[Person].recovery_or_death_time; InfectiousDay++) // loop over days adding to force of infection, probability that other household members will be infected. + for (int InfectiousDay = 0; InfectiousDay < (int)Hosts[person].recovery_or_death_time; InfectiousDay++) // loop over days adding to force of infection, probability that other household members will be infected. { double ProbSurviveToday = 1.0 - Household_Infectiousness * P.infectiousness[InfectiousDay]; ProbSurvive *= ((ProbSurviveToday < 0) ? 0 : ProbSurviveToday); @@ -713,21 +715,21 @@ void InitTransmissionCoeffs(void) // loop over people in households. If household member susceptible (they will be unless already infected in this code block), // and ensuring person doesn't infect themselves, add to household infections, taking account of their age and whether they're a care home resident, - // Person.e. the usual stuff in CalcInfSusc.cpp, but without interventions - for (int HouseholdMember = Households[Hosts[Person].hh].FirstPerson; HouseholdMember < Households[Hosts[Person].hh].FirstPerson + Households[Hosts[Person].hh].nh; HouseholdMember++) - if ((Hosts[HouseholdMember].is_susceptible()) && (HouseholdMember != Person)) - HH_Infections += (1 - ProbSurvive) * P.AgeSusceptibility[HOST_AGE_GROUP(Person)] * ((Hosts[HouseholdMember].care_home_resident) ? P.CareHomeResidentHouseholdScaling : 1.0); - HH_SAR_Denom += (double)(Households[Hosts[Person].hh].nhr - 1); // add to household denominator + // person.e. the usual stuff in CalcInfSusc.cpp, but without interventions + for (int HouseholdMember = Households[Hosts[person].hh].FirstPerson; HouseholdMember < Households[Hosts[person].hh].FirstPerson + Households[Hosts[person].hh].nh; HouseholdMember++) + if ((Hosts[HouseholdMember].is_susceptible()) && (HouseholdMember != person)) + HH_Infections += (1 - ProbSurvive) * P.AgeSusceptibility[HOST_AGE_GROUP(person)] * ((Hosts[HouseholdMember].care_home_resident) ? P.CareHomeResidentHouseholdScaling : 1.0); + HH_SAR_Denom += (double)(Households[Hosts[person].hh].nhr - 1); // add to household denominator } // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** Spatial Infections // Sum over number of days until recovery time, in two parts: entire infection and after symptoms occur, as spatial contact rate differs between these periods. - LatentToSympDelay = (P.LatentToSymptDelay > Hosts[Person].recovery_or_death_time * P.ModelTimeStep) ? Hosts[Person].recovery_or_death_time * P.ModelTimeStep : P.LatentToSymptDelay; + LatentToSympDelay = (P.LatentToSymptDelay > Hosts[person].recovery_or_death_time * P.ModelTimeStep) ? Hosts[person].recovery_or_death_time * P.ModelTimeStep : P.LatentToSymptDelay; // Care home residents less likely to infect via "spatial" contacts. This doesn't correct for non care home residents being less likely to infect care home residents, // but since the latter are a small proportion of the population, this is a minor issue - Spatial_Infectiousness = fabs(Hosts[Person].infectiousness) * P.RelativeSpatialContact[HOST_AGE_GROUP(Person)] * ((Hosts[Person].care_home_resident) ? P.CareHomeResidentSpatialScaling : 1.0) * P.ModelTimeStep; + Spatial_Infectiousness = fabs(Hosts[person].infectiousness) * P.RelativeSpatialContact[HOST_AGE_GROUP(person)] * ((Hosts[person].care_home_resident) ? P.CareHomeResidentSpatialScaling : 1.0) * P.ModelTimeStep; if (P.Got_WAIFW_Matrix_Spatial) { // if doing contact matrices, then need to scale spatial infections of this person (infector) to various infectees. @@ -735,8 +737,8 @@ void InitTransmissionCoeffs(void) double AvContactRate_Infector = 0; // sum over "infectee" age groups to scale infectiousness of this infector. for (int InfecteeAge = 0; InfecteeAge < NUM_AGE_GROUPS; InfecteeAge++) - //AvContactRate_Infector += P.PropAgeGroup[0][InfecteeAge] * P.WAIFW_Matrix_SpatialOnly[InfecteeAge][HOST_AGE_GROUP(Person)]; // use index 0 for admin unit if doing whole country. - AvContactRate_Infector += P.PropAgeGroup[Mcells[Hosts[Person].mcell].adunit][InfecteeAge] * P.WAIFW_Matrix_SpatialOnly[InfecteeAge][HOST_AGE_GROUP(Person)]; + //AvContactRate_Infector += P.PropAgeGroup[0][InfecteeAge] * P.WAIFW_Matrix_SpatialOnly[InfecteeAge][HOST_AGE_GROUP(person)]; // use index 0 for admin unit if doing whole country. + AvContactRate_Infector += P.PropAgeGroup[Mcells[Hosts[person].mcell].adunit][InfecteeAge] * P.WAIFW_Matrix_SpatialOnly[InfecteeAge][HOST_AGE_GROUP(person)]; // scale spatial infectiousness by weighted average. Spatial_Infectiousness *= AvContactRate_Infector; } @@ -744,9 +746,9 @@ void InitTransmissionCoeffs(void) int InfectiousDay; /// Add to spatial infections from all days where latent but not symptomatic for (InfectiousDay = 0; InfectiousDay < NumDaysInfectiousNotSymptomatic; InfectiousDay++) SpatialInfections += Spatial_Infectiousness * P.infectiousness[InfectiousDay]; - Spatial_Infectiousness *= ((Hosts[Person].infectiousness < 0) ? P.SymptSpatialContactRate : 1); + Spatial_Infectiousness *= ((Hosts[person].infectiousness < 0) ? P.SymptSpatialContactRate : 1); /// Add to spatial infections when symptomatic (note do not initialize InfectiousDay again.) - for (; InfectiousDay < (int)Hosts[Person].recovery_or_death_time; InfectiousDay++) SpatialInfections += Spatial_Infectiousness * P.infectiousness[InfectiousDay]; + for (; InfectiousDay < (int)Hosts[person].recovery_or_death_time; InfectiousDay++) SpatialInfections += Spatial_Infectiousness * P.infectiousness[InfectiousDay]; } } // Divide total spatial infections by PopSize to get Spatial R0. @@ -765,13 +767,13 @@ void InitTransmissionCoeffs(void) { double PreviousTotalPlaceInfections = PlaceInfections; #pragma omp parallel for private(ProbSurvive,LatentToSympDelay,Place_Infectiousness,NumPeopleInPlaceGroup,ProbSurviveNonCareHome) schedule(static,1000) reduction(+:PlaceInfections) default(none) shared(P, Hosts, Places, PlaceType) - for (int Person = 0; Person < P.PopSize; Person++) + for (int person = 0; person < P.PopSize; person++) { - int PlaceNum = Hosts[Person].PlaceLinks[PlaceType]; + int PlaceNum = Hosts[person].PlaceLinks[PlaceType]; if (PlaceNum >= 0) //// i.e. if person has a link to a particular Place of this PlaceType. { - LatentToSympDelay = (P.LatentToSymptDelay > Hosts[Person].recovery_or_death_time * P.ModelTimeStep) ? Hosts[Person].recovery_or_death_time * P.ModelTimeStep : P.LatentToSymptDelay; - Place_Infectiousness = fabs(Hosts[Person].infectiousness) * P.ModelTimeStep * P.PlaceTypeTrans[PlaceType]; + LatentToSympDelay = (P.LatentToSymptDelay > Hosts[person].recovery_or_death_time * P.ModelTimeStep) ? Hosts[person].recovery_or_death_time * P.ModelTimeStep : P.LatentToSymptDelay; + Place_Infectiousness = fabs(Hosts[person].infectiousness) * P.ModelTimeStep * P.PlaceTypeTrans[PlaceType]; double PlaceInf_Scaled = Place_Infectiousness / P.PlaceTypeGroupSizeParam1[PlaceType]; ProbSurvive = 1.0; int NumDaysInfectiousNotSymptomatic = (int)(LatentToSympDelay / P.ModelTimeStep); @@ -781,10 +783,10 @@ void InitTransmissionCoeffs(void) double ProbSurviveToday = 1.0 - PlaceInf_Scaled * P.infectiousness[InfectiousDay]; ProbSurvive *= ((ProbSurviveToday < 0) ? 0 : ProbSurviveToday); } - NumPeopleInPlaceGroup = ((double)(_I64(Places[PlaceType][PlaceNum].group_size[Hosts[Person].PlaceGroupLinks[PlaceType]]) - 1)); - PlaceInf_Scaled *= (((Hosts[Person].infectiousness < 0) && (!Hosts[Person].care_home_resident)) ? // if person symptomatic and not a care home resident + NumPeopleInPlaceGroup = ((double)(_I64(Places[PlaceType][PlaceNum].group_size[Hosts[person].PlaceGroupLinks[PlaceType]]) - 1)); + PlaceInf_Scaled *= (((Hosts[person].infectiousness < 0) && (!Hosts[person].care_home_resident)) ? // if person symptomatic and not a care home resident (P.SymptPlaceTypeContactRate[PlaceType] * (1 - P.SymptPlaceTypeWithdrawalProp[PlaceType])) : 1); - for (; InfectiousDay < (int)Hosts[Person].recovery_or_death_time; InfectiousDay++) + for (; InfectiousDay < (int)Hosts[person].recovery_or_death_time; InfectiousDay++) { double ProbSurviveToday = 1.0 - PlaceInf_Scaled * P.infectiousness[InfectiousDay]; ProbSurvive *= ((ProbSurviveToday < 0) ? 0 : ProbSurviveToday); @@ -795,7 +797,7 @@ void InitTransmissionCoeffs(void) // use group structure to model multiple care homes with shared staff - in which case residents of one "group" don't mix with those in another, only staff do. // calculation uses average proportion of care home "members" who are residents. - if (Hosts[Person].care_home_resident) + if (Hosts[person].care_home_resident) PlaceInf_Scaled *= (1.0 - P.CareHomePropResidents) + P.CareHomePropResidents * (P.CareHomeWorkerGroupScaling * (((double)Places[PlaceType][PlaceNum].n - 1) - NumPeopleInPlaceGroup) + NumPeopleInPlaceGroup) / ((double)Places[PlaceType][PlaceNum].n - 1); ProbSurvive = 1.0; for (InfectiousDay = 0; InfectiousDay < NumDaysInfectiousNotSymptomatic; InfectiousDay++) @@ -803,8 +805,8 @@ void InitTransmissionCoeffs(void) double ProbSurviveToday = 1.0 - PlaceInf_Scaled * P.infectiousness[InfectiousDay]; ProbSurvive *= ((ProbSurviveToday < 0) ? 0 : ProbSurviveToday); } - PlaceInf_Scaled *= (((Hosts[Person].infectiousness < 0) && (!Hosts[Person].care_home_resident)) ? (P.SymptPlaceTypeContactRate[PlaceType] * (1 - P.SymptPlaceTypeWithdrawalProp[PlaceType])) : 1); - for (; InfectiousDay < (int)Hosts[Person].recovery_or_death_time; InfectiousDay++) + PlaceInf_Scaled *= (((Hosts[person].infectiousness < 0) && (!Hosts[person].care_home_resident)) ? (P.SymptPlaceTypeContactRate[PlaceType] * (1 - P.SymptPlaceTypeWithdrawalProp[PlaceType])) : 1); + for (; InfectiousDay < (int)Hosts[person].recovery_or_death_time; InfectiousDay++) { double ProbSurviveToday = 1.0 - PlaceInf_Scaled * P.infectiousness[InfectiousDay]; ProbSurvive *= ((ProbSurviveToday < 0) ? 0 : ProbSurviveToday); @@ -818,11 +820,11 @@ void InitTransmissionCoeffs(void) double recovery_time_days = 0; double recovery_time_timesteps = 0; #pragma omp parallel for schedule(static,500) reduction(+:recovery_time_days,recovery_time_timesteps) default(none) shared(P, Hosts) - for (int Person = 0; Person < P.PopSize; Person++) + for (int person = 0; person < P.PopSize; person++) { - recovery_time_days += Hosts[Person].recovery_or_death_time * P.ModelTimeStep; - recovery_time_timesteps += Hosts[Person].recovery_or_death_time; - Hosts[Person].recovery_or_death_time = 0; // reset everybody's recovery_or_death_time + recovery_time_days += Hosts[person].recovery_or_death_time * P.ModelTimeStep; + recovery_time_timesteps += Hosts[person].recovery_or_death_time; + Hosts[person].recovery_or_death_time = 0; // reset everybody's recovery_or_death_time } // Divide total number of place infections by PopSize to get "place" R0. From 495640f9e717e6b5f7bba30348fc76f724a19e0b Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Fri, 1 Jul 2022 16:17:44 +0100 Subject: [PATCH 07/20] moved InitTransmisionCoeffs back to where it was --- src/CovidSim.cpp | 17 ++++++++--------- src/SetupModel.cpp | 4 ++-- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index d74c5caa8..f8418a634 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -174,7 +174,7 @@ int main(int argc, char* argv[]) parse_read_file(input.substr(sep + 1), snapshot_save_file); }; - double cl = (double) clock(); + double StartTime = (double) clock(); // Default bitmap format is platform dependent. #if defined(IMAGE_MAGICK) || defined(_WIN32) @@ -309,6 +309,11 @@ int main(int argc, char* argv[]) ///// initialize model (for all realisations). SetupModel(density_file, out_density_file, load_network_file, save_network_file, school_file, reg_demog_file, output_file_base); + InitTransmissionCoeffs(); + for (int i = 0; i < MAX_ADUNITS; i++) AdUnits[i].NI = 0; + for (auto const& int_file : InterventionFiles) + ReadInterventions(int_file); + Files::xfprintf_stderr("Model setup in %lf seconds\n", ((double)clock() - StartTime) / CLOCKS_PER_SEC); // Allocate memory for Efficacies array P.NumInfectionSettings = MAX_NUM_PLACE_TYPES + 2; // Maximum number of place types, plus household, and spatial @@ -319,12 +324,6 @@ int main(int argc, char* argv[]) // Allocate memory for Betas array (dynamic allocation must be done after P.NumAdunits set in SetupModel.cpp::SetupPopulation) AllocateMemForBetasArray(); - InitTransmissionCoeffs(); - for (int i = 0; i < MAX_ADUNITS; i++) AdUnits[i].NI = 0; - for (auto const& int_file : InterventionFiles) - ReadInterventions(int_file); - Files::xfprintf_stderr("Model setup in %lf seconds\n", ((double) clock() - cl) / CLOCKS_PER_SEC); - //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** @@ -365,7 +364,7 @@ int main(int argc, char* argv[]) if (P.NumRealisations > 1) { output_file = output_file_base + "." + std::to_string(Realisation); - Files::xfprintf_stderr("Realisation %i of %i (time=%lf nr_ne=%i)\n", Realisation + 1, P.NumRealisations, ((double)(clock() - cl)) / CLOCKS_PER_SEC, P.NRactNE); + Files::xfprintf_stderr("Realisation %i of %i (time=%lf nr_ne=%i)\n", Realisation + 1, P.NumRealisations, ((double)(clock() - StartTime)) / CLOCKS_PER_SEC, P.NRactNE); } ///// Set and save seeds if (((Realisation == 0) && (P.FitIter == 1)) || (P.ResetSeeds && P.KeepSameSeeds)) @@ -445,7 +444,7 @@ int main(int argc, char* argv[]) Bitmap_Finalise(); Files::xfprintf_stderr("Extinction in %i out of %i runs\n", P.NRactE, P.NRactNE + P.NRactE); - Files::xfprintf_stderr("Model ran in %lf seconds\n", ((double)clock() - cl) / CLOCKS_PER_SEC); + Files::xfprintf_stderr("Model ran in %lf seconds\n", ((double)clock() - StartTime) / CLOCKS_PER_SEC); Files::xfprintf_stderr("Model finished\n"); } } diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp index 9e0e2aff8..5e08cbb7f 100644 --- a/src/SetupModel.cpp +++ b/src/SetupModel.cpp @@ -668,7 +668,7 @@ void InitTransmissionCoeffs(void) { for (int person = Thread; person < P.PopSize; person += P.NumThreads) // loop over people { - //int AdUnit = Mcells[person->mcell].adunit; + int AdUnit = Mcells[Hosts[person].mcell].adunit; // assign susceptibility of each host. if (P.SusceptibilitySD == 0) @@ -738,7 +738,7 @@ void InitTransmissionCoeffs(void) // sum over "infectee" age groups to scale infectiousness of this infector. for (int InfecteeAge = 0; InfecteeAge < NUM_AGE_GROUPS; InfecteeAge++) //AvContactRate_Infector += P.PropAgeGroup[0][InfecteeAge] * P.WAIFW_Matrix_SpatialOnly[InfecteeAge][HOST_AGE_GROUP(person)]; // use index 0 for admin unit if doing whole country. - AvContactRate_Infector += P.PropAgeGroup[Mcells[Hosts[person].mcell].adunit][InfecteeAge] * P.WAIFW_Matrix_SpatialOnly[InfecteeAge][HOST_AGE_GROUP(person)]; + AvContactRate_Infector += P.PropAgeGroup[AdUnit][InfecteeAge] * P.WAIFW_Matrix_SpatialOnly[InfecteeAge][HOST_AGE_GROUP(person)]; // scale spatial infectiousness by weighted average. Spatial_Infectiousness *= AvContactRate_Infector; } From 0a526db5760f4c9f9679101f4e171d97448af9bf Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Fri, 1 Jul 2022 16:29:42 +0100 Subject: [PATCH 08/20] changed nim to NumImmune in InitModel --- src/CovidSim.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index f8418a634..dd672b1d6 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -896,7 +896,7 @@ void UpdateEfficacyArray() void InitModel(int run) // passing run number so we can save run number in the infection event log: ggilani - 15/10/2014 { - int nim; + int NumImmune = 0; if (P.OutputBitmap) { @@ -1020,7 +1020,6 @@ void InitModel(int run) // passing run number so we can save run number in the i for (int i = 0; i < MAX_CONTACTS+1; i++) StateT[j].contact_dist[i] = 0; } - nim = 0; if (P.DoAdUnits && P.OutputAdUnitAge) for (int Thread = 0; Thread < P.NumThreads; Thread++) @@ -1068,7 +1067,7 @@ void InitModel(int run) // passing run number so we can save run number in the i } } -#pragma omp parallel for reduction(+:nim) schedule(static,1) default(none) \ +#pragma omp parallel for reduction(+:NumImmune) schedule(static,1) default(none) \ shared(P, Cells, Hosts, Households) for (int tn = 0; tn < P.NumThreads; tn++) { @@ -1102,7 +1101,7 @@ void InitModel(int run) // passing run number so we can save run number in the i { if ((P.InitialImmunity[0] == 1) || (ranf_mt(tn) < P.InitialImmunity[0])) { - nim += Households[Hosts[k].hh].nh; + NumImmune += Households[Hosts[k].hh].nh; for (int m = Households[Hosts[k].hh].nh - 1; m >= 0; m--) DoImmune(k + m); } @@ -1114,7 +1113,7 @@ void InitModel(int run) // passing run number so we can save run number in the i int m = HOST_AGE_GROUP(k); if ((P.InitialImmunity[m] == 1) || ((P.InitialImmunity[m] > 0) && (ranf_mt(tn) < P.InitialImmunity[m]))) { - DoImmune(k); nim += 1; + DoImmune(k); NumImmune += 1; } } } From e7463d34ea929b207475228aecfc48e804c6ba5d Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Fri, 1 Jul 2022 16:34:47 +0100 Subject: [PATCH 09/20] whitespace change to rerun GitHub regression tests --- src/CovidSim.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index dd672b1d6..e3378c893 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -350,7 +350,7 @@ int main(int argc, char* argv[]) { Params::ReadParams(param_file, pre_param_file, ad_unit_file, &P, AdUnits); if (!P.FixLocalBeta) InitTransmissionCoeffs(); - output_file_base = output_file_base_f + ".f" + std::to_string(P.FitIter); + output_file_base = output_file_base_f + ".f" + std::to_string(P.FitIter); } } else StopFit = 1; From e9fe1e16dcbc1c6a191b5e5d4ac909b1bd8780ac Mon Sep 17 00:00:00 2001 From: Daniel Laydon Date: Fri, 1 Jul 2022 16:41:57 +0100 Subject: [PATCH 10/20] added pedantic sharing attribute of 'Day' in InfectSweep --- src/Sweep.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Sweep.cpp b/src/Sweep.cpp index 6ed6a4098..29da816ab 100644 --- a/src/Sweep.cpp +++ b/src/Sweep.cpp @@ -302,7 +302,7 @@ void InfectSweep(double t, int run) // added run number as argument in order to FILE* stderr_shared = stderr; #pragma omp parallel for private(CellQueue) schedule(static,1) default(none) \ - shared(t, P, CellLookup, Hosts, AdUnits, Households, Places, SamplingQueue, Cells, Mcells, StateT, seasonality, TimeStepNow, fp, BlanketMoveRestrInPlace, stderr_shared) + shared(t, P, CellLookup, Hosts, AdUnits, Households, Places, SamplingQueue, Cells, Mcells, StateT, seasonality, TimeStepNow, fp, BlanketMoveRestrInPlace, stderr_shared, Day) for (int ThreadNum = 0; ThreadNum < P.NumThreads; ThreadNum++) { Cell* ThisCell; From c8191fd9bf70b379e82f800061a93a966a2e506a Mon Sep 17 00:00:00 2001 From: dlaydon Date: Tue, 24 Jan 2023 17:31:54 +0000 Subject: [PATCH 11/20] added (dummy) mobility data to muliply betas - fitting code not done yet --- covid-sim.vcxproj | 4 +- data/param_files/p_PC7_CI_HQ_SD.txt | 632 ++++++++++++++-------------- src/Constants.h | 8 +- src/CovidSim.cpp | 214 ++++++---- src/Param.h | 2 +- src/Sweep.cpp | 2 - 6 files changed, 453 insertions(+), 409 deletions(-) diff --git a/covid-sim.vcxproj b/covid-sim.vcxproj index a9c0fb53a..a2e8014e6 100644 --- a/covid-sim.vcxproj +++ b/covid-sim.vcxproj @@ -21,13 +21,13 @@ Application true - v142 + v143 MultiByte Application false - v142 + v143 true MultiByte diff --git a/data/param_files/p_PC7_CI_HQ_SD.txt b/data/param_files/p_PC7_CI_HQ_SD.txt index d2129ad40..b5c14b686 100644 --- a/data/param_files/p_PC7_CI_HQ_SD.txt +++ b/data/param_files/p_PC7_CI_HQ_SD.txt @@ -1,317 +1,317 @@ -[Include intervention delays by admin unit] -0 - -[Vary efficacies over time] -0 - -======== PLACE CLOSURE PARAMETERS - -[Place closure start time] -7 - -[Place closure second start time] -100000 - -[Delay to place closure by admin unit] -1 1 1 - -[Duration of place closure by admin unit] -364 364 364 - -[Place closure in administrative units rather than rings] -0 - -[Administrative unit divisor for place closure] -1 - -[Place types to close for admin unit closure (0/1 array)] -1 1 1 0 - -[Cumulative proportion of place members needing to become sick for admin unit closure] -1 - -[Proportion of places in admin unit needing to pass threshold for place closure] -1 - -[Delay to start place closure] -1 - -[Duration of place closure] -364 - -[Proportion of places remaining open after closure by place type] -0 0 0.25 1 - -[Relative household contact rate after closure] -1.5 - -[Relative spatial contact rate after closure] -1.25 - -[Minimum radius for place closure] -1 - -[Place closure incidence threshold] -0 - ^^ needs to be 0 for global triggers - -[Place closure fractional incidence threshold] -0 -^^ needs to be 0 for global triggers or if abs incidence threshold used - -[Trigger incidence per cell for place closure] -0 -*** ^^^ change this for global too *** - -[Number of change times for levels of place closure] -4 - -//// Note: numbers here must match "Number of change times for levels of place closure"; that any times listed here that are before "Place closure start time" and after "Duration of place closure" are irrelevant. -[Change times for levels of place closure] -0 50 100 150 - -//// Example below gives schools closing, then opening etc. -[Proportion of places remaining open after closure by place type over time] -0 0 0.25 1 -1 1 0.25 1 -0 0 0.25 1 -1 1 0.25 1 - -[Relative household contact rates over time after place closure] -1.5 1.5 1.5 1.5 - -[Relative spatial contact rates over time after place closure] -1.25 1.25 1.25 1.25 - -[Place closure incidence threshold over time] -0 0 0 0 -^^ needs to be 0 for global triggers - -[Place closure fractional incidence threshold over time] -0 0 0 0 -^^ needs to be 0 for global triggers or if abs incidence threshold used - -[Trigger incidence per cell for place closure over time] -0 0 0 0 -*** ^^^ change these for global too *** - -//// Note: closure durations longer than interval between change times will be truncated -[Duration of place closure over time] -50 50 50 50 - - -================================== HOUSEHOLD QUARANTINE PARAMETERS - -[Household quarantine start time] -6 - -[Delay to start household quarantine] -1 - -[Delay to household quarantine by admin unit] -1 1 1 - -[Duration of household quarantine by admin unit] -364 364 364 - -[Household quarantine trigger incidence per cell] -0 - -[Length of time households are quarantined] -14 - -[Duration of household quarantine policy] -364 - -[Relative household contact rate after quarantine] -1.5 - -[Residual place contacts after household quarantine by place type] -0.25 0.25 0.25 0.25 - -[Residual spatial contacts after household quarantine] -0.25 - -[Household level compliance with quarantine] -0.75 - -[Individual level compliance with quarantine] -1.0 - -[Number of change times for levels of household quarantine] -3 - -//// Note: numbers here must match "Number of change times for levels of household quarantine"; that any times listed here that are before "Household quarantine start time" and after "Duration of household quarantine policy" are irrelevant. -[Change times for levels of household quarantine] -0 31 121 - -[Relative household contact rates over time after quarantine] -1.5 1.5 1.5 - -[Residual place contacts over time after household quarantine by place type] -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 - -[Residual spatial contacts over time after household quarantine] -0.25 0.25 0.25 - -[Household level compliance with quarantine over time] -0.75 0.75 0.75 - -[Individual level compliance with quarantine over time] -1.0 1.0 1.0 - -[Household quarantine trigger incidence per cell over time] -0 0 0 - - -=================================== CASE ISOLATION PARAMETERS - -[Case isolation start time] -6 - -[Delay to case isolation by admin unit] -1 1 1 - -[Duration of case isolation by admin unit] -364 364 364 - -[Case isolation trigger incidence per cell] -0 - -[Proportion of detected cases isolated] -0.9 - -[Delay to start case isolation] -1 - -///// i.e. for how many days will given case self-isolate? Different from "Duration of case isolation policy" -[Duration of case isolation] -7 - -[Duration of case isolation policy] -364 - -[Residual contacts after case isolation] -0.25 - -[Residual household contacts after case isolation] -0.5 - -[Number of change times for levels of case isolation] -5 - -//// Note: numbers here must match "Number of change times for levels of case isolation"; that any times listed here that are before "Case isolation start time" and after "Duration of case isolation policy" are irrelevant. -[Change times for levels of case isolation] -0 31 61 91 121 - -[Residual contacts after case isolation over time] -0.25 0.25 0.25 0.25 0.25 - -[Residual household contacts after case isolation over time] -0.5 0.5 0.5 0.5 0.5 - -[Proportion of detected cases isolated over time] -0.9 0.9 0.9 0.9 0.9 - -[Case isolation trigger incidence per cell over time] -0 0 0 0 0 - -=================================== SOCIAL DISTANCING PARAMETERS - -[Social distancing start time] -6 - -[Delay to social distancing by admin unit] -0 0 0 - -[Duration of social distancing] -364 - -[Duration of social distancing by admin unit] -364 364 364 - -[Trigger incidence per cell for social distancing] -0 - -[Relative place contact rate given social distancing by place type] -1 1 0.5 0.5 - -[Relative household contact rate given social distancing] -1.25 - -[Relative spatial contact rate given social distancing] -0.1 - -[Minimum radius for social distancing] -1 - -[Proportion compliant with enhanced social distancing] -0.0 - -[Proportion compliant with enhanced social distancing by age group] +[Include intervention delays by admin unit] +0 + +[Vary efficacies over time] +0 + +======== PLACE CLOSURE PARAMETERS + +[Place closure start time] +7 + +[Place closure second start time] +100000 + +[Delay to place closure by admin unit] +1 1 1 + +[Duration of place closure by admin unit] +364 364 364 + +[Place closure in administrative units rather than rings] +0 + +[Administrative unit divisor for place closure] +1 + +[Place types to close for admin unit closure (0/1 array)] +1 1 1 0 + +[Cumulative proportion of place members needing to become sick for admin unit closure] +1 + +[Proportion of places in admin unit needing to pass threshold for place closure] +1 + +[Delay to start place closure] +1 + +[Duration of place closure] +364 + +[Proportion of places remaining open after closure by place type] +0 0 0.25 1 + +[Relative household contact rate after closure] +1.5 + +[Relative spatial contact rate after closure] +1.25 + +[Minimum radius for place closure] +1 + +[Place closure incidence threshold] +0 + ^^ needs to be 0 for global triggers + +[Place closure fractional incidence threshold] +0 +^^ needs to be 0 for global triggers or if abs incidence threshold used + +[Trigger incidence per cell for place closure] +0 +*** ^^^ change this for global too *** + +[Number of change times for levels of place closure] +4 + +"//// Note: numbers here must match ""Number of change times for levels of place closure""; that any times listed here that are before ""Place closure start time"" and after ""Duration of place closure"" are irrelevant." +[Change times for levels of place closure] +0 50 100 150 + +"//// Example below gives schools closing, then opening etc." +[Proportion of places remaining open after closure by place type over time] +0 0 0.25 1 +1 1 0.25 1 +0 0 0.25 1 +1 1 0.25 1 + +[Relative household contact rates over time after place closure] +1.5 1.5 1.5 1.5 + +[Relative spatial contact rates over time after place closure] +1.25 1.25 1.25 1.25 + +[Place closure incidence threshold over time] +0 0 0 0 +^^ needs to be 0 for global triggers + +[Place closure fractional incidence threshold over time] +0 0 0 0 +^^ needs to be 0 for global triggers or if abs incidence threshold used + +[Trigger incidence per cell for place closure over time] +0 0 0 0 +*** ^^^ change these for global too *** + +//// Note: closure durations longer than interval between change times will be truncated +[Duration of place closure over time] +50 50 50 50 + + +================================== HOUSEHOLD QUARANTINE PARAMETERS + +[Household quarantine start time] +6 + +[Delay to start household quarantine] +1 + +[Delay to household quarantine by admin unit] +1 1 1 + +[Duration of household quarantine by admin unit] +364 364 364 + +[Household quarantine trigger incidence per cell] +0 + +[Length of time households are quarantined] +14 + +[Duration of household quarantine policy] +364 + +[Relative household contact rate after quarantine] +1.5 + +[Residual place contacts after household quarantine by place type] +0.25 0.25 0.25 0.25 + +[Residual spatial contacts after household quarantine] +0.25 + +[Household level compliance with quarantine] +0.75 + +[Individual level compliance with quarantine] +1 + +[Number of change times for levels of household quarantine] +3 + +"//// Note: numbers here must match ""Number of change times for levels of household quarantine""; that any times listed here that are before ""Household quarantine start time"" and after ""Duration of household quarantine policy"" are irrelevant." +[Change times for levels of household quarantine] +0 31 121 + +[Relative household contact rates over time after quarantine] +1.5 1.5 1.5 + +[Residual place contacts over time after household quarantine by place type] +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 + +[Residual spatial contacts over time after household quarantine] +0.25 0.25 0.25 + +[Household level compliance with quarantine over time] +0.75 0.75 0.75 + +[Individual level compliance with quarantine over time] +1.0 1.0 1.0 + +[Household quarantine trigger incidence per cell over time] +0 0 0 + + +=================================== CASE ISOLATION PARAMETERS + +[Case isolation start time] +6 + +[Delay to case isolation by admin unit] +1 1 1 + +[Duration of case isolation by admin unit] +364 364 364 + +[Case isolation trigger incidence per cell] +0 + +[Proportion of detected cases isolated] +0.9 + +[Delay to start case isolation] +1 + +"///// i.e. for how many days will given case self-isolate? Different from ""Duration of case isolation policy""" +[Duration of case isolation] +7 + +[Duration of case isolation policy] +364 + +[Residual contacts after case isolation] +0.25 + +[Residual household contacts after case isolation] +0.5 + +[Number of change times for levels of case isolation] +5 + +"//// Note: numbers here must match ""Number of change times for levels of case isolation""; that any times listed here that are before ""Case isolation start time"" and after ""Duration of case isolation policy"" are irrelevant." +[Change times for levels of case isolation] +0 31 61 91 121 + +[Residual contacts after case isolation over time] +0.25 0.25 0.25 0.25 0.25 + +[Residual household contacts after case isolation over time] +0.5 0.5 0.5 0.5 0.5 + +[Proportion of detected cases isolated over time] +0.9 0.9 0.9 0.9 0.9 + +[Case isolation trigger incidence per cell over time] +0 0 0 0 0 + +=================================== SOCIAL DISTANCING PARAMETERS + +[Social distancing start time] +6 + +[Delay to social distancing by admin unit] +0 0 0 + +[Duration of social distancing] +364 + +[Duration of social distancing by admin unit] +364 364 364 + +[Trigger incidence per cell for social distancing] +0 + +[Relative place contact rate given social distancing by place type] +1 1 0.5 0.5 + +[Relative household contact rate given social distancing] +1.25 + +[Relative spatial contact rate given social distancing] +0.1 + +[Minimum radius for social distancing] +1 + +[Proportion compliant with enhanced social distancing] +0 + +[Proportion compliant with enhanced social distancing by age group] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - -[Relative place contact rate given enhanced social distancing by place type] -0.25 0.25 0.25 0.25 - -[Relative household contact rate given enhanced social distancing] -1 - -[Relative spatial contact rate given enhanced social distancing] -0.25 - -[Delay for change in effectiveness of social distancing] -1000 - -[Relative place contact rate given social distancing by place type after change] -1 1 0.75 0.75 - -[Relative household contact rate given social distancing after change] -1.25 - -[Relative spatial contact rate given social distancing after change] -0.25 - -//// Must match "Change times for levels of social distancing" -[Number of change times for levels of social distancing] -6 - -//// Note: numbers here must match "Number of change times for levels of social distancing"; that any times listed here that are before "Social distancing start time" and after "Duration of social distancing" are irrelevant. -[Change times for levels of social distancing] -0 31 61 91 121 151 - -//// Again, want this to supercede "Relative place contact rate given social distancing by place type". Should be matrix of dimension "Number of change times for levels of social distancing" by Number of place types. -[Relative place contact rates over time given social distancing by place type] -1 1 0.5 0.5 -1 1 0.5 0.5 -1 1 0.5 0.5 -1 1 0.5 0.5 -1 1 0.5 0.5 -1 1 0.5 0.5 - -//// Ideally want this to supercede "Relative household contact rate given social distancing", but need to preserved backwards compatibility for now. -[Relative household contact rates over time given social distancing] -1.25 1.25 1.25 1.25 1.25 1.25 - -//// Ideally want this to supercede "Relative spatial contact rate given social distancing", but need to preserved backwards compatibility for now. -[Relative spatial contact rates over time given social distancing] -0.1 0.1 0.1 0.1 0.1 0.1 - -[Relative household contact rates over time given enhanced social distancing] -1 1 1 1 1 1 - -[Relative spatial contact rates over time given enhanced social distancing] -0.25 0.25 0.25 0.25 0.25 0.25 - -[Relative place contact rates over time given enhanced social distancing by place type] -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 - -[Trigger incidence per cell for social distancing over time] -100 100 100 100 100 100 + +[Relative place contact rate given enhanced social distancing by place type] +0.25 0.25 0.25 0.25 + +[Relative household contact rate given enhanced social distancing] +1 + +[Relative spatial contact rate given enhanced social distancing] +0.25 + +[Delay for change in effectiveness of social distancing] +1000 + +[Relative place contact rate given social distancing by place type after change] +1 1 0.75 0.75 + +[Relative household contact rate given social distancing after change] +1.25 + +[Relative spatial contact rate given social distancing after change] +0.25 + +"//// Must match ""Change times for levels of social distancing""" +[Number of change times for levels of social distancing] +6 + +"//// Note: numbers here must match ""Number of change times for levels of social distancing""; that any times listed here that are before ""Social distancing start time"" and after ""Duration of social distancing"" are irrelevant." +[Change times for levels of social distancing] +0 31 61 91 121 151 + +"//// Again, want this to supercede ""Relative place contact rate given social distancing by place type"". Should be matrix of dimension ""Number of change times for levels of social distancing"" by Number of place types." +[Relative place contact rates over time given social distancing by place type] +1 1 0.5 0.5 +1 1 0.5 0.5 +1 1 0.5 0.5 +1 1 0.5 0.5 +1 1 0.5 0.5 +1 1 0.5 0.5 + +"//// Ideally want this to supercede ""Relative household contact rate given social distancing"", but need to preserved backwards compatibility for now." +[Relative household contact rates over time given social distancing] +1.25 1.25 1.25 1.25 1.25 1.25 + +"//// Ideally want this to supercede ""Relative spatial contact rate given social distancing"", but need to preserved backwards compatibility for now." +[Relative spatial contact rates over time given social distancing] +0.1 0.1 0.1 0.1 0.1 0.1 + +[Relative household contact rates over time given enhanced social distancing] +1 1 1 1 1 1 + +[Relative spatial contact rates over time given enhanced social distancing] +0.25 0.25 0.25 0.25 0.25 0.25 + +[Relative place contact rates over time given enhanced social distancing by place type] +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 + +[Trigger incidence per cell for social distancing over time] +100 100 100 100 100 100 diff --git a/src/Constants.h b/src/Constants.h index 5c887f7f3..7ff415a87 100644 --- a/src/Constants.h +++ b/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(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 diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index e3378c893..71fa9de92 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -67,32 +67,57 @@ void CalibrationThresholdCheck(double, int); void CalcLikelihood(int, std::string const&, std::string const&); void CalcOriginDestMatrix_adunit(void); //added function to calculate origin destination matrix: ggilani 28/01/15 -void AllocateMemForBetasArray() // called in main (once per fitting run) +bool FileExists(std::string FileName) +{ + std::ifstream infile(FileName); + bool FileOkay = infile.good(); + infile.close(); + return FileOkay; +} +void AllocateMemForBetasArray() // called in main (once per fitting run). Can you allocate only once but clear values? { P.Betas = new double** [(int)P.SimulationDuration + 1](); for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) { P.Betas[Day] = new double* [P.NumAdunits](); - for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) P.Betas[Day][AdUnit] = new double[P.NumInfectionSettings](); + for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) + P.Betas[Day][AdUnit] = new double[P.NumInfectionSettings](); } + // initialize all values to 1; + for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) + for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) + for (int Setting = 0; Setting < P.NumInfectionSettings; Setting++) + P.Betas[Day][AdUnit][Setting] = 1; } -void InitBetasArray() // called in InitModel (every realistaion/parameter guess iteration). +void ImportMobilityDataToBetasArray(std::string MobilityFilename, int SettingNumber) // assumes all P.Betas values initialized to 1. { - //if (P.VaryBetasOverTimeByRegion == false) // while you're setting up, leave this condition out. Put back in later once mobility data included. - //{ - //// if not varying by region or time, assign them to be the same for every simulation day and every admin unit. - for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) + if (MobilityFilename != "NoFileFound" && FileExists(MobilityFilename)) + { + std::ifstream MobFileStream(MobilityFilename.c_str()); + double ProportionalChange = 0.0; + for (int Day = 0; Day < P.SimulationDuration; Day++) for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) { - // place (by type) - for (int PlaceType = 0; PlaceType < P.NumPlaceTypes; PlaceType++) - P.Betas[Day][AdUnit][PlaceType] = P.PlaceTypeTrans[PlaceType]; - // Household - P.Betas[Day][AdUnit][House] = P.HouseholdTrans; - // Spatial/Community - P.Betas[Day][AdUnit][Spatial] = P.LocalBeta; + MobFileStream >> ProportionalChange; + P.Betas[Day][AdUnit][SettingNumber] += ProportionalChange; } - //} + MobFileStream.close(); + } + else ERR_CRITICAL(("mobility_file: " + MobilityFilename + " does not exist").c_str()); +} +void MultiplyMobDataWithTransmissionParams() +{ + for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) + for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) + { + // place (by type) + for (int PlaceType = 0; PlaceType < P.NumPlaceTypes; PlaceType++) + P.Betas[Day][AdUnit][PlaceType] *= P.PlaceTypeTrans[PlaceType]; + // Household + P.Betas[Day][AdUnit][House] *= P.HouseholdTrans; + // Spatial/Community + P.Betas[Day][AdUnit][Spatial] *= P.LocalBeta; + } } ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// @@ -153,6 +178,8 @@ int main(int argc, char* argv[]) std::string reg_demog_file, fit_file, data_file; std::string ad_unit_file, out_density_file, output_file_base; std::string snapshot_load_file, snapshot_save_file; + // mobility data filenames + std::string spatial_mob_file = "NoFileFound", household_mob_file = "NoFileFound", workplace_mob_file = "NoFileFound"; // set default to "NoFileFound" int StopFit = 0; ///// Flags to ensure various parameters have been read; set to false as default. @@ -196,7 +223,7 @@ int main(int argc, char* argv[]) //// **** PARSE COMMAND-LINE ARGS //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** - + CmdLineArgs args; args.add_string_option("A", parse_read_file, ad_unit_file, "Administrative Division"); args.add_string_option("AP", parse_read_file, air_travel_file, "Air travel data file"); @@ -204,6 +231,11 @@ int main(int argc, char* argv[]) args.add_integer_option("c", P.MaxNumThreads, "Number of threads to use"); args.add_integer_option("C", P.PlaceCloseIndepThresh, "Sets the P.PlaceCloseIndepThresh parameter"); + // Mobility data file names (CL args) + args.add_string_option("MS", parse_read_file, spatial_mob_file , "Spatial mobility data file"); + args.add_string_option("MH", parse_read_file, household_mob_file, "Household mobility data file"); + args.add_string_option("MW", parse_read_file, workplace_mob_file, "Workplace mobility data file"); + /* Wes: Need to allow /CLPxx up to 99. I'll do this naively for now and prevent the help text from looking overly verybose. To satisfiy all behaviour, /CLP0: to /CLP9: should do the same as /CLP00: to /CLP09: - both valid, as are /CLP10: to /CLP99: @@ -321,10 +353,19 @@ int main(int argc, char* argv[]) P.Efficacies = new double* [P.NumInterventionClasses](); for (int intervention = 0; intervention < P.NumInterventionClasses; intervention++) P.Efficacies[intervention] = new double[P.NumInfectionSettings](); - // Allocate memory for Betas array (dynamic allocation must be done after P.NumAdunits set in SetupModel.cpp::SetupPopulation) - AllocateMemForBetasArray(); - - + // Allocate memory for Betas array (dynamic allocation - must be done after P.NumAdunits set in SetupModel.cpp::SetupPopulation) + // By default, make all multipliers of betas equal to 1. If P.VaryBetasOverTimeByRegion, these multipliers will be mobility data + AllocateMemForBetasArray(); + if (P.VaryBetasOverTimeByRegion) + { + ImportMobilityDataToBetasArray(workplace_mob_file , Workplace); + ImportMobilityDataToBetasArray(household_mob_file , House); + ImportMobilityDataToBetasArray(spatial_mob_file , Spatial); + } + //for (int Setting = 0; Setting < P.NumInfectionSettings; Setting++) + // for (int Day = 0; Day < (int)P.SimulationDuration; Day++) + // for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) + // std::cout << "Day " << Day << ", AdUnit " << AdUnit << ", Setting " << Setting << ", Beta " << P.Betas[Day][AdUnit][Setting] << std::endl; ; //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** RUN MODEL @@ -378,6 +419,7 @@ int main(int argc, char* argv[]) int ContCalib, ModelCalibLoop = 0; P.StopCalibration = P.ModelCalibIteration = ModelCalibLoop = 0; + //// Run model (possibly repeatedly if calibrating). do { // has been interrupted to reset holiday time. Note that this currently only happens in the first run, regardless of how many realisations are being run. if ((P.ModelCalibIteration % 14 == 0) && (ModelCalibLoop < 4)) @@ -1209,7 +1251,7 @@ void InitModel(int run) // passing run number so we can save run number in the i //// Add all of the above to P.Efficacies array. UpdateEfficacyArray(); //// InitializeBetasArray (either to be same for all days and regions), or otherwise depending on how we are modelling them). Memory allocated in main using - InitBetasArray(); + MultiplyMobDataWithTransmissionParams(); // Initialize CFR scalings P.CFR_Critical_Scale_Current = P.CFR_TimeScaling_Critical[0]; @@ -3718,73 +3760,75 @@ void CalibrationThresholdCheck(double t,int n) } } -void CalcLikelihood(int run, std::string const& DataFile, std::string const& OutFileBase) +void CalcLikelihood(int RealisationNum, std::string const& DataFile, std::string const& OutFileBase) { - FILE* dat; + FILE* DataFile_dat; - static int DataAlreadyRead = 0, ncols, nrows, * ColTypes; // static i.e. won't be reset upon subsequent calls to CalcLikelihood - static double** Data, NegBinK, sumL; + // static i.e. won't be destroyed at end of function. So in subsequent calls to this function, Data won't have to be read in again, and can add to sum_LogLik. + static int DataAlreadyRead = 0, ncols, nrows, * ColTypes; + static double** Data, NegBinK, sum_LogLik; + // Read in data (if not already read during previous call to this function). if (!DataAlreadyRead) { char FieldName[1024]; - dat = Files::xfopen(DataFile.c_str(), "r"); - // Extract numbers of rows and columns, and overdispersion parameter of negative bionomial distribution, from Data file - Files::xfscanf(dat, 3, "%i %i %lg", &nrows, &ncols, &NegBinK); + DataFile_dat = Files::xfopen(DataFile.c_str(), "r"); + + // From first row of DataFile, extract numbers of rows and columns, and overdispersion parameter of negative bionomial distribution. + Files::xfscanf(DataFile_dat, 3, "%i %i %lg", &nrows, &ncols, &NegBinK); // allocate memory - ColTypes = (int*) Memory::xcalloc(ncols, sizeof(int)); - Data = (double**) Memory::xcalloc(nrows, sizeof(double *)); - for (int i = 0; i < nrows; i++) - Data[i] = (double*) Memory::xcalloc(ncols, sizeof(double)); + ColTypes = (int*) Memory::xcalloc(ncols, sizeof(int)); + Data = (double**) Memory::xcalloc(nrows, sizeof(double *)); + for (int row = 0; row < nrows; row++) + Data[row] = (double*) Memory::xcalloc(ncols, sizeof(double)); - // cycle through columns assigning an int label to each data/column type in data file. Essentially renaming column names to integers. - for (int i = 0; i < ncols; i++) + // From next row of DataFile, cycle through columns assigning an int label to each data/column type in data file. Essentially renaming column names to integers. + for (int col = 0; col < ncols; col++) { - ColTypes[i] = -100; - Files::xfscanf(dat, 1, "%s", FieldName); + ColTypes[col] = -100; // by default assign to negative number as this won't be picked up in loop below. + + // assign entry of this row and column of DataFile to be the variable name + Files::xfscanf(DataFile_dat, 1, "%s", FieldName); + + //// cycle through possibilities if (!strcmp(FieldName, "day")) { - ColTypes[i] = -1; - if (i != 0) ERR_CRITICAL("'day' must be first column in data file\n"); + ColTypes[col] = -1; + if (col != 0) ERR_CRITICAL("'day' must be first column in data file\n"); } - if (!strcmp(FieldName, "all_deaths")) - ColTypes[i] = 0; - else if (!strcmp(FieldName, "hospital_deaths")) - ColTypes[i] = 1; - else if (!strcmp(FieldName, "care_home_deaths")) - ColTypes[i] = 2; - else if (!strcmp(FieldName, "seroprevalence_numerator")) - ColTypes[i] = 3; + if (!strcmp(FieldName, "all_deaths")) ColTypes[col] = 0; + else if (!strcmp(FieldName, "hospital_deaths")) ColTypes[col] = 1; + else if (!strcmp(FieldName, "care_home_deaths")) ColTypes[col] = 2; + else if (!strcmp(FieldName, "seroprevalence_numerator")) ColTypes[col] = 3; else if (!strcmp(FieldName, "seroprevalence_denominator")) { - ColTypes[i] = 4; - if (ColTypes[i - 1] != 3) ERR_CRITICAL("Seroprevalence denominator must be next column after numerator in data file\n"); + ColTypes[col] = 4; + if (ColTypes[col - 1] != 3) ERR_CRITICAL("Seroprevalence denominator must be next column after numerator in data file\n"); } - else if (!strcmp(FieldName, "infection_prevalence_numerator")) - ColTypes[i] = 5; + else if (!strcmp(FieldName, "infection_prevalence_numerator")) ColTypes[col] = 5; else if (!strcmp(FieldName, "infection_prevalence_denominator")) { - ColTypes[i] = 6; - if (ColTypes[i - 1] != 5) ERR_CRITICAL("Infection prevalence denominator must be next column after numerator in data file\n"); + ColTypes[col] = 6; + if (ColTypes[col - 1] != 5) ERR_CRITICAL("Infection prevalence denominator must be next column after numerator in data file\n"); } } // extract data into Data array. - for (int i = 0; i < nrows; i++) - for (int j = 0; j < ncols; j++) - Files::xfscanf(dat, 1, "%lg", &(Data[i][j])); - Files::xfclose(dat); + for (int row = 0; row < nrows; row++) + for (int col = 0; col < ncols; col++) + Files::xfscanf(DataFile_dat, 1, "%lg", &(Data[row][col])); + Files::xfclose(DataFile_dat); DataAlreadyRead = 1; } // calculate likelihood function - double c, LL = 0.0; + double c, LogLik = 0.0; double kp = (P.clP[99] > 0) ? P.clP[99] : NegBinK; // clP[99] reserved for fitting overdispersion. If not positive-definite assign to NegBinK extracted above. c = 1.0; // 1 / ((double)(P.NRactE + P.NRactNE)); int offset = (P.Interventions_StartDate_CalTime > 0) ? ((int)(P.Interventions_StartDate_CalTime - P.DateTriggerReached_SimTime)) : 0; - for (int col = 1; col < ncols; col++) /// cycle through columns (different sources of data contributing to likelihood), and add to log likelihood (LL) accordingly. + for (int col = 1; col < ncols; col++) /// cycle through columns (different sources of data contributing to likelihood), and add to log likelihood (LL) accordingly. Note starts from 1 as 0th column is day of year. { if ((ColTypes[col] >= 0) && (ColTypes[col] <= 2)) // i.e. "all deaths", "hospital deaths", "care home deaths" { @@ -3795,12 +3839,14 @@ void CalcLikelihood(int run, std::string const& DataFile, std::string const& Out if ((Data[row][col] >= -1) && (day < P.NumOutputTimeSteps)) // data is not NA (-ve) and within time range of model run { double ModelValue; - if (ColTypes[col] == 0) - ModelValue = c * TimeSeries[day - offset].incD; // all deaths by date of death - else if (ColTypes[col] == 1) - ModelValue = c * (TimeSeries[day - offset].incDeath_Critical + TimeSeries[day - offset].incDeath_SARI); // hospital deaths (SARI and Critical) by date of death - else if (ColTypes[col] == 2) - ModelValue = c * TimeSeries[day - offset].incDeath_ILI; // care home deaths (ILI) by date of death + + if (ColTypes[col] == 0) // all deaths by date of death + ModelValue = c * TimeSeries[day - offset].incD; + else if (ColTypes[col] == 1) // hospital deaths (SARI and Critical) by date of death + ModelValue = c * (TimeSeries[day - offset].incDeath_Critical + TimeSeries[day - offset].incDeath_SARI); + else if (ColTypes[col] == 2) // care home deaths (ILI) by date of death + ModelValue = c * TimeSeries[day - offset].incDeath_ILI; + ModelValueSum += ModelValue; if (Data[row][col] >= 0) { @@ -3811,13 +3857,13 @@ void CalcLikelihood(int run, std::string const& DataFile, std::string const& Out } if (NegBinK >= 10000) //prob model and data from same underlying poisson - LL += lgamma(2 * (Data[row][col] + ModelValue) + 1) - lgamma(Data[row][col] + ModelValue + 1) - lgamma(Data[row][col] + 1) - lgamma(ModelValue + 1) - (3 * (Data[row][col] + ModelValue) + 1) * log(2); + LogLik += lgamma(2 * (Data[row][col] + ModelValue) + 1) - lgamma(Data[row][col] + ModelValue + 1) - lgamma(Data[row][col] + 1) - lgamma(ModelValue + 1) - (3 * (Data[row][col] + ModelValue) + 1) * log(2); else { //neg bin LL (NegBinK=1 implies no over-dispersion. >1 implies more) double knb = 1.0 + ModelValue / kp; double pnb = kp / (1.0 + kp); - LL += lgamma(Data[row][col] + knb) - lgamma(Data[row][col] + 1) - lgamma(knb) + knb * log(1.0 - pnb) + Data[row][col] * log(pnb); + LogLik += lgamma(Data[row][col] + knb) - lgamma(Data[row][col] + 1) - lgamma(knb) + knb * log(1.0 - pnb) + Data[row][col] * log(pnb); } } } @@ -3830,9 +3876,9 @@ void CalcLikelihood(int run, std::string const& DataFile, std::string const& Out int day = (int)Data[row][0]; // day is day of year - directly indexes TimeSeries[] if ((Data[row][col] >= 0) && (day < P.NumOutputTimeSteps)) // data is not NA (-ve) and within time range of model run { - double m = Data[row][col]; // numerator - double N = Data[row][col + 1]; // denominator - double ModelValue = 0.0; + double m = Data[row][col]; // numerator + double N = Data[row][col + 1]; // denominator + double ModelValue = 0.0; // (re)-initialize ModelValue for (int k = offset; k < day; k++) // loop over all days of infection up to day of sample { double prob_seroconvert = P.SeroConvMaxSens * (1.0 - 0.5 * ((exp(-((double)(_I64(day) - k)) * P.SeroConvP1) + 1.0) * exp(-((double)(_I64(day) - k)) * P.SeroConvP2))); // add P1 to P2 to prevent degeneracy @@ -3840,7 +3886,7 @@ void CalcLikelihood(int run, std::string const& DataFile, std::string const& Out } ModelValue += c * TimeSeries[day - offset].S * (1.0 - P.SeroConvSpec); ModelValue /= ((double)P.PopSize); - LL += m * log((ModelValue + 1e-20) / (m / N + 1e-20)) + (N - m) * log((1.0 - ModelValue + 1e-20) / (1.0 - m / N + 1e-20)); // subtract saturated likelihood + LogLik += m * log((ModelValue + 1e-20) / (m / N + 1e-20)) + (N - m) * log((1.0 - ModelValue + 1e-20) / (1.0 - m / N + 1e-20)); // subtract saturated likelihood } } } @@ -3851,32 +3897,32 @@ void CalcLikelihood(int run, std::string const& DataFile, std::string const& Out int day = (int)Data[row][0]; // day is day of year - directly indexes TimeSeries[] if ((Data[row][col] >= 0) && (day < P.NumOutputTimeSteps)) // data is not NA (-ve) and within time range of model run { - double m = Data[row][col]; // numerator - double N = Data[row][col + 1]; // denominator + double m = Data[row][col]; // numerator + double N = Data[row][col + 1]; // denominator double ModelValue = P.InfPrevSurveyScale * c * TimeSeries[day - offset].I / ((double)P.PopSize); - LL += m * log(ModelValue + 1e-20) + (N - m) * log(1.0 - ModelValue); + LogLik += m * log(ModelValue + 1e-20) + (N - m) * log(1.0 - ModelValue); } } } } - Files::xfprintf_stderr("Log-likelihood = %lg\n", LL); - if (run == 0) - sumL = LL; + Files::xfprintf_stderr("Log-likelihood = %lg\n", LogLik); + if (RealisationNum == 0) sum_LogLik = LogLik; // initialize sum to be likelihood for this realisation. else { - double maxLL = LL; - if (sumL > maxLL) maxLL = sumL; - sumL = maxLL + log(exp(sumL - maxLL) + exp(LL - maxLL)); + // maxLL + double maxLL = LogLik; + if (sum_LogLik > maxLL) maxLL = sum_LogLik; + sum_LogLik = maxLL + log(exp(sum_LogLik - maxLL) + exp(LogLik - maxLL)); } - if (run + 1 == P.NumRealisations) // at final realisation, output log-likelihood + if (RealisationNum + 1 == P.NumRealisations) // at final realisation, output log-likelihood { - LL = sumL - log((double)P.NumRealisations); + double AverageLogLik_AcrossRealisations = sum_LogLik - log((double)P.NumRealisations); /// i.e. average likelihood = sum_LogLik / NumRealisations std::string TmpFile = OutFileBase + ".ll.tmp"; std::string OutFile = OutFileBase + ".ll.txt"; - dat = Files::xfopen(TmpFile.c_str(), "w"); - Files::xfprintf(dat, "%i\t%.8lg\n", P.FitIter, LL); - Files::xfclose(dat); + DataFile_dat = Files::xfopen(TmpFile.c_str(), "w"); + Files::xfprintf(DataFile_dat, "%i\t%.8lg\n", P.FitIter, AverageLogLik_AcrossRealisations); + Files::xfclose(DataFile_dat); Files::xrename(TmpFile.c_str(), OutFile.c_str()); // rename only when file is complete and closed } } diff --git a/src/Param.h b/src/Param.h index 4b02b9d12..4cf85c76f 100644 --- a/src/Param.h +++ b/src/Param.h @@ -120,7 +120,7 @@ struct Param // setting refers to either place, household and spatial/community */ double*** Betas; //// indexed by i) time; ii) region / admin unit; iii) setting. - bool VaryBetasOverTimeByRegion = false; + bool VaryBetasOverTimeByRegion = true; double R0, R0scale; diff --git a/src/Sweep.cpp b/src/Sweep.cpp index 29da816ab..802446a44 100644 --- a/src/Sweep.cpp +++ b/src/Sweep.cpp @@ -343,8 +343,6 @@ void InfectSweep(double t, int run) // added run number as argument in order to (Hosts[InfectiousPersonIndex].digitalContactTracingUser == 1)); // && (TimeStepNow <= (Hosts[InfectiousPersonIndex].detected_time + P.usCaseIsolationDelay))); // BEGIN HOUSEHOLD INFECTIONS - - Household_Beta = (P.DoHouseholds) ? P.Betas[Day][AdUnit_ThisPerson][House] * seasonality * fp : 0; if (Household_Beta > 0) { From 050203e890cba52bc8b3235e4a46f05129a33e26 Mon Sep 17 00:00:00 2001 From: dlaydon Date: Wed, 25 Jan 2023 10:11:38 +0000 Subject: [PATCH 12/20] turning off variable contact rates for betas for github workflow runs --- src/Param.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Param.h b/src/Param.h index 4cf85c76f..4b02b9d12 100644 --- a/src/Param.h +++ b/src/Param.h @@ -120,7 +120,7 @@ struct Param // setting refers to either place, household and spatial/community */ double*** Betas; //// indexed by i) time; ii) region / admin unit; iii) setting. - bool VaryBetasOverTimeByRegion = true; + bool VaryBetasOverTimeByRegion = false; double R0, R0scale; From 0be027786d4aa9db51f783e200bc61af49077ffd Mon Sep 17 00:00:00 2001 From: dlaydon Date: Wed, 25 Jan 2023 14:26:17 +0000 Subject: [PATCH 13/20] old p_PC7 etc. param file back --- data/param_files/p_PC7_CI_HQ_SD.txt | 632 ++++++++++++++-------------- src/CovidSim.cpp | 15 +- src/Param.h | 3 +- 3 files changed, 327 insertions(+), 323 deletions(-) diff --git a/data/param_files/p_PC7_CI_HQ_SD.txt b/data/param_files/p_PC7_CI_HQ_SD.txt index b5c14b686..d2129ad40 100644 --- a/data/param_files/p_PC7_CI_HQ_SD.txt +++ b/data/param_files/p_PC7_CI_HQ_SD.txt @@ -1,317 +1,317 @@ -[Include intervention delays by admin unit] -0 - -[Vary efficacies over time] -0 - -======== PLACE CLOSURE PARAMETERS - -[Place closure start time] -7 - -[Place closure second start time] -100000 - -[Delay to place closure by admin unit] -1 1 1 - -[Duration of place closure by admin unit] -364 364 364 - -[Place closure in administrative units rather than rings] -0 - -[Administrative unit divisor for place closure] -1 - -[Place types to close for admin unit closure (0/1 array)] -1 1 1 0 - -[Cumulative proportion of place members needing to become sick for admin unit closure] -1 - -[Proportion of places in admin unit needing to pass threshold for place closure] -1 - -[Delay to start place closure] -1 - -[Duration of place closure] -364 - -[Proportion of places remaining open after closure by place type] -0 0 0.25 1 - -[Relative household contact rate after closure] -1.5 - -[Relative spatial contact rate after closure] -1.25 - -[Minimum radius for place closure] -1 - -[Place closure incidence threshold] -0 - ^^ needs to be 0 for global triggers - -[Place closure fractional incidence threshold] -0 -^^ needs to be 0 for global triggers or if abs incidence threshold used - -[Trigger incidence per cell for place closure] -0 -*** ^^^ change this for global too *** - -[Number of change times for levels of place closure] -4 - -"//// Note: numbers here must match ""Number of change times for levels of place closure""; that any times listed here that are before ""Place closure start time"" and after ""Duration of place closure"" are irrelevant." -[Change times for levels of place closure] -0 50 100 150 - -"//// Example below gives schools closing, then opening etc." -[Proportion of places remaining open after closure by place type over time] -0 0 0.25 1 -1 1 0.25 1 -0 0 0.25 1 -1 1 0.25 1 - -[Relative household contact rates over time after place closure] -1.5 1.5 1.5 1.5 - -[Relative spatial contact rates over time after place closure] -1.25 1.25 1.25 1.25 - -[Place closure incidence threshold over time] -0 0 0 0 -^^ needs to be 0 for global triggers - -[Place closure fractional incidence threshold over time] -0 0 0 0 -^^ needs to be 0 for global triggers or if abs incidence threshold used - -[Trigger incidence per cell for place closure over time] -0 0 0 0 -*** ^^^ change these for global too *** - -//// Note: closure durations longer than interval between change times will be truncated -[Duration of place closure over time] -50 50 50 50 - - -================================== HOUSEHOLD QUARANTINE PARAMETERS - -[Household quarantine start time] -6 - -[Delay to start household quarantine] -1 - -[Delay to household quarantine by admin unit] -1 1 1 - -[Duration of household quarantine by admin unit] -364 364 364 - -[Household quarantine trigger incidence per cell] -0 - -[Length of time households are quarantined] -14 - -[Duration of household quarantine policy] -364 - -[Relative household contact rate after quarantine] -1.5 - -[Residual place contacts after household quarantine by place type] -0.25 0.25 0.25 0.25 - -[Residual spatial contacts after household quarantine] -0.25 - -[Household level compliance with quarantine] -0.75 - -[Individual level compliance with quarantine] -1 - -[Number of change times for levels of household quarantine] -3 - -"//// Note: numbers here must match ""Number of change times for levels of household quarantine""; that any times listed here that are before ""Household quarantine start time"" and after ""Duration of household quarantine policy"" are irrelevant." -[Change times for levels of household quarantine] -0 31 121 - -[Relative household contact rates over time after quarantine] -1.5 1.5 1.5 - -[Residual place contacts over time after household quarantine by place type] -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 - -[Residual spatial contacts over time after household quarantine] -0.25 0.25 0.25 - -[Household level compliance with quarantine over time] -0.75 0.75 0.75 - -[Individual level compliance with quarantine over time] -1.0 1.0 1.0 - -[Household quarantine trigger incidence per cell over time] -0 0 0 - - -=================================== CASE ISOLATION PARAMETERS - -[Case isolation start time] -6 - -[Delay to case isolation by admin unit] -1 1 1 - -[Duration of case isolation by admin unit] -364 364 364 - -[Case isolation trigger incidence per cell] -0 - -[Proportion of detected cases isolated] -0.9 - -[Delay to start case isolation] -1 - -"///// i.e. for how many days will given case self-isolate? Different from ""Duration of case isolation policy""" -[Duration of case isolation] -7 - -[Duration of case isolation policy] -364 - -[Residual contacts after case isolation] -0.25 - -[Residual household contacts after case isolation] -0.5 - -[Number of change times for levels of case isolation] -5 - -"//// Note: numbers here must match ""Number of change times for levels of case isolation""; that any times listed here that are before ""Case isolation start time"" and after ""Duration of case isolation policy"" are irrelevant." -[Change times for levels of case isolation] -0 31 61 91 121 - -[Residual contacts after case isolation over time] -0.25 0.25 0.25 0.25 0.25 - -[Residual household contacts after case isolation over time] -0.5 0.5 0.5 0.5 0.5 - -[Proportion of detected cases isolated over time] -0.9 0.9 0.9 0.9 0.9 - -[Case isolation trigger incidence per cell over time] -0 0 0 0 0 - -=================================== SOCIAL DISTANCING PARAMETERS - -[Social distancing start time] -6 - -[Delay to social distancing by admin unit] -0 0 0 - -[Duration of social distancing] -364 - -[Duration of social distancing by admin unit] -364 364 364 - -[Trigger incidence per cell for social distancing] -0 - -[Relative place contact rate given social distancing by place type] -1 1 0.5 0.5 - -[Relative household contact rate given social distancing] -1.25 - -[Relative spatial contact rate given social distancing] -0.1 - -[Minimum radius for social distancing] -1 - -[Proportion compliant with enhanced social distancing] -0 - -[Proportion compliant with enhanced social distancing by age group] +[Include intervention delays by admin unit] +0 + +[Vary efficacies over time] +0 + +======== PLACE CLOSURE PARAMETERS + +[Place closure start time] +7 + +[Place closure second start time] +100000 + +[Delay to place closure by admin unit] +1 1 1 + +[Duration of place closure by admin unit] +364 364 364 + +[Place closure in administrative units rather than rings] +0 + +[Administrative unit divisor for place closure] +1 + +[Place types to close for admin unit closure (0/1 array)] +1 1 1 0 + +[Cumulative proportion of place members needing to become sick for admin unit closure] +1 + +[Proportion of places in admin unit needing to pass threshold for place closure] +1 + +[Delay to start place closure] +1 + +[Duration of place closure] +364 + +[Proportion of places remaining open after closure by place type] +0 0 0.25 1 + +[Relative household contact rate after closure] +1.5 + +[Relative spatial contact rate after closure] +1.25 + +[Minimum radius for place closure] +1 + +[Place closure incidence threshold] +0 + ^^ needs to be 0 for global triggers + +[Place closure fractional incidence threshold] +0 +^^ needs to be 0 for global triggers or if abs incidence threshold used + +[Trigger incidence per cell for place closure] +0 +*** ^^^ change this for global too *** + +[Number of change times for levels of place closure] +4 + +//// Note: numbers here must match "Number of change times for levels of place closure"; that any times listed here that are before "Place closure start time" and after "Duration of place closure" are irrelevant. +[Change times for levels of place closure] +0 50 100 150 + +//// Example below gives schools closing, then opening etc. +[Proportion of places remaining open after closure by place type over time] +0 0 0.25 1 +1 1 0.25 1 +0 0 0.25 1 +1 1 0.25 1 + +[Relative household contact rates over time after place closure] +1.5 1.5 1.5 1.5 + +[Relative spatial contact rates over time after place closure] +1.25 1.25 1.25 1.25 + +[Place closure incidence threshold over time] +0 0 0 0 +^^ needs to be 0 for global triggers + +[Place closure fractional incidence threshold over time] +0 0 0 0 +^^ needs to be 0 for global triggers or if abs incidence threshold used + +[Trigger incidence per cell for place closure over time] +0 0 0 0 +*** ^^^ change these for global too *** + +//// Note: closure durations longer than interval between change times will be truncated +[Duration of place closure over time] +50 50 50 50 + + +================================== HOUSEHOLD QUARANTINE PARAMETERS + +[Household quarantine start time] +6 + +[Delay to start household quarantine] +1 + +[Delay to household quarantine by admin unit] +1 1 1 + +[Duration of household quarantine by admin unit] +364 364 364 + +[Household quarantine trigger incidence per cell] +0 + +[Length of time households are quarantined] +14 + +[Duration of household quarantine policy] +364 + +[Relative household contact rate after quarantine] +1.5 + +[Residual place contacts after household quarantine by place type] +0.25 0.25 0.25 0.25 + +[Residual spatial contacts after household quarantine] +0.25 + +[Household level compliance with quarantine] +0.75 + +[Individual level compliance with quarantine] +1.0 + +[Number of change times for levels of household quarantine] +3 + +//// Note: numbers here must match "Number of change times for levels of household quarantine"; that any times listed here that are before "Household quarantine start time" and after "Duration of household quarantine policy" are irrelevant. +[Change times for levels of household quarantine] +0 31 121 + +[Relative household contact rates over time after quarantine] +1.5 1.5 1.5 + +[Residual place contacts over time after household quarantine by place type] +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 + +[Residual spatial contacts over time after household quarantine] +0.25 0.25 0.25 + +[Household level compliance with quarantine over time] +0.75 0.75 0.75 + +[Individual level compliance with quarantine over time] +1.0 1.0 1.0 + +[Household quarantine trigger incidence per cell over time] +0 0 0 + + +=================================== CASE ISOLATION PARAMETERS + +[Case isolation start time] +6 + +[Delay to case isolation by admin unit] +1 1 1 + +[Duration of case isolation by admin unit] +364 364 364 + +[Case isolation trigger incidence per cell] +0 + +[Proportion of detected cases isolated] +0.9 + +[Delay to start case isolation] +1 + +///// i.e. for how many days will given case self-isolate? Different from "Duration of case isolation policy" +[Duration of case isolation] +7 + +[Duration of case isolation policy] +364 + +[Residual contacts after case isolation] +0.25 + +[Residual household contacts after case isolation] +0.5 + +[Number of change times for levels of case isolation] +5 + +//// Note: numbers here must match "Number of change times for levels of case isolation"; that any times listed here that are before "Case isolation start time" and after "Duration of case isolation policy" are irrelevant. +[Change times for levels of case isolation] +0 31 61 91 121 + +[Residual contacts after case isolation over time] +0.25 0.25 0.25 0.25 0.25 + +[Residual household contacts after case isolation over time] +0.5 0.5 0.5 0.5 0.5 + +[Proportion of detected cases isolated over time] +0.9 0.9 0.9 0.9 0.9 + +[Case isolation trigger incidence per cell over time] +0 0 0 0 0 + +=================================== SOCIAL DISTANCING PARAMETERS + +[Social distancing start time] +6 + +[Delay to social distancing by admin unit] +0 0 0 + +[Duration of social distancing] +364 + +[Duration of social distancing by admin unit] +364 364 364 + +[Trigger incidence per cell for social distancing] +0 + +[Relative place contact rate given social distancing by place type] +1 1 0.5 0.5 + +[Relative household contact rate given social distancing] +1.25 + +[Relative spatial contact rate given social distancing] +0.1 + +[Minimum radius for social distancing] +1 + +[Proportion compliant with enhanced social distancing] +0.0 + +[Proportion compliant with enhanced social distancing by age group] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 - -[Relative place contact rate given enhanced social distancing by place type] -0.25 0.25 0.25 0.25 - -[Relative household contact rate given enhanced social distancing] -1 - -[Relative spatial contact rate given enhanced social distancing] -0.25 - -[Delay for change in effectiveness of social distancing] -1000 - -[Relative place contact rate given social distancing by place type after change] -1 1 0.75 0.75 - -[Relative household contact rate given social distancing after change] -1.25 - -[Relative spatial contact rate given social distancing after change] -0.25 - -"//// Must match ""Change times for levels of social distancing""" -[Number of change times for levels of social distancing] -6 - -"//// Note: numbers here must match ""Number of change times for levels of social distancing""; that any times listed here that are before ""Social distancing start time"" and after ""Duration of social distancing"" are irrelevant." -[Change times for levels of social distancing] -0 31 61 91 121 151 - -"//// Again, want this to supercede ""Relative place contact rate given social distancing by place type"". Should be matrix of dimension ""Number of change times for levels of social distancing"" by Number of place types." -[Relative place contact rates over time given social distancing by place type] -1 1 0.5 0.5 -1 1 0.5 0.5 -1 1 0.5 0.5 -1 1 0.5 0.5 -1 1 0.5 0.5 -1 1 0.5 0.5 - -"//// Ideally want this to supercede ""Relative household contact rate given social distancing"", but need to preserved backwards compatibility for now." -[Relative household contact rates over time given social distancing] -1.25 1.25 1.25 1.25 1.25 1.25 - -"//// Ideally want this to supercede ""Relative spatial contact rate given social distancing"", but need to preserved backwards compatibility for now." -[Relative spatial contact rates over time given social distancing] -0.1 0.1 0.1 0.1 0.1 0.1 - -[Relative household contact rates over time given enhanced social distancing] -1 1 1 1 1 1 - -[Relative spatial contact rates over time given enhanced social distancing] -0.25 0.25 0.25 0.25 0.25 0.25 - -[Relative place contact rates over time given enhanced social distancing by place type] -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 -0.25 0.25 0.25 0.25 - -[Trigger incidence per cell for social distancing over time] -100 100 100 100 100 100 + +[Relative place contact rate given enhanced social distancing by place type] +0.25 0.25 0.25 0.25 + +[Relative household contact rate given enhanced social distancing] +1 + +[Relative spatial contact rate given enhanced social distancing] +0.25 + +[Delay for change in effectiveness of social distancing] +1000 + +[Relative place contact rate given social distancing by place type after change] +1 1 0.75 0.75 + +[Relative household contact rate given social distancing after change] +1.25 + +[Relative spatial contact rate given social distancing after change] +0.25 + +//// Must match "Change times for levels of social distancing" +[Number of change times for levels of social distancing] +6 + +//// Note: numbers here must match "Number of change times for levels of social distancing"; that any times listed here that are before "Social distancing start time" and after "Duration of social distancing" are irrelevant. +[Change times for levels of social distancing] +0 31 61 91 121 151 + +//// Again, want this to supercede "Relative place contact rate given social distancing by place type". Should be matrix of dimension "Number of change times for levels of social distancing" by Number of place types. +[Relative place contact rates over time given social distancing by place type] +1 1 0.5 0.5 +1 1 0.5 0.5 +1 1 0.5 0.5 +1 1 0.5 0.5 +1 1 0.5 0.5 +1 1 0.5 0.5 + +//// Ideally want this to supercede "Relative household contact rate given social distancing", but need to preserved backwards compatibility for now. +[Relative household contact rates over time given social distancing] +1.25 1.25 1.25 1.25 1.25 1.25 + +//// Ideally want this to supercede "Relative spatial contact rate given social distancing", but need to preserved backwards compatibility for now. +[Relative spatial contact rates over time given social distancing] +0.1 0.1 0.1 0.1 0.1 0.1 + +[Relative household contact rates over time given enhanced social distancing] +1 1 1 1 1 1 + +[Relative spatial contact rates over time given enhanced social distancing] +0.25 0.25 0.25 0.25 0.25 0.25 + +[Relative place contact rates over time given enhanced social distancing by place type] +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 +0.25 0.25 0.25 0.25 + +[Trigger incidence per cell for social distancing over time] +100 100 100 100 100 100 diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index 71fa9de92..37e001b83 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -95,7 +95,7 @@ void ImportMobilityDataToBetasArray(std::string MobilityFilename, int SettingNum { std::ifstream MobFileStream(MobilityFilename.c_str()); double ProportionalChange = 0.0; - for (int Day = 0; Day < P.SimulationDuration; Day++) + for (int Day = 0; Day < P.SimulationDuration + 1; Day++) for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) { MobFileStream >> ProportionalChange; @@ -472,10 +472,7 @@ int main(int argc, char* argv[]) P.NRactual = P.NRactNE; TSMean = TSMeanNE; TSVar = TSVarNE; - if ((P.DoRecordInfEvents) && (P.RecordInfEventsPerRun == 0)) - { - SaveEvents(output_file); - } + if ((P.DoRecordInfEvents) && (P.RecordInfEventsPerRun == 0)) SaveEvents(output_file); SaveSummaryResults(output_file); P.NRactual = P.NRactE; @@ -1251,7 +1248,13 @@ void InitModel(int run) // passing run number so we can save run number in the i //// Add all of the above to P.Efficacies array. UpdateEfficacyArray(); //// InitializeBetasArray (either to be same for all days and regions), or otherwise depending on how we are modelling them). Memory allocated in main using - MultiplyMobDataWithTransmissionParams(); + MultiplyMobDataWithTransmissionParams(); + + //for (int Setting = 0; Setting < P.NumInfectionSettings; Setting++) + // for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) + // for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) + // std::cout << "Day " << Day << ", AdUnit " << AdUnit << ", Setting " << Setting << ", Beta " << P.Betas[Day][AdUnit][Setting] << std::endl; ; + // Initialize CFR scalings P.CFR_Critical_Scale_Current = P.CFR_TimeScaling_Critical[0]; diff --git a/src/Param.h b/src/Param.h index 4b02b9d12..5ae4a6dd2 100644 --- a/src/Param.h +++ b/src/Param.h @@ -118,11 +118,12 @@ struct Param /**< 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]; From ab39a1ac7f18e655cb9930aaf5188537e28208c9 Mon Sep 17 00:00:00 2001 From: dlaydon Date: Thu, 2 Feb 2023 11:55:39 +0000 Subject: [PATCH 14/20] changed 'run' arg to 'realisation' in InitModel to match main, split up memory allocation and initialization of betas array --- src/CovidSim.cpp | 37 ++++++++++++++++++++++++++++--------- 1 file changed, 28 insertions(+), 9 deletions(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index 37e001b83..64ea374ef 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -83,11 +83,14 @@ void AllocateMemForBetasArray() // called in main (once per fitting run). Can yo for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) P.Betas[Day][AdUnit] = new double[P.NumInfectionSettings](); } +} +void InitializeBetasToOne() +{ // initialize all values to 1; for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) for (int Setting = 0; Setting < P.NumInfectionSettings; Setting++) - P.Betas[Day][AdUnit][Setting] = 1; + P.Betas[Day][AdUnit][Setting] = 1.0; } void ImportMobilityDataToBetasArray(std::string MobilityFilename, int SettingNumber) // assumes all P.Betas values initialized to 1. { @@ -105,7 +108,7 @@ void ImportMobilityDataToBetasArray(std::string MobilityFilename, int SettingNum } else ERR_CRITICAL(("mobility_file: " + MobilityFilename + " does not exist").c_str()); } -void MultiplyMobDataWithTransmissionParams() +void MultiplyBetasArrayWithTransmissionParams() { for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) @@ -119,6 +122,20 @@ void MultiplyMobDataWithTransmissionParams() P.Betas[Day][AdUnit][Spatial] *= P.LocalBeta; } } +void DivideBetasArrayWithTransmissionParams() +{ + for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) + for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) + { + // place (by type) + for (int PlaceType = 0; PlaceType < P.NumPlaceTypes; PlaceType++) + P.Betas[Day][AdUnit][PlaceType] /= P.PlaceTypeTrans[PlaceType]; + // Household + P.Betas[Day][AdUnit][House] /= P.HouseholdTrans; + // Spatial/Community + P.Betas[Day][AdUnit][Spatial] /= P.LocalBeta; + } +} ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** GLOBAL VARIABLES (some structures in CovidSim.h file and some containers) - memory allocated later. @@ -355,13 +372,15 @@ int main(int argc, char* argv[]) // Allocate memory for Betas array (dynamic allocation - must be done after P.NumAdunits set in SetupModel.cpp::SetupPopulation) // By default, make all multipliers of betas equal to 1. If P.VaryBetasOverTimeByRegion, these multipliers will be mobility data - AllocateMemForBetasArray(); + AllocateMemForBetasArray(); + InitializeBetasToOne(); if (P.VaryBetasOverTimeByRegion) { ImportMobilityDataToBetasArray(workplace_mob_file , Workplace); ImportMobilityDataToBetasArray(household_mob_file , House); ImportMobilityDataToBetasArray(spatial_mob_file , Spatial); } + MultiplyBetasArrayWithTransmissionParams(); //for (int Setting = 0; Setting < P.NumInfectionSettings; Setting++) // for (int Day = 0; Day < (int)P.SimulationDuration; Day++) // for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) @@ -933,7 +952,7 @@ void UpdateEfficacyArray() P.Efficacies[DigContactTracing][Spatial ] = P.DCTCaseIsolationEffectiveness; } -void InitModel(int run) // passing run number so we can save run number in the infection event log: ggilani - 15/10/2014 +void InitModel(int Realisation) // passing run number so we can save run number in the infection event log: ggilani - 15/10/2014 { int NumImmune = 0; @@ -1006,7 +1025,7 @@ void InitModel(int run) // passing run number so we can save run number in the i State.cumDC_adunit[i] = State.cumCT_adunit[i] = State.cumCC_adunit[i] = State.trigDC_adunit[i] = State.DCT_adunit[i] = State.cumDCT_adunit[i] = 0; //added hospitalisation, added detected cases, contact tracing per adunit, cases who are contacts: ggilani 03/02/15, 15/06/17 AdUnits[i].place_close_trig = 0; AdUnits[i].CaseIsolationTimeStart = AdUnits[i].HQuarantineTimeStart = AdUnits[i].DigitalContactTracingTimeStart = AdUnits[i].SocialDistanceTimeStart = AdUnits[i].PlaceCloseTimeStart = 1e10; - AdUnits[i].ndct = 0; //noone being digitally contact traced at beginning of run + AdUnits[i].ndct = 0; //noone being digitally contact traced at beginning of Realisation } //update state variables for storing contact distribution @@ -1248,7 +1267,7 @@ void InitModel(int run) // passing run number so we can save run number in the i //// Add all of the above to P.Efficacies array. UpdateEfficacyArray(); //// InitializeBetasArray (either to be same for all days and regions), or otherwise depending on how we are modelling them). Memory allocated in main using - MultiplyMobDataWithTransmissionParams(); + if (Realisation == 0) MultiplyBetasArrayWithTransmissionParams(); //for (int Setting = 0; Setting < P.NumInfectionSettings; Setting++) // for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) @@ -1272,7 +1291,7 @@ void InitModel(int run) // passing run number so we can save run number in the i UpdateProbs(0); DoInitUpdateProbs = 0; } - //initialise event log to zero at the beginning of every run: ggilani - 10/10/2014. UPDATE: 15/10/14 - we are now going to store all events from all realisations in one file + //initialise event log to zero at the beginning of every Realisation: ggilani - 10/10/2014. UPDATE: 15/10/14 - we are now going to store all events from all realisations in one file if ((P.DoRecordInfEvents) && (P.RecordInfEventsPerRun)) { nEvents = 0; @@ -1286,7 +1305,7 @@ void InitModel(int run) // passing run number so we can save run number in the i int* NumSeedingInfections_byLocation = new int[P.NumSeedLocations]; for (int i = 0; i < P.NumSeedLocations; i++) NumSeedingInfections_byLocation[i] = (int) (((double) P.NumInitialInfections[i]) * P.InitialInfectionsAdminUnitWeight[i]* P.SeedingScaling +0.5); - SeedInfection(0, NumSeedingInfections_byLocation, 0, run); + SeedInfection(0, NumSeedingInfections_byLocation, 0, Realisation); delete[] NumSeedingInfections_byLocation; P.ControlPropCasesId = P.PreAlertControlPropCasesId; P.TreatTimeStart = 1e10; @@ -1306,7 +1325,7 @@ void InitModel(int run) // passing run number so we can save run number in the i P.PlaceCloseIncTrig = P.PlaceCloseIncTrig1; P.PlaceCloseTimeStartPrevious = 1e10; P.PlaceCloseCellIncThresh = P.PlaceCloseCellIncThresh1; - P.ResetSeedsFlag = 0; //added this to allow resetting seeds part way through run: ggilani 27/11/2019 + P.ResetSeedsFlag = 0; //added this to allow resetting seeds part way through Realisation: ggilani 27/11/2019 if (!P.StopCalibration) P.DateTriggerReached_SimTime = 0; if (P.InitialInfectionCalTime > 0) { From 4166937717a30f5ecaa05833956d7a52797effb9 Mon Sep 17 00:00:00 2001 From: dlaydon Date: Thu, 2 Feb 2023 12:34:57 +0000 Subject: [PATCH 15/20] moved multiplying betas array outside INitMOdel --- src/CovidSim.cpp | 12 ++++-------- src/SetupModel.cpp | 2 ++ 2 files changed, 6 insertions(+), 8 deletions(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index 64ea374ef..85a9d66da 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -106,7 +106,7 @@ void ImportMobilityDataToBetasArray(std::string MobilityFilename, int SettingNum } MobFileStream.close(); } - else ERR_CRITICAL(("mobility_file: " + MobilityFilename + " does not exist").c_str()); + else ERR_CRITICAL(("mobility_file = " + MobilityFilename + " - file does not exist").c_str()); } void MultiplyBetasArrayWithTransmissionParams() { @@ -381,6 +381,7 @@ int main(int argc, char* argv[]) ImportMobilityDataToBetasArray(spatial_mob_file , Spatial); } MultiplyBetasArrayWithTransmissionParams(); + //for (int Setting = 0; Setting < P.NumInfectionSettings; Setting++) // for (int Day = 0; Day < (int)P.SimulationDuration; Day++) // for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) @@ -408,8 +409,10 @@ int main(int argc, char* argv[]) StopFit = ReadFitIter(fit_file); if (!StopFit) { + DivideBetasArrayWithTransmissionParams(); Params::ReadParams(param_file, pre_param_file, ad_unit_file, &P, AdUnits); if (!P.FixLocalBeta) InitTransmissionCoeffs(); + MultiplyBetasArrayWithTransmissionParams(); output_file_base = output_file_base_f + ".f" + std::to_string(P.FitIter); } } @@ -1266,13 +1269,6 @@ void InitModel(int Realisation) // passing run number so we can save run number //// Add all of the above to P.Efficacies array. UpdateEfficacyArray(); - //// InitializeBetasArray (either to be same for all days and regions), or otherwise depending on how we are modelling them). Memory allocated in main using - if (Realisation == 0) MultiplyBetasArrayWithTransmissionParams(); - - //for (int Setting = 0; Setting < P.NumInfectionSettings; Setting++) - // for (int Day = 0; Day < (int)P.SimulationDuration + 1; Day++) - // for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++) - // std::cout << "Day " << Day << ", AdUnit " << AdUnit << ", Setting " << Setting << ", Beta " << P.Betas[Day][AdUnit][Setting] << std::endl; ; // Initialize CFR scalings diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp index 5e08cbb7f..1e1f5049b 100644 --- a/src/SetupModel.cpp +++ b/src/SetupModel.cpp @@ -661,6 +661,7 @@ void InitTransmissionCoeffs(void) } Files::xfprintf_stderr("Household mean size = %lg\n", HouseholdMeanSize); + //// Loops below sum household and spatial infections double HH_SAR_Denom = 0.0; // household secondary-attack rate denominator. Will sum over following #pragma loop #pragma omp parallel for private(Household_Infectiousness,Spatial_Infectiousness,quantile,LatentToSympDelay,ProbSurvive) schedule(static,1) reduction(+:HH_Infections,SpatialInfections,HH_SAR_Denom) default(none) shared(P, Households, Hosts, Mcells) @@ -759,6 +760,7 @@ void InitTransmissionCoeffs(void) P.R0household = HH_Infections / ((double)P.PopSize); Files::xfprintf_stderr("Household R0 = %lg\n", P.R0household); + // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** // ** Place Infections if (P.DoPlaces) From 8ab0ff555a06e37dbb870dcea71784019a451d6e Mon Sep 17 00:00:00 2001 From: Wes Hinsley Date: Thu, 2 Feb 2023 14:49:30 +0000 Subject: [PATCH 16/20] Attempt to fix MacOS CI --- ci/install_macos_dependencies.sh | 1 + 1 file changed, 1 insertion(+) diff --git a/ci/install_macos_dependencies.sh b/ci/install_macos_dependencies.sh index 6f33f619c..eb7bbd5e1 100755 --- a/ci/install_macos_dependencies.sh +++ b/ci/install_macos_dependencies.sh @@ -1,3 +1,4 @@ #!/bin/sh -e brew install llvm cmake libomp python3 coreutils +export CMAKE_CXX_FLAGS="$CMAKE_CXX_FLAGS -Xpreprocessor -fopenmp From 92ca177cd5f0357be90dee36ac8489e1b3770135 Mon Sep 17 00:00:00 2001 From: Wes Hinsley Date: Thu, 2 Feb 2023 14:51:25 +0000 Subject: [PATCH 17/20] Attempt to fix MacOS CI --- ci/install_macos_dependencies.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ci/install_macos_dependencies.sh b/ci/install_macos_dependencies.sh index eb7bbd5e1..89cf3014d 100755 --- a/ci/install_macos_dependencies.sh +++ b/ci/install_macos_dependencies.sh @@ -1,4 +1,4 @@ #!/bin/sh -e brew install llvm cmake libomp python3 coreutils -export CMAKE_CXX_FLAGS="$CMAKE_CXX_FLAGS -Xpreprocessor -fopenmp +export CMAKE_CXX_FLAGS="$CMAKE_CXX_FLAGS -Xpreprocessor -fopenmp" From 299378abc9b870bf8fd43d1818fa76f9f3824e3d Mon Sep 17 00:00:00 2001 From: Wes Hinsley Date: Thu, 2 Feb 2023 14:55:35 +0000 Subject: [PATCH 18/20] Attempt to fix MacOS CI --- ci/install_macos_dependencies.sh | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ci/install_macos_dependencies.sh b/ci/install_macos_dependencies.sh index 89cf3014d..898be07fe 100755 --- a/ci/install_macos_dependencies.sh +++ b/ci/install_macos_dependencies.sh @@ -2,3 +2,5 @@ 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" From ab98db75f7dd0e43f556807f6e937e143b47a3ad Mon Sep 17 00:00:00 2001 From: Wes Hinsley Date: Thu, 2 Feb 2023 15:04:05 +0000 Subject: [PATCH 19/20] Attempt to fix MacOS CI --- CMakeLists.txt | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 2218a8b31..37e64b8a7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -91,3 +91,18 @@ if(DOXYGEN_FOUND) elseif(DOXYGEN_FOUND) message("Doxygen needs to be installed to generate the doxygen documentation") endif(DOXYGEN_FOUND) + +if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin") + set(OpenMP_C "${CMAKE_C_COMPILER}") + set(OpenMP_C_FLAGS "-Xclang -fopenmp -I/usr/local/Cellar/libomp/14.0.0/include") + set(OpenMP_C_LIB_NAMES "libomp") + set(OpenMP_libomp_LIBRARY "omp") + + set(OpenMP_CXX "${CMAKE_CXX_COMPILER}") + set(OpenMP_CXX_FLAGS "-Xclang -fopenmp -I/usr/local/Cellar/libomp/14.0.0/include") + set(OpenMP_CXX_LIB_NAMES "libomp") + set(OpenMP_libomp_LIBRARY "omp") + + link_directories("/usr/local/Cellar/libomp/14.0.0/lib/") + target_link_directories(${PROJECT_NAME} PRIVATE "/usr/local/Cellar/libomp/14.0.0/lib/") +endif() \ No newline at end of file From 4297928a284cbfb25cf50a6666df3946f330d486 Mon Sep 17 00:00:00 2001 From: Wes Hinsley Date: Thu, 2 Feb 2023 15:08:30 +0000 Subject: [PATCH 20/20] Reverting --- CMakeLists.txt | 15 --------------- 1 file changed, 15 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 37e64b8a7..2218a8b31 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -91,18 +91,3 @@ if(DOXYGEN_FOUND) elseif(DOXYGEN_FOUND) message("Doxygen needs to be installed to generate the doxygen documentation") endif(DOXYGEN_FOUND) - -if (${CMAKE_SYSTEM_NAME} MATCHES "Darwin") - set(OpenMP_C "${CMAKE_C_COMPILER}") - set(OpenMP_C_FLAGS "-Xclang -fopenmp -I/usr/local/Cellar/libomp/14.0.0/include") - set(OpenMP_C_LIB_NAMES "libomp") - set(OpenMP_libomp_LIBRARY "omp") - - set(OpenMP_CXX "${CMAKE_CXX_COMPILER}") - set(OpenMP_CXX_FLAGS "-Xclang -fopenmp -I/usr/local/Cellar/libomp/14.0.0/include") - set(OpenMP_CXX_LIB_NAMES "libomp") - set(OpenMP_libomp_LIBRARY "omp") - - link_directories("/usr/local/Cellar/libomp/14.0.0/lib/") - target_link_directories(${PROJECT_NAME} PRIVATE "/usr/local/Cellar/libomp/14.0.0/lib/") -endif() \ No newline at end of file