.. 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_fields_metamodels_plot_viscous_fall_metamodel.py:
Viscous free fall: metamodel of a field function
================================================
In this example, we present how to create the metamodel of a field function. This examples considers the :ref:`free fall model `. We first compute the Karhunen-Loève decomposition of a sample of trajectories. Then we create a create a polynomial chaos which takes the inputs and returns the KL decomposition modes as outputs. Finally, we create a metamodel by combining the KL decomposition and the polynomial chaos.
Define the model
----------------
.. code-block:: default
from __future__ import print_function
import openturns as ot
import numpy as np
import openturns.viewer as viewer
from matplotlib import pylab as plt
ot.Log.Show(ot.Log.NONE)
We first define the time grid associated with the model.
.. code-block:: default
tmin=0.0 # Minimum time
tmax=12. # Maximum time
gridsize=100 # Number of time steps
mesh = ot.IntervalMesher([gridsize-1]).build(ot.Interval(tmin, tmax))
.. code-block:: default
vertices = mesh.getVertices()
Creation of the input distribution.
.. code-block:: default
distZ0 = ot.Uniform(100.0, 150.0)
distV0 = ot.Normal(55.0, 10.0)
distM = ot.Normal(80.0, 8.0)
distC = ot.Uniform(0.0, 30.0)
distribution = ot.ComposedDistribution([distZ0, distV0, distM, distC])
.. code-block:: default
dimension = distribution.getDimension()
dimension
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
4
Then we define the Python function which computes the altitude at each time value. In order to compute all altitudes with a vectorized evaluation, we first convert the vertices into a Numpy `array` and use the Numpy functions`exp` and `maximum`: this increases the evaluation performance of the script.
.. code-block:: default
def AltiFunc(X):
g = 9.81
z0 = X[0]
v0 = X[1]
m = X[2]
c = X[3]
tau = m / c
vinf = - m * g / c
t = np.array(vertices)
z = z0 + vinf * t + tau * (v0 - vinf) * (1 - np.exp( - t / tau))
z = np.maximum(z,0.)
return [[zeta[0]] for zeta in z]
In order to create a `Function` from this Python function, we use the `PythonPointToFieldFunction` class. Since the altitude is the only output field, the third argument `outputDimension` is equal to `1`. If we had computed the speed as an extra output field, we would have set `2` instead.
.. code-block:: default
outputDimension = 1
alti = ot.PythonPointToFieldFunction(dimension, mesh, outputDimension, AltiFunc)
Compute a training sample.
.. code-block:: default
size = 2000
ot.RandomGenerator.SetSeed(0)
inputSample = distribution.getSample(size)
outputSample = alti(inputSample)
Compute the KL decomposition of the output
------------------------------------------
.. code-block:: default
algo = ot.KarhunenLoeveSVDAlgorithm(outputSample, 1.0e-6)
algo.run()
KLResult = algo.getResult()
scaledModes = KLResult.getScaledModesAsProcessSample()
.. code-block:: default
graph = scaledModes.drawMarginal(0)
graph.setTitle('KL modes')
graph.setXTitle(r'$t$')
graph.setYTitle(r'$z$')
view = viewer.View(graph)
.. image:: /auto_meta_modeling/fields_metamodels/images/sphx_glr_plot_viscous_fall_metamodel_001.png
:alt: KL modes
:class: sphx-glr-single-img
We create the `postProcessingKL` function which takes coefficients of the the K.-L. modes as inputs and returns the trajectories.
.. code-block:: default
karhunenLoeveLiftingFunction = ot.KarhunenLoeveLifting(KLResult)
The `project` method computes the projection of the output sample (i.e. the trajectories) onto the K.-L. modes.
.. code-block:: default
outputSampleChaos = KLResult.project(outputSample)
We limit the sampling size of the Lilliefors selection in order to reduce the computational burden.
.. code-block:: default
ot.ResourceMap.SetAsUnsignedInteger("FittingTest-LillieforsMaximumSamplingSize", 1)
We create a polynomial chaos metamodel which takes the input sample and returns the K.-L. modes.
.. code-block:: default
algo = ot.FunctionalChaosAlgorithm(inputSample, outputSampleChaos)
algo.run()
chaosMetamodel = algo.getResult().getMetaModel()
The final metamodel is a composition of the KL lifting function and the polynomial chaos metamodel. In order to combine these two functions, we use the `PointToFieldConnection` class.
.. code-block:: default
metaModel = ot.PointToFieldConnection(karhunenLoeveLiftingFunction, chaosMetamodel)
Validate the metamodel
----------------------
Create a validation sample.
.. code-block:: default
size = 10
validationInputSample = distribution.getSample(size)
validationOutputSample = alti(validationInputSample)
.. code-block:: default
graph = validationOutputSample.drawMarginal(0)
graph.setColors(['red'])
graph2 = metaModel(validationInputSample).drawMarginal(0)
graph2.setColors(['blue'])
graph.add(graph2)
graph.setTitle('Model/metamodel comparison')
graph.setXTitle(r'$t$')
graph.setYTitle(r'$z$')
view = viewer.View(graph)
plt.show()
.. image:: /auto_meta_modeling/fields_metamodels/images/sphx_glr_plot_viscous_fall_metamodel_002.png
:alt: Model/metamodel comparison
:class: sphx-glr-single-img
We see that the blue trajectories (i.e. the metamodel) are close to the red trajectories (i.e. the validation sample). This shows that the metamodel is quite accurate. However, we observe that the trajectory singularity that occurs when the object touches the ground (i.e. when :math:`z` is equal to zero), makes the metamodel less accurate.
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 0 minutes 27.279 seconds)
.. _sphx_glr_download_auto_meta_modeling_fields_metamodels_plot_viscous_fall_metamodel.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_viscous_fall_metamodel.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_viscous_fall_metamodel.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_