rm quadfailover (#276)

libKriging · Dec 3, 2023 · 372ce0f · 372ce0f
1 parent 927f576
commit 372ce0f
Show file tree

Hide file tree

Showing 6 changed files with 2 additions and 96 deletions.
diff --git a/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogLik.R b/bindings/R/rlibkriging/tests/testthat/bench-KrigingLogLik.R
@@ -45,9 +45,9 @@ rlibkriging:::covariance_use_approx_singular(TRUE)
 
 plot(xx,Vectorize(function(x)DiceKriging::logLikFun(rep(x,3),k))(xx))
 points(xx,Vectorize(function(x)logLikelihoodFun(r,rep(x,3))$logLikelihood)(xx),col='red')
-rlibkriging:::optim_use_quadfailover(FALSE)
+rlibkriging:::covariance_use_approx_singular(FALSE)
 points(xx,Vectorize(function(x)logLikelihoodFun(r_nf,rep(x,3))$logLikelihood)(xx),col='orange')
-rlibkriging:::optim_use_quadfailover(TRUE)
+rlibkriging:::covariance_use_approx_singular(TRUE)
 
 plot(xx,log(times$R_ll),ylim=c(log(min(min(times$R_ll,na.rm = T),min(times$cpp_ll,na.rm = T))),log(max(max(times$R_ll,na.rm = T),max(times$cpp_ll,na.rm = T)))),xlab="x",ylab="log(user_time (s))", panel.first=grid())
 text(20,-1,"DiceKriging")

diff --git a/src/lib/Kriging.cpp b/src/lib/Kriging.cpp
@@ -131,21 +131,6 @@ double Kriging::_logLikelihood(const arma::vec& _theta,
   }
   t0 = Bench::toc(bench, "R = Cov(dX)", t0);
 
-  // Sly turnaround for too long range: use shorter range penalized, and force gradient to point at shorter range
-  // (assuming a Newton like method for wrapping optim)
-  if (Optim::quadfailover)
-    if (arma::rcond(R) < R.n_rows * arma::datum::eps) {
-      // throw std::runtime_error("Covariance matrix is singular");
-      // Try use midpoint of theta and
-      // arma::cout << "Covariance matrix is singular, try use midpoint of theta" << std::endl;
-      double ll_2 = _logLikelihood(_theta / 2, grad_out, hess_out, okm_data, bench);
-      if (grad_out)
-        *grad_out *= 2;
-      if (hess_out)
-        *hess_out *= 4;
-      return ll_2 - log(2);  // emulates likelihood/2
-    }
-
   // Cholesky decompostion of covariance matrix
   fd->T = LinearAlgebra::safe_chol_lower(R);  // Do NOT trimatl T (slower because copy): trimatl(chol(R, "lower"));
   t0 = Bench::toc(bench, "T = Chol(R)", t0);
@@ -455,19 +440,6 @@ double Kriging::_leaveOneOut(const arma::vec& _theta,
   }
   t0 = Bench::toc(bench, "R = Cov(dX)", t0);
 
-  // Sly turnaround for too long range: use shorter range penalized, and force gradient to point at shorter range
-  // (assuming a Newton like method for wrapping optim)
-  if (Optim::quadfailover)
-    if (arma::rcond(R) < R.n_rows * arma::datum::eps) {
-      // throw std::runtime_error("Covariance matrix is singular");
-      // Try use midpoint of theta and
-      // arma::cout << "Covariance matrix is singular, try use midpoint of theta" << std::endl;
-      double loo_2 = _leaveOneOut(_theta / 2, grad_out, yhat_out, okm_data, bench);
-      if (grad_out)
-        *grad_out *= 2;
-      return loo_2 + log(2);  // emulates likelihood/2
-    }
-
   // Cholesky decompostion of covariance matrix
   fd->T = LinearAlgebra::safe_chol_lower(R);
   t0 = Bench::toc(bench, "T = Chol(R)", t0);
@@ -717,19 +689,6 @@ double Kriging::_logMargPost(const arma::vec& _theta,
   }
   t0 = Bench::toc(bench, "R = Cov(dX)", t0);
 
-  // Sly turnaround for too long range: use shorter range penalized, and force gradient to point at shorter range
-  // (assuming a Newton like method for wrapping optim)
-  if (Optim::quadfailover)
-    if (arma::rcond(R) < R.n_rows * arma::datum::eps) {
-      // throw std::runtime_error("Covariance matrix is singular");
-      // Try use midpoint of theta and
-      // arma::cout << "Covariance matrix is singular, try use midpoint of theta" << std::endl;
-      double lmp_2 = _logMargPost(_theta / 2, grad_out, okm_data, bench);
-      if (grad_out)
-        *grad_out *= 2;
-      return lmp_2 - log(2);  // emulates likelihood/2
-    }
-
   // Cholesky decompostion of covariance matrix
   fd->T = LinearAlgebra::safe_chol_lower(R);
   t0 = Bench::toc(bench, "T = Chol(R)", t0);

diff --git a/src/lib/NoiseKriging.cpp b/src/lib/NoiseKriging.cpp
@@ -140,19 +140,6 @@ double NoiseKriging::_logLikelihood(const arma::vec& _theta_sigma2,
   }
   t0 = Bench::toc(bench, "R = Cov(dX)", t0);
 
-  // Sly turnaround for too long range: use shorter range penalized, and force gradient to point at shorter range
-  // (assuming a Newton like method for wrapping optim)
-  if (Optim::quadfailover)
-    if (arma::rcond(R) < R.n_rows * arma::datum::eps) {
-      // throw std::runtime_error("Covariance matrix is singular");
-      // Try use midpoint of theta and
-      // arma::cout << "Covariance matrix is singular, try use midpoint of theta" << std::endl;
-      double ll_2 = _logLikelihood(_theta_sigma2 / 2, grad_out, okm_data, bench);
-      if (grad_out)
-        *grad_out *= 2;
-      return ll_2 - log(2);  // emulates likelihood/2
-    }
-
   // Cholesky decompostion of covariance matrix
   fd->T = LinearAlgebra::safe_chol_lower(R);  // Do NOT trimatl T (slower because copy): trimatl(chol(R, "lower"));
   t0 = Bench::toc(bench, "T = Chol(R)", t0);

diff --git a/src/lib/NuggetKriging.cpp b/src/lib/NuggetKriging.cpp
@@ -136,19 +136,6 @@ double NuggetKriging::_logLikelihood(const arma::vec& _theta_alpha,
   }
   t0 = Bench::toc(bench, "R = Cov(dX)", t0);
 
-  // Sly turnaround for too long range: use shorter range penalized, and force gradient to point at shorter range
-  // (assuming a Newton like method for wrapping optim)
-  if (Optim::quadfailover)
-    if (arma::rcond(R) < R.n_rows * arma::datum::eps) {
-      // throw std::runtime_error("Covariance matrix is singular");
-      // Try use midpoint of theta and
-      // arma::cout << "Covariance matrix is singular, try use midpoint of theta" << std::endl;
-      double ll_2 = _logLikelihood(_theta_alpha / 2, grad_out, okm_data, bench);
-      if (grad_out)
-        *grad_out *= 2;
-      return ll_2 - log(2);  // emulates likelihood/2
-    }
-
   // Cholesky decompostion of covariance matrix
   fd->T = LinearAlgebra::safe_chol_lower(R);  // Do NOT trimatl T (slower because copy): trimatl(chol(R, "lower"));
   t0 = Bench::toc(bench, "T = Chol(R)", t0);
@@ -391,19 +378,6 @@ double NuggetKriging::_logMargPost(const arma::vec& _theta_alpha,
   }
   t0 = Bench::toc(bench, "R = Cov(dX)", t0);
 
-  // Sly turnaround for too long range: use shorter range penalized, and force gradient to point at shorter range
-  // (assuming a Newton like method for wrapping optim)
-  if (Optim::quadfailover)
-    if (arma::rcond(R) < R.n_rows * arma::datum::eps) {
-      // throw std::runtime_error("Covariance matrix is singular");
-      // Try use midpoint of theta and
-      // arma::cout << "Covariance matrix is singular, try use midpoint of theta" << std::endl;
-      double lmp_2 = _logMargPost(_theta_alpha / 2, grad_out, okm_data, bench);
-      if (grad_out)
-        *grad_out *= 2;
-      return lmp_2 - log(2);  // emulates likelihood/2
-    }
-
   // Cholesky decompostion of covariance matrix
   fd->T = LinearAlgebra::safe_chol_lower(R);
   t0 = Bench::toc(bench, "T = Chol(R)", t0);

diff --git a/src/lib/Optim.cpp b/src/lib/Optim.cpp
@@ -98,13 +98,3 @@ LIBKRIGING_EXPORT void Optim::set_objective_rel_tolerance(double objective_rel_t
 LIBKRIGING_EXPORT double Optim::get_objective_rel_tolerance() {
   return Optim::objective_rel_tolerance;
 };
-
-bool Optim::quadfailover = true;
-
-LIBKRIGING_EXPORT void Optim::use_quadfailover(bool do_quadfailover) {
-  Optim::quadfailover = do_quadfailover;
-};
-
-LIBKRIGING_EXPORT bool Optim::is_quadfailover() {
-  return Optim::quadfailover;
-};
diff --git a/src/lib/include/libKriging/Optim.hpp b/src/lib/include/libKriging/Optim.hpp
@@ -46,10 +46,6 @@ class Optim {
   static double objective_rel_tolerance;
   LIBKRIGING_EXPORT static void set_objective_rel_tolerance(double objective_rel_tolerance_val);
   LIBKRIGING_EXPORT static double get_objective_rel_tolerance();
-
-  static bool quadfailover;
-  LIBKRIGING_EXPORT static void use_quadfailover(bool do_quadfailover);
-  LIBKRIGING_EXPORT static bool is_quadfailover();
 };
 
 #endif  // LIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_LINLIBKRIGING_SRC_LIB_INCLUDE_LIBKRIGING_OPTIM_HPPEARALGEBRA_HPP