# Multi-block Parameter Calibration 

# #########################################################################
# Description: L2 Calibration 
# Date: 2021.7.20
# #########################################################################

rm(list=ls(all=TRUE))

library(CVST)
library(kernlab)
library(lhs)
library(laGP)

# #########################################################################
# 1. Problem instance
# #########################################################################

# true process
zeta <- function(x) {
  z = exp(x/10) + ifelse(x<pi, exp(x/10)*sin(x), exp(x/10)*cos(x))  
  return(z) 
}

# computer model 
eta <- function(x,theta) {
  e = zeta(x) - 5*sqrt(theta[1]^2-theta[1]+1)*(sin(theta[1]*x/10)+cos(theta[1]*x/10)) - ifelse(x<pi, (1/2)*sqrt(theta[2]^2-theta[2]+1)*(sin(theta[2]*x/5)+cos(theta[2]*x/5)), (1/2)*sqrt(theta[3]^2-theta[3]+1)*(sin(theta[3]*x/10)+cos(theta[3]*x/10)))
  return(e)
}

# #########################################################################
# 2. Algorithm
# #########################################################################

niter = 200
dim_param = 3 
start = c(-1,-7,7,-1,-7,9,-1,-5,7,-1,-5,9,1,-7,7,1,-7,9,1,-5,7,1,-5,9)
start = matrix(start,8,dim_param,byrow=TRUE)
res_param_mean = matrix(0,nrow(start),dim_param)
res_param_sd = matrix(0,nrow(start),dim_param)
res_mse_mean = matrix(0,nrow(start),dim_param)
res_mse_sd = matrix(0,nrow(start),dim_param)
eps <- sqrt(.Machine$double.eps)

tic <- proc.time()[3]
for (j in 3:3 ) { #nrow(start)

  start_param = matrix(0,niter,dim_param)
  start_mse = matrix(0,niter,dim_param)
  
  for (iter in 1:niter) {
    
    cat("starting point:", j,"\n")
    cat("iteration number:", iter,"\n")
    
    # training sets
    set.seed(iter)
    sgma <- 0.1
    n <- 200
    x <- (1:n) * (2*pi)/n #sort(runif(n)*2*pi)
    y <- zeta(x)+sgma*rnorm(n)
    theta = randomLHS(n,dim_param)

    theta[,1] = (theta[,1]-0.5)*6
    theta[,2] = (theta[,2]-0.5)*6-6
    theta[,3] = (theta[,3]-0.5)*6+8
    
    # training data for emulation 
    x_emulate = cbind(x,theta)
    y_emulate = numeric(0)
    for (i in 1:n) {
      y_emulate[i] = eta(x[i],theta[i,])
    }
    
    # gaussian process 
    gpi <- newGPsep(x_emulate, y_emulate, d=0.1, g=0.1*var(y_emulate), dK=TRUE)
    mle <- mleGPsep(gpi, param="both", tmin=c(eps, eps), tmax=c(10, var(y_emulate)))
   
    # test sets
    xtest = sort(runif(n)*2*pi)
    ytest <- zeta(xtest)+sgma*rnorm(n)
    
    idx_test = (xtest<pi)
    xtest2 = xtest[idx_test]
    xtest3 = xtest[!idx_test]
  
    krr <- constructKRRLearner()
    
    # traning and test sets
    dat <- constructData(x,y)
    dat_tst <- constructData(xtest,0)
    
    # tuning parameter 
    lambdas = 10^(-8:2)
    sigmas = 10^((-8:8)/3)
    
    params <- constructParams(kernel="rbfdot", sigma=sigmas, lambda=lambdas)
    opt <- CV(dat, krr, params, fold=5, verbose=FALSE)
    
    param <- list(kernel="rbfdot", sigma=opt[[1]]$sigma, lambda=opt[[1]]$lambda)
    krr.model <- krr$learn(dat,param)
    pred <- krr$predict(krr.model,dat)
    plot(x,y, main=paste("selected values: sigma=",signif(param$sigma,digits=3),", lambda=",signif(param$lambda,digits=3)))
    lines(x,pred,col='blue')
    
    l2curve <- function(theta.)
    {
      rslt <- sqrt( sum( (as.vector(pred) - predGPsep(gpi, cbind(x,t(replicate(n,theta.))))$mean  )^2 ) )
      return(rslt)
    }
   
    # optim for theta 
    start_param[iter,] = optim(start[j,],l2curve)$par

    start_mse[iter,1] = (1/n)*sum( (ytest - eta(x=xtest,start_param[iter,]))^2 )
    start_mse[iter,2] = (1/sum(idx_test))*sum( (ytest[idx_test] - eta(x=xtest2,start_param[iter,]))^2 )
    start_mse[iter,3] = (1/sum(!idx_test))*sum( (ytest[!idx_test] - eta(x=xtest3,start_param[iter,]))^2 )

    deleteGPsep(gpi)
  }
  
  mean_param = apply(start_param,2,mean)
  sd_param = apply(start_param,2,sd)
  mean_mse = apply(start_mse,2,mean)
  sd_mse = apply(start_mse,2,sd)
  
  res_param_mean[j,] = mean_param
  res_param_sd[j,] = sd_param
  res_mse_mean[j,] = mean_mse
  res_mse_sd[j,] = sd_mse
}
toc <- proc.time()[3]
time = toc - tic 
idx = which.min(res_mse_mean[,1])

cat("calibrate parameter value (mean): ", res_param_mean[idx,], "\n")
cat("calibrate parameter value (sd): ", res_param_sd[idx,], "\n")
cat("mse (mean): ", res_mse_mean[idx,], "\n")
cat("mse (sd): ", res_mse_sd[idx,], "\n")

# #########################################################################
# 3. Save file
# #########################################################################

# save.image(file = paste0(rslt_path,"/L2_result.RData"))
# load(file = paste0(rslt_path,"/L2_result.RData"))

save.image(file = paste0("/root/capsule/results/L2_result.RData"))
# load(file = paste0("/root/capsule/results/L2_result.RData"))

