Running 'test-KrigingFit.R' [174s]
  Running 'test-KrigingLeaveOneOut.R' [172s]
  Running 'test-KrigingLeaveOneOut_3d.R' [2s]
  Running 'test-KrigingLogLik.R' [158s]
  Running 'test-KrigingLogLikGradHess.R' [10s]
  Running 'test-KrigingMethods.R' [6s]
  Running 'test-KrigingPredict.R' [62s]
  Running 'test-NoiseKrigingFit.R' [148s]
  Running 'test-NoiseKrigingLogLik.R' [168s]
  Running 'test-NoiseKrigingMethods.R' [173s]
  Running 'test-NoiseKrigingPredict.R' [20s]
  Running 'test-NuggetKrigingFit.R' [5s]
  Running 'test-NuggetKrigingLogLik.R' [174s]
  Running 'test-NuggetKrigingLogMargPost.R' [10s]
  Running 'test-NuggetKrigingMethods.R' [6s]
  Running 'test-NuggetKrigingPredict.R' [181s]
  Running 'test-RobustGaSP-Nugget.R' [2s]
  Running 'test-RobustGaSP.R' [3s]
  Running 'test-RobustGaSPvsKrigingLMP.R' [171s]
  Running 'test-RobustGaSPvsNuggetKrigingLMP.R' [2s]
  Running 'test-SaveLoad.R' [177s]
  Running 'test-asDiceKriging.R' [13s]
  Running 'test-estimnone.R' [3s]
  Running 'test-normalize.R' [4s]
  Running 'test-rlibkriging-demo.R' [170s]
Running the tests in 'tests/test-KrigingFit.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > context("Fit: 1D")
  > 
  > 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))
  > y = f(X)
  > k = NULL
  > r = NULL
  > k = DiceKriging::km(design=X,response=y,covtype = "gauss",control = list(trace=F))
  > r <- Kriging(y, X, "gauss")
  > 
  > ll = Vectorize(function(x) logLikelihoodFun(r,x)$logLikelihood)
  > plot(ll,xlim=c(0.000001,10))
  >   for (x in seq(0.000001,10,,11)){
  +     envx = new.env()
  +     ll2x = logLikelihoodFun(r,x)$logLikelihood
  +     gll2x = logLikelihoodFun(r,x,grad = T)$logLikelihoodGrad
  +     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
  +   }
  > 
  > theta_ref = optimize(ll,interval=c(0.001,10),maximum=T)$maximum
  > abline(v=theta_ref,col='black')
  > abline(v=as.list(r)$theta,col='red')
  > abline(v=k@covariance@range.val,col='blue')
  > 
  > test_that(desc="Fit: 1D / fit of theta by DiceKriging is right",
  +           expect_equal(theta_ref, k@covariance@range.val, tol= 1e-3))
  Test passed 😸
  > 
  > test_that(desc="Fit: 1D / fit of theta by libKriging is right",
  +           expect_equal(array(theta_ref), array(as.list(r)$theta), tol= 0.01))
  Test passed 🥳
  > 
  > #############################################################
  > 
  > context("Fit: 2D (Branin)")
  > 
  > f = function(X) apply(X,1,DiceKriging::branin)
  > n <- 15
  > set.seed(1234)
  > X <- cbind(runif(n),runif(n))
  > y = f(X)
  > k = NULL
  > r = NULL
  > k = DiceKriging::km(design=X,response=y,covtype = "gauss",control = list(trace=F),parinit = c(.2,.5))
  > r <- Kriging(y, X, "gauss", parameters=list(theta=matrix(c(.2,.5),ncol=2)))
  > 
  > ll = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
  +                   # print(dim(X));
  +                   apply(X,1,
  +                     function(x) {
  +                       # print(dim(x))
  +                       #print(matrix(unlist(x),ncol=2));
  +                       y=-logLikelihoodFun(r,matrix(unlist(x),ncol=2))$logLikelihood
  +                       #print(y);
  +                       y})}
  > #DiceView::contourview(ll,xlim=c(0.01,2),ylim=c(0.01,2))
  > x=seq(0.01,2,,51)
  > contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  > 
  > theta_ref = optim(par=matrix(c(.2,.5),ncol=2),ll,lower=c(0.01,0.01),upper=c(2,2),method="L-BFGS-B")$par
  > points(theta_ref,col='black')
  > points(as.list(r)$theta[1],as.list(r)$theta[2],col='red')
  > points(k@covariance@range.val[1],k@covariance@range.val[2],col='blue')
  > 
  > test_that(desc="Fit: 2D (Branin) / fit of theta 2D is _quite_ the same that DiceKriging one",
  +           expect_equal(ll(array(as.list(r)$theta)), ll(k@covariance@range.val), tol=1e-1))
  Test passed 😸
  > 
  > 
  > 
  > #############################################################
  > 
  > context("Fit: 2D (Branin) multistart")
  > 
  > f = function(X) apply(X,1,DiceKriging::branin)
  > n <- 15
  > set.seed(1234)
  > X <- cbind(runif(n),runif(n))
  > y = f(X)
  > k = NULL
  > r = NULL
  > 
  > parinit = matrix(runif(10*ncol(X)),ncol=ncol(X))
  > k <- tryCatch( # needed to catch warning due to %dopar% usage when using multistart
  +     withCallingHandlers(
  +       {
  +         error_text <- "No error."
  +         DiceKriging::km(design=X,response=y,covtype = "gauss",multistart = 10, parinit=parinit,control = list(trace=F))
  +       },
  +       warning = function(e) {
  +         error_text <<- trimws(paste0("WARNING: ", e))
  +         invokeRestart("muffleWarning")
  +       }
  +     ),
  +     error = function(e) {
  +       return(list(value = NA, error_text = trimws(paste0("ERROR: ", e))))
  +     },
  +     finally = {
  +     }
  +   )
  > r <- Kriging(y, X, "gauss", parameters=list(theta=parinit))
  > l = as.list(r)
  > 
  > # save(list=ls(),file="fit-2d-multistart.Rdata")
  > 
  > ll = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
  + # print(dim(X));
  + apply(X,1,
  +       function(x) {
  +         # print(dim(x))
  +         #print(matrix(unlist(x),ncol=2));
  +         y=-logLikelihoodFun(r,matrix(unlist(x),ncol=2))$logLikelihood
  +         #print(y);
  +         y})}
  > #DiceView::contourview(ll,xlim=c(0.01,2),ylim=c(0.01,2))
  > x=seq(0.01,2,,51)
  > contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  > 
  > theta_ref = optim(par=matrix(c(.2,.5),ncol=2),ll,lower=c(0.01,0.01),upper=c(2,2),method="L-BFGS-B")$par
  > points(theta_ref,col='black')
  > points(as.list(r)$theta[1],as.list(r)$theta[2],col='red')
  > points(k@covariance@range.val[1],k@covariance@range.val[2],col='blue')
  > 
  > test_that(desc="Fit: 2D (Branin) multistart / fit of theta 2D is _quite_ the same that DiceKriging one",
  +           expect_equal(ll(array(as.list(r)$theta)), ll(k@covariance@range.val), tol= 1e-3))
  Test passed 🌈
  > 
  > 
  > 
  > #############################################################
  > 
  > context("Fit: 2D")
  > 
  > f <- function(X) apply(X, 1,
  +                        function(x)
  +                          prod(
  +                            sin(2*pi*
  +                                  ( x * (seq(0,1,l=1+length(x))[-1])^2 )
  +                            )))
  > logn <- 1 #seq(1, 2.5, by=.1)
  > n <- floor(10^logn)
  > d <- 2
  > set.seed(1234)
  > X <- matrix(runif(n*d),ncol=d)
  > y <- f(X)
  > k = NULL
  > r = NULL
  > k = DiceKriging::km(design=X,response=y,covtype = "gauss",control = list(trace=F))
  > 
  > x=seq(0,2,,51)
  > mll_fun <- function(x) -apply(x,1,
  +                               function(theta)
  +                                 DiceKriging::logLikFun(theta,k)
  + )
  > contour(x,x,matrix(mll_fun(expand.grid(x,x)),nrow=length(x)),nlevels = 30)
  > 
  > # use same startup point for convergence
  > r <- Kriging(y, X, "gauss","constant",FALSE,"BFGS","LL",
  +              parameters=list(theta=matrix(k@parinit,ncol=2)))
  > #mll2_fun <- function(x) -apply(x,1,
  > #                              function(theta)
  > #                                r$logLikelihoodFun(theta)$logLikelihood
  > #)
  > #contour(x,x,matrix(mll2_fun(expand.grid(x,x)),nrow=length(x)),nlevels = 30)
  > 
  > l = as.list(r)
  > 
  > # save(list=ls(),file="fit-2d.Rdata")
  > 
  > points(as.list(r)$theta[1],as.list(r)$theta[2],col='red')
  > points(k@covariance@range.val[1],k@covariance@range.val[2],col='blue')
  > 
  > test_that(desc="Fit: 2D / fit of theta 2D is the same that DiceKriging one",
  +           expect_equal(array(as.list(r)$theta),array(k@covariance@range.val),tol=  5e-2))
  Test passed 🎉
  > 
  > ################################################################################
  > 
  > context("Fit: 2D _not_ in [0,1]^2")
  > 
  > # "unnormed" version of Branin: [0,1]x[0,15] -> ...
  > branin_15 <- function (x) {
  +   x1 <- x[1] * 15 - 5
  +   x2 <- x[2] #* 15
  +   (x2 - 5/(4 * pi^2) * (x1^2) + 5/pi * x1 - 6)^2 + 10 * (1 - 1/(8 * pi)) * cos(x1) + 10
  + }
  > 
  > f = function(X) apply(X,1,branin_15)
  > n <- 15
  > set.seed(1234)
  > X <- cbind(runif(n,0,1),runif(n,0,15))
  > y = f(X)
  > k = NULL
  > r = NULL
  > k = DiceKriging::km(design=X,response=y,covtype = "gauss",control = list(trace=F),parinit = c(0.25,10))
  > r <- Kriging(y, X, "gauss",parameters=list(theta=matrix(c(0.25,10),ncol=2)))
  > l = as.list(r)
  > 
  > # save(list=ls(),file="fit-2d-not01.Rdata")
  > 
  > ll_r = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
  + # print(dim(X));
  + apply(X,1,
  +       function(x) {
  +         # print(dim(x))
  +         #print(matrix(unlist(x),ncol=2));
  +         -logLikelihoodFun(r,matrix(unlist(x),ncol=2))$logLikelihood
  +         #print(y);
  +         })}
  > #DiceView::contourview(ll,xlim=c(0.01,2),ylim=c(0.01,2))
  > x1=seq(0.001,2,,51)
  > x2=seq(0.001,30,,51)
  > contour(x1,x2,matrix(ll_r(as.matrix(expand.grid(x1,x2))),nrow=length(x1)),nlevels = 30,col='red')
  > points(as.list(r)$theta[1],as.list(r)$theta[2],col='red')
  > ll_r(t(as.list(r)$theta))
  [1] 66.31244
  > 
  > ll_k = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
  + apply(X,1,function(x) {-DiceKriging::logLikFun(x,k)})}
  > contour(x1,x2,matrix(ll_k(as.matrix(expand.grid(x1,x2))),nrow=length(x1)),nlevels = 30,add=T)
  > points(k@covariance@range.val[1],k@covariance@range.val[2])
  > ll_k(k@covariance@range.val)
  [1] 66.31255
  > 
  > theta_ref = optim(par=matrix(c(.25,10),ncol=2),ll_r,lower=c(0.001,0.001),upper=c(2,30),method="L-BFGS-B")$par
  > points(theta_ref,col='black')
  > 
  > test_that(desc="Fit: 2D _not_ in [0,1]^2 / fit of theta 2D is _quite_ the same that DiceKriging one",
  +           expect_equal(ll_r(array(as.list(r)$theta)), ll_k(k@covariance@range.val), tol=1e-1))
  Test passed 😸
  > 
  > proc.time()
     user  system elapsed 
     4.65    0.42    5.07 
Running the tests in 'tests/test-KrigingLeaveOneOut.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > for (kernel in c("gauss","exp","matern3_2","matern5_2")) {
  + # kernel = "gauss"
  +   context(paste0("Check leaveOneOut for kernel ",kernel))
  +   
  + 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 = DiceKriging::km(design=X,response=y,covtype = kernel,control = list(trace=F), formula = as.formula("~."))
  + ll = function(theta) DiceKriging::leaveOneOutFun(theta,k)
  + 
  + plot(Vectorize(ll),ylab="LL",xlab="theta",xlim=c(0.01,1))
  + for (x in seq(0.01,1,,11)){
  +   envx = new.env()
  +   llx = DiceKriging::leaveOneOutFun(x,k,envx)
  +   gllx = DiceKriging::leaveOneOutGrad(x,k,envx)
  +   arrows(x,llx,x+.1,llx+.1*gllx)
  + }
  + 
  + r <- Kriging(y, X, kernel, regmodel="linear")
  + ll2 = function(theta) leaveOneOutFun(r,theta)$leaveOneOut[1]
  + # plot(Vectorize(ll2),col='red',add=T) # FIXME fails with "error: chol(): decomposition failed"
  + for (x in seq(0.01,1,,11)){
  +   envx = new.env()
  +   ll2x = leaveOneOutFun(r,x)$leaveOneOut[1]
  +   gll2x = leaveOneOutFun(r,x,grad=T)$leaveOneOutGrad
  +   arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
  + }
  + 
  + precision <- 1e-8  # the following tests should work with it, since the computations are analytical
  + x=.5
  + xenv=new.env()
  + test_that(desc="leaveOneOut is the same that DiceKriging one",
  +          expect_equal(leaveOneOutFun(r,x)$leaveOneOut[1],DiceKriging::leaveOneOutFun(x,k,xenv),tolerance = precision))
  + 
  + test_that(desc="leaveOneOut Grad is the same that DiceKriging one",
  +           expect_equal(leaveOneOutFun(r,x,grad=T)$leaveOneOutGrad,DiceKriging::leaveOneOutGrad(x,k,xenv),tolerance= precision))
  + }
  Test passed 😸
  Test passed 🥳
  Test passed 😸
  Test passed 🥳
  Test passed 😸
  Test passed 🥳
  Test passed 😸
  Test passed 🥳
  > 
  > 
  > 
  > context(paste0("Check leaveOneOutVec"))
  >   
  > 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))
  > y = f(X)
  > r <- Kriging(y, X, kernel="gauss", regmodel="linear")
  > 
  > loo = list(mean=matrix(NA,n,1),stdev=matrix(NA,n,1))
  > for (i in 1:n) {
  +   ri <- Kriging(y[-i], X[-i,,drop=FALSE], kernel="gauss", regmodel="linear", optim="none", parameters=list(theta=r$theta(),sigma2=r$sigma2()))
  +   predi = predict(ri,X[i])
  +   loo$mean[i,1] = predi$mean
  +   loo$stdev[i,1] = predi$stdev
  + }
  > 
  > loovec = r$leaveOneOutVec(r$theta())
  > 
  > test_that(desc="leaveOneOut vector mean is equal to k-fold loo",
  +           expect_equal(loovec$mean,loo$mean,tolerance= precision))
  Test passed 🥳
  > test_that(desc="leaveOneOut vector stdev is equal to k-fold loo",
  +           expect_equal(loovec$stdev,loo$stdev,tolerance= precision))
  Test passed 🥳
  > 
  > k = DiceKriging::km(design=X,response=y,covtype = "gauss",control = list(trace=F), formula = as.formula("~."), 
  +                     coef.var=r$sigma2(), coef.cov=r$theta()[1])
  > loovec_km = DiceKriging::leaveOneOut.km(k,type="UK",trend.reestim=TRUE)
  > 
  > test_that(desc="leaveOneOut vector mean is equal to DiceKriging one",
  +           expect_equal(loovec_km$mean,loovec$mean[,1],tolerance= precision))
  Test passed 🎉
  > test_that(desc="leaveOneOut vector stdev is equal to DiceKriging one",
  +           expect_equal(loovec_km$sd,loovec$stdev[,1],tolerance= precision))
  Test passed 😸
  > 
  > proc.time()
     user  system elapsed 
     2.65    0.43    3.17 
Running the tests in 'tests/test-KrigingLogLik.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > for (kernel in c("exp","matern3_2","matern5_2","gauss")) {
  +   context(paste0("Check LogLikelihood for kernel ",kernel))
  +   
  +   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 = DiceKriging::km(design=X,response=y,covtype = kernel,control = list(trace=F))
  +   ll = function(theta) DiceKriging::logLikFun(theta,k)
  +   plot(Vectorize(ll),ylab="LL",xlab="theta",xlim=c(0.01,1))
  +   for (x in seq(0.01,1,,11)){
  +     envx = new.env()
  +     llx = DiceKriging::logLikFun(x,k,envx)
  +     gllx = DiceKriging::logLikGrad(x,k,envx)
  +     arrows(x,llx,x+.1,llx+.1*gllx)
  +   }
  +   
  +   #library(rlibkriging)
  +   r <- Kriging(y, X, kernel)
  +   ll2 = function(theta) logLikelihoodFun(r,theta)$logLikelihood
  +   # plot(Vectorize(ll2),col='red',add=T) # FIXME fails with "error: chol(): decomposition failed"
  +   for (x in seq(0.01,1,,11)){
  +     envx = new.env()
  +     ll2x = logLikelihoodFun(r,x)$logLikelihood
  +     gll2x = logLikelihoodFun(r,x,grad = T)$logLikelihoodGrad
  +     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
  +   }
  +   
  +   precision <- 1e-8  # the following tests should work with it, since the computations are analytical
  +   x=.5
  +   xenv=new.env()
  +   test_that(desc="logLik is the same that DiceKriging one", 
  +             expect_equal(logLikelihoodFun(r,x)$logLikelihood[1],DiceKriging::logLikFun(x,k,xenv),tolerance = precision))
  +   
  +   test_that(desc="logLik Grad is the same that DiceKriging one", 
  +             expect_equal(logLikelihoodFun(r,x,grad=T)$logLikelihoodGrad,DiceKriging::logLikGrad(x,k,xenv),tolerance= precision))
  + }
  Test passed 😸
  Test passed 🥳
  Test passed 😸
  Test passed 🥳
  Test passed 😸
  Test passed 🥳
  Test passed 😸
  Test passed 🥳
  > 
  > 
  > ########################## 2D
  > 
  > for (kernel in c("matern3_2","matern5_2","gauss","exp")) {
  +   context(paste0("Check LogLikelihood for kernel ",kernel))
  +   
  +   f <- function(X) apply(X, 1, function(x) prod(sin((x-.5)^2)))
  +   n <- 100
  +   set.seed(123)
  +   X <- cbind(runif(n),runif(n),runif(n))
  +   y <- f(X)
  + 
  +   k = DiceKriging::km(design=X,response=y,covtype = kernel,control = list(trace=F))
  +   
  +   ##library(rlibkriging)
  +   r <- Kriging(y, X, kernel)
  +   
  +   precision <- 1e-8  # the following tests should work with it, since the computations are analytical
  +   x=c(.2,.5,.7)
  +   xenv=new.env()
  +   test_that(desc="logLik is the same that DiceKriging one", 
  +             expect_equal(logLikelihoodFun(r,x)$logLikelihood[1],DiceKriging::logLikFun(x,k,xenv),tolerance = precision))
  +   
  +   test_that(desc="logLik Grad is the same that DiceKriging one", 
  +             expect_equal(logLikelihoodFun(r,x,grad=T)$logLikelihoodGrad[1,],t(DiceKriging::logLikGrad(x,k,xenv))[1,],tolerance= precision))
  + }
  Test passed 🥳
  Test passed 🥇
  Test passed 🥳
  Test passed 🥇
  Test passed 🥳
  Test passed 🥇
  Test passed 🥳
  Test passed 🥇
  > 
  > proc.time()
     user  system elapsed 
     3.93    0.57    4.53 
Running the tests in 'tests/test-NoiseKrigingFit.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > context("Fit: 1D")
  > 
  > 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))
  > y = f(X) + 0.1*rnorm(nrow(X))
  > k = NULL
  > r = NULL
  > k = DiceKriging::km(design=X,response=y,noise.var=rep(0.1^2,nrow(X)),covtype = "gauss",control = list(trace=F),nugget.estim=F,optim.method='BFGS',multistart = 20)
  Warning message:
  executing %dopar% sequentially: no parallel backend registered 
  > r <- NoiseKriging(y,rep(0.1^2,nrow(X)), X, "gauss", optim = "BFGS")
  > l = as.list(r)
  > 
  > ll = Vectorize(function(x) logLikelihoodFun(r,c(x,k@covariance@sd2))$logLikelihood)
  > plot(ll,xlim=c(0.000001,1))
  >   for (x in seq(0.000001,1,,11)){
  +     envx = new.env()
  +     ll2x = logLikelihoodFun(r,c(x,k@covariance@sd2))$logLikelihood
  +     gll2x = logLikelihoodFun(r,c(x,k@covariance@sd2),grad = T)$logLikelihoodGrad[1]
  +     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
  +   }
  > 
  > theta_ref = optimize(ll,interval=c(0.001,1),maximum=T)$maximum
  > abline(v=theta_ref,col='black')
  > abline(v=as.list(r)$theta,col='red')
  > abline(v=k@covariance@range.val,col='blue')
  > 
  > test_that(desc="Noise / Fit: 1D / fit of theta by DiceKriging is right",
  +           expect_equal(theta_ref, k@covariance@range.val, tol= 1e-3))
  Test passed 🎉
  > 
  > test_that(desc="Noise / Fit: 1D / fit of theta by libKriging is right",
  +           expect_equal(array(theta_ref), array(as.list(r)$theta), tol= 0.01))
  Test passed 😀
  > 
  > #############################################################
  > 
  > context("Fit: 2D (Branin)")
  > 
  > f = function(X) apply(X,1,DiceKriging::branin)
  > n <- 15
  > set.seed(1234)
  > X <- cbind(runif(n),runif(n))
  > y = f(X)+ 10*rnorm(nrow(X))
  > k = NULL
  > r = NULL
  > k = DiceKriging::km(design=X,response=y,noise.var=rep(10^2,nrow(X)),covtype = "gauss",control = list(trace=F),nugget.estim=F,optim.method='BFGS',multistart = 20)
  > r <- NoiseKriging(y, noise=rep(10^2,nrow(X)),X, "gauss", optim = "BFGS")
  > #plot(Vectorize(function(a) r$logLikelihoodFun(c(r$theta(),a))$logLikelihood))
  > l = as.list(r)
  > 
  > # save(list=ls(),file="fit-nugget-2d.Rdata")
  > 
  > sigma2_k = k@covariance@sd2
  > sigma2_r = as.list(r)$sigma2
  > test_that(desc="Noise / Fit: 2D (Branin) / fit of LL by DiceKriging is same that libKriging",
  +           expect_equal(k@logLik,r$logLikelihood(), tol= 1e-2))
  Test passed 😀
  > 
  > ll = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
  +                   # print(dim(X));
  +                   apply(X,1,
  +                     function(x) {
  +                       y=-logLikelihoodFun(r,c(unlist(x),sigma2_k))$logLikelihood
  +                       #print(y);
  +                       y})}
  > #DiceView::contourview(ll,xlim=c(0.1,2),ylim=c(0.1,2))
  > x=seq(0.1,1,,51)
  > contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  > 
  > theta_ref = optim(par=matrix(c(.2,.5),ncol=2),ll,lower=c(0.1,0.1),upper=c(2,2),method="L-BFGS-B")$par
  > points(theta_ref,col='black')
  > points(as.list(r)$theta[1],as.list(r)$theta[2],col='red')
  > points(k@covariance@range.val[1],k@covariance@range.val[2],col='blue')
  > 
  > test_that(desc="Noise / Fit: 2D (Branin) / fit of theta 2D is _quite_ the same that DiceKriging one",
  +           expect_equal(ll(array(as.list(r)$theta)), ll(k@covariance@range.val), tol=1e-1))
  Test passed 🎉
  > 
  > 
  > 
  > #############################################################
  > 
  > context("Fit: 2D (Branin) multistart")
  > 
  > f = function(X) apply(X,1,DiceKriging::branin)
  > n <- 15
  > set.seed(1234)
  > X <- cbind(runif(n),runif(n))
  > y = f(X) + 10*rnorm(nrow(X))
  > k = NULL
  > r = NULL
  > 
  > parinit = matrix(runif(10*ncol(X)),ncol=ncol(X))
  > k <- tryCatch( # needed to catch warning due to %dopar% usage when using multistart
  +     withCallingHandlers(
  +       {
  +         error_text <- "No error."
  +         DiceKriging::km(design=X,response=y,noise.var=rep(10^2,nrow(X)),covtype = "gauss", parinit=parinit,control = list(trace=F),nugget.estim=F,optim.method='BFGS',multistart = 20)
  +       },
  +       warning = function(e) {
  +         error_text <<- trimws(paste0("WARNING: ", e))
  +         invokeRestart("muffleWarning")
  +       }
  +     ),
  +     error = function(e) {
  +       return(list(value = NA, error_text = trimws(paste0("ERROR: ", e))))
  +     },
  +     finally = {
  +     }
  +   )
  > r <- NoiseKriging(y,noise=rep(10^2,nrow(X)), X, "gauss", parameters=list(theta=parinit))
  > l = as.list(r)
  > 
  > # save(list=ls(),file="fit-nugget-multistart.Rdata")
  > 
  > sigma2_k = k@covariance@sd2
  > sigma2_r = as.list(r)$sigma2
  > test_that(desc="Noise / Fit: 2D (Branin) multistart / fit of LL by DiceKriging is same that libKriging",
  +           expect_equal(k@logLik,r$logLikelihood(), tol= 0.01))
  Test passed 🌈
  > 
  > ll = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
  + # print(dim(X));
  + apply(X,1,
  +       function(x) {
  +         # print(dim(x))
  +         #print(matrix(unlist(x),ncol=2));
  +         y=-logLikelihoodFun(r,c(unlist(x),sigma2_k))$logLikelihood
  +         #print(y);
  +         y})}
  > #DiceView::contourview(ll,xlim=c(0.1,2),ylim=c(0.1,2))
  > x=seq(0.1,2,,51)
  > contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  > 
  > theta_ref = optim(par=matrix(c(.2,.5),ncol=2),ll,lower=c(0.1,0.1),upper=c(2,2),method="L-BFGS-B")$par
  > points(theta_ref,col='black')
  > points(as.list(r)$theta[1],as.list(r)$theta[2],col='red')
  > points(k@covariance@range.val[1],k@covariance@range.val[2],col='blue')
  > 
  > test_that(desc="Noise / Fit: 2D (Branin) multistart / fit of theta 2D is _quite_ the same that DiceKriging one",
  +           expect_equal(ll(array(as.list(r)$theta)), ll(k@covariance@range.val), tol= 1e-1))
  Test passed 😀
  > 
  > 
  > ################################################################################
  > 
  > context("Fit: 2D _not_ in [0,1]^2")
  > 
  > # "unnormed" version of Branin: [0,1]x[0,15] -> ...
  > branin_15 <- function (x) {
  +   x1 <- x[1] * 15 - 5
  +   x2 <- x[2] #* 15
  +   (x2 - 5/(4 * pi^2) * (x1^2) + 5/pi * x1 - 6)^2 + 10 * (1 - 1/(8 * pi)) * cos(x1) + 10
  + }
  > 
  > f = function(X) apply(X,1,branin_15)
  > n <- 15
  > set.seed(1234)
  > X <- cbind(runif(n,0,1),runif(n,0,15))
  > y = f(X) + 10*rnorm(nrow(X))
  > k = NULL
  > r = NULL
  > k = DiceKriging::km(design=X,response=y,noise.var=rep(10^2,nrow(X)),covtype = "gauss",control = list(trace=F),nugget.estim=FALSE,optim="BFGS",multistart=20)#,parinit = c(0.5,5))
  > r <- NoiseKriging(y,noise=rep(10^2,nrow(X)), X, "gauss",, optim = "BFGS")#, parameters=list(theta=matrix(c(0.5,5),ncol=2)))
  > l = as.list(r)
  > 
  > # save(list=ls(),file="fit-nugget-2d-not01.Rdata")
  > 
  > sigma2_k = k@covariance@sd2
  > sigma2_r = as.list(r)$sigma2
  > test_that(desc="Noise / Fit: 2D _not_ in [0,1]^2 / fit of LL by DiceKriging is same that libKriging",
  +           expect_equal(k@logLik,r$logLikelihood(), tol= 0.01))
  Test passed 😀
  > 
  > ll_r = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
  + # print(dim(X));
  + apply(X,1,
  +       function(x) {
  +         # print(dim(x))
  +         #print(matrix(unlist(x),ncol=2));
  +         -logLikelihoodFun(r,c(unlist(x),sigma2_k))$logLikelihood
  +         #print(y);
  +         })}
  > #DiceView::contourview(ll,xlim=c(0.1,2),ylim=c(0.1,2))
  > x1=seq(0.001,2,,51)
  > x2=seq(0.001,30,,51)
  > contour(x1,x2,matrix(ll_r(as.matrix(expand.grid(x1,x2))),nrow=length(x1)),nlevels = 30,col='red')
  > points(as.list(r)$theta[1],as.list(r)$theta[2],col='red')
  > ll_r(t(as.list(r)$theta))
  [1] 71.67573
  > 
  > ll_k = function(X) {if (!is.matrix(X)) X = matrix(X,ncol=2);
  + apply(X,1,function(x) {-DiceKriging::logLikFun(c(x,sigma2_k),k)})}
  > contour(x1,x2,matrix(ll_k(as.matrix(expand.grid(x1,x2))),nrow=length(x1)),nlevels = 30,add=T)
  > points(k@covariance@range.val[1],k@covariance@range.val[2])
  > ll_k(k@covariance@range.val)
  [1] 71.62632
  > 
  > theta_ref = optim(par=matrix(c(.2,10),ncol=2),ll_r,lower=c(0.001,0.001),upper=c(2,30),method="L-BFGS-B")$par
  > points(theta_ref,col='black')
  > 
  > test_that(desc="Noise / Fit: 2D _not_ in [0,1]^2 / fit of theta 2D is _quite_ the same that DiceKriging one",
  +           expect_equal(ll_r(array(as.list(r)$theta)), ll_k(k@covariance@range.val), tol=1e-1))
  Test passed 🎉
  > 
  > proc.time()
     user  system elapsed 
     4.90    0.31    5.21 
Running the tests in 'tests/test-NoiseKrigingLogLik.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > for (kernel in c("exp","matern3_2","matern5_2","gauss")) {
  +   context(paste0("Check LogLikelihood for kernel ",kernel))
  +   
  + 
  + ##library(rlibkriging, lib.loc="bindings/R/Rlibs")
  + #rlibkriging:::optim_log(3)
  + #kernel="exp"
  + 
  +   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))
  +   y = f(X) + 0.1*rnorm(nrow(X))
  + 
  + tmin=0.01
  + tmax=1
  + 
  +   k = DiceKriging::km(design=X,response=y,noise.var=rep(0.1^2,nrow(X)),covtype = kernel,control = list(trace=F))
  +   ll_k = function(theta_sigma2) apply(theta_sigma2,1,function(...)DiceKriging::logLikFun(...,k))
  +   x=seq(tmin,tmax,,51)
  +   contour(x,x,matrix(ll_k(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  +   for (x1 in seq(tmin,tmax,,11)){
  +   for (x2 in seq(tmin,tmax,,11)){
  +     envx = new.env()
  +     llx = DiceKriging::logLikFun(c(x1,x2),k,envx)
  +     gllx = DiceKriging::logLikGrad(c(x1,x2),k,envx)
  +     arrows(x1,x2,x1+0.1*gllx[1],x2+.01*gllx[2])
  +   }}
  +   
  +   #library(rlibkriging)
  +   r <- NoiseKriging(y,noise=rep(0.1^2,nrow(X)), X, kernel)
  +   ll_r = function(theta_sigma2) logLikelihoodFun(r,theta_sigma2)$logLikelihood
  +   x=seq(tmin,tmax,,51)
  +   contour(x,x,matrix(ll_r(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  +   for (x1 in seq(tmin,tmax,,11)){
  +   for (x2 in seq(tmin,tmax,,11)){
  +     envx = new.env()
  +     llx = logLikelihoodFun(r,c(x1,x2))$logLikelihood
  +     gllx = logLikelihoodFun(r,c(x1,x2),grad = T)$logLikelihoodGrad
  +     arrows(x1,x2,x1+0.1*gllx[1],x2+.01*gllx[2])
  +   }}
  +   
  + 
  +   precision <- 1e-8  # the following tests should work with it, since the computations are analytical
  +   x=c(.5,.5)
  +   xenv=new.env()
  +   test_that(desc="logLik is the same that DiceKriging one", 
  +             expect_equal(logLikelihoodFun(r,x)$logLikelihood,DiceKriging::logLikFun(x,k,xenv),tolerance = precision))
  +   
  +   test_that(desc="logLik Grad is the same that DiceKriging one", 
  +             expect_equal(logLikelihoodFun(r,x,grad=T)$logLikelihoodGrad,t(DiceKriging::logLikGrad(x,k,xenv)),tolerance= precision))
  + }
  Test passed 🎉
  Test passed 😀
  Test passed 🎉
  Test passed 😀
  Test passed 🎉
  Test passed 😀
  Test passed 🎉
  Test passed 😀
  > 
  > proc.time()
     user  system elapsed 
     6.12    0.34    6.45 
Running the tests in 'tests/test-NoiseKrigingMethods.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > # f <- function(X) apply(X, 1, function(x) prod(sin((x-.5)^2)))
  > f <- function(X) apply(X, 1, function(x)
  +   prod(sin(2*pi*( x * (seq(0,1,l=1+length(x))[-1])^2 )))
  + )
  > n <- 20
  > set.seed(123)
  > X <- cbind(runif(n),runif(n))
  > y <- f(X) + rnorm(n,0,0.1)
  > d = ncol(X)
  > 
  > x=seq(0,1,,51)
  > contour(x,x,matrix(f(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  > points(X)
  > 
  > r <- NoiseKriging(y, noise=rep(0.1^2,nrow(X)), X,"gauss",parameters = list(theta=matrix(runif(40),ncol=2),sigma2=1))
  > l= as.list(r)
  > # ll = function(X) {
  > #   logLikelihoodFun(r,X,grad=F)$logLikelihood
  > # }
  > # contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  > # gll = function(X) {
  > #   logLikelihoodFun(r,X,grad=T)$logLikelihoodGrad
  > # }
  > # for (ix in 1:21) {
  > # for (iy in 1:21) {
  > #   xx = c(ix/21,iy/21)
  > #   g = gll(xx)
  > #   arrows(xx[1],xx[2],xx[1]+g[1]/1000,xx[2]+g[2]/1000,col='grey',length=0.05)
  > # }
  > # }
  > 
  > context("noise / print")
  > 
  > p = capture.output(print(r))
  > 
  > 
  > context("noise / logLikelihood")
  > 
  > ll = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r,c(theta,l$sigma2))$logLikelihood)}
  > t=seq(0.01,2,,51)
  > contour(t,t,matrix(ll(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30)
  > points(as.list(r)$theta[1],as.list(r)$theta[2])
  > # logLikelihoodFun(r,t(as.list(r)$theta))
  > 
  > test_that("noise / logLikelihoodFun returned",
  +           expect_equal(names(logLikelihoodFun(r,runif(d+1))),c("logLikelihood")))
  Test passed 😸
  > test_that("noise / logLikelihoodFun logLikelihoodGrad returned",
  +           expect_equal(names(logLikelihoodFun(r,runif(d+1),grad=TRUE)),c("logLikelihood","logLikelihoodGrad")))
  Test passed 🎊
  > 
  > test_that("noise / logLikelihoodFun dim",
  +           expect_equal(dim(logLikelihoodFun(r,rbind(runif(d+1),runif(d+1)))$logLikelihood),c(2,1)))
  Test passed 🥳
  > test_that("noise / logLikelihoodGrad dim",
  +           expect_equal(dim(logLikelihoodFun(r,rbind(runif(d+1),runif(d+1)),grad=TRUE)$logLikelihoodGrad),c(2,d+1)))
  Test passed 🌈
  > 
  > 
  > context("noise / predict")
  > 
  > test_that("noise / predict mean stdev returned",
  +           expect_equal(names(predict(r,runif(d))),c("mean","stdev")))
  Test passed 🎉
  > test_that("noise / predict mean returned",
  +           expect_equal(names(predict(r,runif(d),stdev=FALSE)),c("mean")))
  Test passed 😸
  > test_that("noise / predict mean stdev cov returned",
  +           expect_equal(names(predict(r,runif(d),cov=TRUE)),c("mean","stdev","cov")))
  Test passed 😸
  > 
  > test_that("noise / predict mean dim",
  +           expect_equal(dim(predict(r,rbind(runif(d),runif(d)))$mean),c(2,1)))
  Test passed 🎊
  > test_that("noise / predict stdev dim",
  +           expect_equal(dim(predict(r,rbind(runif(d),runif(d)))$stdev),c(2,1)))
  Test passed 🎊
  > test_that("noise / predict cov dim",
  +           expect_equal(dim(predict(r,rbind(runif(d),runif(d)),cov=TRUE)$cov),c(2,2)))
  Test passed 🎊
  > 
  > 
  > context("noise / simulate")
  > 
  > test_that("noise / simulate dim",
  +           expect_equal(dim(simulate(r,x=rbind(runif(d),runif(d)))),c(2,1)))
  Test passed 🥇
  > test_that("noise / simulate nsim dim",
  +           expect_equal(dim(simulate(r,x=rbind(runif(d),runif(d)),nsim=10)),c(2,10)))
  Test passed 🎉
  > 
  > expect_not_equal = function(x,y,...)
  +   expect_false(isTRUE(all.equal(x, y)))
  > test_that("noise / simulate mean X",
  +           expect_not_equal(mean(simulate(r,nsim = 100, x=X[1,]+0.00005)),y[1],tolerance = 0.01))
  Test passed 🌈
  > set.seed(12345)
  > x = runif(d)
  > test_that("noise / simulate mean",
  +           expect_equal(mean(simulate(r,nsim = 100, x=x)),predict(r,x)$mean[1],tolerance = 0.01))
  Test passed 🥳
  > test_that("noise / simulate sd",
  +           expect_equal(sd(simulate(r,nsim = 100, x=x)),predict(r,x)$stdev[1],tolerance = 0.01))
  Test passed 🌈
  > 
  > 
  > context("noise / update")
  > 
  > set.seed(1234)
  > X2 = matrix(runif(d*10),ncol=d)
  > y2 = f(X2) + rnorm(nrow(X2),0,0.1)
  > x=seq(0,1,,51)
  > contour(x,x,matrix(f(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  > points(X)
  > points(X2,col='red')
  > 
  > #r20 <- Kriging(c(y,y2), rbind(X,X2),"gauss")
  > #ll = function(X) {
  > #  logLikelihoodFun(r20,X,grad=F)$logLikelihood
  > #}
  > #contour(x,x,matrix(ll(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  > #gll = function(X) {
  > #  logLikelihoodFun(r20,X,grad=T)$logLikelihoodGrad
  > #}
  > #for (ix in 1:21) {
  > #for (iy in 1:21) {
  > #  xx = c(ix/21,iy/21)
  > #  g = gll(xx)
  > #  arrows(xx[1],xx[2],xx[1]+g[1]/1000,xx[2]+g[2]/1000,col='grey',length=0.05)
  > #}
  > #}
  > 
  > 
  > 
  > t=seq(0.01,2,,51)
  > contour(t,t,matrix(ll(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30)
  > points(as.list(r)$theta[1],as.list(r)$theta[2],pch=20)
  > 
  > #cat("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!","\n")
  > #cat(paste0(collapse="\n",p),"\n")
  > #cat(r$logLikelihood(),"\n")
  > #
  > #rlibkriging:::optim_log(1)
  > 
  > set.seed(1234)
  > r2 <- NoiseKriging(c(y,y2), noise=rep(0.1^2,nrow(X)+nrow(X2)),rbind(X,X2),"gauss", parameters = list(theta=matrix(as.list(r)$theta,ncol=2),sigma2=1))
  > ll2 = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r2,c(theta,as.list(r2)$sigma2))$logLikelihood)}
  > contour(t,t,matrix(ll2(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30,add=T,col='red')
  > points(as.list(r2)$theta[1],as.list(r2)$theta[2],col='red',pch=20)
  > 
  > p2 = capture.output(print(r2))
  > 
  > update(object=r,y2,newnoise=rep(0.1^2,nrow(X2)),X2)
  > llu = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r,  c(theta,as.list(r2)$sigma2) )$logLikelihood)}
  > contour(t,t,matrix(llu(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30,add=T,col='blue')
  > points(as.list(r)$theta[1],as.list(r)$theta[2],col='blue',pch=20)
  > 
  > pu = capture.output(print(r))
  > 
  > test_that("noise / update",
  +           expect_false(all(p == pu)))
  Test passed 🌈
  > 
  > # cat(paste0(collapse="\n",p2),"\n")
  > # cat(r2$logLikelihood(),"\n")
  > # cat(paste0(collapse="\n",capture.output(print( r2$logLikelihoodFun(c(r2$theta(),r2$sigma2()/(r2$sigma2()+r2$noise())),TRUE) ))),"\n")
  > # cat(paste0(collapse="\n",pu),"\n")
  > # cat(r$logLikelihood(),"\n")
  > # cat(paste0(collapse="\n",capture.output(print( r$logLikelihoodFun(c(r$theta(),r$sigma2()/(r$sigma2()+r$noise())),TRUE) ))),"\n")
  > 
  > test_that("noise / update almost converge",
  +           expect_equal(as.list(r2)$theta, as.list(r)$theta, tolerance = 2E-2))
  Test passed 🥳
  > 
  > proc.time()
     user  system elapsed 
     4.62    0.45    5.09 
Running the tests in 'tests/test-NuggetKrigingLogLik.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > kernel="gauss"
  > for (kernel in c("exp","matern3_2","matern5_2","gauss")) {
  +   context(paste0("Check LogLikelihood for kernel ",kernel))
  +   
  +   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 = DiceKriging::km(design=X,response=y,covtype = kernel,control = list(trace=F),nugget=0, nugget.estim=TRUE)
  +   alpha0 = k@covariance@sd2/(k@covariance@sd2+k@covariance@nugget)
  +   ll_theta = function(theta) DiceKriging::logLikFun(c(theta,alpha0),k)
  +   plot(Vectorize(ll_theta),ylab="LL",xlab="theta",xlim=c(0.01,1))
  +   for (x in seq(0.01,1,,21)){
  +     envx = new.env()
  +     llx = DiceKriging::logLikFun(c(x,alpha0),k,envx)
  +     gllx = DiceKriging::logLikGrad(c(x,alpha0),k,envx)[1,]
  +     arrows(x,llx,x+.1,llx+.1*gllx)
  +   }
  +   
  +   #library(rlibkriging)
  +   r <- NuggetKriging(y, X, kernel, parameters=list(nugget=0,is_nugget_estim=TRUE))
  +   ll2_theta = function(theta) logLikelihoodFun(r,c(theta,alpha0))$logLikelihood
  +   # second arg is alpha=1 for nugget=0
  +   # plot(Vectorize(ll2),col='red'), add=T) 
  +   for (x in seq(0.01,1,,21)){
  +     envx = new.env()
  +     ll2x = logLikelihoodFun(r,c(x,alpha0))$logLikelihood
  +     gll2x = logLikelihoodFun(r,c(x,alpha0),grad = T)$logLikelihoodGrad[,1]
  +     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
  +   }
  +   
  +   theta0 = k@covariance@range.val
  +   ll_alpha = function(alpha) DiceKriging::logLikFun(c(theta0,alpha),k)
  +   plot(Vectorize(ll_alpha),ylab="LL",xlab="alpha",xlim=c(0.01,1))
  +   for (x in seq(0.01,1,,21)){
  +     envx = new.env()
  +     llx = DiceKriging::logLikFun(c(theta0,x),k,envx)
  +     gllx = DiceKriging::logLikGrad(c(theta0,x),k,envx)[2,]
  +     arrows(x,llx,x+.1,llx+.1*gllx)
  +   }
  +   ll2_alpha = function(alpha) logLikelihoodFun(r,c(theta0,alpha))$logLikelihood
  +   #plot(Vectorize(ll2_alpha),col='red',add=T)
  +   for (x in seq(0.01,1,,21)){
  +     envx = new.env()
  +     ll2x = logLikelihoodFun(r,c(theta0,x))$logLikelihood
  +     gll2x = logLikelihoodFun(r,c(theta0,x),grad = T)$logLikelihoodGrad[,2]
  +     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
  +   }
  + 
  +   precision <- 1e-8  # the following tests should work with it, since the computations are analytical
  +   x=.25
  +   xenv=new.env()
  +   test_that(desc="logLik is the same that DiceKriging one /alpha", 
  +             expect_equal(
  +               logLikelihoodFun(r,c(theta0,x))$logLikelihood[1]
  +               ,
  +               DiceKriging::logLikFun(c(theta0,x),k,xenv)
  +               ,tolerance = precision))
  +   
  +   test_that(desc="logLik Grad is just good /alpha", 
  +             expect_equal(
  +               logLikelihoodFun(r,c(theta0,x),grad=T)$logLikelihoodGrad[,2]
  +               ,
  +               -(logLikelihoodFun(r,c(theta0,x-1e-5))$logLikelihood[1]-logLikelihoodFun(r,c(theta0,x))$logLikelihood[1])/1e-5,
  +               ,tolerance= 1e-5))
  +               
  +   test_that(desc="logLik Grad is the same that DiceKriging one /alpha", 
  +             expect_equal(
  +               logLikelihoodFun(r,c(theta0,x),grad=T)$logLikelihoodGrad[,2]
  +               ,
  +               DiceKriging::logLikGrad(c(theta0,x),k,xenv)[2,]
  +               ,tolerance= precision))
  +               
  +   xenv=new.env()
  +   test_that(desc="logLik is the same that DiceKriging one /theta", 
  +             expect_equal(
  +               logLikelihoodFun(r,c(x,alpha0))$logLikelihood[1]
  +               ,
  +               DiceKriging::logLikFun(c(x,alpha0),k,xenv)
  +               ,tolerance = precision))
  + 
  +   test_that(desc="logLik Grad is just good /theta", 
  +             expect_equal(
  +               logLikelihoodFun(r,c(x,alpha0),grad=T)$logLikelihoodGrad[,1]
  +               ,
  +               (logLikelihoodFun(r,c(x+1e-5,alpha0))$logLikelihood[1]-logLikelihoodFun(r,c(x,alpha0))$logLikelihood[1])/1e-5
  +               ,tolerance= 1e-3))
  + 
  +   test_that(desc="logLik Grad is the same that DiceKriging one /theta", 
  +             expect_equal(
  +               logLikelihoodFun(r,c(x,alpha0),grad=T)$logLikelihoodGrad[,1]
  +               ,
  +               DiceKriging::logLikGrad(c(x,alpha0),k,xenv)[1,]
  +               ,tolerance= precision))
  + 
  + }
  Test passed 🥳
  Test passed 😀
  Test passed 😀
  Test passed 😸
  Test passed 😀
  Test passed 🎉
  Test passed 🥳
  Test passed 😀
  Test passed 😀
  Test passed 😸
  Test passed 😀
  Test passed 🎉
  Test passed 🥳
  Test passed 😀
  Test passed 😀
  Test passed 😸
  Test passed 😀
  Test passed 🎉
  Test passed 🥳
  Test passed 😀
  Test passed 😀
  Test passed 😸
  Test passed 😀
  Test passed 🎉
  > 
  > 
  > ########################## 2D
  > 
  > 
  > 
  > for (kernel in c("matern3_2","matern5_2","gauss","exp")) {
  +   context(paste0("Check LogLikelihood for kernel ",kernel))
  +   
  +   f <- function(X) apply(X, 1, function(x) prod(sin((x-.5)^2)))
  +   n <- 100
  +   set.seed(123)
  +   X <- cbind(runif(n),runif(n),runif(n))
  +   y <- f(X)
  + 
  +   k = DiceKriging::km(design=X,response=y,covtype = kernel,control = list(trace=F),nugget=0, nugget.estim=TRUE)
  +   
  +   ##library(rlibkriging)
  +   r <- NuggetKriging(y, X, kernel, parameters=list(nugget=0,is_nugget_estim=TRUE))
  +   
  +   precision <- 1e-8  # the following tests should work with it, since the computations are analytical
  +   #x = c(.2,.5,.7,0.01)
  +   #x=c(.5,.5,.5,.9999995)
  +   x = c(r$theta(),r$sigma2()/(r$sigma2()+r$nugget()))
  +   #x = c(k@covariance@range.val,k@covariance@sd2/(k@covariance@sd2+k@covariance@nugget))
  +   xenv=new.env()
  +   test_that(desc="logLik is the same that DiceKriging one", 
  +             expect_equal(
  +               logLikelihoodFun(r,x)$logLikelihood[1]
  +               ,
  +               DiceKriging::logLikFun(x,k,xenv)
  +               ,tolerance = precision))
  +   
  +   test_that(desc="logLik Grad is the same that DiceKriging one", 
  +             expect_equal(
  +               logLikelihoodFun(r,x,grad=T)$logLikelihoodGrad[1,]
  +               ,
  +               t(DiceKriging::logLikGrad(x,k,xenv))[1,]
  +               ,tolerance= precision))
  + 
  +   eps=0.000001
  +   test_that(desc="logLik Grad is just good", 
  +             expect_equal(
  +               -logLikelihoodFun(r,x,grad=T)$logLikelihoodGrad[1,]
  +               ,
  +               c( # finite-diff grad
  +                 (logLikelihoodFun(r,x-c(eps,0,0,0))$logLikelihood[1]-logLikelihoodFun(r,x)$logLikelihood[1])/eps,
  +                 (logLikelihoodFun(r,x-c(0,eps,0,0))$logLikelihood[1]-logLikelihoodFun(r,x)$logLikelihood[1])/eps,
  +                 (logLikelihoodFun(r,x-c(0,0,eps,0))$logLikelihood[1]-logLikelihoodFun(r,x)$logLikelihood[1])/eps,
  +                 (logLikelihoodFun(r,x-c(0,0,0,eps))$logLikelihood[1]-logLikelihoodFun(r,x)$logLikelihood[1])/eps
  +               )
  +               ,tolerance= 1e-2))
  + }
  Test passed 🎊
  Test passed 🎊
  Test passed 😸
  Test passed 🎊
  Test passed 🎊
  Test passed 😸
  Test passed 🎊
  Test passed 🎊
  Test passed 😸
  Test passed 🎊
  Test passed 🎊
  Test passed 😸
  > 
  > proc.time()
     user  system elapsed 
     6.67    0.71    7.36 
Running the tests in 'tests/test-NuggetKrigingPredict.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > for (kernel in c("gauss","exp","matern3_2","matern5_2")) {
  +   context(paste0("Check predict 1D for kernel ",kernel))
  + 
  +   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 = DiceKriging::km(design=X,response=y,covtype = "gauss",control = list(trace=F), nugget=0,nugget.estim = TRUE)
  +   #library(rlibkriging)
  +   r <- NuggetKriging(y,X,"gauss","constant",FALSE,"none","LL",
  +                parameters=list(sigma2=k@covariance@sd2,has_sigma2=TRUE, is_sigma2_estim=FALSE,
  +                theta=matrix(k@covariance@range.val),has_theta=TRUE, is_theta_estim=FALSE,
  +                nugget=k@covariance@nugget,has_nugget=TRUE, is_nugget_estim=FALSE))
  +   # m = as.list(r)
  +   
  +   ntest <- 100
  +   Xtest <- as.matrix(runif(ntest))
  +   ptest <- DiceKriging::predict(k,Xtest,type="UK",cov.compute = TRUE,checkNames=F)
  +   Yktest <- ptest$mean
  +   sktest <- ptest$sd
  +   cktest <- c(ptest$cov)
  +   Ytest <- predict(r,Xtest,TRUE,TRUE)
  +   
  +   plot(f)
  +   points(X,y)
  +   points(Xtest,Yktest,col='blue')
  +   points(Xtest,Ytest$mean,col='red')
  +   
  +   precision <- 1e-5
  +   test_that(desc=paste0("pred mean is the same that DiceKriging one:\n ",paste0(collapse=",",Yktest),"\n ",paste0(collapse=",",Ytest$mean)),
  +             expect_equal(array(Yktest),array(Ytest$mean),tol = precision))
  +   
  +   test_that(desc="pred sd is the same that DiceKriging one", 
  +             expect_equal(array(sktest),array(Ytest$stdev) ,tol = precision))
  +               
  +   test_that(desc="pred cov is the same that DiceKriging one", 
  +             expect_equal(cktest,c(Ytest$cov) ,tol = precision))
  +   
  +   plot(f)
  +   points(X,y)
  +   x=seq(0,1,,101)
  +   lines(x,predict(r,x)$mean,lty=2)
  +   polygon(
  +       c(x,rev(x)),
  +       c(predict(r,x)$mean-predict(r,x)$stdev,
  +         predict(r,rev(x))$mean+predict(r,rev(x))$stdev),
  +         col=rgb(0,0,0,.1),border=NA)
  +   #
  +   #for (i in 1:10)
  +   #    lines(x,simulate(r,x=x,seed=i),col='grey')
  +   
  +   .x = seq(0.1,0.9,,11)
  +   p_allx = predict(r,.x,TRUE, cov=FALSE, deriv=TRUE)
  +   for (i in 1:length(.x)) {
  +     # ref from DiceOptim::EI.grad
  +     newdata = .x[i]
  +     model = k
  +     T <- model@T
  +     X <- model@X
  +     z <- model@z
  +     u <- model@M
  +     covStruct <- model@covariance
  +     predx <- predict(object=model, newdata=newdata, type="UK", checkNames = FALSE,se.compute=TRUE,cov.compute=FALSE)
  +     kriging.mean <- predx$mean
  +     kriging.sd <- predx$sd
  +     v <- predx$Tinv.c
  +     c <- predx$c
  +     dc <- DiceKriging::covVector.dx(x=newdata, X=X, object=covStruct, c=c)  
  +     #d = model@d
  +     #h=sqrt(.Machine$double.eps)
  +     #A <- matrix(newdata, nrow=d, ncol=d, byrow=TRUE)
  +     #Apos <- A+h*diag(d)
  +     #Aneg <- A-h*diag(d)
  +     #newpoints <- data.frame(rbind(Apos, Aneg))
  +     #f.newdata <- model.matrix(model@trend.formula, data = newpoints)
  +     #f.deltax <- (f.newdata[1:d,]-f.newdata[(d+1):(2*d),])/(2*h)
  +     f.deltax <- DiceKriging::trend.deltax(x=newdata, model=model)
  +     W <- backsolve(t(T), dc, upper.tri=FALSE)
  +     kriging.mean.grad <- t(W)%*%z + t(model@trend.coef%*%f.deltax)
  +     tuuinv <- solve(t(u)%*%u)
  +     F.newdata <- model.matrix(model@trend.formula, data=as.data.frame(newdata))
  +     kriging.sd2.grad <-  t( -2*t(v)%*%W +
  +                                   2*(F.newdata - t(v)%*%u )%*% tuuinv %*%
  +                                   (f.deltax - t(t(W)%*%u) ))
  +     kriging.sd.grad <- kriging.sd2.grad / (2*kriging.sd)
  +                        
  +     p = predict(r,.x[i],TRUE, cov=FALSE, deriv=TRUE)
  +   
  +     test_that(desc=paste0("vect pred mean deriv is ok:\n ",paste0(collapse=",",p_allx$mean_deriv[i]),"\n ",paste0(collapse=",",p$mean_deriv)),
  +              expect_equal(array(p_allx$mean_deriv[i]),array(p$mean_deriv),tol = precision))
  +     test_that(desc=paste0("vect pred sd deriv is ok:\n ",paste0(collapse=",",p_allx$stdev_deriv[i]),"\n ",paste0(collapse=",",p$stdev_deriv)),
  +              expect_equal(array(p_allx$stdev_deriv[i]),array(p$stdev_deriv),tol = precision))
  + 
  +   arrows(.x[i],p$mean, .x[i]+0.1, p$mean+0.1*p$mean_deriv)
  +   arrows(.x[i],p$mean+p$stdev, .x[i]+0.1, p$mean+p$stdev+0.1*p$mean_deriv+0.1*p$stdev_deriv, col='darkgrey')
  +   
  +     test_that(desc=paste0("pred mean deriv is the same that DiceKriging one:\n ",paste0(collapse=",",kriging.mean.grad),"\n ",paste0(collapse=",",p$mean_deriv)),
  +               expect_equal(array(kriging.mean.grad),array(p$mean_deriv),tol = precision))
  +     
  +     test_that(desc=paste0("pred sd deriv is the same that DiceKriging one:\n ",paste0(collapse=",",kriging.sd.grad),"\n ",paste0(collapse=",",p$stdev_deriv)),
  +               expect_equal(array(kriging.sd.grad),array(p$stdev_deriv),tol = precision))
  +               
  +   }
  + }
  Test passed 🎊
  Test passed 😸
  Test passed 🥳
  Test passed 🥇
  Test passed 🎉
  Test passed 🥇
  Test passed 🎊
  Test passed 😸
  Test passed 🌈
  Test passed 🌈
  Test passed 🎉
  Test passed 🌈
  Test passed 🥇
  Test passed 😸
  Test passed 😀
  Test passed 🎊
  Test passed 🎉
  Test passed 🥳
  Test passed 😀
  Test passed 🎊
  Test passed 🥳
  Test passed 🥇
  Test passed 😸
  Test passed 🌈
  Test passed 🎊
  Test passed 🎊
  Test passed 🌈
  Test passed 😀
  Test passed 🎉
  Test passed 🎉
  Test passed 🌈
  Test passed 😀
  Test passed 🎉
  Test passed 😸
  Test passed 🌈
  Test passed 😸
  Test passed 🥇
  Test passed 🌈
  Test passed 🌈
  Test passed 🌈
  Test passed 🎉
  Test passed 😀
  Test passed 🥳
  Test passed 😸
  Test passed 🌈
  Test passed 😸
  Test passed 😀
  Test passed 🎊
  Test passed 😸
  Test passed 🥳
  Test passed 🥇
  Test passed 🎉
  Test passed 🥇
  Test passed 🎊
  Test passed 😸
  Test passed 🌈
  Test passed 🌈
  Test passed 🎉
  Test passed 🌈
  Test passed 🥇
  Test passed 😸
  Test passed 😀
  Test passed 🎊
  Test passed 🎉
  Test passed 🥳
  Test passed 😀
  Test passed 🎊
  Test passed 🥳
  Test passed 🥇
  Test passed 😸
  Test passed 🌈
  Test passed 🎊
  Test passed 🎊
  Test passed 🌈
  Test passed 😀
  Test passed 🎉
  Test passed 🎉
  Test passed 🌈
  Test passed 😀
  Test passed 🎉
  Test passed 😸
  Test passed 🌈
  Test passed 😸
  Test passed 🥇
  Test passed 🌈
  Test passed 🌈
  Test passed 🌈
  Test passed 🎉
  Test passed 😀
  Test passed 🥳
  Test passed 😸
  Test passed 🌈
  Test passed 😸
  Test passed 😀
  Test passed 🎊
  Test passed 😸
  Test passed 🥳
  Test passed 🥇
  Test passed 🎉
  Test passed 🥇
  Test passed 🎊
  Test passed 😸
  Test passed 🌈
  Test passed 🌈
  Test passed 🎉
  Test passed 🌈
  Test passed 🥇
  Test passed 😸
  Test passed 😀
  Test passed 🎊
  Test passed 🎉
  Test passed 🥳
  Test passed 😀
  Test passed 🎊
  Test passed 🥳
  Test passed 🥇
  Test passed 😸
  Test passed 🌈
  Test passed 🎊
  Test passed 🎊
  Test passed 🌈
  Test passed 😀
  Test passed 🎉
  Test passed 🎉
  Test passed 🌈
  Test passed 😀
  Test passed 🎉
  Test passed 😸
  Test passed 🌈
  Test passed 😸
  Test passed 🥇
  Test passed 🌈
  Test passed 🌈
  Test passed 🌈
  Test passed 🎉
  Test passed 😀
  Test passed 🥳
  Test passed 😸
  Test passed 🌈
  Test passed 😸
  Test passed 😀
  Test passed 🎊
  Test passed 😸
  Test passed 🥳
  Test passed 🥇
  Test passed 🎉
  Test passed 🥇
  Test passed 🎊
  Test passed 😸
  Test passed 🌈
  Test passed 🌈
  Test passed 🎉
  Test passed 🌈
  Test passed 🥇
  Test passed 😸
  Test passed 😀
  Test passed 🎊
  Test passed 🎉
  Test passed 🥳
  Test passed 😀
  Test passed 🎊
  Test passed 🥳
  Test passed 🥇
  Test passed 😸
  Test passed 🌈
  Test passed 🎊
  Test passed 🎊
  Test passed 🌈
  Test passed 😀
  Test passed 🎉
  Test passed 🎉
  Test passed 🌈
  Test passed 😀
  Test passed 🎉
  Test passed 😸
  Test passed 🌈
  Test passed 😸
  Test passed 🥇
  Test passed 🌈
  Test passed 🌈
  Test passed 🌈
  Test passed 🎉
  Test passed 😀
  Test passed 🥳
  Test passed 😸
  Test passed 🌈
  Test passed 😸
  Test passed 😀
  > 
  > test_that(desc="pred cov is the same that DiceKriging one", 
  +           expect_equal(cktest,c(Ytest$cov) ,tol = precision))
  Test passed 🎊
  > 
  > proc.time()
     user  system elapsed 
    19.46    2.21   21.78 
Running the tests in 'tests/test-RobustGaSPvsKrigingLMP.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > kernel_type = function(kernel) {
  +   if (kernel=="matern3_2") return("matern_3_2")
  +   if (kernel=="matern5_2") return("matern_5_2")
  +   stop(paste0("Cannot use ",kernel))
  + }
  > kernel_type_num = function(kernel) {
  +   if (kernel=="matern3_2") return(2)
  +   if (kernel=="matern5_2") return(3)
  +   stop(paste0("Cannot use ",kernel))
  + }
  > 
  > for (kernel in c("matern5_2","matern3_2")) {
  +   context(paste0("Check Marginal Posterior for kernel ",kernel))
  +   
  +   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,kernel_type=kernel_type(kernel))
  + 
  +   lmp = function(theta) {
  +     #cat("theta: ",theta,"\n")
  +     lml = RobustGaSP::log_marginal_lik(param=log(1/theta),nugget=k@nugget,nugget_est=k@nugget.est,
  +       R0=k@R0,X=k@X,zero_mean=k@zero_mean,output=k@output,kernel_type=kernel_type_num(kernel),alpha=k@alpha)
  +     #cat("  lml: ",lml,"\n")
  +     larp = RobustGaSP::log_approx_ref_prior(param=log(1/theta),nugget=k@nugget,nugget_est=k@nugget.est,
  +       CL=k@CL,a=0.2,b=1/(length(y))^{1/dim(as.matrix(X))[2]}*(0.2+dim(as.matrix(X))[2]))
  +     #cat("  larp: ",larp,"\n")
  +     return(lml+larp)
  +   }
  + 
  +   plot(Vectorize(lmp),ylab="LMP",xlab="theta",xlim=c(0.01,2),ylim=c(-5,5))
  +   abline(v=1/k@beta_hat)
  + 
  +   lmp_deriv = function(theta) {
  +     #cat("theta: ",theta,"\n")
  +     lml_d = RobustGaSP::log_marginal_lik_deriv(param=log(1/theta),nugget=k@nugget,nugget_est=k@nugget.est,
  +       R0=k@R0,X=k@X,zero_mean=k@zero_mean,output=k@output,kernel_type=kernel_type_num(kernel),alpha=k@alpha)
  +     #cat("  lml_d: ",lml_d,"\n")
  +     larp_d = RobustGaSP::log_approx_ref_prior_deriv(param=log(1/theta),nugget=k@nugget,nugget_est=k@nugget.est,
  +       CL=k@CL,a=0.2,b=1/(length(y))^{1/dim(as.matrix(X))[2]}*(0.2+dim(as.matrix(X))[2]))
  +     #cat("  larp_d: ",larp_d,"\n")
  +     return((lml_d + larp_d)* 1/theta * (-1/theta))
  +   }
  + 
  +   for (x in seq(0.01,2,,11)){
  +     arrows(x,lmp(x),x+.1,lmp(x)+.1*lmp_deriv(x))
  +   }
  + 
  +   #library(rlibkriging)
  +   r <- Kriging(y, X, kernel, objective="LMP")
  +   ## Should be equal:
  +   #lmp(1.0); lmp_deriv(1.0);
  +   #logMargPostFun(r,1.0,grad = T)
  +   #lmp(0.1); lmp_deriv(0.1);
  +   #logMargPostFun(r,0.1,grad = T)
  +   #ll2 = function(theta) logMargPostFun(r,theta)$logMargPost
  +   # plot(Vectorize(ll2),col='red',add=T,xlim=c(0.01,2)) # FIXME fails with "error: chol(): decomposition failed"
  +   for (x in seq(0.01,2,,11)){
  +     ll2x = logMargPostFun(r,x)$logMargPost
  +     gll2x = logMargPostFun(r,x,grad = T)$logMargPostGrad
  +     arrows(x,ll2x,x+.1,ll2x+.1*gll2x,col='red')
  +   }
  +   
  +   precision <- 1e-8  # the following tests should work with it, since the computations are analytical
  +   x=.5
  +   test_that(desc="logMargPost is the same that RobustGaSP one", 
  +             expect_equal(logMargPostFun(r,x)$logMargPost[1],lmp(x),tolerance = precision))
  +   
  +   test_that(desc="logMargPost Grad is the same that RobustGaSP one", 
  +             expect_equal(logMargPostFun(r,x,grad = T)$logMargPostGrad[1],lmp_deriv(x),tolerance= precision))
  + }
  The upper bounds of the range parameters are 184.9743 
  The initial values of range parameters are 3.699485 
  Start of the optimization  1  : 
  The number of iterations is  10 
   The value of the  marginal posterior  function is  2.497978 
   Optimized range parameters are 0.1921691 
   Optimized nugget parameter is 0 
   Convergence:  TRUE 
  The initial values of range parameters are 0.05223118 
  Start of the optimization  2  : 
  The number of iterations is  6 
   The value of the  marginal posterior  function is  2.497978 
   Optimized range parameters are 0.1921691 
   Optimized nugget parameter is 0 
   Convergence:  TRUE 
  Test passed 🥳
  Test passed 🥳
  The upper bounds of the range parameters are 12784.39 
  The initial values of range parameters are 255.6879 
  Start of the optimization  1  : 
  The number of iterations is  11 
   The value of the  marginal posterior  function is  2.654437 
   Optimized range parameters are 0.2772546 
   Optimized nugget parameter is 0 
   Convergence:  TRUE 
  The initial values of range parameters are 0.05223118 
  Start of the optimization  2  : 
  The number of iterations is  7 
   The value of the  marginal posterior  function is  2.654437 
   Optimized range parameters are 0.2772546 
   Optimized nugget parameter is 0 
   Convergence:  TRUE 
  Test passed 🥳
  Test passed 🥳
  > 
  > proc.time()
     user  system elapsed 
     1.92    0.21    2.12 
Running the tests in 'tests/test-SaveLoad.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > #library(testthat)
  > 
  > ##library(rlibkriging, lib.loc="bindings/R/Rlibs")
  > ##library(rlibkriging)
  > 
  > # f <- function(X) apply(X, 1, function(x) prod(sin((x-.5)^2)))
  > f <- function(X) apply(X, 1, function(x)
  +   prod(sin(2*pi*( x * (seq(0,1,l=1+length(x))[-1])^2 )))
  + )
  > n <- 20
  > set.seed(123)
  > X <- cbind(runif(n),runif(n))
  > y <- f(X)
  > d = ncol(X)
  > 
  > x=seq(0,1,,51)
  > contour(x,x,matrix(f(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
  > points(X)
  > 
  > k <- Kriging(y, X,"gauss",parameters = list(theta=matrix(runif(40),ncol=2)))
  > print(k)
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.254277
  * variance (est.): 0.480222
  * covariance:
    * kernel: gauss
    * range (est.): 0.92171, 0.286891
    * fit:
      * objective: LL
      * optim: BFGS
  > 
  > rlibkriging::save(k, filename="k.h5")
  > 
  > k2 <- rlibkriging::load(filename="k.h5")
  > print(k2)
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.254277
  * variance (est.): 0.480222
  * covariance:
    * kernel: gauss
    * range (est.): 0.92171, 0.286891
    * fit:
      * objective: LL
      * optim: BFGS
  > 
  > test_that("Save/Load NuggetKriging", expect_true( print(k) == print(k2)))
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.254277
  * variance (est.): 0.480222
  * covariance:
    * kernel: gauss
    * range (est.): 0.92171, 0.286891
    * fit:
      * objective: LL
      * optim: BFGS
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.254277
  * variance (est.): 0.480222
  * covariance:
    * kernel: gauss
    * range (est.): 0.92171, 0.286891
    * fit:
      * objective: LL
      * optim: BFGS
  Test passed 🎉
  > 
  > nuk <- NuggetKriging(y, X,"gauss",parameters = list(theta=matrix(runif(40),ncol=2)))
  > print(nuk)
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.254277
  * variance (est.): 0.480222
  * covariance:
    * kernel: gauss
    * range (est.): 0.92171, 0.286891
    * nugget (est.): 0
    * fit:
      * objective: LL
      * optim: BFGS
  > 
  > rlibkriging::save(nuk, filename="nuk.h5")
  > 
  > nuk2 <- rlibkriging::load(filename="nuk.h5")
  > print(nuk2)
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.254277
  * variance (est.): 0.480222
  * covariance:
    * kernel: gauss
    * range (est.): 0.92171, 0.286891
    * nugget (est.): 0
    * fit:
      * objective: LL
      * optim: BFGS
  > 
  > test_that("Save/Load NuggetKriging", expect_true( print(nuk) == print(nuk2)))
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.254277
  * variance (est.): 0.480222
  * covariance:
    * kernel: gauss
    * range (est.): 0.92171, 0.286891
    * nugget (est.): 0
    * fit:
      * objective: LL
      * optim: BFGS
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.254277
  * variance (est.): 0.480222
  * covariance:
    * kernel: gauss
    * range (est.): 0.92171, 0.286891
    * nugget (est.): 0
    * fit:
      * objective: LL
      * optim: BFGS
  Test passed 🎊
  > 
  > nok <- NoiseKriging(y, rep(0.1^2,nrow(X)), X,"gauss",parameters = list(theta=matrix(runif(40),ncol=2)))
  > print(nok)
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.146854
  * variance (est.): 0.305096
  * covariance:
    * kernel: gauss
    * range (est.): 0.749122, 0.241536
    * noise: 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01
    * fit:
      * objective: LL
      * optim: BFGS
  > 
  > rlibkriging::save(nok, filename="nok.h5")
  > 
  > nok2 <- rlibkriging::load(filename="nok.h5")
  > print(nok2)
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.146854
  * variance (est.): 0.305096
  * covariance:
    * kernel: gauss
    * range (est.): 0.749122, 0.241536
    * noise: 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01
    * fit:
      * objective: LL
      * optim: BFGS
  > 
  > test_that("Save/Load NoiseKriging", expect_true( print(nok) == print(nok2)))
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.146854
  * variance (est.): 0.305096
  * covariance:
    * kernel: gauss
    * range (est.): 0.749122, 0.241536
    * noise: 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01
    * fit:
      * objective: LL
      * optim: BFGS
  * data: 20x[0.0420595,0.956833],[0.0246137,0.99427] -> 20x[-0.884835,0.990808]
  * trend constant (est.): 0.146854
  * variance (est.): 0.305096
  * covariance:
    * kernel: gauss
    * range (est.): 0.749122, 0.241536
    * noise: 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01
    * fit:
      * objective: LL
      * optim: BFGS
  Test passed 🎉
  > 
  > proc.time()
     user  system elapsed 
     1.85    0.23    2.09 
Running the tests in 'tests/test-rlibkriging-demo.R' failed.
Complete output:
  > library(testthat)
  >  Sys.setenv('OMP_THREAD_LIMIT'=2)
  >  library(rlibkriging)
  
  Attaching package: 'rlibkriging'
  
  The following objects are masked from 'package:base':
  
      load, save
  
  > 
  > X <- as.matrix(c(0.0, 0.2, 0.5, 0.8, 1.0))
  > f <- function(x) 1 - 1 / 2 * (sin(12 * x) / (1 + x) + 2 * cos(7 * x) * x^5 + 0.7)
  > y <- f(X)
  > 
  > #library(rlibkriging)
  > k_R <- Kriging(y, X, "gauss")
  > 
  > print(k_R)
  * data: 5x[0,1] -> 5x[0.030241,0.772403]
  * trend constant (est.): 0.453087
  * variance (est.): 0.0653536
  * covariance:
    * kernel: gauss
    * range (est.): 0.0400794
    * fit:
      * objective: LL
      * optim: BFGS
  > 
  > x <- as.matrix(seq(0, 1, , 100))
  > p <- predict(k_R, x, TRUE, FALSE)
  > 
  > pdf("rplot.pdf") # plot to file
  > plot(f)
  > points(X, y)
  > 
  > lines(x, p$mean, col = 'blue')
  > polygon(c(x, rev(x)), c(p$mean - 2 * p$stdev, rev(p$mean + 2 * p$stdev)), border = NA, col = rgb(0, 0, 1, 0.2))
  > 
  > s <- simulate(k_R,nsim = 10, seed = 123, x=x)
  > plot(f)
  > points(X,y)
  > matplot(x,s,col=rgb(0,0,1,0.2),type='l',lty=1,add=T)
  > 
  > Xn <- as.matrix(c(0.3,0.4))
  > yn <- f(Xn)
  > print(k_R)
  * data: 5x[0,1] -> 5x[0.030241,0.772403]
  * trend constant (est.): 0.453087
  * variance (est.): 0.0653536
  * covariance:
    * kernel: gauss
    * range (est.): 0.0400794
    * fit:
      * objective: LL
      * optim: BFGS
  > update(k_R, yn, Xn)
  > print(k_R)
  * data: 7x[0,1] -> 7x[0.030241,1.01542]
  * trend constant (est.): 0.50458
  * variance (est.): 0.0979076
  * covariance:
    * kernel: gauss
    * range (est.): 0.117859
    * fit:
      * objective: LL
      * optim: BFGS
  > 
  > proc.time()
     user  system elapsed 
     0.70    0.12    0.79