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

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

.. _sphx_glr_auto_meta_modeling_kriging_metamodel_plot_kriging_advanced.py:


Advanced kriging
================

.. GENERATED FROM PYTHON SOURCE LINES 6-10

In this example we will build a metamodel using gaussian process regression of the :math:`x\sin(x)` function.

We will choose the number of learning points, the basis and the covariance model.


.. GENERATED FROM PYTHON SOURCE LINES 13-21

.. code-block:: Python

    import openturns as ot
    from openturns.viewer import View
    import numpy as np
    import matplotlib.pyplot as plt
    import openturns.viewer as viewer

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








.. GENERATED FROM PYTHON SOURCE LINES 22-27

Generate design of experiment
-----------------------------

We create training samples from the function :math:`x\sin(x)`. We can change their number and distribution in the :math:`[0; 10]` range.
If the `with_error` boolean is `True`, then the data is computed by adding a gaussian noise to the function values.

.. GENERATED FROM PYTHON SOURCE LINES 29-35

.. code-block:: Python

    dim = 1
    xmin = 0
    xmax = 10
    n_pt = 20  # number of initial points
    with_error = True  # whether to use generation with error








.. GENERATED FROM PYTHON SOURCE LINES 36-49

.. code-block:: Python

    ref_func_with_error = ot.SymbolicFunction(["x", "eps"], ["x * sin(x) + eps"])
    ref_func = ot.ParametricFunction(ref_func_with_error, [1], [0.0])
    x = np.vstack(np.linspace(xmin, xmax, n_pt))
    ot.RandomGenerator.SetSeed(1235)
    eps = ot.Normal(0, 1.5).getSample(n_pt)
    X = ot.Sample(n_pt, 2)
    X[:, 0] = x
    X[:, 1] = eps
    if with_error:
        y = np.array(ref_func_with_error(X))
    else:
        y = np.array(ref_func(x))








.. GENERATED FROM PYTHON SOURCE LINES 50-60

.. code-block:: Python

    graph = ref_func.draw(xmin, xmax, 200)
    cloud = ot.Cloud(x, y)
    cloud.setColor("red")
    cloud.setPointStyle("bullet")
    graph.add(cloud)
    graph.setLegends(["Function", "Data"])
    graph.setLegendPosition("upper left")
    graph.setTitle("Sample size = %d" % (n_pt))
    view = viewer.View(graph)




.. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_001.png
   :alt: Sample size = 20
   :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 61-63

Create the kriging algorithm
----------------------------

.. GENERATED FROM PYTHON SOURCE LINES 63-93

.. code-block:: Python


    # 1. basis
    ot.ResourceMap.SetAsBool(
        "GeneralLinearModelAlgorithm-UseAnalyticalAmplitudeEstimate", True
    )
    basis = ot.ConstantBasisFactory(dim).build()
    print(basis)

    # 2. covariance model
    cov = ot.MaternModel([1.0], [2.5], 1.5)
    print(cov)

    # 3. kriging algorithm
    algokriging = ot.KrigingAlgorithm(x, y, cov, basis)

    # error measure
    # algokriging.setNoise([5*1e-1]*n_pt)

    # 4. Optimization
    # algokriging.setOptimizationAlgorithm(ot.NLopt('GN_DIRECT'))
    lhsExperiment = ot.LHSExperiment(ot.Uniform(1e-1, 1e2), 50)
    algokriging.setOptimizationAlgorithm(ot.MultiStart(ot.TNC(), lhsExperiment.generate()))
    algokriging.setOptimizationBounds(ot.Interval([0.1], [1e2]))

    # if we choose not to optimize parameters
    # algokriging.setOptimizeParameters(False)

    # 5. run the algorithm
    algokriging.run()





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

 .. code-block:: none

    Basis( [class=LinearEvaluation name=Unnamed center=[0] constant=[1] linear=[[ 0 ]]] )
    MaternModel(scale=[1], amplitude=[2.5], nu=1.5)




.. GENERATED FROM PYTHON SOURCE LINES 94-96

Results
-------

.. GENERATED FROM PYTHON SOURCE LINES 98-99

get some results

.. GENERATED FROM PYTHON SOURCE LINES 99-108

.. code-block:: Python

    krigingResult = algokriging.getResult()
    print("residual = ", krigingResult.getResiduals())
    print("R2 = ", krigingResult.getRelativeErrors())
    print("Optimal scale= {}".format(krigingResult.getCovarianceModel().getScale()))
    print(
        "Optimal amplitude = {}".format(krigingResult.getCovarianceModel().getAmplitude())
    )
    print("Optimal trend coefficients = {}".format(krigingResult.getTrendCoefficients()))





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

 .. code-block:: none

    residual =  [5.39875e-16]
    R2 =  [2.99965e-31]
    Optimal scale= [0.818671]
    Optimal amplitude = [4.51225]
    Optimal trend coefficients = [-0.115697]




.. GENERATED FROM PYTHON SOURCE LINES 109-110

get the metamodel

.. GENERATED FROM PYTHON SOURCE LINES 110-150

.. code-block:: Python

    krigingMeta = krigingResult.getMetaModel()

    n_pts_plot = 1000
    x_plot = np.vstack(np.linspace(xmin, xmax, n_pts_plot))
    fig, [ax1, ax2] = plt.subplots(1, 2, figsize=(12, 6))

    # On the left, the function
    graph = ref_func.draw(xmin, xmax, n_pts_plot)
    graph.setLegends(["Function"])
    graphKriging = krigingMeta.draw(xmin, xmax, n_pts_plot)
    graphKriging.setColors(["green"])
    graphKriging.setLegends(["Kriging"])
    graph.add(graphKriging)
    cloud = ot.Cloud(x, y)
    cloud.setColor("red")
    cloud.setLegend("Data")
    graph.add(cloud)
    graph.setLegendPosition("upper left")
    View(graph, axes=[ax1])

    # On the right, the conditional kriging variance
    graph = ot.Graph("", "x", "Conditional kriging variance", True, "")
    # Sample for the data
    sample = ot.Sample(n_pt, 2)
    sample[:, 0] = x
    cloud = ot.Cloud(sample)
    cloud.setColor("red")
    graph.add(cloud)
    # Sample for the variance
    sample = ot.Sample(n_pts_plot, 2)
    sample[:, 0] = x_plot
    variance = [[krigingResult.getConditionalCovariance(xx)[0, 0]] for xx in x_plot]
    sample[:, 1] = variance
    curve = ot.Curve(sample)
    curve.setColor("green")
    graph.add(curve)
    View(graph, axes=[ax2])

    fig.suptitle("Kriging result")




.. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_002.png
   :alt: Kriging result
   :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_002.png
   :class: sphx-glr-single-img


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

 .. code-block:: none


    Text(0.5, 0.98, 'Kriging result')



.. GENERATED FROM PYTHON SOURCE LINES 151-153

Display the confidence interval
-------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 155-156

sphinx_gallery_thumbnail_number = 3

.. GENERATED FROM PYTHON SOURCE LINES 156-174

.. code-block:: Python

    level = 0.95
    quantile = ot.Normal().computeQuantile((1 - level) / 2)[0]
    borne_sup = krigingMeta(x_plot) + quantile * np.sqrt(variance)
    borne_inf = krigingMeta(x_plot) - quantile * np.sqrt(variance)

    fig, ax = plt.subplots(figsize=(8, 8))
    ax.plot(x, y, ("ro"))
    ax.plot(x_plot, borne_sup, "--", color="orange", label="Confidence interval")
    ax.plot(x_plot, borne_inf, "--", color="orange")
    View(ref_func.draw(xmin, xmax, n_pts_plot), axes=[ax], plot_kw={"label": "$x sin(x)$"})
    View(
        krigingMeta.draw(xmin, xmax, n_pts_plot),
        plot_kw={"color": "green", "label": "prediction"},
        axes=[ax],
    )
    legend = ax.legend()
    ax.autoscale()




.. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_003.png
   :alt: v0 as a function of v0
   :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_003.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 175-177

Generate conditional trajectories
---------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 179-180

support for trajectories with training samples removed

.. GENERATED FROM PYTHON SOURCE LINES 180-185

.. code-block:: Python

    values = np.linspace(0, 10, 500)
    for xx in x:
        if len(np.argwhere(values == xx)) == 1:
            values = np.delete(values, np.argwhere(values == xx)[0, 0])








.. GENERATED FROM PYTHON SOURCE LINES 186-187

Conditional Gaussian process

.. GENERATED FROM PYTHON SOURCE LINES 187-190

.. code-block:: Python

    krv = ot.KrigingRandomVector(krigingResult, np.vstack(values))
    krv_sample = krv.getSample(5)








.. GENERATED FROM PYTHON SOURCE LINES 191-214

.. code-block:: Python

    x_plot = np.vstack(np.linspace(xmin, xmax, n_pts_plot))
    fig, ax = plt.subplots(figsize=(8, 6))
    ax.plot(x, y, ("ro"))
    for i in range(krv_sample.getSize()):
        if i == 0:
            ax.plot(
                values, krv_sample[i, :], "--", alpha=0.8, label="Conditional trajectories"
            )
        else:
            ax.plot(values, krv_sample[i, :], "--", alpha=0.8)
    View(
        ref_func.draw(xmin, xmax, n_pts_plot),
        axes=[ax],
        plot_kw={"color": "black", "label": "$x*sin(x)$"},
    )
    View(
        krigingMeta.draw(xmin, xmax, n_pts_plot),
        axes=[ax],
        plot_kw={"color": "green", "label": "prediction"},
    )
    legend = ax.legend()
    ax.autoscale()




.. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_004.png
   :alt: v0 as a function of v0
   :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_004.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 215-217

Validation
----------

.. GENERATED FROM PYTHON SOURCE LINES 219-228

.. code-block:: Python

    n_valid = 10
    x_valid = ot.Uniform(xmin, xmax).getSample(n_valid)
    if with_error:
        X_valid = ot.Sample(x_valid)
        X_valid.stack(ot.Normal(0.0, 1.5).getSample(n_valid))
        y_valid = np.array(ref_func_with_error(X_valid))
    else:
        y_valid = np.array(ref_func(X_valid))








.. GENERATED FROM PYTHON SOURCE LINES 229-232

.. code-block:: Python

    validation = ot.MetaModelValidation(x_valid, y_valid, krigingMeta)
    validation.computePredictivityFactor()






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    class=Point name=Unnamed dimension=1 values=[0.861246]
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 233-236

.. code-block:: Python

    graph = validation.drawValidation()
    view = viewer.View(graph)




.. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_005.png
   :alt: Metamodel validation - n = 10
   :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_005.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 237-242

.. code-block:: Python

    graph = validation.getResidualDistribution().drawPDF()
    graph.setXTitle("Residuals")
    view = viewer.View(graph)
    plt.show()




.. image-sg:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_006.png
   :alt: plot kriging advanced
   :srcset: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_advanced_006.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 243-244

Reset default settings

.. GENERATED FROM PYTHON SOURCE LINES 244-245

.. code-block:: Python

    ot.ResourceMap.Reload()








.. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_plot_kriging_advanced.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_advanced.ipynb <plot_kriging_advanced.ipynb>`

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

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