.. 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_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_meta_modeling_kriging_metamodel_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-24 .. 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 25-26 Define the function, the threshold above which the system is considered in failure, and the input probability distribution. .. GENERATED FROM PYTHON SOURCE LINES 26-30 .. 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 31-32 Create the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 32-40 .. 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 41-42 We plot the limit state function, the initial design of experiments and the failure threshold. .. GENERATED FROM PYTHON SOURCE LINES 42-67 .. 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 = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_001.svg :alt: plot gpr active learning :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_001.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 68-70 Define the reliability analysis ------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 72-73 We define the event and estimate the reference failure probability with Monte-Carlo algorithm. .. GENERATED FROM PYTHON SOURCE LINES 75-91 .. 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 92-94 Create the algorithms --------------------- .. GENERATED FROM PYTHON SOURCE LINES 97-112 .. 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 113-122 .. 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 123-206 .. 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 207-209 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 211-216 .. code-block:: Python gprResult = createMyBasicGPR(X, Y) graph = plotMyBasicGPR(gprResult, xMin, xMax, X, Y, event, sampleX, probability) view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_002.svg :alt: Estimated probability = 0.017606, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_002.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 217-219 Active learning Gaussian Process Regressor to sequentially add new points ------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 221-225 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 228-246 .. 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 247-248 We first call `getNewPoint` to get a point to add to the design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 251-254 .. 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 255-256 Then we evaluate the function on the new point and add it to the training design of experiments. .. GENERATED FROM PYTHON SOURCE LINES 258-262 .. code-block:: Python yNew = g(xNew) X.add(xNew) Y.add(yNew) .. GENERATED FROM PYTHON SOURCE LINES 263-264 We now plot the updated Gaussian Process Regressor. .. GENERATED FROM PYTHON SOURCE LINES 266-270 .. code-block:: Python gprResult = createMyBasicGPR(X, Y) graph = plotMyBasicGPR(gprResult, xMin, xMax, X, Y, event, sampleX, probability) view = viewer.View(graph) .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_003.svg :alt: Estimated probability = 0.000000, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_003.svg :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 271-272 The algorithm added in the domain. .. GENERATED FROM PYTHON SOURCE LINES 274-283 .. 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) View(graph) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_004.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_004.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_005.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_005.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_006.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_006.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_007.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_007.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_008.svg :alt: Estimated probability = 0.000339, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_008.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_009.svg :alt: Estimated probability = 0.005248, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_009.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_010.svg :alt: Estimated probability = 0.012020, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_010.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_011.svg :alt: Estimated probability = 0.016083, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_011.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_012.svg :alt: Estimated probability = 0.016760, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_012.svg :class: sphx-glr-multi-img * .. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_013.svg :alt: Estimated probability = 0.016760, Reference probability = 0.016760 :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_gpr_active_learning_013.svg :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 284-288 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 290-300 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 300-302 .. code-block:: Python View.ShowAll() .. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_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 `