fixed generation of random numbers, to really ensure reproducibility; minor fixes to plots

ENDA.py

@@ -281,7 +281,7 @@ class cycle:
                  localization_cutoff,\
                  inflate_rtps,\
                  inflation_rtps_alpha,\
-                 RNG):
+                 rnseed):
 
         # How many cycles do you need to run?
         ncycles = int(trun/assimilation_interval)
@@ -515,8 +515,13 @@ class experiment:
                     elif self.obs_kinds[j]=='v':
                         variable = nr.history['v'][:,time_index]
                     truths[i,j] = observation_operator(variable,state_coordinates,obs_coordinate)
+                    # Change the seed every time you go through this,
+                    # otherwise the perturbations are always the same;
+                    # and you have bias instead of randomness
+                    seed = self.rnseed+j*1000+i*100000
+                    RNG = np.random.default_rng(seed=seed)
                     observations[i,j] = truths[i,j]+\
-                        self.RNG.normal(0,self.obs_error_sdev_generate[j])
+                        RNG.normal(0,self.obs_error_sdev_generate[j])
 
             # Store truths and observations
             self.truths = truths
@@ -540,7 +545,7 @@ class experiment:
                         self.localization_cutoff,
                         self.inflate_rtps,
                         self.inflation_rtps_alpha,
-                        self.RNG)
+                        self.rnseed)
     
         # Compute diagnostics
         if not hasattr(self,'dg'):

graphics.py

@@ -421,7 +421,7 @@ def plot_diagnostics(experiments_pe,experiments_nope,labels,filename):
     fig.set_size_inches(6,4)
     z = experiments_pe[0].obs_coordinates
     z_pbl = z*1.
-    z_pbl[z>1000] = np.nan
+    z_pbl[z>1500] = np.nan
     for i in range(len(experiments_pe)):
         i1 = experiments_pe[i].dg
         i2 = experiments_nope[i].dg

models.py

@@ -97,7 +97,8 @@ class CBL:
         # Overwrite the pre-existing unperturbed values
         # This is needed for safety, because parameter transformations may be nonlinear
         if self.nens>1:
-            if self.do_parameter_estimation and hasattr(self,"initial_perturbed_parameters"):
+            #if self.do_parameter_estimation and hasattr(self,"initial_perturbed_parameters"):
+            if hasattr(self,"initial_perturbed_parameters"):
                 x0[-self.parameter_number:,:] = self.initial_perturbed_parameters
             else:
                 if self.perturb_ensemble_parameters or self.do_parameter_estimation:
@@ -158,7 +159,8 @@ class CBL:
         pp = np.zeros((self.parameter_number,self.nens))
 
         for k in range(-self.parameter_number,0):
-            dum = self.RNG.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
+            RNG = np.random.default_rng(seed=self.rnseed)
+            dum = RNG.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
             pp[k,:] = self.parameter_transform[k](dum , kind='dir')
 
         return pp
@@ -172,10 +174,12 @@ class CBL:
             randomsize = 1
         
         if self.perturbations_type == "random" or self.perturbations_type == "uniform":
-            ppt = self.RNG.standard_normal((randomsize,self.nens))
+            RNG = np.random.default_rng(seed=self.rnseed)
+            ppt = RNG.standard_normal((randomsize,self.nens))
             if self.is_bwind:
-                ppu = self.RNG.standard_normal((self.nens,randomsize))
-                ppv = self.RNG.standard_normal((self.nens,randomsize))
+                RNG = np.random.default_rng(seed=self.rnseed)
+                ppu = RNG.standard_normal((self.nens,randomsize))
+                ppv = RNG.standard_normal((self.nens,randomsize))
 
         # Smooth perturbations are slightly more complicated
         if self.perturbations_type == "smooth":
@@ -188,15 +192,17 @@ class CBL:
             # Draw random perturbations, then interpolate
             randomsize = self.perturbations_smooth_number
             ppt = np.zeros((self.nz,self.nens))+np.nan
-            pert_t = self.RNG.standard_normal((randomsize,self.nens))
+            RNG = np.random.default_rng(seed=self.rnseed)
+            pert_t = RNG.standard_normal((randomsize,self.nens))
             for n in range(self.nens):
                 f = CubicSpline(ipert,pert_t[:,n])
                 ppt[:,n] = f(np.arange(self.nz))
             if self.is_bwind:
                 ppu = np.zeros((self.nz,self.nens))+np.nan
                 ppv = np.zeros((self.nz,self.nens))+np.nan
-                pert_u = self.RNG.standard_normal((randomsize,self.nens))
-                pert_v = self.RNG.standard_normal((randomsize,self.nens))
+                RNG = np.random.default_rng(seed=self.rnseed)
+                pert_u = RNG.standard_normal((randomsize,self.nens))
+                pert_v = RNG.standard_normal((randomsize,self.nens))
                 for n in range(self.nens):
                     f = CubicSpline(ipert,pert_u[:,n])
                     ppu[:,n] = f(np.arange(self.nz))
@@ -311,7 +317,11 @@ class CBL:
                 H0 = Hmax
             
             # Add random perturbations to the initial value
-            H0 += self.RNG.normal(scale=H0_perturbation_ampl_init)
+            # Change the seed every time you go through this,
+            # otherwise the perturbations are always the same
+            seed = self.rnseed+j*200000
+            RNG = np.random.default_rng(seed=seed)
+            H0 += RNG.normal(scale=H0_perturbation_ampl_init)
 
             # Set the surface momentum flux (ustar)
             ustar = 0
@@ -360,8 +370,13 @@ class CBL:
             rdz = 1./self.dz
 
             # Add time-dependent surface perturbations
+            # Change the seed every time you go through this,
+            # otherwise the perturbations are always the same
+            seed = self.rnseed+j*200000
+            RNG = np.random.default_rng(seed=seed)
+
             # Then compute sensible heat flux and integrate T equation
-            H[0,j]     = H0 + self.RNG.normal(scale=H0_perturbation_ampl_time)
+            H[0,j]     = H0 + RNG.normal(scale=H0_perturbation_ampl_time)
             H[1:-1,j]  = -K[1:-1,j]*( (thetap[1:]-thetap[:-1])*rdz - gammac)
             H[nz,j]    = 2*H[nz-1,j]-H[nz-2,j]
             theta[:,j] = thetap[:]-dt*rdz*(H[1:,j]-H[:-1,j])