diff --git a/ci/install_macos_dependencies.sh b/ci/install_macos_dependencies.sh
index 6f33f619c..898be07fe 100755
--- a/ci/install_macos_dependencies.sh
+++ b/ci/install_macos_dependencies.sh
@@ -1,3 +1,6 @@
#!/bin/sh -e
brew install llvm cmake libomp python3 coreutils
+export CMAKE_CXX_FLAGS="$CMAKE_CXX_FLAGS -Xpreprocessor -fopenmp"
+export OpenMP_C_LIB_NAMES="libomp"
+export OpenMP_C_FLAGS="-fopenmp"
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/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/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 b480596a0..85a9d66da 100644
--- a/src/CovidSim.cpp
+++ b/src/CovidSim.cpp
@@ -67,6 +67,76 @@ 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
+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]();
+ }
+}
+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.0;
+}
+void ImportMobilityDataToBetasArray(std::string MobilityFilename, int SettingNumber) // assumes all P.Betas values initialized to 1.
+{
+ if (MobilityFilename != "NoFileFound" && FileExists(MobilityFilename))
+ {
+ std::ifstream MobFileStream(MobilityFilename.c_str());
+ double ProportionalChange = 0.0;
+ for (int Day = 0; Day < P.SimulationDuration + 1; Day++)
+ for (int AdUnit = 0; AdUnit < P.NumAdunits; AdUnit++)
+ {
+ MobFileStream >> ProportionalChange;
+ P.Betas[Day][AdUnit][SettingNumber] += ProportionalChange;
+ }
+ MobFileStream.close();
+ }
+ else ERR_CRITICAL(("mobility_file = " + MobilityFilename + " - file does not exist").c_str());
+}
+void MultiplyBetasArrayWithTransmissionParams()
+{
+ 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;
+ }
+}
+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.
///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** ///// ***** /////
@@ -125,6 +195,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.
@@ -146,7 +218,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)
@@ -168,7 +240,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");
@@ -176,6 +248,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:
@@ -261,9 +338,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 +362,30 @@ 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);
+ 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
- P.NumInterventionClasses = 6; // CI, HQ, PC, SD, Enhanced Social Distancing, DCT
- P.Efficacies = new double*[P.NumInterventionClasses]();
+ 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)
+ // By default, make all multipliers of betas equal to 1. If P.VaryBetasOverTimeByRegion, these multipliers will be mobility data
+ 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++)
+ // std::cout << "Day " << Day << ", AdUnit " << AdUnit << ", Setting " << Setting << ", Beta " << P.Betas[Day][AdUnit][Setting] << std::endl; ;
//// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// ****
//// **** RUN MODEL
@@ -319,9 +409,11 @@ 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();
- output_file_base = output_file_base_f + ".f" + std::to_string(P.FitIter);
+ MultiplyBetasArrayWithTransmissionParams();
+ output_file_base = output_file_base_f + ".f" + std::to_string(P.FitIter);
}
}
else StopFit = 1;
@@ -335,7 +427,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))
@@ -349,6 +441,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))
@@ -401,10 +494,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;
@@ -415,7 +505,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");
}
}
@@ -865,9 +955,9 @@ 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 nim;
+ int NumImmune = 0;
if (P.OutputBitmap)
{
@@ -938,7 +1028,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
@@ -991,7 +1081,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++)
@@ -1039,7 +1128,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++)
{
@@ -1073,7 +1162,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);
}
@@ -1085,7 +1174,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;
}
}
}
@@ -1181,6 +1270,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();
+
// Initialize CFR scalings
P.CFR_Critical_Scale_Current = P.CFR_TimeScaling_Critical[0];
P.CFR_SARI_Scale_Current = P.CFR_TimeScaling_SARI [0];
@@ -1197,7 +1287,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;
@@ -1211,7 +1301,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;
@@ -1231,7 +1321,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)
{
@@ -3688,73 +3778,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"
{
@@ -3765,12 +3857,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)
{
@@ -3781,13 +3875,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);
}
}
}
@@ -3800,9 +3894,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
@@ -3810,7 +3904,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
}
}
}
@@ -3821,32 +3915,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/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..5ae4a6dd2 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,13 +109,23 @@ struct Param
CovidSim::Geometry::BoundingBox2d SpatialBoundingBox;
double** LocationInitialInfection;
- double InitialInfectionsAdminUnitWeight[MAX_NUM_SEED_LOCATIONS], InitialInfectionCalTime, TimeStepsPerDay;
+ double InitialInfectionsAdminUnitWeight[MAX_NUM_SEED_LOCATIONS], InitialInfectionCalTime;
double FalsePositiveRate, FalsePositivePerCapitaIncidence, FalsePositiveAgeRate[NUM_AGE_GROUPS];
double SeroConvMaxSens, SeroConvP1, SeroConvP2, SeroConvSpec, InfPrevSurveyScale;
double AirportTrafficScale;
- double R0, R0scale, LocalBeta;
+ /**< Betas array (average number of infectious contacts per day at time t in given region for particular setting)
+ // by default, betas do not vary over time or between regions, unless P.VaryBetasOverTimeByRegion == true
+ // setting refers to either place, household and spatial/community
+ // Flow is to read in mobility data (proportional change in mobility type by region over time).
+ // If VaryBetasOverTimeByRegion == false, rather than read in mobility data, first set all betas to 1, then multipy each day and regions value by P.LocalBeta, P.HouseholdTrans etc.
+ */
+ double*** Betas; //// indexed by i) time; ii) region / admin unit; iii) setting.
+ bool VaryBetasOverTimeByRegion = false;
+
+ double R0, R0scale;
+ double LocalBeta; //// Spatial/community beta: (roughly the average number of infectious spatial/community contacts per person per unit time)
double infectious_prof[INFPROF_RES + 1], infectiousness[MAX_INFECTIOUS_STEPS];
double InfectiousPeriod; // In days. Mean of icdf (inverse cumulative distribution function).
double R0household, R0places, R0spatial;
@@ -140,7 +151,8 @@ struct Param
int AbsenteeismPlaceClosure; // Default 0 (off), 1 (on) track place closures in more detail
int MaxAbsentTime; // In days. Max number of days absent, range [0, MAX_ABSENT_TIME]. Default 0 if !P.AbsenteeismPlaceClosure, otherwise MAX_ABSENT_TIME
int InvJourneyDurationDistrib[1025], InvLocalJourneyDurationDistrib[1025];
- double HouseholdTrans, HouseholdTransPow;
+ double HouseholdTrans; //// household beta (roughly the average number of infectious household contacts per person per unit time)
+ double HouseholdTransPow;
double** HouseholdSizeDistrib; // [MAX_ADUNITS] [MAX_HOUSEHOLD_SIZE]
double HouseholdDenomLookup[MAX_HOUSEHOLD_SIZE];
int PlaceTypeAgeMin[MAX_NUM_PLACE_TYPES], PlaceTypeAgeMax[MAX_NUM_PLACE_TYPES], PlaceTypeMaxAgeRead[MAX_NUM_PLACE_TYPES];
@@ -149,7 +161,8 @@ struct Param
int PlaceTypeNearestNeighb[MAX_NUM_PLACE_TYPES], PlaceTypeKernelType[MAX_NUM_PLACE_TYPES];
double PlaceTypePropAgeGroup[MAX_NUM_PLACE_TYPES], PlaceTypePropAgeGroup2[MAX_NUM_PLACE_TYPES];
double PlaceTypePropAgeGroup3[MAX_NUM_PLACE_TYPES], PlaceTypeKernelShape[MAX_NUM_PLACE_TYPES], PlaceTypeKernelScale[MAX_NUM_PLACE_TYPES];
- double PlaceTypeKernelP3[MAX_NUM_PLACE_TYPES], PlaceTypeKernelP4[MAX_NUM_PLACE_TYPES], PlaceTypeTrans[MAX_NUM_PLACE_TYPES];
+ double PlaceTypeKernelP3[MAX_NUM_PLACE_TYPES], PlaceTypeKernelP4[MAX_NUM_PLACE_TYPES];
+ double PlaceTypeTrans[MAX_NUM_PLACE_TYPES]; //// indexed by place time place type betas (roughly average number of infectious place contacts per person per unit time).
double PlaceTypeMeanSize[MAX_NUM_PLACE_TYPES], PlaceTypePropBetweenGroupLinks[MAX_NUM_PLACE_TYPES], PlaceTypeSizeSD[MAX_NUM_PLACE_TYPES]; //added PlaceTypeSizeSD for lognormal distribution - ggilani 09/02/17
double PlaceTypeSizePower[MAX_NUM_PLACE_TYPES], PlaceTypeSizeOffset[MAX_NUM_PLACE_TYPES], PlaceTypeSizeMax[MAX_NUM_PLACE_TYPES], PlaceTypeSizeMin[MAX_NUM_PLACE_TYPES];
double PlaceTypeGroupSizeParam1[MAX_NUM_PLACE_TYPES], PlaceExclusivityMatrix[MAX_NUM_PLACE_TYPES * MAX_NUM_PLACE_TYPES]; //changed PlaceExclusivityMatrix from [MAX_NUM_PLACE_TYPES][MAX_NUM_PLACE_TYPES]
@@ -263,6 +276,7 @@ struct Param
int NumInterventionClasses = 0, NumInfectionSettings = 0; // initialized in CovidSim.cpp::main.
+
///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****
///// **** VARIABLE EFFICACIES OVER TIME
///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// **** ///// ****
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");
diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp
index 3372a6bef..1e1f5049b 100644
--- a/src/SetupModel.cpp
+++ b/src/SetupModel.cpp
@@ -661,35 +661,38 @@ 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)
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[Hosts[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 +701,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 +716,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 +738,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[AdUnit][InfecteeAge] * P.WAIFW_Matrix_SpatialOnly[InfecteeAge][HOST_AGE_GROUP(person)];
// scale spatial infectiousness by weighted average.
Spatial_Infectiousness *= AvContactRate_Infector;
}
@@ -744,9 +747,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.
@@ -757,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)
@@ -765,13 +769,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 +785,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 +799,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 +807,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 +822,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.
diff --git a/src/Sweep.cpp b/src/Sweep.cpp
index 9aef2fd73..802446a44 100644
--- a/src/Sweep.cpp
+++ b/src/Sweep.cpp
@@ -292,10 +292,9 @@ 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 +302,33 @@ 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, seasonality, TimeStepNow, fp, BlanketMoveRestrInPlace, stderr_shared, Day)
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)
{
- 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;
+ 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:
@@ -325,11 +337,13 @@ 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)
- && (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
-
+ Household_Beta = (P.DoHouseholds) ? P.Betas[Day][AdUnit_ThisPerson][House] * seasonality * fp : 0;
if (Household_Beta > 0)
{
// For InfectiousPerson's household (InfectiousPerson->hh),
@@ -353,14 +367,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 +389,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
{
@@ -390,7 +404,8 @@ 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.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
@@ -405,32 +420,31 @@ 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 +460,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 +476,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,26 +497,24 @@ 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 +561,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 +583,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 +616,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 +634,44 @@ 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)
+
+ SpatialSeasonal_Beta = seasonality * fp * P.Betas[Day][AdUnit_ThisPerson][Spatial];
+ 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 +703,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 +712,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 +720,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 +739,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 +763,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 +780,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 +828,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 +858,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 +869,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) \