Skip to content

Commit

Permalink
Merge 9209049 into 2fd5c5a
Browse files Browse the repository at this point in the history
  • Loading branch information
@Paul-Saves
Paul-Saves committed Feb 17, 2024
2 parents 2fd5c5a + 9209049 commit 03d9ade
Show file tree
Hide file tree
Showing 3 changed files with 275 additions and 6 deletions.
147 changes: 147 additions & 0 deletions smt/applications/tests/test_mixed_integer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,153 @@


class TestMixedInteger(unittest.TestCase):
def test_periodic_mixed(self):
# Objective function
def f_obj(X):
"""
s01 objective
Parameters
----------
point: array_like
point to evaluate
"""
design_space = DesignSpace(
[
FloatVariable(0.0, 1.0),
CategoricalVariable(
values=[
"1",
"2",
"3",
"4",
"5",
"6",
"7",
"8",
"9",
"10",
"11",
"12",
"13",
]
),
]
)

PI = 3.14159265358979323846
fail = False
x = X[0]
# caté 1
c = str(design_space.decode_values(X, i_dv=1)[1])
x = np.abs(x)
c1 = c == "1"
c2 = c == "2"
c3 = c == "3"
c4 = c == "4"
c5 = c == "5"
c6 = c == "6"
c7 = c == "7"
c8 = c == "8"
c9 = c == "9"
c10 = c == "10"
c11 = c == "11"
c12 = c == "12"
c13 = c == "13"

if np.size(c1) == (
np.sum(c1)
+ np.sum(c2)
+ np.sum(c3)
+ np.sum(c4)
+ np.sum(c5)
+ np.sum(c6)
+ np.sum(c7)
+ np.sum(c8)
+ np.sum(c9)
+ np.sum(c10)
+ np.sum(c11)
+ np.sum(c12)
+ np.sum(c13)
):
u = (
1 * c1
+ 2 * c2
+ 3 * c3
+ 4 * c4
+ 5 * c5
+ 6 * c6
+ 7 * c7
+ 8 * c8
+ 9 * c9
+ 10 * c10
+ 11 * c11
+ 12 * c12
+ 13 * c13
)
hu = (PI * (0.4 + u / 15)) * (c10 + c11 + c12 + c13) - u / 20
y = np.cos(3.5 * PI * x + hu)
else:
print("type error")
print(X)
fail = True
return (y, fail)

n_doe = 98

design_space = DesignSpace(
[
FloatVariable(0.0, 1.0),
CategoricalVariable(
values=[
"1",
"2",
"3",
"4",
"5",
"6",
"7",
"8",
"9",
"10",
"11",
"12",
"13",
]
),
]
)

sampling = MixedIntegerSamplingMethod(
LHS,
design_space,
criterion="ese",
random_state=42,
output_in_folded_space=True,
)
xdoe = sampling(n_doe)

y_doe = [f_obj(xdoe[i])[0] for i in range(len(xdoe))]
# Surrogate
for m in MixIntKernelType:
sm = MixedIntegerKrigingModel(
surrogate=KRG(
design_space=design_space,
categorical_kernel=m,
hyper_opt="Cobyla",
theta0=[0.01],
corr="squar_sin_exp",
n_start=5,
),
)

sm.set_training_values(xdoe, np.array(y_doe))
if m in [MixIntKernelType.EXP_HOMO_HSPHERE, MixIntKernelType.CONT_RELAX]:
with self.assertRaises(ValueError):
sm.train()
else:
sm.train()

def test_krg_mixed_3D(self):
design_space = DesignSpace(
[
Expand Down
56 changes: 50 additions & 6 deletions smt/surrogate_models/krg_based.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
quadratic,
pow_exp,
squar_exp,
squar_sin_exp,
abs_exp,
act_exp,
cross_distances,
Expand Down Expand Up @@ -61,6 +62,7 @@ class KrgBased(SurrogateModel):
"pow_exp": pow_exp,
"abs_exp": abs_exp,
"squar_exp": squar_exp,
"squar_sin_exp": squar_sin_exp,
"act_exp": act_exp,
"matern52": matern52,
"matern32": matern32,
Expand Down Expand Up @@ -217,7 +219,12 @@ def _final_initialize(self):
# initialize default power values
if self.options["corr"] == "squar_exp":
self.options["pow_exp_power"] = 2.0
elif self.options["corr"] in ["abs_exp", "matern32", "matern52"]:
elif self.options["corr"] in [
"abs_exp",
"squar_sin_exp",
"matern32",
"matern52",
]:
self.options["pow_exp_power"] = 1.0

# Check the pow_exp_power is >0 and <=2
Expand Down Expand Up @@ -631,6 +638,7 @@ def _matrix_data_corr(
"pow_exp": pow_exp,
"abs_exp": abs_exp,
"squar_exp": squar_exp,
"squar_sin_exp": squar_sin_exp,
"act_exp": act_exp,
"matern52": matern52,
"matern32": matern32,
Expand Down Expand Up @@ -706,6 +714,18 @@ def _matrix_data_corr(
return r
else:
d_cont = d[:, np.logical_not(cat_features)]
if self.options["corr"] == "squar_sin_exp":
if self.options["categorical_kernel"] != MixIntKernelType.GOWER:
theta_cont_features[-len([self.design_space.is_cat_mask == True]) :] = (
np.atleast_2d(
np.array([True] * len([self.design_space.is_cat_mask == True]))
).T
)
theta_cat_features[1][
-len([self.design_space.is_cat_mask == True]) :
] = np.atleast_2d(
np.array([False] * len([self.design_space.is_cat_mask == True]))
).T

theta_cont = theta[theta_cont_features[:, 0]]
r_cont = _correlation_types[corr](theta_cont, d_cont)
Expand Down Expand Up @@ -2096,7 +2116,30 @@ def _check_param(self):
mat_dim,
)

self.options["theta0"] *= np.ones(n_param)
if self.options["corr"] == "squar_sin_exp":
if self.options["categorical_kernel"] == MixIntKernelType.GOWER:
self.options["theta0"] *= np.ones(2 * n_param)
else:
n_param += len([self.design_space.is_cat_mask == True])
self.options["theta0"] *= np.ones(n_param)

else:
self.options["theta0"] *= np.ones(n_param)
if (
not (self.options["corr"] in ["squar_exp", "abs_exp", "pow_exp"])
and not (self.is_continuous)
and not (
self.options["categorical_kernel"]
in [
MixIntKernelType.GOWER,
MixIntKernelType.COMPOUND_SYMMETRY,
MixIntKernelType.HOMO_HSPHERE,
]
)
):
raise ValueError(
"Categorical kernels should be matrix or exponential based."
)

if len(self.options["theta0"]) != d and (
self.options["categorical_kernel"]
Expand All @@ -2106,10 +2149,11 @@ def _check_param(self):
if len(self.options["theta0"]) == 1:
self.options["theta0"] *= np.ones(d)
else:
raise ValueError(
"the length of theta0 (%s) should be equal to the number of dim (%s)."
% (len(self.options["theta0"]), d)
)
if self.options["corr"] != "squar_sin_exp":
raise ValueError(
"the length of theta0 (%s) should be equal to the number of dim (%s)."
% (len(self.options["theta0"]), d)
)
if self.options["eval_noise"] or np.max(self.options["noise0"]) > 1e-12:
self.options["hyper_opt"] = "Cobyla"
warnings.warn(
Expand Down
78 changes: 78 additions & 0 deletions smt/utils/kriging.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -663,6 +663,84 @@ def pow_exp(theta, d, grad_ind=None, hess_ind=None, derivative_params=None):
return r


def squar_sin_exp(theta, d, grad_ind=None, hess_ind=None, derivative_params=None):
"""
Generative exponential autocorrelation model.
Parameters
----------
theta : list[small_d * n_comp]
Hyperparameters of the correlation model
d: np.ndarray[n_obs * (n_obs - 1) / 2, n_comp]
d_i otherwise
grad_ind : int, optional
Indice for which component the gradient dr/dtheta must be computed. The default is None.
hess_ind : int, optional
Indice for which component the hessian d²r/d²(theta) must be computed. The default is None.
derivative_paramas : dict, optional
List of arguments mandatory to compute the gradient dr/dx. The default is None.
Raises
------
Exception
Assure that theta is of the good length
Returns
-------
r: np.ndarray[n_obs * (n_obs - 1) / 2,1]
An array containing the values of the autocorrelation model.
"""

r = np.zeros((d.shape[0], 1))
n_components = d.shape[1]

# Construct/split the correlation matrix
i, nb_limit = 0, int(1e4)
while i * nb_limit <= d.shape[0]:
theta_array = theta.reshape(1, len(theta))
r[i * nb_limit : (i + 1) * nb_limit, 0] = np.exp(
-np.sum(
np.atleast_2d(theta_array[0][0 : int(len(theta) / 2)])
* np.sin(
np.atleast_2d(theta_array[0][int(len(theta) / 2) : int(len(theta))])
* d[i * nb_limit : (i + 1) * nb_limit, :]
)
** 2,
axis=1,
)
)
i += 1

# =============================================================================
# i = 0
# if grad_ind is not None:
# while i * nb_limit <= d.shape[0]:
# r[i * nb_limit : (i + 1) * nb_limit, 0] = (
# -d[i * nb_limit : (i + 1) * nb_limit, grad_ind]
# * r[i * nb_limit : (i + 1) * nb_limit, 0]
# )
# i += 1
#
# i = 0
# if hess_ind is not None:
# while i * nb_limit <= d.shape[0]:
# r[i * nb_limit : (i + 1) * nb_limit, 0] = (
# -d[i * nb_limit : (i + 1) * nb_limit, hess_ind]
# * r[i * nb_limit : (i + 1) * nb_limit, 0]
# )
# i += 1
#
# if derivative_params is not None:
# dd = derivative_params["dd"]
# r = r.T
# dr = -np.einsum("i,ij->ij", r[0], dd)
# return r.T, dr
#
# =============================================================================
return r


def matern52(theta, d, grad_ind=None, hess_ind=None, derivative_params=None):
"""
Matern 5/2 correlation model.
Expand Down

0 comments on commit 03d9ade

Please sign in to comment.