Fix normalize (#217)

* add test to check normalize behavior. Must fail. * also fix for nugget * clang-format
libKriging · Sep 5, 2022 · 4bcb532 · 4bcb532
1 parent 3f09287
commit 4bcb532
Show file tree

Hide file tree

Showing 5 changed files with 162 additions and 36 deletions.
diff --git a/bindings/R/rlibkriging/tests/testthat/test-normalize.R b/bindings/R/rlibkriging/tests/testthat/test-normalize.R
@@ -0,0 +1,88 @@
+library(testthat)
+#library(rlibkriging, lib.loc="bindings/R/Rlibs")
+#rlibkriging:::optim_log(2)
+#rlibkriging:::optim_use_reparametrize(FALSE)
+#rlibkriging:::optim_set_theta_lower_factor(0.02)
+
+
+f = function(x) 1-1/2*(sin(12*x)/(1+x)+2*cos(7*x)*x^5+0.7)
+n <- 5
+set.seed(123)
+X <- as.matrix(runif(n))
+X10 = 10*X
+
+y = f(X)
+y50 = 50*y
+
+library(rlibkriging)
+
+
+context("no normalize")
+
+r_nonorm <- Kriging(y, X, "gauss", normalize=F)
+r1050_nonorm <- Kriging(y50, X10, "gauss", normalize=F)
+
+test_that(desc="theta nonorm",
+          expect_equal( r_nonorm$theta()*10 , r1050_nonorm$theta() ,tol=0.01))
+test_that(desc="beta nonorm",
+          expect_equal(r_nonorm$beta()*50 , r1050_nonorm$beta(),tol=0.01))
+test_that(desc="sigma2 nonorm",
+          expect_equal(r_nonorm$sigma2()*50*50 , r1050_nonorm$sigma2(),tol=0.01))
+
+
+test_that(desc="predict nonorm",
+          expect_equal(lapply(r_nonorm$predict(0.5),function(...){50*...}), r1050_nonorm$predict(10*0.5),tol=0.01))
+
+test_that(desc="simulate nonorm",
+          expect_equal(50*r_nonorm$simulate(1,x=0.5), r1050_nonorm$simulate(1,x=10*0.5),tol=0.01))
+
+
+r_nonorm$update(newX=0.5,newy=f(0.5))
+r1050_nonorm$update(newX=10*0.5,newy=50*f(0.5))
+
+
+test_that(desc="update theta nonorm",
+          expect_equal(r_nonorm$theta()*10 , r1050_nonorm$theta(),tol=0.01))
+test_that(desc="update beta nonorm",
+          expect_equal(r_nonorm$beta()*50 , r1050_nonorm$beta(),tol=0.01))
+test_that(desc="update sigma2 nonorm",
+          expect_equal(r_nonorm$sigma2()*50*50 , r1050_nonorm$sigma2(),tol=0.01))
+
+
+
+context("normalize")
+
+r_norm <- Kriging(y, X, "gauss", normalize=T)
+r1050_norm <- Kriging(y50, X10, "gauss", normalize=T)
+
+test_that(desc="theta norm",
+          expect_equal(r_norm$theta() , r1050_norm$theta(),tol=0.01))
+test_that(desc="beta norm",
+          expect_equal(r_norm$beta() , r1050_norm$beta(),tol=0.01))
+test_that(desc="sigma2 norm",
+          expect_equal(r_norm$sigma2() , r1050_norm$sigma2(),tol=0.01))
+
+
+test_that(desc="predict norm",
+          expect_equal(lapply(r_norm$predict(0.5),function(...){50*...}), r1050_norm$predict(10*0.5),tol=0.01))
+
+test_that(desc="simulate norm",
+          expect_equal(50*r_norm$simulate(1,x=0.5), r1050_norm$simulate(1,x=10*0.5),tol=0.01))
+
+
+plot(seq(0,1,,101),r_norm$simulate(1,seed=123,x=seq(0,1,,101)))
+points(X,y,col='red')
+plot(seq(0,10,,101),r1050_norm$simulate(1,seed=123,x=seq(0,10,,101)))
+points(X10,y50,col='red')
+
+r_norm$update(newX=0.5,newy=f(0.5))
+r1050_norm$update(newX=10*0.5,newy=50*f(0.5))
+
+test_that(desc="update theta norm",
+          expect_equal(r_norm$theta() , r1050_norm$theta(),tol=0.01))
+test_that(desc="update beta norm",
+          expect_equal(r_norm$beta() , r1050_norm$beta(),tol=0.01))
+test_that(desc="update sigma2 norm",
+          expect_equal(r_norm$sigma2() , r1050_norm$sigma2(),tol=0.01))
+
+
diff --git a/src/lib/Kriging.cpp b/src/lib/Kriging.cpp
@@ -968,7 +968,8 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
   double centerY;
   double scaleY;
   // Normalization of inputs and output
-  if (normalize) {
+  m_normalize = normalize;
+  if (m_normalize) {
     centerX = min(X, 0);
     scaleX = max(X, 0) - min(X, 0);
     centerY = min(y);
@@ -1052,8 +1053,8 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
     m_est_sigma2 = parameters.is_sigma2_estim;
 
   } else {
-    arma::vec theta_lower = Optim::theta_lower_factor * trans(max(X, 0) - min(X, 0));
-    arma::vec theta_upper = Optim::theta_upper_factor * trans(max(X, 0) - min(X, 0));
+    arma::vec theta_lower = Optim::theta_lower_factor * trans(max(m_X, 0) - min(m_X, 0));
+    arma::vec theta_upper = Optim::theta_upper_factor * trans(max(m_X, 0) - min(m_X, 0));
 
     if (Optim::variogram_bounds_heuristic) {
       arma::vec dy2 = arma::zeros(n * n);
@@ -1135,6 +1136,7 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
           arma::cout << "  null gradient tolerance: " << Optim::gradient_tolerance << arma::endl;
           arma::cout << "  constant objective tolerance: " << Optim::objective_rel_tolerance << arma::endl;
           arma::cout << "  reparametrize: " << Optim::reparametrize << arma::endl;
+          arma::cout << "  normalize: " << m_normalize << arma::endl;
           arma::cout << "  lower_bounds: " << theta_lower.t() << " ";
           arma::cout << "  upper_bounds: " << theta_upper.t() << " ";
           arma::cout << "  start_point: " << theta0.row(i).t() << " ";
@@ -1259,6 +1261,7 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
           arma::cout << "  null gradient tolerance: " << Optim::gradient_tolerance << arma::endl;
           arma::cout << "  constant objective tolerance: " << Optim::objective_rel_tolerance << arma::endl;
           arma::cout << "  reparametrize: " << Optim::reparametrize << arma::endl;
+          arma::cout << "  normalize: " << m_normalize << arma::endl;
           arma::cout << "  lower_bounds: " << theta_lower.t() << " ";
           arma::cout << "  upper_bounds: " << theta_upper.t() << " ";
           arma::cout << "  start_point: " << theta0.row(i).t() << " ";
@@ -1308,12 +1311,9 @@ LIBKRIGING_EXPORT void Kriging::fit(const arma::colvec& y,
         }
       }
     } else
-      throw std::runtime_error("Unsupported optim: " + optim + " (supported are: none, BFGS*, Newton*)");
+      throw std::runtime_error("Unsupported optim: " + optim + " (supported are: none, BFGS[#], Newton[#])");
   }
 
-  if (!parameters.sigma2.has_value())
-    m_sigma2 *= scaleY * scaleY;
-
   // arma::cout << "theta:" << m_theta << arma::endl;
 }
 
@@ -1405,6 +1405,8 @@ Kriging::predict(const arma::mat& Xp, bool withStd, bool withCov, bool withDeriv
     arma::mat s2_predict = total_sd2 * (1.0 - s2_predict_1 + s2_predict_2);
     s2_predict.transform([](double val) { return (std::isnan(val) || val < 0 ? 0.0 : val); });
     pred_stdev = sqrt(s2_predict);
+
+    pred_stdev *= m_scaleY;
   }
 
   if (withCov) {
@@ -1417,6 +1419,8 @@ Kriging::predict(const arma::mat& Xp, bool withStd, bool withCov, bool withDeriv
     }
 
     pred_cov = total_sd2 * (R_predpred - trans(Tinv_pred) * Tinv_pred + trans(s2_predict_mat) * s2_predict_mat);
+
+    pred_cov *= m_scaleY;
   }
 
   if (withDeriv) {
@@ -1593,17 +1597,22 @@ LIBKRIGING_EXPORT void Kriging::update(const arma::vec& newy, const arma::mat& n
                              + std::to_string(newX.n_cols) + "), y: (" + std::to_string(newy.n_elem) + ")");
 
   // rebuild starting parameters
-  Parameters parameters{std::make_optional(this->m_sigma2),
+  Parameters parameters{std::make_optional(this->m_sigma2 * this->m_scaleY * this->m_scaleY),
                         this->m_est_sigma2,
-                        std::make_optional(trans(this->m_theta)),
+                        std::make_optional(trans(this->m_theta) % this->m_scaleX),
                         this->m_est_theta,
-                        std::make_optional(trans(this->m_beta)),
+                        std::make_optional(trans(this->m_beta) * this->m_scaleY),
                         this->m_est_beta};
   // re-fit
   // TODO refit() method which will use Shurr forms to fast update matrix (R, ...)
-  bool normalize = (m_centerY == 0.0 && m_scaleY == 1) && (m_centerX.is_zero() && m_scaleX.is_zero());
-  this->fit(
-      arma::join_cols(m_y, newy), arma::join_cols(m_X, newX), m_regmodel, normalize, m_optim, m_objective, parameters);
+  this->fit(arma::join_cols(m_y * this->m_scaleY + this->m_centerY,
+                            newy),  // de-normalize previous data according to suite unnormed new data
+            arma::join_cols((m_X.each_row() % this->m_scaleX).each_row() + this->m_centerX, newX),
+            m_regmodel,
+            m_normalize,
+            m_optim,
+            m_objective,
+            parameters);
 }
 
 LIBKRIGING_EXPORT std::string Kriging::summary() const {
@@ -1616,7 +1625,16 @@ LIBKRIGING_EXPORT std::string Kriging::summary() const {
     });
   };
 
-  oss << "* data: " << m_X.n_rows << " x " << m_X.n_cols << " -> " << m_y.n_rows << " x " << m_y.n_cols << "\n";
+  oss << "* data";
+  oss << ((m_normalize) ? " (normalized): " : ": ");
+  arma::rowvec Xmins = arma::min(m_X, 0);
+  arma::rowvec Xmaxs = arma::max(m_X, 0);
+  for (arma::uword i = 0; i < m_X.n_cols; i++) {
+    oss << "[" << Xmins[i] << "," << Xmaxs[i] << "]";
+    if (i < m_X.n_cols - 1)
+      oss << "x";
+  }
+  oss << " -> [" << arma::min(m_y) << "," << arma::max(m_y) << "]\n";
   oss << "* trend " << Trend::toString(m_regmodel);
   oss << ((m_est_beta) ? " (est.): " : ": ");
   colvec_printer(m_beta);

diff --git a/src/lib/NuggetKriging.cpp b/src/lib/NuggetKriging.cpp
@@ -573,7 +573,8 @@ LIBKRIGING_EXPORT void NuggetKriging::fit(const arma::colvec& y,
   double centerY;
   double scaleY;
   // Normalization of inputs and output
-  if (normalize) {
+  m_normalize = normalize;
+  if (m_normalize) {
     centerX = min(X, 0);
     scaleX = max(X, 0) - min(X, 0);
     centerY = min(y);
@@ -679,8 +680,8 @@ LIBKRIGING_EXPORT void NuggetKriging::fit(const arma::colvec& y,
   } else if (optim.rfind("BFGS", 0) == 0) {
     Random::init();
 
-    arma::vec theta_lower = Optim::theta_lower_factor * trans(max(X, 0) - min(X, 0));
-    arma::vec theta_upper = Optim::theta_upper_factor * trans(max(X, 0) - min(X, 0));
+    arma::vec theta_lower = Optim::theta_lower_factor * trans(max(m_X, 0) - min(m_X, 0));
+    arma::vec theta_upper = Optim::theta_upper_factor * trans(max(m_X, 0) - min(m_X, 0));
 
     if (Optim::variogram_bounds_heuristic) {
       arma::vec dy2 = arma::zeros(n * n);
@@ -778,6 +779,7 @@ LIBKRIGING_EXPORT void NuggetKriging::fit(const arma::colvec& y,
         arma::cout << "  null gradient tolerance: " << Optim::gradient_tolerance << arma::endl;
         arma::cout << "  constant objective tolerance: " << Optim::objective_rel_tolerance << arma::endl;
         arma::cout << "  reparametrize: " << Optim::reparametrize << arma::endl;
+        arma::cout << "  normalize: " << m_normalize << arma::endl;
         arma::cout << "  lower_bounds: " << theta_lower.t() << "";
         arma::cout << "                " << alpha_lower << arma::endl;
         arma::cout << "  upper_bounds: " << theta_upper.t() << "";
@@ -881,12 +883,7 @@ LIBKRIGING_EXPORT void NuggetKriging::fit(const arma::colvec& y,
       }
     }
   } else
-    throw std::runtime_error("Unsupported optim: " + optim + " (supported are: none, BFGS**)");
-
-  if (!parameters.sigma2.has_value())
-    m_sigma2 *= scaleY * scaleY;
-  if (!parameters.nugget.has_value())
-    m_nugget *= scaleY * scaleY;
+    throw std::runtime_error("Unsupported optim: " + optim + " (supported are: none, BFGS[#])");
 
   // arma::cout << "theta:" << m_theta << arma::endl;
 }
@@ -968,6 +965,8 @@ NuggetKriging::predict(const arma::mat& Xp, bool withStd, bool withCov, bool wit
     arma::mat s2_predict = total_sd2 * (1.0 - s2_predict_1 + s2_predict_2);
     s2_predict.transform([](double val) { return (std::isnan(val) || val < 0 ? 0.0 : val); });
     pred_stdev = sqrt(s2_predict);
+
+    pred_stdev *= m_scaleY;
   }
 
   if (withCov) {
@@ -985,6 +984,8 @@ NuggetKriging::predict(const arma::mat& Xp, bool withStd, bool withCov, bool wit
     // cond.cov <- cond.cov + crossprod(s2.predict.mat)
 
     pred_cov = total_sd2 * (R_predpred - trans(Tinv_pred) * Tinv_pred + trans(s2_predict_mat) * s2_predict_mat);
+
+    pred_cov *= m_scaleY;
   }
 
   if (withDeriv) {
@@ -1138,19 +1139,25 @@ LIBKRIGING_EXPORT void NuggetKriging::update(const arma::vec& newy, const arma::
                              + std::to_string(newX.n_cols) + "), y: (" + std::to_string(newy.n_elem) + ")");
 
   // rebuild starting parameters
-  Parameters parameters{std::make_optional(arma::vec(1, arma::fill::value(this->m_nugget))),
-                        this->m_est_nugget,
-                        std::make_optional(arma::vec(1, arma::fill::value(this->m_sigma2))),
-                        this->m_est_sigma2,
-                        std::make_optional(trans(this->m_theta)),
-                        this->m_est_theta,
-                        std::make_optional(trans(this->m_beta)),
-                        this->m_est_beta};
+  Parameters parameters{
+      std::make_optional(arma::vec(1, arma::fill::value(this->m_nugget * this->m_scaleY * this->m_scaleY))),
+      this->m_est_nugget,
+      std::make_optional(arma::vec(1, arma::fill::value(this->m_sigma2 * this->m_scaleY * this->m_scaleY))),
+      this->m_est_sigma2,
+      std::make_optional(trans(this->m_theta) % this->m_scaleX),
+      this->m_est_theta,
+      std::make_optional(trans(this->m_beta) * this->m_scaleY),
+      this->m_est_beta};
   // re-fit
   // TODO refit() method which will use Shurr forms to fast update matrix (R, ...)
-  bool normalize = (m_centerY == 0.0 && m_scaleY == 1) && (m_centerX.is_zero() && m_scaleX.is_zero());
-  this->fit(
-      arma::join_cols(m_y, newy), arma::join_cols(m_X, newX), m_regmodel, normalize, m_optim, m_objective, parameters);
+  this->fit(arma::join_cols(m_y * this->m_scaleY + this->m_centerY,
+                            newy),  // de-normalize previous data according to suite unnormed new data
+            arma::join_cols((m_X.each_row() % this->m_scaleX).each_row() + this->m_centerX, newX),
+            m_regmodel,
+            m_normalize,
+            m_optim,
+            m_objective,
+            parameters);
 }
 
 LIBKRIGING_EXPORT std::string NuggetKriging::summary() const {
@@ -1163,7 +1170,16 @@ LIBKRIGING_EXPORT std::string NuggetKriging::summary() const {
     });
   };
 
-  oss << "* data: " << m_X.n_rows << " x " << m_X.n_cols << " -> " << m_y.n_rows << " x " << m_y.n_cols << "\n";
+  oss << "* data";
+  oss << ((m_normalize) ? " (normalized): " : ": ");
+  arma::rowvec Xmins = arma::min(m_X, 0);
+  arma::rowvec Xmaxs = arma::max(m_X, 0);
+  for (arma::uword i = 0; i < m_X.n_cols; i++) {
+    oss << "[" << Xmins[i] << "," << Xmaxs[i] << "]";
+    if (i < m_X.n_cols - 1)
+      oss << "x";
+  }
+  oss << " -> [" << arma::min(m_y) << "," << arma::max(m_y) << "]\n";
   oss << "* trend " << Trend::toString(m_regmodel);
   oss << ((m_est_beta) ? " (est.): " : ": ");
   colvec_printer(m_beta);
@@ -1186,4 +1202,4 @@ LIBKRIGING_EXPORT std::string NuggetKriging::summary() const {
   oss << "    * objective: " << m_objective << "\n";
   oss << "    * optim: " << m_optim << "\n";
   return oss.str();
-}
+}
diff --git a/src/lib/include/libKriging/Kriging.hpp b/src/lib/include/libKriging/Kriging.hpp
@@ -38,6 +38,7 @@ class Kriging {
   [[nodiscard]] const arma::colvec& y() const { return m_y; };
   [[nodiscard]] const double& centerY() const { return m_centerY; };
   [[nodiscard]] const double& scaleY() const { return m_scaleY; };
+  [[nodiscard]] const bool& normalize() const { return m_normalize; };
   [[nodiscard]] const Trend::RegressionModel& regmodel() const { return m_regmodel; };
   [[nodiscard]] const arma::mat& F() const { return m_F; };
   [[nodiscard]] const arma::mat& T() const { return m_T; };
@@ -58,6 +59,7 @@ class Kriging {
   arma::colvec m_y;
   double m_centerY;
   double m_scaleY;
+  bool m_normalize;
   Trend::RegressionModel m_regmodel;
   std::string m_optim;
   std::string m_objective;

diff --git a/src/lib/include/libKriging/NuggetKriging.hpp b/src/lib/include/libKriging/NuggetKriging.hpp
@@ -40,6 +40,7 @@ class NuggetKriging {
   [[nodiscard]] const arma::colvec& y() const { return m_y; };
   [[nodiscard]] const double& centerY() const { return m_centerY; };
   [[nodiscard]] const double& scaleY() const { return m_scaleY; };
+  [[nodiscard]] const bool& normalize() const { return m_normalize; };
   [[nodiscard]] const Trend::RegressionModel& regmodel() const { return m_regmodel; };
   [[nodiscard]] const arma::mat& F() const { return m_F; };
   [[nodiscard]] const arma::mat& T() const { return m_T; };
@@ -62,6 +63,7 @@ class NuggetKriging {
   arma::colvec m_y;
   double m_centerY;
   double m_scaleY;
+  bool m_normalize;
   Trend::RegressionModel m_regmodel;
   std::string m_optim;
   std::string m_objective;