From f834672a466ea94ba8d0e263cd51849ab6177194 Mon Sep 17 00:00:00 2001
From: Stefano Serafin <serafin@jet02.img.univie.ac.at>
Date: Thu, 18 Jul 2024 16:24:07 +0200
Subject: [PATCH] error-growth state perturbation now follows same strategy as
initial; random number generator object is inherited from settings
---
ENDA.py | 23 +++++++++++------------
1 file changed, 11 insertions(+), 12 deletions(-)
diff --git a/ENDA.py b/ENDA.py
index 5545032..89c69be 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'):
--
GitLab