.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_reliability_sensitivity/reliability/plot_line_sampling.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_reliability_sensitivity_reliability_plot_line_sampling.py: Estimate a probability using Line Sampling ========================================== .. GENERATED FROM PYTHON SOURCE LINES 7-8 In this example, we estimate the probability that the output of a function exceeds a given threshold with the Line Sampling method. .. GENERATED FROM PYTHON SOURCE LINES 10-15 .. code-block:: Python import openturns as ot import openturns.experimental as otexp import openturns.viewer as otv import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 16-18 Define the limit state function and the random vector ----------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 18-27 .. code-block:: Python dim = 2 g_twoBranch = ot.SymbolicFunction( ["x1", "x2"], [ "min(4 + 0.1 * (x1 - x2)^2 - (x1 + x2) / sqrt(2), 4 + 0.1 * (x1 - x2)^2 + (x1 + x2) / sqrt(2))" ], ) X = ot.RandomVector(ot.Normal(dim)) .. GENERATED FROM PYTHON SOURCE LINES 28-30 Define the event ---------------- .. GENERATED FROM PYTHON SOURCE LINES 30-34 .. code-block:: Python Y_twoBranch = ot.CompositeRandomVector(g_twoBranch, X) threshold = 1.5 event_twoBranch = ot.ThresholdEvent(Y_twoBranch, ot.Less(), threshold) .. GENERATED FROM PYTHON SOURCE LINES 35-37 Run FORM approximation ---------------------- .. GENERATED FROM PYTHON SOURCE LINES 37-44 .. code-block:: Python optimAlgo = ot.Cobyla() optimAlgo.setStartingPoint(X.getMean()) algo = ot.FORM(optimAlgo, event_twoBranch) algo.run() resultFORM = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 45-47 Run Line Sampling ----------------- .. GENERATED FROM PYTHON SOURCE LINES 49-50 The design point will define the initial important direction .. GENERATED FROM PYTHON SOURCE LINES 50-52 .. code-block:: Python alpha_twoBranch = resultFORM.getStandardSpaceDesignPoint() .. GENERATED FROM PYTHON SOURCE LINES 53-54 Define the root solver .. GENERATED FROM PYTHON SOURCE LINES 54-57 .. code-block:: Python solver = ot.Brent(1e-3, 1e-3, 1e-3, 5) rootStrategy = ot.MediumSafe(solver) .. GENERATED FROM PYTHON SOURCE LINES 58-59 Create the algorithm .. GENERATED FROM PYTHON SOURCE LINES 59-65 .. code-block:: Python algo = otexp.LineSampling(event_twoBranch, alpha_twoBranch, rootStrategy) algo.setMaximumOuterSampling(1000) algo.setMaximumCoefficientOfVariation(5e-2) # disable adaptive important direction to make plots easier to interpret algo.setAdaptiveImportantDirection(False) .. GENERATED FROM PYTHON SOURCE LINES 66-67 Save the important direction history, and especially all root points .. GENERATED FROM PYTHON SOURCE LINES 67-69 .. code-block:: Python algo.setStoreHistory(True) .. GENERATED FROM PYTHON SOURCE LINES 70-71 Run the simulation .. GENERATED FROM PYTHON SOURCE LINES 71-75 .. code-block:: Python algo.run() result = algo.getResult() print(result) .. rst-class:: sphx-glr-script-out .. code-block:: none probabilityEstimate=8.191688e-03 varianceEstimate=1.672832e-07 standard deviation=4.09e-04 coefficient of variation=4.99e-02 confidenceLength(0.95)=1.60e-03 outerSampling=81 blockSize=1 .. GENERATED FROM PYTHON SOURCE LINES 76-77 Retrieve important directions .. GENERATED FROM PYTHON SOURCE LINES 77-79 .. code-block:: Python alphas = algo.getAlphaHistory() .. GENERATED FROM PYTHON SOURCE LINES 80-81 Retrieve root points and root values .. GENERATED FROM PYTHON SOURCE LINES 81-85 .. code-block:: Python rootPoints = algo.getRootPointsHistory() rootValues = algo.getRootValuesHistory() .. GENERATED FROM PYTHON SOURCE LINES 86-119 .. code-block:: Python def drawLines(algo, n=10): """Draw sampling lines and the corresponding roots.""" rootPoints = algo.getRootPointsHistory() alphas = algo.getAlphaHistory() n = min(n, len(rootPoints)) for i in range(n): # there can be several roots per sample n_roots = len(rootPoints[i]) alpha = alphas[i] for j in range(n_roots): if i + j == 0: print(f"n_roots={n_roots}") dp = rootPoints[i][j] # project design point on the hyperplane orthogonal to alpha to get origin point uPoint = dp uDot = uPoint.dot(alpha) uPoint = uPoint - alpha * uDot # draw segment origin -> design point plt.plot( (uPoint[0], dp[0]), (uPoint[1], dp[1]), color="blue", linestyle="--", linewidth=0.75, ) # draw origin plt.plot(uPoint[0], uPoint[1], "ro", markersize=3) # draw design point plt.plot(dp[0], dp[1], "bo", markersize=3, zorder=3) .. GENERATED FROM PYTHON SOURCE LINES 120-144 .. code-block:: Python def drawLSDesign(algo): """Draw sampling lines, roots, and the limit state function.""" dmin = [-4.0] * 2 dmax = [4.0] * 2 graph1 = g_twoBranch.draw(dmin, dmax) contour_g = graph1.getDrawable(0).getImplementation() contour_g.setColorBarPosition("") contour_g.setColorMap("gray") graph1.setDrawable(contour_g, 0) view1 = otv.View(graph1, square_axes=True) # now draw the the limit state on same figure graph2 = g_twoBranch.draw(dmin, dmax) contour_g = graph2.getDrawable(0).getImplementation() contour_g.setLevels([threshold]) contour_g.setColor("red") graph2.setDrawable(contour_g, 0) plt.axline([-1, 1], [1, -1], linestyle="dotted", color="gray") drawLines(algo) graph2.setTitle("Line Sampling") _ = otv.View(graph2, figure=view1.getFigure()) .. GENERATED FROM PYTHON SOURCE LINES 145-146 Plot the limit state, a few design points and their origin .. GENERATED FROM PYTHON SOURCE LINES 146-148 .. code-block:: Python drawLSDesign(algo) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_line_sampling_001.svg :alt: Line Sampling :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_line_sampling_001.svg :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none n_roots=2 .. GENERATED FROM PYTHON SOURCE LINES 149-151 Now we disable the opposite direction search to try to save function evaluations. This practice is incorrect in this case, however, as the algorithm misses half the probability of the event. .. GENERATED FROM PYTHON SOURCE LINES 151-157 .. code-block:: Python ot.RandomGenerator.SetSeed(0) algo.setSearchOppositeDirection(False) algo.run() result = algo.getResult() print(result) .. rst-class:: sphx-glr-script-out .. code-block:: none probabilityEstimate=4.095844e-03 varianceEstimate=4.182085e-08 standard deviation=2.05e-04 coefficient of variation=4.99e-02 confidenceLength(0.95)=8.02e-04 outerSampling=81 blockSize=1 .. GENERATED FROM PYTHON SOURCE LINES 158-160 We plot the limit state, a few design points and their origin. This time we see each origin sample point in the orthogonal hyperplane yields only one design point. .. GENERATED FROM PYTHON SOURCE LINES 160-162 .. code-block:: Python drawLSDesign(algo) .. image-sg:: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_line_sampling_002.svg :alt: Line Sampling :srcset: /auto_reliability_sensitivity/reliability/images/sphx_glr_plot_line_sampling_002.svg :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none n_roots=1 .. GENERATED FROM PYTHON SOURCE LINES 163-164 Now re-enable adaptive important direction search .. GENERATED FROM PYTHON SOURCE LINES 164-171 .. code-block:: Python ot.RandomGenerator.SetSeed(0) algo.setAdaptiveImportantDirection(True) algo.setSearchOppositeDirection(True) algo.run() result = algo.getResult() print(result) .. rst-class:: sphx-glr-script-out .. code-block:: none probabilityEstimate=8.253432e-03 varianceEstimate=1.683218e-07 standard deviation=4.10e-04 coefficient of variation=4.97e-02 confidenceLength(0.95)=1.61e-03 outerSampling=80 blockSize=1 .. GENERATED FROM PYTHON SOURCE LINES 172-173 Inspect important directions (without duplicates from history) .. GENERATED FROM PYTHON SOURCE LINES 173-180 .. code-block:: Python unique_alphas = ot.Sample(0, dim) for alpha in algo.getAlphaHistory(): if len(unique_alphas) == 0 or alpha != unique_alphas[-1]: unique_alphas.add(alpha) print("unique alphas:") print(unique_alphas) .. rst-class:: sphx-glr-script-out .. code-block:: none unique alphas: 0 : [ 0.707104 0.707109 ] 1 : [ 0.939276 0.343163 ] 2 : [ 0.634607 0.772835 ] 3 : [ 0.671819 0.740716 ] 4 : [ -0.689475 -0.724309 ] 5 : [ 0.71885 0.695165 ] 6 : [ -0.695049 -0.718963 ] 7 : [ 0.69544 0.718584 ] .. GENERATED FROM PYTHON SOURCE LINES 181-182 .. code-block:: Python otv.View.ShowAll() .. _sphx_glr_download_auto_reliability_sensitivity_reliability_plot_line_sampling.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_line_sampling.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_line_sampling.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_line_sampling.zip `