library(DiceKriging)

N <- 200 #nb of points in the trajectory
x <- seq(0,1,length.out=N)
dist_all <- as.matrix(dist(x))
n <- 10 # simulation of n trajectories
Z1 <- matrix(rnorm(N*n),N,n)  # simulation of a gaussian centred vector 
# with covariance matrix identity
f0 = rep(0,N)


# gaussian covariance 
sig2 <- 0.5
theta <- 0.05
## theta <- 0.2
lambda <- 10^(-10)
cov = sig2*(exp(-(dist_all/theta)^2) + lambda * diag(1,N))  # cov between x
Z2 <- matrix(0,N,n)
for (i in 1:n){
	Z2[,i] <- f0+t(chol(cov))%*%Z1[,i]
}

plot(x,Z2[,1], main=paste('Simulation of GP trajectories 
with gaussian covariance and theta =',theta),type='l',ylim=c(min(Z2),max(Z2)))
for (i in 2:n)
	lines(x,Z2[,i],col=i)

# exponential covariance
sig2 <- 0.5
theta <- 0.05
cov <- sig2*exp(-(dist_all/theta))
Z2 <- Z1
for (i in 1:n)
  Z2[,i] <- f0+t(chol(cov))%*%Z1[,i]
plot(x,Z2[,1], main=paste('Simulation of GP trajectories with exponential covariance and theta =',theta),type='l',ylim=c(min(Z2),max(Z2)))
for (i in 2:n)
	lines(x,Z2[,i],col=i)

# gaussian covariance with nugget effect
sig2 <- 0.5
theta <- 0.2
lambda <- 0.05
cov <- sig2*(exp(-(dist_all/theta)^2) + lambda * diag(1,N))
Z2 <- Z1
for (i in 1:n)
  Z2[,i] <- f0+t(chol(cov))%*%Z1[,i]
plot(x,Z2[,1], main=paste('Simulation of GP trajectories with gaussian covariance and theta =',theta),type='l',ylim=c(min(Z2),max(Z2)))
for (i in 2:n)
  lines(x,Z2[,i],col=i)

##### or with DiceKriging
type <- "exp"
coef <- c(theta = 0.05)
sigma <- 0.5
model <- km(design=x, response=rep(0,N),coef.trend = 0,covtype=type, coef.cov=coef, coef.var=sigma)
Z2 <- simulate(model, nsim=10, newdata=NULL)
plot(x,Z2[1,], main=paste('Simulation of GP trajectories with exponential covariance and theta =',theta),type='l',ylim=c(min(Z2),max(Z2)))
for (i in 2:n)
  lines(x,Z2[i,], main=paste('Simulation of GP trajectories with exponential covariance and theta =',theta))


type <- "gauss"
lambda <- 10^(-10) # diagonal term to outperform the conditioning of the covariance matrix 
model <- km(design=x, response=rep(0,N),coef.trend = 0,covtype=type, coef.cov=coef, coef.var=sigma,nugget=lambda)
Z2 <- simulate(model, nsim=10, newdata=NULL)
plot(x,Z2[1,], main=paste('Simulation of GP trajectories with exponential covariance and theta =',theta),type='l',ylim=c(min(Z2),max(Z2)))
for (i in 2:n)
  lines(x,Z2[i,], main=paste('Simulation of GP trajectories with exponential covariance and theta =',theta))

d = 2 # inputs space dimension
N_dim = 50 # nb of points in each dimension of the inputs space
N_tot = N_dim^d 
x <- as.matrix(expand.grid(seq(0,1,length.out=N_dim),seq(0,1,length.out=N_dim)))

mu <- 0
sig2 <- 0.5
dist1 <- as.matrix(dist( x[,1] ))
dist2 <- as.matrix(dist( x[,2] ))
# isotropic Gaussian covariance
theta = 0.1
lambda <- 10^(-10)  # diagonal term to outperform the conditioning of the covariance matrix
cov = sig2*(exp(-((dist1/theta)^2 +(dist2/theta)^2))+ lambda * diag(1,N_tot))
f0 = x %*% c(0,0)

n <- 1 # simulation of n trajectories
Z1 <- matrix(rnorm(N_tot*n),N_tot,n) # simulation of a centred gaussian vector with covariance matrix identity 
Z2=f0+t(chol(cov))%*%Z1  # then linear transformation

par(mfrow=c(1,1))
filled.contour(matrix(Z2,N_dim,N_dim),main=paste('Simulation trajectoires PG de covariance gaussienne isotrope avec theta =', theta),color.palette=heat.colors)

# anisotropic gaussian covariance
theta1 = 0.1
theta2 = 0.03
cov = sig2*(exp(-((dist1/theta1)^2 +(dist2/theta2)^2)) + lambda * diag(1,N_tot))
f0 = x %*% c(0,0)

Z2=f0+t(chol(cov))%*%Z1  # then linear transformation 
filled.contour(matrix(Z2,N_dim,N_dim),main=paste('Simulation of GP trajectories with anisotropic gaussian covariance and theta1 =', theta1,'& theta2=',theta2),color.palette=heat.colors)

# anisotropic exponential covariance 
cov = sig2*exp(-(dist1/theta1 + dist2/theta2)) 
f0 = x %*% c(0,0)

Z2=f0+t(chol(cov))%*%Z1  # then linear transformation 
filled.contour(matrix(Z2,N_dim,N_dim),main=paste('Simulation of GP trajectories with anisotropic exponential covariance and theta1 =', theta1,'& theta2=',theta2),color.palette=heat.colors)

rm(list=ls(all.names=TRUE))
library(lhs)
library(DiceKriging)
library(DiceView)

######################## Learning basis and test basis ##########################
N_BA = 15                ## Number of points in the learning sample (BA)
d = 1                     ## Dimension of the inputs space
S_BA = optimumLHS(N_BA,d) ## Simulation os the learning basis BA by optimized LHS  => use the function optimumLHS
N_BT = 100                 ### Number of points in the test basis
S_BT = seq(0,1,length.out=N_BT)    ### Simulation of an initial test basis
                           ## S_BT : inputs for the test basis
                           ## Y_BT : output of the emulator

fx_exo2 <- function(X)
sin(30*(X-0.9)^4)*cos(2*(X-0.9)) + (X-0.9)/2

Y_BA <- fx_exo2(S_BA)   ## Y_BA : output of the emulator on the learning basis BA
Y_BT <- fx_exo2(S_BT)   ## Y_BT :output of the emulator on the test basis BT

plot(S_BT,Y_BT,type='l')
points(S_BA,Y_BA,pch=3,col=2)
legend("topright",y=c(0.17,0.33),c('Theoretical function','Pts in the learning basis'),lty=c(1,0),
pch=c(-1,3),col=c(1,2),cex=0.8)


library(DiceKriging)
metamodel <- km(formula=~1,design=S_BA,response=Y_BA,covtype="matern5_2") # Construction of the Gaussian Process emulator 
prediction <- predict.km(metamodel,S_BT,'UK',checkNames=FALSE) # Prediction of the metamodel on the test basis => use of the function predict.km

par(mfrow=c(2,1))
plot(S_BT,Y_BT,type='l',main='Estimation of f by the GP metamodel')
points(S_BA,Y_BA,pch=3,col='red')
lines(S_BT,prediction$mean,col='blue')
lines(S_BT,prediction$lower95,col='blue',lty=2)
lines(S_BT,prediction$upper95,col='blue',lty=2)
plot(S_BT,prediction$sd,type='l',main="sqrt(MSE) of the estimation")


library(DiceView)
# with DiceView
par(mfrow=c(2,1))
sectionview(metamodel)
plot(S_BT,prediction$sd,type='l',main="sqrt(MSE) of the estimation")

#### Computation of Q2 on BT ####
# Building of the function allowing to compute Q2
Q2 <- function(y,yhat)
  1-mean((y-yhat)^2)/var(y)
Q2_BT <- Q2(Y_BT,prediction$mean) # Computation of Q2 on the test basis


#### Display of the results ####
print(paste('Number of points in the learning sample :',N_BA))
print(paste('Inputs space dimension :',d))
print(paste('Q2 on the test basis :',Q2_BT))


#### Optional part ####
## Computation of the likelihood and optimization of hyperparameters
x <- seq(0.01,1,length=1000)
y <- sapply(x,logLikFun,model=metamodel)
plot(x,y,type='l',main='Log-likelihood with respect to theta')

theta_optim <- x[which.max(y)] # Identification of optimal theta 


theta_optim

## Step 3 : adaptive construction of the design of experiments
par(mfrow=c(2,1))
sectionview(metamodel,xlim=c(0,1))
plot(S_BT,prediction$sd,type='l',main="sqrt(MSE) of the estimation")

## Adaptive design : one iteration
maxi <- which.max(prediction$sd)  # Research of the point where MSE is maximal
points(S_BT[maxi],prediction$sd[maxi],pch=16,col=3)
S_BA2 <- rbind(S_BA,S_BT[maxi])   # add the new point to the learning basis BA
Y_BA2 <- c(Y_BA,fx_exo2(S_BT[maxi])) # add the corresponding output Y 
metamodel2 <- km(formula=~1,design=S_BA2,response=Y_BA2,covtype="matern5_2") # Update the GP metamodel
prediction2 <- predict.km(metamodel2,S_BT,'UK',checkNames=FALSE) # Update the predictions and the MSE of the metamodel

# Display of the results
par(mfrow=c(2,1))
sectionview(metamodel2)
plot(S_BT,prediction2$sd,type='l',main="sqrt(MSE) of the estimation when adding one point")
points(S_BT[maxi],prediction2$sd[maxi],pch=16,col=3)

## Adaptive design : Automation of the process
nIte <- 15
for (i in 1:nIte){ 
	maxi <- which.max(prediction2$sd)   # Research of the point where MSE is maximal
	S_BA2 <- rbind(S_BA2,S_BT[maxi])    # add the new point to the learning basis BA
	Y_BA2 <- c(Y_BA2,fx_exo2(S_BT[maxi]))
	capture.output(metamodel2 <- km(formula=~1,design=S_BA2,response=Y_BA2,covtype="matern5_2")) # Update the GP metamodel
	prediction2 <- predict.km(metamodel2,S_BT,'UK',checkNames=FALSE) # Update predictions and MSE 
  
  # Affichage
  par(mfrow=c(2,1))
  sectionview(metamodel2)
  plot(S_BT,prediction2$sd,type='l',main="sqrt(MSE) wahen adding one point")
  points(S_BT[maxi],prediction2$sd[maxi],pch=16,col=3)
  Sys.sleep(2)
}
Q2_BTnew <- Q2(Y_BT,prediction2$mean)

library(lhs)

######################## learning basis and test basis ##########################
N_BA = 80;                ## Number of points in the learning basis
d = 2;     

      
                   ## Inputs space dimension 
S_BA = optimumLHS(N_BA,d);     ## Construction of the learning basis


                           ## Y_BA : outputs of the emulator
N_BT_dim = 70;
N_BT = N_BT_dim^2;               ### Number of points in the test basis
S_BT_dim <- seq(0,1,length.out=N_BT_dim)
S_BT <- expand.grid(S_BT_dim,S_BT_dim) ### Construction of an initial test basis
                           ## S_BT : inputs in the test basis
                           ## Y_BT : corresponding emulator's outputs


 schwefel <- function(x){
x2 <- x*400-200
rowSums(-x*sin(sqrt(abs(x2))))
}

Y_BA <- schwefel(S_BA)
Y_BT <- schwefel(S_BT)

par(mfrow=c(1,1))
filled.contour(S_BT_dim,S_BT_dim,matrix(Y_BT,N_BT_dim,N_BT_dim),main='Function Schwefel',xlab='X1',ylab='X2',color.palette=heat.colors,
plot.axes={points(S_BA[,1],S_BA[,2],pch=3)})

library(DiceKriging)
metamodel <- km(design=S_BA,response= Y_BA, covtype="matern3_2")

# ###### Prediction on BT with GP metamodel  #######
predictions_BT<- predict(metamodel,S_BT,'UK',checkNames=F)

filled.contour(S_BT_dim,S_BT_dim,matrix(predictions_BT$mean,N_BT_dim,N_BT_dim),xlab='X1',ylab='X2',main='Prediction du metamodele',color.palette=heat.colors) 
filled.contour(S_BT_dim,S_BT_dim,matrix(abs(Y_BT - predictions_BT$mean),N_BT_dim,N_BT_dim),xlab='X1',ylab='X2',main='Erreur du metamodele PG (valeur absolue)',color.palette=heat.colors) 


# ###### Computation of Q2 on BT #########
Q2 <- function(y,yhat)
  1-mean((y-yhat)^2)/var(y)
Q2_BT = Q2(Y_BT,predictions_BT$mean) # Computation of Q2 on the test basis 

#### Results' display ####
print(paste('Number of points in the test basis:',N_BA))
print(paste('Inputs space dimension:',d))
print(paste('Q2 on the test basis:',Q2_BT))








library(lhs)



#### Learning basis and test basis ####

N_BA <- 150;                 ## Number of points in the learning basis
d <- 3;                     ## Inputs space dimension
S_BA <- randomLHS(N_BA,d)   ## Construction of the learning basis
                            ## S_BA : inputs for the learning basis
                            ## Y_BA : corresponding outputs of the emulator
S_BT <- randomLHS(100,3)    ## Construction of the test basis


ishigami <- function(X,coeff){
  ## Analytical Function of Ishigami for scaled inputs => inputs in [0:1]
	## Ishigami function : Y = sinX1+a*(sinX2).^2 + b*X3.^4*sinX1 for Xi ~ U[-pi;+pi]
	## with a = coeff(1); b=coeff(2) 
	## Generally : coeff = [7 0.1];
	X <- 2*pi*X-pi
	sin(X[,1]) + coeff[1]*sin(X[,2])^2 + coeff[2]*X[,3]^4*sin(X[,1])
}

coeff_fx_test = c(7, 0.1) 
Y_BA <- ishigami(S_BA,coeff_fx_test)
Y_BT <- ishigami(S_BT,coeff_fx_test)

#### METAMODELING ####
library(DiceKriging)
metamodel <- km(formula=~1,design=S_BA,response=Y_BA,covtype="matern5_2")
plot(metamodel) # Affichage des resultats

## Computation of Q2 on the test basis
Q2 <- function(y,yhat)
  1-mean((y-yhat)^2)/var(y)

prediction <- predict(metamodel, S_BT, 'UK', checkNames=FALSE)$mean
Q2_BT<-Q2(Y_BT,prediction)
print(paste('Q2 sur Base de test :',Q2_BT))

## Computation of Q2 by cross validation
prediction_LOO <- leaveOneOut.km(metamodel,'UK')$mean
Q2_LOO<- Q2(Y_BA,prediction_LOO)
print(paste('Q2 sur Base de test :',Q2_LOO))


library(sensitivity)
library(boot)
library(numbers)
#### Estimation of first-order and total Sobol' indices ####
# Building of the gaussian process metamodel used to compute Sobol' indices
f <- function(x) 
  predict(metamodel,x,type='UK',checkNames=FALSE)$mean

## Computation of Sobol' indices by SobolEff Janon et al. (2014)
n <- 10000
sample1 <- data.frame(X1=runif(n,0,1),X2=runif(n,0,1),X3=runif(n,0,1))
sample2 <- data.frame(X1=runif(n,0,1),X2=runif(n,0,1),X3=runif(n,0,1))
indices.sobol <- sobolEff(f,sample1,sample2,order=1)
indices.sobol.totaux <- sobolEff(f,sample1,sample2,order=0)
indices.sobol
indices.sobol.totaux
## Computation of first- and closed second-order Sobol' indices with replicated designs
indices.sobol.rep.o1 <- sobolroalhs(model = f, factors = 3, N = 10000, order = 1, nboot=100)
print(indices.sobol.rep.o1)
plot(indices.sobol.rep.o1)
indices.sobol.rep.o2 <- sobolroalhs(model = f, factors = 3, N = 10000, order = 2, nboot=100)
print(indices.sobol.rep.o2)
plot(indices.sobol.rep.o2)