Example 40: set inversion on a fixed grid

Script: examples/example40_setinversion_gridset.py

Purpose

The script estimates the inverse image of an output-space box on a fixed input grid. The target set is

\[\Gamma_B = \{x : f(x) \in B\}.\]

At each step, posterior box-membership probabilities and a box weighted-MSE criterion are evaluated on the grid. The box-membership probability is the same type of Gaussian probability used for constrained Bayesian optimization [3].

What is computed

• posterior mean and variance on the fixed grid.

• box-membership probabilities P(Y(x) in B | observations).

• box_wMSE values used for selecting the next grid point.

• one new evaluation per sequential step.

Main objects

• gpmpcontrib.optim.setinversion.SetInversionGridSearch

• gpmpcontrib.samplingcriteria.box_probability

• gpmpcontrib.samplingcriteria.box_wMSE

Outputs

Run python examples/example40_setinversion_gridset.py from the repository root to execute the example. Regenerate the static figure with cd docs && python make_example_results.py.

Top panel: target output box, observations, and GP posterior. Middle panel: criterion used to select new evaluations. Lower panel: posterior probability that the function belongs to the target box at each input location.

Source excerpt

# -- define a mono-objective problem


problem = gpc.ComputerExperiment(
    1,  # dim of search space
    [[-2.0], [2.0]],  # search box
    single_function=test_function,  # test function
)


# -- create initial dataset

nt = 2000
xt = gp.misc.designs.regulargrid(problem.input_dim, nt, problem.input_box)
zt = problem(xt)

ni = 4
ind = [i * 20 for i in [20, 46, 60, 80]]
xi = xt[ind]

box_target = gnp.array([[0.99], [1.02]])

# -- initialize a model and the ei algorithm
model = gpc.Model_ConstantMean_Maternp_REMAP(
    "GP1d",
    output_dim=problem.output_dim,
    mean_specification={"type": "constant"},
    covariance_specification={"p": 2},
)

beta = 1.0
algo = si.SetInversionGridSearch(problem, model, xt, box_target, options={"beta": beta})

algo.set_initial_design(xi)

# -- visualization


def plot(show=True, x=None, z=None):
    zpm, zpv = algo.predict(xt, convert_out=False)
    crit0, _ = sampcrit.box_wMSE(algo.box_target, zpm, zpv, beta=beta)
    # crit1, _ = sampcrit.box_wMSE(algo.box_target, zpm, zpv, beta=0.5)
    pe, _ = sampcrit.box_probability(algo.box_target, zpm, zpv)
    
    fig = gp.plot.Figure(nrows=3, ncols=1, isinteractive=True)
    fig.subplot(1)
    fig.plot(xt, zt, "k", linewidth=0.5)
    fig.axhline(box_target[0], color="k", linewidth=0.5)
    fig.axhline(box_target[1], color="k", linewidth=0.5)
    if z is not None: