From 616a814ee94497c540c9d7b94895a8f17584eb35 Mon Sep 17 00:00:00 2001 From: Ian Lynagh Date: Thu, 23 Apr 2020 14:54:13 +0100 Subject: [PATCH] Refactor setting of RNG seeds There were a couple of places where we use P.newseed1 = (int)(ranf() * 1e8); to make new seeds, but Rand's setall already implements RNG splitting, so we now use that instead. This needed a few knock-on changes in the places where we seed the RNG, as if we wanted to reuse the seed values then we now need to store them. Finally, we now also split the seeds used for the setup phase, and reinitialise the RNG with them after the network has been created or loaded. It looks like historically the second pair of seeds had been used at this point, to make the runs identical regardless of how the network was made, but that this had been changed when seed-resetting was implemented. Now the determinism is back, fixing #116. --- src/CovidSim.cpp | 68 +++++++++++++++++----------------------------- src/Param.h | 6 ++-- src/Rand.cpp | 16 +++++++---- src/Rand.h | 2 +- src/SetupModel.cpp | 14 ++++++---- 5 files changed, 49 insertions(+), 57 deletions(-) diff --git a/src/CovidSim.cpp b/src/CovidSim.cpp index abe79761e..f3b247554 100644 --- a/src/CovidSim.cpp +++ b/src/CovidSim.cpp @@ -133,10 +133,10 @@ int main(int argc, char* argv[]) { ///// Get seeds. i = argc - 4; - sscanf(argv[i], "%li", &P.seed1); - sscanf(argv[i + 1], "%li", &P.seed2); - sscanf(argv[i + 2], "%li", &P.seed3); - sscanf(argv[i + 3], "%li", &P.seed4); + sscanf(argv[i], "%li", &P.setupSeed1); + sscanf(argv[i + 1], "%li", &P.setupSeed2); + sscanf(argv[i + 2], "%li", &P.runSeed1); + sscanf(argv[i + 3], "%li", &P.runSeed2); ///// Set parameter defaults - read them in after P.PlaceCloseIndepThresh = P.LoadSaveNetwork = P.DoHeteroDensity = P.DoPeriodicBoundaries = P.DoSchoolFile = P.DoAdunitDemog = P.OutputDensFile = P.MaxNumThreads = P.DoInterventionFile = 0; @@ -347,14 +347,6 @@ int main(int argc, char* argv[]) //print out number of calls to random number generator - ///// Set seeds - if (!P.ResetSeeds) - { - P.newseed1 = P.seed3; - P.newseed2 = P.seed4; - setall(P.newseed1, P.newseed2); - } - //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** RUN MODEL //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** @@ -370,31 +362,21 @@ int main(int argc, char* argv[]) } ///// Set and save seeds - if (P.ResetSeeds) + if (i == 0 || P.ResetSeeds && P.KeepSameSeeds) { - if (P.KeepSameSeeds) - { - P.newseed1 = P.seed3; - P.newseed2 = P.seed4; - } - else - { - if (i == 0) - { - P.newseed1 = P.seed3; - P.newseed2 = P.seed4; - } - else - { - //sample 2 new seeds for random number generators - P.newseed1 = (int)(ranf() * 1e8); - P.newseed2 = (int)(ranf() * 1e8); - } - } + P.nextRunSeed1 = P.runSeed1; + P.nextRunSeed2 = P.runSeed2; + } + if (P.ResetSeeds) { //save these seeds to file SaveRandomSeeds(); - //reset seeds - setall(P.newseed1, P.newseed2); + } + // Now that we have set P.nextRunSeed* ready for the run, we need to save the values in case + // we need to reinitialise the RNG after the run is interrupted. + long thisRunSeed1 = P.nextRunSeed1; + long thisRunSeed2 = P.nextRunSeed2; + if (i == 0 || P.ResetSeeds) { + setall(&P.nextRunSeed1, &P.nextRunSeed2); //fprintf(stderr, "%i, %i\n", P.newseed1,P.newseed2); //fprintf(stderr, "%f\n", ranf()); } @@ -404,7 +386,9 @@ int main(int argc, char* argv[]) if (P.DoLoadSnapshot) LoadSnapshot(); while (RunModel(i)) { // 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. - setall(P.newseed1, P.newseed2); // reset random number seeds to generate same run again after calibration. + long tmp1 = thisRunSeed1; + long tmp2 = thisRunSeed2; + setall(&tmp1, &tmp2); // reset random number seeds to generate same run again after calibration. InitModel(i); } if (P.OutputNonSummaryResults) @@ -2928,9 +2912,7 @@ int RunModel(int run) //added run number as parameter if (P.ResetSeedsPostIntervention) if ((P.ResetSeedsFlag == 0) && (ts >= (P.TimeToResetSeeds * P.TimeStepsPerDay))) { - P.newseed1 = (int)(ranf() * 1e8); - P.newseed2 = (int)(ranf() * 1e8); - setall(P.newseed1, P.newseed2); + setall(&P.nextRunSeed1, &P.nextRunSeed2); P.ResetSeedsFlag = 1; } @@ -3934,7 +3916,7 @@ void SaveRandomSeeds(void) sprintf(outname, "%s.seeds.xls", OutFile); if (!(dat = fopen(outname, "wb"))) ERR_CRITICAL("Unable to open output file\n"); - fprintf(dat, "%li\t%li\n", P.newseed1, P.newseed2); + fprintf(dat, "%li\t%li\n", P.nextRunSeed1, P.nextRunSeed2); fclose(dat); } @@ -3993,8 +3975,8 @@ void LoadSnapshot(void) fread_big((void*)& i, sizeof(int), 1, dat); if (i != P.NCP) ERR_CRITICAL("Incorrect NCP in snapshot file.\n"); fread_big((void*)& i, sizeof(int), 1, dat); if (i != P.ncw) ERR_CRITICAL("Incorrect ncw in snapshot file.\n"); fread_big((void*)& i, sizeof(int), 1, dat); if (i != P.nch) ERR_CRITICAL("Incorrect nch in snapshot file.\n"); - fread_big((void*)& l, sizeof(long), 1, dat); if (l != P.seed1) ERR_CRITICAL("Incorrect seed1 in snapshot file.\n"); - fread_big((void*)& l, sizeof(long), 1, dat); if (l != P.seed2) ERR_CRITICAL("Incorrect seed2 in snapshot file.\n"); + fread_big((void*)& l, sizeof(long), 1, dat); if (l != P.setupSeed1) ERR_CRITICAL("Incorrect setupSeed1 in snapshot file.\n"); + fread_big((void*)& l, sizeof(long), 1, dat); if (l != P.setupSeed2) ERR_CRITICAL("Incorrect setupSeed2 in snapshot file.\n"); fread_big((void*)& t, sizeof(double), 1, dat); if (t != P.TimeStep) ERR_CRITICAL("Incorrect TimeStep in snapshot file.\n"); fread_big((void*) & (P.SnapshotLoadTime), sizeof(double), 1, dat); P.NumSamples = 1 + (int)ceil((P.SampleTime - P.SnapshotLoadTime) / P.SampleStep); @@ -4067,9 +4049,9 @@ void SaveSnapshot(void) fprintf(stderr, "## %i\n", i++); fwrite_big((void*) & (P.nch), sizeof(int), 1, dat); fprintf(stderr, "## %i\n", i++); - fwrite_big((void*) & (P.seed1), sizeof(long), 1, dat); + fwrite_big((void*) & (P.setupSeed1), sizeof(long), 1, dat); fprintf(stderr, "## %i\n", i++); - fwrite_big((void*) & (P.seed2), sizeof(long), 1, dat); + fwrite_big((void*) & (P.setupSeed2), sizeof(long), 1, dat); fprintf(stderr, "## %i\n", i++); fwrite_big((void*) & (P.TimeStep), sizeof(double), 1, dat); fprintf(stderr, "## %i\n", i++); diff --git a/src/Param.h b/src/Param.h index cb2c56d85..c037368c0 100644 --- a/src/Param.h +++ b/src/Param.h @@ -50,8 +50,10 @@ typedef struct PARAM { int DoAirports, Nairports, Air_popscale, DoSchoolFile, DoRealSymptWithdrawal, CaseAbsentChildAgeCutoff, DoEarlyCaseDiagnosis, DoInterventionFile; int PlaceTypeNoAirNum; // If DoAirports then this is the number of non-airport place types (< PlaceTypeNum), else == PlaceTypeNum (~ no airport places). int HotelPlaceType; // If DoAirports then this is place type for hotel (>= PlaceTypeNoAirNum, < PlaceTypeNum), else == PlaceTypeNum (~ unused). - long seed1, seed2, seed3, seed4; - long newseed1, newseed2; //added these to allow for seeds to be reset - ggilani 09/03/17 + long setupSeed1, setupSeed2; // RNG seeds from the command line, used to initialise the RNG for setup + long runSeed1, runSeed2; // RNG seeds from the command line, used to initialise the RNG for running the model + long nextSetupSeed1, nextSetupSeed2; // The next RNG seeds to use when we need to reinitialise the RNG for setup + long nextRunSeed1, nextRunSeed2; // The next RNG seeds to use when we need to reinitialise the RNG for the model int ResetSeeds,KeepSameSeeds, ResetSeedsPostIntervention, ResetSeedsFlag, TimeToResetSeeds; double SpatialBoundingBox[4], LocationInitialInfection[MAX_NUM_SEED_LOCATIONS][2], InitialInfectionsAdminUnitWeight[MAX_NUM_SEED_LOCATIONS], TimeStepsPerDay; double FalsePositiveRate, FalsePositivePerCapitaIncidence, FalsePositiveAgeRate[NUM_AGE_GROUPS]; diff --git a/src/Rand.cpp b/src/Rand.cpp index a4faf8a3b..c0fa05e22 100644 --- a/src/Rand.cpp +++ b/src/Rand.cpp @@ -66,7 +66,7 @@ double ranf_mt(int tn) return ((double)z) / Xm1; } -void setall(long iseed1, long iseed2) +void setall(long *pseed1, long *pseed2) /* ********************************************************************** void setall(long iseed1,long iseed2) @@ -88,12 +88,16 @@ void setall(long iseed1, long iseed2) { int g; - *Xcg1 = iseed1; - *Xcg2 = iseed2; - for (g = 1; g < MAX_NUM_THREADS; g++) { - *(Xcg1 + g * CACHE_LINE_SIZE) = mltmod(Xa1vw, *(Xcg1 + (g - 1) * CACHE_LINE_SIZE), Xm1); - *(Xcg2 + g * CACHE_LINE_SIZE) = mltmod(Xa2vw, *(Xcg2 + (g - 1) * CACHE_LINE_SIZE), Xm2); + long iseed1 = *pseed1; + long iseed2 = *pseed2; + + for (g = 0; g < MAX_NUM_THREADS; g++) { + *(Xcg1 + g * CACHE_LINE_SIZE) = iseed1 = mltmod(Xa1vw, iseed1, Xm1); + *(Xcg2 + g * CACHE_LINE_SIZE) = iseed2 = mltmod(Xa2vw, iseed2, Xm2); } + + *pseed1 = iseed1; + *pseed2 = iseed2; } long mltmod(long a, long s, long m) /* diff --git a/src/Rand.h b/src/Rand.h index 6a9079d92..927e48fac 100644 --- a/src/Rand.h +++ b/src/Rand.h @@ -19,7 +19,7 @@ long ignbin_mt(long, double, int); long ignpoi_mt(double, int); double ranf(void); double ranf_mt(int); -void setall(long, long); +void setall(long *, long *); double sexpo_mt(int); double sexpo(void); long mltmod(long, long, long); diff --git a/src/SetupModel.cpp b/src/SetupModel.cpp index 5c673b593..8bd5971da 100644 --- a/src/SetupModel.cpp +++ b/src/SetupModel.cpp @@ -34,7 +34,9 @@ void SetupModel(char* DensityFile, char* NetworkFile, char* SchoolFile, char* Re if (!(Xcg1 = (long*)malloc(MAX_NUM_THREADS * CACHE_LINE_SIZE * sizeof(long)))) ERR_CRITICAL("Unable to allocate ranf storage\n"); if (!(Xcg2 = (long*)malloc(MAX_NUM_THREADS * CACHE_LINE_SIZE * sizeof(long)))) ERR_CRITICAL("Unable to allocate ranf storage\n"); - setall(P.seed1, P.seed2); + P.nextSetupSeed1 = P.setupSeed1; + P.nextSetupSeed2 = P.setupSeed2; + setall(&P.nextSetupSeed1, &P.nextSetupSeed2); P.DoBin = -1; if (P.DoHeteroDensity) @@ -280,7 +282,9 @@ void SetupModel(char* DensityFile, char* NetworkFile, char* SchoolFile, char* Re SavePeopleToPlaces(NetworkFile); //SaveDistribs(); - //setall(P.seed3,P.seed4); + // From here on, we want the same random numbers regardless of whether we used the RNG to make the network, + // or loaded the network from a file. Therefore we need to reseed the RNG. + setall(&P.nextSetupSeed1, &P.nextSetupSeed2); StratifyPlaces(); for (i = 0; i < P.NC; i++) @@ -2544,7 +2548,7 @@ void LoadPeopleToPlaces(char* NetworkFile) fread_big(&s2, sizeof(long), 1, dat); if (i != npt) ERR_CRITICAL("Number of place types does not match saved value\n"); if (j != P.N) ERR_CRITICAL("Population size does not match saved value\n"); - if ((s1 != P.seed1) || (s2 != P.seed2)) ERR_CRITICAL("Random number seeds do not match saved values\n"); + if ((s1 != P.setupSeed1) || (s2 != P.setupSeed2)) ERR_CRITICAL("Random number seeds do not match saved values\n"); k = (P.N + 999999) / 1000000; for (i = 0; i < P.N; i++) for (j = 0; j < P.PlaceTypeNum; j++) @@ -2596,8 +2600,8 @@ void SavePeopleToPlaces(char* NetworkFile) { fwrite_big(&npt, sizeof(int), 1, dat); fwrite_big(&(P.N), sizeof(int), 1, dat); - fwrite_big(&P.seed1, sizeof(long), 1, dat); - fwrite_big(&P.seed2, sizeof(long), 1, dat); + fwrite_big(&P.setupSeed1, sizeof(long), 1, dat); + fwrite_big(&P.setupSeed2, sizeof(long), 1, dat); for (i = 0; i < P.N; i++) { if ((i + 1) % 100000 == 0) fprintf(stderr, "%i saved \r", i + 1);