From 1103982c39cc7af0fd1f9180a67aef385f9bf244 Mon Sep 17 00:00:00 2001 From: Saves Paul Date: Thu, 19 Sep 2024 14:51:25 +0200 Subject: [PATCH] Fix is_ri computational time (#652) * fix is_ri * remove c noisy variance * fix test is ri --- smt/applications/ego.py | 2 + smt/surrogate_models/krg_based.py | 150 +++++++++++++------ smt/surrogate_models/tests/test_krg_noise.py | 2 +- 3 files changed, 111 insertions(+), 43 deletions(-) diff --git a/smt/applications/ego.py b/smt/applications/ego.py index 9fb2a04ad..b00966d22 100644 --- a/smt/applications/ego.py +++ b/smt/applications/ego.py @@ -274,6 +274,8 @@ def _setup_optimizer(self, fun): """ # Set the model self.gpr = self.options["surrogate"] + if "is_ri" in self.gpr.options: + self.gpr.options["is_ri"] = self.options["is_ri"] self.design_space: BaseDesignSpace = self.gpr.design_space if isinstance(self.design_space, DesignSpace): self.design_space.seed = self.options["random_state"] diff --git a/smt/surrogate_models/krg_based.py b/smt/surrogate_models/krg_based.py index 93fc0dc1d..87556538d 100644 --- a/smt/surrogate_models/krg_based.py +++ b/smt/surrogate_models/krg_based.py @@ -196,6 +196,9 @@ def _initialize(self): desc="definition of the (hierarchical) design space: " "use `smt.utils.design_space.DesignSpace` as the main API. Also accepts list of float variable bounds", ) + declare( + "is_ri", False, types=bool, desc="activate reinterpolation for noisy cases" + ) self.options.declare( "random_state", default=41, @@ -913,8 +916,10 @@ def _reduced_likelihood_function(self, theta): nugget = self.options["nugget"] # Nugget to ensure that the Cholesky decomposition can be performed if self.options["eval_noise"]: - nugget = 100.0 * np.finfo(np.double).eps - + if self.options["is_ri"]: + nugget = 100.0 * np.finfo(np.double).eps + else: + nugget = 0 noise = self.noise0 tmp_var = theta if self.options["use_het_noise"]: @@ -988,10 +993,14 @@ def _reduced_likelihood_function(self, theta): ) return reduced_likelihood_function_value, par - R_noisy = np.eye(self.nt) * (1.0 + nugget + noise) - R_noisy[self.ij[:, 0], self.ij[:, 1]] = r[:, 0] - R_noisy[self.ij[:, 1], self.ij[:, 0]] = r[:, 0] - R = np.eye(self.nt) * (1.0 + nugget) + if self.options["is_ri"]: + R_noisy = np.eye(self.nt) * (1.0 + nugget + noise) + R_noisy[self.ij[:, 0], self.ij[:, 1]] = r[:, 0] + R_noisy[self.ij[:, 1], self.ij[:, 0]] = r[:, 0] + R = np.eye(self.nt) * (1.0 + nugget) + else: + 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] @@ -1010,45 +1019,105 @@ def _reduced_likelihood_function(self, theta): print(np.linalg.eig(R)[0]) return reduced_likelihood_function_value, par - # Computation of R_ri for the reinterpolation case - C_inv = np.linalg.inv(C) - R_inv = np.dot(C_inv.T, C_inv) - R_ri = R_noisy @ R_inv @ R_noisy - par["C"] = C + if self.options["is_ri"]: + # Computation of R_ri for the reinterpolation case + C_inv = np.linalg.inv(C) + R_inv = np.dot(C_inv.T, C_inv) + R_ri = R_noisy @ R_inv @ R_noisy + par["C"] = C + + par["sigma2"] = None + par["sigma2_ri"] = None + _, _, sigma2_ri = self._compute_sigma2( + R_ri, reduced_likelihood_function_value, par, p, q, is_ri=True + ) + if sigma2_ri is not None: + par["sigma2_ri"] = sigma2_ri * self.y_std**2.0 - par["sigma2"] = None - par["sigma2_ri"] = None - _, _, sigma2_ri = self._compute_sigma2( - R_ri, reduced_likelihood_function_value, par, p, q, is_ri=True - ) - if sigma2_ri is not None: - par["sigma2_ri"] = sigma2_ri * self.y_std**2.0 + reduced_likelihood_function_value, par, sigma2 = self._compute_sigma2( + R_noisy, reduced_likelihood_function_value, par, p, q, is_ri=False + ) + if sigma2 is not None: + par["sigma2"] = sigma2 * self.y_std**2.0 - reduced_likelihood_function_value, par, sigma2 = self._compute_sigma2( - R_noisy, reduced_likelihood_function_value, par, p, q, is_ri=False - ) - if sigma2 is not None: - par["sigma2"] = sigma2 * self.y_std**2.0 + if self.name in ["MGP"]: + reduced_likelihood_function_value += self._reduced_log_prior(theta) - if self.name in ["MGP"]: - reduced_likelihood_function_value += self._reduced_log_prior(theta) + # A particular case when f_min_cobyla fail + if (self.best_iteration_fail is not None) and ( + not np.isinf(reduced_likelihood_function_value) + ): + if reduced_likelihood_function_value > self.best_iteration_fail: + self.best_iteration_fail = reduced_likelihood_function_value + self._thetaMemory = np.array(tmp_var) - # A particular case when f_min_cobyla fail - if (self.best_iteration_fail is not None) and ( - not np.isinf(reduced_likelihood_function_value) - ): - if reduced_likelihood_function_value > self.best_iteration_fail: + elif (self.best_iteration_fail is None) and ( + not np.isinf(reduced_likelihood_function_value) + ): self.best_iteration_fail = reduced_likelihood_function_value self._thetaMemory = np.array(tmp_var) + if reduced_likelihood_function_value > 1e15: + reduced_likelihood_function_value = 1e15 + return reduced_likelihood_function_value, par + else: + # Get generalized least squared solution + Ft = linalg.solve_triangular(C, self.F, lower=True) + Q, G = linalg.qr(Ft, mode="economic") + sv = linalg.svd(G, compute_uv=False) + rcondG = sv[-1] / sv[0] + if rcondG < 1e-10: + # Check F + sv = linalg.svd(self.F, compute_uv=False) + condF = sv[0] / sv[-1] + if condF > 1e15: + raise Exception( + "F is too ill conditioned. Poor combination " + "of regression model and observations." + ) - elif (self.best_iteration_fail is None) and ( - not np.isinf(reduced_likelihood_function_value) - ): - self.best_iteration_fail = reduced_likelihood_function_value - self._thetaMemory = np.array(tmp_var) - if reduced_likelihood_function_value > 1e15: - reduced_likelihood_function_value = 1e15 - return reduced_likelihood_function_value, par + else: + # Ft is too ill conditioned, get out (try different theta) + return reduced_likelihood_function_value, par + + Yt = linalg.solve_triangular(C, self.y_norma, lower=True) + beta = linalg.solve_triangular(G, np.dot(Q.T, Yt)) + rho = Yt - np.dot(Ft, beta) + + # The determinant of R is equal to the squared product of the diagonal + # elements of its Cholesky decomposition C + detR = (np.diag(C) ** (2.0 / self.nt)).prod() + + sigma2 = (rho**2.0).sum(axis=0) / (self.nt - p - q) + reduced_likelihood_function_value = -(self.nt - p - q) * np.log10( + sigma2.sum() + ) - self.nt * np.log10(detR) + par["sigma2"] = sigma2 * self.y_std**2.0 + par["beta"] = beta + par["gamma"] = linalg.solve_triangular(C.T, rho) + par["C"] = C + par["Ft"] = Ft + par["G"] = G + par["Q"] = Q + + if self.name in ["MGP"]: + reduced_likelihood_function_value += self._reduced_log_prior(theta) + + # A particular case when f_min_cobyla fail + if (self.best_iteration_fail is not None) and ( + not np.isinf(reduced_likelihood_function_value) + ): + if reduced_likelihood_function_value > self.best_iteration_fail: + self.best_iteration_fail = reduced_likelihood_function_value + self._thetaMemory = np.array(tmp_var) + + elif (self.best_iteration_fail is None) and ( + not np.isinf(reduced_likelihood_function_value) + ): + self.best_iteration_fail = reduced_likelihood_function_value + self._thetaMemory = np.array(tmp_var) + if reduced_likelihood_function_value > 1e15: + reduced_likelihood_function_value = 1e15 + return reduced_likelihood_function_value, par def _compute_sigma2( self, R, reduced_likelihood_function_value, par, p, q, is_ri=False @@ -1726,10 +1795,7 @@ def _predict_variances( # Compute the correlation function r = self.corr(d).reshape(n_eval, self.nt) X_cont = (X_cont - self.X_offset) / self.X_scale - if is_ri: - C = self.optimal_par["C"] - else: - C = self.optimal_par["C_noisy"] + C = self.optimal_par["C"] rt = linalg.solve_triangular(C, r.T, lower=True) u = linalg.solve_triangular( diff --git a/smt/surrogate_models/tests/test_krg_noise.py b/smt/surrogate_models/tests/test_krg_noise.py index 304a9be47..aba8757d3 100644 --- a/smt/surrogate_models/tests/test_krg_noise.py +++ b/smt/surrogate_models/tests/test_krg_noise.py @@ -86,7 +86,7 @@ def test_null_predict_variance_ri_at_data_points(self): yt = np.array([0.0, 1.0, 1.5, 1.1, 1.0]) # defining the models - sm_noisy = KRG(noise0=[1e-1], print_global=False) + sm_noisy = KRG(noise0=[1e-1], print_global=False, is_ri=True) # training the models sm_noisy.set_training_values(xt, yt) sm_noisy.train()