library(MASS)      # for stepAIC
library(sm)        # for sm.density in diagnostics
library(akima)     # for interp to smooth for plotting
library(fields)    # for surface to plot surface plot
library(FNN)       # for k nearest neighbor to snap grids together
library(locfit)    # for use of local polynomial regression
library(MPV)
library(prodlim)


######## READING IN DATA
data=read.table(file='http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/climatol-ann.txt')

Y=data$V4      #ann avg precip at 246 locations
Y = Y/10		#convert to cm
X=data[,1:3]   #Parameters (log,lat,elev)
names(X)=c("lon","lat","elev")

N = length(Y)  #length of the data
nvar=3
porder=1
minalpha=2*(nvar*porder+1)/N
alpha1=seq(minalpha,1.0,by=0.05)
n=length(alpha1)

porder=2
minalpha=2*(nvar*porder+1)/N
alpha2=seq(minalpha,1.0,by=0.05)
alpha=c(alpha1,alpha2)

#get the GCV values for all the alpha values in alpha for order of
# polynomial = 1. kern="bisq" argument is to use the bisquare kernel
# in computing the weights of the neighbors, which are then used in
# the weighted least squares..

zz=gcvplot(Y ~lon+lat+elev,data=X, alpha=alpha1,deg=1,kern="bisq",ev=dat(),scale=TRUE,maxk=5000)
z1=gcvplot(Y ~lon+lat+elev,data = X, alpha=alpha2,deg=2,kern="bisq",ev=dat(),scale=TRUE,maxk=5000)

# pick the best alpha and the degree of the polynomial that
# gives the least GCV
z2=order(c(zz$values,z1$values))
deg1=1
if(z2[1] > n)deg1=2
bestalpha=alpha[z2[1]]
bestdeg=deg1

##### just order 1
z2=order(zz$values)
bestdeg=1
bestalpha=alpha[z2[1]]

zz=locfit(Y~lon+lat+elev,data=X,alpha=bestalpha,deg=bestdeg,kern="bisq",scale=TRUE,ev=dat(),maxk=5000)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~ii~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
##Test Goodness of Fit  (F-Test)
Yp = predict(zz)
RSS1 = sum((Y-Yp)^2)
#nu1 = sum(fitted(zz,what="infl"))     # trace(L)
#nu2 = sum(fitted(zz,what="vari"))     # trace(L^T L)

nu1 = zz$dp[6]
nu2 = zz$dp[7]

nu11 = N-2*nu1 + nu2

#linear regression ###
zzl=locfit(Y~V1+V2+V3,data=X,alpha=1,deg=1,kern="bisq",scale=TRUE,ev=none(),maxk=5000)

Yp = predict(zzl)
RSS0 = sum((Y-Yp)^2)

hatm = fitted(zzl,what="infl")
II = diag(N)
delta0 = t(II-hatm)%*%(II-hatm)    #Equation 9.2
nu00 = sum(diag(delta0))

Fdata = (RSS0 - RSS1)/(nu00 - nu11)
Fdata = (Fdata / (RSS1 / nu11))
Ftheor = qf(0.95,(nu00-nu11), nu11)    #95% confidence level..

## Fdata > Ftheor   - reject null - i.e., data is otherwise (local polynomial)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~v~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#
## read te India DEM and estimate on it
### Create Rajeevan grid with elevation
## read the Rajeevan grid (lat and lon)
par(mfrow=c(1,1))

### India and Central Asia topo
test1=matrix(scan("india-grid-topo.txt"),ncol=3,byrow=T)

### Rajeevan grid
test2=matrix(scan("http://civil.colorado.edu/~balajir/CVEN6833/HWs/HW-1/Rajeevan-grid.txt"),ncol=2,byrow=T)
### put longitude first
test3=cbind(test2[,2],test2[,1])

xv1=data.frame(test1[,1:2])
xv2=data.frame(test3[,1:2])

## find common points in the top grid that contains
## Rajeevan grid points and the missing points
##

zzint=row.match(xv2,xv1)
indexnotna = which(!is.na(zzint))
indexna = which( is.na(zzint) )

## create the new rajeevan grid
rajeevgridel = cbind(test3[indexnotna,1:2],test1[zzint[indexnotna],3])  

X=data[,1:3]
names(X)=c("lon","lat","elev")
zz=locfit(Y~lon+lat+elev,data=X,alpha=bestalpha, deg=bestdeg, kern="bisq", scale=TRUE,maxk=5000)

X = as.data.frame(rajeevgridel)
names(X)=c("lon","lat","elev")
ypred <- predict(zz,newdata=X,se.fit=T)

############ Simple plotting..
quilt.plot(X[,1],X[,2],ypred$fit)

##### Splines
X=data[,1:3]
names(X)=c("lon","lat","elev")
ztps = Tps(X,Y)

X = as.data.frame(rajeevgridel)
names(X)=c("lon","lat","elev")
ypred <- predict(ztps,X)

sefit =  predictSE.Krig(ztps,X)

############# Simple plotting ########

quilt.plot(X[,1],X[,2],ypred)


########### plotting

nsta = 357   # number of grid points over India
ygrid=seq(6.5,38.5,by=1)
ny=length(ygrid)
ncols=ny

xgrid=seq(66.5,100.5,by=1)
nx = length(xgrid)
nrows=nx

## Rectangular area cover India
index=1:(nx*ny)

## grid points of India
index1 = scan("india-grid-index.txt")  

zfull = rep(NaN,length(index))

## replace the negative model values to zero
ypredpos = ypred$fit
ypredpos[ypredpos < 0]=0

zfull[index1]=ypredpos
zfull[index1[indexnotna]]=ypredpos

zmat = matrix(zfull,nrow=nrows,ncol=ncols)

image.plot(xgrid,ygrid,zmat,ylim=range(6,40),xlim=range(60,100),xlab="Long",ylab="Lat")
world(ylim=range(6,35),xlim=range(60,95),add=TRUE)

### plot elevation
zfull = rep(NaN,length(index))
zfull[index1]=rajeevgridel[,3]

zmat = matrix(zfull,nrow=nrows,ncol=ncols)

image.plot(xgrid,ygrid,zmat,ylim=range(6,40),xlim=range(60,100),xlab="Long",ylab="Lat")
world(ylim=range(6,35),xlim=range(60,95),add=TRUE)
##########################################################