.. 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 6-17

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 19-21

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

.. GENERATED FROM PYTHON SOURCE LINES 23-30

.. 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 31-32

We load the use case :

.. GENERATED FROM PYTHON SOURCE LINES 32-34

.. code-block:: Python

    cb = cantilever_beam.CantileverBeam()








.. GENERATED FROM PYTHON SOURCE LINES 35-36

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

.. GENERATED FROM PYTHON SOURCE LINES 36-38

.. code-block:: Python

    model = cb.model








.. GENERATED FROM PYTHON SOURCE LINES 39-40

Then we define the distribution of the input random vector.

.. GENERATED FROM PYTHON SOURCE LINES 40-43

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 44-46

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

.. GENERATED FROM PYTHON SOURCE LINES 48-51

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 53-57

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 58-60

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

.. GENERATED FROM PYTHON SOURCE LINES 62-66

In order to create the Kriging metamodel, we first select a constant trend with the `ConstantBasisFactory` class.
Then we use a squared exponential covariance kernel.
The `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 68-71

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 72-74

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

.. GENERATED FROM PYTHON SOURCE LINES 76-79

.. 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 80-82

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 84-88

.. 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 89-93

Finally, we use the `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 95-99

.. 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 100-101

Run the algorithm and get the result.

.. GENERATED FROM PYTHON SOURCE LINES 103-107

.. code-block:: Python

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








.. GENERATED FROM PYTHON SOURCE LINES 108-109

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

.. GENERATED FROM PYTHON SOURCE LINES 111-113

.. code-block:: Python

    print(result.getTrendCoefficients())





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

 .. code-block:: none

    [0.108219]




.. GENERATED FROM PYTHON SOURCE LINES 114-115

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

.. GENERATED FROM PYTHON SOURCE LINES 117-119

.. 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 120-122

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

.. GENERATED FROM PYTHON SOURCE LINES 124-133

.. 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 134-147

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 149-160

.. 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 161-174

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 (dim+1)}{2}=\frac{4\times 5}{2}=10


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

.. GENERATED FROM PYTHON SOURCE LINES 176-178

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

.. GENERATED FROM PYTHON SOURCE LINES 180-181

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 183-188

.. code-block:: Python

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









.. GENERATED FROM PYTHON SOURCE LINES 189-199

.. code-block:: Python

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









.. GENERATED FROM PYTHON SOURCE LINES 200-214

.. 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 215-227

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>`