From 914f0e18153bd3c0123ce152ebb9b3f2da9ae95e Mon Sep 17 00:00:00 2001 From: Federico Blasi <33967974+blasif@users.noreply.github.com> Date: Thu, 21 Aug 2025 09:25:37 +0000 Subject: [PATCH 1/2] 0.1.5 : fix erorr in sparse density function, added unit tests, fixed example cocoOptim --- R/getFunctions.R | 16 +++--- R/optim.R | 6 +-- man/cocoOptim.Rd | 6 +-- tests/coco_test.R | 129 +++++++++++++++++++++++++++++++++++++++++++++- 4 files changed, 140 insertions(+), 17 deletions(-) diff --git a/R/getFunctions.R b/R/getFunctions.R index f863b31..fd6d94a 100644 --- a/R/getFunctions.R +++ b/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")} diff --git a/R/optim.R b/R/optim.R index d720cee..3d9a922 100644 --- a/R/optim.R +++ b/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) #' diff --git a/man/cocoOptim.Rd b/man/cocoOptim.Rd index 89fe87f..7fa5d48 100644 --- a/man/cocoOptim.Rd +++ b/man/cocoOptim.Rd @@ -40,7 +40,7 @@ model.list <- list('mean' = 0, '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, ], @@ -49,7 +49,7 @@ coco_object <- coco(type = 'dense', model.list = model.list) optim_coco <- cocoOptim(coco_object, -boundaries = getBoundaries(coco_object) +boundaries = getBoundaries(coco_object,-3,3)) plotOptimInfo(optim_coco) @@ -59,7 +59,7 @@ plot(optim_coco, type = 'ellipse') plot(optim_coco, type = 'correlations', index = c(2,3,5)) -summary(optim_coco) +summary(optim_coco,inv.hess = getHessian(optim_coco)) getEstims(optim_coco) diff --git a/tests/coco_test.R b/tests/coco_test.R index 81ed535..02c4a51 100644 --- a/tests/coco_test.R +++ b/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).") +} \ No newline at end of file From 74229283c7528f194f9a778e49cb3643e4e988a9 Mon Sep 17 00:00:00 2001 From: Federico Blasi <33967974+blasif@users.noreply.github.com> Date: Thu, 6 Nov 2025 21:42:07 +0100 Subject: [PATCH 2/2] Update article link to published work in Environmetrics Updated the article link to reflect publication in Environmetrics. --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 7bd16ce..0667f50 100644 --- a/README.md +++ b/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).