version 1.0.2; change vecchia_cov package dependency from spGNNP to GpGp

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: bspme
 Type: Package
 Title: Bayesian Spatial Measurement Error Models
-Version: 1.0.1
+Version: 1.0.2
 Authors@R: c(person("Changwoo", "Lee", role=c("aut", "cre"), email="[email protected]"), person("Eun Sug", "Park", role = c("aut")))
 Author: Changwoo Lee[aut, cre], Eun Sug Park[aut]
 Maintainer: Changwoo Lee <[email protected]>
@@ -16,7 +16,7 @@ Imports:
     coda,
     fields,
     spam,
-    spNNGP
+    GpGp
 Depends: 
     Matrix,
     R (>= 2.10)

R/vecchia_cov.R

@@ -4,19 +4,22 @@
 #' this function finds a sparse precision matrix (inverse covariance) \eqn{Q} based on the Vecchia approximation (Vecchia 1988, Katzfuss and Guinness 2021),
 #' where \eqn{N(\mu, Q^{-1})} is the sparse MVN that approximates the original MVN \eqn{N(\mu, \Sigma)}.
 #' The algorithm is based on the pseudocode 2 of Finley et al. (2019).
+#' Nearest neighbor finding algorithm is based on GpGp::find_ordered_nn.
 #'
 #' @param Sigma *matrix<num>*, n by n covariance matrix
 #' @param coords *matrix<num>*, n by 2 coordinate matrix for nearest neighborhood search
 #' @param n.neighbors *integer*, the number of nearest neighbors (k) to determine conditioning set of Vecchia approximation
 #' @param ord *vector<int>*, length n vector, ordering of data. If NULL, ordering based on the first coordinate will be used.
 #' @param KLdiv *logical*, If TRUE, return KL divergence \eqn{D_{KL}(p || \tilde{p})} where \eqn{p} is multivariate normal with original covariance matrix and \eqn{\tilde{p}} is the approximated multivariate normal with sparse precision matrix.
+#' @param lonlat *logical*, (default FALSE) If TRUE, the coordinates are in longitude and latitude. If FALSE, the coordinates are in Euclidean space.
 #'
 #' @return list of the following:
 #' \describe{
 #'   \item{Q}{n by n sparse precision matrix in \link{Matrix} format}
 #'   \item{ord}{ordering used for Vecchia approximation}
 #'   \item{cputime}{time taken to run Vecchia approximation}
 #'   \item{KLdiv}{ (if \code{KLdiv = TRUE}) KL divergence \eqn{D_{KL}(p || \tilde{p})} where \eqn{p} is the multivariate normal with original covariance matrix and \eqn{\tilde{p}} is the approximated multivariate normal with a sparse precision matrix.}
+#'   \item{Linv}{ n by n sparse reverse Cholesky factor of Q, a lower triangular matrix such that Q = t(Linv)%*%Linv (before ordering changes). In other words, Linv = chol(Q[n:1,n:1])[n:1,n:1] (before ordering changes).}
 #' }
 #'
 #' @references Vecchia, A. V. (1988). Estimation and model identification for continuous spatial processes. Journal of the Royal Statistical Society Series B: Statistical Methodology, 50(2), 297-312.
@@ -25,23 +28,30 @@
 #'
 #' Finley, A. O., Datta, A., Cook, B. D., Morton, D. C., Andersen, H. E., & Banerjee, S. (2019). Efficient algorithms for Bayesian nearest neighbor Gaussian processes. Journal of Computational and Graphical Statistics, 28(2), 401-414.
 #'
-#' Zhang, L., (2020), public Github repository <https://github.com/LuZhangstat/NNGP_STAN>.
-#'
 #' @importFrom Matrix Diagonal
 #'
 #' @export
 #'
 #' @examples
 #'
+#' library(bspme)
 #' n = 1000
 #' coords = cbind(runif(n), runif(n))
+#' library(GpGp)
+#' ord <- GpGp::order_maxmin(coords) # from GpGp
 #' Sigma = fields::Exp.cov(coords, aRange = 1)
-#' fit5 = vecchia_cov(Sigma, coords, n.neighbors = 5, KLdiv = TRUE)
+#' fit5 = vecchia_cov(Sigma, coords, n.neighbors = 5, ord = ord, KLdiv = TRUE)
 #' fit5$KLdiv
-#' fit10 = vecchia_cov(Sigma, coords, n.neighbors = 10, KLdiv = TRUE)
+#' fit10 = vecchia_cov(Sigma, coords, n.neighbors = 10, ord = ord, KLdiv = TRUE)
 #' fit10$KLdiv
 #'
-vecchia_cov = function(Sigma, coords, n.neighbors, ord = NULL, KLdiv = FALSE){
+#' # Check Linv
+#' Q = fit5$Q
+#' Q_reordered = Q[ord,ord]
+#' all.equal(fit5$Linv, chol(Q_reordered[n:1,n:1])[n:1,n:1])
+#' all.equal(Q_reordered, t(fit5$Linv) %*% fit5$Linv)
+#'
+vecchia_cov = function(Sigma, coords, n.neighbors, ord = NULL, KLdiv = FALSE, lonlat = FALSE){
   start.time = Sys.time()
 
   n = ncol(Sigma)
@@ -55,16 +65,7 @@ vecchia_cov = function(Sigma, coords, n.neighbors, ord = NULL, KLdiv = FALSE){
       stop("error: supplied order vector ord must be of length n")
     }
   }
-
-  indx <- spNNGP:::mkNNIndxCB(coords, n.neighbors, n.omp.threads = 1) # search.type == "cb"
-  nn.indx <- indx$nnIndx
-  n.indx = spNNGP:::mk.n.indx.list(nn.indx, n, n.neighbors)
-
-  if(n.neighbors > 1){
-    NN_ind <- t(sapply(1: (n - 1), get_NN_ind, n.indx[-1], n.neighbors))
-  }else if(n.neighbors == 1){ # n.neighbors == 1
-    NN_ind <- as.matrix(sapply(1: (n - 1), get_NN_ind, n.indx[-1], n.neighbors))
-  }else{ # n.neighbors == 0
+  if(n.neighbors == 0){
     Q = Matrix::Diagonal(n, x=1/diag(Sigma))
     Q = as(Q, "dgCMatrix")
     out = list(Q = Q, ord = ord, cputime =  difftime(Sys.time(),start.time, units = "secs"))
@@ -73,17 +74,18 @@ vecchia_cov = function(Sigma, coords, n.neighbors, ord = NULL, KLdiv = FALSE){
     }
     return(out)
   }
+
   Sigma = Sigma[ord,ord] # Sigma is now reordered, overwrite to reduce memory usage
 
+  # GpGp::find_ordered_nn allows lon,lat input;
+  NNarray <- GpGp::find_ordered_nn(coords[ord,],n.neighbors, lonlat = lonlat)
+  NNarray = NNarray[,-1] # remove first column
   Nlist = list()
   Nlist[[1]] = c()
   for(i in 1:(n-1)){
-    if(i<=n.neighbors){
-      Nlist[[i+1]] <- NN_ind[i,1:i]
-    }else{
-      Nlist[[i+1]] <- NN_ind[i,]
-    }
+    Nlist[[i+1]] <- NNarray[i+1,!is.na(NNarray[i+1,])]
   }
+
   # sparse matrix
   A = Matrix::Matrix(0, n, n);
   D = numeric(n) #D = Matrix(0, n, n);
@@ -95,10 +97,13 @@ vecchia_cov = function(Sigma, coords, n.neighbors, ord = NULL, KLdiv = FALSE){
     A[i+1, Nlist[[i+1]]]= temp
     D[i+1]= Sigma[i+1, i+1] - sum(Sigma[i+1, Nlist[[i+1]]]*temp)
   }
-  D = Matrix::Diagonal(n, x = D)
-  Q_reordered = Matrix::t((Matrix::Diagonal(n) - A))%*%Matrix::solve(D)%*%(Matrix::Diagonal(n) - A)
+  Dmat = Matrix::Diagonal(n, x = D)
+  Q_reordered = Matrix::t((Matrix::Diagonal(n) - A))%*%Matrix::solve(Dmat)%*%(Matrix::Diagonal(n) - A)
+  Linv =  Matrix::Diagonal(n, x = 1/sqrt(D))%*%(Matrix::Diagonal(n) - A)
+
   Q = Q_reordered[order(ord),order(ord)]
   Q = as(Q, "dgCMatrix")
+
   end.time = Sys.time()
 
   out = list()
@@ -108,23 +113,8 @@ vecchia_cov = function(Sigma, coords, n.neighbors, ord = NULL, KLdiv = FALSE){
   if(KLdiv){
     out$KLdiv = as.numeric(0.5*(sum(Q*Sigma[order(ord),order(ord)]) - n - Matrix::determinant(Q)$modulus - determinant(Sigma)$modulus))
   }
+  out$Linv = Linv
   return(out)
 }
 
 
-
-
-
-
-#' helper function, get nearest neighbor indices
-#' source: https://github.com/LuZhangstat/NNGP_STAN/blob/master/NNMatrix.R
-#' @noRd
-get_NN_ind <- function (ind, ind_distM_i, M) {
-  if (ind < M ){l = ind } else {l = M};
-  D_i <- rep(0, M);
-  D_i[1:l] <- c(ind_distM_i[[ind]])[1:l]
-  return(D_i)
-}
-
-
-