.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_meta_modeling/kriging_metamodel/plot_kriging_sequential.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_meta_modeling_kriging_metamodel_plot_kriging_sequential.py: Sequentially adding new points to a Gaussian Process metamodel ============================================================== .. GENERATED FROM PYTHON SOURCE LINES 7-9 In this example, we show how to sequentially add new points to a Gaussian Process fitter (Kriging) in order to improve the predictivity of the metamodel. In order to create simple graphics, we consider a 1-d function. .. GENERATED FROM PYTHON SOURCE LINES 11-13 Create the function and the design of experiments ------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 15-22 .. code-block:: Python import openturns as ot import openturns.experimental as otexp from openturns.viewer import View import numpy as np from openturns import viewer .. GENERATED FROM PYTHON SOURCE LINES 23-26 .. code-block:: Python sampleSize = 4 dimension = 1 .. GENERATED FROM PYTHON SOURCE LINES 27-28 Define the function. .. GENERATED FROM PYTHON SOURCE LINES 30-32 .. code-block:: Python g = ot.SymbolicFunction(["x"], ["0.5*x^2 + sin(2.5*x)"]) .. GENERATED FROM PYTHON SOURCE LINES 33-34 Create the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 36-42 .. code-block:: Python xMin = -0.9 xMax = 1.9 X_distr = ot.Uniform(xMin, xMax) X = ot.LHSExperiment(X_distr, sampleSize, False, False).generate() Y = g(X) .. GENERATED FROM PYTHON SOURCE LINES 43-50 .. code-block:: Python graph = g.draw(xMin, xMax) data = ot.Cloud(X, Y) data.setColor("red") graph.add(data) view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_001.svg :alt: y0 as a function of x :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_001.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 51-53 Create the algorithms --------------------- .. GENERATED FROM PYTHON SOURCE LINES 56-71 .. code-block:: Python def createMyBasicGPfitter(X, Y): """ Create a Gaussian Process model from a pair of X and Y samples. We use a 3/2 Matérn covariance model and a constant trend. """ basis = ot.ConstantBasisFactory(dimension).build() covarianceModel = ot.MaternModel([1.0], 1.5) fitter = otexp.GaussianProcessFitter(X, Y, covarianceModel, basis) fitter.run() algo = otexp.GaussianProcessRegression(fitter.getResult()) algo.run() gprResult = algo.getResult() return gprResult .. GENERATED FROM PYTHON SOURCE LINES 72-81 .. code-block:: Python def linearSample(xmin, xmax, npoints): """Returns a sample created from a regular grid from xmin to xmax with npoints points.""" step = (xmax - xmin) / (npoints - 1) rg = ot.RegularGrid(xmin, step, npoints) vertices = rg.getVertices() return vertices .. GENERATED FROM PYTHON SOURCE LINES 82-83 The following `sqrt` function will be used later to compute the standard deviation from the variance. .. GENERATED FROM PYTHON SOURCE LINES 85-88 .. code-block:: Python sqrt = ot.SymbolicFunction(["x"], ["sqrt(x)"]) .. GENERATED FROM PYTHON SOURCE LINES 89-152 .. code-block:: Python def plotMyBasicGPfitter(gprResult, xMin, xMax, X, Y, level=0.95): """ Given a metamodel result, plot the data, the GP fitter metamodel and a confidence interval. """ samplesize = X.getSize() meta = gprResult.getMetaModel() graphKriging = meta.draw(xMin, xMax) graphKriging.setLegends(["Gaussian Process Regression"]) # Create a grid of points and evaluate the function and the metamodel nbpoints = 50 xGrid = linearSample(xMin, xMax, nbpoints) yFunction = g(xGrid) yKrig = meta(xGrid) # Compute the conditional covariance gpcc = otexp.GaussianProcessConditionalCovariance(gprResult) epsilon = ot.Sample(nbpoints, [1.0e-8]) conditionalVariance = gpcc.getConditionalMarginalVariance(xGrid) + epsilon conditionalSigma = sqrt(conditionalVariance) # Compute the quantile of the Normal distribution alpha = 1 - (1 - level) / 2 quantileAlpha = ot.DistFunc.qNormal(alpha) # Graphics of the bounds epsilon = 1.0e-8 dataLower = [ yKrig[i, 0] - quantileAlpha * conditionalSigma[i, 0] for i in range(nbpoints) ] dataUpper = [ yKrig[i, 0] + quantileAlpha * conditionalSigma[i, 0] for i in range(nbpoints) ] # Compute the Polygon graphics boundsPoly = ot.Polygon.FillBetween(xGrid.asPoint(), dataLower, dataUpper) boundsPoly.setLegend("95% bounds") # Validate the metamodel metamodelPredictions = meta(xGrid) mmv = ot.MetaModelValidation(yFunction, metamodelPredictions) r2Score = mmv.computeR2Score()[0] # Plot the function graphFonction = ot.Curve(xGrid, yFunction) graphFonction.setLineStyle("dashed") graphFonction.setColor("magenta") graphFonction.setLineWidth(2) graphFonction.setLegend("Function") # Draw the X and Y observed cloudDOE = ot.Cloud(X, Y) cloudDOE.setPointStyle("circle") cloudDOE.setColor("red") cloudDOE.setLegend("Data") # Assemble the graphics graph = ot.Graph() graph.add(boundsPoly) graph.add(graphFonction) graph.add(cloudDOE) graph.add(graphKriging) graph.setLegendPosition("lower right") graph.setAxes(True) graph.setGrid(True) graph.setTitle("Size = %d, R2=%.2f%%" % (samplesize, 100 * r2Score)) graph.setXTitle("X") graph.setYTitle("Y") return graph .. GENERATED FROM PYTHON SOURCE LINES 153-154 We start by creating the initial GP fitter metamodel on the 4 points in the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 156-161 .. code-block:: Python gprResult = createMyBasicGPfitter(X, Y) graph = plotMyBasicGPfitter(gprResult, xMin, xMax, X, Y) view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_002.svg :alt: Size = 4, R2=96.34% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_002.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 162-164 Sequentially add new points --------------------------- .. GENERATED FROM PYTHON SOURCE LINES 166-167 The following function is the building block of the algorithm. It returns a new point which maximizes the conditional variance. .. GENERATED FROM PYTHON SOURCE LINES 170-185 .. code-block:: Python def getNewPoint(xMin, xMax, gprResult): """ Returns a new point to be added to the design of experiments. This point maximizes the conditional variance of the metamodel. """ nbpoints = 50 xGrid = linearSample(xMin, xMax, nbpoints) gpcc = otexp.GaussianProcessConditionalCovariance(gprResult) conditionalVariance = gpcc.getConditionalMarginalVariance(xGrid) iMaxVar = int(np.argmax(conditionalVariance)) xNew = xGrid[iMaxVar, 0] xNew = ot.Point([xNew]) return xNew .. GENERATED FROM PYTHON SOURCE LINES 186-187 We first call `getNewPoint` to get a point to add to the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 190-193 .. code-block:: Python xNew = getNewPoint(xMin, xMax, gprResult) xNew .. raw:: html
class=Point name=Unnamed dimension=1 values=[-0.9]


.. GENERATED FROM PYTHON SOURCE LINES 194-195 Then we evaluate the function on the new point and add it to the training design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 197-201 .. code-block:: Python yNew = g(xNew) X.add(xNew) Y.add(yNew) .. GENERATED FROM PYTHON SOURCE LINES 202-203 We now plot the updated GP fitter. .. GENERATED FROM PYTHON SOURCE LINES 205-206 sphinx_gallery_thumbnail_number = 3 .. GENERATED FROM PYTHON SOURCE LINES 206-211 .. code-block:: Python gprResult = createMyBasicGPfitter(X, Y) graph = plotMyBasicGPfitter(gprResult, xMin, xMax, X, Y) graph.setTitle("GP fitter #0") view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_003.svg :alt: GP fitter #0 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_003.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 212-213 The algorithm added a point to the right bound of the domain. .. GENERATED FROM PYTHON SOURCE LINES 215-225 .. code-block:: Python for krigingStep in range(5): xNew = getNewPoint(xMin, xMax, gprResult) yNew = g(xNew) X.add(xNew) Y.add(yNew) gprResult = createMyBasicGPfitter(X, Y) graph = plotMyBasicGPfitter(gprResult, xMin, xMax, X, Y) graph.setTitle("GP fitter #%d " % (krigingStep + 1) + graph.getTitle()) View(graph) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_004.svg :alt: GP fitter #1 Size = 6, R2=98.80% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_004.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_005.svg :alt: GP fitter #2 Size = 7, R2=99.74% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_005.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_006.svg :alt: GP fitter #3 Size = 8, R2=99.87% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_006.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_007.svg :alt: GP fitter #4 Size = 9, R2=99.89% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_007.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_008.svg :alt: GP fitter #5 Size = 10, R2=99.90% :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_sequential_008.svg :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 226-230 We observe that the second added point is the left bound of the domain. The remaining points were added strictly inside the domain where the accuracy was drastically improved. With only 10 points, the metamodel accuracy is already very good with a :math:`Q^2` which is equal to 99.9%. .. GENERATED FROM PYTHON SOURCE LINES 232-237 Conclusion ---------- The current example presents the naive implementation on the creation of a sequential design of experiments based on Gaussian Process metamodel. More practical algorithms are presented in [ginsbourger2018]_. .. GENERATED FROM PYTHON SOURCE LINES 237-239 .. code-block:: Python View.ShowAll() .. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_plot_kriging_sequential.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_kriging_sequential.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_kriging_sequential.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_kriging_sequential.zip `