History Matching

import torch

import matplotlib.pyplot as plt

from autoemulate.emulators import GaussianProcess
from autoemulate.core.compare import AutoEmulate
from autoemulate.simulations.projectile import Projectile
from autoemulate.calibration.history_matching import HistoryMatching

In this tutorial we demonstrate the use of History Matching to determine which points in the input space are plausible given a set of observations.

Performing History Matching requires:

  • a fit emulator that predicts uncertainty (e.g., Gaussian Process) and

  • an observation associated with the simulator output.

The emulator is used to efficiently estimate the simulator output, accounting for all uncertainties. The emulated output is then compared with the observation and parameter inputs that are unlikely to produce the observation can then be “ruled out” as implausible, reducing the input space.

History Matching#

1. Simulate data and train an emulator#

In this example we use the Projectile simulator. It returns distance travelled by a projectile given a drag coefficient c and initial velocity v0.

simulator = Projectile(log_level="error", parameters_range={'c': (-5.0, 1.0), 'v0': (0.0, 10)})
x = simulator.sample_inputs(1000)
y, _ = simulator.forward_batch(x)

plt.scatter(x[:, 0], x[:, 1], c=y, cmap='viridis')
plt.xlabel('Drag coefficient')
plt.ylabel('Initial velocity')
plt.colorbar(label="Distance")
plt.show

<function matplotlib.pyplot.show(close=None, block=None)>
../../_images/85de18fec89a41a830185994cd74dc9508f948b0781ab5d1240a9069902f1056.png

Since we need an uncertainty aware emulator, we restrict AutoEmulate to only train a Gaussian Process. For the purposes of this tutorial, we use default parameters.

ae = AutoEmulate(x, y, models=[GaussianProcess], model_params={}, log_level="error")

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.010147 0.999988 0.000002 0.999989 7.142013e-07 
model = ae.best_result().model

2. Calibrate#

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

true_drag = -2.0
true_velocity = 7.0
distance = simulator.forward(torch.tensor([[true_drag, true_velocity]]))
stdev = distance.item() * 0.05
observation = {"distance": (distance.item(), stdev**2)}

distance

tensor([[6.2728]], dtype=torch.float64)

#  observed data: (mean, variance)
hm = HistoryMatching(observations=observation)

We also need predictions for a set of query points using the trained emulator. These must have uncertainty estimates.

x_new = simulator.sample_inputs(1000)
mean, variance = model.predict_mean_and_variance(x_new)

Primary use of the HistoryMatching class is the calculate_implausibility method, which returns the implausibility metric (a number of standard deviations between the emulator mean and the observation) for the queried points.

implausability = hm.calculate_implausibility(mean, variance)

# first 10 implausability scores
implausability[:10]

tensor([[16.0332],
        [15.9673],
        [ 0.6845],
        [ 7.4618],
        [ 5.2313],
        [19.1409],
        [16.5915],
        [ 0.7034],
        [19.5580],
        [16.8082]])

The get_nroy method can be used to filter the queried points given the implausability scores and only retain those that have not been ruled out yet (NROY).

nroy_params = hm.get_nroy(implausability, x_new)

Below we plot the NROY parameters against the true parameters used to generate the observed data. History matching has successfully ruled out a portion of the input space.

fig, ax = plt.subplots(1, 2, figsize=(10, 4))

ax[0].boxplot(nroy_params[:, 0], vert=True, patch_artist=True, boxprops=dict(facecolor='lightblue'))
ax[0].axhline(true_drag, color="k", linestyle="--", label="True drag")
ax[0].set_ylabel("Drag coefficient")
ax[0].legend()

ax[1].boxplot(nroy_params[:, 1], vert=True, patch_artist=True, boxprops=dict(facecolor='lightgreen'))
ax[1].axhline(true_velocity, color="k", linestyle="--", label="True initial velocity")
ax[1].legend()
ax[1].set_ylabel("Initial velocity")

plt.suptitle("Box plot of NROY parameters after history matching")
plt.show()
../../_images/b8f8d550ea6e82bd5c1df8b4cd6c00c9bea0ea4a1d46b2694ebccfd1bba9546d.png
Contents