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

source("dace.sim.r")




# Example in 2D with Gaussian Correlation
library(rgl); rgl.clear()
set.seed(66)

n=25
x1=runif(n)
x2=runif(n)
X=as.matrix(cbind(x1,x2))
l1=list(m1=abs(outer(X[,1],X[,1],"-")))
l2=list(m2=abs(outer(X[,2],X[,2],"-")))
l.dez=list(l1=l1,l2=l2)

rho=c(0.9,0.1)
R=rhogeodacecormat(l.dez,rho)$R

L=t(chol(R+diag(n)*.Machine$double.eps*100))
Z=rnorm(n)
Y=L%*%Z

rho1=seq(.001,1,length=20)
rho2=seq(.001,1,length=20)
rho.grid=as.matrix(expand.grid(rho1,rho2))
l.vals=rep(0,nrow(rho.grid))

for(i in 1:nrow(rho.grid))
    l.vals[i]=logl(rho.grid[i,],Y,l.dez,
                   conditioning=.Machine$double.eps*1000)

persp3d(matrix(l.vals,20,20),col="grey",xlab="rho1",
        ylab="rho2",zlab="logl",xlim=range(rho1),
        ylim=range(rho2))





# Example in 2D with Gaussian Correlation
rgl.clear()
set.seed(66)

n=15
x1=runif(n)
x2=runif(n)
X=as.matrix(cbind(x1,x2))
l1=list(m1=abs(outer(X[,1],X[,1],"-")))
l2=list(m2=abs(outer(X[,2],X[,2],"-")))
l.dez=list(l1=l1,l2=l2)

rho=c(0.9,0.1)
R=rhogeodacecormat(l.dez,rho)$R

L=t(chol(R+diag(n)*.Machine$double.eps*100))
Z=rnorm(n)
Y=L%*%Z

rho1=seq(.001,1,length=20)
rho2=seq(.001,1,length=20)
rho.grid=as.matrix(expand.grid(rho1,rho2))
l.vals=rep(0,nrow(rho.grid))

for(i in 1:nrow(rho.grid))
    l.vals[i]=logl(rho.grid[i,],Y,l.dez,conditioning=.Machine$double.eps*1000)

persp3d(matrix(l.vals,20,20),col="grey",xlab="rho1",ylab="rho2",zlab="logl",xlim=range(rho1),ylim=range(rho2))




# Example in 2D with Gaussian Correlation
rgl.clear()
set.seed(66)

n=8
x1=runif(n)
x2=runif(n)
X=as.matrix(cbind(x1,x2))
l1=list(m1=abs(outer(X[,1],X[,1],"-")))
l2=list(m2=abs(outer(X[,2],X[,2],"-")))
l.dez=list(l1=l1,l2=l2)

rho=c(0.9,0.1)
R=rhogeodacecormat(l.dez,rho)$R

L=t(chol(R+diag(n)*.Machine$double.eps*100))
Z=rnorm(n)
Y=L%*%Z

rho1=seq(.001,1,length=20)
rho2=seq(.001,1,length=20)
rho.grid=as.matrix(expand.grid(rho1,rho2))
l.vals=rep(0,nrow(rho.grid))

for(i in 1:nrow(rho.grid))
    l.vals[i]=logl(rho.grid[i,],Y,l.dez,conditioning=.Machine$double.eps*1000)

persp3d(matrix(l.vals,20,20),col="grey",xlab="rho1",ylab="rho2",zlab="logl",xlim=range(rho1),ylim=range(rho2))





# Same example with penalized likelihood
l.vals=rep(0,nrow(rho.grid))
lambda=0.1

for(i in 1:nrow(rho.grid))
    l.vals[i]=logl(rho.grid[i,],Y,l.dez,conditioning=.Machine$double.eps*1000)+n*lambda*(rho.grid[i,1]^2*(1-rho.grid[i,1])^2+rho.grid[i,2]^2*(1-rho.grid[i,2])^2)

rgl.clear()
persp3d(matrix(l.vals,20,20),col="grey",xlab="rho1",ylab="rho2",zlab="logl",xlim=range(rho1),ylim=range(rho2))





# Same example with penalized likelihood
rgl.clear()
l.vals=rep(0,nrow(rho.grid))
lambda=5

for(i in 1:nrow(rho.grid))
    l.vals[i]=logl(rho.grid[i,],Y,l.dez,conditioning=.Machine$double.eps*1000)+n*lambda*(rho.grid[i,1]^2*(1-rho.grid[i,1])^2+rho.grid[i,2]^2*(1-rho.grid[i,2])^2)

persp3d(matrix(l.vals,20,20),col="grey",xlab="rho1",ylab="rho2",zlab="logl",xlim=range(rho1),ylim=range(rho2))




# Same example with penalized likelihood
rgl.clear()
l.vals=rep(0,nrow(rho.grid))
lambda=10

for(i in 1:nrow(rho.grid))
    l.vals[i]=logl(rho.grid[i,],Y,l.dez,conditioning=.Machine$double.eps*1000)+n*lambda*(rho.grid[i,1]^2*(1-rho.grid[i,1])^2+rho.grid[i,2]^2*(1-rho.grid[i,2])^2)

persp3d(matrix(l.vals,20,20),col="grey",xlab="rho1",ylab="rho2",zlab="logl",xlim=range(rho1),ylim=range(rho2))






# Prediction and Uncertainty Quantification
# 1D Example - generate our data
set.seed(77)

n=4
x1=seq(0.1,0.9,length=n)+runif(n,-.05,.05)
l1=list(m1=abs(outer(x1,x1,"-")))
l.dez=list(l1=l1)

rho=c(0.001)
R=rhogeodacecormat(l.dez,rho)$R

L=t(chol(R+diag(n)*.Machine$double.eps*100))
Z=rnorm(n)
Y=L%*%Z


# 1D Example - estimate via MLE
rho.hat=optimize(logl,interval=c(0,1),Y,l.dez,lower=0,
        upper=0.9,maximum=TRUE)$maximum
rho.hat

R=rhogeodacecormat(l.dez,rho.hat,2)$R
cR=chol(R+diag(n)*.Machine$double.eps*0)
Rinv=chol2inv(cR)
s2.hat=(1/n)*t(Y)%*%Rinv%*%Y
s2.hat


# 1D Example - predict and pointwise uncertainty intervals
ngrid=100
pred.grid=seq(0,1,length=100)
X=c(pred.grid,x1)
l1=list(m1=abs(outer(X,X,"-")))
l.dez=list(l1=l1)

Rall=rhogeodacecormat(l.dez,rho.hat)$R
R0=Rall[1:ngrid,(ngrid+1):(ngrid+n)]
rm(Rall)

yhat=R0%*%Rinv%*%Y
s2hat=s2.hat*diag(1-R0%*%Rinv%*%t(R0))


# 1D Example - plot
plot(x1,Y,pch=20,cex=2,xlim=c(0,1),ylim=c(-2,2),
     xlab="x",ylab="y(x)")
lines(pred.grid,yhat,lwd=2,col="black")
lines(pred.grid,yhat-1.96*sqrt(s2hat),col="grey")
lines(pred.grid,yhat+1.96*sqrt(s2hat),col="grey")






# 2D Example - generate our data
set.seed(99)

n=10
x1=runif(n)
x2=runif(n)
X=cbind(x1,x2)
l1=list(m1=abs(outer(X[,1],X[,1],"-")))
l2=list(m2=abs(outer(X[,2],X[,2],"-")))
l.dez=list(l1=l1,l2=l2)

rho=c(0.001,0.5)
R=rhogeodacecormat(l.dez,rho)$R

L=t(chol(R+diag(n)*.Machine$double.eps*100))
Z=rnorm(n)
Y=L%*%Z


# 2D Example - estimate via MLE
rho.hat=optim(c(0.5,0.5),logl,gr=NULL,lower=0,
        upper=0.9,method="L-BFGS-B",
        control=list(fnscale=-1),Y,l.dez)$par
rho.hat

R=rhogeodacecormat(l.dez,rho.hat,2)$R
cR=chol(R+diag(n)*.Machine$double.eps*0)
Rinv=chol2inv(cR)
s2.hat=(1/n)*t(Y)%*%Rinv%*%Y
s2.hat


# 2D Example - predict and pointwise uncertainty intervals
ngrid=10
pred.grid=as.matrix(expand.grid(seq(0,1,length=ngrid),
                    seq(0,1,length=ngrid)))
Xall=rbind(pred.grid,X)
l1=list(m1=abs(outer(Xall[,1],Xall[,1],"-")))
l2=list(m2=abs(outer(Xall[,2],Xall[,2],"-")))
l.dez=list(l1=l1,l2=l2)

Rall=rhogeodacecormat(l.dez,rho.hat)$R
R0=Rall[1:(ngrid^2),(ngrid^2+1):(ngrid^2+n)]
rm(Rall)

yhat=R0%*%Rinv%*%Y
se.pred=sqrt(s2.hat*diag(1-R0%*%Rinv%*%t(R0)))


# 2D Example - plot
rgl.clear()
plot3d(X[,1],X[,2],Y,type='s',radius=0.075,col="red",
       xlim=c(0,1),ylim=c(0,1),zlim=c(-3,3))
persp3d(seq(0,1,length=ngrid),seq(0,1,length=ngrid),
        yhat,col="blue",add=TRUE)
persp3d(seq(0,1,length=ngrid),seq(0,1,length=ngrid),
       yhat-1.96*se.pred,col="grey",add=TRUE)
persp3d(seq(0,1,length=ngrid),seq(0,1,length=ngrid),
       yhat+1.96*se.pred,col="grey",add=TRUE)






# Dealing with Computational Limitations
N=seq(10,1010,by=50)
times=rep(0,length(N))
for(i in 1:length(N))
{
    n=N[i]; x1=runif(n); x2=runif(n); X=cbind(x1,x2)
    l1=list(m1=abs(outer(X[,1],X[,1],"-")))
    l2=list(m2=abs(outer(X[,2],X[,2],"-")))
    l.dez=list(l1=l1,l2=l2)

    rho=c(0.01,0.01)
    elapt=system.time({
    R=rhogeodacecormat(l.dez,rho)$R;
    Ri=chol2inv(chol(R+diag(n)*.Machine$double.eps*100));
    rm(R); rm(Ri)
    })
    times[i]=elapt[[1]]
}

par(mfrow=c(1,2))
plot(N,times,pch=20,xlab="sample size",ylab="time (s)")
plot(N,times^(1/3),pch=20,xlab="sample size",ylab="time^(1/3)")




# Compact Correlation Functions - Wendland and friends
x=seq(0,1,length=100)
l1=list(m1=abs(outer(x,x,"-")))
l.dez=list(l1=l1)
h=sqrt((0-x)^2)

R.gauss=rhogeodacecormat(l.dez,rho=0.001,alpha=2)$R
R.exp=rhogeodacecormat(l.dez,rho=0.001,alpha=1)$R
R.wend1=wendland1(l.dez,0.5)$R
R.wend2=wendland2(l.dez,0.5)$R
R.gw1=generalized.wendland(l.dez,0.5,1)$R
R.gw2=generalized.wendland(l.dez,0.5,2)$R
R.gw3=generalized.wendland(l.dez,0.5,3)$R
R.gw4=generalized.wendland(l.dez,0.5,4)$R

plot(h,R.gauss[1,],type='l',ylim=c(0,1),col="grey",lwd=3,ylab="R(h)")
lines(h,R.exp[1,],col="grey",lwd=3)
lines(h,R.wend1[1,],lwd=1,lty=2)
lines(h,R.wend2[1,],lwd=1,lty=2)
lines(h,R.gw1[1,],lwd=1,lty=3)
lines(h,R.gw2[1,],lwd=1,lty=3)
lines(h,R.gw3[1,],lwd=1,lty=3)
lines(h,R.gw4[1,],lwd=1,lty=3)