Skip to content

Commit

Permalink
fix bug for evaluating noise automatically for redundants points with…
Browse files Browse the repository at this point in the history 
… noisy output when the inputs are mixed categorical (SMTorg#509)

* cleaning X,y and standardization management in krg_based
  • Loading branch information
@Paul-Saves
Paul-Saves authored Feb 15, 2024
1 parent d978ccd commit 0c4b243
Show file tree
Hide file tree
Showing 2 changed files with 120 additions and 36 deletions.
76 changes: 76 additions & 0 deletions smt/applications/tests/test_mixed_integer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1553,6 +1553,82 @@ def test_mixed_homo_hyp_3D_PLS_cate(self):
self.assertTrue((np.abs(np.sum(np.array(sm.predict_values(xt) - yt)))) < 1e-6)
self.assertTrue((np.abs(np.sum(np.array(sm.predict_variances(xt) - 0)))) < 1e-6)

def test_compound_hetero_noise_auto(self):
xt1 = np.array([[0, 0.0], [0, 2.0], [0, 4.0]])
xt2 = np.array([[1, 0.0], [1, 2.0], [1, 3.0]])
xt3 = np.array([[2, 1.0], [2, 2.0], [2, 4.0]])

xt = np.concatenate((xt1, xt2, xt3), axis=0)
xt[:, 1] = xt[:, 1].astype(np.float64)
yt1 = np.array([0.0, 9.0, 16.0])
yt2 = np.array([0.0, -4, -13.0])
yt3 = np.array([-10, 3, 11.0])
yt = np.concatenate((yt1, yt2, yt3), axis=0)

design_space = DesignSpace(
[
CategoricalVariable(["Blue", "Red", "Green"]),
# OrdinalVariable([0,1,2]),
FloatVariable(0, 4),
]
)

# Surrogate
sm = MixedIntegerKrigingModel(
surrogate=KRG(
design_space=design_space,
categorical_kernel=MixIntKernelType.COMPOUND_SYMMETRY,
theta0=[1e-1],
corr="squar_exp",
n_start=20,
eval_noise=True,
# noise0 = [1.0,0.1,10,0.01,1.0,0.1,10,0.01,10],
use_het_noise=True,
hyper_opt="Cobyla",
),
)
sm.set_training_values(xt, yt)
sm.train()

# DOE for validation
n = 100
x_cat1 = []
x_cat2 = []
x_cat3 = []

for i in range(n):
x_cat1.append(0)
x_cat2.append(1)
x_cat3.append(2)

x_cont = np.linspace(0.0, 4.0, n)
x1 = np.concatenate(
(np.asarray(x_cat1).reshape(-1, 1), x_cont.reshape(-1, 1)), axis=1
)
x2 = np.concatenate(
(np.asarray(x_cat2).reshape(-1, 1), x_cont.reshape(-1, 1)), axis=1
)
x3 = np.concatenate(
(np.asarray(x_cat3).reshape(-1, 1), x_cont.reshape(-1, 1)), axis=1
)

y1 = sm.predict_values(x1)
y2 = sm.predict_values(x2)
y3 = sm.predict_values(x3)

# estimated variance
s2_1 = sm.predict_variances(x1)
s2_2 = sm.predict_variances(x2)
s2_3 = sm.predict_variances(x3)

self.assertEqual(
np.array_equal(
np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]),
sm._surrogate.optimal_noise,
),
True,
)

def test_mixed_homo_gaussian_3D_ord_cate(self):
xt = np.array([[0.5, 0, 5], [2, 3, 4], [5, 2, -1], [-2, 4, 0.5]])
yt = np.array([[0.0], [3], [1.0], [1.5]])
Expand Down
80 changes: 44 additions & 36 deletions smt/surrogate_models/krg_based.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -375,8 +375,42 @@ def _new_train(self):
self._corr_params = None
_, self.cat_features = compute_X_cont(self.X_train, self.design_space)
D = None # For SGP, D is not computed at all
# Center and scale X and y
(
self.X_norma,
self.y_norma,
self.X_offset,
self.y_mean,
self.X_scale,
self.y_std,
) = standardization(X.copy(), y.copy())

if not self.options["eval_noise"]:
self.optimal_noise = np.array(self.options["noise0"])
elif self.options["use_het_noise"]:
# hetGP works with unique design variables when noise variance are not given
(
self.X_norma,
index_unique,
nt_reps,
) = np.unique(self.X_norma, return_inverse=True, return_counts=True, axis=0)
self.nt = self.X_norma.shape[0]

# computing the mean of the output per unique design variable (see Binois et al., 2018)
y_norma_unique = []
for i in range(self.nt):
y_norma_unique.append(np.mean(self.y_norma[index_unique == i]))
# pointwise sensible estimates of the noise variances (see Ankenman et al., 2010)
self.optimal_noise = self.options["noise0"] * np.ones(self.nt)
for i in range(self.nt):
diff = self.y_norma[index_unique == i] - y_norma_unique[i]
if np.sum(diff**2) != 0.0:
self.optimal_noise[i] = np.std(diff, ddof=1) ** 2
self.optimal_noise = self.optimal_noise / nt_reps
self.y_norma = y_norma_unique

if not (self.is_continuous):
D, self.ij, X = gower_componentwise_distances(
D, self.ij, X_cont = gower_componentwise_distances(
X=X,
x_is_acting=is_acting,
design_space=self.design_space,
Expand All @@ -397,15 +431,15 @@ def _new_train(self):
mixint_type=MixIntKernelType.GOWER,
)
if self.options["categorical_kernel"] == MixIntKernelType.CONT_RELAX:
X2, _ = self.design_space.unfold_x(self.training_points[None][0][0])
X2, _ = self.design_space.unfold_x(X)
(
self.X2_norma,
_,
self.X2_offset,
_,
self.X2_scale,
_,
) = standardization(X2, self.training_points[None][0][1])
) = standardization(X2.copy(), y.copy())
D, _ = cross_distances(self.X2_norma)
self.Lij, self.n_levels = cross_levels(
X=self.X_train, ij=self.ij, design_space=self.design_space
Expand All @@ -422,39 +456,15 @@ def _new_train(self):
mixint_type=MixIntKernelType.CONT_RELAX,
)

# Center and scale X and y
(
self.X_norma,
self.y_norma,
self.X_offset,
self.y_mean,
self.X_scale,
self.y_std,
) = standardization(X, y)

if not self.options["eval_noise"]:
self.optimal_noise = np.array(self.options["noise0"])
elif self.options["use_het_noise"]:
# hetGP works with unique design variables when noise variance are not given
# Center and scale X_cont and y
(
self.X_norma,
index_unique,
nt_reps,
) = np.unique(self.X_norma, return_inverse=True, return_counts=True, axis=0)
self.nt = self.X_norma.shape[0]

# computing the mean of the output per unique design variable (see Binois et al., 2018)
y_norma_unique = []
for i in range(self.nt):
y_norma_unique.append(np.mean(self.y_norma[index_unique == i]))
# pointwise sensible estimates of the noise variances (see Ankenman et al., 2010)
self.optimal_noise = self.options["noise0"] * np.ones(self.nt)
for i in range(self.nt):
diff = self.y_norma[index_unique == i] - y_norma_unique[i]
if np.sum(diff**2) != 0.0:
self.optimal_noise[i] = np.std(diff, ddof=1) ** 2
self.optimal_noise = self.optimal_noise / nt_reps
self.y_norma = y_norma_unique
self.y_norma,
self.X_offset,
self.y_mean,
self.X_scale,
self.y_std,
) = standardization(X_cont.copy(), y.copy())

if self.name not in ["SGP"]:
if self.is_continuous:
Expand Down Expand Up @@ -871,7 +881,6 @@ def _reduced_likelihood_function(self, theta):
noise = tmp_var[-1]
if not (self.is_continuous):
dx = self.D

if self.options["categorical_kernel"] == MixIntKernelType.CONT_RELAX:
if "MFK" in self.name:
if (
Expand Down Expand Up @@ -920,7 +929,6 @@ def _reduced_likelihood_function(self, theta):
r = self._correlation_types[self.options["corr"]](theta, self.D).reshape(
-1, 1
)

R = np.eye(self.nt) * (1.0 + nugget + noise)
R[self.ij[:, 0], self.ij[:, 1]] = r[:, 0]
R[self.ij[:, 1], self.ij[:, 0]] = r[:, 0]
Expand Down

0 comments on commit 0c4b243

Please sign in to comment.