Skip to content

Commit

Permalink
Merge branch 'master' into master
Browse files Browse the repository at this point in the history
  • Loading branch information
@relf
relf authored Dec 3, 2024
2 parents cd12c86 + 4121b9e commit cf83bb7
Show file tree
Hide file tree
Showing 5 changed files with 2,681 additions and 2,522 deletions.
1 change: 1 addition & 0 deletions .github/workflows/tests_coverage.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ jobs:
run: |
ruff --version
ruff check .
ruff format --check
- name: Test with pytest and coverage
run: |
Expand Down
206 changes: 129 additions & 77 deletions smt/applications/mfck.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,15 @@
Conference on Neural Information Processing Systems.
"""

#import warnings
# import warnings
import numpy as np
from scipy.linalg import solve_triangular
from scipy import optimize
from smt.sampling_methods import LHS
from smt.utils.kriging import (differences, componentwise_distance)
from smt.utils.kriging import differences, componentwise_distance
from smt.surrogate_models.krg_based import KrgBased


class MFCK(KrgBased):
def _initialize(self):
super()._initialize()
Expand All @@ -37,7 +38,7 @@ def _initialize(self):
)
declare(
"rho_bounds",
[-5., 5.],
[-5.0, 5.0],
types=(list, np.ndarray),
desc="Bounds for the rho parameter used in the autoregressive model",
)
Expand All @@ -54,8 +55,12 @@ def _initialize(self):
desc="Bounds for the variance parameter",
)

self.options["nugget"] = 1e-9 # Incresing the nugget for numerical stability reasons
self.options["hyper_opt"] = "Cobyla" # MFCK doesn't support gradient-based optimizers
self.options["nugget"] = (
1e-9 # Incresing the nugget for numerical stability reasons
)
self.options["hyper_opt"] = (
"Cobyla" # MFCK doesn't support gradient-based optimizers
)

def train(self):
"""
Expand All @@ -74,20 +79,23 @@ def train(self):
i = i + 1
xt.append(self.training_points[None][0][0])
yt.append(self.training_points[None][0][1])
self.lvl = i+1
self.lvl = i + 1
self.X = xt
self.y = np.vstack(yt)
self._check_param()
self.nx = 1 # Forcing this in order to consider isotropic kernels (i.e., only one lengthscale)
self.nx = 1 # Forcing this in order to consider isotropic kernels (i.e., only one lengthscale)
if self.lvl == 1:
# For a single level, initialize theta_ini, lower_bounds, and
# upper_bounds with consistent shapes
theta_ini = np.hstack((self.options["sigma0"],
self.options["theta0"][0])) # Kernel variance + theta0
lower_bounds = np.hstack((self.options["sigma_bounds"][0],
self.options["theta_bounds"][0]))
upper_bounds = np.hstack((self.options["sigma_bounds"][1],
self.options["theta_bounds"][1]))
theta_ini = np.hstack(
(self.options["sigma0"], self.options["theta0"][0])
) # Kernel variance + theta0
lower_bounds = np.hstack(
(self.options["sigma_bounds"][0], self.options["theta_bounds"][0])
)
upper_bounds = np.hstack(
(self.options["sigma_bounds"][1], self.options["theta_bounds"][1])
)
theta_ini = np.log10(theta_ini)
lower_bounds = np.log10(lower_bounds)
upper_bounds = np.log10(upper_bounds)
Expand All @@ -96,47 +104,73 @@ def train(self):
for lvl in range(self.lvl):
if lvl == 0:
# Initialize theta_ini for level 0
theta_ini = np.hstack((self.options["sigma0"],
self.options["theta0"][0])) # Variance + initial theta values
lower_bounds = np.hstack((self.options["sigma_bounds"][0],
np.full(self.nx, self.options["theta_bounds"][0])))
upper_bounds = np.hstack((self.options["sigma_bounds"][1],
np.full(self.nx, self.options["theta_bounds"][1])))
theta_ini = np.hstack(
(self.options["sigma0"], self.options["theta0"][0])
) # Variance + initial theta values
lower_bounds = np.hstack(
(
self.options["sigma_bounds"][0],
np.full(self.nx, self.options["theta_bounds"][0]),
)
)
upper_bounds = np.hstack(
(
self.options["sigma_bounds"][1],
np.full(self.nx, self.options["theta_bounds"][1]),
)
)
# Apply log10 to theta_ini and bounds
nb_params = len(self.options["theta0"])
theta_ini[:nb_params+1] = np.log10(theta_ini[:nb_params+1])
lower_bounds[:nb_params+1] = np.log10(lower_bounds[:nb_params+1])
upper_bounds[:nb_params+1] = np.log10(upper_bounds[:nb_params+1])
theta_ini[: nb_params + 1] = np.log10(theta_ini[: nb_params + 1])
lower_bounds[: nb_params + 1] = np.log10(
lower_bounds[: nb_params + 1]
)
upper_bounds[: nb_params + 1] = np.log10(
upper_bounds[: nb_params + 1]
)

elif lvl > 0:
# For additional levels, append to theta_ini, lower_bounds, and upper_bounds
thetat = np.hstack((self.options["sigma0"],
self.options["theta0"][0]))
lower_boundst = np.hstack((self.options["sigma_bounds"][0],
np.full(self.nx, self.options["theta_bounds"][0])))
upper_boundst = np.hstack((self.options["sigma_bounds"][1],
np.full(self.nx, self.options["theta_bounds"][1])))
thetat = np.hstack(
(self.options["sigma0"], self.options["theta0"][0])
)
lower_boundst = np.hstack(
(
self.options["sigma_bounds"][0],
np.full(self.nx, self.options["theta_bounds"][0]),
)
)
upper_boundst = np.hstack(
(
self.options["sigma_bounds"][1],
np.full(self.nx, self.options["theta_bounds"][1]),
)
)
# Apply log10 to the newly added values
thetat = np.log10(thetat)
lower_boundst = np.log10(lower_boundst)
upper_boundst = np.log10(upper_boundst)
# Append to theta_ini, lower_bounds, and upper_bounds
theta_ini = np.hstack([theta_ini, thetat,self.options["rho0"]])
theta_ini = np.hstack([theta_ini, thetat, self.options["rho0"]])
lower_bounds = np.hstack([lower_bounds, lower_boundst])
upper_bounds = np.hstack([upper_bounds, upper_boundst])
# Finally, append the rho bounds
lower_bounds = np.hstack([lower_bounds, self.options["rho_bounds"][0]])
upper_bounds = np.hstack([upper_bounds, self.options["rho_bounds"][1]])
lower_bounds = np.hstack(
[lower_bounds, self.options["rho_bounds"][0]]
)
upper_bounds = np.hstack(
[upper_bounds, self.options["rho_bounds"][1]]
)

theta_ini = theta_ini[:].T
x_opt = theta_ini

if self.options["hyper_opt"] == "Cobyla":
if self.options["n_start"] > 1:
sampling = LHS(
xlimits = np.stack((lower_bounds, upper_bounds), axis=1),
criterion = "ese",
random_state = 0,
xlimits=np.stack((lower_bounds, upper_bounds), axis=1),
criterion="ese",
random_state=0,
)
theta_lhs_loops = sampling(self.options["n_start"])
theta0 = np.vstack((theta_ini, theta_lhs_loops))
Expand All @@ -151,11 +185,9 @@ def train(self):
optimal_theta_res_loop = optimize.minimize(
self.neg_log_likelihood_scipy,
theta0[j, :],
method = "COBYLA",
constraints = [
{"fun": con, "type": "ineq"} for con in constraints
],
options = {
method="COBYLA",
constraints=[{"fun": con, "type": "ineq"} for con in constraints],
options={
"rhobeg": 0.1,
"tol": 1e-4,
"maxiter": 200,
Expand Down Expand Up @@ -190,9 +222,9 @@ def train(self):
f"The optimizer {self.options['hyper_opt']} is not available"
)

x_opt[0:2] = 10**(x_opt[0:2])
x_opt[2::3] = 10**(x_opt[2::3])
x_opt[3::3] = 10**(x_opt[3::3])
x_opt[0:2] = 10 ** (x_opt[0:2])
x_opt[2::3] = 10 ** (x_opt[2::3])
x_opt[3::3] = 10 ** (x_opt[3::3])
self.optimal_theta = x_opt

def eta(self, j, jp, rho):
Expand Down Expand Up @@ -236,8 +268,9 @@ def compute_cross_K(self, x, xp, L, Lp, param):
cov_gamma_j = self._compute_K(x, xp, param0)
else:
# Cov(γ_j(x), γ_j(x')) using the kernel
cov_gamma_j = self._compute_K(x, xp,
[sigmas_gamma[j-1], ls_gamma[j-1]])
cov_gamma_j = self._compute_K(
x, xp, [sigmas_gamma[j - 1], ls_gamma[j - 1]]
)
# Add to the value of the covariance
cov_value += eta_j_l * eta_j_lp * cov_gamma_j

Expand All @@ -261,34 +294,46 @@ def predict_all_levels(self, x):
means = []
covariances = []
if self.lvl == 1:
k_XX = self._compute_K(self.X[0], self.X[0],self.optimal_theta[0:2])
k_XX = self._compute_K(self.X[0], self.X[0], self.optimal_theta[0:2])
k_xX = self._compute_K(x, self.X[0], self.optimal_theta[0:2])
k_xx = self._compute_K(x,x,self.optimal_theta[0:2])
k_xx = self._compute_K(x, x, self.optimal_theta[0:2])
# To be adapted using the Cholesky decomposition
k_XX_inv = np.linalg.inv(k_XX + self.options["nugget"]*np.eye(k_XX.shape[0]))
means.append( np.dot(k_xX, np.matmul(k_XX_inv, self.y)))
covariances.append(k_xx - np.matmul(k_xX,
np.matmul(k_XX_inv,
k_xX.transpose())))
k_XX_inv = np.linalg.inv(
k_XX + self.options["nugget"] * np.eye(k_XX.shape[0])
)
means.append(np.dot(k_xX, np.matmul(k_XX_inv, self.y)))
covariances.append(
k_xx - np.matmul(k_xX, np.matmul(k_XX_inv, k_xX.transpose()))
)

else:
L = np.linalg.cholesky(self.K + self.options["nugget"] * np.eye(self.K.shape[0]))
L = np.linalg.cholesky(
self.K + self.options["nugget"] * np.eye(self.K.shape[0])
)
k_xX = []
for ind in range(self.lvl):
k_xx = self.compute_cross_K(x, x, ind, ind,self.optimal_theta)
k_xx = self.compute_cross_K(x, x, ind, ind, self.optimal_theta)
for j in range(self.lvl):
if ind >= j:
k_xX.append(self.compute_cross_K(self.X[j], x, ind, j,self.optimal_theta))
k_xX.append(
self.compute_cross_K(
self.X[j], x, ind, j, self.optimal_theta
)
)
else:
k_xX.append(self.compute_cross_K(self.X[j], x, j, ind,self.optimal_theta))
k_xX.append(
self.compute_cross_K(
self.X[j], x, j, ind, self.optimal_theta
)
)

beta1 = solve_triangular(L, np.vstack(k_xX), lower=True)
alpha1 = solve_triangular(L, self.y, lower=True)
means.append(np.dot(beta1.T, alpha1))
covariances.append(k_xx - np.dot(beta1.T, beta1))
k_xX.clear()

return means,covariances
return means, covariances

def _predict(self, x):
"""
Expand All @@ -306,61 +351,65 @@ def _predict(self, x):
"""
means, covariances = self.predict_all_levels(x)

return means[self.lvl-1], covariances[self.lvl-1]
return means[self.lvl - 1], covariances[self.lvl - 1]

def neg_log_likelihood(self, param, grad=None):

if self.lvl == 1:
self.K = self._compute_K(self.X[0], self.X[0], param[0:2])
else:
self.K = self.compute_blockwise_K(param)

L = np.linalg.cholesky(self.K + self.options["nugget"]*np.eye(self.K.shape[0]))
L = np.linalg.cholesky(
self.K + self.options["nugget"] * np.eye(self.K.shape[0])
)
beta = solve_triangular(L, self.y, lower=True)
NMLL = 1/2*(2*np.sum(np.log(np.diag(L))) + np.dot(beta.T,beta))
nmll, = NMLL[0]
NMLL = 1 / 2 * (2 * np.sum(np.log(np.diag(L))) + np.dot(beta.T, beta))
(nmll,) = NMLL[0]
return nmll

def neg_log_likelihood_scipy(self, param):
"""
Likelihood for Cobyla-scipy (SMT) optimizer
"""
param = np.array(param, copy=True)
param[0:2] = 10**(param[0:2])
param[2::3] = 10**(param[2::3])
param[3::3] = 10**(param[3::3])
param[0:2] = 10 ** (param[0:2])
param[2::3] = 10 ** (param[2::3])
param[3::3] = 10 ** (param[3::3])
return self.neg_log_likelihood(param)

def neg_log_likelihood_nlopt(self, param, grad=None):
"""
Likelihood for nlopt optimizers
"""
param = np.array(param, copy=True)
param[0:2] = 10**(param[0:2])
param[2::3] = 10**(param[2::3])
param[3::3] = 10**(param[3::3])
param[0:2] = 10 ** (param[0:2])
param[2::3] = 10 ** (param[2::3])
param[3::3] = 10 ** (param[3::3])
return self.neg_log_likelihood(param, grad)


def compute_blockwise_K(self, param):
K_block = {}
n = self.y.shape[0]
for jp in range(self.lvl):
for j in range(self.lvl):
if jp >= j:
K_block[str(jp)+str(j)] = self.compute_cross_K(self.X[j],
self.X[jp],
jp, j, param)
K_block[str(jp) + str(j)] = self.compute_cross_K(
self.X[j], self.X[jp], jp, j, param
)

K = np.zeros((n, n))
row_init, col_init = 0, 0
for j in range(self.lvl):
col_init = row_init
for jp in range(j, self.lvl):
r, c = K_block[str(jp)+str(j)].shape
K[row_init:row_init+r, col_init:col_init+c] = K_block[str(jp)+str(j)]
r, c = K_block[str(jp) + str(j)].shape
K[row_init : row_init + r, col_init : col_init + c] = K_block[
str(jp) + str(j)
]
if j != jp:
K[col_init:col_init+c, row_init:row_init+r] = K_block[str(jp)+str(j)].T
K[col_init : col_init + c, row_init : row_init + r] = K_block[
str(jp) + str(j)
].T
col_init += c
row_init += r
return K
Expand All @@ -372,9 +421,12 @@ def _compute_K(self, A: np.ndarray, B: np.ndarray, param):
"""
# Compute pairwise componentwise L1-distances between A and B
dx = differences(A, B)
d = componentwise_distance(dx, self.options["corr"],
self.X[0].shape[1],
power=self.options["pow_exp_power"])
d = componentwise_distance(
dx,
self.options["corr"],
self.X[0].shape[1],
power=self.options["pow_exp_power"],
)
self.corr.theta = np.full(self.X[0].shape[1], param[1])
r = self.corr(d)
R = r.reshape(A.shape[0], B.shape[0])
Expand Down
Loading

0 comments on commit cf83bb7

Please sign in to comment.