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