.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_surrogate_modeling/gaussian_process_regression/plot_gpr_active_learning.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_surrogate_modeling_gaussian_process_regression_plot_gpr_active_learning.py: Gaussian Process-based active learning for reliability ====================================================== .. GENERATED FROM PYTHON SOURCE LINES 5-7 .. code-block:: Python # sphinx_gallery_thumbnail_number = 13 .. GENERATED FROM PYTHON SOURCE LINES 8-12 In this example, we show how to sequentially add new points to a Gaussian Progress Regression model (GPR). The goal is to improve the predictivity of the surrogate model for reliability estimation. This kind of strategy is called "active learning". In order to create simple graphs, we consider a 1-d function. .. GENERATED FROM PYTHON SOURCE LINES 14-16 Create the function and the design of experiments ------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 18-23 .. code-block:: Python import openturns as ot import openturns.experimental as otexp import openturns.viewer as otv import numpy as np .. GENERATED FROM PYTHON SOURCE LINES 24-25 Define the function, the threshold above which the system is considered in failure, and the input probability distribution. .. GENERATED FROM PYTHON SOURCE LINES 25-29 .. code-block:: Python g = ot.SymbolicFunction(["x"], ["0.5*x^2 + sin(5*x)"]) threshold = 1.25 distribution = ot.Normal(0, 0.4) .. GENERATED FROM PYTHON SOURCE LINES 30-31 Create the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 31-39 .. code-block:: Python dimension = 1 DoESize = 4 xMin = -2.0 xMax = 2.0 X_distr = ot.Uniform(xMin, xMax) X = ot.LHSExperiment(X_distr, DoESize, False, False).generate() Y = g(X) .. GENERATED FROM PYTHON SOURCE LINES 40-41 We plot the limit state function, the initial design of experiments and the failure threshold. .. GENERATED FROM PYTHON SOURCE LINES 41-66 .. code-block:: Python thresholdFunction = ot.Curve([xMin, xMax], [threshold] * 2) thresholdFunction.setLineStyle("dashed") thresholdFunction.setColor("red") thresholdFunction.setLineWidth(2) thresholdFunction.setLegend("Failure threshold") data = ot.Cloud(X, Y) data.setColor("red") data.setPointStyle("circle") data.setLegend("Design of Experiments") graphFunction = g.draw(xMin, xMax) graphFunction.setColors(["magenta"]) graphFunction.setLegends(["Limit state function"]) graph = ot.Graph() graph.add(graphFunction) graph.add(thresholdFunction) graph.add(data) graph.setLegendPosition("lower right") graph.setAxes(True) graph.setGrid(True) view = otv.View(graph) .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_001.svg :alt: plot gpr active learning :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_001.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 67-69 Define the reliability analysis ------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 71-72 We define the event and estimate the reference failure probability with Monte-Carlo algorithm. .. GENERATED FROM PYTHON SOURCE LINES 74-90 .. code-block:: Python vect = ot.RandomVector(distribution) g = ot.MemoizeFunction(g) G = ot.CompositeRandomVector(g, vect) event = ot.ThresholdEvent(G, ot.Greater(), threshold) experiment = ot.MonteCarloExperiment() algo = ot.ProbabilitySimulationAlgorithm(event, experiment) algo.setMaximumCoefficientOfVariation(0.1) algo.setMaximumOuterSampling(int(1e5)) algo.run() result = algo.getResult() probability = result.getProbabilityEstimate() sampleX = g.getInputHistory() print("Reference probability on the real function =", probability) .. rst-class:: sphx-glr-script-out .. code-block:: none Reference probability on the real function = 0.016759776536312866 .. GENERATED FROM PYTHON SOURCE LINES 91-93 Create the algorithms --------------------- .. GENERATED FROM PYTHON SOURCE LINES 96-111 .. code-block:: Python def createMyBasicGPR(X, Y): """ Create a Gaussian Process 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 112-121 .. 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 122-205 .. code-block:: Python def plotMyBasicGPR( gprResult, xMin, xMax, X, Y, event, sampleX, refProbability, level=0.95 ): """ Given a gaussian process result, plot the data, the GP metamodel and a confidence interval. """ meta = gprResult.getMetaModel() graphKriging = meta.draw(xMin, xMax) graphKriging.setLegends(["GPR"]) # Create a grid of points and evaluate the function and the kriging nbpoints = 50 xGrid = linearSample(xMin, xMax, nbpoints) yFunction = g(xGrid) yKrig = meta(xGrid) # Compute the conditional covariance if event.getOperator().compare(1, 2): proba = ( np.sum(np.array(gprResult.getMetaModel()(sampleX)) < event.getThreshold()) / sampleX.getSize() ) else: proba = ( np.sum(np.array(gprResult.getMetaModel()(sampleX)) > event.getThreshold()) / sampleX.getSize() ) gpcc = otexp.GaussianProcessConditionalCovariance(gprResult) conditionalVariance = gpcc.getConditionalMarginalVariance(xGrid) conditionalSigma = np.sqrt(conditionalVariance) # Compute the quantile of the Normal distribution alpha = 1 - (1 - level) / 2 quantileAlpha = ot.DistFunc.qNormal(alpha) # Draw the bounds 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") # Plot the function graphFonction = ot.Curve(xGrid, yFunction) graphFonction.setLineStyle("dashed") graphFonction.setColor("magenta") graphFonction.setLineWidth(2) graphFonction.setLegend("Function") thresholdFunction = ot.Curve([xMin, xMax], [event.getThreshold()] * 2) thresholdFunction.setLineStyle("dashed") thresholdFunction.setColor("red") thresholdFunction.setLineWidth(2) thresholdFunction.setLegend("Threshold") # Draw the X and Y observed cloudDOE = ot.Cloud(X, Y) cloudDOE.setPointStyle("circle") cloudDOE.setColor("red") cloudDOE.setLegend("Data") # Assemble the graph graph = ot.Graph() graph.add(boundsPoly) graph.add(graphFonction) graph.add(thresholdFunction) graph.add(cloudDOE) graph.add(graphKriging) graph.setLegendPosition("lower right") graph.setAxes(True) graph.setGrid(True) graph.setTitle( "Estimated probability = %f, Reference probability = %f" % (proba, refProbability) ) graph.setXTitle("X") graph.setYTitle("Y") return graph .. GENERATED FROM PYTHON SOURCE LINES 206-208 We start by creating the initial Gaussian Process Regressor model :math:`\hat{\mathcal{M}}` on the 4 points in the design of experiments. We estimate the probability on this surrogate model and compare with the reference probability computed on the real limit state function. .. GENERATED FROM PYTHON SOURCE LINES 210-215 .. code-block:: Python gprResult = createMyBasicGPR(X, Y) graph = plotMyBasicGPR(gprResult, xMin, xMax, X, Y, event, sampleX, probability) view = otv.View(graph) .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_002.svg :alt: Estimated probability = 0.017606, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_002.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 216-218 Active learning Gaussian Process Regressor to sequentially add new points ------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 220-224 To sequentially add the new points, the "U criterion" is used. It consists in finding the new point as the sample :math:`\mathbf{x}` in the Monte-Carlo experiment that minimizes the following expression: :math:`\frac{ \left| T - \hat{\mathcal{M}} ( \mathbf{x}) \right|}{\hat{\sigma}(\mathbf{x})}` with :math:`\hat{\sigma}(\mathbf{x})` the square root of the marginal covariance of the Gaussian Process evaluated on :math:`\mathbf{x}`, and :math:`T` the event threshold (here 1.5) .. GENERATED FROM PYTHON SOURCE LINES 227-245 .. code-block:: Python def getNewPoint(X, gprResult, threshold): """ Returns a new point to be added to the design of experiments. This point maximizes the U criterion. """ gpcc = otexp.GaussianProcessConditionalCovariance(gprResult) response = gprResult.getMetaModel()(X) conditionalVariance = gpcc.getConditionalMarginalVariance(X) criterion = np.abs( ot.Sample([ot.Point([event.getThreshold()])] * X.getSize()) - response ) / np.sqrt(conditionalVariance + 1e-12) iMaxU = int(np.argmin(criterion)) xNew = X[iMaxU] return xNew .. GENERATED FROM PYTHON SOURCE LINES 246-247 We first call `getNewPoint` to get a point to add to the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 250-253 .. code-block:: Python xNew = getNewPoint(sampleX, gprResult, event.getThreshold()) print(xNew) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.821871] .. GENERATED FROM PYTHON SOURCE LINES 254-255 Then we evaluate the function on the new point and add it to the training design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 257-261 .. code-block:: Python yNew = g(xNew) X.add(xNew) Y.add(yNew) .. GENERATED FROM PYTHON SOURCE LINES 262-263 We now plot the updated Gaussian Process Regressor. .. GENERATED FROM PYTHON SOURCE LINES 265-269 .. code-block:: Python gprResult = createMyBasicGPR(X, Y) graph = plotMyBasicGPR(gprResult, xMin, xMax, X, Y, event, sampleX, probability) view = otv.View(graph) .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_003.svg :alt: Estimated probability = 0.000000, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_003.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 270-271 The algorithm added in the domain. .. GENERATED FROM PYTHON SOURCE LINES 273-282 .. code-block:: Python for GPRStep in range(10): xNew = getNewPoint(sampleX, gprResult, event) yNew = g(xNew) X.add(xNew) Y.add(yNew) gprResult = createMyBasicGPR(X, Y) graph = plotMyBasicGPR(gprResult, xMin, xMax, X, Y, event, sampleX, probability) otv.View(graph) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_004.svg :alt: Estimated probability = 0.000677, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_004.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_005.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_005.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_006.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_006.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_007.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_007.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_008.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_008.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_009.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_009.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_010.svg :alt: Estimated probability = 0.005756, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_010.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_011.svg :alt: Estimated probability = 0.013035, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_011.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_012.svg :alt: Estimated probability = 0.016421, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_012.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_013.svg :alt: Estimated probability = 0.016760, Reference probability = 0.016760 :srcset: /auto_surrogate_modeling/gaussian_process_regression/images/sphx_glr_plot_gpr_active_learning_013.svg :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 283-287 We can see that the metamodel only needs to be accurate near the event threshold to retrieve a precise estimation probability of failure. With only 10 points evaluated on the real limit state function, the metamodel accuracy is sufficient to estimate the failure probability. Indeed, the estimated probability is very close to the reference probability. This kind of active learning strategies allows one to save a large number of simulations. .. GENERATED FROM PYTHON SOURCE LINES 289-299 Conclusion ---------- The current example presents the naive implementation on the creation of a sequential design of experiments (active learning) based on GPR for failure probability estimation. See `Modules `_ for module `ot-ak` that implements active learning algorithms for reliability. More practical algorithms are presented in the following references. * Echard, B., Gayton, N., & Lemaire, M. (2011). AK-MCS: an active learning reliability method combining Kriging and Monte Carlo simulation. Structural Safety, 33(2), 145-154. * Echard, B. (2012). Assessment by kriging of the reliability of structures subjected to fatigue stress, Université Blaise Pascal, PhD thesis .. GENERATED FROM PYTHON SOURCE LINES 299-301 .. code-block:: Python otv.View.ShowAll() .. _sphx_glr_download_auto_surrogate_modeling_gaussian_process_regression_plot_gpr_active_learning.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_gpr_active_learning.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_gpr_active_learning.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_gpr_active_learning.zip `