Repository: equinor/GraphSPME (language: python)

Issue:
R crashes when `Z` is specifies independence through being identity matrix
This should be the easiest theoretical example


Discussion / hints:
This fails:
```
library(GraphSPME)
library(Matrix)

p <- 10
n <- 100
x <- matrix(rnorm(p*n), nrow=n, ncol=p)
Z <- bandSparse(p,p, 0:0, list(rep(1,p)))
prec <- prec_sparse(x,Z,T)
```
This fails because Z becomes object of type `dtCMatrix` but it needs to be of type `dgCMatrix`
This works:
```
Z <- bandSparse(p, p, 0, list(rep(1,p)))
Z <- as(Z, "dgCMatrix") # ensure correct sparsity class
Z
# 5 x 5 sparse Matrix of class "dgCMatrix"
# 
# [1,] 1 . . . .
# [2,] . 1 . . .
# [3,] . . 1 . .
# [4,] . . . 1 .
# [5,] . . . . 1
prec_est <- prec_sparse(x, Z)
prec_est
# 5 x 5 sparse Matrix of class "dgCMatrix"
# 
# [1,] 0.8224125 .        .        .        .       
# [2,] .         1.117901 .        .        .       
# [3,] .         .        1.132829 .        .       
# [4,] .         .        .        1.006042 .       
# [5,] .         .        .        .        1.065029
```
Thus the C++ code is okay.
R (and likely Python as well @sondreso ) should verify correct sparsity type on input for `prec_sparse`

Fix patch:
diff --git a/R-package/R/cov_shrink_spd.R b/R-package/R/cov_shrink_spd.R
index 5a03a8a..4cb8a79 100644
--- a/R-package/R/cov_shrink_spd.R
+++ b/R-package/R/cov_shrink_spd.R
@@ -1,6 +1,6 @@
 ## License: GPL-3
 
-# Stein-type covariance shrinkage estimate 
+# Stein-type covariance shrinkage estimate
 #
 
 #' Asymptotic Stein-type shrinkage of SPD covariance estimate.
@@ -10,9 +10,9 @@
 #' @param x an \eqn{nxp} matrix with \eqn{n-}observations of a \eqn{p-}vector
 #'
 #' @details
-#' 
+#'
 #' The standard frequentist covariance matrix estimate is calculated.
-#' Using techniques of Touloumis (2015) the matrix is shrunk towards a 
+#' Using techniques of Touloumis (2015) the matrix is shrunk towards a
 #' sparse target matrix, set to be the diagonal matrix of the frequentist
 #' estimate. The amount of shrinkage is found using asymptotic results
 #' found in the reference paper.
@@ -25,9 +25,9 @@
 #'
 #' @rdname cov_shrink_spd
 #' @export
-cov_shrink_spd <- function(x){
-    # check input
-    
-    # calculate cov
-    return(.cov_shrink_spd(x))
+cov_shrink_spd <- function(x) {
+  # check input
+
+  # calculate cov
+  return(.cov_shrink_spd(x))
 }
diff --git a/R-package/R/graphspme-package.R b/R-package/R/graphspme-package.R
index f0c39b5..3eab23c 100644
--- a/R-package/R/graphspme-package.R
+++ b/R-package/R/graphspme-package.R
@@ -3,14 +3,14 @@
 #' \code{GraphSPME} is a package for sparse precision matrix estimation through
 #' a graph.
 #' It employs the method of Lunde 2022, imported through the \code{graph_spme} function.
-#' 
+#'
 #' Important functions:
-#' 
+#'
 #' \itemize{
 #' \item \code{\link{prec_sparse}}: Sparse precision matrix estimation through a graph
 #' \item \code{\link{cov_shrink_spd}}: Asymptotic Stein-type shrinkage of SPD covariance estimate
 #' }
-#' 
+#'
 #' See individual function documentation for usage.
 #'
 #' @docType package
@@ -21,4 +21,4 @@
 #' @useDynLib GraphSPME, .registration = TRUE
 #' @import Rcpp RcppEigen
 #' @importFrom Rcpp evalCpp
-NULL
\ No newline at end of file
+NULL
diff --git a/R-package/R/prec_sparse.R b/R-package/R/prec_sparse.R
index fd88c27..f03b828 100644
--- a/R-package/R/prec_sparse.R
+++ b/R-package/R/prec_sparse.R
@@ -7,18 +7,18 @@
 #'
 #' \code{prec_sparse} is the precision matrix estimator in \pkg{GraphSPME}.
 #'
-#' @param x an \eqn{nxp} matrix with \eqn{n-}observations of a \eqn{p-}vector
+#' @param X an \eqn{nxp} matrix with \eqn{n-}observations of a \eqn{p-}vector
 #' @param Z a sparse matrix indicating non-zero elements of the precision matrix
 #' @param shrinkage if shrinkage should be applied to block-frequentist-covariance see details for more
 #'
 #' @details
-#' 
+#'
 #' The algorithm utilizes the knowledge of the sparsity provided through
 #' \code{Z}, using the algorithm of Le (2021). If \code{shrinkage=TRUE} then
-#' the method of Le (2021) is modified by the method in 
+#' the method of Le (2021) is modified by the method in
 #' Lunde (2022) to iteratively employ asymptotic shrinkage when building
 #' the sparse precision matrix estimator.
-#' 
+#'
 #' The algorithm should be robust in the face of extremely high dimensions.
 #'
 #' @return
@@ -29,9 +29,20 @@
 #'
 #' @rdname prec_sparse
 #' @export
-prec_sparse <- function(x, Z, shrinkage=TRUE){
-    # check input
-    
-    # calculate cov
-    return(.prec_sparse(x, Z, shrinkage))
+prec_sparse <- function(X, Z, shrinkage = TRUE) {
+  # check xtype
+  if (!is.matrix(X)) {
+    stop("X must be a matrix")
+  }
+  # dimensions
+  if (!all(ncol(X) == dim(Z))) {
+    stop("ncol(X) must equal both dim(Z)")
+  }
+  # sparsity type
+  if (!(class(Z) == "dgCMatrix")) {
+    stop("Z must be a sparse matrix of type dgCMatrix")
+  }
+
+  # calculate precision
+  return(.prec_sparse(X, Z, shrinkage))
 }


Test patch:
diff --git a/R-package/tests/testthat/test_graphspme.R b/R-package/tests/testthat/test_graphspme.R
index bf83e2a..2b595ab 100644
--- a/R-package/tests/testthat/test_graphspme.R
+++ b/R-package/tests/testthat/test_graphspme.R
@@ -1,6 +1,6 @@
 library(GraphSPME)
 
 test_that("cov_ml(X) is the same as R's cov(X)", {
-    X = matrix(rnorm(36), nrow=6)
-    expect_equal(.cov_ml(X), cov(X))
-})
\ No newline at end of file
+  X <- matrix(rnorm(36), nrow = 6)
+  expect_equal(.cov_ml(X), cov(X))
+})
diff --git a/R-package/tests/testthat/test_prec_sparse.R b/R-package/tests/testthat/test_prec_sparse.R
new file mode 100644
index 0000000..a19a075
--- /dev/null
+++ b/R-package/tests/testthat/test_prec_sparse.R
@@ -0,0 +1,34 @@
+# Tests for prec_sparse()
+
+library(GraphSPME)
+library(Matrix)
+
+p <- 10
+n <- 100
+X <- matrix(rnorm(p * n), nrow = n, ncol = p)
+Z <- bandSparse(p, p, 0:0, list(rep(1, p)))
+Z_wrong_sparsity <- as(Z, "dtCMatrix") # ensure wrong sparsity class
+Z_correct_sparsity <- as(Z, "dgCMatrix") # ensure correct sparsity class
+
+test_that("R throws an error for wrong sparsity type", {
+  expect_error(
+    prec_sparse(X, Z_wrong_sparsity),
+    "Z must be a sparse matrix of type dgCMatrix"
+  )
+})
+
+test_that("R throws an error for wrong data input type", {
+  X_df <- data.frame(x)
+  expect_error(
+    prec_sparse(X_df, Z_correct_sparsity),
+    "X must be a matrix"
+  )
+})
+
+test_that("R throws an error when data and neighbourhood structure does not match", {
+  X_extended <- cbind(X, rnorm(n))
+  expect_error(prec_sparse(X_extended, Z_correct_sparsity),
+    "ncol(X) must equal both dim(Z)",
+    fixed = TRUE
+  )
+})