Observations

Summary

File: src/observations.py

Standalone class (no inheritance). Generates or loads ground-truth data, applies a configurable manual bias, adds noise, and exposes observation time indices for the DA loop.

Key attributes:

Attribute	Description
y_true	Clean biased truth signal (Nt, Nq, 1)
y_raw	Noisy observed signal (Nt, Nq, 1)
b_true	Applied bias (Nt, Nq, 1)
t_true	Full time vector
y_obs, t_obs	Observations at assimilation times
obs_idx	Indices into t_true at which observations are taken
Nt_obs	Subsampling rate (every Nt_obs steps)

Noise options (noise_type): 'gauss, add', 'gauss, mult', coloured noise variants.

Manual bias options (manual_bias): 'linear', 'periodic', 'time', 'cosine', or any callable f(y_true, t_true) -> (b, name).

Key method: plot_truth(case) — five-panel figure (raw, truth, PDF, PSD, difference).

---

romda.observations.Observations(model=None, **kwargs)

Reference truth and observations for (twin) data assimilation experiments.

The truth can be generated by integrating a Model (twin experiments), loaded from a file (experimental data), or provided directly as arrays. A manual model bias and measurement noise can be added on top, and the observations are sampled from the raw data every Nt_obs time steps between t_start and t_stop.

Parameters:

Name	Type	Description	Default
model	(Model, type[Model], str or None)	Source of the truth: a model instance/class to integrate, a filename to load, or None to provide y_true/y_raw/t_true directly.	None
**kwargs		Time windows (t_start, t_stop, t_max, t_min), sampling (Nt_obs), noise options (add_noise, noise_level, noise_type), bias options (manual_bias: 'linear', 'periodic', 'time', 'cosine' or a callable f(y, t) -> (b, name)), and model parameters.	{}

Attributes:

Name	Type	Description
y_raw	ndarray	Raw (measured) data — biased truth plus noise — shape \((N_t, N_q, L)\), where \(L\) is the number of independent realizations.
y_true	ndarray	Biased truth (without measurement noise), shape \((N_t, N_q, L)\).
b_true	ndarray	Bias added to the truth (zero if no manual bias), shape \((N_t, N_q, L)\).
t_true	ndarray	Time points of the truth, shape \((N_t,)\).
y_obs	ndarray	Observations to assimilate (subset of y_raw), shape \((N_\mathrm{obs}, N_q)\).
t_obs	ndarray	Observation times (subset of t_true), shape \((N_\mathrm{obs},)\).

Notes

Measurement noise (kwargs add_noise, noise_level, noise_type) is drawn either from a Gaussian, \(\epsilon_q(t) \sim \mathcal{N}(0, \texttt{noise\_level}^2)\), or with a prescribed spectral colour (see colour_noise), and combined with the (biased) truth \(\mathbf{y}\) as

\[ \mathbf{y}_\mathrm{raw} = \mathbf{y} + \boldsymbol{\epsilon}\, \max_t|\mathbf{y}| \quad\text{(additive, default)}, \qquad \mathbf{y}_\mathrm{raw} = \mathbf{y}\,(1 + \boldsymbol{\epsilon}) \quad\text{(multiplicative)}, \]

selected via the 'add'/'mult' substring of noise_type.

Initializes the Observations object, loading or creating truth data, applying bias, adding noise, and interpolating to observation times.

Source code in src/observations.py

@typechecked
def __init__(self, model: Union[Model, type[Model], str, None]=None, **kwargs):
    """
    Initializes the Observations object, loading or creating truth data,
    applying bias, adding noise, and interpolating to observation times.
    """
    # 1. Update instance attributes with any passed kwargs
    model_dict = kwargs.copy()
    for key in kwargs.keys():
        if hasattr(Observations, key):
            setattr(self, key, model_dict.pop(key))

    # 2. Generate or Load Truth Data
    if model is None:

        assert ('y_raw' in kwargs or 'y_true' in kwargs) and 't_true' in kwargs, "If model is None, y_raw, y_true, and t_true must be provided as kwargs."

        self.y_raw = kwargs.get('y_raw', None)
        if self.y_raw is None:
            self.y_raw = kwargs.get('y_true')

        self.y_true = kwargs.get('y_true', self.y_raw)

        for key in ['y_raw', 'y_true']:
            val = getattr(self, key)
            if val.ndim == 1:
                val = val[:, np.newaxis, np.newaxis]
            elif val.ndim == 2:
                val = val[:, :, np.newaxis]
            assert val.ndim == 3
            setattr(self, key, val)

        self.t_true = kwargs['t_true']
        self.t_start = kwargs.get('t_start', self.t_true[0]) # type: float
        self.t_stop = kwargs.get('t_stop', self.t_true[-1]) # type: float
        self.name_truth =  kwargs.get('name_truth', 'Truth_Provided')

    else:
        self.y_raw, self.y_true, self.t_true, self.name_truth = self._create_observations(model, **model_dict)

    self.dt = self.t_true[1] - self.t_true[0]

    # 3. Add noise and bias if requested (in this order and only if y_raw is None)
    self._set_bias()
    self._apply_noise()


    # 4. Compute Observation Times
    # Adjust all times by t_min if t_min > 0 (to start t_true[0] at 0)
    if self.t_min > 0:
        self.t_true -= self.t_min
        if self.t_start is not None:
            self.t_start -= self.t_min
        if self.t_stop is not None:
            self.t_stop -= self.t_min

    self.update_obs_idx(self.t_start, self.t_stop, self.Nt_obs)

    # Include washout period if requested
    if kwargs.get('include_washout', False):
        t_wash_0 = self.t_start - 20 * self.dt_obs
        t_wash_end = self.t_start - self.dt_obs

        self.wash_idx = np.arange(np.searchsorted(self.t_true, t_wash_0), np.searchsorted(self.t_true, t_wash_end))

    # Calculate indices
    self._frozen = True  # Freeze attributes to prevent further modification

obs_idx property

Indices into t_true at which observations are sampled.

Raises:

Type	Description
AttributeError	If update_obs_idx has not been called yet.

dt_obs property

Time between consecutive observations, Nt_obs * dt.

y_wash property

Raw data at the washout indices, shape \((N_\mathrm{wash}, N_q)\).

None unless include_washout=True was passed at construction.

t_wash property

Time points of the washout period, shape \((N_\mathrm{wash},)\).

None unless include_washout=True was passed at construction.

y_obs property

Observations to assimilate: y_raw sampled at obs_idx, shape \((N_\mathrm{obs}, N_q)\).

t_obs property

Observation times: t_true sampled at obs_idx, shape \((N_\mathrm{obs},)\).

y_raw property writable

Raw (measured) data — biased truth plus noise. See class docstring for shape.

y_true property writable

Biased truth (without measurement noise). See class docstring for shape.

t_true property writable

Time points of the truth. See class docstring for shape.

name_truth property writable

Descriptive name of the truth data (e.g. used to build result filenames).

update_obs_idx(t_start=None, t_stop=None, Nt_obs=None)

Recompute obs_idx from a time window and sampling stride.

Any argument left as None falls back to the corresponding current attribute (self.t_start, self.t_stop, self.Nt_obs); passing a value updates that attribute in place before recomputing obs_idx.

Parameters:

Name	Type	Description	Default
t_start	float	Start time of the observation window.	None
t_stop	float	End time of the observation window.	None
Nt_obs	int	Number of raw time steps between consecutive observations.	None

Source code in src/observations.py

def update_obs_idx(self, t_start=None, t_stop=None, Nt_obs=None):
    """Recompute `obs_idx` from a time window and sampling stride.

    Any argument left as ``None`` falls back to the corresponding current
    attribute (``self.t_start``, ``self.t_stop``, ``self.Nt_obs``); passing a
    value updates that attribute in place before recomputing `obs_idx`.

    Parameters
    ----------
    t_start : float, optional
        Start time of the observation window.
    t_stop : float, optional
        End time of the observation window.
    Nt_obs : int, optional
        Number of raw time steps between consecutive observations.
    """

    # Use existing attributes if parameters are not provided
    for key, val in zip(['t_start', 't_stop', 'Nt_obs'], [t_start, t_stop, Nt_obs]):
        if val is None:
            val = getattr(self, key)
        else:
            setattr(self, key, val)

    assert self.t_start is not None and self.t_stop is not None and self.Nt_obs is not None, "t_start, t_stop, and Nt_obs must be defined to update obs_idx."
    start_idx = np.searchsorted(self.t_true, self.t_start)
    stop_idx = np.searchsorted(self.t_true, self.t_stop, side='right') - 1

    self._obs_idx =  np.arange(start_idx, stop_idx + 1, self.Nt_obs, dtype=int)

plot_truth(case, Nq=None, fig_width=12, window=None, f_max=None) staticmethod

Plot raw vs. true time series, PDFs, PSDs and their difference.

Produces a grid with one row per observable and 5 columns: raw signal, true signal, PDF, PSD, and the raw-minus-true difference (bias overlaid if present).

Parameters:

Name	Type	Description	Default
case	Observations	Instance providing y_raw, y_true, t_true, y_obs, t_obs, b_true and (optionally) y_wash/t_wash.	required
Nq	int	Number of observables to plot. Defaults to all observables in case.	None
fig_width	float	Figure width in inches. Default 12.	12
window	float	Length of the time window shown in the time-domain plots. Defaults to a fixed fraction of the available history.	None
f_max	float	Maximum frequency shown on the PSD plots.	None

Source code in src/observations.py

@staticmethod
def plot_truth(case, Nq=None, fig_width=12, window=None, f_max=None):
    """Plot raw vs. true time series, PDFs, PSDs and their difference.

    Produces a grid with one row per observable and 5 columns: raw signal, true
    signal, PDF, PSD, and the raw-minus-true difference (bias overlaid if
    present).

    Parameters
    ----------
    case : Observations
        Instance providing `y_raw`, `y_true`, `t_true`, `y_obs`, `t_obs`,
        `b_true` and (optionally) `y_wash`/`t_wash`.
    Nq : int, optional
        Number of observables to plot. Defaults to all observables in `case`.
    fig_width : float, optional
        Figure width in inches. Default 12.
    window : float, optional
        Length of the time window shown in the time-domain plots. Defaults to a
        fixed fraction of the available history.
    f_max : float, optional
        Maximum frequency shown on the PSD plots.
    """

    # 1. Data Extraction and Setup
    keys = ['y_raw', 'y_true', 't_true', 'y_obs', 't_obs', 'b_true', 'y_wash', 't_wash']

    y_raw, y_true, t_true, y_obs, t_obs, b, y_wash, t_wash = tuple((val.squeeze() if val is not None else None)
                                                                    for key in keys
                                                                    for val in [getattr(case, key)])

    assert isinstance(y_true, np.ndarray), "y_true is required for plotting but is not available in the case data."
    assert isinstance(y_raw, np.ndarray), "y_raw is required for plotting but is not available in the case data."
    assert isinstance(t_true, np.ndarray), "t_true is required for plotting but is not available in the case data."
    assert isinstance(b, np.ndarray), "b_true is required for plotting but is not available in the case data."

    if y_true.ndim == 1:
        y_true = y_true[:, np.newaxis]
        y_raw = y_raw[:, np.newaxis]
        if y_obs is not None:
            y_obs = y_obs[:, np.newaxis]
        if y_wash is not None:
            y_wash = y_wash[:, np.newaxis]


    dt = t_true[1] - t_true[0]
    # Calculate noise: Difference between raw and true signal
    noise = y_raw - y_true


    if Nq is None:
        Nq = y_true.shape[1]

    # Compute PSDs
    # find first index for t_obs
    if hasattr(case, 'wash_idx'):
        t0 = case.t_wash[0] - case.dt_obs * 2
    else:
        t0 = case.t_obs[0] - case.dt_obs * 10

    t0_idx = np.argmin(np.abs(case.t_true -  t0))  # Start a bit before the first observation to capture initial conditions in PSD

    nt_PSD = int((len(t_true) - t0_idx) // 2)
    f_raw, PSD_raw = fun_PSD(dt, y_raw[t0_idx:nt_PSD + t0_idx])
    _, PSD_true = fun_PSD(dt, y_true[t0_idx:nt_PSD + t0_idx])

    # Determine plotting time window (simplified: use the first X data points if no window is given)
    # Using a simplified window selection for demonstration:
    if window is None:
        # Use a fixed fraction of the data for the time plots, e.g., 20%
        t1_idx = len(t_true) // 5 if len(t_true) > 100 else len(t_true)
    else:
        # Index corresponding to the window time
        t1_idx = int(window // dt)


    # Trim data for time-domain plots
    t_plot = t_true[t0_idx:t0_idx+t1_idx]
    y_raw_plot = y_raw[t0_idx:t0_idx+t1_idx]
    y_true_plot = y_true[t0_idx:t0_idx+t1_idx]
    bias_plot = b[t0_idx:t0_idx+t1_idx]
    if np.sum(abs(bias_plot)) < 1e-10:
        bias_plot = None
    else:
        if bias_plot.ndim == 1:
            bias_plot = bias_plot[:, np.newaxis]
        _, PSD_bias = fun_PSD(dt, b[t0_idx:nt_PSD + t0_idx])



    # X-limits for time plots
    xlim_time = [t_plot[0], t_plot[-1]]


    # 2. Figure Setup
    _, axes = plt.subplots(
        Nq, 5,
        figsize=(fig_width, 2. * Nq),
        layout='constrained',
        gridspec_kw={'width_ratios': [1, 1, 0.5, 1, 1], 'wspace': 0.1, 'hspace': 0.1}
    )

    # If Nq=1, axes will be a 1D array; ensure it's 2D for consistent indexing
    if Nq == 1:
        axes = axes.reshape(1, 5)

    titles = ['Raw', 'Truth', 'PDF', 'PSD', 'Difference']
    xlabels = ['$t$', '$t$', '$p$', '$f$', '$t$']


    c_raw = '#20b2aae5'
    c_true = '#000080ff'
    c_unbiased = "#8362caff"
    c_diff = "#6b256fff"
    c_bias = "#db76deff"


    # 3. Plotting Loop

    # Plotting column by column (more readable than the original's structure)
    for q_i in range(Nq):
        # Column 0: Raw Time Series (y_raw)
        ax = axes[q_i, 0]
        ax.plot(t_plot, y_raw_plot[:, q_i], color=c_raw, label=f'$y_{q_i}$')
        if y_obs is not None:
            ax.plot(t_obs, y_obs[:, q_i], 'ro', ms=3, mec='k', lw=.1)
        if y_wash is not None:
            ax.plot(t_wash, y_wash[:, q_i], 'rx', ms=3)

        ax.legend(fontsize='x-small', )
        ax.set(xlim=xlim_time)
        y_lim_base = ax.get_ylim()


        # Column 1: True Time Series (y_true)
        ax = axes[q_i, 1]
        ax.plot(t_plot, y_true_plot[:, q_i], color=c_true, label=f'$y^t_{q_i}$')
        if bias_plot is not None:
            ax.plot(t_plot, y_true_plot[:, q_i]-bias_plot[:, q_i], color=c_unbiased, label=f'$y^t_{q_i}-b^t_{q_i}$')
        # if q_i == 0:
        ax.legend(fontsize='x-small', ncol=2)

        y_lim_2 = ax.get_ylim()
        y_lim_base = [min(y_lim_base[0], y_lim_2[0]),
                      max(y_lim_base[1], y_lim_2[1])]
        # reset ax0 if changed
        axes[q_i, 0].set_ylim(y_lim_base)
        ax.set(xlim=xlim_time, ylim=y_lim_base)

        # Column 2: PDF (uses full data)
        ax = axes[q_i, 2]
        # Raw and true PDF
        for ds, c in zip([y_true, y_raw], [c_true, c_raw]):
            ax.hist(ds[:, q_i], bins=20, density=True, orientation='horizontal', color=c, histtype='stepfilled', alpha = 0.7)

        if bias_plot is not None:
            ax.hist(y_true_plot[:, q_i]-bias_plot[:, q_i], bins=20, density=True, orientation='horizontal', alpha = 0.7,
                    color=c_unbiased, label=f'$y^t_{q_i}-b^t_{q_i}$')

        if y_obs is not None:
            ax.hist(y_obs[:, q_i], bins=20, color='r', lw=1, histtype='step', density=True, orientation='horizontal')


        ax.set(ylim=y_lim_base)

        # Column 3: PSD (uses full data)
        ax = axes[q_i, 3]
        for ds, c, a in zip([PSD_true, PSD_raw], [c_true, c_raw], [1., .8]):
            ax.semilogy(f_raw, ds[q_i], color=c, alpha=a)
        if bias_plot is not None:
            ax.semilogy(f_raw, PSD_bias[q_i], color=c_unbiased, alpha=.8) #type: ignore

        if q_i == 0:
            ylims_PSD = [np.min(PSD_raw) * 0.1, np.max(PSD_raw) * 10]
        ax.set_xlim([0, f_max])
        ax.set_ylim(ylims_PSD)#type: ignore

        # Column 4: Difference Time Series (Noise)
        ax = axes[q_i, 4]

        noise = y_true_plot[:, q_i] - y_raw_plot[:, q_i]
        ax.plot(t_plot, noise, color=c_diff, label=f'$y^t - y_{q_i}$')
        ax.axhline(np.mean(noise), color='k', lw=.5, ls='--')

        if bias_plot is not None:
            bias = bias_plot[:, q_i]
            ax.plot(t_plot, bias, color=c_bias, label=f'$b^t_{q_i}$')

            # if q_i == 0:
        ax.legend(fontsize='x-small', ncol=2)

        ax.set(xlim=xlim_time)

        if q_i < Nq:
            for jj, ax in enumerate(axes[q_i, :]):
                if q_i != Nq-1:
                    ax.set_xticklabels([]) # No x-axis labels
                if jj in [1,2]:
                    ax.set_yticklabels([]) # No y-axis labels


    # Set titles and xlabels
    for i, (title, xlbl) in enumerate(zip(titles, xlabels)):
        axes[0, i].set_title(title)
        axes[-1, i].set_xlabel(xlbl)