R version 4.0.0 (2020-04-24) -- "Arbor Day"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>  # 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")
Testing:  Krig mKrig d coef
PASSED test at tolerance  1e-08
> test.for.zero( look$c, test$c, tag="Krig mKrig c coef")
Testing:  Krig mKrig c coef
PASSED test at tolerance  1e-08
> 
> # 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")
Testing:  Monte Carlo eff.df
PASSED test at tolerance  0.01
> 
> 
> # 
> 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")
Testing:  Test of wendland and great circle
PASSED test at tolerance  1e-08
> 
> test.for.zero(look$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D)
+               ,tol=.01, tag="eff.df")
Testing:  eff.df
PASSED test at tolerance  0.01
> 
> # 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")
Testing:  Test of sparse wendland and great circle
PASSED test at tolerance  8e-07
> test.for.zero(look4$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D),
+                         tol=.01, tag="sparse eff.df")
Testing:  sparse eff.df
PASSED test at tolerance  0.01
> 
> # 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")
Testing:  exact sparse eff.df
PASSED test at tolerance  0.01
> 
> # 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)
PASSED test at tolerance  5e-07
> test.for.zero( look3$d, lookKrig$d, tol=2e-8)
PASSED test at tolerance  2e-08
> test.for.zero( look3$fitted.values, lookKrig$fitted.values, tol=1e-7)
PASSED test at tolerance  1e-07
> 
> test.for.zero( predict( look3, xnew= look3$x), predict( lookKrig, xnew= lookKrig$x),
+                tol=5e-7)
PASSED test at tolerance  5e-07
> 
> test.for.zero( hold[,1], hold2[,1], tol=1e-7, relative=FALSE)
PASSED test at tolerance  1e-07
> 
> test.for.zero(diag(hold),diag(hold2), tol=2E-7,
+               relative=FALSE, tag="exact sparse eff.df by predict -- fastTps")
Testing:  exact sparse eff.df by predict -- fastTps
PASSED test at tolerance  2e-07
> #plot( diag(hold), ( 1- diag(hold2)/ diag(hold))  )
> 
> test.for.zero(look3$eff.df,sum( diag(hold)) , tag="fastTps ef.df exact" )
Testing:  fastTps ef.df exact
PASSED test at tolerance  1e-08
> 
> test.for.zero(look3$eff.df, Krig.ftrace( lookKrig$lambda, lookKrig$matrices$D),
+                    tol=2e-7, tag="exact sparse eff.df through mKrig-- fastTps")
Testing:  exact sparse eff.df through mKrig-- fastTps
PASSED test at tolerance  2e-07
> 
> # 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")
Testing:  initial check on d for likelihood
PASSED test at tolerance  1e-08
> 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")
Testing:  rho hat from likelihood
PASSED test at tolerance  1e-08
> 
> test.for.zero( look, out2$rhohat, tag="rho hat from likelihood compared to Krig")
Testing:  rho hat from likelihood compared to Krig
PASSED test at tolerance  1e-08
> 
> 
> 
> # 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)
PASSED test at tolerance  1e-08
> test.for.zero( out$lnDetCov, determinant(M, log=TRUE)$modulus)
PASSED test at tolerance  1e-08
> 
> # 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)
PASSED test at tolerance  1e-08
> test.for.zero(  look2, determinant(M, log=TRUE)$modulus)
PASSED test at tolerance  1e-08
> test.for.zero( out$lnDetCov, determinant(M, log=TRUE)$modulus)
PASSED test at tolerance  1e-08
> 
> 
> 
> # 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)
> 
> 
> 
>             
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> 
> proc.time()
   user  system elapsed 
  1.633   0.362   1.969