Skip to content

Commit

Permalink
Refactor setting of RNG seeds
Browse files Browse the repository at this point in the history 
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.
  • Loading branch information
@igfoo
igfoo committed Apr 23, 2020
1 parent 40da459 commit 616a814
Show file tree
Hide file tree
Showing 5 changed files with 49 additions and 57 deletions.
68 changes: 25 additions & 43 deletions src/CovidSim.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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
//// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// ****
Expand All @@ -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());
}
Expand All @@ -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)
Expand Down Expand Up @@ -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;
}

Expand Down Expand Up @@ -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);
}

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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++);
Expand Down
6 changes: 4 additions & 2 deletions src/Param.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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];
Expand Down
16 changes: 10 additions & 6 deletions src/Rand.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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)
/*
Expand Down
2 changes: 1 addition & 1 deletion src/Rand.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
14 changes: 9 additions & 5 deletions src/SetupModel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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++)
Expand Down Expand Up @@ -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++)
Expand Down Expand Up @@ -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);
Expand Down

0 comments on commit 616a814

Please sign in to comment.