# Examples from slides8
# For STAT8810
# Fall, 2017
# M.T. Pratola



# Example - Minimax Distance Design
library(fields)
cands=as.matrix(expand.grid(seq(0,1,length=10),seq(0,1,length=10)))
nd=9
design=cover.design(cands,nd,nruns=10)$design
plot(design,pch=20,xlab="x1",ylab="x2")
points(cands,col="grey")



# Latin Hypercube Designs (LHS)
library(lhs)
set.seed(66) # only to replicate this output
n=9
p=2
design1=randomLHS(n,p)  # default algorithm
set.seed(66)
design2=optimumLHS(n,p) # maximize mean distance 
                        # between design points
set.seed(66)
design3=maximinLHS(n,p) # maximize the min distance 
                        # between design points

plot(design1,pch='1',xlab="x1",ylab="x2",xlim=c(0,1),ylim=c(0,1))
points(design2,pch='2',col="red")
points(design3,pch='3',col="blue")






# Sequential Design using Expected Improvement
library(DiceOptim)
library(rgl)

# get our initial starting design
#set.seed(7)
#design=optimumLHS(9,2)
design=as.matrix(expand.grid(seq(0,1,length=3),seq(0,1,length=3)))
colnames(design)=c("x1","x2")
y.branin=apply(design,1,branin)

# Fit the GP model - built into the DiceOptim package
fit.gp=km(~1, design = data.frame(x = design), 
          response = y.branin,covtype = "gauss")

# Plot surface
X=as.matrix(expand.grid(seq(0,1,length=10),
            seq(0,1,length=10)))
colnames(X)=c("x1","x2")
yhat=predict(fit.gp, newdata = data.frame(x = X), 
             type = "UK")
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$mean,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$lower95),col="blue",add=TRUE)
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$upper95),col="blue",add=TRUE)
plot3d(design[,1],design[,2],y.branin,type="s",
        radius=7,col="red",add=TRUE)


# Calculate EI
ego=apply(as.matrix(X),1,EI,fit.gp,type="UK",
          minimization=TRUE)
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(ego,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")


# Get next x that maximizes the EI
x.new=max_EI(fit.gp,lower=rep(0,2),upper=rep(1,2),
             parinit = 0.5,minimization=TRUE)$par
yhat.xnew=predict(fit.gp,newdata=x.new,type="UK")$mean


# Calculate EI
ego=apply(as.matrix(X),1,EI,fit.gp,type="UK",
          minimization=TRUE)
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(ego,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")


# Get next x that maximizes the EI
x.new=max_EI(fit.gp,lower=rep(0,2),upper=rep(1,2),
             parinit = 0.5,minimization=TRUE)$par
yhat.xnew=predict(fit.gp,newdata=x.new,type="UK")$mean

persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$mean,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")
plot3d(design[,1],design[,2],y.branin,type="s",
       radius=7,col="red",add=TRUE)
plot3d(x.new[1],x.new[2],yhat.xnew,type="s",
       radius=7,col="green",add=TRUE)



# Update by evaluating our expensive function
y.new=apply(x.new,1,branin)
y.branin=c(y.branin,y.new)
design=rbind(design,x.new)


# Refit the GP model
fit.gp=km(~1, design = data.frame(x = design), 
          response = y.branin, covtype = "gauss")
yhat=predict(fit.gp, newdata = data.frame(x = X), 
          type = "UK")


# Plot
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$mean,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")
plot3d(design[,1],design[,2],y.branin,type="s",
        radius=7,col="red",add=TRUE)
plot3d(x.new[1],x.new[2],yhat.xnew,type="s",
        radius=7,col="green",add=TRUE)


# Plot with uncertainties
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$mean,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")
plot3d(design[,1],design[,2],y.branin,type="s",
        radius=7,col="red",add=TRUE)
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$lower95),col="blue",add=TRUE)
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$upper95),col="blue",add=TRUE)



# Calculate EI
ego=apply(as.matrix(X),1,EI,fit.gp,type="UK",
          minimization=TRUE)
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(ego,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")


# Get next x that maximizes the EI
x.new=max_EI(fit.gp,lower=rep(0,2),upper=rep(1,2),
             parinit = 0.5,minimization=TRUE)$par
yhat.xnew=predict(fit.gp,newdata=x.new,type="UK")$mean


persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$mean,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")
plot3d(design[,1],design[,2],y.branin,type="s",
        radius=7,col="red",add=TRUE)
plot3d(x.new[1],x.new[2],yhat.xnew,type="s",
        radius=7,col="green",add=TRUE)



# Update by evaluating our expensive function
y.new=apply(x.new,1,branin)
y.branin=c(y.branin,y.new)
design=rbind(design,x.new)

# Refit the GP model
fit.gp=km(~1, design = data.frame(x = design), 
          response = y.branin,covtype = "gauss")
yhat=predict(fit.gp, newdata = data.frame(x = X),
          type = "UK")



# Plot
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$mean,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")
plot3d(design[,1],design[,2],y.branin,type="s",
        radius=7,col="red",add=TRUE)
plot3d(x.new[1],x.new[2],yhat.xnew,type="s",
        radius=7,col="green",add=TRUE)


# Plot with uncertainties
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$mean,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")
plot3d(design[,1],design[,2],y.branin,type="s",
        radius=7,col="red",add=TRUE)
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$lower95),col="blue",add=TRUE)
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$upper95),col="blue",add=TRUE)



# Calculate EI
ego=apply(as.matrix(X),1,EI,fit.gp,type="UK",
          minimization=TRUE)
persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(ego,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")


# Get next x that maximizes the EI
x.new=max_EI(fit.gp,lower=rep(0,2),upper=rep(1,2),
             parinit = 0.5,minimization=TRUE)$par
yhat.xnew=predict(fit.gp,newdata=x.new,type="UK")$mean


persp3d(seq(0,1,length=10),seq(0,1,length=10),
        matrix(yhat$mean,10,10),col="grey",
        xlab="x1",ylab="x2",zlab="y")
plot3d(design[,1],design[,2],y.branin,type="s",
        radius=7,col="red",add=TRUE)
plot3d(x.new[1],x.new[2],yhat.xnew,type="s",
        radius=7,col="green",add=TRUE)


# and so on...






# Sensitivity Analysis Example
# Generate data
set.seed(88) # just to replicate this example
n=500
X=matrix(runif(n*5),ncol=5)
f=10*sin(2*pi*X[,1]*X[,2])+(X[,3]-0.5)^2+X[,4]+X[,5]
y=f+rnorm(n,sd=1)

# true 1-way marginal effects
f1=-10/(2*pi*X[,1])*cos(2*pi*X[,1])+10/(2*pi*X[,1])+13/12
f2=-10/(2*pi*X[,2])*cos(2*pi*X[,2])+10/(2*pi*X[,2])+13/12
f3=3.87964+(X[,3]-0.5)^2+1
f4=3.87964+7/12+X[,4]
f5=3.87964+7/12+X[,5]


# plot
par(mfrow=c(2,3))
plot(X[,1],y,xlab="X1",ylab="Y",pch=20,xlim=c(0,1))
ix=sort(X[,1],index.return=TRUE)$ix
lines(X[ix,1],f1[ix],lwd=4,col="blue")
plot(X[,2],y,xlab="X1",ylab="Y",pch=20,xlim=c(0,1))
ix=sort(X[,2],index.return=TRUE)$ix
lines(X[ix,2],f2[ix],lwd=4,col="blue")
plot(X[,3],y,xlab="X1",ylab="Y",pch=20,xlim=c(0,1))
ix=sort(X[,3],index.return=TRUE)$ix
lines(X[ix,3],f3[ix],lwd=4,col="blue")
plot(X[,4],y,xlab="X1",ylab="Y",pch=20,xlim=c(0,1))
ix=sort(X[,4],index.return=TRUE)$ix
lines(X[ix,4],f4[ix],lwd=4,col="blue")
plot(X[,5],y,xlab="X1",ylab="Y",pch=20,xlim=c(0,1))
ix=sort(X[,5],index.return=TRUE)$ix
lines(X[ix,5],f5[ix],lwd=4,col="blue")


# Calculate the Sobol Indicies
library(sensitivity)
N=10000
X1=data.frame(matrix(runif(N*5),ncol=5))
X2=data.frame(matrix(runif(N*5),ncol=5))
f.test <-function(X) { 
    10*sin(2*pi*X[,1]*X[,2])+(X[,3]-0.5)^2+X[,4]+X[,5] 
}
si.S=sobolEff(model=f.test,X1=X1,X2=X2,order=1,nboot=0)
si.TS=sobolEff(model=f.test,X1=X1,X2=X2,order=0,nboot=0)

# First-order sensitivity indices:
si.S$S

# Total sensitivity indices:
si.TS$S