Repository: libKriging/libKriging (language: python)

Issue:
Error in LMP/linear trend
```
# 1D function, small design
X <- as.matrix(c(0.0, 0.25, 0.5, 0.75, 1.0))
f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
y <- f(X) # + 10*rnorm(nrow(X))

library(rlibkriging)
k <- Kriging(y, X, kernel="matern5_2", objective = "LMP", optim = "BFGS", regmodel = "linear")
plot( function(t) k$logMargPostFun(t)$logMargPost ,xlab=expression(theta))
abline(v=k$theta(), col="red")

# LL
obj_fun = function(theta,...) logLikelihoodFun(k,theta)$logLikelihood
obj_grad = function(theta,...) logLikelihoodFun(k,theta,grad = T)$logLikelihoodGrad
obj_grad_fd = function(theta,...) (obj_fun(theta+0.01)-obj_fun(theta-0.01))/0.02

plot(Vectorize(obj_fun),col='blue')
for (x in seq(0.01,1,,11)){
  fx = obj_fun(x)
  gx = obj_grad(x)
  gx_ref = obj_grad_fd(x)
  arrows(x,fx,x+.1,fx+.1*gx,col='red')
  arrows(x,fx,x+.11,fx+.11*gx_ref,col='blue')
}

# LMP
obj_fun = function(theta,...) logMargPostFun(k,theta)$logMargPost
obj_grad = function(theta,...) logMargPostFun(k,theta,grad = T)$logMargPostGrad
obj_grad_fd = function(theta,...) (obj_fun(theta+0.01)-obj_fun(theta-0.01))/0.02

plot(Vectorize(obj_fun),col='blue')
for (x in seq(0.01,1,,11)){
  fx = obj_fun(x)
  gx = obj_grad(x)
  gx_ref = obj_grad_fd(x)
  arrows(x,fx,x+.1,fx+.1*gx,col='red')
  arrows(x,fx,x+.11,fx+.11*gx_ref,col='blue')
}
```


Discussion / hints:
derivative of lmp / trend linear is ok indeed.
... but the lmp itself is not

```
# 1D function, small design
X <- as.matrix(c(0.0, 0.25, 0.5, 0.75, 1.0))
f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
y <- f(X) # + 10*rnorm(nrow(X))

library(rlibkriging)
k <- Kriging(y, X, kernel="matern5_2", objective = "LMP", optim = "BFGS", regmodel = "linear")
plot( function(t) k$logMargPostFun(t)$logMargPost ,xlab=expression(theta))
abline(v=k$theta(), col="red")

# LL
obj_fun = function(theta,...) logLikelihoodFun(k,theta)$logLikelihood
obj_grad = function(theta,...) logLikelihoodFun(k,theta,grad = T)$logLikelihoodGrad
obj_grad_fd = function(theta,...) (obj_fun(theta+0.01)-obj_fun(theta-0.01))/0.02

plot(Vectorize(obj_fun),col='blue')
for (x in seq(0.01,1,,11)){
  fx = obj_fun(x)
  gx = obj_grad(x)
  gx_ref = obj_grad_fd(x)
  arrows(x,fx,x+.1,fx+.1*gx,col='red')
  arrows(x,fx,x+.11,fx+.11*gx_ref,col='blue')
}

# LMP
obj_fun = function(theta,...) logMargPostFun(k,theta)$logMargPost
obj_grad = function(theta,...) logMargPostFun(k,theta,grad = T)$logMargPostGrad
obj_grad_fd = function(theta,...) (obj_fun(theta+0.01)-obj_fun(theta-0.01))/0.02

plot(Vectorize(obj_fun),col='blue')
for (x in seq(0.01,1,,11)){
  fx = obj_fun(x)
  gx = obj_grad(x)
  gx_ref = obj_grad_fd(x)
  arrows(x,fx,x+.1,fx+.1*gx,col='red', lty=2)
  arrows(x,fx,x+.11,fx+.11*gx_ref,col='blue')
}

library(RobustGaSP)
r = rgasp(design=X,response = y,kernel_type="matern_5_2",optimization = "lbfgs",trend=cbind(matrix(1,length(y),1),X))

obj_fun_rgasp = function(theta,...) 
  -RobustGaSP:::neg_log_marginal_post_approx_ref(param=log(1/theta),nugget=0, nugget.est=r@nugget.est, 
                                              R0=r@R0,X=r@X, zero_mean=r@zero_mean,output=r@output, 
                                              CL=r@CL, 
                                              a=0.2,
                                              b=1/(length(r@output))^{1/dim(as.matrix(r@input))[2]}*(0.2+dim(as.matrix(r@input))[2]),
                                              kernel_type=rep(as.integer(3),ncol(X)),alpha=r@alpha)

obj_grad_rgasp = function(theta,...) 
  -RobustGaSP:::neg_log_marginal_post_approx_ref_deriv(param=log(1/theta),nugget=0, nugget.est=r@nugget.est, 
                                                 R0=r@R0,X=r@X, zero_mean=r@zero_mean,output=r@output, 
                                                 CL=r@CL, 
                                                 a=0.2,
                                                 b=1/(length(r@output))^{1/dim(as.matrix(r@input))[2]}*(0.2+dim(as.matrix(r@input))[2]),
                                                 kernel_type=rep(as.integer(3),ncol(X)),alpha=r@alpha) / -theta

lines(seq(0,1,,101),(Vectorize(obj_fun_rgasp)(seq(0,1,,101))),col='red')
abline(v=1/r@beta_hat,col='red')
for (x in seq(0.01,1,,11)){
  fx = obj_fun_rgasp(x)
  gx = obj_grad_rgasp(x)
  arrows(x,fx,x+.1,fx+.1*gx,col='red', lty=1)
}
```


Fix patch:
diff --git a/src/lib/Kriging.cpp b/src/lib/Kriging.cpp
index 72728640..4bbb102b 100644
--- a/src/lib/Kriging.cpp
+++ b/src/lib/Kriging.cpp
@@ -721,10 +721,10 @@ double Kriging::_logMargPost(const arma::vec& _theta,
   t0 = Bench::toc(bench, "S2 = YtRi * y - yt * RiFFtRiFiFtRi * y", t0);
 
   if (fd->is_sigma2_estim)  // means no sigma2 provided
-    fd->sigma2 = S_2(0, 0) / (n - d);
+    fd->sigma2 = S_2(0, 0) / (n - m_F.n_cols);
 
   double log_S_2 = log(S_2(0, 0));
-  double log_marginal_lik = -sum(log(L.diag())) - sum(log(LX.diag())) - (n - d) / 2.0 * log_S_2;
+  double log_marginal_lik = -sum(log(L.diag())) - sum(log(LX.diag())) - (n - m_F.n_cols) / 2.0 * log_S_2;
   t0 = Bench::toc(bench, "lml = -Sum(log(diag(T))) - Sum(log(diag(TF)))...", t0);
   // arma::cout << " log_marginal_lik:" << log_marginal_lik << arma::endl;
 
diff --git a/src/lib/NuggetKriging.cpp b/src/lib/NuggetKriging.cpp
index cbcd6b9c..76721c63 100644
--- a/src/lib/NuggetKriging.cpp
+++ b/src/lib/NuggetKriging.cpp
@@ -426,15 +426,15 @@ double NuggetKriging::_logMargPost(const arma::vec& _theta_alpha,
   // arma::cout << " S_2:" << S_2 << arma::endl;
 
   if (fd->is_sigma2_estim)
-    fd->sigma2 = S_2(0, 0) / (n - d);
+    fd->sigma2 = S_2(0, 0) / (n - m_F.n_cols);
 
   if (fd->is_nugget_estim)
-    fd->nugget = S_2(0, 0) / (n - d) * (1 / _alpha - 1);
+    fd->nugget = S_2(0, 0) / (n - m_F.n_cols) * (1 / _alpha - 1);
 
-  fd->var = S_2(0, 0) / (n - d) / _alpha;
+  fd->var = S_2(0, 0) / (n - m_F.n_cols) / _alpha;
 
   double log_S_2 = log(S_2(0, 0));
-  double log_marginal_lik = -sum(log(L.diag())) - sum(log(LX.diag())) - (n - d) / 2.0 * log_S_2;
+  double log_marginal_lik = -sum(log(L.diag())) - sum(log(LX.diag())) - (n - m_F.n_cols) / 2.0 * log_S_2;
   t0 = Bench::toc(bench, "lml = -Sum(log(diag(T))) - Sum(log(diag(TF)))...", t0);
   // arma::cout << " log_marginal_lik:" << log_marginal_lik << arma::endl;
 
@@ -501,7 +501,7 @@ double NuggetKriging::_logMargPost(const arma::vec& _theta_alpha,
              - gradR_k * R_inv_X_Xt_R_inv_X_inv_Xt_R_inv;
       t0 = Bench::toc(bench, "Wb_k = gradR_k \\ T \\ Tt - gradR_k * RiFFtRiFiFtRi", t0);
 
-      ans[k] = -0.5 * sum(Wb_k.diag()) + (n - d) / 2.0 * (trans(m_y) * trans(Wb_k) * Q_output / S_2(0, 0))[0];
+      ans[k] = -0.5 * sum(Wb_k.diag()) + (n - m_F.n_cols) / 2.0 * (trans(m_y) * trans(Wb_k) * Q_output / S_2(0, 0))[0];
       t0 = Bench::toc(bench, "ans[k] = Sum(diag(Wb_k)) + yt * Wb_kt * Qo / S2...", t0);
     }
     // arma::cout << " log_marginal_lik_deriv:" << -ans * pow(_theta,2) << arma::endl;
@@ -513,7 +513,8 @@ double NuggetKriging::_logMargPost(const arma::vec& _theta_alpha,
     Wb_k = trans(
                solve(trans(L), solve(L, gradR_d, LinearAlgebra::default_solve_opts), LinearAlgebra::default_solve_opts))
            - gradR_d * R_inv_X_Xt_R_inv_X_inv_Xt_R_inv;
-    double ans_d = -0.5 * sum(Wb_k.diag()) + (n - d) / 2.0 * (trans(m_y) * trans(Wb_k) * Q_output / S_2(0, 0))[0];
+    double ans_d
+        = -0.5 * sum(Wb_k.diag()) + (n - m_F.n_cols) / 2.0 * (trans(m_y) * trans(Wb_k) * Q_output / S_2(0, 0))[0];
 
     (*grad_out).at(d) = ans_d - (a / t - b) / pow(_alpha, 2.0);
 


Test patch:
diff --git a/bindings/R/rlibkriging/tests/testthat/test-RobustGaSPtrendlinear.R b/bindings/R/rlibkriging/tests/testthat/test-RobustGaSPtrendlinear.R
new file mode 100644
index 00000000..28def647
--- /dev/null
+++ b/bindings/R/rlibkriging/tests/testthat/test-RobustGaSPtrendlinear.R
@@ -0,0 +1,102 @@
+library(testthat)
+ Sys.setenv('OMP_THREAD_LIMIT'=2)
+library(rlibkriging)
+
+library(RobustGaSP)
+
+context("RobustGaSP / Fit: 1D")
+
+f = function(x) 1-1/2*(sin(12*x)/(1+x)+2*cos(7*x)*x^5+0.7)
+#plot(f)
+n <- 5
+set.seed(123)
+X <- as.matrix(runif(n))
+y = f(X)
+#points(X,y)
+k = RobustGaSP::rgasp(design=X,response=y, trend=cbind(matrix(1,length(y),1),X))
+#library(rlibkriging)
+r <- Kriging(y, X,
+  kernel="matern5_2",
+  regmodel = "linear", normalize = FALSE,
+  optim = "BFGS",
+  objective = "LMP")
+# m = as.list(r)
+
+# Check lmp function
+
+lmp_rgasp = function(X, model=k) {if (!is.matrix(X)) X = matrix(X,ncol=1);
+                  # print(dim(X));
+                  apply(X,1,
+                    function(x) {
+                      #y=-logMargPostFun(r,matrix(unlist(x),ncol=2))$logMargPost
+                      y=RobustGaSP:::neg_log_marginal_post_approx_ref(param=(x),nugget=0, nugget.est=model@nugget.est, 
+                                        R0=model@R0,X=model@X, zero_mean=model@zero_mean,output=model@output, 
+                                        CL=model@CL, 
+                                        a=0.2,
+                                        b=1/(length(model@output))^{1/dim(as.matrix(model@input))[2]}*(0.2+dim(as.matrix(model@input))[2]),
+                                        kernel_type=rep(as.integer(3),ncol(X)),alpha=model@alpha
+                                        )
+                      y})}
+lmp_rgasp(1)
+
+plot(lmp_rgasp,xlim=c(0.01,6))
+abline(v=(log(k@beta_hat)))
+
+lmp_lk = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=1);
+                  # print(dim(X));
+                  apply(X,1,
+                    function(x) {
+                      y=-logMargPostFun(r,matrix(unlist(exp(-(x))),ncol=1))$logMargPost
+                      y})}
+lmp_lk(1)
+
+lines(seq(0.1,6,,101),lmp_lk(seq(0.1,6,,101)),col='red')
+abline(v=(log(1/as.list(r)$theta)),col='red')
+
+precision <- 1e-3
+test_that(desc=paste0("RobustGaSP / Fit: 1D / rgasp/lmp is the same that lk/lmp one"),
+          expect_equal(lmp_rgasp(1),lmp_lk(1),tol = precision))
+test_that(desc=paste0("RobustGaSP / Fit: 1D / fitted theta is the same that RobustGaSP one"),
+          expect_equal(as.list(r)$theta[1],1/k@beta_hat,tol = precision))
+
+
+
+dlmp_rgasp = function(X, model=k) {if (!is.matrix(X)) X = matrix(X,ncol=1);
+                  # print(dim(X));
+                  apply(X,1,
+                    function(x) {
+
+#    print(RobustGaSP:::log_marginal_lik_deriv(param=(x),nugget=0,nugget_est=model@nugget.est, 
+#                                        R0=model@R0,X=model@X, zero_mean=model@zero_mean,
+#                                        output=model@output, 
+#                                        kernel_type=rep(as.integer(3),ncol(X)),alpha=model@alpha))
+#
+#    print(RobustGaSP:::log_approx_ref_prior_deriv(param=(x),nugget=0, nugget_est=model@nugget.est, 
+#                                        CL=model@CL, 
+#                                        a=0.2,
+#                                        b=1/(length(model@output))^{1/dim(as.matrix(model@input))[2]}*(0.2+dim(as.matrix(model@input))[2])))
+
+
+                      #y=-logMargPostFun(r,matrix(unlist(x),ncol=2))$logMargPost
+                      y=RobustGaSP:::neg_log_marginal_post_approx_ref_deriv(param=(x),nugget=0, nugget.est=model@nugget.est, 
+                                        R0=model@R0,X=model@X, zero_mean=model@zero_mean,output=model@output, 
+                                        CL=model@CL, 
+                                        a=0.2,
+                                        b=1/(length(model@output))^{1/dim(as.matrix(model@input))[2]}*(0.2+dim(as.matrix(model@input))[2]),
+                                        kernel_type=rep(as.integer(3),ncol(X)),alpha=model@alpha
+                                        )
+                      y})}
+dlmp_rgasp(1)
+
+dlmp_lk = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=1);
+                  apply(X,1,
+                    function(x) {
+                      y=-logMargPostFun(r,matrix(unlist(exp(-(x))),ncol=1),TRUE)$logMargPostGrad
+                      y})}
+-exp(-1)*dlmp_lk(1)
+
+precision <- 1e-3
+test_that(desc=paste0("RobustGaSP / Fit: 1D / rgasp/lmp deriv is the same that lk/lmp deriv"),
+          expect_equal(dlmp_rgasp(1),-exp(-1)*dlmp_lk(1),tol = precision))
+
+