diff --git a/ENDA.py b/ENDA.py index 5545032abd366f29a47468a48b80ce38a1cf5c8c..89c69be23a61e010e7cb480f149027faa369cdd2 100644 --- a/ENDA.py +++ b/ENDA.py @@ -20,9 +20,6 @@ from models import CBL verbose = True -# Initiate a random number generator with preset seed, to ensure reproducibility -RNG = np.random.default_rng(seed=240718) - def transform_none(x,kind): if kind=='dir': return x @@ -245,7 +242,8 @@ class cycle: obs_error_sdev_assimilate,\ localization_cutoff,\ inflate_rtps,\ - inflation_rtps_alpha): + inflation_rtps_alpha,\ + RNG): # How many cycles do you need to run? ncycles = int(trun/assimilation_interval) @@ -311,6 +309,12 @@ class cycle: # Store prior backgrounds[k,:,:] = da.x[:,:] + # If you want to add/simulate some error growth in the run, before the analysis, + # perturb the prior state here. + if isinstance(da,CBL): + if cbl_settings["simulate_error_growth"]: + backgrounds[k,:,:] = CBL.perturb_state(da,backgrounds[k,:,:]) + # Compute prior spread (necessary for variance inflation) if inflate_rtps: sigp = np.std(backgrounds[k,:,:],axis=1) @@ -348,12 +352,6 @@ class cycle: # Store posterior analyses[k,:,:] = x0 - # If you want to add/simulate some error growth in the next cycle, - # perturb the state here. - if isinstance(nr,CBL): - if cbl_settings["simulate_error_growth"]: - x0[:-cbl_settings["parameter_number"],:] += RNG.normal(scale=cbl_settings["error_growth_perturbations_amplitude"],size=nens)[None,:] - # Save initial conditions for next cycle da.update(x0) @@ -445,7 +443,7 @@ class experiment: variable = nr.history['v'][:,time_index] truths[i,j] = observation_operator(variable,state_coordinates,obs_coordinate) observations[i,j] = truths[i,j]+\ - RNG.normal(0,self.obs_error_sdev_generate[j]) + self.RNG.normal(0,self.obs_error_sdev_generate[j]) # Store truths and observations self.truths = truths @@ -468,7 +466,8 @@ class experiment: self.obs_error_sdev_assimilate, self.localization_cutoff, self.inflate_rtps, - self.inflation_rtps_alpha) + self.inflation_rtps_alpha, + self.RNG) # Compute diagnostics if not hasattr(self,'dg'):