Running ‘test-AllKrigingConcistency.R’ [8s/10s]
  Running ‘test-KrigingCopy.R’
  Running ‘test-KrigingFit.R’
  Running ‘test-KrigingLeaveOneOut.R’
  Running ‘test-KrigingLeaveOneOut_3d.R’
  Running ‘test-KrigingLogLik.R’
  Running ‘test-KrigingLogLikGradHess.R’ [27s/28s]
  Running ‘test-KrigingMethods.R’
  Running ‘test-KrigingPredict.R’ [25s/29s]
  Running ‘test-KrigingSimulate.R’
  Running ‘test-KrigingUpdate.R’
  Running ‘test-KrigingUpdateSimulate.R’
  Running ‘test-LinearAlgebra.R’
  Running ‘test-NoiseKrigingFit.R’ [12s/14s]
  Running ‘test-NoiseKrigingLogLik.R’
  Running ‘test-NoiseKrigingMethods.R’
  Running ‘test-NoiseKrigingPredict.R’ [17s/21s]
  Running ‘test-NoiseKrigingSimulate.R’ [10s/12s]
  Running ‘test-NoiseKrigingUpdate.R’
  Running ‘test-NoiseKrigingUpdateSimulate.R’
  Running ‘test-NuggetKrigingFit.R’ [15s/16s]
  Running ‘test-NuggetKrigingLogLik.R’ [11s/13s]
  Running ‘test-NuggetKrigingLogMargPost.R’ [18s/21s]
  Running ‘test-NuggetKrigingMethods.R’
  Running ‘test-NuggetKrigingPredict.R’ [17s/22s]
  Running ‘test-NuggetKrigingSimulate.R’
  Running ‘test-NuggetKrigingUpdate.R’
  Running ‘test-NuggetKrigingUpdateSimulate.R’
  Running ‘test-RobustGaSP-Nugget.R’
  Running ‘test-RobustGaSP.R’
  Running ‘test-RobustGaSPtrendlinear.R’
  Running ‘test-RobustGaSPvsKrigingLMP.R’
  Running ‘test-RobustGaSPvsNuggetKrigingLMP.R’
  Running ‘test-SaveLoad.R’
  Running ‘test-asDiceKriging.R’ [22s/26s]
  Running ‘test-estimnone.R’
  Running ‘test-normalize.R’
  Running ‘test-rlibkriging-demo.R’
  Running ‘test-unstableLL.R’
Running the tests in ‘tests/test-RobustGaSP.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(rlibkriging, lib.loc="bindings/R/Rlibs")
  > ##library(testthat)
  > 
  > library(RobustGaSP)
  #########
  ##
  ## Robust Gaussian Stochastic Process, RobustGaSP Package
  ## Copyright (C) 2016-2025 Mengyang Gu, Jesus Palomo and James O. Berger
  #########
  
  Attaching package: 'RobustGaSP'
  
  The following object is masked from 'package:rlibkriging':
  
      simulate
  
  The following object is masked from 'package:stats':
  
      simulate
  
  > 
  > 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)
  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  30 
   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  30 
   The value of the  marginal posterior  function is  1.035387 
   Optimized range parameters are 0.05296527 
   Optimized nugget parameter is 0 
   Convergence:  TRUE 
  > #library(rlibkriging)
  > r <- Kriging(y, X,
  +   kernel="matern5_2",
  +   regmodel = "constant", 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)
  [1] -1.901254
  > 
  > 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)
  [1] -1.901254
  > 
  > lines(seq(0.1,6,,5),lmp_lk(seq(0.1,6,,5)),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 passed 🥳
  > 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))
  Test passed 🥳
  > 
  > 
  > 
  > 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)
  [1] -1.703845
  > 
  > 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)
  [1] -1.703845
  > 
  > 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))
  Test passed 🎉
  > 
  > 
  > # Check predict
  > 
  > ntest <- 10
  > Xtest <- seq(0,1,,ntest)
  > Ytest_rgasp <- predict(k,matrix(Xtest,ncol=1))
  > Ytest_libK <- predict(r,Xtest)
  > 
  > plot(f)
  > points(X,y)
  > lines(Xtest,Ytest_rgasp$mean,col='blue')
  > polygon(c(Xtest,rev(Xtest)),
  +         c(Ytest_rgasp$mean+2*Ytest_rgasp$sd,rev(Ytest_rgasp$mean-2*Ytest_rgasp$sd)),
  +         col=rgb(0,0,1,0.1), border=NA)
  > 
  > lines(Xtest,Ytest_libK$mean,col='red')
  > polygon(c(Xtest,rev(Xtest)),
  +         c(Ytest_libK$mean+2*Ytest_libK$stdev,rev(Ytest_libK$mean-2*Ytest_libK$stdev)),
  +         col=rgb(1,0,0,0.1), border=NA)
  > 
  > precision <- 1e-3
  > test_that(desc=paste0("pred mean is the same that RobustGaSP one"),
  +           expect_equal(predict(r,0.7)$mean[1],predict(k,matrix(0.7))$mean,tol = precision))
  Test passed 😸
  > test_that(desc=paste0("pred sd is the same that RobustGaSP one"),
  +           expect_equal(predict(r,0.7)$stdev[1],predict(k,matrix(0.7))$sd,tol = precision))
  Test passed 🥇
  > 
  > 
  > ## RobustGaSP examples
  > 
  >   #---------------------------------------
  >   # a 1 dimensional example 
  >   #---------------------------------------
  > context("RobustGaSP / 1 dimensional example")
  > 
  > 
  >   input=10*seq(0,1,1/14)
  >   output<-higdon.1.data(input)
  >   #the following code fit a GaSP with zero mean by setting zero.mean="Yes"
  >   model<- rgasp(design = input, response = output, zero.mean="No")
  The upper bounds of the range parameters are 670.0756 
  The initial values of range parameters are 13.40151 
  Start of the optimization  1  : 
  The number of iterations is  30 
   The value of the  marginal posterior  function is  -10.48964 
   Optimized range parameters are 13.2106 
   Optimized nugget parameter is 0 
   Convergence:  TRUE 
  The initial values of range parameters are 0.08888889 
  Start of the optimization  2  : 
  The number of iterations is  30 
   The value of the  marginal posterior  function is  -15.24592 
   Optimized range parameters are 0.1706386 
   Optimized nugget parameter is 0 
   Convergence:  TRUE 
  >   model
  
  Call:
  rgasp(design = input, response = output, zero.mean = "No")
  Mean parameters:  2.174187e-10 
  Variance parameter:  4249.587 
  Range parameters:  13.2106 
  Noise parameter:  0 
  >   
  >   testing_input = as.matrix(seq(0,10,1/100))
  >   model.predict<-predict(model,testing_input)
  >   names(model.predict)
  [1] "mean"    "lower95" "upper95" "sd"     
  >   
  >   #########plot predictive distribution
  >   testing_output=higdon.1.data(testing_input)
  >   plot(testing_input,model.predict$mean,type='l',col='blue',
  +        xlab='input',ylab='output')
  >   polygon( c(testing_input,rev(testing_input)),c(model.predict$lower95,
  +         rev(model.predict$upper95)),col =  "grey80", border = FALSE)
  >   lines(testing_input, testing_output)
  >   lines(testing_input,model.predict$mean,type='l',col='blue')
  >   lines(input, output,type='p')
  >   
  >   ## mean square erros
  >   mean((model.predict$mean-testing_output)^2)
  [1] 4.63608e-05
  > 
  > model_libK = Kriging(matrix(output,ncol=1), matrix(input,ncol=1), 
  +   kernel="matern5_2", 
  +   regmodel = "constant", normalize = FALSE, 
  +   optim = "BFGS", 
  +   objective = "LMP", parameters = NULL)
  > 
  >     lines(testing_input,predict(model_libK,testing_input)$mean,type='l',col='red')
  >     polygon( 
  +       c(testing_input,rev(testing_input)),
  +       c(
  +         predict(model_libK,testing_input)$mean+2*predict(model_libK,testing_input)$stdev,
  +         rev(predict(model_libK,testing_input)$mean-2*predict(model_libK,testing_input)$stdev)),
  +       col = rgb(1,0,0,0.1), border = FALSE)
  > 
  > precision <- 1e-3
  > test_that(desc=paste0("RobustGaSP / 1 dimensional example / pred mean is the same that RobustGaSP one"),
  +           expect_equal(predict(model_libK,0.7)$mean[1],predict(model,matrix(0.7))$mean,tol = precision))
  ── Failure: RobustGaSP / 1 dimensional example / pred mean is the same that RobustGaSP one ──
  predict(model_libK, 0.7)$mean[1] not equal to predict(model, matrix(0.7))$mean.
  1/1 mismatches
  [1] 0.621 - 0.623 == -0.00162
  
  Error:
  ! Test failed
  Backtrace:
      ▆
   1. ├─testthat::test_that(...)
   2. │ └─withr (local) `<fn>`()
   3. └─reporter$stop_if_needed()
   4.   └─rlang::abort("Test failed", call = NULL)
  Execution halted