Visualization & Diagnostics¶

History matching produces a lot to look at — emulator fits, the shrinking plausible region, how close the simulated outputs sit to their targets. This notebook tours the built-in plotting and display API on the stochastic SIR model used throughout these tutorials.

Every plot_* method returns Matplotlib axes and accepts an ax= argument, so the figures render inline and compose into your own layouts. The same figures are written to the engine's output_dir after each wave.

See the Diagnostics & plotting guide for how to interpret each figure.

In [1]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

import historymatching as hm
from model import SIR, generate_observed_data

%matplotlib inline
np.random.seed(0)

import numpy as np import pandas as pd import matplotlib.pyplot as plt import historymatching as hm from model import SIR, generate_observed_data %matplotlib inline np.random.seed(0)

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1780534098.224547 1042889 port.cc:153] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.
I0000 00:00:1780534098.287259 1042889 cpu_feature_guard.cc:227] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
To enable the following instructions: AVX2 AVX_VNNI FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.

WARNING: All log messages before absl::InitializeLog() is called are written to STDERR
I0000 00:00:1780534099.532026 1042889 port.cc:153] oneDNN custom operations are on. You may see slightly different numerical results due to floating-point round-off errors from different computation orders. To turn them off, set the environment variable `TF_ENABLE_ONEDNN_OPTS=0`.

/software/conda/envs/hm/lib/python3.13/site-packages/gpflow/versions.py:1: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
  import pkg_resources

Set up and run a calibration¶

We recover the transmission/recovery rates of a synthetic outbreak, then use the result to demonstrate each plot. (See 01_basic_workflow.ipynb for the workflow itself in detail.)

In [2]:

beta_true, gamma_true = 1.3, 0.5
population_size, n_seed = 10_000, 100
incidence_obs, _ = generate_observed_data(beta_true=beta_true, gamma_true=gamma_true,
                                          population_size=population_size,
                                          n_seed_infections=n_seed)

def sir_sim(samples: pd.DataFrame) -> pd.DataFrame:
    df = samples.copy()
    if "rand_seed" not in df.columns:
        df["rand_seed"] = np.random.default_rng(0).integers(0, 2**31, size=len(df))
    rows = []
    for _, row in df.iterrows():
        model = SIR(beta=row["beta"], gamma=row["gamma"],
                    s0=population_size - n_seed, i0=n_seed, seed=int(row["rand_seed"]))
        inc = model.get_incidence()
        rec = {f"incidence_{i}": inc[i] for i in range(len(inc))}
        rec["peak_incidence"] = inc.max()
        rec["total_cases"] = inc.sum()
        rows.append(rec)
    return pd.DataFrame(rows)

observations = {
    "peak_incidence": (incidence_obs.max(), 50),
    "total_cases":    (incidence_obs.sum(), 200),
    "incidence_10":   (incidence_obs.iloc[10], 40),
}

engine = hm.HistoryMatching(
    bounds={"beta": (0.5, 3.0), "gamma": (0.1, 1.0)},
    observations=observations,
    function=sir_sim,
    emulator_type="gpr",
    n_samples=300,
    max_iterations=3,
    random_seed=123,
    feature_selection=["peak_incidence", "incidence_10"],
    output_dir=None,  # focus on interactive plots; no files written here
)
results = engine.run()
print(f"Completed {len(results)} waves.")

beta_true, gamma_true = 1.3, 0.5 population_size, n_seed = 10_000, 100 incidence_obs, _ = generate_observed_data(beta_true=beta_true, gamma_true=gamma_true, population_size=population_size, n_seed_infections=n_seed) def sir_sim(samples: pd.DataFrame) -> pd.DataFrame: df = samples.copy() if "rand_seed" not in df.columns: df["rand_seed"] = np.random.default_rng(0).integers(0, 2**31, size=len(df)) rows = [] for _, row in df.iterrows(): model = SIR(beta=row["beta"], gamma=row["gamma"], s0=population_size - n_seed, i0=n_seed, seed=int(row["rand_seed"])) inc = model.get_incidence() rec = {f"incidence_{i}": inc[i] for i in range(len(inc))} rec["peak_incidence"] = inc.max() rec["total_cases"] = inc.sum() rows.append(rec) return pd.DataFrame(rows) observations = { "peak_incidence": (incidence_obs.max(), 50), "total_cases": (incidence_obs.sum(), 200), "incidence_10": (incidence_obs.iloc[10], 40), } engine = hm.HistoryMatching( bounds={"beta": (0.5, 3.0), "gamma": (0.1, 1.0)}, observations=observations, function=sir_sim, emulator_type="gpr", n_samples=300, max_iterations=3, random_seed=123, feature_selection=["peak_incidence", "incidence_10"], output_dir=None, # focus on interactive plots; no files written here ) results = engine.run() print(f"Completed {len(results)} waves.")

I0000 00:00:1780534103.638570 1042889 gpu_device.cc:2043] Created device /job:localhost/replica:0/task:0/device:GPU:0 with 2175 MB memory:  -> device: 0, name: NVIDIA GeForce RTX 3050 Laptop GPU, pci bus id: 0000:01:00.0, compute capability: 8.6

I0000 00:00:1780534108.354130 1043120 cuda_solvers.cc:175] Creating GpuSolver handles for stream 0x59e111230a40

Completed 3 waves.

Text summaries¶

engine.get_status_summary() folds the per-wave convergence and the surviving parameter ranges into one report. engine.nroy_summary() returns the same per-parameter information as a tidy DataFrame.

In [3]:

print(engine.get_status_summary())

print(engine.get_status_summary())

=== History Matching Status ===
State: completed
Progress: 3/3 waves
Acceptance rate: 8.7%
Total samples generated: 17,236
Total samples accepted: 1,800
Emulators trained: 6
[OK] Simulator function configured
[OK] All waves completed successfully

In [4]:

engine.nroy_summary()

engine.nroy_summary()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[4], line 1
----> 1 engine.nroy_summary()

AttributeError: 'HistoryMatching' object has no attribute 'nroy_summary'

Convergence¶

The NROY fraction is the share of fresh prior samples passing all emulator constraints so far. A steadily falling fraction means the calibration is successfully ruling out parameter space.

In [5]:

engine.plot_convergence();

engine.plot_convergence();

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[5], line 1
----> 1 engine.plot_convergence();

AttributeError: 'HistoryMatching' object has no attribute 'plot_convergence'

The NROY parameter cloud¶

The headline result: a corner plot of the non-implausible region — marginals on the diagonal, pairwise scatter below. We overlay the known true values with truth=.

In [6]:

truth = {"beta": beta_true, "gamma": gamma_true}
engine.plot_nroy(truth=truth);

truth = {"beta": beta_true, "gamma": gamma_true} engine.plot_nroy(truth=truth);

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[6], line 2
      1 truth = {"beta": beta_true, "gamma": gamma_true}
----> 2 engine.plot_nroy(truth=truth);

AttributeError: 'HistoryMatching' object has no attribute 'plot_nroy'

plot_marginals shows just the one-dimensional posteriors, with the sample median and the true value marked.

In [7]:

engine.plot_marginals(truth=truth);

engine.plot_marginals(truth=truth);

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[7], line 1
----> 1 engine.plot_marginals(truth=truth);

AttributeError: 'HistoryMatching' object has no attribute 'plot_marginals'

Z-scores vs targets¶

For every target, the distribution of (simulated - target) / target_std across the NROY samples, by wave. Bands inside the green ±threshold region and centred on zero are consistent with the data.

In [8]:

engine.plot_zscores();

engine.plot_zscores();

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[8], line 1
----> 1 engine.plot_zscores();

AttributeError: 'HistoryMatching' object has no attribute 'plot_zscores'

Constrained directions¶

A PCA view of which combinations of parameters the data constrained most — useful for spotting non-identifiable directions that marginal plots miss.

In [9]:

engine.plot_constrained_dims(n_top=2);

engine.plot_constrained_dims(n_top=2);

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[9], line 1
----> 1 engine.plot_constrained_dims(n_top=2);

AttributeError: 'HistoryMatching' object has no attribute 'plot_constrained_dims'

Emulator quality¶

The emulators are only trustworthy if they reproduce the simulator. Inspect the fit per feature for the final wave.

In [10]:

final = results[-1]
final.quality_table()

final = results[-1] final.quality_table()

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[10], line 2
      1 final = results[-1]
----> 2 final.quality_table()

AttributeError: 'IterationResult' object has no attribute 'quality_table'

In [11]:

final.plot_emulator_quality();

final.plot_emulator_quality();

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[11], line 1
----> 1 final.plot_emulator_quality();

AttributeError: 'IterationResult' object has no attribute 'plot_emulator_quality'

In [12]:

final.plot_predicted_vs_actual("peak_incidence");

final.plot_predicted_vs_actual("peak_incidence");

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[12], line 1
----> 1 final.plot_predicted_vs_actual("peak_incidence");

AttributeError: 'IterationResult' object has no attribute 'plot_predicted_vs_actual'

Ensemble fan plot vs observed data¶

History matching constrains parameters; to check the outputs, re-run the simulator at NROY parameter sets and compare the trajectories to the observed data. hm.HistoryMatching.plot_ensemble_fan is model-agnostic — give it a 2-D array of trajectories and an observed series.

In [13]:

nroy = engine.get_nroy_samples(40, method="lhs")  # unbiased draw
trajectories = np.array([
    SIR(beta=r["beta"], gamma=r["gamma"], s0=population_size - n_seed, i0=n_seed).get_incidence()
    for _, r in nroy.iterrows()
])
hm.HistoryMatching.plot_ensemble_fan(trajectories, observed=incidence_obs.values,
                     xlabel="Day", ylabel="Daily incidence",
                     title="NROY trajectories vs observed");

nroy = engine.get_nroy_samples(40, method="lhs") # unbiased draw trajectories = np.array([ SIR(beta=r["beta"], gamma=r["gamma"], s0=population_size - n_seed, i0=n_seed).get_incidence() for _, r in nroy.iterrows() ]) hm.HistoryMatching.plot_ensemble_fan(trajectories, observed=incidence_obs.values, xlabel="Day", ylabel="Daily incidence", title="NROY trajectories vs observed");

Domain-object views¶

ParameterSpace, ObservationData, and EmulatorBank each have a summary() and (where it makes sense) a plot.

In [14]:

print(engine.parameter_space.summary())
print()
print(engine.observations.summary())
print()
print(engine.emulator_bank.summary())

print(engine.parameter_space.summary()) print() print(engine.observations.summary()) print() print(engine.emulator_bank.summary())

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[14], line 1
----> 1 print(engine.parameter_space.summary())
      2 print()
      3 print(engine.observations.summary())
      4 print()

AttributeError: 'ParameterSpace' object has no attribute 'summary'

In [15]:

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
engine.parameter_space.plot_bounds(reference=engine.parameter_space, ax=axes[0])
engine.observations.plot_targets(ax=axes[1])
plt.tight_layout();

fig, axes = plt.subplots(1, 2, figsize=(12, 4)) engine.parameter_space.plot_bounds(reference=engine.parameter_space, ax=axes[0]) engine.observations.plot_targets(ax=axes[1]) plt.tight_layout();

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Cell In[15], line 2
      1 fig, axes = plt.subplots(1, 2, figsize=(12, 4))
----> 2 engine.parameter_space.plot_bounds(reference=engine.parameter_space, ax=axes[0])
      3 engine.observations.plot_targets(ax=axes[1])
      4 plt.tight_layout();

AttributeError: 'ParameterSpace' object has no attribute 'plot_bounds'

All of these methods return Matplotlib axes, so you can drop them into your own multi-panel figures (as above) or save them with fig.savefig(...). For the full catalogue and reading guidance, see the Diagnostics & plotting guide.