a

R/getFunctions.R

@@ -38,12 +38,13 @@ getCovMatrix <- function(coco.object, type = "global", index = NULL){
 
   par.pos <- getDesignMatrix(coco.object@model.list, data = coco.object@data)$par.pos
 
+  theta_list <- cocons::getModelLists(coco.object@output$par,par.pos = par.pos, 
+                                      type = "diff")
+
   if(coco.object@type == "dense"){
 
     if(type == "global"){
-      theta_list <- cocons::getModelLists(coco.object@output$par,par.pos = par.pos, 
-                                         type = "diff")
-
+
       return(cocons::cov_rns(theta = theta_list,
                             locs = coco.object@locs,
                             x_covariates = x_covs,
@@ -59,10 +60,7 @@ getCovMatrix <- function(coco.object, type = "global", index = NULL){
   if(coco.object@type == "sparse"){
 
     if(type == "global"){
-
-      theta_list <- cocons::getModelLists(coco.object@output$par,par.pos = par.pos, 
-                                         type = "diff")
-
+
       # taper
       ref_taper <- coco.object@info$taper(
         spam::nearest.dist(coco.object@locs, delta = coco.object@info$delta, upper = NULL),
@@ -143,7 +141,7 @@ getDensityFromDelta <- function(coco.object, delta){
 
   coco.object@info$delta <- delta
 
-  return(summary(getCovMatrix(coco.object))$density)
+  return(summary(getCovMatrix(coco.object))$density / 100)
 
 }
 
@@ -698,7 +696,7 @@ getBoundariesV4 <- function(coco.object, lower.bound = 2, upper.bound = 2){
 #' @returns (\code{list}) a list with boundaries and simple init values for the optim L-BFGS-B routine
 #' @author Federico Blasi
 #' 
-getBoundaries <- function(x, lower.value, upper.value){
+getBoundaries <- function(x, lower.value = -2, upper.value = 2){
 
   if(upper.value < lower.value){stop("upper.value lower than lower.value")}
 

R/optim.R

@@ -37,7 +37,7 @@
 #'                    'smooth' = 3/2,
 #'                    'nugget' = -Inf)
 #'                    
-#' sample_index <- sample(1:dim(holes[[1]])[1],100)
+#' sample_index <- sample(1:dim(holes[[1]])[1],200)
 #'                    
 #' coco_object <- coco(type = 'dense',
 #'                     data = holes[[1]][sample_index, ],
@@ -46,7 +46,7 @@
 #'                     model.list = model.list)
 #'                     
 #' optim_coco <- cocoOptim(coco_object,
-#' boundaries = getBoundaries(coco_object)
+#' boundaries = getBoundaries(coco_object,-3,3))
 #' 
 #' plotOptimInfo(optim_coco)
 #' 
@@ -56,7 +56,7 @@
 #' 
 #' plot(optim_coco, type = 'correlations', index = c(2,3,5))
 #' 
-#' summary(optim_coco)
+#' summary(optim_coco,inv.hess = getHessian(optim_coco))
 #'  
 #' getEstims(optim_coco)
 #' 

README.md

@@ -23,6 +23,6 @@ For a brief introduction:
 
 or you can check the vignette directly [here](https://github.com/blasif/cocons/blob/main/inst/doc/cocons.pdf).
 
-For a description of the class of modular covariate-based covariance functions implemented in this package, please check the following [article](https://arxiv.org/abs/2410.16716#).
+**UPDATE!** Our article has been published in Environmetrics! Learn more about the type of modular and flexible Gaussian Process models you can fit here [article](https://doi.org/10.1002/env.70038)
 
 **Bugs and/or errors?:** feel free to add new issues on [https://github.com/blasif/cocons/issues](https://github.com/blasif/cocons/issues).

tests/coco_test.R

@@ -1,7 +1,17 @@
 # tests/coco_test.R
 
-# Load your package
-library(cocons)
+## --- Load package ----------------------------------------------------------
+
+suppressPackageStartupMessages({
+  library(cocons)
+})
+
+## Optional: check version at least 0.1.5
+if (packageVersion("cocons") < "0.1.5") {
+  stop("cocons >= 0.1.5 required for these checks.", call. = FALSE)
+}
+
+## --------
 
 data_test <- cocons::holes$training[1:50,]
 
@@ -190,3 +200,118 @@ stopifnot(all(eigen(cmat)$values > 0))
 
 predss <- cocoPredict(test_optim, newdataset = cocons::holes$test[1:200,1:4],newlocs = as.matrix(holes$test[1:200,1:2]))
 stopifnot(all(!is.na(predss$stochastic)))
+
+aeq <- function(a, b, tol = 1e-8) {
+  if (!isTRUE(all.equal(a, b, tolerance = tol))) {
+    stop(sprintf("Numeric mismatch:\nexpected: %s\ngot:      %s",
+                 paste0(head(b,10), collapse=", "),
+                 paste0(head(a,10), collapse=", ")),
+         call. = FALSE)
+  }
+  invisible(TRUE)
+}
+
+expect_error <- function(expr, pattern = NULL) {
+  err <- NULL
+  tryCatch(eval.parent(substitute(expr)), error = function(e) err <<- e)
+  if (is.null(err)) stop("Expected an error, but none was thrown.", call. = FALSE)
+  if (!is.null(pattern) && !grepl(pattern, conditionMessage(err))) {
+    stop(sprintf("Error did not match /%s/; got: %s",
+                 pattern, conditionMessage(err)), call. = FALSE)
+  }
+  invisible(TRUE)
+}
+
+set.seed(1)
+locs <- expand.grid(seq(0, 1, length.out = 5),
+                    seq(0, 1, length.out = 5))
+toydata <- data.frame(x = locs[, 1])
+z <- rnorm(nrow(locs))
+
+model.list <- list(
+  mean    = 0,
+  std.dev = formula(~ 1),
+  scale   = formula(~ 1 + x),
+  aniso   = 0,
+  tilt    = 0,
+  smooth  = 3/2,
+  nugget  = -Inf
+)
+
+coco_dense <- coco(
+  type = "dense",
+  data = toydata,
+  locs = as.matrix(locs),
+  z = z,
+  model.list = model.list
+)
+
+## S4 class + slots
+stopifnot(isS4(coco_dense))
+stopifnot(methods::is(coco_dense, "coco"))
+needed_slots <- c("type","data","locs","z","model.list","info","output")
+stopifnot(all(needed_slots %in% slotNames(coco_dense)))
+stopifnot(identical(coco_dense@type, "dense"))
+stopifnot(nrow(coco_dense@locs) == nrow(toydata))
+stopifnot(nrow(coco_dense@locs) == length(z))
+
+## --- 2) invalid 'type' should error ---------------------------------------
+
+expect_error(
+  coco(type = "weird",
+       data = toydata,
+       locs = as.matrix(locs),
+       z = z,
+       model.list = model.list),
+  pattern = "type|dense|sparse"
+)
+
+## --- 3) is.formula() helper ------------------------------------------------
+
+stopifnot(isTRUE( cocons::is.formula(~ 1) ))
+stopifnot(isTRUE(!cocons::is.formula(1)))
+stopifnot(isTRUE(!cocons::is.formula("~ 1 + x")))
+
+## --- 4) Scoring rules: getLogScore / getCRPS -------------------------------
+
+z.pred    <- c(0.0,  1.0, -0.5)
+mean.pred <- c(0.1,  0.8, -0.4)
+sd.pred   <- c(0.5,  1.2,  2.0)
+
+## reference formulas (match package code)
+logscore_ref <- (log(2*pi) + ((z.pred - mean.pred)/sd.pred)^2)/2 + log(sd.pred)
+v <- (mean.pred - z.pred)/sd.pred
+crps_ref <- sd.pred * ( v * (2 * stats::pnorm(v) - 1) +
+                          2 * stats::dnorm(v) - 1 / sqrt(pi) )
+
+aeq(cocons::getLogScore(z.pred, mean.pred, sd.pred), logscore_ref, tol = 1e-10)
+aeq(cocons::getCRPS(z.pred, mean.pred, sd.pred),    crps_ref,    tol = 1e-8)
+
+## --- 5) getDensityFromDelta() behavior ------------------------------------
+
+## Dense objects should error
+expect_error(cocons::getDensityFromDelta(coco_dense, delta = 0.2),
+             pattern = "only for sparse coco objects")
+
+## Optional: sparse path (skipped if 'spam' not available)
+if (requireNamespace("spam", quietly = TRUE)) {
+  coco_sparse <- coco(
+    type = "sparse",
+    data = toydata,
+    locs = as.matrix(locs),
+    z = z,
+    model.list = model.list,
+    info = list(
+      taper = spam::cov.wend1,
+      delta = 0.15,
+      smooth.limits = c(0.5, 2.5)
+    )
+  )
+
+  dens <- cocons::getDensityFromDelta(coco_sparse, delta = 0.20)
+  stopifnot(is.numeric(dens), length(dens) == 1L, is.finite(dens))
+  ## For a 25x25 grid with short taper, density should be in (0, 1)
+  stopifnot(dens > 0, dens < 1)
+} else {
+  message("Skipping sparse/taper density check (package 'spam' not available).")
+}