.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_reliability_sensitivity/reliability_processes/plot_field_fca_sobol.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_processes_plot_field_fca_sobol.py: Estimate Sobol indices on a field to point function =================================================== .. GENERATED FROM PYTHON SOURCE LINES 7-26 In this example, we are going to perform sensitivity analysis of an application that takes fields as input and vectors as output from a sample of data: .. math:: h: \left| \begin{array}{ccl} \cM_N \times (\Rset^d)^N & \rightarrow & \Rset^p \\ \mat{X} & \mapsto & \vect{Y} \end{array} \right. This involves these steps: - Generate some input/output data matching the application :math:`h` - Run the :class:`~openturns.experimental.FieldToPointFunctionalChaosAlgorithm` class - Validate the Karhunen-Loeve decompositions of the inputs - Validate the chaos metamodel between the KL coefficients and the outputs - Retrieve the Sobol' indices from :class:`~openturns.FieldFunctionalChaosSobolIndices` .. GENERATED FROM PYTHON SOURCE LINES 28-32 .. code-block:: Python import openturns as ot import openturns.experimental as otexp from openturns.viewer import View .. GENERATED FROM PYTHON SOURCE LINES 33-35 First build a process to generate the input data. We assemble a 4-d process from functional and Gaussian processes. .. GENERATED FROM PYTHON SOURCE LINES 35-50 .. code-block:: Python T = 3.0 NT = 32 tg = ot.RegularGrid(0.0, T / NT, NT) f1 = ot.SymbolicFunction(["t"], ["sin(t)"]) f2 = ot.SymbolicFunction(["t"], ["cos(t)^2"]) coeff1_dist = ot.Normal([1.0] * 2, [0.6] * 2, ot.CorrelationMatrix(2)) p1 = ot.FunctionalBasisProcess(coeff1_dist, ot.Basis([f1, f2]), tg) p2 = ot.GaussianProcess(ot.SquaredExponential([1.0], [T / 4.0]), tg) coeff3_dist = ot.ComposedDistribution([ot.Uniform(), ot.Normal()]) f1 = ot.SymbolicFunction(["t"], ["1", "0"]) f2 = ot.SymbolicFunction(["t"], ["0", "1"]) p3 = ot.FunctionalBasisProcess(coeff3_dist, ot.Basis([f1, f2])) X = ot.AggregatedProcess([p1, p2, p3]) X.setMesh(tg) .. GENERATED FROM PYTHON SOURCE LINES 51-52 Draw some input trajectories from our process .. GENERATED FROM PYTHON SOURCE LINES 52-58 .. code-block:: Python ot.RandomGenerator.SetSeed(0) x = X.getSample(10) graph = x.drawMarginal(0) graph.setTitle(f"{x.getSize()} input trajectories") _ = View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_001.png :alt: 10 input trajectories :srcset: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 59-60 Generate input realizations and the corresponding output from a Field->Point function .. GENERATED FROM PYTHON SOURCE LINES 60-80 .. code-block:: Python class pyf2p(ot.OpenTURNSPythonFieldToPointFunction): def __init__(self, mesh): super(pyf2p, self).__init__(mesh, 4, 1) self.setInputDescription(["x1", "x2", "x3", "x4"]) self.setOutputDescription(["y"]) def _exec(self, X): Xs = ot.Sample(X) x1, x2, x3, x4 = Xs.computeMean() y = x1 + x2 + x3 - x4 + x1 * x2 - x3 * x4 - 0.1 * x1 * x2 * x3 return [y] f = ot.FieldToPointFunction(pyf2p(tg)) N = 1000 x = X.getSample(N) y = f(x) .. GENERATED FROM PYTHON SOURCE LINES 81-83 Run the field-vector algorithm that performs KL-decomposition of the inputs and chaos learning between the KL coefficients and the output vectors .. GENERATED FROM PYTHON SOURCE LINES 83-101 .. code-block:: Python algo = otexp.FieldToPointFunctionalChaosAlgorithm(x, y) # 1. KL parameters algo.setCenteredSample(False) # our input sample is not centered (default) algo.setThreshold(4e-2) # we expect to explain 96% of variance algo.setRecompress( False ) # whether to re-truncate modes according to a global eigen value threshold across inputs (default) algo.setNbModes(10) # max KL modes (default=unlimited) # 2. chaos parameters: bs = ot.ResourceMap.GetAsUnsignedInteger("FunctionalChaosAlgorithm-BasisSize") ot.ResourceMap.SetAsUnsignedInteger( "FunctionalChaosAlgorithm-BasisSize", N ) # chaos basis size ot.ResourceMap.SetAsBool("FunctionalChaosAlgorithm-Sparse", True) algo.run() ot.ResourceMap.SetAsUnsignedInteger("FunctionalChaosAlgorithm-BasisSize", bs) result = algo.getResult() .. GENERATED FROM PYTHON SOURCE LINES 102-104 Retrieve the eigen values of each KL decomposition: we observe that each input process is represented by a different number of modes. .. GENERATED FROM PYTHON SOURCE LINES 104-108 .. code-block:: Python kl_results = result.getInputKLResultCollection() n_modes = [len(res.getEigenvalues()) for res in kl_results] print(f"n_modes={n_modes}") .. rst-class:: sphx-glr-script-out .. code-block:: none n_modes=[2, 3, 1, 1] .. GENERATED FROM PYTHON SOURCE LINES 109-111 Retrieve the ratios of selected variance over cumulated variance: we see that all 3 inputs are perfectly represented, and the 2nd input almost perfectly. .. GENERATED FROM PYTHON SOURCE LINES 111-114 .. code-block:: Python ratios = [res.getSelectionRatio() for res in kl_results] print(f"ratios={ratios}") .. rst-class:: sphx-glr-script-out .. code-block:: none ratios=[1.0, 0.9851877006609377, 1.0, 1.0] .. GENERATED FROM PYTHON SOURCE LINES 115-118 Graphically validate the KL decompositions: we also see that the 2nd input appear to be less well represented than the others. Note however that the 0.98 selected/cumulated variance ratio actually means it is very good. .. GENERATED FROM PYTHON SOURCE LINES 118-126 .. code-block:: Python graphs = [] for i in range(x.getDimension()): validation = ot.KarhunenLoeveValidation(x.getMarginal(i), kl_results[i]) graph = validation.drawValidation().getGraph(0, 0) graph.setTitle(f"KL validation - marginal #{i} ratio={100.0 * ratios[i]:.2f} %") View(graph) graphs.append(graph) .. rst-class:: sphx-glr-horizontal * .. image-sg:: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_002.png :alt: KL validation - marginal #0 ratio=100.00 % :srcset: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_002.png :class: sphx-glr-multi-img * .. image-sg:: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_003.png :alt: KL validation - marginal #1 ratio=98.52 % :srcset: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_003.png :class: sphx-glr-multi-img * .. image-sg:: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_004.png :alt: KL validation - marginal #2 ratio=100.00 % :srcset: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_004.png :class: sphx-glr-multi-img * .. image-sg:: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_005.png :alt: KL validation - marginal #3 ratio=100.00 % :srcset: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_005.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 127-129 On the 2nd marginal we can filter out the points inside the 99% level-set to see that actually only a few points out of N are actually outliers. .. GENERATED FROM PYTHON SOURCE LINES 129-150 .. code-block:: Python graph = graphs[1] data = graph.getDrawable(1).getData() normal = ot.NormalFactory().build(data) log_pdf = normal.computeLogPDF(data).asPoint() l_pair = [(log_pdf[i], data[i]) for i in range(len(data))] l_pair.sort(key=lambda t: t[0]) index_bad = int(0.01 * len(data)) # here 0.01 = (100-99)% beta = l_pair[index_bad][0] gnorm = normal.drawLogPDF(data.getMin(), data.getMax()) bad = [l_pair[i][1] for i in range(index_bad + 1)] c = ot.Cloud(bad) c.setPointStyle("bullet") c.setColor("blue") graph.setDrawable(c, 1) dr = gnorm.getDrawable(0) dr.setLevels([beta]) dr.setColor("red") dr.setLegend("99% level-set") graph.add(dr) _ = View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_006.png :alt: KL validation - marginal #1 ratio=98.52 % :srcset: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 151-153 Inspect the chaos quality: residuals and relative errors. The relative error is very low; that means the chaos decomposition performs very well. .. GENERATED FROM PYTHON SOURCE LINES 153-156 .. code-block:: Python print(f"residuals={result.getFCEResult().getResiduals()}") print(f"relative errors={result.getFCEResult().getRelativeErrors()}") .. rst-class:: sphx-glr-script-out .. code-block:: none residuals=[0.000223921] relative errors=[1.44099e-08] .. GENERATED FROM PYTHON SOURCE LINES 157-160 Graphically validate the chaos result: we can see the points are very close to the diagonal; this means approximated points are very close to the learning points. .. GENERATED FROM PYTHON SOURCE LINES 160-170 .. code-block:: Python modes = result.getModesSample() metamodel = result.getFCEResult().getMetaModel() output = result.getOutputSample() validation = ot.MetaModelValidation(modes, output, metamodel) q2 = validation.computePredictivityFactor() print(f"q2={q2}") graph = validation.drawValidation() graph.setTitle(f"Chaos validation - q2={q2}") _ = View(graph) .. image-sg:: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_007.png :alt: Chaos validation - q2=[0.999986] :srcset: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none q2=[0.999986] .. GENERATED FROM PYTHON SOURCE LINES 171-173 Perform an evaluation on a new realization and ensure the output is close to the evaluation with the reference function .. GENERATED FROM PYTHON SOURCE LINES 173-179 .. code-block:: Python metamodel = result.getFieldToPointMetamodel() x0 = X.getRealization() y0 = f(x0) y0hat = metamodel(x0) print(f"y0={y0} y0^={y0hat}") .. rst-class:: sphx-glr-script-out .. code-block:: none y0=[6.01011] y0^=[6.00689] .. GENERATED FROM PYTHON SOURCE LINES 180-182 Retrieve the first order Sobol' indices The preponderant variables are x2, x4 whereas x1, x3 have a low influence on the output .. GENERATED FROM PYTHON SOURCE LINES 182-186 .. code-block:: Python sensitivity = otexp.FieldFunctionalChaosSobolIndices(result) sobol_0 = sensitivity.getFirstOrderIndices() print(f"first order={sobol_0}") .. rst-class:: sphx-glr-script-out .. code-block:: none first order=[0.066632,0.441135,0.0953934,0.275428] .. GENERATED FROM PYTHON SOURCE LINES 187-190 Retrieve the total Sorder obol' indices The x3,x4 variables have total order indices significantly different than their first order indices counterpart meaning they interact with other variables .. GENERATED FROM PYTHON SOURCE LINES 190-193 .. code-block:: Python sobol_0t = sensitivity.getTotalOrderIndices() print(f"total order={sobol_0t}") .. rst-class:: sphx-glr-script-out .. code-block:: none total order=[0.0902712,0.465189,0.193238,0.372786] .. GENERATED FROM PYTHON SOURCE LINES 194-195 Draw the Sobol' indices .. GENERATED FROM PYTHON SOURCE LINES 195-200 .. code-block:: Python graph = sensitivity.draw() view = View(graph) View.ShowAll() .. image-sg:: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_008.png :alt: Sobol' indices :srcset: /auto_reliability_sensitivity/reliability_processes/images/sphx_glr_plot_field_fca_sobol_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 201-202 Reset default settings .. GENERATED FROM PYTHON SOURCE LINES 202-203 .. code-block:: Python ot.ResourceMap.Reload() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 6.487 seconds) .. _sphx_glr_download_auto_reliability_sensitivity_reliability_processes_plot_field_fca_sobol.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_field_fca_sobol.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_field_fca_sobol.py `