Skip to content
Snippets Groups Projects
Commit f8559de7 authored by Stefano Serafin's avatar Stefano Serafin
Browse files

for reproducibility, replaced all np.rand occurrences with a generator object, with preset seed

parent c7fd8e17
Branches
Tags
No related merge requests found
......@@ -20,6 +20,9 @@ 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
......@@ -349,7 +352,8 @@ class cycle:
# perturb the state here.
if isinstance(nr,CBL):
if cbl_settings["simulate_error_growth"]:
x0[:-cbl_settings["parameter_number"],:] += np.random.normal(scale=cbl_settings["error_growth_perturbations_amplitude"],size=nens)[None,:]
x0[:-cbl_settings["parameter_number"],:] += RNG.normal(scale=cbl_settings["error_growth_perturbations_amplitude"],size=nens)[None,:]
#x0[:-cbl_settings["parameter_number"],:] += np.random.normal(scale=cbl_settings["error_growth_perturbations_amplitude"],size=nens)[None,:]
# Save initial conditions for next cycle
da.update(x0)
......@@ -442,7 +446,8 @@ class experiment:
variable = nr.history['v'][:,time_index]
truths[i,j] = observation_operator(variable,state_coordinates,obs_coordinate)
observations[i,j] = truths[i,j]+\
np.random.normal(0,self.obs_error_sdev_generate[j])
RNG.normal(0,self.obs_error_sdev_generate[j])
#np.random.normal(0,self.obs_error_sdev_generate[j])
# Store truths and observations
self.truths = truths
......
......@@ -5,6 +5,9 @@ from scipy.interpolate import CubicSpline
verbose = True
# Initiate a random number generator with preset seed, to ensure reproducibility
RNG = np.random.default_rng(seed=240718)
class CBL:
def __init__(self,settings):
for k, v in settings.items():
......@@ -98,10 +101,13 @@ class CBL:
randomsize = 1
if self.perturbations_type == "random" or self.perturbations_type == "uniform":
ppt = np.random.randn(randomsize,self.nens)
ppt = RNG.standard_normal((randomsize,self.nens))
#ppt = np.random.randn(randomsize,self.nens)
if self.is_bwind:
ppu = np.random.randn(self.nens,randomsize)
ppv = np.random.randn(self.nens,randomsize)
ppu = RNG.standard_normal((self.nens,randomsize))
ppv = RNG.standard_normal((self.nens,randomsize))
#ppu = np.random.randn(self.nens,randomsize)
#ppv = np.random.randn(self.nens,randomsize)
# Smooth perturbations are slightly more complicated
if self.perturbations_type == "smooth":
......@@ -114,15 +120,18 @@ class CBL:
# Draw random perturbations, then interpolate
randomsize = self.perturbations_smooth_number
ppt = np.zeros((self.nz,self.nens))+np.nan
pert_t = np.random.randn(randomsize,self.nens)
pert_t = RNG.standard_normal((randomsize,self.nens))
#pert_t = np.random.randn(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 = np.random.randn(randomsize,self.nens)
pert_v = np.random.randn(randomsize,self.nens)
pert_u = RNG.standard_normal((randomsize,self.nens))
pert_v = RNG.standard_normal((randomsize,self.nens))
#pert_u = np.random.randn(randomsize,self.nens)
#pert_v = np.random.randn(randomsize,self.nens)
for n in range(self.nens):
f = CubicSpline(ipert,pert_u[:,n])
ppu[:,n] = f(np.arange(self.nz))
......@@ -212,7 +221,8 @@ class CBL:
pp = np.zeros((self.parameter_number,self.nens))
for k in range(-self.parameter_number,0):
dum = np.random.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
dum = RNG.uniform(self.parameter_ensemble_min[k], self.parameter_ensemble_max[k], size=self.nens)
#dum = np.random.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
......@@ -307,7 +317,8 @@ class CBL:
H0 = Hmax
# Add random perturbations to the initial value
H0 += np.random.normal(scale=H0_perturbation_ampl_init)
H0 += RNG.normal(scale=H0_perturbation_ampl_init)
#H0 += np.random.normal(scale=H0_perturbation_ampl_init)
# Set the surface momentum flux (ustar)
ustar = 0
......@@ -357,7 +368,8 @@ class CBL:
# Add time-dependent surface perturbations
# Then compute sensible heat flux and integrate T equation
H[0,j] = H0 + np.random.normal(scale=H0_perturbation_ampl_time)
H[0,j] = H0 + RNG.normal(scale=H0_perturbation_ampl_time)
#H[0,j] = H0 + np.random.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])
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment