.. 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_kriging_beam_trend.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_kriging_metamodel_plot_kriging_beam_trend.py>`
        to download the full example code.

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

.. _sphx_glr_auto_meta_modeling_kriging_metamodel_plot_kriging_beam_trend.py:


Kriging: choose a polynomial trend on the beam model
====================================================

.. GENERATED FROM PYTHON SOURCE LINES 7-18

The goal of this example is to show how to configure the trend in a Kriging metamodel.
This example focuses on three polynomial trends:

* :class:`~openturns.ConstantBasisFactory`,
* :class:`~openturns.LinearBasisFactory`,
* :class:`~openturns.QuadraticBasisFactory`.

In the :doc:`/auto_meta_modeling/kriging_metamodel/plot_kriging_chose_trend` example,
we give another example of this procedure.

For this purpose, we use the :ref:`cantilever beam <use-case-cantilever-beam>` example.

.. GENERATED FROM PYTHON SOURCE LINES 20-22

Definition of the model
-----------------------

.. GENERATED FROM PYTHON SOURCE LINES 24-31

.. code-block:: Python

    from openturns.usecases import cantilever_beam
    import openturns as ot
    from openturns.viewer import View

    ot.RandomGenerator.SetSeed(0)
    ot.Log.Show(ot.Log.NONE)








.. GENERATED FROM PYTHON SOURCE LINES 32-33

We load the use case :

.. GENERATED FROM PYTHON SOURCE LINES 33-35

.. code-block:: Python

    cb = cantilever_beam.CantileverBeam()








.. GENERATED FROM PYTHON SOURCE LINES 36-37

We define the function which evaluates the output depending on the inputs.

.. GENERATED FROM PYTHON SOURCE LINES 37-39

.. code-block:: Python

    model = cb.model








.. GENERATED FROM PYTHON SOURCE LINES 40-41

Then we define the distribution of the input random vector.

.. GENERATED FROM PYTHON SOURCE LINES 41-44

.. code-block:: Python

    dimension = cb.dim  # number of inputs
    myDistribution = cb.distribution








.. GENERATED FROM PYTHON SOURCE LINES 45-47

Create the design of experiments
--------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 49-52

We consider a simple Monte-Carlo sampling as a design of experiments.
This is why we generate an input sample using the :meth:`~openturns.Distribution.getSample` method of the distribution.
Then we evaluate the output using the `model` function.

.. GENERATED FROM PYTHON SOURCE LINES 54-58

.. code-block:: Python

    sampleSize_train = 10
    X_train = myDistribution.getSample(sampleSize_train)
    Y_train = model(X_train)








.. GENERATED FROM PYTHON SOURCE LINES 59-61

Create the metamodel
--------------------

.. GENERATED FROM PYTHON SOURCE LINES 63-67

In order to create the Kriging metamodel, we first select a constant trend with the :class:`~openturns.ConstantBasisFactory` class.
Then we use a squared exponential covariance kernel.
The :class:`~openturns.SquaredExponential` kernel has one amplitude coefficient and 4 scale coefficients.
This is because this covariance kernel is anisotropic : each of the 4 input variables is associated with its own scale coefficient.

.. GENERATED FROM PYTHON SOURCE LINES 69-72

.. code-block:: Python

    basis = ot.ConstantBasisFactory(dimension).build()
    covarianceModel = ot.SquaredExponential(dimension)








.. GENERATED FROM PYTHON SOURCE LINES 73-75

Typically, the optimization algorithm is quite good at setting optimization bounds.
In this case, however, the range of the input domain is extreme.

.. GENERATED FROM PYTHON SOURCE LINES 77-80

.. code-block:: Python

    print("Lower and upper bounds of X_train:")
    print(X_train.getMin(), X_train.getMax())





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

 .. code-block:: none

    Lower and upper bounds of X_train:
    [6.50185e+10,262.654,2.50948,1.40294e-07] [6.88439e+10,323.088,2.59143,1.5807e-07]




.. GENERATED FROM PYTHON SOURCE LINES 81-83

We need to manually define sensible optimization bounds.
Note that since the amplitude parameter is computed analytically (this is possible when the output dimension is 1), we only need to set bounds on the scale parameter.

.. GENERATED FROM PYTHON SOURCE LINES 85-89

.. code-block:: Python

    scaleOptimizationBounds = ot.Interval(
        [1.0, 1.0, 1.0, 1.0e-10], [1.0e11, 1.0e3, 1.0e1, 1.0e-5]
    )








.. GENERATED FROM PYTHON SOURCE LINES 90-94

Finally, we use the :class:`~openturns.KrigingAlgorithm` class to create the Kriging metamodel.
It requires a training sample, a covariance kernel and a trend basis as input arguments.
We need to set the initial scale parameter for the optimization. The upper bound of the input domain is a sensible choice here.
We must not forget to actually set the optimization bounds defined above.

.. GENERATED FROM PYTHON SOURCE LINES 96-100

.. code-block:: Python

    covarianceModel.setScale(X_train.getMax())
    algo = ot.KrigingAlgorithm(X_train, Y_train, covarianceModel, basis)
    algo.setOptimizationBounds(scaleOptimizationBounds)








.. GENERATED FROM PYTHON SOURCE LINES 101-102

Run the algorithm and get the result.

.. GENERATED FROM PYTHON SOURCE LINES 104-108

.. code-block:: Python

    algo.run()
    result = algo.getResult()
    krigingWithConstantTrend = result.getMetaModel()








.. GENERATED FROM PYTHON SOURCE LINES 109-110

The `getTrendCoefficients` method returns the coefficients of the trend.

.. GENERATED FROM PYTHON SOURCE LINES 112-114

.. code-block:: Python

    print(result.getTrendCoefficients())





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

 .. code-block:: none

    [0.108219]




.. GENERATED FROM PYTHON SOURCE LINES 115-116

The constant trend always has only one coefficient (if there is one single output).

.. GENERATED FROM PYTHON SOURCE LINES 118-120

.. code-block:: Python

    print(result.getCovarianceModel())





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

 .. code-block:: none

    SquaredExponential(scale=[6.88439e+10,323.088,2.59143,1.5807e-07], amplitude=[0.247034])




.. GENERATED FROM PYTHON SOURCE LINES 121-123

Setting the trend
-----------------

.. GENERATED FROM PYTHON SOURCE LINES 125-134

.. code-block:: Python

    covarianceModel.setScale(X_train.getMax())
    basis = ot.LinearBasisFactory(dimension).build()
    algo = ot.KrigingAlgorithm(X_train, Y_train, covarianceModel, basis)
    algo.setOptimizationBounds(scaleOptimizationBounds)
    algo.run()
    result = algo.getResult()
    krigingWithLinearTrend = result.getMetaModel()
    result.getTrendCoefficients()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Unnamed dimension=5 values=[-8.81529e-25,-1.00576e-13,-2.65632e-22,-2.22154e-24,-1.3279e-31]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 135-148

The number of coefficients in the linear and quadratic trends depends on the number of inputs, which is
equal to

.. math::
   dim = 4


in the cantilever beam case.

We observe that the number of coefficients in the trend is 5, which corresponds to:

* 1 coefficient for the constant part,
* dim = 4 coefficients for the linear part.

.. GENERATED FROM PYTHON SOURCE LINES 150-161

.. code-block:: Python

    covarianceModel.setScale(X_train.getMax())
    basis = ot.QuadraticBasisFactory(dimension).build()
    algo = ot.KrigingAlgorithm(X_train, Y_train, covarianceModel, basis)
    algo.setOptimizationBounds(scaleOptimizationBounds)
    algo.run()
    result = algo.getResult()
    krigingWithQuadraticTrend = result.getMetaModel()
    result.getTrendCoefficients()
    print(algo.getOptimizationBounds())
    print(result.getCovarianceModel())





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

 .. code-block:: none

    [1, 1e+11]
    [1, 1000]
    [1, 10]
    [1e-10, 1e-05]
    SquaredExponential(scale=[6.88439e+10,323.088,2.59143,1.5807e-07], amplitude=[1.74417e-10])




.. GENERATED FROM PYTHON SOURCE LINES 162-175

This time, the trend has 15 coefficients:

* 1 coefficient for the constant part,
* 4 coefficients for the linear part,
* 10 coefficients for the quadratic part.

This is because the number of coefficients in the quadratic part has

.. math::
   \frac{dim \times (dim+1)}{2}=\frac{4\times 5}{2}=10


coefficients, associated with the symmetric matrix of the quadratic function.

.. GENERATED FROM PYTHON SOURCE LINES 177-179

Validate the metamodel
----------------------

.. GENERATED FROM PYTHON SOURCE LINES 181-182

We finally want to validate the Kriging metamodel. This is why we generate a validation sample with size 100 and we evaluate the output of the model on this sample.

.. GENERATED FROM PYTHON SOURCE LINES 184-189

.. code-block:: Python

    sampleSize_test = 100
    X_test = myDistribution.getSample(sampleSize_test)
    Y_test = model(X_test)









.. GENERATED FROM PYTHON SOURCE LINES 190-201

.. code-block:: Python

    def drawMetaModelValidation(X_test, Y_test, krigingMetamodel, title):
        metamodelPredictions = krigingMetamodel(X_test)
        val = ot.MetaModelValidation(Y_test, metamodelPredictions)
        r2Score = val.computeR2Score()[0]
        graph = val.drawValidation().getGraph(0, 0)
        graph.setLegends([""])
        graph.setLegends(["%s, R2 = %.2f%%" % (title, 100 * r2Score), ""])
        graph.setLegendPosition("upper left")
        return graph









.. GENERATED FROM PYTHON SOURCE LINES 202-216

.. code-block:: Python

    grid = ot.GridLayout(1, 3)
    grid.setTitle("Different trends")
    graphConstant = drawMetaModelValidation(
        X_test, Y_test, krigingWithConstantTrend, "Constant"
    )
    graphLinear = drawMetaModelValidation(X_test, Y_test, krigingWithLinearTrend, "Linear")
    graphQuadratic = drawMetaModelValidation(
        X_test, Y_test, krigingWithQuadraticTrend, "Quadratic"
    )
    grid.setGraph(0, 0, graphConstant)
    grid.setGraph(0, 1, graphLinear)
    grid.setGraph(0, 2, graphQuadratic)
    _ = View(grid, figure_kw={"figsize": (13, 4)})




.. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_beam_trend_001.png
   :alt: Different trends
   :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_beam_trend_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 217-229

We observe that the three trends perform very well in this case.
With more coefficients, the Kriging metamodel is more flexibile and can adjust better to the training sample.
This does not mean, however, that the trend coefficients will provide a good fit for the validation sample.

The number of parameters in each Kriging metamodel is the following :

* the constant trend Kriging has 6 coefficients to estimate: 5 coefficients for the covariance matrix and 1 coefficient for the trend,
* the linear trend Kriging has 10 coefficients to estimate: 5 coefficients for the covariance matrix and 5 coefficients for the trend,
* the quadratic trend Kriging has 20 coefficients to estimate: 5 coefficients for the covariance matrix and 15 coefficients for the trend.

In the cantilever beam example, fitting the metamodel to a training sample with only 10 points is made much easier because the function to approximate is almost linear.
In this case, a quadratic trend is not advisable because it can interpolate all points in the training sample.


.. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_plot_kriging_beam_trend.py:

.. only:: html

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

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

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

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

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

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_kriging_beam_trend.zip <plot_kriging_beam_trend.zip>`