Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

fix bug for evaluating noise automatically for redundants points with noisy output when the inputs are mixed categorical #509

Merged
merged 4 commits into from
Feb 15, 2024
Merged

fix bug for evaluating noise automatically for redundants points with noisy output when the inputs are mixed categorical #509

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
005b6a0
fix
Paul-Saves Feb 10, 2024
e31ab31
cleaning X,y and standardization management in krg_based
Paul-Saves Feb 14, 2024
6610be2
fix
Paul-Saves Feb 14, 2024
909f3e9
fix
Paul-Saves Feb 14, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
76 changes: 76 additions & 0 deletions smt/applications/tests/test_mixed_integer.py
View file
Open in desktop
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
Open in desktop
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