R Under development (unstable) (2024-08-21 r87038 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> 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,return_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,return_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.25    0.40    2.64