Dataloader-based parameter selection

This example shows how observations can be passed through Dataset and DataLoader objects during parameter selection. Use this form when the selection criterion should be evaluated on batches instead of a single full array at every call.

What this example does

The script creates a test problem, stores observations in a Dataset, wraps it with a DataLoader, and calls select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior with the dataloader argument. The selected model is then used for prediction and the preview plots predicted values against reference values.

Mathematical description

The model and REMAP criterion are the same as in the array-based examples. The dataloader changes how the criterion is evaluated. If batches \(b_1,\ldots,b_B\) have sizes \(n_1,\ldots,n_B\), the batch wrapper evaluates a weighted scalar objective of the form

\[\overline J(\theta) = \frac{\sum_{\ell=1}^B n_\ell J_\ell(\theta)} {\sum_{\ell=1}^B n_\ell},\]

where \(J_\ell(\theta)\) is the selection criterion evaluated on batch \(b_\ell\). With batches_per_eval=0, one criterion call uses the full loader. With a positive batches_per_eval, one criterion call uses only that many successive batches, cycling through the loader.

Outputs

The displayed quantities are GP predictions and reference values at test points. Points near the diagonal indicate accurate predictions. Systematic deviations from the diagonal suggest bias, poor covariance parameters, or insufficient observations.

Functions used

  • Selection functions accept either explicit xi, zi arrays or a dataloader. Do not pass both.

  • DataLoader controls batching. The selection criterion still returns a scalar objective for optimization.

  • The selected model.covparam is used normally by model.predict after batched selection.

../_images/dataloader_0_0.png

Script: examples/gpmp_example40_dataloader.py

  1"""
  2Prediction of some classical test functions in dimension > 2
  3
  4An anisotropic Matern covariance function is used for the Gaussian
  5Process (GP) prior. The parameters of this covariance function
  6(variance and ranges) are estimated using the Restricted Maximum
  7A Posteriori (ReMAP).
  8
  9This script uses a DataLoader object to structure the input dataset
 10into fixed-size batches for parameter estimation.
 11
 12----
 13Author: Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr>
 14Copyright (c) 2022-2026, CentraleSupelec
 15License: GPLv3 (see LICENSE)
 16"""
 17import gpmp.num as gnp
 18import gpmp as gp
 19from gpmp.dataloader import Dataset, DataLoader
 20
 21
 22def choose_test_case(problem, ni=2000):
 23    if problem == 1:
 24        problem_name = "Hartmann4"
 25        f = gp.misc.testfunctions.hartmann4
 26        dim = 4
 27        box = [[0.0] * 4, [1.0] * 4]
 28        xi = gp.misc.designs.ldrandunif(dim, ni, box)
 29        nt = 1000
 30        xt = gp.misc.designs.ldrandunif(dim, nt, box)
 31
 32    elif problem == 2:
 33        problem_name = "Hartmann6"
 34        f = gp.misc.testfunctions.hartmann6
 35        dim = 6
 36        box = [[0.0] * 6, [1.0] * 6]
 37        xi = gp.misc.designs.ldrandunif(dim, ni, box)
 38        nt = 1000
 39        xt = gp.misc.designs.ldrandunif(dim, nt, box)
 40
 41    elif problem == 3:
 42        problem_name = "Borehole"
 43        f = gp.misc.testfunctions.borehole
 44        dim = 8
 45        box = [
 46            [0.05, 100.0, 63070.0, 990.0, 63.1, 700.0, 1120.0, 9855.0],
 47            [0.15, 50000.0, 115600.0, 1110.0, 116.0, 820.0, 1680.0, 12045.0],
 48        ]
 49        xi = gp.misc.designs.maximinldlhs(dim, ni, box)
 50        nt = 1000
 51        xt = gp.misc.designs.ldrandunif(dim, nt, box)
 52
 53    elif problem == 4:
 54        problem_name = "detpep8d"
 55        f = gp.misc.testfunctions.detpep8d
 56        dim = 8
 57        box = [[0.0] * 8, [1.0] * 8]
 58        xi = gp.misc.designs.maximinldlhs(dim, ni, box)
 59        nt = 1000
 60        xt = gp.misc.designs.ldrandunif(dim, nt, box)
 61
 62    return problem_name, f, dim, box, ni, xi, nt, xt
 63
 64
 65def constant_mean(x, param):
 66    return gnp.ones((x.shape[0], 1))
 67
 68
 69def kernel(x, y, covparam, pairwise=False):
 70    p = 10
 71    return gp.kernel.maternp_covariance(x, y, p, covparam, pairwise)
 72
 73
 74def visualize_predictions(problem_name, zt, zpm):
 75    fig = gp.plot.Figure()
 76    fig.plot(zt, zpm, "ko", markersize=3)
 77    (xmin, xmax), (ymin, ymax) = fig.ax.get_xlim(), fig.ax.get_ylim()
 78    xmin = min(xmin, ymin)
 79    xmax = max(xmax, ymax)
 80    fig.plot([xmin, xmax], [xmin, xmax], "--", linewidth=1)
 81    fig.xylabels("reference values", "posterior mean")
 82    fig.title(f"{problem_name}: test predictions")
 83    fig.grid()
 84    fig.show()
 85
 86
 87def main():
 88    problem = 2
 89    problem_name, f, dim, box, ni, xi, nt, xt = choose_test_case(problem, ni=1000)
 90
 91    zi = f(xi)
 92    zt = f(xt)
 93    dataset = Dataset(xi, zi)
 94    loader = DataLoader(dataset, batch_size=200, shuffle=False)
 95
 96    model = gp.core.Model(constant_mean, kernel)
 97
 98    model, info = (
 99        gp.kernel.select_parameters_sigma2_rho_with_remap_logsigma2_logrho_prior(
100            model, dataloader=loader, info=True
101        )
102    )
103
104    gp.modeldiagnosis.diag(
105        model, "linear_mean_maternp_anisotropic", info, xi, zi
106    )
107
108    gp.modeldiagnosis.plot_selection_criterion_crosssections(
109        info=info,
110        selection_criterion=info.selection_criterion_nograd,
111        covparam=None,
112        n_points=100,
113        param_names=None,
114        criterion_name="selection criterion",
115        criterion_name_full="Cross sections for negative log restricted likelihood",
116        ind=[0, 1, 2],
117        ind_pooled=None,
118        param_box=None,
119        param_box_pooled=None,
120        delta=5.0,
121    )
122    
123    (zpm, zpv) = model.predict(xi, zi, xt)
124
125    visualize_predictions(problem_name, zt, zpm)
126
127    zloom, zloov, eloo = model.loo(xi, zi)
128    gp.plot.plot_loo(zi, zloom, zloov)
129
130    gp.modeldiagnosis.perf(
131        model,
132        xi,
133        zi,
134        loo=True,
135        loo_res=(zloom, zloov, eloo),
136        xtzt=(xt, zt),
137        zpmzpv=(zpm, zpv),
138    )
139
140
141if __name__ == "__main__":
142    main()