.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_meta_modeling/polynomial_chaos_metamodel/plot_functional_chaos.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_meta_modeling_polynomial_chaos_metamodel_plot_functional_chaos.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_meta_modeling_polynomial_chaos_metamodel_plot_functional_chaos.py:


Create a polynomial chaos metamodel from a data set
===================================================

.. GENERATED FROM PYTHON SOURCE LINES 6-11

In this example, we create a polynomial chaos expansion (PCE) using
a data set.
More precisely, given a pair of input and output samples,
we create a PCE without the knowledge of the input distribution.
In this example, we use a relatively small sample size equal to 80.

.. GENERATED FROM PYTHON SOURCE LINES 13-34

In this example we create a global approximation of a model response using
polynomial chaos expansion.

Let :math:`h` be the function defined by:

.. math::
   g(\mathbf{x}) = \left[\cos(x_1 + x_2), (x_2 + 1) e^{x_1}\right]


for any :math:`\mathbf{x}\in\mathbb{R}^2`.

We assume that

.. math::
   X_1 \sim \mathcal{N}(0,1) \textrm{ and } X_2 \sim \mathcal{N}(0,1)

and that :math:`X_1` and :math:`X_2` are independent.

An interesting point in this example is that the output is multivariate.
This is why we are going to use the `getMarginal` method in the script
in order to select the output marginal that we want to manage.

.. GENERATED FROM PYTHON SOURCE LINES 36-38

Simulate input and output samples
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. GENERATED FROM PYTHON SOURCE LINES 40-46

.. code-block:: Python

    import openturns as ot
    import openturns.viewer as viewer
    from matplotlib import pylab as plt

    ot.Log.Show(ot.Log.NONE)








.. GENERATED FROM PYTHON SOURCE LINES 47-48

We first create the function `model`.

.. GENERATED FROM PYTHON SOURCE LINES 50-56

.. code-block:: Python

    ot.RandomGenerator.SetSeed(2)
    dimension = 2
    input_names = ["x1", "x2"]
    formulas = ["cos(x1 + x2)", "(x2 + 1) * exp(x1)"]
    model = ot.SymbolicFunction(input_names, formulas)








.. GENERATED FROM PYTHON SOURCE LINES 57-59

Then we create a sample `inputSample` and compute the corresponding output
sample `outputSample`.

.. GENERATED FROM PYTHON SOURCE LINES 61-66

.. code-block:: Python

    distribution = ot.Normal(dimension)
    samplesize = 80
    inputSample = distribution.getSample(samplesize)
    outputSample = model(inputSample)








.. GENERATED FROM PYTHON SOURCE LINES 67-69

Create the PCE
~~~~~~~~~~~~~~

.. GENERATED FROM PYTHON SOURCE LINES 71-81

Create a functional chaos model.
The algorithm needs to fit a distribution on the input sample.
To do this, the algorithm in :class:`~openturns.FunctionalChaosAlgorithm`
uses the :class:`~openturns.FunctionalChaosAlgorithm.BuildDistribution`
static method to fit the distribution to the input sample.
Please read :doc:`Fit a distribution from an input sample </auto_meta_modeling/polynomial_chaos_metamodel/plot_chaos_build_distribution>`
for an example of this method.
The algorithm does this automatically using the Lilliefors test.
In order to make the algorithm a little faster, we reduce the
value of the sample size used in the Lilliefors test.

.. GENERATED FROM PYTHON SOURCE LINES 83-85

.. code-block:: Python

    ot.ResourceMap.SetAsUnsignedInteger("FittingTest-LillieforsMaximumSamplingSize", 50)








.. GENERATED FROM PYTHON SOURCE LINES 86-89

The main topic of this example is to introduce the next constructor of
:class:`~openturns.FunctionalChaosAlgorithm`.
Notice that the only input arguments are the input and output sample.

.. GENERATED FROM PYTHON SOURCE LINES 89-94

.. code-block:: Python

    algo = ot.FunctionalChaosAlgorithm(inputSample, outputSample)
    algo.run()
    result = algo.getResult()
    result






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    FunctionalChaosResult
    <ul>
      <li>input dimension: 2</li>
      <li>output dimension: 2</li>
      <li>distribution dimension: 2</li>
      <li>transformation: 2 -> 2</li>
      <li>inverse transformation: 2 -> 2</li>
      <li>orthogonal basis dimension: 2</li>
      <li>indices size: 33</li>
      <li>relative errors: [0.000423845,3.97704e-08]</li>
      <li>residuals: [0.00156354,0.000147243]</li>
    </ul>
    <table>
      <tr>
        <th>Index</th>
        <th>Multi-index</th>
        <th>Coeff.#0</th>
        <th>Coeff.#1</th>
      </tr>
      <tr>
        <th>0</th>
        <td>[0,0]</td>
        <td>-0.9778362</td>
        <td>2.494928</td>
      </tr>
      <tr>
        <th>1</th>
        <td>[1,0]</td>
        <td>3.946319</td>
        <td>1.310505</td>
      </tr>
      <tr>
        <th>2</th>
        <td>[0,1]</td>
        <td>-0.5668731</td>
        <td>2.058245</td>
      </tr>
      <tr>
        <th>3</th>
        <td>[2,0]</td>
        <td>-6.595425</td>
        <td>3.31883</td>
      </tr>
      <tr>
        <th>4</th>
        <td>[0,2]</td>
        <td>0.05754212</td>
        <td>0.1086024</td>
      </tr>
      <tr>
        <th>5</th>
        <td>[3,0]</td>
        <td>7.262131</td>
        <td>-0.841792</td>
      </tr>
      <tr>
        <th>6</th>
        <td>[0,3]</td>
        <td>-0.9758324</td>
        <td>-0.12278</td>
      </tr>
      <tr>
        <th>7</th>
        <td>[1,1]</td>
        <td>-0.7887584</td>
        <td>2.418963</td>
      </tr>
      <tr>
        <th>8</th>
        <td>[4,0]</td>
        <td>-6.489266</td>
        <td>1.935301</td>
      </tr>
      <tr>
        <th>9</th>
        <td>[0,4]</td>
        <td>0.2844189</td>
        <td>0.1285879</td>
      </tr>
      <tr>
        <th>10</th>
        <td>[5,0]</td>
        <td>5.352722</td>
        <td>-1.168838</td>
      </tr>
      <tr>
        <th>11</th>
        <td>[0,5]</td>
        <td>-0.8491851</td>
        <td>-0.1064447</td>
      </tr>
      <tr>
        <th>12</th>
        <td>[2,1]</td>
        <td>0.3063163</td>
        <td>1.643328</td>
      </tr>
      <tr>
        <th>13</th>
        <td>[1,2]</td>
        <td>-0.0281957</td>
        <td>0.01273302</td>
      </tr>
      <tr>
        <th>14</th>
        <td>[6,0]</td>
        <td>-3.800862</td>
        <td>0.9235032</td>
      </tr>
      <tr>
        <th>15</th>
        <td>[0,6]</td>
        <td>0.004040826</td>
        <td>0.1013839</td>
      </tr>
      <tr>
        <th>16</th>
        <td>[7,0]</td>
        <td>2.402248</td>
        <td>-0.5406551</td>
      </tr>
      <tr>
        <th>17</th>
        <td>[0,7]</td>
        <td>-0.441342</td>
        <td>-0.0634571</td>
      </tr>
      <tr>
        <th>18</th>
        <td>[3,1]</td>
        <td>-0.003325215</td>
        <td>1.022703</td>
      </tr>
      <tr>
        <th>19</th>
        <td>[1,3]</td>
        <td>-0.2267971</td>
        <td>0.002919773</td>
      </tr>
      <tr>
        <th>20</th>
        <td>[2,2]</td>
        <td>0.4659559</td>
        <td>-0.001704627</td>
      </tr>
      <tr>
        <th>21</th>
        <td>[8,0]</td>
        <td>-1.272689</td>
        <td>0.2902903</td>
      </tr>
      <tr>
        <th>22</th>
        <td>[0,8]</td>
        <td>-0.03341318</td>
        <td>0.05202352</td>
      </tr>
      <tr>
        <th>23</th>
        <td>[4,1]</td>
        <td>0.1240007</td>
        <td>0.2562661</td>
      </tr>
      <tr>
        <th>24</th>
        <td>[1,4]</td>
        <td>-0.04681154</td>
        <td>-0.006325388</td>
      </tr>
      <tr>
        <th>25</th>
        <td>[9,0]</td>
        <td>0.5059943</td>
        <td>-0.1094002</td>
      </tr>
      <tr>
        <th>26</th>
        <td>[0,9]</td>
        <td>-0.1155697</td>
        <td>-0.01848551</td>
      </tr>
      <tr>
        <th>27</th>
        <td>[3,2]</td>
        <td>-0.04218489</td>
        <td>0.02369182</td>
      </tr>
      <tr>
        <th>28</th>
        <td>[2,3]</td>
        <td>-0.05014101</td>
        <td>-0.01482048</td>
      </tr>
      <tr>
        <th>29</th>
        <td>[10,0]</td>
        <td>-0.1825852</td>
        <td>0.03945734</td>
      </tr>
      <tr>
        <th>30</th>
        <td>[0,10]</td>
        <td>-0.01467434</td>
        <td>0.01355966</td>
      </tr>
      <tr>
        <th>31</th>
        <td>[5,1]</td>
        <td>-0.08889085</td>
        <td>0.1317192</td>
      </tr>
      <tr>
        <th>32</th>
        <td>[1,5]</td>
        <td>-0.2461735</td>
        <td>0.001139069</td>
      </tr>
    </table>

    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 95-100

Not all coefficients are selected in this PCE.
Indeed, the default constructor of :class:`~openturns.FunctionalChaosAlgorithm`
creates a sparse PCE.
Please read :doc:`Create a full or sparse polynomial chaos expansion </auto_meta_modeling/polynomial_chaos_metamodel/plot_functional_chaos_database>`
for more details on this topic.

.. GENERATED FROM PYTHON SOURCE LINES 102-103

Get the metamodel.

.. GENERATED FROM PYTHON SOURCE LINES 103-105

.. code-block:: Python

    metamodel = result.getMetaModel()








.. GENERATED FROM PYTHON SOURCE LINES 106-109

Plot the second output of our model depending on :math:`x_2` with :math:`x_1=0.5`.
In order to do this, we create a `ParametricFunction` and set the value of :math:`x_1`.
Then we use the `getMarginal` method to extract the second output (which index is equal to 1).

.. GENERATED FROM PYTHON SOURCE LINES 111-129

.. code-block:: Python

    x1index = 0
    x1value = 0.5
    x2min = -3.0
    x2max = 3.0
    outputIndex = 1
    metamodelParametric = ot.ParametricFunction(metamodel, [x1index], [x1value])
    graph = metamodelParametric.getMarginal(outputIndex).draw(x2min, x2max)
    graph.setLegends(["Metamodel"])
    modelParametric = ot.ParametricFunction(model, [x1index], [x1value])
    curve = modelParametric.getMarginal(outputIndex).draw(x2min, x2max).getDrawable(0)
    curve.setColor("red")
    curve.setLegend("Model")
    graph.add(curve)
    graph.setLegendPosition("lower right")
    graph.setXTitle("X2")
    graph.setTitle("Metamodel Validation, output #%d" % (outputIndex))
    view = viewer.View(graph)




.. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_functional_chaos_001.png
   :alt: Metamodel Validation, output #1
   :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_functional_chaos_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 130-134

We see that the metamodel fits approximately to the model, except
perhaps for extreme values of :math:`x_2`.
However, there is a better way of globally validating the metamodel,
using the :class:`~openturns.MetaModelValidation` on a validation design of experiment.

.. GENERATED FROM PYTHON SOURCE LINES 136-141

.. code-block:: Python

    n_valid = 100
    inputTest = distribution.getSample(n_valid)
    outputTest = model(inputTest)









.. GENERATED FROM PYTHON SOURCE LINES 142-143

Plot the corresponding validation graphics.

.. GENERATED FROM PYTHON SOURCE LINES 145-151

.. code-block:: Python

    val = ot.MetaModelValidation(outputTest, metamodel(inputTest))
    R2 = val.computeR2Score()
    graph = val.drawValidation()
    graph.setTitle("Metamodel validation R2=" + str(R2))
    view = viewer.View(graph)




.. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_functional_chaos_002.png
   :alt: Metamodel validation R2=[0.486703,0.999523]
   :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_functional_chaos_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 152-156

The coefficient of predictivity is not extremely satisfactory for the
first output, but is would be sufficient for a central dispersion study.
The second output has a much more satisfactory R2: only one single
extreme point is far from the diagonal of the graphics.

.. GENERATED FROM PYTHON SOURCE LINES 158-160

Compute and print Sobol' indices
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. GENERATED FROM PYTHON SOURCE LINES 162-165

.. code-block:: Python

    chaosSI = ot.FunctionalChaosSobolIndices(result)
    chaosSI






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    FunctionalChaosSobolIndices
    <ul>
      <li>input dimension: 2</li>
      <li>output dimension: 2</li>
      <li>basis size: 33</li>
      <li>mean: [-0.977836,2.49493]</li>
      <li>std-dev: [14.4244,5.81255]</li>
    </ul>
    Output marginal: 0
    <table>
      <tr>
        <th>Input</th>
        <th>Variable</th>
        <th>Sobol' index</th>
        <th>Total index</th>
      </tr>
      <tr>
        <td>0</td>
        <td>X0</td>
        <td>0.983831</td>
        <td>0.989001</td>
      </tr>
      <tr>
        <td>1</td>
        <td>X1</td>
        <td>0.010999</td>
        <td>0.016169</td>
      </tr>
    </table>
    <table>
      <tr>
        <th>Index</th>
        <th>Multi-index</th>
        <th>Part of variance</th>
      </tr>
      <tr>
        <td>5</td>
        <td>[3,0]</td>
        <td>0.253472</td>
      </tr>
      <tr>
        <td>3</td>
        <td>[2,0]</td>
        <td>0.209068</td>
      </tr>
      <tr>
        <td>8</td>
        <td>[4,0]</td>
        <td>0.202392</td>
      </tr>
      <tr>
        <td>10</td>
        <td>[5,0]</td>
        <td>0.137706</td>
      </tr>
      <tr>
        <td>1</td>
        <td>[1,0]</td>
        <td>0.074849</td>
      </tr>
      <tr>
        <td>14</td>
        <td>[6,0]</td>
        <td>0.069433</td>
      </tr>
      <tr>
        <td>16</td>
        <td>[7,0]</td>
        <td>0.027736</td>
      </tr>
    </table>
    Output marginal: 1
    <table>
      <tr>
        <th>Input</th>
        <th>Variable</th>
        <th>Sobol' index</th>
        <th>Total index</th>
      </tr>
      <tr>
        <td>0</td>
        <td>X0</td>
        <td>0.585905</td>
        <td>0.872471</td>
      </tr>
      <tr>
        <td>1</td>
        <td>X1</td>
        <td>0.127529</td>
        <td>0.414095</td>
      </tr>
    </table>
    <table>
      <tr>
        <th>Index</th>
        <th>Multi-index</th>
        <th>Part of variance</th>
      </tr>
      <tr>
        <td>3</td>
        <td>[2,0]</td>
        <td>0.326015</td>
      </tr>
      <tr>
        <td>7</td>
        <td>[1,1]</td>
        <td>0.173191</td>
      </tr>
      <tr>
        <td>2</td>
        <td>[0,1]</td>
        <td>0.125390</td>
      </tr>
      <tr>
        <td>8</td>
        <td>[4,0]</td>
        <td>0.110857</td>
      </tr>
      <tr>
        <td>12</td>
        <td>[2,1]</td>
        <td>0.079931</td>
      </tr>
      <tr>
        <td>1</td>
        <td>[1,0]</td>
        <td>0.050833</td>
      </tr>
      <tr>
        <td>10</td>
        <td>[5,0]</td>
        <td>0.040437</td>
      </tr>
      <tr>
        <td>18</td>
        <td>[3,1]</td>
        <td>0.030958</td>
      </tr>
      <tr>
        <td>14</td>
        <td>[6,0]</td>
        <td>0.025243</td>
      </tr>
      <tr>
        <td>5</td>
        <td>[3,0]</td>
        <td>0.020974</td>
      </tr>
    </table>

    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 166-174

Let us analyse the results of this global sensitivity analysis.

* We see that the first output involves significant multi-indices with
  higher marginal degree.
* For the second output, the first variable is especially significant,
  with a significant contribution of the interactions.
  The contribution of the interactions are very
  significant in this model.

.. GENERATED FROM PYTHON SOURCE LINES 176-177

Draw Sobol' indices.

.. GENERATED FROM PYTHON SOURCE LINES 179-183

.. code-block:: Python

    sensitivityAnalysis = ot.FunctionalChaosSobolIndices(result)
    first_order = [sensitivityAnalysis.getSobolIndex(i) for i in range(dimension)]
    total_order = [sensitivityAnalysis.getSobolTotalIndex(i) for i in range(dimension)]








.. GENERATED FROM PYTHON SOURCE LINES 184-189

.. code-block:: Python

    input_names = model.getInputDescription()
    graph = ot.SobolIndicesAlgorithm.DrawSobolIndices(input_names, first_order, total_order)
    graph.setLegendPosition("center")
    view = viewer.View(graph)




.. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_functional_chaos_003.png
   :alt: Sobol' indices
   :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_functional_chaos_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 190-205

Testing the sensitivity to the degree
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

With the specific constructor of `FunctionalChaosAlgorithm` that
we use, the `FunctionalChaosAlgorithm-MaximumTotalDegree`
in the `ResourceMap` configure the maximum degree explored by
the algorithm. This degree is a trade-off.

* If the maximum degree is too low, the polynomial may miss some
  coefficients so that the quality is lower than possible.
* If the maximum degree is too large, the number of coefficients
  to explore is too large, so that the coefficients might be poorly estimated.

This is why the following `for` loop explores various degrees to see
the sensitivity of the metamodel predictivity depending on the degree.

.. GENERATED FROM PYTHON SOURCE LINES 207-208

The default value of this parameter is 10.

.. GENERATED FROM PYTHON SOURCE LINES 210-212

.. code-block:: Python

    ot.ResourceMap.GetAsUnsignedInteger("FunctionalChaosAlgorithm-MaximumTotalDegree")





.. rst-class:: sphx-glr-script-out

 .. code-block:: none


    10



.. GENERATED FROM PYTHON SOURCE LINES 213-214

This is why we explore the values from 1 to 15.

.. GENERATED FROM PYTHON SOURCE LINES 216-234

.. code-block:: Python

    degrees = range(1, 12)
    r2 = ot.Sample(len(degrees), 2)
    for maximumDegree in degrees:
        ot.ResourceMap.SetAsUnsignedInteger(
            "FunctionalChaosAlgorithm-MaximumTotalDegree", maximumDegree
        )
        print("Maximum total degree =", maximumDegree)
        algo = ot.FunctionalChaosAlgorithm(inputSample, outputSample)
        algo.run()
        result = algo.getResult()
        metamodel = result.getMetaModel()
        for outputIndex in range(2):
            val = ot.MetaModelValidation(
                outputTest[:, outputIndex], metamodel.getMarginal(outputIndex)(inputTest)
            )
            r2Value = min(1.0, max(0.0, val.computeR2Score()[0]))  # Get lucky.
            r2[maximumDegree - degrees[0], outputIndex] = r2Value





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Maximum total degree = 1
    Maximum total degree = 2
    Maximum total degree = 3
    Maximum total degree = 4
    Maximum total degree = 5
    Maximum total degree = 6
    Maximum total degree = 7
    Maximum total degree = 8
    Maximum total degree = 9
    Maximum total degree = 10
    Maximum total degree = 11




.. GENERATED FROM PYTHON SOURCE LINES 235-251

.. code-block:: Python

    graph = ot.Graph("Predictivity", "Total degree", "R2", True)
    cloud = ot.Cloud([[d] for d in degrees], r2[:, 0])
    cloud.setLegend("Output #0")
    cloud.setPointStyle("bullet")
    graph.add(cloud)
    cloud = ot.Cloud([[d] for d in degrees], r2[:, 1])
    cloud.setLegend("Output #1")
    cloud.setColor("red")
    cloud.setPointStyle("bullet")
    graph.add(cloud)
    graph.setLegendPosition("upper right")
    view = viewer.View(graph, legend_kw={"bbox_to_anchor": (1.0, 1.0), "loc": "upper left"})
    plt.subplots_adjust(right=0.7)

    plt.show()




.. image-sg:: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_functional_chaos_004.png
   :alt: Predictivity
   :srcset: /auto_meta_modeling/polynomial_chaos_metamodel/images/sphx_glr_plot_functional_chaos_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 252-275

We see that the R2 score increases then gets constant or decreases.
A low total polynomial degree is not sufficient to describe
the first output with good predictivity.
However, the coefficient of predictivity can decrease when the total degree
gets greater.
The predictivity of the second output seems to be much less
satisfactory: a little more work would be required to improve the metamodel.

In this situation, the following methods may be used.

* Since the distribution of the input is known, we may want to give
  this information to the `FunctionalChaosAlgorithm`.
  This prevents the algorithm from trying to fit the input distribution
  which best fit to the data.
* We may want to customize the `adaptiveStrategy` by selecting an enumerate
  function which best fit to this particular example.
  In this specific example, however, the interactions plays a great role so that the
  linear enumerate function may provide better results than the hyperbolic rule.
* We may want to customize the `projectionStrategy` by selecting a method
  to compute the coefficient which improves the estimation.
  For example, it might be interesting to
  try an integration rule instead of the least squares method.
  Notice that a specific design of experiment is required in this case.

.. GENERATED FROM PYTHON SOURCE LINES 277-278

Reset default settings

.. GENERATED FROM PYTHON SOURCE LINES 278-279

.. code-block:: Python

    ot.ResourceMap.Reload()








.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 13.119 seconds)


.. _sphx_glr_download_auto_meta_modeling_polynomial_chaos_metamodel_plot_functional_chaos.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_functional_chaos.ipynb <plot_functional_chaos.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_functional_chaos.py <plot_functional_chaos.py>`