from autoemulate.emulators import GaussianProcess
from autoemulate.core.compare import AutoEmulate
from autoemulate.simulations.epidemic import Epidemic
from autoemulate.calibration.history_matching import HistoryMatchingWorkflow
from autoemulate.calibration.history_matching_dashboard import HistoryMatchingDashboard
random_seed = 42

History Matching workflow#

The HistoryMatching calibration method can be a useful way to iteratively decide which simulations to run to generate data refine the emulator on. The HistoryMatchingWorkflow implements this iterative sample-predict-refit workflow. Each time it is run:

  • parameters are sampled from the not ruled out yet (NROY) space

  • an emulator is used in combination with HistoryMatching to score the implausability of the samples

  • simulations are run for a subset of the NROY samples

  • the emulator is refit given the newly simulated data

In this tutorial, we demonstrate how to implement this simulator in the loop workflow.

1. Simulate data and train an emulator#

In this example, we’ll use the Epidemic simulator, which returns the peak infection rate given two input parameters, beta(the transimission rate per day) and gamma (the recovery rate per day).

simulator = Epidemic(log_level="error")
x = simulator.sample_inputs(500, random_seed=random_seed)
y = simulator.forward_batch(x)

For the purposes of this tutorial, we will restrict the model choice to GaussianProcess.

ae = AutoEmulate(
    x,
    y,
    models=[GaussianProcess],
    model_tuning=False,
    log_level="error",
    random_seed=random_seed
)

We can verify that the fitted emulator performs well on both the train and test data.

ae.summarise()
model_name x_transforms y_transforms params rmse_test r2_test r2_test_std r2_train r2_train_std
0 GaussianProcess [StandardizeTransform()] [StandardizeTransform()] {} 0.038385 0.999956 0.000011 0.999965 0.000004
model = ae.best_result().model

2. Calibrate#

To instantiate the HistoryMatchingWorkflow object, we need an observed mean and, optionally, variance for each simulator output.

observations = {"infection_rate": (0.3, 0.05)}

We also pass the fitted emulator and the simulator to the HistoryMatchingWorkflow.

hmw = HistoryMatchingWorkflow(
    simulator=simulator,
    emulator=model,
    observations=observations,
    threshold=3.0,
    train_x=x,
    train_y=y, 
    random_seed=random_seed
)

The run method implements the iterative sample-predict-refit workflow:

  • sample n_test_samples to test from the not ruled out yet (NROY) space

  • use emulator to filter out implausible samples and update the NROY space

  • run n_simulations predictions for the sampled parameters using the simulator

  • refit the emulator using the simulated data

The HistoryMatchingWorkflow object maintains and updates the internal state each time run() is called so the full workflow can be run over a number of iterations.

test_parameters, impl_scores = hmw.run(n_simulations=20, n_test_samples=100)

3. Visualise results#

For visualising results, you can use the HistoryMatchingDashboard which implements a number of interactive plots. To initialize it just pass in outputs of the run() method and information about the simulation.

dashboard = HistoryMatchingDashboard(
    samples=test_parameters,
    impl_scores=impl_scores,
    param_names=simulator.param_names,  
    output_names=simulator.output_names, 
)

The plots can be viewed using the display() method. Below you can see an example image of what the dashboard looks like.

# dashboard.display()
Work Flow