.. only:: html .. note:: :class: sphx-glr-download-link-note Click :ref:here  to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_meta_modeling_kriging_metamodel_plot_kriging_simulate.py: Simulate new trajectories from a kriging metamodel ================================================== The main goal of this example is to show how to simulate new trajectories from a kriging metamodel. Introduction ------------ We consider the sine function: .. math:: y = \sin(x) for any :math:x\in[0,12]. We want to create a metamodel of this function. This is why we create a sample of :math:n observations of the function: .. math:: y_i=\sin(x_i) for :math:i=1,...,7, where :math:x_i is the i-th input and :math:y_i is the corresponding output. We consider the seven following inputs : ============ === === === === ===== ==== ====== :math:i 1 2 3 4 5 6 7 :math:x_i 1 3 4 6 7.9 11 11.5 ============ === === === === ===== ==== ====== We are going to consider a kriging metamodel with a * constant trend, * a Matern covariance model. Creation of the metamodel ------------------------- We begin by defining the function g as a symbolic function. Then we define the x_train variable which contains the inputs of the design of experiments of the training step. Then we compute the y_train corresponding outputs. The variable n_train is the size of the training design of experiments. .. code-block:: default import numpy as np import openturns as ot import openturns.viewer as viewer from matplotlib import pylab as plt ot.Log.Show(ot.Log.NONE) .. code-block:: default g = ot.SymbolicFunction(['x'], ['sin(x)']) .. code-block:: default x_train = ot.Sample([[x] for x in [1.,3.,4.,6.,7.9,11., 11.5]]) y_train = g(x_train) n_train = x_train.getSize() n_train .. rst-class:: sphx-glr-script-out Out: .. code-block:: none 7 In order to compare the function and its metamodel, we use a test (i.e. validation) design of experiments made of a regular grid of 100 points from 0 to 12. Then we convert this grid into a Sample and we compute the outputs of the function on this sample. .. code-block:: default xmin = 0. xmax = 12. n_test = 100 step = (xmax-xmin)/(n_test-1) myRegularGrid = ot.RegularGrid(xmin, step, n_test) x_test_coord = myRegularGrid.getValues() x_test = ot.Sample([[x] for x in x_test_coord]) y_test = g(x_test) In order to observe the function and the location of the points in the input design of experiments, we define the following functions which plots the data. .. code-block:: default def plot_data_train(x_train,y_train): '''Plot the data (x_train,y_train) as a Cloud, in red''' graph_train = ot.Cloud(x_train,y_train) graph_train.setColor("red") graph_train.setLegend("Data") return graph_train .. code-block:: default def plot_data_test(x_test,y_test): '''Plot the data (x_test,y_test) as a Curve, in dashed black''' graphF = ot.Curve(x_test,y_test) graphF.setLegend("Exact") graphF.setColor("black") graphF.setLineStyle("dashed") return graphF .. code-block:: default graph = ot.Graph() graph.add(plot_data_test(x_test,y_test)) graph.add(plot_data_train(x_train,y_train)) graph.setAxes(True) graph.setXTitle("X") graph.setYTitle("Y") graph.setLegendPosition("topright") view = viewer.View(graph) .. image:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_simulate_001.png :alt: plot kriging simulate :class: sphx-glr-single-img We use the ConstantBasisFactory class to define the trend and the MaternModel class to define the covariance model. This Matérn model is based on the regularity parameter :math:\nu=3/2. .. code-block:: default dimension = 1 basis = ot.ConstantBasisFactory(dimension).build() covarianceModel = ot.MaternModel([1.]*dimension, 1.5) algo = ot.KrigingAlgorithm(x_train, y_train, covarianceModel, basis) algo.run() krigingResult = algo.getResult() krigingResult .. raw:: html

KrigingResult(covariance models=MaternModel(scale=[1.26541], amplitude=[0.818994], nu=1.5), covariance coefficients=0 : [ 1.14764 ]
1 : [ 1.006 ]
2 : [ -1.75286 ]
3 : [ -0.55873 ]
4 : [ 1.78675 ]
5 : [ -1.59663 ]
6 : [ -0.032173 ], basis=[Basis( [class=LinearEvaluation name=Unnamed center= constant= linear=[[ 0 ]]] )], trend coefficients=[[0.00671564]])

We observe that the scale and amplitude hyper-parameters have been optimized by the run method. Then we get the metamodel with getMetaModel and evaluate the outputs of the metamodel on the test design of experiments. .. code-block:: default krigeageMM = krigingResult.getMetaModel() y_test_MM = krigeageMM(x_test) The following function plots the kriging data. .. code-block:: default def plot_data_kriging(x_test,y_test_MM): '''Plots (x_test,y_test_MM) from the metamodel as a Curve, in blue''' graphK = ot.Curve(x_test,y_test_MM) graphK.setColor("blue") graphK.setLegend("Kriging") return graphK .. code-block:: default graph = ot.Graph() graph.add(plot_data_test(x_test,y_test)) graph.add(plot_data_train(x_train,y_train)) graph.add(plot_data_kriging(x_test,y_test_MM)) graph.setAxes(True) graph.setXTitle("X") graph.setYTitle("Y") graph.setLegendPosition("topright") view = viewer.View(graph) .. image:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_simulate_002.png :alt: plot kriging simulate :class: sphx-glr-single-img Simulate new trajectories ------------------------- In order to generate new trajectories of the conditioned gaussian process, we couild technically use the KrigingRandomVector class, because it provides the getSample method that we need. However, the KrigingRandomVector class was more specifically designed to create a RandomVector so that it can feed, for example, a function which has a field as input argument. This is why we use the ConditionedGaussianProcess, which provides a Process. .. code-block:: default n_test = 100 step = (xmax-xmin)/(n_test-1) myRegularGrid = ot.RegularGrid(xmin, step, n_test) vertices = myRegularGrid.getVertices() If we directly use the vertices values, we get: RuntimeError: InternalException : Error: the matrix is not definite positive. Indeed, some points in vertices are also in x_train. This is why the conditioned covariance matrix is singular at these points. This is why we define the following function which deletes points in vertices which are also found in x_train. .. code-block:: default def deleteCommonValues(x_train,x_test): ''' Delete from x_test the values which are in x_train so that values in x_test have no interect with x_train. ''' x_test_filtered = x_test # Initialize for x_train_value in x_train: print("Checking %s" % (x_train_value)) indices = np.argwhere(x_test==x_train_value) if len(indices) == 1: print(" Delete %s" % (x_train_value)) x_test_filtered = np.delete(x_test_filtered, indices[0, 0]) else: print(" OK") return x_test_filtered .. code-block:: default vertices_filtered = deleteCommonValues(np.array(x_train.asPoint()),np.array(vertices.asPoint())) .. rst-class:: sphx-glr-script-out Out: .. code-block:: none Checking 1.0 OK Checking 3.0 OK Checking 4.0 Delete 4.0 Checking 6.0 OK Checking 7.9 OK Checking 11.0 OK Checking 11.5 OK .. code-block:: default evaluationMesh = ot.Mesh(ot.Sample([[vf] for vf in vertices_filtered])) .. code-block:: default process = ot.ConditionedGaussianProcess(krigingResult, evaluationMesh) .. code-block:: default trajectories = process.getSample(10) type(trajectories) The getSample method returns a ProcessSample. By comparison, the getSample method of a KrigingRandomVector would return a Sample. .. code-block:: default graph = trajectories.drawMarginal() graph.add(plot_data_test(x_test,y_test)) graph.add(plot_data_train(x_train,y_train)) graph.setAxes(True) graph.setXTitle("X") graph.setYTitle("Y") graph.setLegendPosition("topright") graph.setTitle("10 simulated trajectories") view = viewer.View(graph) plt.show() .. image:: /auto_meta_modeling/kriging_metamodel/images/sphx_glr_plot_kriging_simulate_003.png :alt: 10 simulated trajectories :class: sphx-glr-single-img References ---------- * Metamodeling with Gaussian processes, Bertrand Iooss, EDF R&D, 2014, www.gdr-mascotnum.fr/media/sssamo14_iooss.pdf .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 0.265 seconds) .. _sphx_glr_download_auto_meta_modeling_kriging_metamodel_plot_kriging_simulate.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:Download Python source code: plot_kriging_simulate.py  .. container:: sphx-glr-download sphx-glr-download-jupyter :download:Download Jupyter notebook: plot_kriging_simulate.ipynb  .. only:: html .. rst-class:: sphx-glr-signature Gallery generated by Sphinx-Gallery _