v524 - EnKF random number generator bug fix

bug fix to GaussRandom() in ModelEnsemble.cpp - had one in 33,000 chance of throwing a -inf

Author: James Craig

src/EnKF.cpp

@@ -188,6 +188,10 @@ EnKF_mode CEnKFEnsemble::GetEnKFMode() const
 
 //////////////////////////////////////////////////////////////////
 /// \brief initializes EnKF assimilation run
+/// - performs QA/QC on inputs 
+/// - gets #s and names of state variables adjusted during assimilation
+/// - allocate & populate _output_matrix, _obs_matrix, _noise_matrix
+/// - open blank EnKFOutput.csv file 
 /// \param &Options [out] Global model options information
 //
 void CEnKFEnsemble::Initialize(const CModel* pModel,const optStruct &Options)
@@ -197,8 +201,8 @@ void CEnKFEnsemble::Initialize(const CModel* pModel,const optStruct &Options)
 
   // QA/QC
   //-----------------------------------------------
-  if(_nMembers==0) {
-    ExitGracefully("CEnKFEnsemble::Initialize: number of EnKF ensemble members must be >0",BAD_DATA);
+  if(_nMembers<=1) {
+    ExitGracefully("CEnKFEnsemble::Initialize: number of EnKF ensemble members must be >1",BAD_DATA);
   }
   if ((pModel->GetNumForcingPerturbations() == 0) && (_EnKF_mode != ENKF_FORECAST) && (_EnKF_mode != ENKF_OPEN_FORECAST)) {
     ExitGracefully("CEnKFEnsemble::Initialize: at least one forcing perturbation must be set using :ForcingPerturbation command",BAD_DATA);
@@ -358,7 +362,7 @@ void CEnKFEnsemble::Initialize(const CModel* pModel,const optStruct &Options)
     ExitGracefullyIf(_noise_matrix[e]==NULL,"EnKF Initialization: noise matrix",OUT_OF_MEMORY);
   }
 
-  //populate _output_matrix, _obs_matrix
+  //populate _output_matrix, _obs_matrix, _noise_matrix
   //-----------------------------------------------
   int j;
   double eps;
@@ -427,7 +431,7 @@ void CEnKFEnsemble::StartTimeStepOps(CModel* pModel,optStruct& Options,const tim
   }
 }
 //////////////////////////////////////////////////////////////////
-/// \brief called at end of each time step
+/// \brief called at end of each time step - updates state matrix and builds output matrix
 /// \param pModel [out] pointer to global model instance
 /// \param &Options [out] Global model options information
 /// \params tt [in] time structure
@@ -439,7 +443,7 @@ void CEnKFEnsemble::CloseTimeStepOps(CModel* pModel,optStruct& Options,const tim
     AddToStateMatrix(pModel,Options,e);
   }
 
-  //Build output Matrix
+  //Build output Matrix (typically final timestep of simulation)
   //----------------------------------------------------------
   int curr_nn=(int)((tt.model_time+TIME_CORRECTION)/Options.timestep);//current timestep index
 
@@ -680,7 +684,7 @@ void CEnKFEnsemble::AddToStateMatrix(CModel* pModel,optStruct& Options,const int
   for(int i=0;i<_nAssimStates;i++)
   {
     kk=_aAssimGroupID[i];
-    if(_aAssimStates[i]==STREAMFLOW)
+    if     (_aAssimStates[i]==STREAMFLOW)
     {
       //cout<<"RETRIEVE STATES : # subbasins: "<<pModel->GetSubBasinGroup(kk)->GetNumSubbasins()<<endl;
       for(int pp=0;pp<pModel->GetSubBasinGroup(kk)->GetNumSubbasins();pp++)
@@ -808,62 +812,66 @@ void CEnKFEnsemble::FinishEnsembleRun(CModel *pModel,optStruct &Options,const ti
 
   if(e==_nEnKFMembers-1) //After all ensemble members have run
   {
-    _ENKFOUT<<"PRE-ASSIMILATION STATE MATRIX:"<<endl;
-    _ENKFOUT<<"member"<<",";
-    for (int i = 0; i<_nStateVars; i++) {
-      _ENKFOUT<<_state_names[i]<<",";
-    }
-    _ENKFOUT<<endl;
-    for(int ee=0;ee<_nEnKFMembers;ee++) {
-      _ENKFOUT<<ee+1<<",";
-      for(int i=0;i<_nStateVars;i++) {
-        _ENKFOUT<<_state_matrix[ee][i]<<",";
+    {
+      _ENKFOUT<<"PRE-ASSIMILATION STATE MATRIX:"<<endl;
+      _ENKFOUT<<"member"<<",";
+      for (int i = 0; i<_nStateVars; i++) {
+        _ENKFOUT<<_state_names[i]<<",";
       }
       _ENKFOUT<<endl;
+      for(int ee=0;ee<_nEnKFMembers;ee++) {
+        _ENKFOUT<<ee+1<<",";
+        for(int i=0;i<_nStateVars;i++) {
+          _ENKFOUT<<_state_matrix[ee][i]<<",";
+        }
+        _ENKFOUT<<endl;
+      }
     }
 
     cout<<endl<<"ENKF: Performing Assimilation Calculations with "<<_nObsDatapoints<<" observation datapoints..."<<endl;
     AssimilationCalcs();
 
-    _ENKFOUT<<"POST-ASSIMILATION STATE MATRIX:"<<endl;
-    _ENKFOUT<<"member"<<",";
-    for (int i = 0; i<_nStateVars; i++) {
-      _ENKFOUT<<_state_names[i]<<",";
-    }
-    _ENKFOUT<<endl;
-    for(int ee=0;ee<_nEnKFMembers;ee++) {
-      _ENKFOUT<<ee+1<<",";
-      for(int i=0;i<_nStateVars;i++) {
-        _ENKFOUT<<_state_matrix[ee][i]<<",";
+    {
+      _ENKFOUT<<"POST-ASSIMILATION STATE MATRIX:"<<endl;
+      _ENKFOUT<<"member"<<",";
+      for (int i = 0; i<_nStateVars; i++) {
+        _ENKFOUT<<_state_names[i]<<",";
       }
       _ENKFOUT<<endl;
-    }
+      for(int ee=0;ee<_nEnKFMembers;ee++) {
+        _ENKFOUT<<ee+1<<",";
+        for(int i=0;i<_nStateVars;i++) {
+          _ENKFOUT<<_state_matrix[ee][i]<<",";
+        }
+        _ENKFOUT<<endl;
+      }
 
-    _ENKFOUT<<"NOISE MATRIX:"<<endl;
-    for(int ee=0;ee<_nEnKFMembers;ee++) {
-      _ENKFOUT<<ee+1<<",";
-      for(int i=0;i<_nObsDatapoints;i++) {
-        _ENKFOUT<<_noise_matrix[ee][i]<<",";
+      _ENKFOUT<<"NOISE MATRIX:"<<endl;
+      for(int ee=0;ee<_nEnKFMembers;ee++) {
+        _ENKFOUT<<ee+1<<",";
+        for(int i=0;i<_nObsDatapoints;i++) {
+          _ENKFOUT<<_noise_matrix[ee][i]<<",";
+        }
+        _ENKFOUT<<endl;
       }
-      _ENKFOUT<<endl;
-    }
-    _ENKFOUT<<"OBSERVATION MATRIX:"<<endl;
-    for(int ee=0;ee<_nEnKFMembers;ee++) {
-      _ENKFOUT<<ee+1<<",";
-      for(int i=0;i<_nObsDatapoints;i++) {
-        _ENKFOUT<<_obs_matrix[ee][i]<<",";
+      _ENKFOUT<<"OBSERVATION MATRIX:"<<endl;
+      for(int ee=0;ee<_nEnKFMembers;ee++) {
+        _ENKFOUT<<ee+1<<",";
+        for(int i=0;i<_nObsDatapoints;i++) {
+          _ENKFOUT<<_obs_matrix[ee][i]<<",";
+        }
+        _ENKFOUT<<endl;
       }
-      _ENKFOUT<<endl;
-    }
-    _ENKFOUT<<"SIMULATED OUTPUT MATRIX:"<<endl;
-    for(int ee=0;ee<_nEnKFMembers;ee++) {
-      _ENKFOUT<<ee+1<<",";
-      for(int i=0;i<_nObsDatapoints;i++) {
-        _ENKFOUT<<_output_matrix[ee][i]<<",";
+      _ENKFOUT<<"SIMULATED OUTPUT MATRIX:"<<endl;
+      for(int ee=0;ee<_nEnKFMembers;ee++) {
+        _ENKFOUT<<ee+1<<",";
+        for(int i=0;i<_nObsDatapoints;i++) {
+          _ENKFOUT<<_output_matrix[ee][i]<<",";
+        }
+        _ENKFOUT<<endl;
       }
-      _ENKFOUT<<endl;
+      _ENKFOUT.close();
     }
-    _ENKFOUT.close();
 
     //Write EnKF-updated solution files
     // requires full re-reading of ensemble member solution files

src/Model.cpp

@@ -2726,15 +2726,17 @@ void CModel::RecalculateHRUDerivedParams(const optStruct    &Options,
 //
 void CModel::PrepareForcingPerturbation(const optStruct &Options, const time_struct &tt)
 {
+  if (_nPerturbations==0){return;}
+  
+  double partday      = Options.julian_start_day-floor(Options.julian_start_day+TIME_CORRECTION);
+  int    nn           = (int)(rvn_round((tt.model_time+partday-floor(tt.model_time+partday+TIME_CORRECTION))/Options.timestep));
+  bool   start_of_day = ((nn==0) || tt.day_changed); //nn==0 corresponds to midnight
 
-  for(int i=0;i<_nPerturbations;i++)
-  {
+  if (start_of_day)  { //get all random perturbation samples for the day
     int    nStepsPerDay = (int)(rvn_round(1.0/Options.timestep));
-    double partday      = Options.julian_start_day-floor(Options.julian_start_day+TIME_CORRECTION);
-    int    nn           = (int)(rvn_round((tt.model_time+partday-floor(tt.model_time+partday+TIME_CORRECTION))/Options.timestep));
-    bool   start_of_day = ((nn==0) || tt.day_changed); //nn==0 corresponds to midnight
 
-    if (start_of_day)  { //get all random perturbation samples for the day
+    for(int i=0;i<_nPerturbations;i++)
+    {
       for (int n=0;n<nStepsPerDay;n++){
         _pPerturbations[i]->eps[n]=SampleFromDistribution(_pPerturbations[i]->distribution,_pPerturbations[i]->distpar);
       }

src/ModelEnsemble.cpp

@@ -23,7 +23,7 @@ double UniformRandom()
 //
 double GaussRandom()
 {
-  double u1=(double)(rand())/RAND_MAX;
+  double u1=(double)(rand()+1)/(RAND_MAX+1.0); //avoid zero
   double u2=(double)(rand())/RAND_MAX;
   return sqrt(-2.0*log(u1))*cos(2.0*PI*u2);
 }

src/ParseInput.cpp

@@ -3619,6 +3619,8 @@ bool ParseMainInputFile (CModel     *&pModel,
                    "ParseMainInputFile::Must have a postitive time step",BAD_DATA);
   ExitGracefullyIf(Options.duration<0,
                    "ParseMainInputFile::Model duration less than zero. Make sure :EndDate is after :StartDate.",BAD_DATA_WARN);
+  ExitGracefullyIf(Options.julian_start_year==1666,
+                   "ParseMainInputFile:: no :StartDate supplied in .rvi or runinfo.nc file.",BAD_DATA_WARN);
 
   if((Options.nNetCDFattribs>0) && (Options.output_format!=OUTPUT_NETCDF)){
     WriteAdvisory("ParseMainInputFile: NetCDF attributes were specified but output format is not NetCDF.",Options.noisy);

src/RavenMain.cpp

@@ -189,7 +189,7 @@ int main(int argc, char* argv[])
 }
 
 //////////////////////////////////////////////////////////////////
-/// \param argc [in] number of arguments to executable
+/// \param argc [in] number of command line arguments to executable
 /// \param argv[] [in] executable arguments; Raven.exe [filebase] [-p rvp_file] [-h hru_file] [-t rvt_file] [-c rvc_file] [-o output_dir]
 /// \details initializes input files and output directory
 /// \details filebase has no extension, all others require .rv* extension

src/UnitTesting.cpp

@@ -36,6 +36,7 @@ void RavenUnitTesting(const optStruct &Options)
   //TestWetBulbTemps();
   //TestDateStrings();
   //IncGammaLimit();
+  //TestRandomSampling();
 }
 /////////////////////////////////////////////////////////////////
 /// \brief Tests DateStringToTimeStruct() function
@@ -614,4 +615,13 @@ void IncGammaLimit() {
     cout<<x<<" "<<IncompleteGamma(xx,1.0)<<" "<<IncompleteGamma(xx,2.0)<<" "<<IncompleteGamma(xx,3.0)<<" "<<IncompleteGamma(xx,4.0)<<" "<<IncompleteGamma(xx,5.0)<<endl;
   }
   ExitGracefully("UnitTesting:: IncGammaLimit",SIMULATION_DONE);
-}
+}
+double GaussRandom();
+void TestRandomSampling() {
+  ofstream TEST;
+  TEST.open("GaussSample.csv");
+  for(int m=0; m<10000;m++) {
+    TEST<<GaussRandom()<<endl;
+  }
+  ExitGracefully("UnitTesting:: TestRandomSampling",SIMULATION_DONE);
+}

src/UnitTesting.h

@@ -26,4 +26,5 @@ void TestGammaSampling();
 void TestWetBulbTemps();
 void TestDateStrings();
 void IncGammaLimit();
+void TestRandomSampling();
 #endif