# fields is a package for analysis of spatial data written for
  # the R software environment .
  # Copyright (C) 2018
  # University Corporation for Atmospheric Research (UCAR)
  # Contact: Douglas Nychka, nychka@ucar.edu,
  # National Center for Atmospheric Research,
  # PO Box 3000, Boulder, CO 80307-3000
  #
  # This program is free software; you can redistribute it and/or modify
  # it under the terms of the GNU General Public License as published by
  # the Free Software Foundation; either version 2 of the License, or
  # (at your option) any later version.
  # This program is distributed in the hope that it will be useful,
  # but WITHOUT ANY WARRANTY; without even the implied warranty of
  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  # GNU General Public License for more details.


suppressMessages(library(fields))
#options( echo=FALSE)
test.for.zero.flag<- 1
data(ozone2)
y<- ozone2$y[16,]
x<- ozone2$lon.lat
#
# Omit the NAs
good<- !is.na( y)
x<- x[good,]
y<- y[good]
#source("~/Home/Src/fields/R/mKrig.family.R")

# now look at mKrig w/o sparse matrix 
look<- mKrig( x,y, cov.function="stationary.cov", theta=10, lambda=.3,
                      chol.args=list( pivot=FALSE))


lookKrig<- Krig( x,y, cov.function="stationary.cov",
               theta=10) 

test.df<-Krig.ftrace(look$lambda,lookKrig$matrices$D)

test<- Krig.coef( lookKrig, lambda=look$lambda)

test.for.zero( look$d, test$d, tag="Krig mKrig d coef")
test.for.zero( look$c, test$c, tag="Krig mKrig c coef")

# test of trace calculation

look<- mKrig( x,y, cov.function="stationary.cov", theta=10, lambda=.3,
         
          find.trA=TRUE, NtrA= 1000, iseed=243)

test.for.zero( look$eff.df, test.df,tol=.01, tag="Monte Carlo eff.df")


# 
lookKrig<-Krig( x,y, cov.function="stationary.cov",
               theta=350, Distance="rdist.earth",Covariance="Wendland", 
               cov.args=list( k=2, dimension=2) ) 

look<- mKrig( x,y, cov.function="stationary.cov", 
        theta=350, 
        Distance="rdist.earth",Covariance="Wendland",  
        cov.args=list( k=2, dimension=2),
        lambda=lookKrig$lambda,
        find.trA=TRUE, NtrA= 1000, iseed=243)

test.for.zero( look$c, lookKrig$c, tag="Test of wendland and great circle")

test.for.zero(look$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D)
              ,tol=.01, tag="eff.df")

# same calculation using sparse matrices.

look4<- mKrig( x,y, cov.function="wendland.cov", 
        theta=350, 
        Dist.args=list( method="greatcircle"),  
        cov.args=list( k=2),
        lambda=lookKrig$lambda,
        find.trA=TRUE, NtrA=500, iseed=243)

test.for.zero( look$c, look4$c,tol=8e-7, 
           tag="Test of sparse wendland and great circle")
test.for.zero(look4$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D),
                        tol=.01, tag="sparse eff.df")

# great circle distance switch has been a  big bug -- test some options

look<- mKrig( x,y, cov.function="wendland.cov", 
 theta=350, Dist.args=list( method="greatcircle"),  
 cov.args=list( k=2),lambda=lookKrig$lambda,
 find.trA=TRUE, NtrA=1000, iseed=243)

test.for.zero(look$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D),
                   tol=1e-2, tag="exact sparse eff.df")

# compare to fast Tps 
look3<-  fastTps( x,y,theta=350,lambda=lookKrig$lambda, NtrA=200, iseed=243, 
                lon.lat=TRUE)
#look3$c<- lookKrig$c
#look3$d<-  lookKrig$d
object<- look3
np<- object$np
Ey <- diag(1, np)
NtrA <- np
hold <- predict.mKrig(object, ynew = Ey, collapseFixedEffect=FALSE)
hold2<- matrix( NA, np,np)
for(  k in 1:np){
hold2[,k] <- predict.Krig(lookKrig, y = Ey[,k])
}
#plot( diag(hold), diag(hold2))


test.for.zero( look3$c, lookKrig$c, tol=5e-7)
test.for.zero( look3$d, lookKrig$d, tol=2e-8)
test.for.zero( look3$fitted.values, lookKrig$fitted.values, tol=1e-7)

test.for.zero( predict( look3, xnew= look3$x), predict( lookKrig, xnew= lookKrig$x),
               tol=5e-7)

test.for.zero( hold[,1], hold2[,1], tol=1e-7, relative=FALSE)

test.for.zero(diag(hold),diag(hold2), tol=2E-7,
              relative=FALSE, tag="exact sparse eff.df by predict -- fastTps")
#plot( diag(hold), ( 1- diag(hold2)/ diag(hold))  )

test.for.zero(look3$eff.df,sum( diag(hold)) , tag="fastTps ef.df exact" )

test.for.zero(look3$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D),
                   tol=2e-7, tag="exact sparse eff.df through mKrig-- fastTps")

# calculations of likelihood, rho and sigma

lam<-.2

out<- mKrig( x,y, cov.function =Exp.cov, theta=4, lambda=lam)
out2<- Krig( x,y, cov.function =Exp.cov, theta=4, lambda=lam)

            
Sigma<- Exp.cov( x,x,theta=4)
X<-  cbind( rep(1, nrow(x)), x)

Sinv<- solve( Sigma + lam* diag( 1, nrow( x)))

#checks on  likelihoods            

# quadratic form:
dhat<- c(solve( t(X)%*%Sinv%*%(X) ) %*% t(X) %*%Sinv%*%y)
test.for.zero( dhat, out$d, tag="initial check on d for likelihood")
r<- y -X%*%dhat
N<- nrow(x)
look<-  t( r)%*%(Sinv)%*%r/N



test.for.zero( look, out$rho.MLE, tag="rho hat from likelihood")

test.for.zero( look, out2$rhohat, tag="rho hat from likelihood compared to Krig")



# check determinant
lam<- .2
Sigma<- Exp.cov( x,x,theta=4)
M<- Sigma + lam * diag( 1, nrow(x))
chol( M)-> Mc
look2<- sum( log(diag( Mc)))*2

out<-mKrig( x,y,cov.function =Exp.cov, theta=4, lambda=lam)

test.for.zero( out$lnDetCov, look2)
test.for.zero( out$lnDetCov, determinant(M, log=TRUE)$modulus)

# weighted version 
lam<- .2
Sigma<- Exp.cov( x,x,theta=4)
set.seed( 123)
weights<- runif(nrow( x))
M<- Sigma +  diag(lam/ weights)
chol( M)-> Mc
look2<- sum( log(diag( Mc)))*2

out<-mKrig( x,y,weights=weights, cov.function =Exp.cov, theta=4, lambda=lam)

test.for.zero( out$lnDetCov, look2)
test.for.zero(  look2, determinant(M, log=TRUE)$modulus)
test.for.zero( out$lnDetCov, determinant(M, log=TRUE)$modulus)



# check profile likelihood by estimating MLE
lam.true<- .2
N<- nrow( x)
Sigma<- Exp.cov( x,x,theta=4)
M<- Sigma + lam.true * diag( 1, nrow(x))
chol( M)-> Mc
t(Mc)%*%Mc -> test




##D set.seed( 234)
##D NSIM<- 100
##D hold2<-rep( NA, NSIM)
##D temp.fun<- function(lglam){
##D             out<-mKrig( x,ytemp,
##D                         cov.function =Exp.cov, theta=4, lambda=exp(lglam))
##D             return(-1* out$lnProfileLike)}

##D hold1<-rep( NA, NSIM)
##D yt<- rep( 1, N) 
##D obj<- Krig( x,yt, theta=4)


##D E<- matrix( rnorm( NSIM*N), ncol=NSIM)

##D for ( j in 1:NSIM){
##D cat( j, " ")
##D ytemp <- x%*%c(1,2) +  t(Mc)%*%E[,j] 
##D out<- optim( log(.2), temp.fun, method="BFGS")
##D hold2[j]<- exp(out$par)
##D hold1[j]<-  gcv.Krig(obj, y=ytemp)$lambda.est[6,1]

##D }
##D test.for.zero( median( hold1), .2, tol=.08)
##D test.for.zero( median( hold2), .2, tol=.12)