Merge pull request #327 from DrylandEcology/feature_pcg_seeding

* Random number generators now produce sequences that can be exactly reproduced.
* `RandSeed()` gains arguments "initstate" and "initseq" (and lost "seed") to
  fully seed a `pcg32` random number generator.
* `RandNorm()` is now re-entrant and discards one of the two generated values.
  Compilation with `"RANDNORMSTATIC"` re-produces the old, not re-entrant
  implementation.
* SOILWAT2 gains `rng_seed` as new user input (`"weathsetup.in"`).
* `SW_MKV_construct()` now only seeds `markov_rng` (the random number generator
  of the weather generator) if run as `SOILWAT2` using the new input `rng_seed`;
  `SW_MKV_construct()` does not seed `markov_rng` when run as part of `STEPWAT2`
  or `rSOILWAT2` (both of which use their own `RNG` initialization procedures).
* `SW_WTH_init_run()` now also initializes yesterday's weather values.

NEWS.md

@@ -0,0 +1,14 @@
+
+# SOILWAT2 v6.6.0
+* Random number generators now produce sequences that can be exactly reproduced.
+* `RandSeed()` gains arguments "initstate" and "initseq" (and lost "seed") to
+  fully seed a `pcg32` random number generator.
+* `RandNorm()` is now re-entrant and discards one of the two generated values.
+  Compilation with `"RANDNORMSTATIC"` re-produces the old, not re-entrant
+  implementation.
+* SOILWAT2 gains `rng_seed` as new user input (`"weathsetup.in"`).
+* `SW_MKV_construct()` now only seeds `markov_rng` (the random number generator
+  of the weather generator) if run as `SOILWAT2` using the new input `rng_seed`;
+  `SW_MKV_construct()` does not seed `markov_rng` when run as part of `STEPWAT2`
+  or `rSOILWAT2` (both of which use their own `RNG` initialization procedures).
+* `SW_WTH_init_run()` now also initializes yesterday's weather values.

SW_Markov.c

@@ -206,9 +206,14 @@ void SW_MKV_construct(void) {
 	SW_MARKOV *m = &SW_Markov;
 	size_t s = sizeof(RealD);
 
-	/* STEPWAT2: The markov_rng seed will be reset with `Globals.randseed` by
-		 its `main` at the beginning of each iteration */
-	RandSeed(0, &markov_rng);
+	/* Set seed of `markov_rng`
+	  - SOILWAT2: set seed here
+	  - STEPWAT2: `main()` uses `Globals.randseed` to (re-)set for each iteration
+	  - rSOILWAT2: R API handles RNGs
+	*/
+	#if defined(SOILWAT)
+	RandSeed(SW_Weather.rng_seed, 1u, &markov_rng);
+	#endif
 
 	m->ppt_events = 0;
 

SW_Weather.c

@@ -243,7 +243,9 @@ void SW_WTH_init_run(void) {
 	 * and with ppt=0 there's nothing to freeze.
 	 */
 	SW_Weather.now.temp_max[Today] = SW_Weather.now.temp_min[Today] = 0.;
+	SW_Weather.now.temp_max[Yesterday] = SW_Weather.now.temp_min[Yesterday] = 0.;
 	SW_Weather.now.ppt[Today] = SW_Weather.now.rain[Today] = 0.;
+	SW_Weather.now.ppt[Yesterday] = SW_Weather.now.rain[Yesterday] = 0.;
 	SW_Weather.snow = SW_Weather.snowmelt = SW_Weather.snowloss = 0.;
 	SW_Weather.snowRunoff = 0.;
 	SW_Weather.surfaceRunoff = SW_Weather.surfaceRunon = 0.;
@@ -366,7 +368,7 @@ void SW_WTH_new_day(void) {
 void SW_WTH_read(void) {
 	/* =================================================== */
 	SW_WEATHER *w = &SW_Weather;
-	const int nitems = 17;
+	const int nitems = 18;
 	FILE *f;
 	int lineno = 0, month, x;
 	RealF sppt, stmax, stmin;
@@ -398,12 +400,16 @@ void SW_WTH_read(void) {
 			break;
 
 		case 4:
+			w->rng_seed = atoi(inbuf);
+			break;
+
+		case 5:
 			x = atoi(inbuf);
 			w->yr.first = (x < 0) ? SW_Model.startyr : yearto4digit(x);
 			break;
 
 		default:
-			if (lineno == 5 + MAX_MONTHS)
+			if (lineno == 6 + MAX_MONTHS)
 				break;
 
 			x = sscanf(

SW_Weather.h

@@ -67,6 +67,7 @@ typedef struct {
 			// swTRUE: use weather generator for missing weather input (values/files)
 			// swFALSE: fail if any weather input is missing (values/files)
 		use_snow;
+	int rng_seed; // initial state for `markov_rng`
 	RealD pct_snowdrift, pct_snowRunoff;
 	SW_TIMES yr;
   RealD

doc/SOILWAT2.bib

@@ -332,3 +332,19 @@ @article{huang2018JAMC
   publisher = {American Meteorological Society},
   doi = {10.1175/JAMC-D-17-0334.1}
 }
+
+
+@article{marsaglia1964SR,
+  title = {A Convenient Method for Generating Normal Variables},
+  author = {Marsaglia, G. and Bray, T. A.},
+  year = {1964},
+  journal = {SIAM Review},
+  volume = {6},
+  number = {3},
+  pages = {260--264},
+  publisher = {Society for Industrial and Applied Mathematics},
+  issn = {0036-1445},
+  doi = {10.1137/1006063}
+}
+
+

rands.c

@@ -25,48 +25,54 @@
 /*                  Local Variables                    */
 /* --------------------------------------------------- */
 
-#ifndef RSOILWAT
-  static uint64_t stream = 1u; //stream id. this is given out to a pcg_rng then incremented.
-#endif
-
 
 /* =================================================== */
 /*             Global Function Definitions             */
 /* --------------------------------------------------- */
 
 /*****************************************************/
 /**
-  \brief Sets the random number seed.
+  \brief Set the seed of a random number generator.
 
-  \param seed The initial state of the system; if 0 then use system time.
-  \param[in,out] pcg_rng The random number generator to set.
+  The sequence produced by a `pcg` random number generator
+  (https://github.com/imneme/pcg-c-basic) is determined by two elements:
+  (i) an initial state and (ii) a sequence selector (or stream identification).
 
-  \note If using this function with STEPWAT2, then call RandSeed() only once
-    per iteration.
+  A sequence is exactly reproduced if `pcg_rng` is initialized with both
+  the same initial state and the same sequence selector.
+  Unique `initseq` values guarantee that different `pcg_rng` produce
+  non-coinciding sequences.
 
-  \sideeffect Increment the stream so that no two generators have the same
-    sequence.
+  If `initstate` is `0u`, then the state of a `pcg_rng` is initialized to
+  system time and the sequence selector to `initseq` plus system time.
+  Thus, two calls (with `initstate = 0u`) will produce different
+  random number sequences if `initseq` was different,
+  even if they occurred during the same system time.
+
+  \param initstate The initial state of the system.
+  \param initseq The initial sequence selector.
+  \param[in,out] pcg_rng The random number generator to set.
 */
-void RandSeed(signed long seed, pcg32_random_t* pcg_rng) {
-//we don't need to set a random seed if RSOILWAT is used
+void RandSeed(
+  unsigned long initstate,
+  unsigned long initseq,
+  pcg32_random_t* pcg_rng
+) {
+// R uses its own random number generators
 #ifndef RSOILWAT
 
-  if (seed == 0) {
-    //seed with a random value. Uses the system time to generate
-    //a pseudo-random seed.
-    pcg32_srandom_r(pcg_rng, time(NULL), stream);
-  }
-  else {
-    //Seed with a specific value.
-    pcg32_srandom_r(pcg_rng, (int) seed, stream);
-  }
+  if (initstate == 0u) {
+    // See `pcg32-demo.c`: seed with time
+    pcg32_srandom_r(pcg_rng, time(NULL), initseq + time(NULL));
 
-  //Increment the stream so no two generators have the same sequence.
-  stream++;
+  } else {
+    // Seed with specific values to make rng output sequence reproducible
+    pcg32_srandom_r(pcg_rng, initstate, initseq);
+  }
 
 #else
   // silence compile warnings [-Wunused-parameter]
-  if (pcg_rng == NULL && seed > 0) {}
+  if (pcg_rng == NULL && initstate > 1u && initseq > 1u) {}
 
 #endif
 }
@@ -77,6 +83,9 @@ void RandSeed(signed long seed, pcg32_random_t* pcg_rng) {
 /**
   \brief A pseudo-random number from the uniform distribution.
 
+  If `pcg_rng` was not initialized with `RandSeed()`, then the first call to
+  `RandUni()` returns 0.
+
   \param[in,out] *pcg_rng The random number generator to use.
 
   @return A pseudo-random number between 0 and 1.
@@ -282,8 +291,34 @@ void RandUniList(long count, long first, long last, RandListType list[], pcg32_r
 
 /*****************************************************/
 /**
-	 \brief A pseudo-random number from a normal distribution.
-
+	\brief A pseudo-random number from a normal distribution.
+
+	If compiled with defined `SOILWAT` or `STEPWAT`, then the function
+	implements the Marsaglia polar method
+	(Marsaglia & Bray 1964 \cite marsaglia1964SR).
+	The implementation is re-entrant
+	(but discards one of the two generated random values).
+	The original, more efficient but not re-entrant implementation is used
+	if compiled with defined `RANDNORMSTATIC`.
+
+	The version compiled with `RANDNORMSTATIC` is not re-entrant because
+	of internal static storage.
+	A first call produces two random values one of which is returned immediately;
+	the second value is internally stored in static `gset`
+	and marked by static `set`.
+	The following call will return the stored value of `gset`
+	(whether or not the call used the same `pcg_rng`) and resets `gset` and `set`.
+
+	If compiled with defined `RSOILWAT`, then `rnorm()` from header <Rmath.h>
+	is called and R handles the random number generator.
+
+	\param mean The mean of the distribution.
+	\param stddev Standard deviation of the distribution.
+	\param[in,out] *pcg_rng The random number generator to use.
+*/
+double RandNorm(double mean, double stddev, pcg32_random_t* pcg_rng) {
+/* History:
+	 cwb - 6/20/00
 	 This routine is
 	 adapted from FUNCTION GASDEV in
 	 Press, et al., 1986, Numerical Recipes,
@@ -293,32 +328,26 @@ void RandUniList(long count, long first, long last, RandListType list[], pcg32_r
 	 prior to calling any function, that
 	 depends on RandUni().
 
-	\param mean The mean of the distribution
-	\param stddev Standard deviation of the distribution
-  \param[in,out] *pcg_rng The random number generator to use.
-*/
-double RandNorm(double mean, double stddev, pcg32_random_t* pcg_rng) {
-/* History:
-	 cwb - 6/20/00
 	 cwb - 09-Dec-2002 -- FINALLY noticed that
 	 gasdev and gset have to be static!
 	 might as well set the others.
 */
 	double res;
 
 	#ifndef RSOILWAT
-		static short set = 0;
+		double v1, v2, r;
 
-		static double v1, v2, r, fac, gset, gasdev;
+		#ifdef RANDNORMSTATIC
+		/* original, non-reentrant code: issue #326 */
+		static short set = 0;
+		static double fac, gset, gasdev;
 
 		if (!set) {
-			/* `RandUni(pcg_rng)` returns 0 if `pcg_rng` is not initialized
-				 causing an infinite while-loop */
 			do {
 				v1 = 2.0 * RandUni(pcg_rng) - 1.0;
 				v2 = 2.0 * RandUni(pcg_rng) - 1.0;
 				r = v1 * v1 + v2 * v2;
-			} while (r >= 1.0);
+			} while (r >= 1.0 || r == 0.);
 			fac = sqrt(-2.0 *log(r)/r);
 			gset = v1 * fac;
 			gasdev = v2 * fac;
@@ -330,6 +359,18 @@ double RandNorm(double mean, double stddev, pcg32_random_t* pcg_rng) {
 
 		res = mean + gasdev * stddev;
 
+		#else
+		/* reentrant code: discard one of the two generated random variables */
+		do {
+			v1 = 2.0 * RandUni(pcg_rng) - 1.0;
+			v2 = 2.0 * RandUni(pcg_rng) - 1.0;
+			r = v1 * v1 + v2 * v2;
+		} while (r >= 1.0 || r == 0.);
+
+		res = mean + stddev * v2 * sqrt(-2.0 * log(r) / r);
+		#endif
+
+
 	#else
 		if (pcg_rng == NULL) {} // silence compile warnings [-Wunused-parameter]
 

rands.h

@@ -28,7 +28,11 @@ typedef long RandListType;
 /* =================================================== */
 /*             Global Function Declarations            */
 /* --------------------------------------------------- */
-void RandSeed(signed long seed, pcg32_random_t* pcg_rng);
+void RandSeed(
+  unsigned long initstate,
+  unsigned long initseq,
+  pcg32_random_t* pcg_rng
+);
 double RandUni(pcg32_random_t* pcg_rng);
 int RandUniIntRange(const long first, const long last, pcg32_random_t* pcg_rng);
 float RandUniFloatRange(const float min, const float max, pcg32_random_t* pcg_rng);

test/test_SW_Flow_Lib.cc

@@ -210,7 +210,7 @@ namespace
     double *drain2 = new double[nlyrs];
 
     pcg32_random_t infiltrate_rng;
-    RandSeed(0,&infiltrate_rng);
+    RandSeed(0u, 0u, &infiltrate_rng);
 
     for (i = 0; i < MAX_LAYERS; i++)
     {

test/test_SW_Flow_lib_temp.cc

@@ -95,7 +95,7 @@ namespace {
     unsigned int nlyrs, nRgr = 65;
     Bool ptr_stError = swFALSE;
     pcg32_random_t STInit_rng;
-    RandSeed(0,&STInit_rng);
+    RandSeed(0u, 0u, &STInit_rng);
 
     // *****  Test when nlyrs = 1  ***** //
     unsigned int i =0.;
@@ -176,7 +176,7 @@ namespace {
     double *fc2 = new double[nlyrs];
     double *wp2 = new double[nlyrs];
     pcg32_random_t STInitDeath_rng;
-    RandSeed(0,&STInitDeath_rng);
+    RandSeed(0u, 0u, &STInitDeath_rng);
 
     for (i = 0; i < nlyrs; i++) {
       bDensity2[i] = RandNorm(1.,0.5,&STInitDeath_rng);
@@ -213,7 +213,7 @@ namespace {
     Bool ptr_stError = swFALSE;
 
     pcg32_random_t SLIF_rng;
-    RandSeed(0,&SLIF_rng);
+    RandSeed(0u, 0u, &SLIF_rng);
 
     // *****  Test when nlyrs = 1  ***** //
     unsigned int i = 0.;
@@ -370,7 +370,7 @@ namespace {
     Bool ptr_stError = swFALSE;
 
     pcg32_random_t STTF_rng;
-    RandSeed(0,&STTF_rng);
+    RandSeed(0u, 0u, &STTF_rng);
 
     // declare input in for loop for non-error causing conditions;
     /// don't use RandNorm for fcR, wpR, vwcR, and bDensityR because will trigger
@@ -444,7 +444,7 @@ namespace {
     unsigned int k, i;
 
     pcg32_random_t MSTF_Lyer1_rng;
-    RandSeed(0,&MSTF_Lyer1_rng);
+    RandSeed(0u, 0u, &MSTF_Lyer1_rng);
 
     // *****  Test when nlyrs = 1  ***** //
     unsigned int nlyrs = 1, nRgr = 65;
@@ -568,7 +568,7 @@ namespace {
   TEST(SWFlowTempTest, MainSoilTemperatureFunction_LyrMAX) {
     // *****  Test when nlyrs = MAX_LAYERS  ***** //
     pcg32_random_t soilTemp_rng;
-    RandSeed(0, &soilTemp_rng);
+    RandSeed(0u, 0u, &soilTemp_rng);
 
 
     unsigned int i, k;