test_loic.py

# -*- coding: utf-8 -*-
"""
Created on Mon Nov 23 14:14:19 2020

@author: lbrevaul
"""

import openturns as ot
import math as m
import numpy as np

##### function mixte
def fun_mixte(X):
    xx = X[0]
    z = X[1]
    
    if z==0:
        y = np.sin(7*xx)
    else:
        y = 2*np.sin(7*xx)
    return y


XX_input = np.array([[0.1,0],
          [0.32,0],
          [0.6,0],
          [0.9,0],
          [0.07,1],
          [0.1,1],
          [0.4,1],
          [0.5,1],
          [0.85,1]])
y_output = np.zeros((len(XX_input),1))
for i in range(len(XX_input)):
    y_output[i] = fun_mixte(XX_input[i,:])



def C(s, t):
    return m.exp( -4.0 * abs(s - t) / (1 + (s * s + t * t)))


N = 32
a = 4.0
myMesh = ot.IntervalMesher([N]).build(ot.Interval(-a, a))


myCovariance = ot.CovarianceMatrix(myMesh.getVerticesNumber())
for k in range(myMesh.getVerticesNumber()):
    t = myMesh.getVertices()[k]
    for l in range(k + 1):
        s = myMesh.getVertices()[l]
        myCovariance[k, l] = C(s[0], t[0])
        
        
        
#covModel_discrete = ot.UserDefinedCovarianceModel(myMesh, myCovariance)

def Gower_fun(X):
    tau = X[0]
    sigma = X[1]
    theta = X[2]
    if tau == 0:
        output = [sigma**2*np.exp(-tau/theta)]
    else:
        tau_ = 1
        output = [sigma**2*np.exp(-tau_/theta)]
    return output


def Gower_fun(X):
    tau = X[0]
    sigma = 1
    theta = 1
    if tau == 0:
        output = [sigma**2*np.exp(-tau/theta)]
    else:
        tau_ = 1
        output = [sigma**2*np.exp(-tau_/theta)]
    return output

f_ = ot.PythonFunction(1, 1, Gower_fun)
f_ = ot.SymbolicFunction(["tau", "theta", "sigma"], ["(tau!=0) * exp(-1/theta) * sigma * sigma +  (tau==0) * exp(0) * sigma * sigma"])
rho = ot.ParametricFunction(f_, [1, 2], [0.2, 0.3])
covModel_discrete = ot.StationaryFunctionalCovarianceModel([1.0], [1.0], rho)

inputDimension = 1
covModel_continuous = ot.SquaredExponential([1.0],[1.0])

covarianceModel = ot.ProductCovarianceModel([covModel_continuous, covModel_discrete])



dimension = 2
basis = ot.ConstantBasisFactory(dimension).build()

ot.ResourceMap.SetAsBool('GeneralLinearModelAlgorithm-UseAnalyticalAmplitudeEstimate', True)
algo = ot.KrigingAlgorithm(XX_input, y_output, covarianceModel, basis)

#solver_kriging = ot.Cobyla()
#algo.setOptimizationAlgorithm(solver_kriging)
#algo.setOptimizationBounds(ot.Interval([0.1,0.1,0.1], [1e2,1e2,1e2]))
algo.setOptimizationBounds(ot.Interval([0.1,0.1], [1e2,1e2]))
#ot.Log.Show(ot.Log.ALL)
#ot.TBB.Disable()
algo.run()