R Under development (unstable) (2024-08-23 r87049 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

> 
> ##library(rlibkriging, lib.loc="bindings/R/Rlibs")
> ##library(testthat)
> 
> # 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.01)
> d = ncol(X)
> 
> x=seq(0,1,,5)
> contour(x,x,matrix(f(as.matrix(expand.grid(x,x))),nrow=length(x)),nlevels = 30)
> points(X)
> 
> r <- NuggetKriging(y, X,"gauss",parameters = list(theta=matrix(runif(40),ncol=2),nugget=0,sigma2=1))
> l= as.list(r)
> # ll = function(X) {
> #   logLikelihoodFun(r,X,return_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,return_grad=TRUE)$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("nugget / print")
> 
> p = capture.output(print(r))
> 
> 
> context("nugget / logLikelihood")
> 
> ll = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r,c(theta,l$sigma2/(l$nugget+l$sigma2)))$logLikelihood)}
> t=seq(0.01,2,,5)
> 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("nugget / logLikelihoodFun returned",
+           expect_equal(names(logLikelihoodFun(r,runif(d+1))),c("logLikelihood")))
Test passed 😸
> test_that("nugget / logLikelihoodFun logLikelihoodGrad returned",
+           expect_equal(names(logLikelihoodFun(r,runif(d+1),return_grad=TRUE)),c("logLikelihood","logLikelihoodGrad")))
Test passed 🎊
> 
> test_that("nugget / logLikelihoodFun dim",
+           expect_equal(dim(logLikelihoodFun(r,rbind(runif(d+1),runif(d+1)))$logLikelihood),c(2,1)))
Test passed 🥳
> test_that("nugget / logLikelihoodGrad dim",
+           expect_equal(dim(logLikelihoodFun(r,rbind(runif(d+1),runif(d+1)),return_grad=TRUE)$logLikelihoodGrad),c(2,d+1)))
Test passed 🌈
> 
> 
> context("nugget / predict")
> 
> test_that("nugget / predict mean stdev returned",
+           expect_equal(names(predict(r,runif(d))),c("mean","stdev")))
Test passed 🎉
> test_that("nugget / predict mean returned",
+           expect_equal(names(predict(r,runif(d),return_stdev=F)),c("mean")))
Test passed 😸
> test_that("nugget / predict mean stdev cov returned",
+           expect_equal(names(predict(r,runif(d),return_cov=TRUE)),c("mean","stdev","cov")))
Test passed 😸
> 
> test_that("nugget / predict mean dim",
+           expect_equal(dim(predict(r,rbind(runif(d),runif(d)))$mean),c(2,1)))
Test passed 🎊
> test_that("nugget / predict stdev dim",
+           expect_equal(dim(predict(r,rbind(runif(d),runif(d)))$stdev),c(2,1)))
Test passed 🎊
> test_that("nugget / predict cov dim",
+           expect_equal(dim(predict(r,rbind(runif(d),runif(d)),return_cov=TRUE)$cov),c(2,2)))
Test passed 🎊
> 
> 
> context("nugget / simulate")
> 
> test_that("nugget / simulate dim",
+           expect_equal(dim(simulate(r,x=rbind(runif(d),runif(d)))),c(2,1)))
Test passed 🥇
> test_that("nugget / simulate nsim dim",
+           expect_equal(dim(simulate(r,x=rbind(runif(d),runif(d)),nsim=10)),c(2,10)))
Test passed 🎉
> 
> # offset is not legitimate because nugget kriging is interpolating
> #test_that("nugget / simulate mean X",
> #          expect_equal(mean(simulate(r,nsim = 100, x=X[1,]+0.00005)),y[1],tolerance = 0.01))
> test_that("nugget / simulate mean X",
+           expect_equal(mean(simulate(r,nsim = 100, x=X[1,])),y[1],tolerance = 0.01))
Test passed 🌈
> set.seed(12345)
> x = runif(d)
> test_that("nugget / simulate mean",
+           expect_equal(mean(simulate(r,nsim = 10000, x=x)),predict(r,x)$mean[1],tolerance = 0.01))
Test passed 🥳
> test_that("nugget / simulate sd",
+           expect_equal(sd(simulate(r,nsim = 10000, x=x)),predict(r,x)$stdev[1],tolerance = 0.01))
Test passed 🌈
> 
> 
> context("nugget / update")
> 
> set.seed(1234)
> X2 = matrix(runif(d*10),ncol=d)
> y2 = f(X2) + rnorm(nrow(X2),0,0.01)
> x=seq(0,1,,5)
> 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,return_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,return_grad=TRUE)$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,10,,5)
> contour(t,t,matrix(ll(as.matrix(expand.grid(t,t))),nrow=length(t)),xlim=c(0,10),ylim=c(0,1),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 <- NuggetKriging(c(y,y2), rbind(X,X2),"gauss", parameters = list(theta=matrix(as.list(r)$theta,ncol=2),nugget=0,sigma2=1))
> ll2 = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r2,c(theta,as.list(r2)$sigma2/(as.list(r2)$sigma2+as.list(r2)$nugget)))$logLikelihood)}
> contour(t,t,matrix(ll2(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30,add=TRUE,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,X2)
> llu = function(Theta){apply(Theta,1,function(theta) logLikelihoodFun(r,  c(theta,as.list(r2)$sigma2/(as.list(r2)$sigma2+as.list(r2)$nugget)) )$logLikelihood)}
> contour(t,t,matrix(llu(as.matrix(expand.grid(t,t))),nrow=length(t)),nlevels = 30,add=TRUE,col='blue')
> points(as.list(r)$theta[1],as.list(r)$theta[2],col='blue',pch=20)
> 
> pu = capture.output(print(r))
> 
> test_that("nugget / 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$nugget())),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$nugget())),TRUE) ))),"\n")
> 
> test_that("nugget / update almost converge",
+           expect_equal(as.list(r2)$theta, as.list(r)$theta, tolerance = 2E-2))
Test passed 🥳
> 
> proc.time()
   user  system elapsed 
   3.40    0.31    3.65