# coding: utf-8

# # Exemple de krigeage 1D
# 
# # Références
# 
# Metamodeling with Gaussian processes, Bertrand Iooss, EDF R&D, 2014, www.gdr-mascotnum.fr/media/sssamo14_iooss.pdf

#
import numpy as np
import openturns as ot
import pylab as pl

#
n_test = 100
x_test_coord = np.linspace(0,12,n_test)
x_test = ot.Sample(x_test_coord,1)
y_test = np.sin(x_test)

# Base d'apprentissage de 7 points.
#
x_train_coord = [1.,3.,4.,6.,7.9,11., 11.5]
x_train = ot.Sample(x_train_coord,1)
y_train = np.sin(x_train)
n_train = len(x_train)
n_train

#
dimension = 1
basis = ot.ConstantBasisFactory(dimension).build()
covarianceModel = ot.SquaredExponential([1.]*dimension, [5.0])
algo = ot.KrigingAlgorithm(x_train, y_train, covarianceModel, basis)
algo.setOptimizeParameters(False)
algo.run()
result = algo.getResult()
krigeageMM = result.getMetaModel()

#
krigeageMM = result.getMetaModel()
y_test_MM = krigeageMM(x_test)

#
level = 0.95
alpha = (1-level)/2
sampleQinf = ot.Sample(n_test,1)
sampleQsup = ot.Sample(n_test,1)
for i in range(n_test):
    xi = x_test[i]
    var = result.getConditionalCovariance(xi)[0,0]
    yi = y_test_MM[i][0]
    N = ot.Normal(yi,np.sqrt(var))
    q025 = N.computeQuantile(alpha)[0]
    sampleQinf[i,0] = q025
    q975 = N.computeQuantile(1-alpha)[0]
    sampleQsup[i,0] = q975

#
pl.plot(x_test,y_test,"k-",label="Exact")
pl.plot(x_test,y_test_MM,"b--",label="Krigeage")
pl.plot(x_train,y_train,"ro",label="Apprentissage")
pl.plot(x_test,sampleQinf,"g-",label="%.2f%%" % (100*level))
pl.plot(x_test,sampleQsup,"g-")
pl.title("sin")
pl.xlabel("X")
pl.ylabel("Y")
pl.legend()