# psuedo inverse of  a matrix.

svdinv_function(a,tol=1e-6){
 sva_svd(a)
 maxsv_max(sva$d)
 dd_ifelse(sva$d<tol*maxsv,0,1/sva$d)
 t(sva$u %*% (t(sva$v)*dd))
}

# generate multivariate Normals
rmultnorm_function(n,mu,sigma){
   # returns n rows of multinormal mu,sigma random vectors
   # ie: returns a n x length(mu) matrix
  p_length(mu)
  z_matrix(rnorm(n * p),nrow=n)
  ch_chol(sigma,pivot=T)
  piv_attr(ch,"pivot")
  zz_(z%*%ch)
  zzz_0*zz
  zzz[,piv]_zz
  zzz + matrix(mu,nrow=n,ncol=p,byrow=T)
}


# some examples of spatial surfaces for different 
# covariance/precision specifications

# covariance functions
covbm_function(locs,pow=1,scale=1){
  # assume anchored at 1st component
  # browser()
  d1_dist(locs)
  n_dim(locs)[1]
  mat_matrix(0,n,n)
  mat[lower.tri(mat)]_d1
  mat_mat+t(mat)

  vario_(mat/scale)^pow
  covm_matrix(0,n,n)
  vj_matrix(vario[1,],n,n,byrow=T)
  vi_matrix(vario[,1],n,n,byrow=F)
  cc_(vi+vj-vario)/2

  return(cc)
}
covpow_function(locs,pow=2,scale=5){
  # browser()
  d1_dist(locs)
  n_dim(locs)[1]
  mat_matrix(0,n,n)
  mat[lower.tri(mat)]_d1
  mat_mat+t(mat)
  cc_exp(-(mat/scale)^pow)
  return(cc)
}
covsph_function(locs,scale=5){
  # browser()
  d1_dist(locs)
  n_dim(locs)[1]
  mat_matrix(0,n,n)
  mat[lower.tri(mat)]_d1
  mat_mat+t(mat)
  mat_mat/scale
  cc_ifelse(mat>1,0,
     ifelse(mat>0,
     1-2/pi*(sqrt(1-mat^2)*mat + asin(mat)),1))
  return(cc)
}
  
# source("/moses/local/higdon/Sfiles/rmultnorm.s")
# a random walk / brownian motion
locs_expand.grid(1:20,1:20)
locs_cbind(locs[,1],locs[,2])
cbm_covbm(locs)
z_rmultnorm(1,rep(0,20*20),cbm)
z_matrix(z,ncol=20)
persp(z)

# fractal brownian motion
cbm_covbm(locs,pow=1.8)
z_rmultnorm(1,rep(0,20*20),cbm)
z_matrix(z,ncol=20)
persp(z)

# gaussian surface
cga_covpow(locs,pow=2,scale=4)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
persp(z)

# exponentail surface
cga_covpow(locs,pow=1,scale=12)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
persp(z)

# now make some figs for the notes
# exponential
par(mfrow=c(1,2),oma=c(18,0,2,0))
cga_covpow(locs,pow=1,scale=25)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mar=c(0,4,0,0))
persp(z)
par(mar=c(7,4,7,2))
plot(z[,1],ylab="z",xlab="x",ylim=range(z),type="l")
mtext("random realization - exponential covariogram",side=3,line=0,outer=T)

# Gaussian
par(mfrow=c(1,2),oma=c(18,0,2,0))
cga_covpow(locs,pow=2,scale=6)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mar=c(0,4,0,0))
persp(z)
par(mar=c(7,4,7,2))
plot(z[,1],ylab="z",xlab="x",ylim=range(z),type="l")
mtext("random realization - Gaussian covariogram",side=3,line=0,outer=T)

# p = 1.5
par(mfrow=c(1,2),oma=c(18,0,2,0))
cga_covpow(locs,pow=1.5,scale=9)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mar=c(0,4,0,0))
persp(z)
par(mar=c(7,4,7,2))
plot(z[,1],ylab="z",xlab="x",ylim=range(z),type="l")
mtext("random realization - C(d) = exp(-d^1.5)",side=3,line=0,outer=T)

# BM p = 1
par(mfrow=c(1,2),oma=c(18,0,2,0))
cga_covbm(locs,pow=1.0,scale=5)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mar=c(0,4,0,0))
persp(z)
par(mar=c(7,4,7,2))
plot(z[,1],ylab="z",xlab="x",ylim=range(z),type="l")
mtext("random realization - Brownian motion p=1",side=3,line=0,outer=T)

par(mfrow=c(1,2),oma=c(18,0,2,0))
cga_covbm(locs,pow=1.3,scale=1)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mar=c(0,4,0,0))
persp(z)
par(mar=c(7,4,7,2))
plot(z[,1],ylab="z",xlab="x",ylim=range(z),type="l")
mtext("random realization - Brownian motion p=1.3",side=3,line=0,outer=T)

par(mfrow=c(1,2),oma=c(18,0,2,0))
cga_covbm(locs,pow=1.7,scale=1)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mar=c(0,4,0,0))
persp(z)
par(mar=c(7,4,7,2))
plot(z[,1],ylab="z",xlab="x",ylim=range(z),type="l")
mtext("random realization - Brownian motion p=1.7",side=3,line=0,outer=T)

par(mfrow=c(1,2),oma=c(18,0,2,0))
cga_covbm(locs,pow=2.0,scale=1)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mar=c(0,4,0,0))
persp(z)
par(mar=c(7,4,7,2))
plot(z[,1],ylab="z",xlab="x",ylim=range(z),type="l")
mtext("random realization - Brownian motion p=2.0",side=3,line=0,outer=T)

# spherical

par(mfrow=c(1,2),oma=c(18,0,2,0))
cga_covsph(locs,scale=15)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mar=c(0,4,0,0))
persp(z)
par(mar=c(7,4,7,2))
plot(z[,1],ylab="z",xlab="x",ylim=range(z),type="l")
mtext("random realization - spherical",side=3,line=0,outer=T)



# Some interpolation stuff:
# compute the distribution of x2 given x1

source("/moses/local/higdon/Sfiles/svdinv.s")
condgau_function(x,mu,cov,n1){
# computes the conditional distribution of x2 = x[n1+1..n1+n2] given
# the observed values of x1 = x[1..n1].  It is assumed that the
# first n1 values of x hold the observed values.  The function
# returns a vector of length n1+n2, where the first n1 values
# are as were entered in the function.
# browser()
 N_length(x); n2_N-n1
 x1_x[1:n1]
 mu1_mu[1:n1]
 mu2_mu[(n1+1):N]
 s11_cov[1:n1,1:n1]
 s22_cov[(n1+1):N,(n1+1):N]
 s12_cov[1:n1,(n1+1):N]
 s21_t(s12)
 s11i_svdinv(s11)
 #svds11$values_ifelse(svds11$values < 1e-8,0,svds11$values)
 s22o1_s22-s21%*%s11i%*%s12
 mu22o1_mu2+s21%*%s11i%*%(x1-mu1)
 return(mean=mu22o1,cov=(s22o1 + t(s22o1))/2)
}
st_covbm(cbind(1:5,rep(0,5)))
mu_rep(0,5)
x_c(0,1,0,0,0)
n1_2

# Show conditioning on an edge
# Gaussian
cga_covpow(locs,pow=2,scale=6)
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mfrow=c(2,2),mar=c(0,4,0,0))
persp(z)
mtext("a realization",side=3,line=-2)
a_condgau(as.vector(z),rep(0,20*20),cga,20)
cm_c(as.vector(z)[1:20],a$mean)
cm_matrix(cm,ncol=20)
persp(cm,zlim=range(z))
mtext("mean conditional on Y=1 points",side=3,line=-2)
cz_rmultnorm(1,a$mean,a$cov)
cz_c(z[,1],cz)
cz_matrix(cz,ncol=20)
persp(cz,zlim=range(z))
mtext("realization conditional on Y=1 points",side=3,line=-2)

# Brownian motion
locs_expand.grid(1:20,1:20)
locs_cbind(locs[,1],locs[,2])
cga_covbm(locs,pow=1.5,scale=3)+1
z_rmultnorm(1,rep(0,20*20),cga)
z_matrix(z,ncol=20)
par(mfrow=c(2,2),mar=c(0,4,0,0))
persp(z,zlim=range(z)*1.5)
mtext("a realization",side=3,line=-2)
a_condgau(as.vector(z),rep(0,20*20),cga,20)
cm_c(as.vector(z)[1:20],a$mean)
cm_matrix(cm,ncol=20)
persp(cm,zlim=range(z)*1.5)
mtext("mean conditional on Y=1 points",side=3,line=-2)
cz_rmultnorm(1,a$mean,a$cov)
cz_c(z[,1],cz)
cz_matrix(cz,ncol=20)
persp(cz,zlim=range(z)*1.5)
mtext("realization conditional on Y=1 points",side=3,line=-2)

# some simple 1-d examples
zx_c(3,6,9,12,seq(-5,20,length=101))
zy_c(0,0,0,0,rep(0,101))
locs_cbind(zx,zy)
z_c(0,1,1,0,rep(0,101))
#z_c(-1.5,0,0,1.5,rep(0,101))

par(mfcol=c(3,3),mar=c(5,5,3,1))
zl_c(-3.0,3.0)
 # Gaussian
cv_covpow(locs,pow=2,scale=1)
a_condgau(z,rep(0,length(z)),cv,4)
zr_rmultnorm(4,a$mean,a$cov)
matplot(zx[-(1:4)],t(zr),type="l",xlab="x",ylab="z",ylim=zl)
points(zx[1:4],z[1:4],cex=1.6)
lines(zx[-(1:4)],a$mean,lwd=2)
mtext("Gaussian C(r), scale = 2",side=3,line=1)

cv_covpow(locs,pow=2,scale=3)
a_condgau(z,rep(0,length(z)),cv,4)
zr_rmultnorm(4,a$mean,a$cov)
matplot(zx[-(1:4)],t(zr),type="l",xlab="x",ylab="z",ylim=zl)
points(zx[1:4],z[1:4],cex=1.6)
lines(zx[-(1:4)],a$mean,lwd=2)
mtext("Gaussian C(r), scale = 3",side=3,line=1)

cv_covpow(locs,pow=2,scale=5)
a_condgau(z,rep(0,length(z)),cv,4)
zr_rmultnorm(4,a$mean,a$cov)
matplot(zx[-(1:4)],t(zr),type="l",xlab="x",ylab="z",ylim=zl)
points(zx[1:4],z[1:4],cex=1.6)
lines(zx[-(1:4)],a$mean,lwd=2)
mtext("Gaussian C(r), scale = 5",side=3,line=1)

 # Exponential
cv_covpow(locs,pow=1,scale=1)
a_condgau(z,rep(0,length(z)),cv,4)
zr_rmultnorm(4,a$mean,a$cov)
matplot(zx[-(1:4)],t(zr),type="l",xlab="x",ylab="z",ylim=zl)
points(zx[1:4],z[1:4],cex=1.6)
lines(zx[-(1:4)],a$mean,lwd=2)
mtext("Exponential C(r), scale = 1",side=3,line=1)

cv_covpow(locs,pow=1,scale=10)
a_condgau(z,rep(0,length(z)),cv,4)
zr_rmultnorm(4,a$mean,a$cov)
matplot(zx[-(1:4)],t(zr),type="l",xlab="x",ylab="z",ylim=zl)
points(zx[1:4],z[1:4],cex=1.6)
lines(zx[-(1:4)],a$mean,lwd=2)
mtext("Exponential C(r), scale = 10",side=3,line=1)

cv_covpow(locs,pow=1,scale=20)
a_condgau(z,rep(0,length(z)),cv,4)
zr_rmultnorm(4,a$mean,a$cov)
matplot(zx[-(1:4)],t(zr),type="l",xlab="x",ylab="z",ylim=zl)
points(zx[1:4],z[1:4],cex=1.6)
lines(zx[-(1:4)],a$mean,lwd=2)
mtext("Exponential C(r), scale = 20",side=3,line=1)

 # Brownian motion p=1.5
cv_covbm(locs,pow=1.5,scale=1)+1
a_condgau(z,rep(0,length(z)),cv,4)
zr_rmultnorm(4,a$mean,a$cov)
matplot(zx[-(1:4)],t(zr),type="l",xlab="x",ylab="z",ylim=zl)
points(zx[1:4],z[1:4],cex=1.6)
lines(zx[-(1:4)],a$mean,lwd=2)
mtext("Brownian motion C(r), p = 1.5 scale = 1",side=3,line=1)

cv_covbm(locs,pow=1.5,scale=3)+1
a_condgau(z,rep(0,length(z)),cv,4)
zr_rmultnorm(4,a$mean,a$cov)
matplot(zx[-(1:4)],t(zr),type="l",xlab="x",ylab="z",ylim=zl)
points(zx[1:4],z[1:4],cex=1.6)
lines(zx[-(1:4)],a$mean,lwd=2)
mtext("Brownian motion C(r), p = 1.5 scale = 3",side=3,line=1)

cv_covbm(locs,pow=1.5,scale=5)+1
a_condgau(z,rep(0,length(z)),cv,4)
zr_rmultnorm(4,a$mean,a$cov)
matplot(zx[-(1:4)],t(zr),type="l",xlab="x",ylab="z",ylim=zl)
points(zx[1:4],z[1:4],cex=1.6)
lines(zx[-(1:4)],a$mean,lwd=2)
mtext("Brownian motion C(r), p = 1.5 scale = 5",side=3,line=1)