Note
Click here to download the full example code
Draw fields¶
The objective here is to manipulate a multivariate stochastic process , where is discretized on the mesh and exhibit some of the services exposed by the Process objects:
ask for the dimension, with the method getOutputDimension
ask for the mesh, with the method getMesh
ask for the mesh as regular 1-d mesh, with the getTimeGrid method
ask for a realization, with the method the getRealization method
ask for a continuous realization, with the getContinuousRealization method
ask for a sample of realizations, with the getSample method
ask for the normality of the process with the isNormal method
ask for the stationarity of the process with the isStationary method
from __future__ import print_function
import openturns as ot
import openturns.viewer as viewer
from matplotlib import pylab as plt
import math as m
ot.Log.Show(ot.Log.NONE)
We create a mesh -a time grid- which is a RegularGrid :
tMin = 0.0
timeStep = 0.1
n = 100
time_grid = ot.RegularGrid(tMin, timeStep, n)
time_grid.setName('time')
We create a Normal process in dimension 3 with an exponential covariance model. We define the amplitude and the scale of the ExponentialModel
scale = [4.0]
amplitude = [1.0, 2.0, 3.0]
We define a spatial correlation :
spatialCorrelation = ot.CorrelationMatrix(3)
spatialCorrelation[0, 1] = 0.8
spatialCorrelation[0, 2] = 0.6
spatialCorrelation[1, 2] = 0.1
The covariance model is now created with :
myCovarianceModel = ot.ExponentialModel(scale, amplitude, spatialCorrelation)
Eventually, the process is built with :
process = ot.GaussianProcess(myCovarianceModel, time_grid)
The dimension d of the process may be retrieved by
dim = process.getOutputDimension()
print("Dimension : %d" % dim)
Out:
Dimension : 3
The underlying mesh of the process is obtained with the getMesh method :
mesh = process.getMesh()
We have access to peculiar data of the mesh such as the corners :
minMesh = mesh.getVertices().getMin()[0]
maxMesh = mesh.getVertices().getMax()[0]
We draw it :
graph = mesh.draw()
graph.setTitle("Time grid (mesh)")
graph.setXTitle("t")
graph.setYTitle("")
view = viewer.View(graph)
We can get the time grid of the process when the mesh can be interpreted as a regular time grid :
print(process.getTimeGrid())
Out:
RegularGrid(start=0, step=0.1, n=100)
A typical feature of a process is the generation of a realization of itself :
realization = process.getRealization()
Here it is a sample of size (100 time steps, 3 spatial cooordinates and the time variable). We are able to draw its marginals, for instance the first (index 0) one , against the time with no interpolation :
interpolate = False
graph = realization.drawMarginal(0, interpolate)
graph.setTitle("First marginal of a realization of the process")
graph.setXTitle("t")
graph.setYTitle(r"$X_t^0$")
view = viewer.View(graph)
The same graph, but with interpolated values (default behaviour) :
graph = realization.drawMarginal(0)
graph.setTitle("First marginal of a realization of the process")
graph.setXTitle("t")
graph.setYTitle(r"$X_t^0$")
view = viewer.View(graph)
We can build a function representing the process using P1-Lagrange interpolation (when not defined from a functional model).
continuousRealization = process.getContinuousRealization()
Once again we draw its first marginal :
marginal0 = continuousRealization.getMarginal(0)
graph = marginal0.draw(minMesh, maxMesh)
graph.setTitle(
"First marginal of a P1-Lagrange continuous realization of the process")
graph.setXTitle("t")
view = viewer.View(graph)
Please note that the marginal0 object is a function. Consequently we can evaluate it at any point of the domain, say
t0 = 3.5678
print(t0, marginal0([t0]))
Out:
3.5678 [-0.593188]
Several realizations of the process may be determined at once :
number = 10
fieldSample = process.getSample(number)
Let us draw them the first marginal)
graph = fieldSample.drawMarginal(0, False)
graph.setTitle("First marginal of 10 realizations of the process")
graph.setXTitle("t")
graph.setYTitle(r"$X_t^0$")
view = viewer.View(graph)
Same graph, but with interpolated values :
graph = fieldSample.drawMarginal(0)
graph.setTitle("First marginal of 10 realizations of the process")
graph.setXTitle("t")
graph.setYTitle(r"$X_t^0$")
view = viewer.View(graph)
Miscellanies¶
We can extract any marginal of the process with the getMarginal method. Beware the numerotation begins at 0 ! It may be not implemented yet for some processes. The extracted marginal is a 1D gaussian process :
print(process.getMarginal([1]))
Out:
GaussianProcess(trend=[x0]->[0.0], covariance=ExponentialModel(scale=[4], amplitude=[2], no spatial correlation))
If we extract simultaneously two indices we build a 2D gaussian process :
print(process.getMarginal([0, 2]))
Out:
GaussianProcess(trend=[x0]->[0.0,0.0], covariance=ExponentialModel(scale=[4], amplitude=[1,3], spatial correlation=
[[ 1 0.6 ]
[ 0.6 1 ]]))
We can check whether the process is normal or not :
print("Is normal ? ", process.isNormal())
Out:
Is normal ? True
and the stationarity as well :
print("Is stationary ? ", process.isStationary())
plt.show()
Out:
Is stationary ? True
Total running time of the script: ( 0 minutes 0.697 seconds)