Pull request #28: Final linting fixes

Merge in STAT/prestogp from build-workflow to master

* commit '81b675dda9b57f2a362c34876364adf5bcf45560':
  Fix remaining linting issues
  Fix all indentation errors
  Add indentation_linter configuration to lintr config file

.lintr

@@ -3,6 +3,11 @@ linters: linters_with_defaults(
     commented_code_linter = NULL,
     object_name_linter = NULL,
     cyclocomp_linter = NULL,
-    object_length_linter = NULL
+    object_length_linter = NULL,
+    indentation_linter(
+      indent = 2L,
+      hanging_indent_style = "never",
+      assignment_as_infix = FALSE
+    )
   )
 encoding: "UTF-8"

R/Log_Likelihood.R

@@ -12,8 +12,7 @@
 #'
 #' @examples
 #' @noRd
-negloglik_vecchia_ST <- function(logparms, res, vecchia.approx, param.seq,
-                                 scaling, nscale) {
+negloglik_vecchia_ST <- function(logparms, res, vecchia.approx, param.seq, scaling, nscale) {
   parms <- unlog.params(logparms, param.seq, 1)
   locs.scaled <- vecchia.approx$locsord
   for (j in 1:nscale) {
@@ -141,8 +140,7 @@ mvnegloglik <- function(logparams, vecchia.approx, y, param.seq, P) {
 ##############################################################################
 ### Flexible Spatiotemporal Multivariate Matern Negative Loglikelihood Function ###########
 
-mvnegloglik_ST <- function(logparams, vecchia.approx, y, param.seq, P, scaling,
-                           nscale) {
+mvnegloglik_ST <- function(logparams, vecchia.approx, y, param.seq, P, scaling, nscale) {
   #  Input-
   #  logparams: A numeric vector of length (4*P)+(4*choose(P,2)).
   #             To construct these parameters we unlist a list of the 7 covariance
@@ -235,8 +233,7 @@ mvnegloglik.full <- function(logparams, locs, y, param.seq) {
 }
 
 ##############################################################################
-create.cov.upper.flex <- function(P, marg.var, marg.range, marg.smooth,
-                                  nugget, R.corr) {
+create.cov.upper.flex <- function(P, marg.var, marg.range, marg.smooth, nugget, R.corr) {
   # Create the symmetrical marginal+cross-covariance flexible matern from the
   # given parameters. Output is a list of the 4 Matern parameters as matrices
   sig2.mat <- diag(marg.var, P, P)
@@ -250,15 +247,12 @@ create.cov.upper.flex <- function(P, marg.var, marg.range, marg.smooth,
       j <- combs[iter, 2]
 
       smoothness.mat[i, j] <- (marg.smooth[i] + marg.smooth[j]) / 2
-      range.mat[i, j] <- 1 / sqrt(((1 / marg.range[i])^2 +
-        (1 / marg.range[j])^2) / 2)
+      range.mat[i, j] <- 1 / sqrt(((1 / marg.range[i])^2 + (1 / marg.range[j])^2) / 2)
 
       s1 <- sqrt(marg.var[i] * marg.var[j])
-      s2 <- ((1 / marg.range[i])^marg.smooth[i] *
-        (1 / marg.range[j])^marg.smooth[j]) /
+      s2 <- ((1 / marg.range[i])^marg.smooth[i] * (1 / marg.range[j])^marg.smooth[j]) /
         ((1 / range.mat[i, j])^(2 * smoothness.mat[i, j]))
-      s3 <- gamma(smoothness.mat[i, j]) / (sqrt(gamma(marg.smooth[i])) *
-        sqrt(gamma(marg.smooth[j])))
+      s3 <- gamma(smoothness.mat[i, j]) / (sqrt(gamma(marg.smooth[i])) * sqrt(gamma(marg.smooth[j])))
       s4 <- R.corr[iter]
       sig2.mat[i, j] <- s1 * s2 * s3 * s4
     }
@@ -340,8 +334,7 @@ cat.covariances <- function(locs.list, sig2, range, smoothness, nugget) {
 ##############################################################################
 ### Create the likelihood initial values                    #########
 
-create.initial.values.flex <- function(marg.var, marg.range, marg.smooth,
-                                       nugget, R.corr, P) {
+create.initial.values.flex <- function(marg.var, marg.range, marg.smooth, nugget, R.corr, P) {
   # Log-transform the covariance parameters and arrange in the proper order
   # for the likelihood function
   logparams.init <- c(

R/PrestoGP_CreateU_Multivariate.R

@@ -66,8 +66,7 @@ max_min_ordering <- function(locs, dist_func) {
 #' @param dist_func Any distance function with a signature of dist(query_location, locations_matrix)
 #'
 #' @return A vector containing the indices of the neighbors
-knn_indices <- function(ordered_locs, query, n_neighbors,
-                        dist_func, dist_func_code) {
+knn_indices <- function(ordered_locs, query, n_neighbors, dist_func, dist_func_code) {
   if (dist_func_code == "custom") {
     dists <- dist_func(query, ordered_locs)
     dists_order <- order(dists)
@@ -92,8 +91,7 @@ knn_indices <- function(ordered_locs, query, n_neighbors,
 #'
 #' @return A list containing two matrices, each with one row per location:
 #' an indices matrix with the indices of nearest neighbors for each location, and a distance matrix with the associated distances
-sparseNN <- function(ordered_locs, n_neighbors,
-                     dist_func, dist_func_code, ordered_locs_pred = NULL) {
+sparseNN <- function(ordered_locs, n_neighbors, dist_func, dist_func_code, ordered_locs_pred = NULL) {
   ee <- min(apply(ordered_locs, 2, stats::sd))
   n <- nrow(ordered_locs)
   ordered_locs <- ordered_locs + matrix(
@@ -112,8 +110,7 @@ sparseNN <- function(ordered_locs, n_neighbors,
   for (row in 1:n_neighbors) {
     # for the locations from 1 to n_neighbors, use the entire locs list to find the neighbors
     nn <- knn_indices(
-      ordered_locs[1:
-      (n_neighbors + 1), , drop = FALSE][-row, ,
+      ordered_locs[1:(n_neighbors + 1), , drop = FALSE][-row, ,
         drop = FALSE
       ],
       ordered_locs[row, , drop = FALSE], n_neighbors,
@@ -198,9 +195,9 @@ calc.q <- function(nn.obj, firstind.pred) {
 
 #' @export
 vecchia_Mspecify <- function(locs.list, m, locs.list.pred = NULL,
-                             dist.func = NULL,
-                             ordering.pred = c("obspred", "general"),
-                             pred.cond = c("independent", "general")) {
+  dist.func = NULL,
+  ordering.pred = c("obspred", "general"),
+  pred.cond = c("independent", "general")) {
   ordering.pred <- match.arg(ordering.pred)
   pred.cond <- match.arg(pred.cond)
 
@@ -365,8 +362,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) {
     ajj <- 1 / rangep[ondx[2]]
     aij <- sqrt((aii^2 + ajj^2) / 2)
     K1 <- rho.mat[ondx[1], ondx[2]] * sqrt(sig2[ondx[1]]) * sqrt(sig2[ondx[2]]) *
-      aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) *
-        gamma(vjj))) *
+      aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) *
       cov_func(dist_func(olocs[1, , drop = FALSE], olocs[2, , drop = FALSE], ),
         smoothness = vij, alpha = aij
       )
@@ -421,8 +417,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) {
     ajj <- 1 / rangep[ondx[2]]
     aij <- sqrt((aii^2 + ajj^2) / 2)
     K1 <- rho.mat[ondx[1], ondx[2]] * sqrt(sig2[ondx[1]]) * sqrt(sig2[ondx[2]]) *
-      aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) *
-        gamma(vjj))) *
+      aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) *
       cov_func(dist_func(olocs[1, , drop = FALSE], olocs[2, , drop = FALSE], ),
         smoothness = vij, alpha = aij
       )
@@ -461,8 +456,7 @@ createUMultivariate <- function(vec.approx, params, cov_func = NULL) {
         # positive definite. See equation (9) in Apanasovich (2011).
         K1[j] <- rho.mat[ondx[i], ondx[cur.q[j]]] *
           sqrt(sig2[ondx[i]]) * sqrt(sig2[ondx[cur.q[j]]]) *
-          aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) *
-            sqrt(gamma(vii) * gamma(vjj))) *
+          aii^vii * ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) *
           cov_func(
             dist_func(
               olocs[i, , drop = FALSE],

R/PrestoGP_Model.R

@@ -71,16 +71,16 @@ setGeneric("show_theta", function(object, Y_names) standardGeneric("show_theta")
 setGeneric(
   "prestogp_fit",
   function(model, Y, X, locs, scaling = NULL, apanasovich = FALSE,
-           covparams = NULL, beta.hat = NULL, tol = 0.999999, max_iters = 100, verbose = FALSE,
-           optim.method = "Nelder-Mead", optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000),
-           parallel = FALSE, foldid = NULL) {
+    covparams = NULL, beta.hat = NULL, tol = 0.999999, max_iters = 100, verbose = FALSE,
+    optim.method = "Nelder-Mead", optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000),
+    parallel = FALSE, foldid = NULL) {
     standardGeneric("prestogp_fit")
   }
 )
 setGeneric(
   "prestogp_predict",
   function(model, X = "matrix", locs = "matrix", m = "numeric", ordering.pred = c("obspred", "general"),
-           pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) {
+    pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) {
     standardGeneric("prestogp_predict")
   }
 )
@@ -215,10 +215,10 @@ setMethod(
 setMethod(
   "prestogp_fit", "PrestoGPModel",
   function(model, Y, X, locs, scaling = NULL, apanasovich = NULL,
-           covparams = NULL, beta.hat = NULL, tol = 0.999999,
-           max_iters = 100, verbose = FALSE, optim.method = "Nelder-Mead",
-           optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000),
-           parallel = FALSE, foldid = NULL) {
+    covparams = NULL, beta.hat = NULL, tol = 0.999999,
+    max_iters = 100, verbose = FALSE, optim.method = "Nelder-Mead",
+    optim.control = list(trace = 0, reltol = 1e-3, maxit = 5000),
+    parallel = FALSE, foldid = NULL) {
     model <- check_input(model, Y, X, locs)
     if (!is.null(beta.hat)) {
       if (!is.vector(beta.hat) | !is.numeric(beta.hat)) {
@@ -284,10 +284,7 @@ setMethod(
         parallel = parallel,
         foldid = foldid
       )
-      beta.hat <- as.matrix(predict(beta0.glmnet,
-        type = "coefficients",
-        s = beta0.glmnet$lambda.1se
-      ))
+      beta.hat <- as.matrix(predict(beta0.glmnet, type = "coefficients", s = beta0.glmnet$lambda.1se))
     }
     Y.hat <- beta.hat[1, 1] + model@X_train %*% beta.hat[-1, ]
 
@@ -521,8 +518,7 @@ setMethod("scale_locs", "PrestoGPModel", function(model, locs) {
       for (j in 1:model@nscale) {
         locs.out[[i]][, model@scaling == j] <-
           locs[[i]][, model@scaling == j] /
-            model@covparams[model@param_sequence[2, 1] +
-              model@nscale * (i - 1) + j - 1]
+            model@covparams[model@param_sequence[2, 1] + model@nscale * (i - 1) + j - 1]
       }
     }
     return(locs.out)
@@ -535,25 +531,20 @@ setMethod("transform_covariance_parameters", "PrestoGPModel", function(model) {
     model@covparams <- c(
       exp(model@logparams[1:model@param_sequence[2, 2]]),
       gtools::inv.logit(
-        model@logparams[model@param_sequence[3, 1]:
-        model@param_sequence[3, 2]],
+        model@logparams[model@param_sequence[3, 1]:model@param_sequence[3, 2]],
         0, 2.5
       ),
-      exp(model@logparams[model@param_sequence[4, 1]:
-      model@param_sequence[4, 2]]),
-      tanh(model@logparams[model@param_sequence[5, 1]:
-      model@param_sequence[5, 2]])
+      exp(model@logparams[model@param_sequence[4, 1]:model@param_sequence[4, 2]]),
+      tanh(model@logparams[model@param_sequence[5, 1]:model@param_sequence[5, 2]])
     )
   } else {
     model@covparams <- c(
       exp(model@logparams[1:model@param_sequence[2, 2]]),
       gtools::inv.logit(
-        model@logparams[model@param_sequence[3, 1]:
-        model@param_sequence[3, 2]],
+        model@logparams[model@param_sequence[3, 1]:model@param_sequence[3, 2]],
         0, 2.5
       ),
-      exp(model@logparams[model@param_sequence[4, 1]:
-      model@param_sequence[4, 2]]), 1
+      exp(model@logparams[model@param_sequence[4, 1]:model@param_sequence[4, 2]]), 1
     )
   }
   invisible(model)

R/PrestoGP_Multivariate_Vecchia.R

@@ -45,78 +45,76 @@ setMethod("initialize", "MultivariateVecchiaModel", function(.Object, n_neighbor
 #' prediction <- prestogp_predict(model, X.test, locs.test)
 #' Vec.mean <- prediction[[1]]
 #' Vec.sds <- prediction[[2]]
-setMethod("prestogp_predict",
-  "MultivariateVecchiaModel",
+setMethod("prestogp_predict", "MultivariateVecchiaModel",
   function(model, X, locs, m = NULL, ordering.pred = c("obspred", "general"), pred.cond = c("independent", "general"), return.values = c("mean", "meanvar")) {
-  # validate parameters
-  ordering.pred <- match.arg(ordering.pred)
-  pred.cond <- match.arg(pred.cond)
-  return.values <- match.arg(return.values)
-  X <- check_input_pred(model, X, locs)
-  if (is.null(m)) { # m defaults to the value used for training
-    m <- model@n_neighbors
-  }
-  if (m < model@min_m) {
-    stop(paste("m must be at least ", model@min_m, sep = ""))
-  }
-  if (m >= nrow(model@X_train)) {
-    warning("Conditioning set size m chosen to be >=n. Changing to m=n-1")
-    m <- nrow(model@X_train) - 1
-  }
-
-  # Vecchia prediction at new locations
-  Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx])
-  # Vecchia trend prediction at observed data
-  Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx])
+    # validate parameters
+    ordering.pred <- match.arg(ordering.pred)
+    pred.cond <- match.arg(pred.cond)
+    return.values <- match.arg(return.values)
+    X <- check_input_pred(model, X, locs)
+    if (is.null(m)) { # m defaults to the value used for training
+      m <- model@n_neighbors
+    }
+    if (m < model@min_m) {
+      stop(paste("m must be at least ", model@min_m, sep = ""))
+    }
+    if (m >= nrow(model@X_train)) {
+      warning("Conditioning set size m chosen to be >=n. Changing to m=n-1")
+      m <- nrow(model@X_train) - 1
+    }
 
-  # Test set prediction
-  res <- model@Y_train - Vecchia.hat
+    # Vecchia prediction at new locations
+    Vecchia.Pred <- predict(model@linear_model, newx = X, s = model@linear_model$lambda[model@lambda_1se_idx])
+    # Vecchia trend prediction at observed data
+    Vecchia.hat <- predict(model@linear_model, newx = model@X_train, s = model@linear_model$lambda[model@lambda_1se_idx])
 
-  locs.train.scaled <- scale_locs(model, model@locs_train)
-  locs.scaled <- scale_locs(model, locs)
-  vec.approx.test <- vecchia_Mspecify(locs.train.scaled, m,
-    locs.list.pred = locs.scaled,
-    ordering.pred = ordering.pred,
-    pred.cond = pred.cond
-  )
+    # Test set prediction
+    res <- model@Y_train - Vecchia.hat
 
-  ## carry out prediction
-  if (!model@apanasovich) {
-    params <- model@covparams
-    param.seq <- model@param_sequence
-    pred <- vecchia_Mprediction(res, vec.approx.test,
-      c(
-        params[1:param.seq[1, 2]],
-        rep(1, param.seq[2, 2] - param.seq[2, 1] + 1),
-        params[param.seq[3, 1]:
-        param.seq[5, 2]]
-      ),
-      return.values = return.values
-    )
-  } else {
-    pred <- vecchia_Mprediction(res, vec.approx.test, model@covparams,
-      return.values = return.values
+    locs.train.scaled <- scale_locs(model, model@locs_train)
+    locs.scaled <- scale_locs(model, locs)
+    vec.approx.test <- vecchia_Mspecify(locs.train.scaled, m,
+      locs.list.pred = locs.scaled,
+      ordering.pred = ordering.pred,
+      pred.cond = pred.cond
     )
-  }
 
-  # prediction function can return both mean and sds
-  # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord
-  Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend
-  if (return.values == "mean") {
-    return.list <- list(means = Vec.mean)
-  } else {
-    warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.")
-    vec.var <- pred$var.pred
-    ndx.out <- NULL
-    for (i in seq_along(length(locs))) {
-      vec.sds[ndx.out == i] <- sqrt(vec.var[ndx.out == i] +
-        model@covparams[model@param_sequence[4, i]])
+    ## carry out prediction
+    if (!model@apanasovich) {
+      params <- model@covparams
+      param.seq <- model@param_sequence
+      pred <- vecchia_Mprediction(res, vec.approx.test,
+        c(
+          params[1:param.seq[1, 2]],
+          rep(1, param.seq[2, 2] - param.seq[2, 1] + 1),
+          params[param.seq[3, 1]:param.seq[5, 2]]
+        ),
+        return.values = return.values
+      )
+    } else {
+      pred <- vecchia_Mprediction(res, vec.approx.test, model@covparams,
+        return.values = return.values
+      )
     }
-    return.list <- list(means = Vec.mean, sds = vec.sds)
-  }
 
-  return(return.list)
-})
+    # prediction function can return both mean and sds
+    # returns a list with elements mu.pred,mu.obs,var.pred,var.obs,V.ord
+    Vec.mean <- pred$mu.pred + Vecchia.Pred # residual + mean trend
+    if (return.values == "mean") {
+      return.list <- list(means = Vec.mean)
+    } else {
+      warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.")
+      vec.var <- pred$var.pred
+      ndx.out <- NULL
+      for (i in seq_along(length(locs))) {
+        vec.sds[ndx.out == i] <- sqrt(vec.var[ndx.out == i] + model@covparams[model@param_sequence[4, i]])
+      }
+      return.list <- list(means = Vec.mean, sds = vec.sds)
+    }
+
+    return(return.list)
+  }
+)
 
 setMethod("check_input", "MultivariateVecchiaModel", function(model, Y, X, locs) {
   if (!is.list(locs)) {
@@ -222,8 +220,7 @@ setMethod("specify", "MultivariateVecchiaModel", function(model) {
       for (j in 1:model@nscale) {
         olocs.scaled[model@vecchia_approx$ondx == i, model@scaling == j] <-
           olocs.scaled[model@vecchia_approx$ondx == i, model@scaling == j] *
-            model@covparams[model@param_sequence[2, 1] +
-              model@nscale * (i - 1) + j - 1]
+            model@covparams[model@param_sequence[2, 1] + model@nscale * (i - 1) + j - 1]
       }
     }
     model@vecchia_approx$locsord <- olocs.scaled
@@ -294,12 +291,8 @@ setMethod("transform_data", "MultivariateVecchiaModel", function(model, Y, X) {
       vecchia.approx = vecchia.approx,
       c(
         params[1:param.seq[1, 2]],
-        rep(
-          1,
-          param.seq[2, 2] - param.seq[2, 1] + 1
-        ),
-        params[param.seq[3, 1]:
-        param.seq[5, 2]]
+        rep(1, param.seq[2, 2] - param.seq[2, 1] + 1),
+        params[param.seq[3, 1]:param.seq[5, 2]]
       )
     )
   } else {