Skip to content

Commit

Permalink
Merge branch 'build-workflow' of sciome-bot-git:Spatiotemporal-Exposu…
Browse files Browse the repository at this point in the history
…res-and-Toxicology/PrestoGP into build-workflow
  • Loading branch information
sciome-bot committed Jan 12, 2024
2 parents 245eaa3 + 43ae272 commit 0460202
Show file tree
Hide file tree
Showing 24 changed files with 1,150 additions and 813 deletions.
3 changes: 2 additions & 1 deletion .Rbuildignore
Original file line number Diff line number Diff line change
Expand Up @@ -5,4 +5,5 @@
^\.Rproj\.user$
^doc$
^Meta$
^\.github/
^\.lintr
^\.github
5 changes: 2 additions & 3 deletions .github/workflows/lint.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
on:
push:
branches: [main, master]
branches: [main, master, build-workflow]
pull_request:
branches: [main, master]

Expand All @@ -26,8 +26,7 @@ jobs:
needs: lint

- name: Lint
run: lintr::lint_package(linters = lintr::linters_with_defaults(object_name_linter = NULL,
commented_code_linter = NULL, cyclocomp_linter = NULL))
run: lintr::lint_package()
shell: Rscript {0}
env:
LINTR_ERROR_ON_LINT: true
6 changes: 6 additions & 0 deletions .lintr
Original file line number Diff line number Diff line change
@@ -0,0 +1,6 @@
linters: linters_with_defaults(
line_length_linter = NULL,
commented_code_linter = NULL,
object_name_linter = NULL,
cyclocomp_linter = NULL)
encoding: "UTF-8"
40 changes: 26 additions & 14 deletions R/Log_Likelihood.R
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,13 @@ negloglik_vecchia_ST <- function(logparms, res, vecchia.approx, param.seq,
parms[param.seq[2, 1] + j - 1]
}
vecchia.approx$locsord <- locs.scaled
-vecchia_likelihood(res, vecchia.approx, c(parms[1], 1,
parms[param.seq[3, 1]]),
parms[param.seq[4, 1]])
-vecchia_likelihood(
res, vecchia.approx, c(
parms[1], 1,
parms[param.seq[3, 1]]
),
parms[param.seq[4, 1]]
)
}

#' negloglik_vecchia
Expand All @@ -42,8 +46,10 @@ negloglik_vecchia_ST <- function(logparms, res, vecchia.approx, param.seq,
#' @noRd
negloglik_vecchia <- function(logparms, res, vecchia.approx, param.seq) {
parms <- unlog.params(logparms, param.seq, 1)
-vecchia_likelihood(res, vecchia.approx, c(parms[1], parms[2], parms[3]),
parms[4])
-vecchia_likelihood(
res, vecchia.approx, c(parms[1], parms[2], parms[3]),
parms[4]
)
}

#' negloglik_full_ST
Expand Down Expand Up @@ -95,8 +101,10 @@ negloglik.full <- function(logparams, d, y, param.seq) {
params <- unlog.params(logparams, param.seq, 1)
# d <- fields::rdist(locs)
N <- nrow(d)
cov.mat <- params[1] * fields::Matern(d, range = params[2],
smoothness = params[3]) +
cov.mat <- params[1] * fields::Matern(d,
range = params[2],
smoothness = params[3]
) +
params[4] * diag(N)
return(-1 * mvtnorm::dmvnorm(y, rep(0, N), cov.mat, log = TRUE))
}
Expand Down Expand Up @@ -243,14 +251,14 @@ create.cov.upper.flex <- function(P, marg.var, marg.range, marg.smooth,

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)
(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]) /
(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])))
sqrt(gamma(marg.smooth[j])))
s4 <- R.corr[iter]
sig2.mat[i, j] <- s1 * s2 * s3 * s4
}
Expand Down Expand Up @@ -292,12 +300,16 @@ cat.covariances <- function(locs.list, sig2, range, smoothness, nugget) {
# Calculate the covariance matrix - if/then based on its location in the super-matrix
N <- nrow(d)
if (i == j) { # To accomodate varying size outcomes- the nugget is not included on cross-covariances
cov.mat.ij <- sig2[i, j] * geoR::matern(d, phi = range[i, j], kappa =
smoothness[i, j]) +
cov.mat.ij <- sig2[i, j] * geoR::matern(d,
phi = range[i, j], kappa =
smoothness[i, j]
) +
nugget[i, j] * diag(N)
} else {
cov.mat.ij <- sig2[i, j] * geoR::matern(d, phi = range[i, j], kappa =
smoothness[i, j])
cov.mat.ij <- sig2[i, j] * geoR::matern(d,
phi = range[i, j], kappa =
smoothness[i, j]
)
}


Expand Down
39 changes: 24 additions & 15 deletions R/PrestoGP_CreateU_Multivariate.R
Original file line number Diff line number Diff line change
Expand Up @@ -100,25 +100,34 @@ sparseNN <- function(ordered_locs, n_neighbors,
stats::rnorm(n * ncol(ordered_locs)),
n, ncol(ordered_locs)
)
indices_matrix <- matrix(data = NA, nrow = nrow(ordered_locs),
ncol = n_neighbors)
distances_matrix <- matrix(data = NA, nrow = nrow(ordered_locs),
ncol = n_neighbors)
indices_matrix <- matrix(
data = NA, nrow = nrow(ordered_locs),
ncol = n_neighbors
)
distances_matrix <- matrix(
data = NA, nrow = nrow(ordered_locs),
ncol = 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, ,
drop = FALSE],
ordered_locs[row, , drop = FALSE], n_neighbors,
dist_func, dist_func_code)
nn <- knn_indices(
ordered_locs[1:
(n_neighbors + 1), , drop = FALSE][-row, ,
drop = FALSE
],
ordered_locs[row, , drop = FALSE], n_neighbors,
dist_func, dist_func_code
)
indices_matrix[row, 1:n_neighbors] <- nn$indices[1:n_neighbors]
distances_matrix[row, 1:n_neighbors] <- nn$distances[1:n_neighbors]
}
for (row in (n_neighbors + 1):nrow(ordered_locs)) {
# get the m nearest neighbors from the locs before this one in the max-min order
nn <- knn_indices(ordered_locs[1:(row - 1), , drop = FALSE],
ordered_locs[row, , drop = FALSE], n_neighbors,
dist_func, dist_func_code)
nn <- knn_indices(
ordered_locs[1:(row - 1), , drop = FALSE],
ordered_locs[row, , drop = FALSE], n_neighbors,
dist_func, dist_func_code
)
indices_matrix[row, 1:n_neighbors] <- nn$indices[1:n_neighbors]
distances_matrix[row, 1:n_neighbors] <- nn$distances[1:n_neighbors]
}
Expand Down Expand Up @@ -170,7 +179,7 @@ calc.q <- function(nn.obj, firstind.pred) {
for (j in 2:m) {
cur.k <- cur.q[j]
cur.qy <- intersect(q.y[[cur.k]], cur.q)
if (length(cur.qy) > length(best.qy) & cur.k < firstind.pred) {
if (length(cur.qy) > length(best.qy) && cur.k < firstind.pred) {
best.k <- cur.k
best.qy <- cur.qy
}
Expand Down Expand Up @@ -245,7 +254,7 @@ vecchia_Mspecify <- function(locs.list, m, locs.list.pred = NULL,
loc.order <- max_min_ordering(locs.all, dist.func)
loc.order <- c(unique(loc.order), setdiff(1:n, loc.order))
} else {
if (is.null(locs.list.pred) | ordering.pred == "general") {
if (is.null(locs.list.pred) || ordering.pred == "general") {
loc.order <- GPvecchia::order_maxmin_exact(locs.all)
# I am not sure why the next two lines are here. I added them because
# similar code exists in the GPvecchia package. But I don't know why
Expand Down Expand Up @@ -274,7 +283,7 @@ vecchia_Mspecify <- function(locs.list, m, locs.list.pred = NULL,
# is non-deterministic, so there may be some slight differences
# between the output of this function and the output of createU
# in the GPvecchia package.
if (is.null(locs.list.pred) | pred.cond == "general") {
if (is.null(locs.list.pred) || pred.cond == "general") {
nn.mat <- sparseNN(olocs, m, dist.func, dist.func.code)
} else {
nn.mat <- sparseNN(
Expand Down
9 changes: 4 additions & 5 deletions R/PrestoGP_Multivariate_Vecchia.R
Original file line number Diff line number Diff line change
Expand Up @@ -199,13 +199,12 @@ setMethod("check_input_pred", "MultivariateVecchiaModel", function(model, X, loc
}
}
if (length(X) == 1) {
X <- X[[1]]
}
else {
X <- psych::superMatrix(X)
X <- X[[1]]
} else {
X <- psych::superMatrix(X)
}
if (ncol(X) != ncol(model@X_train)) {
stop("X and X_train must have the same number of predictors")
stop("X and X_train must have the same number of predictors")
}
return(X)
})
Expand Down
4 changes: 2 additions & 2 deletions R/PrestoGP_Util_Functions.R
Original file line number Diff line number Diff line change
Expand Up @@ -233,11 +233,11 @@ vecchia_Mprediction <- function(z, vecchia.approx, covparms, var.exact = NULL,
mu.pred = vecchia.mean$mu.pred, mu.obs = vecchia.mean$mu.obs,
var.pred = NULL, var.obs = NULL, V.ord = NULL, U.obj = NULL
)
if (return.values == "meanmat" | return.values == "all") {
if (return.values == "meanmat" || return.values == "all") {
return.list$V.ord <- V.ord
return.list$U.obj <- U.obj
}
if (return.values == "meanvar" | return.values == "all") {
if (return.values == "meanvar" || return.values == "all") {
if (is.null(var.exact)) {
var.exact <- (sum(!vecchia.approx$obs) < 4 * 10000)
}
Expand Down
28 changes: 14 additions & 14 deletions R/PrestoGP_Vecchia.R
Original file line number Diff line number Diff line change
Expand Up @@ -71,24 +71,24 @@ setMethod("prestogp_predict", "VecchiaModel", function(model, X, locs, m = NULL,

locs.train.scaled <- scale_locs(model, model@locs_train)[[1]]
locs.scaled <- scale_locs(model, list(locs))[[1]]
vec.approx.test <- vecchia_specify(locs.train.scaled, m, locs.pred = locs.scaled, ordering.pred=ordering.pred, pred.cond=pred.cond)
vec.approx.test <- vecchia_specify(locs.train.scaled, m, locs.pred = locs.scaled, ordering.pred = ordering.pred, pred.cond = pred.cond)

## carry out prediction
if (!model@apanasovich) {
pred <- vecchia_prediction(res, vec.approx.test, c(model@covparams[1], 1, model@covparams[3]), model@covparams[4], return.values=return.values)
} else {
pred <- vecchia_prediction(res, vec.approx.test, c(model@covparams[1], model@covparams[2], model@covparams[3]), model@covparams[4], return.values=return.values)
}
## carry out prediction
if (!model@apanasovich) {
pred <- vecchia_prediction(res, vec.approx.test, c(model@covparams[1], 1, model@covparams[3]), model@covparams[4], return.values = return.values)
} else {
pred <- vecchia_prediction(res, vec.approx.test, c(model@covparams[1], model@covparams[2], model@covparams[3]), model@covparams[4], return.values = return.values)
}

# 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)
return.list <- list(means = Vec.mean)
} else {
warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.")
Vec.sds <- sqrt(pred$var.pred + model@covparams[4])
return.list <- list(means = Vec.mean, sds = vec.sds)
warning("Variance estimates do not include model fitting variance and are anticonservative. Use with caution.")
Vec.sds <- sqrt(pred$var.pred + model@covparams[4])
return.list <- list(means = Vec.mean, sds = vec.sds)
}

return(return.list)
Expand Down Expand Up @@ -130,13 +130,13 @@ setMethod("check_input_pred", "VecchiaModel", function(model, X, locs) {
stop("X must be a matrix")
}
if (ncol(locs) != ncol(model@locs_train[[1]])) {
stop("locs must have the same number of columns as locs_train")
stop("locs must have the same number of columns as locs_train")
}
if (nrow(X) != nrow(locs)) {
stop("X must have the same number of rows as locs")
stop("X must have the same number of rows as locs")
}
if (ncol(X) != ncol(model@X_train)) {
stop("X and X_train must have the same number of predictors")
stop("X and X_train must have the same number of predictors")
}
invisible(model)
})
Expand Down
5 changes: 2 additions & 3 deletions R/RcppExports.R
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,9 @@
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

na_omit_c <- function(x) {
.Call('_PrestoGP_na_omit_c', PACKAGE = 'PrestoGP', x)
.Call("_PrestoGP_na_omit_c", PACKAGE = "PrestoGP", x)
}

createU_helper_mat <- function(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning) {
.Call('_PrestoGP_createU_helper_mat', PACKAGE = 'PrestoGP', olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning)
.Call("_PrestoGP_createU_helper_mat", PACKAGE = "PrestoGP", olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning)
}

1 change: 0 additions & 1 deletion inst/doc/test.R
Original file line number Diff line number Diff line change
@@ -1,3 +1,2 @@
## -----------------------------------------------------------------------------
cat("test", "\n")

83 changes: 49 additions & 34 deletions src/RcppExports.cpp
Original file line number Diff line number Diff line change
@@ -1,55 +1,70 @@
// Generated by using Rcpp::compileAttributes() -> do not edit by hand
// Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#include <RcppArmadillo.h>
#include <Rcpp.h>
#include <RcppArmadillo.h>

using namespace Rcpp;

#ifdef RCPP_USE_GLOBAL_ROSTREAM
Rcpp::Rostream<true>& Rcpp::Rcout = Rcpp::Rcpp_cout_get();
Rcpp::Rostream<false>& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
Rcpp::Rostream<true> &Rcpp::Rcout = Rcpp::Rcpp_cout_get();
Rcpp::Rostream<false> &Rcpp::Rcerr = Rcpp::Rcpp_cerr_get();
#endif

// na_omit_c
arma::vec na_omit_c(arma::vec x);
RcppExport SEXP _PrestoGP_na_omit_c(SEXP xSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< arma::vec >::type x(xSEXP);
rcpp_result_gen = Rcpp::wrap(na_omit_c(x));
return rcpp_result_gen;
END_RCPP
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter<arma::vec>::type x(xSEXP);
rcpp_result_gen = Rcpp::wrap(na_omit_c(x));
return rcpp_result_gen;
END_RCPP
}
// createU_helper_mat
arma::sp_mat createU_helper_mat(const arma::mat& olocs, const arma::vec& ondx, const arma::mat& curqys, const arma::mat& curqzs, const arma::mat& vijs, const arma::mat& aijs, const arma::mat& full_const, const arma::vec& nugget, const arma::vec& sig2, const arma::vec& U_beginning);
RcppExport SEXP _PrestoGP_createU_helper_mat(SEXP olocsSEXP, SEXP ondxSEXP, SEXP curqysSEXP, SEXP curqzsSEXP, SEXP vijsSEXP, SEXP aijsSEXP, SEXP full_constSEXP, SEXP nuggetSEXP, SEXP sig2SEXP, SEXP U_beginningSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter< const arma::mat& >::type olocs(olocsSEXP);
Rcpp::traits::input_parameter< const arma::vec& >::type ondx(ondxSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type curqys(curqysSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type curqzs(curqzsSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type vijs(vijsSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type aijs(aijsSEXP);
Rcpp::traits::input_parameter< const arma::mat& >::type full_const(full_constSEXP);
Rcpp::traits::input_parameter< const arma::vec& >::type nugget(nuggetSEXP);
Rcpp::traits::input_parameter< const arma::vec& >::type sig2(sig2SEXP);
Rcpp::traits::input_parameter< const arma::vec& >::type U_beginning(U_beginningSEXP);
rcpp_result_gen = Rcpp::wrap(createU_helper_mat(olocs, ondx, curqys, curqzs, vijs, aijs, full_const, nugget, sig2, U_beginning));
return rcpp_result_gen;
END_RCPP
arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx,
const arma::mat &curqys,
const arma::mat &curqzs, const arma::mat &vijs,
const arma::mat &aijs,
const arma::mat &full_const,
const arma::vec &nugget, const arma::vec &sig2,
const arma::vec &U_beginning);
RcppExport SEXP _PrestoGP_createU_helper_mat(SEXP olocsSEXP, SEXP ondxSEXP,
SEXP curqysSEXP, SEXP curqzsSEXP,
SEXP vijsSEXP, SEXP aijsSEXP,
SEXP full_constSEXP,
SEXP nuggetSEXP, SEXP sig2SEXP,
SEXP U_beginningSEXP) {
BEGIN_RCPP
Rcpp::RObject rcpp_result_gen;
Rcpp::RNGScope rcpp_rngScope_gen;
Rcpp::traits::input_parameter<const arma::mat &>::type olocs(olocsSEXP);
Rcpp::traits::input_parameter<const arma::vec &>::type ondx(ondxSEXP);
Rcpp::traits::input_parameter<const arma::mat &>::type curqys(curqysSEXP);
Rcpp::traits::input_parameter<const arma::mat &>::type curqzs(curqzsSEXP);
Rcpp::traits::input_parameter<const arma::mat &>::type vijs(vijsSEXP);
Rcpp::traits::input_parameter<const arma::mat &>::type aijs(aijsSEXP);
Rcpp::traits::input_parameter<const arma::mat &>::type full_const(
full_constSEXP);
Rcpp::traits::input_parameter<const arma::vec &>::type nugget(nuggetSEXP);
Rcpp::traits::input_parameter<const arma::vec &>::type sig2(sig2SEXP);
Rcpp::traits::input_parameter<const arma::vec &>::type U_beginning(
U_beginningSEXP);
rcpp_result_gen =
Rcpp::wrap(createU_helper_mat(olocs, ondx, curqys, curqzs, vijs, aijs,
full_const, nugget, sig2, U_beginning));
return rcpp_result_gen;
END_RCPP
}

static const R_CallMethodDef CallEntries[] = {
{"_PrestoGP_na_omit_c", (DL_FUNC) &_PrestoGP_na_omit_c, 1},
{"_PrestoGP_createU_helper_mat", (DL_FUNC) &_PrestoGP_createU_helper_mat, 10},
{NULL, NULL, 0}
};
{"_PrestoGP_na_omit_c", (DL_FUNC)&_PrestoGP_na_omit_c, 1},
{"_PrestoGP_createU_helper_mat", (DL_FUNC)&_PrestoGP_createU_helper_mat,
10},
{NULL, NULL, 0}};

RcppExport void R_init_PrestoGP(DllInfo *dll) {
R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
R_useDynamicSymbols(dll, FALSE);
R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
R_useDynamicSymbols(dll, FALSE);
}
Loading

0 comments on commit 0460202

Please sign in to comment.