Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Refactor setting of RNG seeds #121

Merged
merged 6 commits into from Apr 25, 2020
Merged
Changes from all commits
Commits
File filter...
Filter file types
Jump to…
Jump to file
Failed to load files.

Always

Just for now

@@ -32,7 +32,7 @@ CovidSim
/O:OutputFilesPrefix
[/D:PopulationDensityFile]
[/L:NetworkFileToLoad | /S:NetworkFileToSave]
Seed1 Seed2 Seed3 Seed4
SetupSeed1 SetupSeed2 RunSeed1 RunSeed2
```

Explanation of the arguments with examples:
@@ -43,7 +43,8 @@ Explanation of the arguments with examples:
- `/O:./output/NoInt_R0=1` specifies the prefix pathname for a collection of output files that contain simulation data. The output files are tabular `tsv` data (but with the extension `.xls`)
- `[/D:pop_usa_adm2.txt]` a population density file for a specific geography (e.g. a country)
- `[/L:NetworkFileToLoad | /S:NetworkFileToSave]`. For efficiency, we can run and, as a side-effect, generate a [network file](./model-glossary.md#Network\ file) that assigns [people](./model-glossary.md#People) to [places](./model-glossary.md#Places). The [network file](./model-glossary.md#Network\ file) may then be re-used for subsequent runs (with different input parameters for the same geography). The network file is a non-portable `.bin`. Generate this file with the `/S` option and re-use it (in a subsequent run) with the `/L` option.
- `Seed1 Seed2 Seed3 Seed4` Random seeds (larger positive integers).
- `SetupSeed1 SetupSeed2` Random number generator seeds used when initialising the model, including creating the networkfile (large positive integers).
- `RunSeed1 RunSeed2` Random number generator seeds used when running the model. These can be varied to do multiple runs with the same network file (large positive integers).

## Additional command-line arguments

@@ -62,7 +63,7 @@ CovidSim
[/AP:AirTravelFile]
[/s:SchoolFile]
[/R:R0scaling]
Seed1 Seed2 Seed3 Seed4
SetupSeed1 SetupSeed2 RunSeed1 RunSeed2
```
Explanation of additional arguments:
- `[/AP:AirTravelFile]` Air travel data for a specific geography (unused currently)
@@ -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;
@@ -288,7 +288,7 @@ int main(int argc, char* argv[])
sprintf(OutFile, "%s", OutFileBase);

fprintf(stderr, "Param=%s\nOut=%s\nDens=%s\n", ParamFile, OutFile, DensityFile);
if (Perr) ERR_CRITICAL_FMT("Syntax:\n%s /P:ParamFile /O:OutputFile [/AP:AirTravelFile] [/s:SchoolFile] [/D:DensityFile] [/L:NetworkFileToLoad | /S:NetworkFileToSave] [/R:R0scaling] Seed1 Seed2 Seed3 Seed4\n", argv[0]);
if (Perr) ERR_CRITICAL_FMT("Syntax:\n%s /P:ParamFile /O:OutputFile [/AP:AirTravelFile] [/s:SchoolFile] [/D:DensityFile] [/L:NetworkFileToLoad | /S:NetworkFileToSave] [/R:R0scaling] SetupSeed1 SetupSeed2 RunSeed1 RunSeed2\n", argv[0]);

//// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// **** //// ****
//// **** SET UP OMP / THREADS
@@ -327,7 +327,7 @@ int main(int argc, char* argv[])
if (GotScF) P.DoSchoolFile = 1;
if (P.DoAirports)
{
if (!GotAP) ERR_CRITICAL_FMT("Syntax:\n%s /P:ParamFile /O:OutputFile /AP:AirTravelFile [/s:SchoolFile] [/D:DensityFile] [/L:NetworkFileToLoad | /S:NetworkFileToSave] [/R:R0scaling] Seed1 Seed2 Seed3 Seed4\n", argv[0]);
if (!GotAP) ERR_CRITICAL_FMT("Syntax:\n%s /P:ParamFile /O:OutputFile /AP:AirTravelFile [/s:SchoolFile] [/D:DensityFile] [/L:NetworkFileToLoad | /S:NetworkFileToSave] [/R:R0scaling] SetupSeed1 SetupSeed2 RunSeed1 RunSeed2\n", argv[0]);
ReadAirTravel(AirTravelFile);
}

@@ -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)
@@ -2931,9 +2915,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;
}

@@ -3937,7 +3919,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);
}

@@ -3996,8 +3978,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);
@@ -4070,9 +4052,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++);
@@ -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];
@@ -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)
/*
@@ -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);
@@ -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)
@@ -284,7 +286,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++)
@@ -2560,7 +2564,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++)
@@ -2612,8 +2616,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);
@@ -1,7 +1,10 @@
CI_100_91_R0=2.2.00100.bmp f13526a03252773523b04195f6edcfabc1a4b311bd6f8e06e30c837c22128864414eb308754c3faa9ec3b4dea9601a462bedba1a96a83f43394db97e9719fc0a
CI_100_91_R0=2.2.avNE.age.xls 0ed19510c9ffc036e34100f25cd106a055a37ed7b9da565d67c82a416638e5d21dc051c81e29e444dacd2c407ead453ab104f995e5a937221d599c214eac95b8
CI_100_91_R0=2.2.avNE.severity.adunit.xls fac454e05f5f42d9cb5263aba6ce4ed934be49a9a3744a010403ca8ac53d86fa30de08992989567d4d3048810bf11f9f285b4d95ebd3e51fd1d60a874a0c1c9b
CI_100_91_R0=2.2.avNE.severity.xls 74b8850021cbde47c015f3b117735caa7e85a13a7f87ff951d6da1210283c1d3fab67eafad0f1a8ccefd50e1969769261e0ae90fd675fc53102aa119ce80305b
NoInt_R0=2.2.avNE.age.xls 9859e60445c690af59dec5c4a74259c399ba80744e5e6bc1e5beba95c2c15478d311825adf8bc7e2f7b401a88a2620e1daef7f3ebefab46191feff180bae0ece
NoInt_R0=2.2.avNE.severity.adunit.xls db2b23dfc70480940cf2d1528d182549b77834ac0bbdc2b503524a18cf1a9401c2545293478de9fb10f0e9abd07c5862f481fe6493ab9a2a8c384024749409e9
NoInt_R0=2.2.avNE.severity.xls bd1fbc152f0524dbccc9e700fb8ac911f5f542594439bafa89fb773f39666168449013639e2a659c5f85225fd2a424acd682988cfe84de7492ef8339eaaa07ee
CI_100_91_R0=2.2.00100.bmp 253bb53d534929fc5f7ed3e011ba0a3baf2b6b19ff567169ef39be3e8d3625ab1e8aaf31d4c5f659132c296cba7687f7706d26206b0a7be4d834e2a639e3f755
CI_100_91_R0=2.2.avNE.age.xls 604bb4fb66468cda445d26a0c1fad563b9edabe4ca7b50940d612c04b06b639f79ae7e79542c1ee5da956c61ceddc3ff3ac418458f9487fbed12bfb04425909f
CI_100_91_R0=2.2.avNE.severity.adunit.xls 79f1b6adc0ba9bbc050fb13630f7887ced692f3a505fff3b5e13e3bc9363d90037792373b3a7dde2fe1f85a0c508a73203ae320e5d7016aeb824fc023802474c
CI_100_91_R0=2.2.avNE.severity.xls 212f43ef160a640da374b8942624468b35d39471c2d8356221adc8360a481bb0f437068aeb91f0ace2a428114aa2fed3f978f8a15922e69641c81c15cdbb04e6
NoInt_R0=2.2-repeat.avNE.age.xls ae837285aafa50c46395fde69de3fa0f423242c25f9356b947d54ba42068d6401c985e1aea72a6d60cece8321154300cf537e16d76dbcf004c268b3320ca208d
NoInt_R0=2.2-repeat.avNE.severity.adunit.xls 872c7c19e5be46636618c9d975566c02f7172d812e054ee01d295f1a68d2b15a4e87417564d0cbc1fffe830fbf5b51ad2ffc262927569f01e5a6d18313db46ee
NoInt_R0=2.2-repeat.avNE.severity.xls 546ec5241e7cff3a12be90b66bcdda38e3ff34cbf79bcf441bcecb1f21e69a93378b90a46363c8c5a558ca0c9e9efe8cf9b9d599039dfed1ea3ed655be825479
NoInt_R0=2.2.avNE.age.xls ae837285aafa50c46395fde69de3fa0f423242c25f9356b947d54ba42068d6401c985e1aea72a6d60cece8321154300cf537e16d76dbcf004c268b3320ca208d
NoInt_R0=2.2.avNE.severity.adunit.xls 872c7c19e5be46636618c9d975566c02f7172d812e054ee01d295f1a68d2b15a4e87417564d0cbc1fffe830fbf5b51ad2ffc262927569f01e5a6d18313db46ee
NoInt_R0=2.2.avNE.severity.xls 546ec5241e7cff3a12be90b66bcdda38e3ff34cbf79bcf441bcecb1f21e69a93378b90a46363c8c5a558ca0c9e9efe8cf9b9d599039dfed1ea3ed655be825479
@@ -68,6 +68,7 @@
# that you can compare the failing versions of the .xls files against
# good versions.

import glob
import gzip
import os
import sys
@@ -76,6 +77,7 @@

testdir='regressiontest_UK_100th'

failed = False
accept_results = False
if len(sys.argv) > 1 and sys.argv[1] == '--accept':
accept_results = True
@@ -120,6 +122,16 @@
'/S:NetworkUKN_32T_100th.bin', '/R:1.1',
'98798150', '729101', '17389101', '4797132'
])
print('=== Starting running repeat:')
subprocess.check_call(
[covidsim_exe, '/c:1',
'/PP:' + updir1 + 'preUK_R0=2.0.txt',
'/P:' + updir1 + 'p_NoInt.txt', '/CLP1:100000',
'/CLP2:0', '/O:NoInt_R0=2.2-repeat', '/D:' + wpop_bin,
'/A:' + updir1 + 'sample_admin.txt',
'/L:NetworkUKN_32T_100th.bin', '/R:1.1',
'98798150', '729101', '17389101', '4797132'
])
print('=== Starting running:')
subprocess.check_call(
[covidsim_exe, '/c:1',
@@ -133,6 +145,23 @@
])
print('=== Done')

repeat_files_checked = 0
for fn2 in glob.glob('NoInt_R0=2.2-repeat*.xls'):
fn1 = fn2.replace('-repeat', '')
with open(fn1, 'rb') as f:
dat1 = f.read()
with open(fn2, 'rb') as f:
dat2 = f.read()
# We expect multiple reasonably large files, so only count those that are at least 10 bytes long
if len(dat1) > 10:
repeat_files_checked += 1
if dat1 != dat2:
print('FAILURE: Contents of ' + fn1 + ' does not match that of ' + fn2)
failed = True
if repeat_files_checked < 3:
print('FAILURE: Not enough repeat files found')
failed = True

checksums_filename='regressiontest_UK_100th.checksums.txt'

# Compute SHA-512 checksums for the generated files.
@@ -200,4 +229,8 @@
with open(os.pardir + os.sep + checksums_filename, 'wb') as checksums_out:
shutil.copyfileobj(checksums_in, checksums_out)
else:
sys.exit(1)
failed = True

if failed:
print('TEST FAILED.')
sys.exit(1)
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.