1#! /usr/bin/env python
2
3from __future__ import print_function
4import openturns as ot
5import openturns.testing as ott
6from math import sqrt
7
8ot.TESTPREAMBLE()
9
10
11def test_model(myModel, test_partial_grad=True, x1=None, x2=None):
12
13    inputDimension = myModel.getInputDimension()
14    dimension = myModel.getOutputDimension()
15
16    if x1 is None and x2 is None:
17        x1 = ot.Point(inputDimension)
18        x2 = ot.Point(inputDimension)
19        for j in range(inputDimension):
20            x1[j] = -1.0 - j
21            x2[j] = 3.0 + 2.0 * j
22    else:
23        x1 = ot.Point(x1)
24        x2 = ot.Point(x2)
25
26    if myModel.isStationary():
27        ott.assert_almost_equal(
28            myModel(x1 - x2), myModel(x1, x2), 1e-14, 1e-14)
29        ott.assert_almost_equal(
30            myModel(x2 - x1), myModel(x1, x2), 1e-14, 1e-14)
31
32    eps = 1e-3
33
34    mesh = ot.IntervalMesher([7]*inputDimension).build(
35        ot.Interval([-10]*inputDimension, [10]*inputDimension))
36
37    C = myModel.discretize(mesh)
38    if dimension == 1:
39        # Check that discretize & computeAsScalar provide the
40        # same values
41        vertices = mesh.getVertices()
42        for j in range(len(vertices)):
43            for i in range(j, len(vertices)):
44                ott.assert_almost_equal(C[i, j], myModel.computeAsScalar(
45                    vertices[i], vertices[j]), 1e-14, 1e-14)
46    else:
47        # Check that discretize & operator() provide the same values
48        vertices = mesh.getVertices()
49        localMatrix = ot.SquareMatrix(dimension)
50        for j in range(len(vertices)):
51            for i in range(j, len(vertices)):
52                for localJ in range(dimension):
53                    for localI in range(dimension):
54                        localMatrix[localI, localJ] = C[i *
55                                                        dimension + localI, j * dimension + localJ]
56                ott.assert_almost_equal(localMatrix, myModel(
57                    vertices[i], vertices[j]), 1e-14, 1e-14)
58
59    # Now we suppose that discretize is ok
60    # we look at crossCovariance of (vertices, vertices) which should return the same values
61    C.getImplementation().symmetrize()
62    crossCov = myModel.computeCrossCovariance(vertices, vertices)
63    ott.assert_almost_equal(crossCov, C, 1e-14, 1e-14,  "in " +
64                            myModel.getClassName() + "::computeCrossCovariance")
65
66    # Now crossCovariance(sample, sample) is ok
67    # Let us validate crossCovariance(Sample, point) with 1st column(s) of previous calculations
68    crossCovSamplePoint = myModel.computeCrossCovariance(vertices, vertices[0])
69    crossCovCol = crossCov.reshape(crossCov.getNbRows(), dimension)
70    ott.assert_almost_equal(crossCovSamplePoint, crossCovCol, 1e-14, 1e-14,
71                            "in " + myModel.getClassName() + "::computeCrossCovarianceSamplePoint")
72
73    if test_partial_grad:
74        grad = myModel.partialGradient(x1, x2)
75
76        if (dimension == 1):
77            gradfd = ot.Matrix(inputDimension, 1)
78            for j in range(inputDimension):
79                x1_g = ot.Point(x1)
80                x1_d = ot.Point(x1)
81                x1_g[j] = x1_d[j] + eps
82                x1_d[j] = x1_d[j] - eps
83                gradfd[j, 0] = (myModel.computeAsScalar(
84                    x1_g, x2) - myModel.computeAsScalar(x1_d, x2)) / (2 * eps)
85        else:
86            gradfd = ot.Matrix(inputDimension, dimension * dimension)
87            covarianceX1X2 = myModel(x1, x2)
88            centralValue = ot.Point(covarianceX1X2.getImplementation())
89            # Loop over the shifted points
90            for i in range(inputDimension):
91                currentPoint = ot.Point(x1)
92                currentPoint[i] += eps
93                localCovariance = myModel(currentPoint, x2)
94                currentValue = ot.Point(localCovariance.getImplementation())
95                for j in range(currentValue.getSize()):
96                    gradfd[i, j] = (currentValue[j] - centralValue[j]) / eps
97
98        ott.assert_almost_equal(grad, gradfd, 1e-5, 1e-5,
99                                "in " + myModel.getClassName() + " grad")
100
101
102def test_scalar_model(myModel, x1=None, x2=None):
103    if x1 is None and x2 is None:
104        x1 = 2.0
105        x2 = -3.0
106    # Check that computeAsScalar(Scalar) == computeAsScalar(Point)
107    ott.assert_almost_equal(myModel.computeAsScalar(
108        [x1], [x2]), myModel.computeAsScalar(x1, x2), 1.0e-14, 1.0e-14)
109
110    # Gradient testing
111    eps = 1e-5
112
113    grad = myModel.partialGradient([x1], [x2])[0, 0]
114
115    x1_g = x1 + eps
116    x1_d = x1 - eps
117    gradfd = (myModel.computeAsScalar(x1_g, x2) -
118              myModel.computeAsScalar(x1_d, x2)) / (2.0 * eps)
119    ott.assert_almost_equal(gradfd, grad, 1e-5, 1e-5)
120
121
122inputDimension = 2
123
124# 1) SquaredExponential
125myModel = ot.SquaredExponential([2.0], [3.0])
126ott.assert_almost_equal(myModel.getScale(), [2], 0, 0)
127ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
128test_model(myModel)
129
130myModel = ot.SquaredExponential([2.0] * inputDimension, [3.0])
131ott.assert_almost_equal(myModel.getScale(), [2, 2], 0, 0)
132ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
133test_model(myModel)
134
135
136# 2) GeneralizedExponential
137myModel = ot.GeneralizedExponential([2.0], [3.0], 1.5)
138ott.assert_almost_equal(myModel.getScale(), [2], 0, 0)
139ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
140ott.assert_almost_equal(myModel.getP(), 1.5, 0, 0)
141test_model(myModel)
142
143myModel = ot.GeneralizedExponential([2.0] * inputDimension, [3.0], 1.5)
144ott.assert_almost_equal(myModel.getScale(), [2, 2], 0, 0)
145ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
146ott.assert_almost_equal(myModel.getP(), 1.5, 0, 0)
147test_model(myModel)
148
149# 3) AbsoluteExponential
150myModel = ot.AbsoluteExponential([2.0], [3.0])
151ott.assert_almost_equal(myModel.getScale(), [2], 0, 0)
152ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
153test_model(myModel)
154
155myModel = ot.AbsoluteExponential([2.0] * inputDimension, [3.0])
156ott.assert_almost_equal(myModel.getScale(), [2, 2], 0, 0)
157ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
158test_model(myModel)
159
160# 4) MaternModel
161myModel = ot.MaternModel([2.0], [3.0], 1.5)
162ott.assert_almost_equal(myModel.getScale(), [2], 0, 0)
163ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
164ott.assert_almost_equal(myModel.getNu(), 1.5, 0, 0)
165test_model(myModel)
166
167myModel = ot.MaternModel([2.0] * inputDimension, [3.0], 1.5)
168ott.assert_almost_equal(myModel.getScale(), [2, 2], 0, 0)
169ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
170ott.assert_almost_equal(myModel.getNu(), 1.5, 0, 0)
171test_model(myModel)
172
173# 5) ExponentiallyDampedCosineModel
174myModel = ot.ExponentiallyDampedCosineModel([2.0], [3.0], 1)
175ott.assert_almost_equal(myModel.getScale(), [2], 0, 0)
176ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
177ott.assert_almost_equal(myModel.getFrequency(), 1, 0, 0)
178test_model(myModel)
179myModel.setFrequency(3)
180ott.assert_almost_equal(myModel.getFrequency(), 3, 0, 0)
181
182myModel = ot.ExponentiallyDampedCosineModel([2.0] * inputDimension, [3.0], 1)
183ott.assert_almost_equal(myModel.getScale(), [2, 2], 0, 0)
184ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
185ott.assert_almost_equal(myModel.getFrequency(), 1, 0, 0)
186test_model(myModel)
187
188# 6) SphericalModel
189myModel = ot.SphericalModel([2.0], [3.0], 4.5)
190ott.assert_almost_equal(myModel.getScale(), [2], 0, 0)
191ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
192ott.assert_almost_equal(myModel.getRadius(), 4.5, 0, 0)
193test_model(myModel)
194
195myModel = ot.SphericalModel([2.0] * inputDimension, [3.0], 4.5)
196ott.assert_almost_equal(myModel.getScale(), [2, 2], 0, 0)
197ott.assert_almost_equal(myModel.getAmplitude(), [3], 0, 0)
198ott.assert_almost_equal(myModel.getRadius(), 4.5, 0, 0)
199test_model(myModel)
200myModel.setRadius(1.5)
201ott.assert_almost_equal(myModel.getRadius(), 1.5, 0, 0)
202
203# 7) FractionalBrownianMotionModel
204myModel = ot.FractionalBrownianMotionModel(2.0, 3.0, 0.25)
205test_model(myModel)
206
207# 8) DiracCovarianceModel
208myModel = ot.DiracCovarianceModel()
209# Should not check the partialGradient Dirac model
210test_model(myModel, test_partial_grad=False)
211
212amplitude = [1.5 + 2.0 * k for k in range(2)]
213dimension = 2
214spatialCorrelation = ot.CorrelationMatrix(dimension)
215for j in range(dimension):
216    for i in range(j + 1, dimension):
217        spatialCorrelation[i, j] = (
218            i + 1.0) / dimension - (j + 1.0) / dimension
219myModel = ot.DiracCovarianceModel(
220    inputDimension, amplitude, spatialCorrelation)
221ott.assert_almost_equal(myModel.getScale(), [1, 1], 0, 0)
222ott.assert_almost_equal(myModel.getAmplitude(), amplitude, 0, 0)
223test_model(myModel, test_partial_grad=False, x1=[0.5, 0.0], x2=[0.5, 0.0])
224
225# 9) StationaryFunctionalCovarianceModel
226rho = ot.SymbolicFunction(['tau'], ['exp(-abs(tau))*cos(2*pi_*abs(tau))'])
227myModel = ot.StationaryFunctionalCovarianceModel([1.0], [1.0], rho)
228ott.assert_almost_equal(myModel.getScale(), [1], 0, 0)
229ott.assert_almost_equal(myModel.getAmplitude(), [1], 0, 0)
230test_model(myModel)
231
232# 10) ProductCovarianceModel
233myModel = ot.ProductCovarianceModel()
234test_model(myModel)
235
236cov1 = ot.AbsoluteExponential([2.0], [3.0])
237cov2 = ot.SquaredExponential([2.0], [3.0])
238myModel = ot.ProductCovarianceModel([cov1, cov2])
239test_model(myModel)
240ott.assert_almost_equal(myModel.getScale(), [2, 2], 0, 0)
241ott.assert_almost_equal(myModel.getAmplitude(), [9], 0, 0)
242point = [0.50, -6]
243x = [point[0]]
244y = [point[1]]
245ott.assert_almost_equal(myModel.computeAsScalar(point), cov1.computeAsScalar(
246    x) * cov2.computeAsScalar(y), 1.0e-15, 1.0e-15)
247
248# 11) TensorizedCovarianceModel
249
250# Collection ==> add covariance models
251myAbsoluteExponential = ot.AbsoluteExponential([2.0] * inputDimension, [3.0])
252mySquaredExponential = ot.SquaredExponential([2.0] * inputDimension, [3.0])
253myGeneralizedExponential = ot.GeneralizedExponential(
254    [2.0] * inputDimension, [3.0], 1.5)
255# Build TensorizedCovarianceModel with scale = [1,..,1]
256# Tensorized ignore scales
257myModel = ot.TensorizedCovarianceModel(
258    [myAbsoluteExponential, mySquaredExponential, myGeneralizedExponential])
259ott.assert_almost_equal(myModel.getScale(), [1, 1], 0, 0)
260ott.assert_almost_equal(myModel.getAmplitude(), [3, 3, 3], 0, 0)
261test_model(myModel)
262
263# Define new scale
264scale = [2.5, 1.5]
265myModel.setScale(scale)
266ott.assert_almost_equal(myModel.getScale(), [2.5, 1.5], 0, 0)
267ott.assert_almost_equal(myModel.getAmplitude(), [3, 3, 3], 0, 0)
268test_model(myModel)
269
270# 12) Testing 1d in/out dimension & stationary
271
272# Test scalar input models
273coll = [ot.AbsoluteExponential(), ot.SquaredExponential(),
274        ot.GeneralizedExponential()]
275coll += [ot.MaternModel(), ot.SphericalModel(),
276         ot.ExponentiallyDampedCosineModel()]
277for model in coll:
278    test_scalar_model(model)
279
280# 13) Isotropic covariance model
281
282scale = 3.5
283amplitude = 1.5
284myOneDimensionalKernel = ot.SquaredExponential([scale], [amplitude])
285myIsotropicKernel = ot.IsotropicCovarianceModel(
286    myOneDimensionalKernel, inputDimension)
287
288# Test consistency of isotropic model with underlying 1D kernel
289ott.assert_almost_equal(myIsotropicKernel.getAmplitude()[
290                        0], amplitude, 1e-12, 0.0)
291ott.assert_almost_equal(myIsotropicKernel.getScale()[0], scale, 1e-12, 0.0)
292ott.assert_almost_equal(myIsotropicKernel.getKernel().getAmplitude()[
293                        0], amplitude, 1e-12, 0.0)
294ott.assert_almost_equal(myIsotropicKernel.getKernel().getScale()[
295                        0], scale, 1e-12, 0.0)
296
297# Standard tests applied
298test_model(myIsotropicKernel)
299
300# Test consistency of isotropic kernel's discretization
301inputVector = ot.Point([0.3, 1.7])
302inputVectorNorm = ot.Point([inputVector.norm()])
303ott.assert_almost_equal(myOneDimensionalKernel(inputVectorNorm)[
304                        0, 0], 1.992315565746, 1e-12, 0.0)
305ott.assert_almost_equal(myIsotropicKernel(inputVector)[
306                        0, 0], 1.992315565746, 1e-12, 0.0)
307inputSample = ot.Sample([ot.Point(2), inputVector])
308inputSampleNorm = ot.Sample([ot.Point(1), inputVectorNorm])
309oneDimensionalCovMatrix = myOneDimensionalKernel.discretize(inputSampleNorm)
310isotropicCovMatrix = myIsotropicKernel.discretize(inputSample)
311ott.assert_almost_equal(
312    oneDimensionalCovMatrix[0, 0], 2.250000000002, 1e-12, 0.0)
313ott.assert_almost_equal(
314    oneDimensionalCovMatrix[1, 1], 2.250000000002, 1e-12, 0.0)
315ott.assert_almost_equal(isotropicCovMatrix[0, 0], 2.250000000002, 1e-12, 0.0)
316ott.assert_almost_equal(isotropicCovMatrix[1, 1], 2.250000000002, 1e-12, 0.0)
317ott.assert_almost_equal(
318    oneDimensionalCovMatrix[0, 1], 1.992315565746, 1e-12, 0.0)
319ott.assert_almost_equal(isotropicCovMatrix[0, 1], 1.992315565746, 1e-12, 0.0)
320
321# Exponential covariance model
322inputDimension = 2
323scale = [4, 5]
324spatialCovariance = ot.CovarianceMatrix(inputDimension)
325spatialCovariance[0, 0] = 4
326spatialCovariance[1, 1] = 5
327spatialCovariance[1, 0] = 1.2
328myModel = ot.ExponentialModel(scale, spatialCovariance)
329test_model(myModel)
330# assert that spatialCovariance is taken into account
331checkDiag = spatialCovariance.isDiagonal() == myModel.isDiagonal()
332if (not checkDiag):
333    raise Exception(
334        "isDiagonal differ between spatial covariance & covariance model")
335rho = spatialCovariance[1, 0] / \
336    sqrt(spatialCovariance[0, 0] * spatialCovariance[1, 1])
337ott.assert_almost_equal(myModel.getOutputCorrelation(
338)[0, 1], rho, 0, 0, "in ExponentialModel correlation")
339
340# Kronecker covariance model
341
342# rho correlation
343scale = [4, 5]
344rho = ot.GeneralizedExponential(scale, 1)
345
346# Amplitude values
347amplitude = [1, 2]
348myModel = ot.KroneckerCovarianceModel(rho, amplitude)
349test_model(myModel)
350
351outputCorrelation = ot.CorrelationMatrix(2)
352outputCorrelation[0, 1] = 0.8
353myModel = ot.KroneckerCovarianceModel(rho, amplitude, outputCorrelation)
354test_model(myModel)
355
356outputCovariance = ot.CovarianceMatrix(2)
357outputCovariance[0, 0] = 4.0
358outputCovariance[1, 1] = 5.0
359outputCovariance[1, 0] = 1.2
360
361myModel = ot.KroneckerCovarianceModel(rho, outputCovariance)
362test_model(myModel)
363
364# New kronecker model involving isotropic cov model
365rho = ot.IsotropicCovarianceModel(ot.MaternModel(), 3)
366outputCorrelation = ot.CorrelationMatrix(2)
367outputCorrelation[0, 1] = 0.8
368amplitude = [1, 2]
369myModel = ot.KroneckerCovarianceModel(rho, amplitude, outputCorrelation)
370test_model(myModel)
371ott.assert_almost_equal(myModel.getInputDimension(), 3,
372                        0, 0, "in kronecker dimension check")
373ott.assert_almost_equal(myModel.getScale(), [
374                        1], 0, 0, "in kronecker scale check")
375# full param size = 5 (scale(1), amplitude(2), spatialCorrelation(1), Matern nu(1))
376ott.assert_almost_equal(myModel.getFullParameter(), [
377                        1, 1, 2, 0.8, 1.5], 0, 0, "in kronecker full param check")
378ott.assert_almost_equal(myModel.getFullParameter(
379).getSize(), 5, 0, 0, "in kronecker param size check")
380ott.assert_almost_equal(myModel.getFullParameterDescription(
381).getSize(), 5, 0, 0, "in kronecker param description size check")
382ott.assert_almost_equal(myModel.getActiveParameter(), [
383                        0, 1, 2], "in kronecker active param check")
384myModel.setFullParameter([2, 1, 2, .5, 2.5])
385ott.assert_almost_equal(myModel.getFullParameter(), [
386                        2, 1, 2, .5, 2.5], 0, 0, "in kronecker param check")
387myModel.setActiveParameter([0, 1, 2, 4])
388ott.assert_almost_equal(myModel.getActiveParameter(), [
389                        0, 1, 2, 4], "in kronecker active param check")
390# Now we should get all values except correlation
391ott.assert_almost_equal(myModel.getParameter(), [
392                        2, 1, 2, 2.5], 0, 0, "in kronecker param check")
393assert myModel.getFullParameterDescription(
394) == ["scale_0", "amplitude_0", "amplitude_1", "R_1_0", "nu"]
395