Utils

Helper utilities: typed-dict conversion, saving/loading, AGATA analysis and the plotting functions.

Typed-dict helpers

Convert a mapping to a Numba typed dictionary with float64 values.

This helper normalizes a Python dict-like object into a numba.typed.Dict with keys of type unicode_type and values of type float64. It is useful when passing parameter dictionaries into Numba- compiled code or APIs that expect typed containers.

Parameters:

Name	Type	Description	Default
src	mapping	A mapping or mapping-like object to convert. If None, an empty typed dictionary is returned.	: None
strict	bool	If True, raise TypeError when a value cannot be cast to float64. If False, invalid entries are skipped silently.	: False

Returns:

Type	Description
Dict	A typed dictionary with str keys and float64 values.

Raises:

Type	Description
TypeError	If strict is True and a key/value pair cannot be converted to str/float64.

Source code in src/py_replay_bg/utils/numba_dicts.py

def to_typed_f64_dict(src=None, *, strict=False) -> Dict:
    """Convert a mapping to a Numba typed dictionary with float64 values.

    This helper normalizes a Python dict-like object into a
    ``numba.typed.Dict`` with keys of type ``unicode_type`` and values of type
    ``float64``. It is useful when passing parameter dictionaries into Numba-
    compiled code or APIs that expect typed containers.

    Parameters
    ----------
    src : mapping, optional, default : None
        A mapping or mapping-like object to convert. If ``None``, an empty typed
        dictionary is returned.
    strict : bool, optional, default : False
        If ``True``, raise ``TypeError`` when a value cannot be cast to
        ``float64``. If ``False``, invalid entries are skipped silently.

    Returns
    -------
    numba.typed.Dict
        A typed dictionary with ``str`` keys and ``float64`` values.

    Raises
    ------
    TypeError
        If ``strict`` is ``True`` and a key/value pair cannot be converted to
        ``str``/``float64``.
    """
    out = Dict.empty(key_type=types.unicode_type, value_type=float64)
    if src is None:
        return out

    for k, v in src.items():
        try:
            out[str(k)] = np.float64(v)
        except Exception:
            if strict:
                raise TypeError(f"Cannot cast key={k!r}, value={v!r} to float32")
            # non-strict: ignore invalid entries
    return out

Saving and loading results

Pickles results to <path>/<prefix>_<name>.pkl, creating path if needed.

Parameters:

Name	Type	Description	Default
results	dict	The dictionary of objects to persist.	required
path	str	Directory where the file will be written. Created if missing.	required
name	str or None	Base file name (with or without the .pkl extension). If None, defaults to today's date (YYYY_mm_dd).	required
prefix	str	Prefix always prepended to the file name as <prefix>_ (e.g. "twin" or "replay").	required

Returns:

Type	Description
str	The full path of the file written.

Source code in src/py_replay_bg/utils/save_results.py

def save_results(results: dict, path: str, name: str | None, prefix: str) -> str:
    """Pickles ``results`` to ``<path>/<prefix>_<name>.pkl``, creating ``path`` if needed.

    Parameters
    ----------
    results : dict
        The dictionary of objects to persist.
    path : str
        Directory where the file will be written. Created if missing.
    name : str or None
        Base file name (with or without the ``.pkl`` extension). If ``None``,
        defaults to today's date (``YYYY_mm_dd``).
    prefix : str
        Prefix always prepended to the file name as ``<prefix>_`` (e.g.
        ``"twin"`` or ``"replay"``).

    Returns
    -------
    str
        The full path of the file written.
    """
    os.makedirs(path, exist_ok=True)
    if name is None:
        name = f"{datetime.now().strftime('%Y_%m_%d')}"
    if not name.endswith('.pkl'):
        name += '.pkl'
    name = prefix + '_' + name
    file_path = os.path.join(path, name)
    with open(file_path, 'wb') as f:
        pickle.dump(results, f)
    return file_path

Loads results previously pickled by :func:save_results.

Parameters:

Name	Type	Description	Default
path	str	Directory where the file lives.	required
name	str	File name (with or without the .pkl extension), matching the name passed to :func:save_results.	required
prefix	str	Prefix used when the file was saved (e.g. "twin" or "replay").	required

Returns:

Type	Description
dict	The unpickled results dictionary.

Source code in src/py_replay_bg/utils/load_results.py

def load_results(path: str, name: str, prefix: str) -> dict:
    """Loads results previously pickled by :func:`save_results`.

    Parameters
    ----------
    path : str
        Directory where the file lives.
    name : str
        File name (with or without the ``.pkl`` extension), matching the ``name``
        passed to :func:`save_results`.
    prefix : str
        Prefix used when the file was saved (e.g. ``"twin"`` or ``"replay"``).

    Returns
    -------
    dict
        The unpickled results dictionary.
    """
    if not name.endswith('.pkl'):
        name += '.pkl'
    name = prefix + '_' + name
    file_path = os.path.join(path, name)
    with open(file_path, 'rb') as f:
        return pickle.load(f)

AGATA analysis

Analyze the outcome of a :meth:ReplayBG.twin run with AGATA.

Simulates the model forward with the estimated parameters to obtain the fitted (predicted) glucose trace, then:

• runs :meth:py_agata.Agata.analyze_glucose_profile on that trace, and
• computes the comprehensive py_agata.error metric set (rmse, mard, cod, Clarke EGA zones, g_rmse) against the reference CGM that was fitted (rbg_data.y).

Parameters:

Name	Type	Description	Default
twin_results	dict	The dict returned by :meth:ReplayBG.twin (must contain theta).	required
rbg_data	object	The data object passed to :meth:ReplayBG.twin. Provides the reference CGM (y), the input matrix (u), tsteps and the data sampling stride yts.	required
model	object	A model instance implementing reset(theta_dict, g0, u0), step(u, t) and output(t).	required
integration_ts	int	Minutes represented by one integration step (environment.ts). The AGATA cadence is rbg_data.yts * integration_ts.	: 1
glycemic_target	str	Glycemic target set forwarded to :class:py_agata.Agata.	: "diabetes"
verbose	bool	When True, prints a compact metric summary.	: False
path	str or None	When given, the analysis is added under the analysis key and the results (including rbg_data) are re-saved to <path>/twin_*.pkl, preserving the shape twin() persists.	: None
save_name	str or None	File name for the re-saved pickle.	: None

Returns:

Type	Description
dict	A dict with keys profile (AGATA profile of the fitted trace) and error (the fit-error metric set).

Source code in src/py_replay_bg/utils/agata_analysis.py

def analyze_twin(
    twin_results: dict,
    rbg_data,
    model,
    integration_ts: int = 1,
    glycemic_target: str = "diabetes",
    verbose: bool = False,
    path: str | None = None,
    save_name: str | None = None,
) -> dict:
    """Analyze the outcome of a :meth:`ReplayBG.twin` run with AGATA.

    Simulates the model forward with the estimated parameters to obtain the
    fitted (predicted) glucose trace, then:

    - runs :meth:`py_agata.Agata.analyze_glucose_profile` on that trace, and
    - computes the comprehensive ``py_agata.error`` metric set (``rmse``,
      ``mard``, ``cod``, Clarke EGA zones, ``g_rmse``) against the reference
      CGM that was fitted (``rbg_data.y``).

    Parameters
    ----------
    twin_results : dict
        The dict returned by :meth:`ReplayBG.twin` (must contain ``theta``).
    rbg_data : object
        The data object passed to :meth:`ReplayBG.twin`. Provides the reference
        CGM (``y``), the input matrix (``u``), ``tsteps`` and the data sampling
        stride ``yts``.
    model : object
        A model instance implementing ``reset(theta_dict, g0, u0)``,
        ``step(u, t)`` and ``output(t)``.
    integration_ts : int, optional, default : 1
        Minutes represented by one integration step (``environment.ts``). The
        AGATA cadence is ``rbg_data.yts * integration_ts``.
    glycemic_target : str, optional, default : "diabetes"
        Glycemic target set forwarded to :class:`py_agata.Agata`.
    verbose : bool, optional, default : False
        When ``True``, prints a compact metric summary.
    path : str or None, optional, default : None
        When given, the analysis is added under the ``analysis`` key and the
        results (including ``rbg_data``) are re-saved to ``<path>/twin_*.pkl``,
        preserving the shape ``twin()`` persists.
    save_name : str or None, optional, default : None
        File name for the re-saved pickle.

    Returns
    -------
    dict
        A dict with keys ``profile`` (AGATA profile of the fitted trace) and
        ``error`` (the fit-error metric set).
    """
    theta = twin_results["theta"]

    # Simulate the fitted model forward (same reset/step/output pass the twinner
    # uses for the likelihood; mirrors utils.plot_twinning). Re-seed the
    # cold-start glucose and the t=0 input row from data, same as the twinner
    # and plot_twinning, so this reproduces the exact trace that was fitted
    # (see Model.reset() for exactly when g0/u0 take effect).
    g0 = float(rbg_data.y[rbg_data.y_idxs[0]]) if len(rbg_data.y_idxs) > 0 else np.nan
    model.reset(to_typed_f64_dict(theta), g0, rbg_data.u[0])
    inputs = np.asarray(rbg_data.u)
    tsteps = int(rbg_data.tsteps)
    output = np.zeros(tsteps)
    output[0] = model.output(0)
    for k in range(1, tsteps):
        model.step(inputs[k], k)
        output[k] = model.output(k)

    yts = int(rbg_data.yts)
    pred = np.asarray(output[0::yts], dtype=float)
    ref = np.asarray(rbg_data.y, dtype=float)
    # Align both traces on the same homogeneous grid.
    n = min(pred.size, ref.size)
    pred, ref = pred[:n], ref[:n]

    sample_time = int(yts * integration_ts)
    ref_df = glucose_vector_to_dataframe(ref, sample_time)
    pred_df = glucose_vector_to_dataframe(pred, sample_time)

    agata = Agata(glycemic_target=glycemic_target)
    analysis = {
        "profile": agata.analyze_glucose_profile(pred_df),
        "error": {
            "rmse": rmse(ref_df, pred_df),
            "mard": mard(ref_df, pred_df),
            "cod": cod(ref_df, pred_df),
            "clarke": clarke(ref_df, pred_df),
            "g_rmse": g_rmse(ref_df, pred_df),
        },
    }

    if verbose:
        _print_analysis("Twin", analysis)

    if path is not None:
        save_results(
            {**twin_results, "rbg_data": rbg_data, "analysis": analysis},
            path, save_name, prefix="twin",
        )

    return analysis

Analyze the glucose traces of a :meth:ReplayBG.replay run with AGATA.

Runs :meth:py_agata.Agata.analyze_glucose_profile on the replayed model output and, when the replay used a sensor, also on the sensor measurement trace.

Parameters:

Name	Type	Description	Default
replay_results	dict	The dict returned by :meth:ReplayBG.replay (must contain output; measurement/measurement_time are used when present).	required
ts	int	AGATA glucose sampling cadence in minutes. The output trace is downsampled to this cadence before analysis.	: 5
integration_ts	int	Minutes represented by one integration step (environment.ts). Used to convert ts into a downsampling step and to scale the measurement cadence.	: 1
glycemic_target	str	Glycemic target set forwarded to :class:py_agata.Agata.	: "diabetes"
verbose	bool	When True, prints a compact metric summary.	: False
path	str or None	When given, the analysis is added to replay_results under the analysis key and the dict is re-saved to <path>/replay_*.pkl.	: None
save_name	str or None	File name for the re-saved pickle (see :func:utils.save_results.save_results).	: None

Returns:

Type	Description
dict	A dict with key output (the profile of the model output) and, when a sensor was used, measurement (the profile of the sensor trace).

Source code in src/py_replay_bg/utils/agata_analysis.py

def analyze_replay(
    replay_results: dict,
    ts: int = 5,
    integration_ts: int = 1,
    glycemic_target: str = "diabetes",
    verbose: bool = False,
    path: str | None = None,
    save_name: str | None = None,
) -> dict:
    """Analyze the glucose traces of a :meth:`ReplayBG.replay` run with AGATA.

    Runs :meth:`py_agata.Agata.analyze_glucose_profile` on the replayed model
    ``output`` and, when the replay used a sensor, also on the sensor
    ``measurement`` trace.

    Parameters
    ----------
    replay_results : dict
        The dict returned by :meth:`ReplayBG.replay` (must contain ``output``;
        ``measurement``/``measurement_time`` are used when present).
    ts : int, optional, default : 5
        AGATA glucose sampling cadence in minutes. The ``output`` trace is
        downsampled to this cadence before analysis.
    integration_ts : int, optional, default : 1
        Minutes represented by one integration step (``environment.ts``). Used to
        convert ``ts`` into a downsampling step and to scale the measurement
        cadence.
    glycemic_target : str, optional, default : "diabetes"
        Glycemic target set forwarded to :class:`py_agata.Agata`.
    verbose : bool, optional, default : False
        When ``True``, prints a compact metric summary.
    path : str or None, optional, default : None
        When given, the analysis is added to ``replay_results`` under the
        ``analysis`` key and the dict is re-saved to ``<path>/replay_*.pkl``.
    save_name : str or None, optional, default : None
        File name for the re-saved pickle (see
        :func:`utils.save_results.save_results`).

    Returns
    -------
    dict
        A dict with key ``output`` (the profile of the model output) and, when a
        sensor was used, ``measurement`` (the profile of the sensor trace).
    """
    agata = Agata(glycemic_target=glycemic_target)

    # The output lives at integration resolution; subsample it to the AGATA
    # cadence so the DataFrame timestamps are honest.
    step = max(1, int(round(ts / integration_ts)))
    out = np.asarray(replay_results["output"][0::step], dtype=float)
    out_df = glucose_vector_to_dataframe(out, int(ts))
    analysis = {"output": agata.analyze_glucose_profile(out_df)}

    # Sensor measurements are already at sensor cadence; recover their real
    # spacing from the integration-step indices of consecutive samples.
    if "measurement" in replay_results:
        meas = np.asarray(replay_results["measurement"], dtype=float)
        mt = np.asarray(replay_results.get("measurement_time", []))
        meas_ts = int((mt[1] - mt[0]) * integration_ts) if mt.size >= 2 else int(ts)
        meas_df = glucose_vector_to_dataframe(meas, meas_ts)
        analysis["measurement"] = agata.analyze_glucose_profile(meas_df)

    if verbose:
        _print_analysis("Replay", analysis)

    if path is not None:
        replay_results["analysis"] = analysis
        save_results(replay_results, path, save_name, prefix="replay")

    return analysis

Plotting

Plot the outcome of a ReplayBG.twin() run.

Twinning fits the model parameters to the observed data; this utility shows how well the fitted model reproduces them. The model is simulated forward with the estimated theta (the same reset/step/output pass the twinner uses for the likelihood), and the resulting trace is compared against the observations that were fitted.

The figure follows the same rationale as :func:plot_replay: a vertical stack of subplots sharing a common time axis.

1. Fit (first subplot): the simulated model output for the estimated theta as a continuous line and the observed data (rbg_data.y at the non-missing indices rbg_data.y_idxs) as markers. Each value in thresholds is drawn as a dashed horizontal line spanning the whole width, acting as a visual reference level, and the span between the outermost levels is shaded.
2. Inputs (one subplot per input group): the inputs that drove the fit, labelled via data_to_input. By default every channel gets its own subplot; pass input_groups to overlay channels together.

The utility is model-agnostic: nothing about the meaning of the output or inputs is assumed. Channel names come straight from data_to_input.

Parameters:

Name	Type	Description	Default
rbg_data	object	The data object passed to :meth:ReplayBG.twin (e.g. a MultiMealT1DData). Must expose u, data_to_input, tsteps, yts, y and y_idxs.	required
model	object	A model instance implementing the model interface contract (reset(theta0, g0, u0), step(u, t), output(t)).	required
theta	dict or None	The estimated parameters as returned in twin()['theta']. When given, the model is reset with them before simulating. When None, the model is reset with its own theta0 (e.g. already constructed with the estimated parameters).	: None
thresholds	list of float or None	Optional reference levels drawn as dashed horizontal lines on the fit subplot. Two or more levels also shade the span between the outermost.	: None
input_groups	list of list of int or None	Optional list of channel-index groups. Each group becomes one subplot with its channels overlaid. None plots one channel per subplot, in data_to_input order.	: None
mask_inputs	list of int or None	Optional list of channel indices to hide. Masked channels are dropped from the (default or explicit) input_groups so they get no subplot; e.g. pass the t_hour index to skip it.	: None
ts_min	float	Minutes represented by one integration step, used to scale the time axis.	: 1.0
output_label	str	Legend label for the simulated fit line.	: "Fit"
observation_label	str	Legend label for the observed-data markers.	: "Data"
figsize	tuple of float or None	Optional figure size. Defaults to a height that grows with the number of subplots.	: None
hover	bool	Whether to attach the interactive readout: hovering any subplot drops a vertical cursor across the whole stack and reports the values at that minute. Ignored on non-interactive backends.	: True

Returns:

Type	Description
Figure	The assembled figure.

Source code in src/py_replay_bg/utils/plot_twinning.py

def plot_twinning(
    rbg_data,
    model,
    theta: dict | None = None,
    thresholds: list[float] | None = None,
    input_groups: list[list[int]] | None = None,
    mask_inputs: list[int] | None = None,
    ts_min: float = 1.0,
    output_label: str = "Fit",
    observation_label: str = "Data",
    figsize: tuple[float, float] | None = None,
    hover: bool = True,
) -> plt.Figure:
    """Plot the outcome of a ``ReplayBG.twin()`` run.

    Twinning fits the model parameters to the observed data; this utility shows
    how well the fitted model reproduces them. The model is simulated forward
    with the estimated ``theta`` (the same ``reset``/``step``/``output`` pass the
    twinner uses for the likelihood), and the resulting trace is compared against
    the observations that were fitted.

    The figure follows the same rationale as :func:`plot_replay`: a vertical
    stack of subplots sharing a common time axis.

    1. **Fit** (first subplot): the simulated model ``output`` for the estimated
       ``theta`` as a continuous line and the observed data (``rbg_data.y`` at the
       non-missing indices ``rbg_data.y_idxs``) as markers. Each value in
       ``thresholds`` is drawn as a dashed horizontal line spanning the whole
       width, acting as a visual reference level, and the span between the
       outermost levels is shaded.
    2. **Inputs** (one subplot per input *group*): the inputs that drove the
       fit, labelled via ``data_to_input``. By default every channel gets its own
       subplot; pass ``input_groups`` to overlay channels together.

    The utility is model-agnostic: nothing about the meaning of the output or
    inputs is assumed. Channel names come straight from ``data_to_input``.

    Parameters
    ----------
    rbg_data : object
        The data object passed to :meth:`ReplayBG.twin` (e.g. a
        ``MultiMealT1DData``). Must expose ``u``, ``data_to_input``, ``tsteps``,
        ``yts``, ``y`` and ``y_idxs``.
    model : object
        A model instance implementing the model interface contract
        (``reset(theta0, g0, u0)``, ``step(u, t)``, ``output(t)``).
    theta : dict or None, optional, default : None
        The estimated parameters as returned in ``twin()['theta']``. When given,
        the model is reset with them before simulating. When ``None``, the model
        is reset with its own ``theta0`` (e.g. already constructed with the
        estimated parameters).
    thresholds : list of float or None, optional, default : None
        Optional reference levels drawn as dashed horizontal lines on the fit
        subplot. Two or more levels also shade the span between the outermost.
    input_groups : list of list of int or None, optional, default : None
        Optional list of channel-index groups. Each group becomes one subplot
        with its channels overlaid. ``None`` plots one channel per subplot, in
        ``data_to_input`` order.
    mask_inputs : list of int or None, optional, default : None
        Optional list of channel indices to hide. Masked channels are dropped
        from the (default or explicit) ``input_groups`` so they get no subplot;
        e.g. pass the ``t_hour`` index to skip it.
    ts_min : float, optional, default : 1.0
        Minutes represented by one integration step, used to scale the time axis.
    output_label : str, optional, default : "Fit"
        Legend label for the simulated fit line.
    observation_label : str, optional, default : "Data"
        Legend label for the observed-data markers.
    figsize : tuple of float or None, optional, default : None
        Optional figure size. Defaults to a height that grows with the number of
        subplots.
    hover : bool, optional, default : True
        Whether to attach the interactive readout: hovering any subplot drops a
        vertical cursor across the whole stack and reports the values at that
        minute. Ignored on non-interactive backends.

    Returns
    -------
    matplotlib.figure.Figure
        The assembled figure.
    """
    inputs = np.asarray(rbg_data.u)
    data_to_input = rbg_data.data_to_input

    # --- simulate the fitted model forward (matches the twinner's pass) -------
    output = _simulate(model, rbg_data, theta)
    tsteps = output.shape[0]

    t = np.arange(tsteps) * ts_min

    # --- resolve the subplot layout ------------------------------------------
    input_groups, input_colors = resolve_input_layout(data_to_input, input_groups, mask_inputs)

    n_rows = 1 + len(input_groups)
    fig, axes = make_figure(n_rows, figsize)
    hover_series = {}

    # --- 1. fit subplot ------------------------------------------------------
    ax = axes[0]
    draw_thresholds(ax, thresholds)
    ax.plot(t, output, color="tab:blue", linewidth=1.2, label=output_label)

    # Observations live at data resolution; each sample i sits yts integration
    # steps apart. Only the non-missing samples (y_idxs) were used in the fit.
    y = np.asarray(rbg_data.y, dtype=float)
    y_idxs = np.asarray(rbg_data.y_idxs)
    obs_t = y_idxs * rbg_data.yts * ts_min
    ax.plot(
        obs_t, y[y_idxs],
        linestyle="none", marker="o", markersize=3,
        color="tab:red", alpha=0.7, label=observation_label,
    )
    hover_series[ax] = [
        series(output_label, t, output, KIND_DENSE),
        series(observation_label, obs_t, y[y_idxs], KIND_SPARSE),
    ]

    ax.set_ylabel(output_label)
    set_panel_title(ax, "Twinning fit")
    ax.legend(fontsize="small", loc="upper right", framealpha=0.9)

    # --- 2. input subplots ---------------------------------------------------
    for row, group in enumerate(input_groups, start=1):
        hover_series[axes[row]] = draw_input_group(
            axes[row], t, inputs, group, data_to_input, input_colors)

    finalize(fig, axes, "Time (min)", hover_series, hover, x_fmt="t = {:.0f} min")
    return fig

Plot several x0-chained twin() segments as one continuous figure.

The "interval" workflow fits one day, carries its final state into the next (via x0/theta_prev), and repeats. Each day is a separate twin over its own rbg_data and fitted theta; this utility glues them back together — laying the segments end to end on a single time axis — so the multi-day fit reads as one trace. It is the multi-segment counterpart of :func:plot_twinning and produces the same panel layout:

1. Fit (first subplot): every segment's simulated output drawn as one continuous line, with the observed data of all segments as markers.
2. Inputs (one subplot per input group): the glued input channels.

A dotted vertical separator marks each join between consecutive days.

Each segment must expose the same input-channel layout (data_to_input); the layout and the observation cadence (yts) are taken from the first.

Parameters:

Name	Type	Description	Default
segments	list of dict	The per-day segments in chronological order. Each is a mapping with keys rbg_data (the day's prepared data), model (a model instance) and optional theta (the day's estimated parameters). When theta is given the model is reset with it before simulating; otherwise the model is reset with its own theta0.	required
thresholds	list of float or None	Optional reference levels drawn as dashed horizontal lines on the fit subplot. Two or more levels also shade the span between the outermost.	: None
input_groups	list of list of int or None	Optional list of channel-index groups. Each group becomes one subplot with its channels overlaid. None plots one channel per subplot, in data_to_input order.	: None
mask_inputs	list of int or None	Optional list of channel indices to hide, e.g. the t_hour index.	: None
ts_min	float	Minutes represented by one integration step, used to scale the time axis.	: 1.0
output_label	str	Legend label for the simulated fit line.	: "Fit"
observation_label	str	Legend label for the observed-data markers.	: "Data"
figsize	tuple of float or None	Optional figure size. Defaults to a height that grows with the number of subplots.	: None
hover	bool	Whether to attach the interactive readout. Ignored on non-interactive backends.	: True

Returns:

Type	Description
Figure	The assembled figure.

Source code in src/py_replay_bg/utils/plot_twinning.py

def plot_twinning_intervals(
    segments: list[dict],
    thresholds: list[float] | None = None,
    input_groups: list[list[int]] | None = None,
    mask_inputs: list[int] | None = None,
    ts_min: float = 1.0,
    output_label: str = "Fit",
    observation_label: str = "Data",
    figsize: tuple[float, float] | None = None,
    hover: bool = True,
) -> plt.Figure:
    """Plot several x0-chained ``twin()`` segments as one continuous figure.

    The "interval" workflow fits one day, carries its final state into the next
    (via ``x0``/``theta_prev``), and repeats. Each day is a separate twin over
    its own ``rbg_data`` and fitted ``theta``; this utility glues them back
    together — laying the segments end to end on a single time axis — so the
    multi-day fit reads as one trace. It is the multi-segment counterpart of
    :func:`plot_twinning` and produces the same panel layout:

    1. **Fit** (first subplot): every segment's simulated ``output`` drawn as one
       continuous line, with the observed data of all segments as markers.
    2. **Inputs** (one subplot per input *group*): the glued input channels.

    A dotted vertical separator marks each join between consecutive days.

    Each segment must expose the same input-channel layout (``data_to_input``);
    the layout and the observation cadence (``yts``) are taken from the first.

    Parameters
    ----------
    segments : list of dict
        The per-day segments in chronological order. Each is a mapping with keys
        ``rbg_data`` (the day's prepared data), ``model`` (a model instance) and
        optional ``theta`` (the day's estimated parameters). When ``theta`` is
        given the model is reset with it before simulating; otherwise the model
        is reset with its own ``theta0``.
    thresholds : list of float or None, optional, default : None
        Optional reference levels drawn as dashed horizontal lines on the fit
        subplot. Two or more levels also shade the span between the outermost.
    input_groups : list of list of int or None, optional, default : None
        Optional list of channel-index groups. Each group becomes one subplot
        with its channels overlaid. ``None`` plots one channel per subplot, in
        ``data_to_input`` order.
    mask_inputs : list of int or None, optional, default : None
        Optional list of channel indices to hide, e.g. the ``t_hour`` index.
    ts_min : float, optional, default : 1.0
        Minutes represented by one integration step, used to scale the time axis.
    output_label : str, optional, default : "Fit"
        Legend label for the simulated fit line.
    observation_label : str, optional, default : "Data"
        Legend label for the observed-data markers.
    figsize : tuple of float or None, optional, default : None
        Optional figure size. Defaults to a height that grows with the number of
        subplots.
    hover : bool, optional, default : True
        Whether to attach the interactive readout. Ignored on non-interactive
        backends.

    Returns
    -------
    matplotlib.figure.Figure
        The assembled figure.
    """
    if not segments:
        raise ValueError("segments must contain at least one segment")

    data_to_input = segments[0]["rbg_data"].data_to_input

    # --- simulate every segment and glue the traces end to end ---------------
    outputs, inputs_list, tsteps_list = [], [], []
    obs_t_parts, obs_y_parts = [], []
    tsteps_cum = [int(seg["rbg_data"].tsteps) for seg in segments]
    starts, boundaries = segment_offsets(tsteps_cum, ts_min)

    for seg, start in zip(segments, starts):
        rbg_data = seg["rbg_data"]
        output = _simulate(seg["model"], rbg_data, seg.get("theta"))
        outputs.append(output)
        inputs_list.append(np.asarray(rbg_data.u))
        tsteps_list.append(output.shape[0])

        # Observations live at data resolution; only the non-missing samples
        # (y_idxs) were fitted. Offset each segment onto the glued axis.
        y = np.asarray(rbg_data.y, dtype=float)
        y_idxs = np.asarray(rbg_data.y_idxs)
        obs_t_parts.append(start + y_idxs * rbg_data.yts * ts_min)
        obs_y_parts.append(y[y_idxs])

    output = np.concatenate(outputs)
    inputs = np.concatenate(inputs_list, axis=0)
    t = np.concatenate([start + np.arange(n) * ts_min
                        for start, n in zip(starts, tsteps_list)])
    obs_t = np.concatenate(obs_t_parts)
    obs_y = np.concatenate(obs_y_parts)

    # --- resolve the subplot layout ------------------------------------------
    input_groups, input_colors = resolve_input_layout(data_to_input, input_groups, mask_inputs)

    n_rows = 1 + len(input_groups)
    fig, axes = make_figure(n_rows, figsize)
    hover_series = {}

    # --- 1. fit subplot ------------------------------------------------------
    ax = axes[0]
    draw_thresholds(ax, thresholds)
    ax.plot(t, output, color="tab:blue", linewidth=1.2, label=output_label)
    ax.plot(
        obs_t, obs_y,
        linestyle="none", marker="o", markersize=3,
        color="tab:red", alpha=0.7, label=observation_label,
    )
    draw_segment_boundaries(ax, boundaries)
    hover_series[ax] = [
        series(output_label, t, output, KIND_DENSE),
        series(observation_label, obs_t, obs_y, KIND_SPARSE),
    ]

    ax.set_ylabel(output_label)
    set_panel_title(ax, f"Twinning fit ({len(segments)} segments)")
    ax.legend(fontsize="small", loc="upper right", framealpha=0.9)

    # --- 2. input subplots ---------------------------------------------------
    for row, group in enumerate(input_groups, start=1):
        hover_series[axes[row]] = draw_input_group(
            axes[row], t, inputs, group, data_to_input, input_colors)
        draw_segment_boundaries(axes[row], boundaries)

    finalize(fig, axes, "Time (min)", hover_series, hover, x_fmt="t = {:.0f} min")
    return fig

Plot the outcome of a ReplayBG.replay() run.

The figure is a vertical stack of subplots sharing a common time axis:

1. Output (first subplot): the replayed output as a continuous line and, when a sensor was used, the replayed measurement as markers. Each value in thresholds is drawn as a dashed horizontal line spanning the whole width, acting as a visual reference level, and the span between the outermost levels is shaded.
2. Inputs (one subplot per input group): the applied input channels, labelled via data_to_input. By default every channel gets its own subplot; pass input_groups to overlay channels together.
3. Callback actions (one subplot per distinct callback): the values a callback logged, stemmed at the integration minute of each action.

The utility is model-agnostic: nothing about the meaning of the output, inputs, or actions is assumed. Channel names come straight from data_to_input and callback labels from the action log.

Parameters:

Name	Type	Description	Default
replay_results	dict	The dict returned by :meth:ReplayBG.replay. Must contain output, input, data_to_input and actions; measurement and measurement_time are used when present.	required
thresholds	list of float or None	Optional reference levels drawn as dashed horizontal lines on the output subplot. Two or more levels also shade the span between the outermost.	: None
input_groups	list of list of int or None	Optional list of channel-index groups. Each group becomes one subplot with its channels overlaid. None plots one channel per subplot, in data_to_input order.	: None
mask_inputs	list of int or None	Optional list of channel indices to hide. Masked channels are dropped from the (default or explicit) input_groups so they get no subplot; e.g. pass the t_hour index to skip it.	: None
action_field	str or None	Name of the scalar field to plot for each callback action. None resolves it per callback to the first numeric field the callback logged (anything other than k/callback), which by convention is the quantity the action delivered. An explicit name is used for every callback that logs it; callbacks that do not fall back to the automatic choice.	: None
ts_min	float	Minutes represented by one integration step, used to scale the time axis.	: 1.0
output_label	str	Legend label for the replayed output line.	: "Output"
measurement_label	str	Legend label for the replayed measurement markers.	: "Measurement"
figsize	tuple of float or None	Optional figure size. Defaults to a height that grows with the number of subplots.	: None
hover	bool	Whether to attach the interactive readout: hovering any subplot drops a vertical cursor across the whole stack and reports the values at that minute. Ignored on non-interactive backends.	: True

Returns:

Type	Description
Figure	The assembled figure.

Source code in src/py_replay_bg/utils/plot_replay.py

def plot_replay(
    replay_results: dict,
    thresholds: list[float] | None = None,
    input_groups: list[list[int]] | None = None,
    mask_inputs: list[int] | None = None,
    action_field: str | None = None,
    ts_min: float = 1.0,
    output_label: str = "Output",
    measurement_label: str = "Measurement",
    figsize: tuple[float, float] | None = None,
    hover: bool = True,
) -> plt.Figure:
    """Plot the outcome of a ``ReplayBG.replay()`` run.

    The figure is a vertical stack of subplots sharing a common time axis:

    1. **Output** (first subplot): the replayed ``output`` as a continuous line
       and, when a sensor was used, the replayed ``measurement`` as markers.
       Each value in ``thresholds`` is drawn as a dashed horizontal line
       spanning the whole width, acting as a visual reference level, and the span
       between the outermost levels is shaded.
    2. **Inputs** (one subplot per input *group*): the applied ``input``
       channels, labelled via ``data_to_input``. By default every channel gets
       its own subplot; pass ``input_groups`` to overlay channels together.
    3. **Callback actions** (one subplot per distinct callback): the values a
       callback logged, stemmed at the integration minute of each action.

    The utility is model-agnostic: nothing about the meaning of the output,
    inputs, or actions is assumed. Channel names come straight from
    ``data_to_input`` and callback labels from the action log.

    Parameters
    ----------
    replay_results : dict
        The dict returned by :meth:`ReplayBG.replay`. Must contain ``output``,
        ``input``, ``data_to_input`` and ``actions``; ``measurement`` and
        ``measurement_time`` are used when present.
    thresholds : list of float or None, optional, default : None
        Optional reference levels drawn as dashed horizontal lines on the output
        subplot. Two or more levels also shade the span between the outermost.
    input_groups : list of list of int or None, optional, default : None
        Optional list of channel-index groups. Each group becomes one subplot
        with its channels overlaid. ``None`` plots one channel per subplot, in
        ``data_to_input`` order.
    mask_inputs : list of int or None, optional, default : None
        Optional list of channel indices to hide. Masked channels are dropped
        from the (default or explicit) ``input_groups`` so they get no subplot;
        e.g. pass the ``t_hour`` index to skip it.
    action_field : str or None, optional, default : None
        Name of the scalar field to plot for each callback action. ``None``
        resolves it per callback to the first numeric field the callback logged
        (anything other than ``k``/``callback``), which by convention is the
        quantity the action delivered. An explicit name is used for every
        callback that logs it; callbacks that do not fall back to the automatic
        choice.
    ts_min : float, optional, default : 1.0
        Minutes represented by one integration step, used to scale the time axis.
    output_label : str, optional, default : "Output"
        Legend label for the replayed output line.
    measurement_label : str, optional, default : "Measurement"
        Legend label for the replayed measurement markers.
    figsize : tuple of float or None, optional, default : None
        Optional figure size. Defaults to a height that grows with the number of
        subplots.
    hover : bool, optional, default : True
        Whether to attach the interactive readout: hovering any subplot drops a
        vertical cursor across the whole stack and reports the values at that
        minute. Ignored on non-interactive backends.

    Returns
    -------
    matplotlib.figure.Figure
        The assembled figure.
    """
    output = np.asarray(replay_results["output"])
    inputs = np.asarray(replay_results["input"])
    data_to_input = replay_results["data_to_input"]
    actions = replay_results.get("actions", []) or []

    tsteps = output.shape[0]
    t = np.arange(tsteps) * ts_min

    # --- resolve the subplot layout -----------------------------------------
    input_groups, input_colors = resolve_input_layout(data_to_input, input_groups, mask_inputs)

    # Distinct callbacks, in first-seen order, that actually logged something.
    callback_names = list(dict.fromkeys(rec["callback"] for rec in actions))

    n_rows = 1 + len(input_groups) + len(callback_names)
    fig, axes = make_figure(n_rows, figsize)
    hover_series = {}

    # --- 1. output subplot ---------------------------------------------------
    ax = axes[0]
    draw_thresholds(ax, thresholds)
    ax.plot(t, output, color="tab:blue", linewidth=1.2, label=output_label)
    hover_series[ax] = [series(output_label, t, output, KIND_DENSE)]

    measurement = replay_results.get("measurement")
    measurement_time = replay_results.get("measurement_time")
    if measurement is not None and measurement_time is not None:
        measurement = np.asarray(measurement)
        measurement_time = np.asarray(measurement_time)
        mask = np.isfinite(measurement)
        meas_t = measurement_time[mask] * ts_min
        ax.plot(
            meas_t, measurement[mask],
            linestyle="none", marker="o", markersize=3,
            color="tab:red", alpha=0.7, label=measurement_label,
        )
        hover_series[ax].append(series(measurement_label, meas_t, measurement[mask], KIND_SPARSE))

    ax.set_ylabel(output_label)
    set_panel_title(ax, "Output")
    ax.legend(fontsize="small", loc="upper right", framealpha=0.9)

    # --- 2. input subplots ---------------------------------------------------
    for row, group in enumerate(input_groups, start=1):
        hover_series[axes[row]] = draw_input_group(
            axes[row], t, inputs, group, data_to_input, input_colors)

    # --- 3. callback action subplots ----------------------------------------
    base = 1 + len(input_groups)
    for row, cb_name in enumerate(callback_names, start=base):
        ax = axes[row]
        records = [rec for rec in actions if rec["callback"] == cb_name]
        field = _resolve_action_field(records, action_field)
        ks = np.array([rec["k"] * ts_min for rec in records if field in rec])
        values = np.array([rec[field] for rec in records if field in rec], dtype=float)
        # Registered even when empty, so the cursor still tracks across a panel
        # whose callback logged nothing plottable.
        hover_series[ax] = []
        if ks.size:
            markerline, _, _ = ax.stem(ks, values, linefmt="tab:purple",
                                       markerfmt="o", basefmt=" ")
            plt.setp(markerline, color="tab:purple", markersize=4)
            hover_series[ax] = [series(field, ks, values, KIND_SPARSE)]
        ax.set_ylabel(field or "")
        set_panel_title(ax, f"Callback: {cb_name}")

    finalize(fig, axes, "Time (min)", hover_series, hover, x_fmt="t = {:.0f} min")
    return fig

Plot several x0-chained replay() segments as one continuous figure.

The "interval" workflow replays one day, carries its final state into the next (via x0/theta_prev), and repeats. Each day is a separate replay; this utility glues them back together — laying the segments end to end on a single time axis — so the multi-day replay reads as one trace. It is the multi-segment counterpart of :func:plot_replay and produces the same panel layout (output, inputs, one panel per distinct callback), with a dotted vertical separator marking each join between consecutive days.

Each segment must expose the same input-channel layout (data_to_input), taken from the first. Measurement markers are shown only when every segment carries a sensor measurement.

Parameters:

Name	Type	Description	Default
replay_results_list	list of dict	The per-day dicts returned by :meth:ReplayBG.replay, in chronological order. Each must contain output, input, data_to_input and actions; measurement/measurement_time are used when present on all segments.	required
thresholds	list of float or None	Optional reference levels drawn as dashed horizontal lines on the output subplot. Two or more levels also shade the span between the outermost.	: None
input_groups	list of list of int or None	Optional list of channel-index groups, each becoming one overlaid subplot. None plots one channel per subplot, in data_to_input order.	: None
mask_inputs	list of int or None	Optional list of channel indices to hide, e.g. the t_hour index.	: None
action_field	str or None	Name of the scalar field to plot for each callback action, as in :func:plot_replay. None resolves it per callback automatically.	: None
ts_min	float	Minutes represented by one integration step, used to scale the time axis.	: 1.0
output_label	str	Legend label for the replayed output line.	: "Output"
measurement_label	str	Legend label for the replayed measurement markers.	: "Measurement"
figsize	tuple of float or None	Optional figure size. Defaults to a height that grows with the number of subplots.	: None
hover	bool	Whether to attach the interactive readout. Ignored on non-interactive backends.	: True

Returns:

Type	Description
Figure	The assembled figure.

Source code in src/py_replay_bg/utils/plot_replay.py

def plot_replay_intervals(
    replay_results_list: list[dict],
    thresholds: list[float] | None = None,
    input_groups: list[list[int]] | None = None,
    mask_inputs: list[int] | None = None,
    action_field: str | None = None,
    ts_min: float = 1.0,
    output_label: str = "Output",
    measurement_label: str = "Measurement",
    figsize: tuple[float, float] | None = None,
    hover: bool = True,
) -> plt.Figure:
    """Plot several x0-chained ``replay()`` segments as one continuous figure.

    The "interval" workflow replays one day, carries its final state into the
    next (via ``x0``/``theta_prev``), and repeats. Each day is a separate replay;
    this utility glues them back together — laying the segments end to end on a
    single time axis — so the multi-day replay reads as one trace. It is the
    multi-segment counterpart of :func:`plot_replay` and produces the same panel
    layout (output, inputs, one panel per distinct callback), with a dotted
    vertical separator marking each join between consecutive days.

    Each segment must expose the same input-channel layout (``data_to_input``),
    taken from the first. Measurement markers are shown only when *every* segment
    carries a sensor measurement.

    Parameters
    ----------
    replay_results_list : list of dict
        The per-day dicts returned by :meth:`ReplayBG.replay`, in chronological
        order. Each must contain ``output``, ``input``, ``data_to_input`` and
        ``actions``; ``measurement``/``measurement_time`` are used when present
        on all segments.
    thresholds : list of float or None, optional, default : None
        Optional reference levels drawn as dashed horizontal lines on the output
        subplot. Two or more levels also shade the span between the outermost.
    input_groups : list of list of int or None, optional, default : None
        Optional list of channel-index groups, each becoming one overlaid
        subplot. ``None`` plots one channel per subplot, in ``data_to_input``
        order.
    mask_inputs : list of int or None, optional, default : None
        Optional list of channel indices to hide, e.g. the ``t_hour`` index.
    action_field : str or None, optional, default : None
        Name of the scalar field to plot for each callback action, as in
        :func:`plot_replay`. ``None`` resolves it per callback automatically.
    ts_min : float, optional, default : 1.0
        Minutes represented by one integration step, used to scale the time axis.
    output_label : str, optional, default : "Output"
        Legend label for the replayed output line.
    measurement_label : str, optional, default : "Measurement"
        Legend label for the replayed measurement markers.
    figsize : tuple of float or None, optional, default : None
        Optional figure size. Defaults to a height that grows with the number of
        subplots.
    hover : bool, optional, default : True
        Whether to attach the interactive readout. Ignored on non-interactive
        backends.

    Returns
    -------
    matplotlib.figure.Figure
        The assembled figure.
    """
    if not replay_results_list:
        raise ValueError("replay_results_list must contain at least one segment")

    data_to_input = replay_results_list[0]["data_to_input"]

    # --- glue outputs and inputs end to end ----------------------------------
    outputs = [np.asarray(r["output"]) for r in replay_results_list]
    inputs_list = [np.asarray(r["input"]) for r in replay_results_list]
    tsteps_list = [o.shape[0] for o in outputs]
    starts, boundaries = segment_offsets(tsteps_list, ts_min)
    step_starts = [int(round(s / ts_min)) for s in starts]

    output = np.concatenate(outputs)
    inputs = np.concatenate(inputs_list, axis=0)
    t = np.concatenate([start + np.arange(n) * ts_min
                        for start, n in zip(starts, tsteps_list)])

    # Actions carry a step index local to their segment; offset each onto the
    # glued axis so callback panels line up with the trace.
    actions = []
    for r, step_start in zip(replay_results_list, step_starts):
        for rec in (r.get("actions") or []):
            actions.append({**rec, "k": rec["k"] + step_start})

    # Measurements are shown only when every segment produced them.
    have_meas = all(r.get("measurement") is not None
                    and r.get("measurement_time") is not None
                    for r in replay_results_list)
    meas_t, meas_v = None, None
    if have_meas:
        meas_t_parts, meas_v_parts = [], []
        for r, start in zip(replay_results_list, starts):
            m = np.asarray(r["measurement"])
            mt = np.asarray(r["measurement_time"])
            valid = np.isfinite(m)
            meas_t_parts.append(start + mt[valid] * ts_min)
            meas_v_parts.append(m[valid])
        meas_t = np.concatenate(meas_t_parts)
        meas_v = np.concatenate(meas_v_parts)

    # --- resolve the subplot layout ------------------------------------------
    input_groups, input_colors = resolve_input_layout(data_to_input, input_groups, mask_inputs)
    callback_names = list(dict.fromkeys(rec["callback"] for rec in actions))

    n_rows = 1 + len(input_groups) + len(callback_names)
    fig, axes = make_figure(n_rows, figsize)
    hover_series = {}

    # --- 1. output subplot ---------------------------------------------------
    ax = axes[0]
    draw_thresholds(ax, thresholds)
    ax.plot(t, output, color="tab:blue", linewidth=1.2, label=output_label)
    hover_series[ax] = [series(output_label, t, output, KIND_DENSE)]
    if have_meas:
        ax.plot(
            meas_t, meas_v,
            linestyle="none", marker="o", markersize=3,
            color="tab:red", alpha=0.7, label=measurement_label,
        )
        hover_series[ax].append(series(measurement_label, meas_t, meas_v, KIND_SPARSE))
    draw_segment_boundaries(ax, boundaries)
    ax.set_ylabel(output_label)
    set_panel_title(ax, f"Output ({len(replay_results_list)} segments)")
    ax.legend(fontsize="small", loc="upper right", framealpha=0.9)

    # --- 2. input subplots ---------------------------------------------------
    for row, group in enumerate(input_groups, start=1):
        hover_series[axes[row]] = draw_input_group(
            axes[row], t, inputs, group, data_to_input, input_colors)
        draw_segment_boundaries(axes[row], boundaries)

    # --- 3. callback action subplots -----------------------------------------
    base = 1 + len(input_groups)
    for row, cb_name in enumerate(callback_names, start=base):
        ax = axes[row]
        records = [rec for rec in actions if rec["callback"] == cb_name]
        field = _resolve_action_field(records, action_field)
        ks = np.array([rec["k"] * ts_min for rec in records if field in rec])
        values = np.array([rec[field] for rec in records if field in rec], dtype=float)
        hover_series[ax] = []
        if ks.size:
            markerline, _, _ = ax.stem(ks, values, linefmt="tab:purple",
                                       markerfmt="o", basefmt=" ")
            plt.setp(markerline, color="tab:purple", markersize=4)
            hover_series[ax] = [series(field, ks, values, KIND_SPARSE)]
        draw_segment_boundaries(ax, boundaries)
        ax.set_ylabel(field or "")
        set_panel_title(ax, f"Callback: {cb_name}")

    finalize(fig, axes, "Time (min)", hover_series, hover, x_fmt="t = {:.0f} min")
    return fig

Plots the optimization history from a twinning run.

The figure follows the same rationale as :func:plot_replay and :func:plot_twinning: a vertical stack of subplots sharing a common x axis, here the function evaluation rather than time.

1. Parameter trajectory: each unknown parameter's value per evaluation.
2. Log-prior, log-likelihood, log-posterior: the three terms the twinner is optimising, one per subplot.

Parameters:

Name	Type	Description	Default
history	dict	Dict with keys theta, log_prior, log_likelihood and log_posterior. Each value is a list of recorded values per function evaluation. theta entries are numpy arrays.	required
param_names	list of str or None	Optional list of parameter names for the theta legend. Defaults to theta_0, theta_1, ...	: None
figsize	tuple of float or None	Optional figure size. Defaults to a height that grows with the number of subplots.	: None
hover	bool	Whether to attach the interactive readout: hovering any subplot drops a vertical cursor across the whole stack and reports the values at that evaluation. Ignored on non-interactive backends.	: True

Returns:

Type	Description
Figure	The assembled figure.

Source code in src/py_replay_bg/utils/plot_twinning_history.py

def plot_twinning_history(
    history: dict,
    param_names: list[str] | None = None,
    figsize: tuple[float, float] | None = None,
    hover: bool = True,
) -> plt.Figure:
    """Plots the optimization history from a twinning run.

    The figure follows the same rationale as :func:`plot_replay` and
    :func:`plot_twinning`: a vertical stack of subplots sharing a common x axis,
    here the function evaluation rather than time.

    1. **Parameter trajectory**: each unknown parameter's value per evaluation.
    2. **Log-prior**, **log-likelihood**, **log-posterior**: the three terms the
       twinner is optimising, one per subplot.

    Parameters
    ----------
    history : dict
        Dict with keys ``theta``, ``log_prior``, ``log_likelihood`` and
        ``log_posterior``. Each value is a list of recorded values per function
        evaluation. ``theta`` entries are numpy arrays.
    param_names : list of str or None, optional, default : None
        Optional list of parameter names for the theta legend. Defaults to
        ``theta_0``, ``theta_1``, ...
    figsize : tuple of float or None, optional, default : None
        Optional figure size. Defaults to a height that grows with the number of
        subplots.
    hover : bool, optional, default : True
        Whether to attach the interactive readout: hovering any subplot drops a
        vertical cursor across the whole stack and reports the values at that
        evaluation. Ignored on non-interactive backends.

    Returns
    -------
    matplotlib.figure.Figure
        The assembled figure.
    """
    thetas = np.array(history['theta'])          # (n_evals, n_params)
    log_prior = np.array(history['log_prior'])
    log_likelihood = np.array(history['log_likelihood'])
    log_posterior = np.array(history['log_posterior'])

    n_evals = len(log_posterior)
    steps = np.arange(n_evals)

    if thetas.ndim == 1:
        thetas = thetas[:, np.newaxis]

    n_params = thetas.shape[1]
    if param_names is None:
        param_names = [f'theta_{i}' for i in range(n_params)]

    fig, axes = make_figure(4, figsize or (10, 10))
    hover_series = {}

    # --- theta ---
    ax = axes[0]
    hover_series[ax] = []
    for i in range(n_params):
        mask = np.isfinite(thetas[:, i])
        ax.plot(steps[mask], thetas[mask, i], label=param_names[i], linewidth=0.8)
        hover_series[ax].append(series(param_names[i], steps[mask], thetas[mask, i], KIND_DENSE))
    ax.set_ylabel('Parameter value')
    set_panel_title(ax, 'Parameter trajectory')
    if n_params <= 10:
        ax.legend(fontsize='small', ncol=min(n_params, 4), framealpha=0.9)

    # --- the three optimisation terms ---
    terms = [
        ('Log-prior', log_prior, 'tab:blue'),
        ('Log-likelihood', log_likelihood, 'tab:orange'),
        ('Log-posterior', log_posterior, 'tab:green'),
    ]
    for ax, (label, values, color) in zip(axes[1:], terms):
        mask = np.isfinite(values)
        ax.plot(steps[mask], values[mask], color=color, linewidth=0.8)
        hover_series[ax] = [series(label, steps[mask], values[mask], KIND_DENSE)]
        ax.set_ylabel(label)
        set_panel_title(ax, label)

    finalize(fig, axes, 'Function evaluation', hover_series, hover, x_fmt='eval = {:.0f}')
    return fig