unstable LL when theta to high about libkriging HOT 2 CLOSED

# remotes::install_url("https://github.com/libKriging/rlibkriging/releases/download/0.8-0.1/rlibkriging_0.8-0_R_x86_64-pc-linux-gnu.tar.gz", dependencies = TRUE, build=FALSE)

library(rlibkriging)

# 1D function, small design
X <- as.matrix(c(0.0, 0.25, 0.33, 0.45, 0.5, 0.75, 1.0))
f <- function(x) 1 - 1 / 2 * (sin(4 * x) / (1 + x) + 2 * cos(12 * x) * x^6 + 0.7)
y <- f(X)  #+ 0.5*rnorm(nrow(X))

plot(f)
points(X, y)

rlibkriging:::optim_log(4)
# Build Kriging (https://libkriging.readthedocs.io/en/latest/math/KrigingModels.html)
k <- Kriging(y, X, kernel="gauss", optim="BFGS")#, parameters=list(sigma2=var(y), is_sigma2_estim=FALSE))#, objective = "LL", optim = "BFGS", regmodel = "constant")
# kernel: https://libkriging.readthedocs.io/en/latest/math/kernel.html
# regmodel: https://libkriging.readthedocs.io/en/latest/math/trend.html
# parameters: https://libkriging.readthedocs.io/en/latest/math/parameters.html
print(k)

# # Predict
# .x <- sort(c(as.matrix(seq(0, 1, , 101)),X))
# p.x <- predict(k, .x, TRUE, FALSE)
# 
# plot(f)
# points(X, y)
# lines(.x, p.x$mean, col = 'blue')
# polygon(c(.x, rev(.x)), c(p.x$mean - 2 * p.x$stdev, rev(p.x$mean + 2 * p.x$stdev)), border = NA, col = rgb(0, 0, 1, 0.2))
# 
# # log-Likelihood (https://libkriging.readthedocs.io/en/latest/math/likelihood.html)
# print(k$logLikelihood())
# 
# plot( function(t) k$logLikelihoodFun(t)$logLikelihood ,xlab=expression(theta),ylab="LL", xlim=c(0.001,10))
# abline(v=k$theta(),col='blue')
# abline(h=k$logLikelihood(),col='blue')

getR = function(k) {
  # Save / load /inspect
  save(k,"k.h5")
  # install.packages("hdf5r")
  k.h5 = hdf5r::h5file("k.h5")
  # print(k.h5)
  T = k.h5[['T']]$read()
  R = (T) %*% t(T)
  return(R)
}

Rcpp::sourceCpp('rcond.cpp') #library(Matrix)

# leave-one-out (https://libkriging.readthedocs.io/en/latest/math/leaveOneOut.html)
#print(k$leaveOneOut())
plot( function(t) k$logLikelihoodFun(t)$logLikelihood ,xlab=expression(theta),ylab="LL", xlim=c(0.00,10))

for (t in seq(0.00,10,,101)) {
  k_t <- Kriging(y, X, kernel="gauss", optim="none", parameters=list(theta=matrix(t), is_theta_estim=FALSE))
  R_t = getR(k_t)
  if (rcond(R_t) > nrow(R_t) * 1e-13) {
    points(t,k$logLikelihoodFun(t)$logLikelihood,col='blue')
  } else {
    points(t,k$logLikelihoodFun(t)$logLikelihood,col='red')
  }
}

##  # plot( function(t) k$leaveOneOutFun(t)$leaveOneOut ,xlab=expression(theta),ylab="LOO", xlim=c(0.01,10))
##  # abline(v=k$theta(),col='blue')
##  # abline(h=k$leaveOneOut(),col='blue')
##  
##  ## k <- Kriging(y, X, kernel="gauss", optim="none", parameters=list(theta=matrix(9.2), is_theta_estim=FALSE))#, objective = "LL", optim = "BFGS", regmodel = "constant")
##  ## 
##  
##  ## kappa(R)
##  ## library(RcppArmadillo)
##  ## Rcpp::sourceCpp("inv_sympd.cpp")
##  ## inv_sympd(R)
##  
##  
 plot( function(t) k$leaveOneOutFun(t)$leaveOneOut ,xlab=expression(theta),ylab="LOO", xlim=c(0.00,10))

 for (t in seq(0.00,10,,101)) {
   k_t <- Kriging(y, X, kernel="gauss", optim="none", parameters=list(theta=matrix(t), is_theta_estim=FALSE))
   R_t = getR(k_t)
   if (rcond(R_t) > nrow(R_t) * 1e-13) {
     points(t,k$leaveOneOutFun(t)$leaveOneOut,col='blue')
   } else {
     points(t,k$leaveOneOutFun(t)$leaveOneOut,col='red')
   }
 }

 
 plot( function(t) k$logMargPostFun(t)$logMargPost ,xlab=expression(theta),ylab="LMP", xlim=c(0.00,10))
 
 for (t in seq(0.00,10,,101)) {
   k_t <- Kriging(y, X, kernel="gauss", optim="none", objective="LMP",parameters=list(theta=matrix(t), is_theta_estim=FALSE))
   R_t = getR(k_t)
   if (rcond(R_t) > nrow(R_t) * 1e-13) {
     points(t,k$logMargPostFun(t)$logMargPost,col='blue')
   } else {
     points(t,k$logMargPostFun(t)$logMargPost,col='red')
   }
 }