R Under development (unstable) (2024-06-26 r86840 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

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) 2022 Colorado School of Mines
> # 1500 Illinois St., Golden, CO 80401
> # Contact: Douglas Nychka,  douglasnychka@gmail.edu,
> #
> # 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.
> #
> # You should have received a copy of the GNU General Public License
> # along with the R software environment if not, write to the Free Software
> # Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
> # or see http://www.r-project.org/Licenses/GPL-2
> ##END HEADER
> ##END HEADER
> 
> 
> # test of sreg and related functions
> 
> suppressMessages(library(fields))
> #options(echo=FALSE)
> 
> test.for.zero.flag<- 1
> 
> 
> 
> # simple covariance function for implementation
> exp_cov <- function(dist){
+   sigma2<-1 
+   covariance <- sigma2* exp(-dist / 10) # 10 is arbitrary 
+   return(covariance)
+ }
> 
> # -----------------------------
> # Define grid and observations
> # -----------------------------
>  
> m<- 10
> n<- 11
> nx<- m
> ny<- n
> M<- 15
> dx<- 1
> dy<- 1
> sigma2<-2.0 
> np<-3 
> 
> 
> # first a case where all obs in same grid box.
> # addition of "dx" also tests that this works when grid is not just integers
> # set dx=1 for the most basic case
> dx<- .5 
> s0<- rbind( 
+             c(5.1,6.2),
+             c(5.1,6.5),
+             c( 5.85,6.45)
+             )
>  s0<- s0*dx
>  
> test0<-  offGridWeights( s0, list( x= (1:m)*dx, y=(1:n)*dx),
+                                               aRange=10*dx, sigma2=sigma2, 
+                                               Covariance="Exponential", 
+                                               np=2, 
+                                               debug=TRUE) 
[1] 2
Found 1 grid box(es) containing more than 1 obs location
> # explicit nearest neighbors in this case
> sTmp<- cbind( rep(4:7,4), rep(5:8,each=4) )
> sGrid<- sTmp*dx
> 
> # check that same grid being used by function
> test.for.zero(sGrid, cbind(test0$gridX[,1], test0$gridY[,1]) )
PASSED test at tolerance  1e-08
> 
> S21  <- 2.0* exp( -rdist( s0, sGrid)/(10*dx))
> S11  <- 2.0* exp( -rdist( sGrid , sGrid)/(10*dx) )
> S22 <-  2.0* exp( -rdist( s0, s0)/(10*dx))
> # local weights applied for prediction 
> Btest<- S21%*% solve( S11)
> # find indices for neigborhood
> sIndex<- sTmp[,1] + (sTmp[,2]-1)*m
> # Kriging weights
> Bfull<-  spam2full(test0$B[,sIndex])
> test.for.zero( Bfull, Btest)
PASSED test at tolerance  1e-08
> # standard error matrix 
> # note that transpsoe also taken so SEtest%*%t( SEtest) = cov matrix
> SEtest<-  t(chol(S22 - S21%*% solve( S11)%*%t(S21) ))
> SEfull <- spam2full(test0$SE)
> test.for.zero( SEfull, SEtest)
PASSED test at tolerance  1e-08
> 
> 
> # now test several observation locations
> 
>  dx<- .45
>  s<- rbind( 
+    c(5.1,6.2),
+    c(7.1,7.2),
+    c(5.1,6.5),
+    c(8.5,4.4),
+    c( 5.85,6.45),
+    c(7.3,7.4)
+  )
>  s<- s * dx
> # Note s0 from above is  s[c(1,3,5),]
> ind1<- c(1,3,5)
>  
>  sTmp<- cbind( rep(4:7,4), rep(5:8,each=4) )
>  sGrid<- sTmp*dx
>  sIndex<- sTmp[,1] + (sTmp[,2]-1)*m
>  
>  S21<- 2.0* exp( -rdist( s[ind1,], sGrid)/(10*dx) )
>  S11<- 2.0* exp( -rdist( sGrid , sGrid)/(10*dx) )
>  S22<- 2.0* exp( -rdist( s[ind1,], s[ind1,])/(10*dx) )
>  
>  sparseObj<-  offGridWeights( s, list( x= (1:m)*dx, y=(1:n)*dx),
+                                  aRange=(10*dx), sigma2=sigma2, 
+                                  Covariance="Exponential", 
+                                  np=2, 
+                                  debug=TRUE)
[1] 2
Found 2 grid box(es) containing more than 1 obs location
>  
> test.for.zero( sparseObj$Sigma21Star[ind1,], S21 )
PASSED test at tolerance  1e-08
> test.for.zero( sparseObj$Sigma11Inv, solve(S11) )
PASSED test at tolerance  1e-08
> 
> Btest<- S21%*% solve( S11)
> look2<- spam2full( sparseObj$B)
> test.for.zero( Btest,look2[ind1, sIndex] )
PASSED test at tolerance  1e-08
> 
> 
> SEfull<- spam2full( sparseObj$SE)
> SE2full<- (SEfull)%*%t(SEfull)
> test.for.zero(diag( SE2full), sparseObj$predictionVariance )
PASSED test at tolerance  1e-08
> 
> SEtest<- t(chol(S22 - S21%*%solve( S11)%*%t( S21) ))
> test.for.zero(SEtest, SEfull[ind1, ind1] )
PASSED test at tolerance  1e-08
> 
> # check that debug FALSE also works
> 
> sparseObj1<-  offGridWeights( s, list( x= (1:m)*dx, y=(1:n)*dx),
+                                 aRange=(10*dx), sigma2=sigma2, 
+                                 Covariance="Exponential", 
+                                 np=2, 
+                                 debug=FALSE)
[1] 2
Found 2 grid box(es) containing more than 1 obs location
> 
> test.for.zero( sparseObj$B, sparseObj1$B)
PASSED test at tolerance  1e-08
> test.for.zero( sparseObj$SE, sparseObj1$SE)
PASSED test at tolerance  1e-08
> 
> cat("all done with off grid weight tests part 2", fill=TRUE)
all done with off grid weight tests part 2
> 
> 
> proc.time()
   user  system elapsed 
   0.51    0.15    0.60