Hi,
Thanks for your replay and sorry if I was not clear. It seems that I can’t upload attachments as a new user. So I put below the script to obtain the first graph above. To obtain the second one, you just need to change the value of factor to 1000 in the definition of the fun_ref function. Let me know if I can provide the script in a better way.
import numpy as np
import openturns as ot
import matplotlib.pyplot as plt
# %% definitions of functions for kriging calibration
def get_param(X):
# return objects needed for normalization and standardization
vec_min_col = np.amin(X,axis=0)
vec_max_col = np.amax(X,axis=0)
vec_mean_col = np.mean(X, axis =0)
vec_std_col = np.std(X, axis =0)
return vec_min_col, vec_max_col, vec_mean_col, vec_std_col
def normalize(X, params, bool_standard = True):
# return normalization or standardization of X using params
vec_min_col = params[0]
vec_max_col = params[1]
vec_mean_col = params[2]
vec_std_col = params[3]
if not(bool_standard):
return (X - vec_min_col)/(vec_max_col-vec_min_col)
else:
return (X-vec_mean_col) / vec_std_col
def kriging_calibration(DoE, Response, cov_str= 'Gaussian', Hyperparam_old=[ ], return_likelihood = False, bool_noise = False, noise_vec = [ ] ):
#Inputs : DoE : DoE used to build the kriging (vector of input points of the model)
# Response : vector of outputs of the model at the points of the DoE
#Outputs : the calibrated Kriging model
dim_space= DoE.shape[1]
Response = np.vstack(Response)
#ot.RandomGenerator.SetSeed(1235)
# --- 1. basis ---
#basis = [] #trend 0
basis = ot.ConstantBasisFactory(dim_space).build() # trend constante
#basis = ot.LinearBasisFactory(dim_space).build() # trend lineaire
#basis = ot.QuadraticBasisFactory(dim_space).build() # trend quadratic
# --- 2. covariance model ---
if cov_str == 'Matern1.5':
cov = ot.MaternModel(np.ones(dim_space), 1.5)
elif cov_str == 'Matern2.5':
cov = ot.MaternModel(np.ones(dim_space), 2.5)
elif cov_str == 'Gaussian':
cov = ot.SquaredExponential(np.ones(dim_space))
#cov.setNuggetFactor(1e-6) # to add a nugget effect ("effet pepite")
#cov.setScale(np.array([1.57188,0.992836,1.27414,2.05575,5,1.32729])) # to set the hyperparameters
# --- 3. kriging algorithm ---
algokriging = ot.KrigingAlgorithm(DoE, Response, cov, basis)
if bool_noise:
algokriging.setNoise(noise_vec)
# --- 4. Optimization ---
# for kriging calibration
nb_multistart_hyperparam = 100
lower_bound_theta = 1e-5
upper_bound_theta = 10.
if type(Hyperparam_old)== list :
distribution = ot.ComposedDistribution( np.ones(dim_space)*ot.Uniform(lower_bound_theta , upper_bound_theta ) )
else :
distrib = []
Lower_bounds = []
Upper_bounds = []
for theta in Hyperparam_old :
Lower_bounds.append(max(lower_bound_theta,theta - 1.))
Upper_bounds.append(min(theta +1.,upper_bound_theta))
distrib.append(ot.Uniform(max(lower_bound_theta,theta - 1.) , min(theta +1.,upper_bound_theta)))
distribution = ot.ComposedDistribution(distrib)
startingPoint = ot.LHSExperiment(distribution, nb_multistart_hyperparam).generate()
#startingPoint = loguniform.rvs(1e-3,10.,size = (10, dim_space))
algokriging.setOptimizationAlgorithm(ot.MultiStart(ot.TNC(), startingPoint))
#algokriging.setOptimizationAlgorithm(ot.MultiStart(ot.NLopt('LD_LBFGS'), startingPoint))
algokriging.setOptimizationBounds(ot.Interval(np.ones( dim_space)*lower_bound_theta, np.ones(dim_space)* upper_bound_theta ))
#algokriging.setOptimizationBounds(ot.Interval(np.array(Lower_bounds), np.array(Upper_bounds)))
algokriging.setOptimizeParameters(True)
# # if we choose not to optimize parameters
#algokriging.setOptimizeParameters(False)
# --- 5. run the algorithm ---
algokriging.run()
# # get the metamodel
krigingResult = algokriging.getResult()
#krigingMeta = krigingResult.getMetaModel()
if return_likelihood:
F = algokriging.getReducedLogLikelihoodFunction()
cov_model = krigingResult.getCovarianceModel()
print(cov_model)
print(cov_model.getParameter())
likelihood = F(cov_model.getParameter()[:dim_space])
return krigingResult, cov_model, likelihood
else :
return krigingResult
# %% toy case
def fun_ref(x):
factor = 1 # !!! or factor = 1000 !!!
return factor*np.sin(x)
# --- DoE ---
nb_point = 6
x_max = np.pi
np.random.seed(1) ; DoE = np.random.random_sample(nb_point)*x_max
params_DoE = get_param(np.array(DoE))
DoE_normalized = np.vstack(normalize(np.array(DoE), params_DoE))
Y_true = fun_ref(DoE)
# -- noise ---
vec_noise = 0.1 * np.ones(nb_point)
# --- kriging calibration ---
krigingResult = kriging_calibration(DoE_normalized, Y_true, cov_str='Matern2.5', bool_noise = True, noise_vec = vec_noise)
# --- kriging prediction at x_plot ---
x_plot = np.linspace(0, x_max, 1000)
x_plot_normalized = np.vstack(normalize(np.array(x_plot), params_DoE))
means = np.squeeze(krigingResult.getConditionalMean(np.vstack(x_plot_normalized)))
variances = np.array(krigingResult.getConditionalMarginalVariance(np.vstack(x_plot_normalized)))
stds = np.squeeze(np.sqrt(variances.clip(min=0)))
# --- plot ---
lb_predict = np.squeeze(means) - 2 * stds
ub_predict = np.squeeze(means) + 2 * stds
plt.plot(x_plot, np.squeeze(means), label = 'km', linestyle= '--', c = 'black')
plt.fill_between(x_plot, lb_predict , ub_predict , alpha = 0.1, color = 'black')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.scatter(DoE, Y_true)
plt.plot(x_plot, fun_ref(x_plot), label = 'sin(x)')
plt.legend()
#plt.savefig('sin')
plt.show()
