0.1.5 : several updates to cocons

DESCRIPTION

@@ -1,14 +1,14 @@
 Package: cocons
 Type: Package
 Title: Covariate-Based Covariance Functions for Nonstationary Spatial Modeling
-Version: 0.1.4
+Version: 0.1.5
 Author: Federico Blasi [aut, cre] (<https://orcid.org/0000-0001-9337-7154>), Reinhard Furrer [ctb] (<https://orcid.org/0000-0002-6319-2332>)
-Authors@R: c(person("Federico", "Blasi", email = "federico.blasi@gmail.com",
+Authors@R: c(person("Federico", "Blasi", email = "federicoblasi at gmail.com",
   role = c("aut", "cre"), comment = c(ORCID = "0000-0001-9337-7154")),
   person("Reinhard", "Furrer",
   role = c("ctb"), comment = c(ORCID = "0000-0002-6319-2332"))
   )
-Maintainer: Federico Blasi <federico.blasi@gmail.com>
+Maintainer: Federico Blasi <federicoblasi@gmail.com>
 Description: Estimation, prediction, and simulation of nonstationary Gaussian process with modular covariate-based covariance functions. 
   Sources of nonstationarity, such as spatial mean, variance, geometric anisotropy, smoothness, and nugget, can be considered based on spatial characteristics. 
   An induced compact-supported nonstationary covariance function is provided, enabling fast and memory-efficient computations when handling densely sampled domains.

NEWS.md

@@ -1,3 +1,10 @@
+## cocons 0.1.5
+
+- Vignette updated
+- changed default max.iters from 200 to 500.
+- now `coco` objects store information about optim.control arguments for full reproducibility.
+- introduce a two-step optimization approach for automatic model selection. See `help(coco)`.
+
 ## cocons 0.1.4
 
 - Documentation polished

R/RcppExports.R

@@ -1,6 +1,16 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+#' smoothed-L1 penalization over the covariate-driven parameters
+#'
+#' @param x a vector of parameters to penalize
+#' @param lambda penalization hyperparameter
+#' @param alpha hyperparameter to approximate L1 penalization. By default 1e6 provides a decent approximation to L1, while also being a differentiable function
+#' @return penalization
+sumsmoothlone <- function(x, lambda, alpha = 1e6) {
+    .Call(`_cocons_sumsmoothlone`, x, lambda, alpha)
+}
+
 #' Dense covariance function (difference parameterization)
 #'
 #' @param theta vector of parameters

R/checkFunctions.R

@@ -17,6 +17,10 @@
     optim.control$control$factr <- to_update$control$factr
   }
 
+  if(is.null(optim.control$control$ndeps)){
+    optim.control$control$ndeps <- to_update$control$ndeps
+  }
+
   if(is.null(optim.control$control$trace)){
     optim.control$control$trace <- to_update$control$trace
   }
@@ -269,9 +273,21 @@
     stop("smooth limits not specified.")
   }
 
-  if (!is.null(info$lambda)) {
-    if(info$lambda < 0){
-      stop("lambda must be non-negative.")
+  if (!is.null(info$lambda.reg)) {
+    if(info$lambda.reg < 0){
+      stop("lambda.reg must be non-negative.")
+    }
+  }
+
+  if (!is.null(info$lambda.Sigma)) {
+    if(info$lambda.Sigma < 0){
+      stop("lambda.Sigma must be non-negative.")
+    }
+  }
+
+  if (!is.null(info$lambda.betas)) {
+    if(info$lambda.betas < 0){
+      stop("lambda.betas must be non-negative.")
     }
   }
 
@@ -451,11 +467,21 @@
 
 .cocons.getPen <- function(n, lambda, theta_list, smooth.limits){
 
-  return(2 * n * lambda * exp(theta_list$scale[1]) * 
-           sqrt(((smooth.limits[2]-smooth.limits[1])/ 
-                   (1 + exp(-theta_list$smooth[1])) + 
-                   smooth.limits[1]))
-  )
+  if(T){
+    summ <- lambda[3] * exp(theta_list$scale[1]) * 
+      sqrt(((smooth.limits[2]-smooth.limits[1])/ 
+              (1 + exp(-theta_list$smooth[1])) + 
+              smooth.limits[1])) + sumsmoothlone(x = unlist(theta_list[[1]][-1]), lambda[2])    
+  }
+
+  for(ii in 2:6){
+
+    lambda_i <- sumsmoothlone(x = unlist(theta_list[[ii]][-1]), lambda[1])
+    summ <- summ + lambda_i
+
+  }
+
+  return(2 * n * summ)
 
 }
 
@@ -478,3 +504,98 @@
   return(namesToUpdate)
 
 }
+
+
+.cocons.update.coco.first.step <- function(coco.object, output, boundaries){
+
+  # Update model structure
+
+  par_pos <- getDesignMatrix(coco.object@model.list,data = coco.object@data)$par.pos
+
+  number_parameters_models <- lapply(par_pos,sum)
+
+  new_formulas_list <- coco.object@model.list
+
+  coco.object@output <- output
+
+  parss <- getEstims(coco.object)
+
+  # What if one estimate is exactly 0 ?
+  # Check models formulas keeping environment
+
+  for(ii in 1:7){
+
+    if(!is.formula(coco.object@model.list[[ii]])){
+      next
+    }
+
+    to_zero <- which(abs(parss[[ii]][par_pos[[ii]]]) <= coco.object@info$sparse.point)
+
+    if(length(to_zero) == 0){
+      next
+    }
+
+    to_zero_covs <- sum(abs(parss[[ii]][par_pos[[ii]]]) <= coco.object@info$sparse.point)
+
+    if( to_zero_covs == (number_parameters_models[[ii]] - 1) || (to_zero_covs == (number_parameters_models[[ii]]))){
+      new_formulas_list[[ii]] <- stats::as.formula("~1", env = globalenv())
+      next
+    }
+
+    # Update formula
+
+    # If the only parameter to drop is the intercept
+    if(length(to_zero) == 1 && to_zero == 1){
+      next
+    }
+
+    new_terms <- stats::drop.terms(stats::terms(coco.object@model.list[[ii]]), dropx = (to_zero - 1), keep.response = TRUE)
+
+    if(attr(new_terms,"intercept") & (length(attr(new_terms,"term.labels")) > 0)){
+      new_terms <- c("1",attr(new_terms, "term.labels"))
+    } else{
+      new_terms <- attr(new_terms, "term.labels")
+    }
+
+    new_formula <- stats::reformulate(new_terms, response = NULL, env = globalenv())
+    new_formulas_list[[ii]] <- new_formula
+
+  }
+
+  index_where <- 0
+
+  vector_logical <- logical(length = length(output$par))
+
+  for(ii in 1:7){
+
+    if(!is.formula(coco.object@model.list[[ii]])){
+      next
+    }
+
+    to_zero <- which(abs(parss[[ii]][par_pos[[ii]]]) <= coco.object@info$sparse.point)
+
+    if(any(to_zero == 1) && (length(to_zero) > 1)){
+      to_zero <- to_zero[-1]
+    }
+
+    if(length(to_zero) == 1 && to_zero == 1){
+      index_where <- index_where + number_parameters_models[[ii]]
+      next
+    }
+
+    index_to_place <- (index_where + 1 ) : (index_where + number_parameters_models[[ii]])
+
+    vector_logical[index_to_place][to_zero] <- TRUE
+
+    index_where <- index_where + number_parameters_models[[ii]]
+
+  }
+
+  boundaries_new <- lapply(boundaries, function(x) x[!vector_logical])
+
+  coco.object@info$boundaries <- boundaries_new
+
+  coco.object@model.list <- new_formulas_list
+
+  return(coco.object)
+}

R/cocons.R

@@ -26,7 +26,10 @@
 #' \eqn{\nu_{l}}, and \eqn{\nu_{u}} are the lower and upper bounds limiting the range of variation of the spatially-varying smoothness, and where \eqn{\boldsymbol{X}_{\ell}} relates to the design matrix defined by the specific models for each of the source of nonstationarity.
 #' 
 #' Lastly, arguments for the \code{"info"} list argument involve: \itemize{
-#' \item \code{"lambda"}: (\code{numeric}) a positive scalar specifying the regularization parameter. Larger values penalizes highly-smoothed long-tailed covariance functions. 
+#' \item \code{"lambda.reg"}: (\code{numeric}) a positive scalar specifying the regularization parameter. Larger values regularize highly-smoothed long-tailed covariance functions. 
+#' \item \code{"lambda.Sigma"}: (\code{numeric}) a positive scalar specifying the penalization parameter for the covariate-driven covariance parameters. 
+#' \item \code{"lambda.betas"}: (\code{numeric}) a positive scalar specifying the penalization parameter for the covariate-driven spatial mean parameters. 
+#' \item \code{"sparse.point"}: (\code{numeric}) a cutting point for which smaller coefficients in absolute value will be set to zero after the smoothed L1 penalization optimization. Used in combination with \code{lambda.Sigma} and \code{lambda.betas}. By default, it is set to \code{1e-4}.
 #' \item \code{"smooth.limits"}: (\code{numeric vector}) specifying the range of variation for the spatially varying smoothness (e.g. c(0.5, 2.5)).
 #' \item \code{"taper"}: (\code{numeric}) specifying the desired taper function from the spam package (only for "sparse" coco objects).
 #' \item \code{"delta"}: (\code{numeric}) specifying the taper range/scale (only for "sparse" coco objects).
@@ -135,8 +138,20 @@ coco <- function(type,
 
   .cocons.check.output(output)
 
-  if(is.null(info$lambda)){
-    info$lambda <- 0
+  if(is.null(info$lambda.reg)){
+    info$lambda.reg <- 0
+  }
+
+  if(is.null(info$lambda.betas)){
+    info$lambda.betas <- 0
+  }
+
+  if(is.null(info$lambda.Sigma)){
+    info$lambda.Sigma <- 0
+  }
+
+  if(any(c(info$lambda.betas,info$lambda.Sigma) > 0) && (is.null(info$sparse.point))){
+    info$sparse.point <- 1e-4
   }
 
   if(T){

R/getFunctions.R

@@ -617,6 +617,71 @@ getModelLists <- function(theta, par.pos, type = "diff"){
   }
 }
 
+#' Simple build of boundaries
+#' @description provides a set of upper and lower bounds for the L-BFGS-B routine
+#' 
+#' @usage getBoundariesV4(coco.object, lower.bound, upper.bound)
+#' @param coco.object \code{(S4)} a coco.object.
+#' @param lower.bound \code{(numeric vector)} lower bound for covariate-driven parameters of the nonstationary covariance function.
+#' @param upper.bound \code{(numeric vector)} upper bound for covariate-driven parameters of the nonstationary covariance function.
+#' @returns (\code{list}) a list with boundaries and simple init values for the optim L-BFGS-B routine
+#' @author Federico Blasi
+#' 
+getBoundariesV4 <- function(coco.object, lower.bound = 2, upper.bound = 2){
+
+  designMatrix <- getDesignMatrix(model.list = coco.object@model.list, data = coco.object@data)
+
+  term_labels <- paste("z", "~" , paste(coco.object@model.list$mean)[2])
+
+  # Skip scaling
+  if(any(colnames(coco.object@data[,coco.object@info$skip.scale]) %in% 
+         colnames(coco.object@data))){
+    tmp_info <- .cocons.setDesignMatrixCat(coco.object, designMatrix = designMatrix)
+    tmp_values <- tmp_info$tmp_values
+    mod_DM <- tmp_info$mod_DM
+  } else {
+    tmp_values <- cocons::getScale(designMatrix$model.matrix)
+    mod_DM <- tmp_values$std.covs      
+  }
+
+  coefs_lm <- stats::coef(stats::lm(term_labels,cbind.data.frame('z' = coco.object@z,mod_DM)))
+
+  boundaries_B <- getBoundariesV2(coco.object = coco.object,
+                                  mean.limits = c(-Inf, 0, Inf),
+                                  std.dev.limits = c(-lower.bound, 0, upper.bound),
+                                  scale.limits = c(-lower.bound, 0, upper.bound),
+                                  aniso.limits =  c(-lower.bound, 0, upper.bound),
+                                  tilt.limits =  c(-lower.bound, 0, upper.bound),
+                                  smooth.limits = c(-lower.bound, 0, upper.bound),
+                                  nugget.limits = c(-lower.bound, 0, upper.bound))
+
+  boundaries_B$theta_init[1:length(coefs_lm)] <- coefs_lm
+
+  first_sd <- which(names(boundaries_B$theta_init) == "std.dev.limits")[1]
+  n_var <- length(which(names(boundaries_B$theta_init) == "std.dev.limits")) - 1
+
+  first_range <- which(names(boundaries_B$theta_init) == "scale.limits")[1]
+  n_range <- length(which(names(boundaries_B$theta_init) == "scale.limits")) - 1
+
+  first_smooth <- which(names(boundaries_B$theta_init) == "smooth.limits")[1]
+  n_smooth <- length(which(names(boundaries_B$theta_init) == "smooth.limits")) - 1
+
+  boundaries_B$theta_init[first_range] <- (log(stats::sd(coco.object@z)) - log(stats::sd(c(stats::dist(coco.object@locs)))))/2
+  boundaries_B$theta_init[first_sd] <- (log(stats::sd(coco.object@z)) + log(stats::sd(c(stats::dist(coco.object@locs)))))/2
+
+  boundaries_B$theta_upper[c(first_sd, first_range)] <- c(3 * abs(boundaries_B$theta_init[first_sd]), 3 * abs(boundaries_B$theta_init[first_range]))
+  boundaries_B$theta_lower[c(first_sd, first_range)] <- c(-3 * abs(boundaries_B$theta_init[first_sd]), -3 * abs(boundaries_B$theta_init[first_range]))
+
+  boundaries_B$theta_upper[first_smooth] <- 3
+  boundaries_B$theta_lower[first_smooth] <- -3.5
+
+  boundaries_B$theta_init[1] <- mean(coco.object@z)
+  boundaries_B$theta_upper[1] <- boundaries_B$theta_init[1] + 3 * boundaries_B$theta_init[1]
+  boundaries_B$theta_lower[1] <- boundaries_B$theta_init[1] - 3 * boundaries_B$theta_init[1]
+
+  return(boundaries_B)
+}
+
 #' Simple build of boundaries
 #' @description provides a generic set of upper and lower bounds for the L-BFGS-B routine
 #' 
@@ -894,7 +959,7 @@ getHessian <- function(coco.object, ncores = parallel::detectCores() - 1,
                                    smooth.limits = coco.object@info$smooth.limits, 
                                    n = n, 
                                    z = coco.object@z,
-                                   lambda = coco.object@info$lambda)
+                                   lambda = c(coco.object@info$lambda.Sigma, coco.object@info$lambda.betas, coco.object@info$lambda.reg))
 
     } else{f00 <- coco.object@output$value}
 
@@ -925,23 +990,23 @@ getHessian <- function(coco.object, ncores = parallel::detectCores() - 1,
                                                                                     smooth.limits = coco.object@info$smooth.limits, 
                                                                                     n = n, 
                                                                                     z = coco.object@z,
-                                                                                    lambda = coco.object@info$lambda)
+                                                                                    lambda = c(coco.object@info$lambda.Sigma, coco.object@info$lambda.betas, coco.object@info$lambda.reg))
 
                                              f10 <- cocons::GetNeg2loglikelihood(t10, par.pos = par.pos,
                                                                                     locs = coco.object@locs,
                                                                                     x_covariates = x_covariates, 
                                                                                     smooth.limits = coco.object@info$smooth.limits, 
                                                                                     n = n, 
                                                                                     z = coco.object@z,
-                                                                                    lambda = coco.object@info$lambda)
+                                                                                    lambda = c(coco.object@info$lambda.Sigma, coco.object@info$lambda.betas, coco.object@info$lambda.reg))
 
                                              f11 <- cocons::GetNeg2loglikelihood(t11, par.pos = par.pos,
                                                                                     locs = coco.object@locs,
                                                                                     x_covariates = x_covariates, 
                                                                                     smooth.limits = coco.object@info$smooth.limits, 
                                                                                     n = n, 
                                                                                     z = coco.object@z, 
-                                                                                    lambda = coco.object@info$lambda)
+                                                                                    lambda = c(coco.object@info$lambda.Sigma, coco.object@info$lambda.betas, coco.object@info$lambda.reg))
 
                                              0.5 * ((f11 - f01 - f10 + f00) / (eps*eps))
                                            })
@@ -1009,7 +1074,7 @@ getHessian <- function(coco.object, ncores = parallel::detectCores() - 1,
                                                cholS = cholS,
                                                z = coco.object@z,
                                                n = n,
-                                               lambda = coco.object@info$lambda)
+                                               lambda = c(coco.object@info$lambda.Sigma, coco.object@info$lambda.betas, coco.object@info$lambda.reg))
 
     } else{f00 <- coco.object@output$value}
 
@@ -1044,7 +1109,7 @@ getHessian <- function(coco.object, ncores = parallel::detectCores() - 1,
                                                                                           cholS = cholS,
                                                                                           z = coco.object@z,
                                                                                           n = n,
-                                                                                          lambda = coco.object@info$lambda)
+                                                                                          lambda = c(coco.object@info$lambda.Sigma, coco.object@info$lambda.betas, coco.object@info$lambda.reg))
 
                                              f10 <- cocons::GetNeg2loglikelihoodTaper(theta = t10,
                                                                                           par.pos = par.pos,
@@ -1055,7 +1120,7 @@ getHessian <- function(coco.object, ncores = parallel::detectCores() - 1,
                                                                                           cholS = cholS,
                                                                                           z = coco.object@z,
                                                                                           n = n,
-                                                                                          lambda = coco.object@info$lambda)
+                                                                                          lambda = c(coco.object@info$lambda.Sigma, coco.object@info$lambda.betas, coco.object@info$lambda.reg))
 
                                              f11 <- cocons::GetNeg2loglikelihoodTaper(theta = t11,
                                                                                           par.pos = par.pos,
@@ -1066,7 +1131,7 @@ getHessian <- function(coco.object, ncores = parallel::detectCores() - 1,
                                                                                           cholS = cholS,
                                                                                           z = coco.object@z,
                                                                                           n = n,
-                                                                                          lambda = coco.object@info$lambda)
+                                                                                          lambda = c(coco.object@info$lambda.Sigma, coco.object@info$lambda.betas, coco.object@info$lambda.reg))
 
                                              0.5 * ((f11 - f01 - f10 + f00) / (eps*eps))