Bias estimators

Summary

File: src/bias_estimators/bias.py

Base class for observation-bias estimators. Mixes in HistoryTracker. Wraps a forecaster model to produce bias corrections at each assimilation step.

Key attributes: innovation, dt, forecaster, Nq, N_dim, upsample, biased_observations

Subclasses must implement: - init_forecaster(**kwargs) — build/train the internal forecasting model - state_derivative() — return the Jacobian J = db/dy used by bias-aware filters

Bias
├── ESN_bias
└── ConstantBias
    └── NoBias

Class	File	Notes
ESN_bias	src/bias_estimators/esn.py	Correlation-based training; biased_observations = True
ConstantBias	src/bias_estimators/constantbias.py	Persistent bias, reset to the latest innovation each analysis step
NoBias	src/bias_estimators/constantbias.py	ConstantBias fixed at zero; unbiased-limit placeholder

Forecaster

The forecaster attribute of a Bias instance is a Model subclass — typically a data-driven model trained on the residual between observations and model output.

Bias instance
  .forecaster  →  Model instance   (usually ESN_model)

The forecaster is initialised inside init_forecaster() and stepped forward in sync with the main model during Estimator.forecast_step(). Its output is the predicted bias b(t), which enters the analysis step as a correction to the observation.

Basic configuration of the ESN bias estimator: the reservoir forecasts the innovation between observations.

---

romda.bias_estimators.Bias(innovation, t, dt, **kwargs)

Base class for the model-bias estimators used in bias-aware data assimilation.

A bias estimator provides three things to the assimilation loop:

1. a forecast of the bias between analyses (time_integrate), driven by its internal forecaster (an ESN, a constant map, a linear model, ...);
2. the Jacobian of the bias with respect to the observables (state_derivative), \(\mathbf{J} = \mathrm{d}\mathbf{b}/\mathrm{d}\mathbf{q}\), required by the regularized bias-aware EnKF;
3. an update rule from the analysis innovation (update_state_from_innovation), optionally Bayesian (an internal EnSRKF on the bias state).

The estimator state has \(N_\mathrm{dim}\) components: \([\mathbf{b}]\) if the observations are unbiased, or \([\mathbf{b}; \mathbf{i}]\) (bias and innovations, \(N_\mathrm{dim} = 2 N_q\)) if biased_observations is set. Child classes may add hidden components (e.g., the ESN reservoir).

Parameters:

Name	Type	Description	Default
innovation	ndarray	Initial innovation/bias estimate, used to set the observable dimension \(N_q\).	required
t	float	Initial time.	required
dt	float	Time step of the output history.	required
**kwargs		Class-attribute overrides (see Attributes) and forecaster options.	{}

Attributes:

Name	Type	Description
upsample	int	Upsampling factor of the internal forecaster time step relative to dt.
L	int	Number of trajectories in the training dataset (data-driven estimators).
augment_data	bool or int	Whether (and how much) to augment the training data.
bayesian_update	bool	If True, the innovation update is a Bayesian (EnSRKF) update of the full estimator state; otherwise the innovation is assigned directly.
biased_observations	bool	If True, the observations themselves are assumed biased and the estimator tracks bias and innovations separately.
force_retrain	bool	If True, retrain the forecaster even if a cached configuration exists.

Source code in src/bias_estimators/bias.py

def __init__(self, innovation, t, dt, **kwargs):


    # ===================== ASSIGN PROVIDED KWARGS ======================= ##
    keys = list(kwargs.keys())
    [setattr(self, key, kwargs.pop(key)) for key in keys if hasattr(self, key)]

    self.keys_to_print += self.extra_keys_to_print

    # ================== Setup dimensions ================= ##

    self.precision_t = int(-np.log10(dt)) + 2
    self.dt = dt
    self.Nq = self._format_state(innovation).shape[1]


    # ================== Initialize Forecaster & HISTORY ================= ##

    self.init_forecaster(**kwargs)

    bias_state = self.initialize_bias_state

    assert bias_state.shape[-2] == self.N, f"Bias state shape {bias_state.shape} does not match expected (Nt, N = {self.N}, Nens)."
    self.update_history(bias_state, t=t, reset=True)

initialize_bias_state property

Only used at initialization. If the forecaster is a model, this should be handled by the child class. The state has N_dim components: [bias] if the observations are unbiased, or [bias; innovations] if the observations are biased (N_dim = 2 * Nq).

integrator property

This is the integrator used by the model of the bias. E.g., DiscreteIntegrator if using ESN_model as forecaster.

washout_data property writable

Returns the washout data used for initializing the bias model, which is typically obtained from the washout phase using the validation data. This property can be used to access the washout data for further processing or analysis.

Returns:

Type	Description
Tuple of (washout_data, washout_time) where:	washout_data: np.ndarray - Washout data used for initializing the bias model. washout_time: np.ndarray - Time points corresponding to the washout data.

hist property

Returns only the valid (non-empty) portion of the history buffer.

hist_t property

Returns only the valid portion of the time history.

current_state property

Returns the current state (last entry in history).

current_time property

Returns the current time (last entry in time history).

current_bias property

Returns the current (ensemble-mean) bias computed from the current state. Shape: (Nq, 1) -- the bias is defined on the ensemble mean.

current_innovations property

Returns the current (ensemble-mean) innovations computed from the current state. Shape: (Nq, 1).

DA_method property

Cached measurement operator M for the Bayesian innovation update (_ensrkf_update).

init_forecaster(**kwargs)

.....

Source code in src/bias_estimators/bias.py

def init_forecaster(self, **kwargs):
    """
    .....
    """
    raise NotImplementedError('Bias child classes must implement _init_forecaster() method.')

state_derivative()

Returns the derivative of the bias state, which is used for time integration. This is computed by the forecaster model.

Source code in src/bias_estimators/bias.py

def state_derivative(self):
    """
    Returns the derivative of the bias state, which is used for time integration.
    This is computed by the forecaster model.
    """
    raise NotImplementedError('Bias child classes must implement state_derivative property, typically computed by the forecaster model.')

washout_phase(d_wash, t_wash, **kwargs)

Optional method to initialize the bias model if needed, e.g., by running a washout phase with given data. By default, does nothing, but can be implemented in child classes if needed.

Source code in src/bias_estimators/bias.py

def washout_phase(self, d_wash, t_wash, **kwargs) -> Optional[tuple[np.ndarray, np.ndarray]]:
    """
    Optional method to initialize the bias model if needed, e.g., by running a washout phase with given data.
    By default, does nothing, but can be implemented in child classes if needed.
    """
    return None

update_history_aux(**kwargs)

Auxiliary method to update any additional history attributes in child classes if needed.

Source code in src/bias_estimators/bias.py

def update_history_aux(self, **kwargs):
    """Auxiliary method to update any additional history attributes in child classes if needed."""
    pass

update_state_from_innovation(input_innovation)

Optional method to perform a Bayesian update to the state using the bias model. This can be implemented in child classes if needed, e.g., for ESN bias model. By default, does nothing, but can be implemented in child classes if needed.

Parameters:

Name	Type	Description	Default
input_innovation	ndarray	Analysis innovation ensemble, shape (Nq, m), (Nq, 1) or (Nq,).	required

Returns:

Type	Description
ndarray	Updated estimator state, shape (N, N_ens).

Source code in src/bias_estimators/bias.py

def update_state_from_innovation(self, input_innovation):
    """
    Optional method to perform a Bayesian update to the state using the bias model. This can be implemented in child classes if needed, e.g., for ESN bias model.
    By default, does nothing, but can be implemented in child classes if needed.

    Parameters
    ----------
    input_innovation : np.ndarray
        Analysis innovation ensemble, shape ``(Nq, m)``, ``(Nq, 1)`` or ``(Nq,)``.

    Returns
    -------
    np.ndarray
        Updated estimator state, shape ``(N, N_ens)``.
    """
    input_innovation = self._format_state(input_innovation) # Ensure shape (nt, nstate, nens)
    assert input_innovation.shape[0] == 1, "Input innovation must have only one time step (shape[0] == 1) for state_from_innovation method."

    forecast_state = self.current_state


    if self.bayesian_update:
        mean_innovation = np.mean(input_innovation[0], axis=-1)  # Average innovation across ensemble (obs_dim, Nens) -> (obs_dim,)
        # cov_innovation = (inn_uncertainty * np.max(np.abs(mean_innovation)))**2 * np.eye(mean_innovation.shape[0])  # Diagonal covariance of innovation (obs_dim, obs_dim)
        # cov_innovation = np.cov(input_innovation[0], rowvar=True)  # Full covariance of innovation (obs_dim, obs_dim)
        # cov_innovation = np.atleast_2d(cov_innovation)
        conv_inn = input_innovation[0] - mean_innovation[:, np.newaxis]  # Centered innovations (obs_dim, Nens)
        cov_innovation = np.atleast_2d(np.cov(conv_inn, rowvar=True))  # Full covariance of centered innovations (obs_dim, obs_dim)

        updated_state = _ensrkf_update(Af=forecast_state, d=mean_innovation, Cdd=cov_innovation, M=self.DA_method)
    else:
        updated_state = forecast_state.copy()
        mean_innovation = np.mean(input_innovation[0], axis=-1, keepdims=True)  # Average innovation across ensemble (obs_dim, Nens) -> (obs_dim, 1)
        if input_innovation.shape[-1] == self.N_ens:
            # One innovation per bias-ensemble member: assign directly
            updated_state[self.observed_idx, :] = input_innovation[0]

        elif self.N_ens == 1 or input_innovation.shape[-1] == 1:
            # Assign the same mean innovation to all ensemble members
            updated_state[self.observed_idx, :] = np.repeat(mean_innovation, self.N_ens, axis=-1)

        else:  # ensemble sizes differ => resample from the innovation statistics
            cov_innovation = np.cov(input_innovation[0], rowvar=True)
            mean_innovation = mean_innovation.flatten()
            cov_innovation = np.atleast_2d(cov_innovation)
            resampled_innovation = np.random.multivariate_normal(mean_innovation, cov_innovation, size=self.N_ens).T  # Resample innovations for each ensemble member (obs_dim, Nens)
            updated_state[self.observed_idx, :] = resampled_innovation

    return updated_state

close()

Close any resources used by the bias model, e.g., forecaster model.

Source code in src/bias_estimators/bias.py

def close(self):
    """Close any resources used by the bias model, e.g., forecaster model."""
    if hasattr(self, '_forecaster') and hasattr(self._forecaster, 'close'):
        self._forecaster.close()

romda.bias_estimators.ESN_bias(rom, reference_data=None, **kwargs)

Bases: Bias

Echo-state-network bias estimator.

An ESN_model forecasts the model bias (and, if biased_observations, the innovations) in closed loop between analyses, and its open-loop linearization provides the Jacobian \(\mathbf{J} = \mathrm{d}\mathbf{b}/\mathrm{d}\mathbf{q}\) used by the regularized bias-aware EnKF. The network is trained offline on synthetic innovation data generated by perturbing the low-order model (see create_bias_training_dataset); trained configurations are cached on disk and reloaded when the configuration hash matches (unless force_retrain).

Parameters:

Name	Type	Description	Default
rom	Model	The (biased) low-order forecast model whose bias is estimated.	required
reference_data	Observations or list of Observations	Reference data used to build the training dataset if no cached forecaster or dataset is found.	None
**kwargs		ESN hyperparameters (N_units, N_wash, rho_range, ...), training times (t_train, t_val, t_test), dataset options (L, augment_data, correlation_based_training) and Bias options.	{}

References

Nóvoa, Racca & Magri (2023). Inferring unknown unknowns: Regularized bias-aware ensemble Kalman filter. Comput. Methods Appl. Mech. Eng., 418, 116502. DOI: 10.1016/j.cma.2023.116502.

Source code in src/bias_estimators/esn.py

def __init__(self,
             rom: Model,
             reference_data = None,
             **kwargs):


    super().__init__(
                    innovation=kwargs.pop('innovation', np.zeros((rom.Nq,))),
                    dt=kwargs.pop('dt', rom.dt),
                    t=kwargs.pop('t', rom.t),
                    rom=rom,
                    reference_data=reference_data,
                    **kwargs)

    assert isinstance(self.forecaster, ESN_model), "forecaster must be an instance of ESN_model"

initialize_bias_state property

Initialize the ESN reservoir state from the validation data.

Called during Bias.__init__ to set the initial state of the ESN_model forecaster from its validation data, which improves training and closed-loop performance relative to a zero/random initial state.

Returns:

Type	Description
ndarray	Initialized reservoir state for the ESN_model forecaster, for N_ens ensemble members.

washout_phase(d_wash, t_wash, **kwargs)

Run an open-loop washout to initialize the reservoir from real data.

Parameters:

Name	Type	Description	Default
d_wash	ndarray	Washout data used to initialize the bias model, shape (Nt, Nq) or (Nq, Nt).	required
t_wash	ndarray	Time points corresponding to d_wash.	required

Returns:

Name	Type	Description
psi	ndarray	Bias state trajectory after washout (excluding the initial state).
t_wash	ndarray	Time points of psi.

Source code in src/bias_estimators/esn.py

def washout_phase(self, d_wash, t_wash, **kwargs):
    """Run an open-loop washout to initialize the reservoir from real data.

    Parameters
    ----------
    d_wash : ndarray
        Washout data used to initialize the bias model, shape ``(Nt, Nq)`` or
        ``(Nq, Nt)``.
    t_wash : ndarray
        Time points corresponding to `d_wash`.

    Returns
    -------
    psi : ndarray
        Bias state trajectory after washout (excluding the initial state).
    t_wash : ndarray
        Time points of `psi`.
    """
    assert hasattr(self, 'forecaster'), "Forecaster must be initialized before calling washout_phase."


    esn = self.forecaster #type: ESN_model

    # Make sure first dimension is time
    assert d_wash.ndim == 2, f"Washout data must be a 2D array with shape (Nt, Nq) or (Nq, Nt), got {d_wash.shape}."
    if d_wash.shape[0] != len(t_wash):
        u_wash = d_wash.copy().T
    else:
        u_wash = d_wash.copy()

    # apply upsample and cut if needed to match the washout time points
    u_wash, t_wash = [xx[::esn.upsample][:esn.N_wash+1] for xx in [u_wash, t_wash]]


    # store washout data for potential future plotting
    self.washout_data = (u_wash, t_wash)


    # Get current reservoir state and corresponding physical state from the ESN_model forecaster
    Nt = len(t_wash) + 1
    r_open = esn.reservoir_state
    u_out, r_out = np.empty((Nt, self.N_dim, r_open.shape[1])), np.empty((Nt, *r_open.shape))

    r_out[0]= esn.reservoir_state
    u_out[0] = esn.reservoir_to_physical(r_out[0])

    # Open-loop reservoir
    for kk in range(len(t_wash)):
        u_open, r_open = esn.step(u_wash[kk], r_out[kk])
        u_out[kk+1], r_out[kk+1] = u_open, r_open

    #store final state into the initialization arrays

    psi = esn.build_psi(u=u_out, r=r_out)
    return psi[1:], t_wash

init_forecaster(training_data_filename=None, **kwargs)

Load or create the ESN_model forecaster and store it as self.forecaster.

Parameters:

Name	Type	Description	Default
training_data_filename	str	Path to load/save the bias-training dataset (see load_or_create_bias_training_dataset).	None
**kwargs		ESN hyperparameters and other options, forwarded to load_or_create_forecaster. Must include rom (the low-order model whose bias is estimated).	{}

Source code in src/bias_estimators/esn.py

def init_forecaster(self,
                    training_data_filename: Optional[str] = None,
                    **kwargs):
    """Load or create the `ESN_model` forecaster and store it as `self.forecaster`.

    Parameters
    ----------
    training_data_filename : str, optional
        Path to load/save the bias-training dataset (see
        `load_or_create_bias_training_dataset`).
    **kwargs
        ESN hyperparameters and other options, forwarded to
        `load_or_create_forecaster`. Must include `rom` (the low-order model
        whose bias is estimated).
    """
    cfg = self.config.copy()
    cfg.update(kwargs)

    cfg['N_dim'] = self.N_dim
    cfg['training_data_filename'] = training_data_filename
    cfg['dt'] = self.dt
    rom = kwargs.get('rom')

    assert rom is not None, "ROM object must be provided."

    # add training times if not provided in kwargs, with default values based on the ROM time scales
    t_test_default = 5 * rom.t_CR if kwargs.get('perform_test', True) else 0
    for key, default_value in zip(['t_train', 't_val', 't_test'], [rom.t_transient / 2, rom.t_CR, t_test_default]):
        if key not in cfg.keys():
            cfg[key] = kwargs.get(key, default_value)

    cfg['rom'] = rom

    min_training_time = sum([cfg[key] for key in ['t_train', 't_val', 't_test']])
    self.minimum_training_steps = int(np.ceil(min_training_time / self.dt))

    if not hasattr(self, 'L'):
        self.L = rom.m

    # Load or create training dataset for bias model
    self.forecaster = self.load_or_create_forecaster(**cfg) #type: ESN_model

load_or_create_forecaster(hash=None, reference_data=None, rom=None, **kwargs)

Load a cached ESN_model forecaster from disk, or train a new one.

Parameters:

Name	Type	Description	Default
hash	str	Hash identifying the ESN_model configuration to load. If None, it is derived from kwargs via ESNConfig.from_init_params.	None
reference_data	Observations or list of Observations	Reference data used to build the training dataset if no cached forecaster/dataset is found. Required (with rom) when training a new forecaster.	None
rom	Model	The (biased) low-order model to sample training states from. Required (with reference_data) when training a new forecaster.	None
**kwargs		ESN hyperparameters and dataset options, used both to compute the configuration hash and (if training) to construct the new ESN_model.	{}

Returns:

Type	Description
ESN_model	A trained forecaster, either loaded from cache or newly trained.

Source code in src/bias_estimators/esn.py

def load_or_create_forecaster(self, hash=None, reference_data=None, rom=None, **kwargs) -> ESN_model:
    """Load a cached `ESN_model` forecaster from disk, or train a new one.

    Parameters
    ----------
    hash : str, optional
        Hash identifying the `ESN_model` configuration to load. If None, it is
        derived from `kwargs` via `ESNConfig.from_init_params`.
    reference_data : Observations or list of Observations, optional
        Reference data used to build the training dataset if no cached
        forecaster/dataset is found. Required (with `rom`) when training a new
        forecaster.
    rom : Model, optional
        The (biased) low-order model to sample training states from. Required
        (with `reference_data`) when training a new forecaster.
    **kwargs
        ESN hyperparameters and dataset options, used both to compute the
        configuration hash and (if training) to construct the new `ESN_model`.

    Returns
    -------
    ESN_model
        A trained forecaster, either loaded from cache or newly trained.
    """

    cfg = kwargs.copy()

    if hash is None:
        query_config = ESNConfig.from_init_params(**kwargs)
        query_hash = query_config.to_hash()
    else:
        query_hash = hash

    # Try to load forecaster configuration from disk
    loaded_case = load_esn_model_from_config(q=query_hash)


    if loaded_case is not None and not self.force_retrain:
        assert isinstance(loaded_case, ESN_model), f'Loaded case must be an instance of ESN_model, but got {type(loaded_case)}.'
        assert loaded_case.trained is True, f'{loaded_case.name} model must be trained after initialization.'
        return loaded_case

    elif rom is None or reference_data is None:
        raise ValueError('Both rom and reference_data must be provided to create a new ESN_bias model')

    else:
        # Create new forecaster. First, create training dataset for bias model, then use it to create the forecaster.
        train_data_dict = self.load_or_create_bias_training_dataset(training_data_filename=kwargs.get('training_data_filename'),
                                                                    rom=rom,
                                                                    reference_data=reference_data,
                                                                    std_alpha=cfg.get('std_alpha'),
                                                                    std_phi=cfg.get('std_phi'))
        if train_data_dict is None:
            raise RuntimeError('Failed to load or create training data for bias model.')

        # Create a new instance of the ESN_model to use as forecaster
        cfg.update(train_data_dict)

        # print('CONFIGURATION')
        # print('data ', cfg['data'].shape)
        # print('state ', cfg['state'].shape)
        # print('------------------')

        new_esn_model = ESN_model(**cfg) # Note: ESN_model trains itself during initialization using the provided training data, so we don't need a separate training step here. If the ESN_model implementation changes in the future to require a separate training step, this code will need to be updated accordingly.

        # Save under the QUERY hash (as `auto_load_or_create` does): re-deriving the
        # hash from the trained model would use the t_train/t_val/t_test mutated by
        # `_process_initialization_data`, so the next lookup would miss forever.
        _ = save_esn_model_to_config(new_esn_model, name=f"{query_hash}")

        return new_esn_model

load_or_create_bias_training_dataset(training_data_filename=None, rom=None, reference_data=None, std_phi=None, std_alpha=None)

Load the bias-training dataset from disk, or create and cache a new one.

Delegates to load_bias_training_dataset / create_bias_training_dataset (defined in romda.bias_estimators.aux), which handle preprocessing, augmentation, and formatting for training the bias model.

Parameters:

Name	Type	Description	Default
training_data_filename	str	Filename to load/save the training dataset.	None
rom	Model	ROM to sample states from, if the dataset needs to be created.	None
reference_data	Observations or list of Observations	Reference data, if the dataset needs to be created.	None
std_phi	float	Standard deviation for sampling initial conditions, if the dataset needs to be created; see sample_model_states.	None
std_alpha	float or dict	Standard deviation for sampling parameters, if the dataset needs to be created; see sample_model_states.	None

Returns:

Type	Description
dict	Training-data dictionary; see create_bias_training_dataset.

Source code in src/bias_estimators/esn.py

def load_or_create_bias_training_dataset(self,
                                         training_data_filename=None,
                                         rom=None,
                                         reference_data=None,
                                         std_phi=None,
                                         std_alpha=None):
    """Load the bias-training dataset from disk, or create and cache a new one.

    Delegates to `load_bias_training_dataset` / `create_bias_training_dataset`
    (defined in `romda.bias_estimators.aux`), which handle preprocessing,
    augmentation, and formatting for training the bias model.

    Parameters
    ----------
    training_data_filename : str, optional
        Filename to load/save the training dataset.
    rom : Model, optional
        ROM to sample states from, if the dataset needs to be created.
    reference_data : Observations or list of Observations, optional
        Reference data, if the dataset needs to be created.
    std_phi : float, optional
        Standard deviation for sampling initial conditions, if the dataset needs
        to be created; see `sample_model_states`.
    std_alpha : float or dict, optional
        Standard deviation for sampling parameters, if the dataset needs to be
        created; see `sample_model_states`.

    Returns
    -------
    dict
        Training-data dictionary; see `create_bias_training_dataset`.
    """

    expected_Ndim = None
    if rom is not None:
        expected_Ndim = 2 * rom.Nq if self.biased_observations else rom.Nq

    train_data_dict = load_bias_training_dataset(filename=training_data_filename,
                                                necessary_properties=self.config.copy(),
                                                minimum_training_steps=self.minimum_training_steps,
                                                augment_data_length=self.augment_data_length,
                                                L=self.L,
                                                expected_Ndim=expected_Ndim)


    if train_data_dict is not None:
        assert isinstance(train_data_dict, dict), 'ERROR: Loaded training data for bias model must be a dictionary.'
        return train_data_dict

    print('Creating training data for bias model...')

    assert rom is not None and reference_data is not None, 'ROM and reference data must be provided to create bias training dataset.'



    train_data_dict = create_bias_training_dataset(
                            config=self.config,
                            rom=rom,
                            reference_data=reference_data,
                            std_phi=std_phi,
                            std_alpha=std_alpha,
                            L=self.L,
                            augment_data_length=self.augment_data_length,
                            minimum_training_steps=self.minimum_training_steps,
                            correlation_based_training=self.correlation_based_training,
                            biased_observations=self.biased_observations,
                        )

    if training_data_filename is not None:
        save_to_pickle_file(training_data_filename, train_data_dict)

    return train_data_dict

romda.bias_estimators.ConstantBias(innovation, t, dt, k=None, **kwargs)

Bases: Bias

Constant (persistent) bias estimator.

The bias is held constant between analysis steps, i.e., the forecast model of the bias is \(\dot{\mathbf{b}} = \mathbf{0}\). At each analysis step, the bias state is reset to the latest innovation (see Bias.update_state_from_innovation). This is the classic persistent-bias assumption.

Parameters:

Name	Type	Description	Default
innovation	ndarray	Initial innovation/bias estimate, shape (Nq,), (Nq, N_ens) or (1, Nq, N_ens).	required
t	float	Initial time.	required
dt	float	Time step of the output history.	required
k	float or ndarray	If provided, the initial bias state is set to the constant value(s) k instead of the provided innovation.	None

Source code in src/bias_estimators/constantbias.py

def __init__(self, innovation, t, dt, k=None, **kwargs):

    if k is not None:
        innovation = np.ones(np.shape(innovation)) * k  # Constant initial bias

    self._b0 = np.asarray(innovation, dtype=float)

    super().__init__(innovation=innovation, t=t, dt=dt, **kwargs)

init_forecaster(**kwargs)

Constant forecaster: no underlying model, the state is simply held in time.

Source code in src/bias_estimators/constantbias.py

def init_forecaster(self, **kwargs):
    """
    Constant forecaster: no underlying model, the state is simply held in time.
    """
    self._forecaster = SimpleNamespace()

    initial_capacity = kwargs.pop('initial_capacity', 1000)
    self._forecaster.history = HistoryTracker(initial_capacity=initial_capacity)

    #  INITIALISE INTEGRATOR STRATEGY ================== ##
    self._forecaster.integrator = ConstantIntegrator(self)

romda.bias_estimators.NoBias(innovation, t, dt, **kwargs)

Bases: ConstantBias

Placeholder bias estimator that always returns zero bias. Useful to run the bias-aware machinery in its unbiased limit.

Source code in src/bias_estimators/constantbias.py

def __init__(self, innovation, t, dt, **kwargs):
    kwargs.pop('k', None)
    super().__init__(innovation=innovation, t=t, dt=dt, k=0., **kwargs)

romda.bias_estimators.aux.create_bias_training_dataset(config, rom, reference_data, minimum_training_steps, L, correlation_based_training, augment_data_length, biased_observations, std_phi=None, std_alpha=None)

Build a training dataset for the bias model from ROM samples and reference data.

Parameters:

Name	Type	Description	Default
config	dict	Configuration for creating the training dataset; echoed into the returned dictionary.	required
rom	Model	The reduced-order model to sample states from.	required
reference_data	Observations or list of Observations	Reference data to prepare for training.	required
minimum_training_steps	int	Minimum number of time steps required for training the bias model.	required
L	int	Number of samples to generate from the ROM.	required
correlation_based_training	bool	If True, correlate the model-generated data with the raw observations (via correlate_data) rather than pairing them directly at matching time steps.	required
augment_data_length	int	Number of augmented samples per observed variable: 1 uses only the best lag/direct pairing; 2 adds a mid-point (or scaled) sample; 3 also adds the worst lag (or oppositely-scaled) sample.	required
biased_observations	bool	Whether to also include the model bias (true minus model-generated data) in the training dataset, alongside the innovations.	required
std_phi	float	Standard deviation for sampling initial conditions; see sample_model_states.	None
std_alpha	float or dict	Standard deviation for sampling parameters; see sample_model_states.	None

Returns:

Type	Description
dict	Training data for the bias model, plus the entries of config: 'data' : ndarray, shape (L * augment_data_length, minimum_training_steps, Ndim) where Ndim = Nq if not biased_observations, else 2 * Nq. 'y_model' : ndarray, shape (L, minimum_training_steps, Nq) — the model-generated data used for training.

Source code in src/bias_estimators/aux.py

def create_bias_training_dataset(config: dict,
                                rom: Model,
                                reference_data,
                                minimum_training_steps: int,
                                L: int,
                                correlation_based_training: bool,
                                augment_data_length: int,
                                biased_observations: bool,
                                std_phi: Optional[float] = None,
                                std_alpha: Optional[Union[float, dict[str, Union[float, list[float]]]]] = None,
                            ) -> dict:

    """Build a training dataset for the bias model from ROM samples and reference data.

    Parameters
    ----------
    config : dict
        Configuration for creating the training dataset; echoed into the returned
        dictionary.
    rom : Model
        The reduced-order model to sample states from.
    reference_data : Observations or list of Observations
        Reference data to prepare for training.
    minimum_training_steps : int
        Minimum number of time steps required for training the bias model.
    L : int
        Number of samples to generate from the ROM.
    correlation_based_training : bool
        If True, correlate the model-generated data with the raw observations (via
        `correlate_data`) rather than pairing them directly at matching time steps.
    augment_data_length : int
        Number of augmented samples per observed variable: 1 uses only the best
        lag/direct pairing; 2 adds a mid-point (or scaled) sample; 3 also adds the
        worst lag (or oppositely-scaled) sample.
    biased_observations : bool
        Whether to also include the model bias (true minus model-generated data) in
        the training dataset, alongside the innovations.
    std_phi : float, optional
        Standard deviation for sampling initial conditions; see `sample_model_states`.
    std_alpha : float or dict, optional
        Standard deviation for sampling parameters; see `sample_model_states`.

    Returns
    -------
    dict
        Training data for the bias model, plus the entries of `config`:

        - ``'data'`` : ndarray, shape (L * augment_data_length, minimum_training_steps, Ndim)
          where ``Ndim = Nq`` if not `biased_observations`, else ``2 * Nq``.
        - ``'y_model'`` : ndarray, shape (L, minimum_training_steps, Nq) — the
          model-generated data used for training.
    """

    y_model_L = sample_model_states(rom=rom,
                                    L=L,
                                    minimum_training_steps=minimum_training_steps,
                                    std_phi=std_phi,
                                    std_alpha=std_alpha)

    print('\n\n Preparing reference data for training...')
    print('y_model_L shape:', y_model_L.shape)

    y_raw, y_true = prepare_reference_data(reference_data, minimum_training_steps=minimum_training_steps)
    augment = augment_data_length > 1

    if not correlation_based_training:
        innovations_all, model_bias_all = [], []
        for yr, yt in zip(y_raw, y_true):
            innovations = (yr - y_model_L[-minimum_training_steps:]).transpose((2, 0, 1))
            innovations_all.append(innovations)

            if augment:
                innovations_all.append(innovations * 1e-1)
                innovations_all.append(innovations * -1e-2)

            if biased_observations:
                model_bias = (yt - y_model_L[-minimum_training_steps:]).transpose((2, 0, 1))
                model_bias_all.append(model_bias)
                if augment:
                    model_bias_all.append(model_bias * 1e-1)
                    model_bias_all.append(model_bias * -1e-2)
    else:
        innovations_all, model_bias_all = [], []
        y_model = []
        for yr, yt in zip(y_raw, y_true):
            ym_L = correlate_data(y_model_L, yr, augment_data_length, minimum_training_steps)
            y_model.append(ym_L.copy())

            innovations = (yr - ym_L).transpose((2, 0, 1)) #shape (L, minimum_training_steps, Nq)
            innovations_all.append(innovations)


            if biased_observations:
                model_bias = (yt - ym_L).transpose((2, 0, 1))
                model_bias_all.append(model_bias)
        y_model_L = np.concatenate(y_model, axis=0)

    if not biased_observations:
        train_data = np.concatenate(innovations_all, axis=0)
    else:
        innovations_all = np.concatenate(innovations_all, axis=0)
        model_bias_all = np.concatenate(model_bias_all, axis=0)
        train_data = np.concatenate([model_bias_all, innovations_all], axis=2)
    # Save to dictionary
    train_data_dict = {key: val for key, val in config.items()}
    train_data_dict.update(data=train_data,
                           y_model=y_model_L)
    return train_data_dict

romda.bias_estimators.aux.sample_model_states(rom, L, minimum_training_steps, std_phi=None, std_alpha=None)

Sample model states from the ROM to build a training dataset for the bias model.

1. Initializes the ROM with an ensemble of states (sampled initial conditions and, if specified, parameters).
2. Integrates the ROM forward in time to generate model data.
3. Processes this data into a training dataset for the bias model.

Parameters:

Name	Type	Description	Default
rom	Model	The reduced-order model to sample states from.	required
L	int	Number of samples (ensemble members) to generate.	required
minimum_training_steps	int	Minimum number of time steps to integrate the ROM for, to generate enough data for training the bias model.	required
std_phi	float	Standard deviation for sampling the initial conditions. If None (default), the standard deviation of the ROM's current state is used.	None
std_alpha	float or dict	Standard deviation for sampling the parameters. A float is used as a multiplier on the current-state standard deviation, giving a [mean - std_alpha*std, mean + std_alpha*std] range for each estimated parameter; a dict should have parameter names as keys and per-parameter standard deviations as values. If None (default), the range is taken from the current state of the ROM.	None

Returns:

Type	Description
ndarray	Observable history of the (re-)integrated ROM ensemble, used as the model data for training.

Source code in src/bias_estimators/aux.py

@typechecked
def sample_model_states(rom: Model,
                        L: int,
                        minimum_training_steps: int,
                        std_phi: Optional[Union[float, np.ndarray]] = None,
                        std_alpha: Optional[Union[float, dict[str, Union[float, list[float], tuple[float, float]]]]] = None,
                    ) -> np.ndarray:
    """Sample model states from the ROM to build a training dataset for the bias model.

    1. Initializes the ROM with an ensemble of states (sampled initial conditions and,
       if specified, parameters).
    2. Integrates the ROM forward in time to generate model data.
    3. Processes this data into a training dataset for the bias model.

    Parameters
    ----------
    rom : Model
        The reduced-order model to sample states from.
    L : int
        Number of samples (ensemble members) to generate.
    minimum_training_steps : int
        Minimum number of time steps to integrate the ROM for, to generate enough
        data for training the bias model.
    std_phi : float, optional
        Standard deviation for sampling the initial conditions. If None (default),
        the standard deviation of the ROM's current state is used.
    std_alpha : float or dict, optional
        Standard deviation for sampling the parameters. A float is used as a
        multiplier on the current-state standard deviation, giving a
        ``[mean - std_alpha*std, mean + std_alpha*std]`` range for each estimated
        parameter; a dict should have parameter names as keys and per-parameter
        standard deviations as values. If None (default), the range is taken from
        the current state of the ROM.

    Returns
    -------
    ndarray
        Observable history of the (re-)integrated ROM ensemble, used as the model
        data for training.
    """

    model = rom.copy()

    if std_phi is None:
        std_phi = np.std(model.current_state[:model.Nphi, :], axis=-1)

    if std_alpha is None:
        std_alpha = {}
        for i, key in enumerate(model.est_alpha):
            param = model.current_state[model.Nphi + i, :]
            std_alpha[key] = [min(param), max(param)]
    elif isinstance(std_alpha, float):
        # A scalar std_alpha is a multiplier on the ensemble standard deviation:
        # build a [mean - const*std, mean + const*std] range for each estimated
        # parameter from the current state of the ROM.
        const = std_alpha
        std_alpha = {}
        for i, key in enumerate(model.est_alpha):
            param = model.current_state[model.Nphi + i, :]
            mean_param = np.mean(param)
            std_alpha[key] = [mean_param - const * np.std(param), mean_param + const * np.std(param)]

    assert std_phi is not None, "std_phi must be specified or computed from the model state."
    assert std_alpha is not None, "std_alpha must be specified or computed from the model state."

    def sample_ensemble(psi0_mean, ensemble_size):
        new_phi = mean_vector_to_ensemble(
            rng=model.rng,
            mean_vec=psi0_mean[:model.Nphi],
            std=std_phi,
            m=ensemble_size,
            method='uniform',
        )
        if std_alpha:

            new_alpha = mean_vector_to_ensemble(
                rng=model.rng,
                mean_vec=psi0_mean[model.Nphi:model.Nphi + model.Na],
                std=std_alpha,
                m=ensemble_size,
                method='uniform',
            )

            return np.concatenate([new_phi, new_alpha], axis=0)
        else:
            return new_phi


    psi0 = np.mean(model.current_state.copy(), axis=-1)
    Nt = int(np.round(model.t_transient / model.dt, model.precision_t)) - 1

    # Add parameters to the state vector if parameter uncertainty is given
    if std_alpha and psi0.shape[0] == model.Nphi:
        assert isinstance(std_alpha, dict), "std_alpha must be a dict if parameter uncertainty is specified."
        model.ensemble_cfg = dict(Na=len(std_alpha), est_alpha=list(std_alpha.keys()), m=L)
        psi0 = np.hstack([psi0, np.zeros((len(std_alpha),))])
        for ii, key in enumerate(std_alpha.keys()):
            psi0[model.Nphi + ii] = getattr(model, key)
        psi0_ens = sample_ensemble(psi0, L)
        model.update_history(psi=psi0_ens[np.newaxis, :, :], t=0.0, reset=True)

    elif model.m != L:
        psi0_ens = sample_ensemble(psi0, L)
        model.update_history(psi=psi0_ens[np.newaxis, :, :], t=0.0, reset=True)

    psi, t = model.time_integrate(Nt=Nt)
    model.update_history(psi=psi, t=t, reset=True)

    y_L_model = model.get_observable_hist()
    psi_last = psi[-1, :, :]

    tol = 1e-1
    N_CR = int(round(model.t_CR / model.dt))
    range_y = np.max(np.max(y_L_model[-N_CR:], axis=0) - np.min(y_L_model[-N_CR:], axis=0), axis=0)
    idx_fixed = range_y < tol

    if len(np.flatnonzero(idx_fixed)) / len(idx_fixed) >= 0.2:
        idx_fixed[np.flatnonzero(idx_fixed)[0]] = 0
        psi0 = psi_last[:, ~idx_fixed]
        new_psi0 = sample_ensemble(np.mean(psi0, axis=-1), len(np.flatnonzero(idx_fixed)))
        psi0 = np.concatenate([psi0, new_psi0], axis=-1)
    else:
        psi0 = psi_last

    model.update_history(psi=psi0[np.newaxis, :, :], reset=True)
    psi, t = model.time_integrate(Nt=minimum_training_steps + N_CR)
    model.update_history(psi=psi, t=t, reset=True)
    model.close()

    return model.get_observable_hist()

romda.bias_estimators.plot_train_data(truth, bias_data, t_CR)

Plot the observable and bias training samples against the truth.

Shows one window before the first observation, colouring each training sample (of bias_data) by its bias RMS.

Parameters:

Name	Type	Description	Default
truth	Observations	Reference truth used to select the plotting window and overlay the true observable/bias signals.	required
bias_data	dict	Training-data dictionary as returned by create_bias_training_dataset.	required
t_CR	float	Characteristic (e.g. oscillation) time scale, used to set the plotting window length.	required

Source code in src/bias_estimators/aux.py

def plot_train_data(truth, bias_data, t_CR):
    """Plot the observable and bias training samples against the truth.

    Shows one window before the first observation, colouring each training sample
    (of `bias_data`) by its bias RMS.

    Parameters
    ----------
    truth : Observations
        Reference truth used to select the plotting window and overlay the true
        observable/bias signals.
    bias_data : dict
        Training-data dictionary as returned by `create_bias_training_dataset`.
    t_CR : float
        Characteristic (e.g. oscillation) time scale, used to set the plotting
        window length.
    """
    L, _, _ = bias_data['data'].shape

    Nt = int(t_CR / truth.dt)
    i0_t = np.argmin(np.abs(truth.t_true - truth.t_obs[0]))

    # Build a common valid time window and select the segment before first observation.
    n_common = min(
        len(truth.t_true),
        truth.y_true.shape[0],
        truth.b_true.shape[0],
        bias_data['y_model'].shape[0],
        bias_data['data'].shape[1],
    )
    i_end = min(max(i0_t, 1), n_common)
    i_start = max(i_end - Nt, 0)
    if i_end - i_start < 2:
        i_end = n_common
        i_start = max(i_end - Nt, 0)

    yt = truth.y_true[i_start:i_end]
    bt = truth.b_true[i_start:i_end]
    yr = bias_data['y_model'][i_start:i_end].transpose(2, 0, 1)
    Nq = yt.shape[1]
    br = bias_data['data'][:, i_start:i_end, :Nq]
    tt = truth.t_true[i_start:i_end]

    if len(tt) == 0:
        raise ValueError('Selected plotting window is empty. Check t_CR and training_data dimensions.')


    RS = []
    for ii in range(L):
        RS.append(np.linalg.norm(br[ii][:, 0]) / np.sqrt(len(yt)))

    RS = np.asarray(RS, dtype=float)
    true_RMS = np.linalg.norm(bt[:, 0]) / np.sqrt(len(yt))

    # Plot training data (single row) --------------------------
    fig = plt.figure(figsize    =[12, 2.7], layout='constrained')
    axs = fig.subplots(1, 2)

    # Robust color mapping: clip outliers; if RMS are nearly equal, force distinct member colors.
    if np.ptp(RS) < 1e-12:
        color_values = np.linspace(0.0, 1.0, L)
        norm = Normalize(vmin=0.0, vmax=1.0)
        cmap = plt.cm.ScalarMappable(norm=norm, cmap=plt.get_cmap('viridis'))
        cbar_extend = 'neither'
        cbar_title = 'Member'
    else:
        lo, hi = np.percentile(RS, [5, 95])
        if np.isclose(lo, hi):
            lo = float(np.min(RS))
            hi = float(np.max(RS))
        color_values = np.clip(RS, lo, hi)
        norm = Normalize(vmin=float(lo), vmax=float(hi))
        cmap = plt.cm.ScalarMappable(norm=norm, cmap=plt.get_cmap('viridis'))
        cbar_extend = 'both'
        cbar_title = '$\\mathrm{RMS}$'

    xlim = [tt[0], tt[-1]]

    axs[0].plot(tt, yt[:, 0], color='silver', linewidth=6, alpha=.8)
    axs[1].plot(tt, bt[:, 0], color='silver', linewidth=4, alpha=.8)

    for ii in range(L):
        clr = cmap.to_rgba(color_values[ii])
        axs[0].plot(tt, yr[ii][:, 0], color=clr)
        axs[1].plot(tt, br[ii][:, 0], color=clr)

    axs[0].legend(['Truth'], bbox_to_anchor=(0., 0.25), loc='upper left')
    axs[1].legend([f'True RMS $={true_RMS:.3f}$'], bbox_to_anchor=(0., 0.25), loc='upper left')
    axs[0].set(xlabel='$t$', ylabel='$\\eta$', xlim=xlim)
    axs[1].set(xlabel='$t$', ylabel='$b$', xlim=xlim)

    clb = fig.colorbar(cmap, ax=axs, orientation='vertical', extend=cbar_extend)
    clb.ax.set_title(cbar_title)