1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
|
"""
Manipulate stochastic processes
===============================
"""
# %%
# The objective here is to manipulate a multivariate stochastic process :math:`X: \Omega \times \mathcal{D} \rightarrow \mathbb{R}^d`,
# where :math:`\mathcal{D} \in \mathbb{R}^n` is discretized on the mesh :math:`\mathcal{M}`
# and exhibit some of the services exposed by the :class:`~openturns.Process` objects:
#
# - ask for the dimension, with the method :meth:`~openturns.Process.getOutputDimension` ;
# - ask for the mesh, with the method :meth:`~openturns.Process.getMesh` ;
# - ask for the mesh as regular 1-d mesh, with the :meth:`~openturns.Process.getTimeGrid` method ;
# - ask for a realization, with the method the :meth:`~openturns.Process.getRealization` method ;
# - ask for a continuous realization, with the :meth:`~openturns.Process.getContinuousRealization` method ;
# - ask for a sample of realizations, with the :meth:`~openturns.Process.getSample` method ;
# - ask for the normality of the process with the :meth:`~openturns.Process.isNormal` method ;
# - ask for the stationarity of the process with the :meth:`~openturns.Process.isStationary` method.
#
# %%
import openturns as ot
import openturns.viewer as otv
# %%
# We create a mesh -a time grid- which is a :class:`~openturns.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 :math:`X_t = (X_t^0, X_t^1, X_t^2)` in dimension 3 with an exponential covariance model.
# We define the amplitude and the scale of the :class:`~openturns.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)
# %%
# 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 = otv.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())
# %%
# A typical feature of a process is the generation of a realization of itself :
realization = process.getRealization()
# %%
# Here it is a sample of size :math:`100 \times 4` (100 time steps, 3 spatial cooordinates and the time variable).
# We are able to draw its marginals, for instance the first (index 0) one :math:`X_t^0`, 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 = otv.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 = otv.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 = otv.View(graph)
# %%
# Please note that the `marginal0` object is a function. Consequently we can
# evaluate it at any point of the domain, say :math:`t_0=3.5678`
t0 = 3.5678
print(t0, marginal0([t0]))
# %%
# Several realizations of the process may be determined at once :
number = 10
fieldSample = process.getSample(number)
# %%
# Let us draw them the first marginal
# sphinx_gallery_thumbnail_number = 5
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 = otv.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 = otv.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 1-d Gaussian process
print(process.getMarginal([1]))
# %%
# If we extract simultaneously two indices we build a 2-d Gaussian process
print(process.getMarginal([0, 2]))
# %%
# We can check whether the process is Gaussian or not
print("Is normal ? ", process.isNormal())
# %%
# and the stationarity as well
print("Is stationary ? ", process.isStationary())
# %%
# Display all figures
otv.View.ShowAll()
|
