From 63dbb7840329c158463c0e0aa163ffc617b0b206 Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Wed, 3 Jan 2024 13:53:40 -0500 Subject: [PATCH 1/9] enable linting on build-workflow --- .github/workflows/lint.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml index df2a9fb..cf3d710 100644 --- a/.github/workflows/lint.yaml +++ b/.github/workflows/lint.yaml @@ -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] From ea0c663b2d2fbd06413709fae21fe090be821d01 Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Wed, 3 Jan 2024 14:30:51 -0500 Subject: [PATCH 2/9] Run auto-lint on vscode --- .lintr | 5 + tests/testthat/sim_multivariate.R | 129 ++++----- tests/testthat/sim_vecchia.R | 24 +- tests/testthat/sim_vecchia_st.R | 30 ++- tests/testthat/test-Log_Likelihood.R | 34 ++- .../test-PrestoGP_CreateU_Multivariate.R | 54 ++-- ...test-PrestoGP_Multivariate_Vecchia_Model.R | 254 +++++++++++------- .../testthat/test-PrestoGP_Univariate_Model.R | 236 +++++++++------- 8 files changed, 446 insertions(+), 320 deletions(-) create mode 100644 .lintr diff --git a/.lintr b/.lintr new file mode 100644 index 0000000..16c0bd2 --- /dev/null +++ b/.lintr @@ -0,0 +1,5 @@ +linters: linters_with_defaults( + line_length_linter = NULL, + commented_code_linter = NULL +) +encoding: "UTF-8" diff --git a/tests/testthat/sim_multivariate.R b/tests/testthat/sim_multivariate.R index f33725b..f1bf46f 100644 --- a/tests/testthat/sim_multivariate.R +++ b/tests/testthat/sim_multivariate.R @@ -13,35 +13,35 @@ set.seed(1234) # To keep the simulation small and tractable to quickly run, we # simulate 25 spatial locations and 5 temporal locations -n.spatial=25 -n.spatial.xy = sqrt(n.spatial) -n.temporal = 5 +n.spatial <- 25 +n.spatial.xy <- sqrt(n.spatial) +n.temporal <- 5 # Mean trend, X(s): spatiotemporal data, but the trend is spatial only. # Each outcome has a different combination of true betas # p = number of potential covariates -p=10 +p <- 10 -p.nz=4 -beta1= c(rep(0,p-p.nz),rep(1,p.nz)) -p.nz=2 -beta2= c(rep(1,p.nz),rep(0,p-p.nz)) -p.nz=3 -beta3= c(rep(1,p.nz),rep(0,p-p.nz)) +p.nz <- 4 +beta1 <- c(rep(0, p - p.nz), rep(1, p.nz)) +p.nz <- 2 +beta2 <- c(rep(1, p.nz), rep(0, p - p.nz)) +p.nz <- 3 +beta3 <- c(rep(1, p.nz), rep(0, p - p.nz)) # Combine all of the betas into 1 vector -beta.all <- c(beta1,beta2,beta3) +beta.all <- c(beta1, beta2, beta3) # correlated predictors -Sigma.X=exp(-rdist(sample(1:p))/3) -X=mvrnorm(n.spatial,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial, rep(0, p), Sigma.X) # replicate for space-time for a single-outcomes -X.st <- rbind(X,X,X,X,X) +X.st <- rbind(X, X, X, X, X) # Create the multivariate X-matrix, for multiple outcomes # The rows are by outcome (i), space (j), and then time (k) -X.all <- psych::superMatrix(list(X.st,X.st,X.st)) +X.all <- psych::superMatrix(list(X.st, X.st, X.st)) # S-T design matrix with X's repeated is needed @@ -62,12 +62,12 @@ mean.trend.st <- X.all %*% beta.all # rho is the correlation between outcome S-T errors, 1 for with itself rho.vec <- c(0.8, 0.3, 0.5) -rho <- matrix(0,nrow=3,ncol=3) +rho <- matrix(0, nrow = 3, ncol = 3) rho[upper.tri(rho)] <- rho.vec -rho <- rho+t(rho)+diag(1,3) +rho <- rho + t(rho) + diag(1, 3) # nu, marginal smoothness: -marg.smoothness <- c(0.5,1,1.5) +marg.smoothness <- c(0.5, 1, 1.5) # Non-spatial error/noise/nugget # epsilon(s,t) will be simulated based on the nugget parameters, @@ -75,7 +75,7 @@ nuggets <- c(0.5, 0.7, 2) # marginal variances of the Matern # The marginal variances are scaled linearly with the nuggets -x.variance <- runif(3,1.5,4) +x.variance <- runif(3, 1.5, 4) marg.var <- nuggets + x.variance # ranges, the true ranges, we will scale coordinates by this, but then use a @@ -84,16 +84,16 @@ marg.var <- nuggets + x.variance # outcome1 space and time, outcome 2 space and time, outcome 3 space and time # Spatial domain is on the unit [0,1] # Temporal domain is a year [0,365] -ranges <- c(0.8,0.1,0.5) +ranges <- c(0.8, 0.1, 0.5) # set up the true spatial and temporal dimensions -space.dim1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -space.dim2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -space.dim3 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -time.dim1 <- seq(0,3,length.out = n.temporal)+rnorm(n.temporal,0,0.0001) -time.dim2 <- seq(0,3,length.out = n.temporal)+rnorm(n.temporal,0,0.0001) -time.dim3 <- seq(0,3,length.out = n.temporal)+rnorm(n.temporal,0,0.0001) +space.dim1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +space.dim2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +space.dim3 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +time.dim1 <- seq(0, 3, length.out = n.temporal) + rnorm(n.temporal, 0, 0.0001) +time.dim2 <- seq(0, 3, length.out = n.temporal) + rnorm(n.temporal, 0, 0.0001) +time.dim3 <- seq(0, 3, length.out = n.temporal) + rnorm(n.temporal, 0, 0.0001) # scale the space and time dimensions according to the true covariance ranges @@ -101,62 +101,69 @@ d11 <- fields::rdist(expand.grid(space.dim1, space.dim1, time.dim1)) d22 <- fields::rdist(expand.grid(space.dim2, space.dim2, time.dim2)) d33 <- fields::rdist(expand.grid(space.dim3, space.dim3, time.dim3)) -d12 <- fields::rdist(expand.grid(space.dim1, space.dim1, time.dim1), - expand.grid(space.dim2, space.dim2, time.dim2)) +d12 <- fields::rdist( + expand.grid(space.dim1, space.dim1, time.dim1), + expand.grid(space.dim2, space.dim2, time.dim2) +) d21 <- t(d12) -d13 <- fields::rdist(expand.grid(space.dim1, space.dim1, time.dim1), - expand.grid(space.dim3, space.dim3, time.dim3)) +d13 <- fields::rdist( + expand.grid(space.dim1, space.dim1, time.dim1), + expand.grid(space.dim3, space.dim3, time.dim3) +) d31 <- t(d13) -d23 <- fields::rdist(expand.grid(space.dim2, space.dim2, time.dim2), - expand.grid(space.dim3, space.dim3, time.dim3)) +d23 <- fields::rdist( + expand.grid(space.dim2, space.dim2, time.dim2), + expand.grid(space.dim3, space.dim3, time.dim3) +) d32 <- t(d23) ## Create the correlation matrices -Sigma11 <- marg.var[1] * fields::Matern(d11,range = ranges[1], smoothness = marg.smoothness[1]) -Sigma22 <- marg.var[2] * fields::Matern(d22,range = ranges[2], smoothness = marg.smoothness[2]) -Sigma33 <- marg.var[3] * fields::Matern(d33,range = ranges[3], smoothness = marg.smoothness[3]) +Sigma11 <- marg.var[1] * fields::Matern(d11, range = ranges[1], smoothness = marg.smoothness[1]) +Sigma22 <- marg.var[2] * fields::Matern(d22, range = ranges[2], smoothness = marg.smoothness[2]) +Sigma33 <- marg.var[3] * fields::Matern(d33, range = ranges[3], smoothness = marg.smoothness[3]) vii <- marg.smoothness[1] vjj <- marg.smoothness[2] -vij <- (vii+vjj)/2 -aii <- 1/ranges[1] -ajj <- 1/ranges[2] -aij <- sqrt((aii^2+ajj^2)/2) -Sigma12 <- rho[1,2] * sqrt(marg.var[1]) * sqrt(marg.var[2]) * aii^vii * - ajj^vjj * gamma(vij) / (aij^(2*vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(d12, smoothness=vij, alpha=aij) +vij <- (vii + vjj) / 2 +aii <- 1 / ranges[1] +ajj <- 1 / ranges[2] +aij <- sqrt((aii^2 + ajj^2) / 2) +Sigma12 <- rho[1, 2] * sqrt(marg.var[1]) * sqrt(marg.var[2]) * aii^vii * + ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(d12, smoothness = vij, alpha = aij) Sigma21 <- t(Sigma12) vii <- marg.smoothness[1] vjj <- marg.smoothness[3] -vij <- (vii+vjj)/2 -aii <- 1/ranges[1] -ajj <- 1/ranges[3] -aij <- sqrt((aii^2+ajj^2)/2) -Sigma13 <- rho[1,3] * sqrt(marg.var[1]) * sqrt(marg.var[3]) * aii^vii * - ajj^vjj * gamma(vij) / (aij^(2*vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(d13, smoothness=vij, alpha=aij) +vij <- (vii + vjj) / 2 +aii <- 1 / ranges[1] +ajj <- 1 / ranges[3] +aij <- sqrt((aii^2 + ajj^2) / 2) +Sigma13 <- rho[1, 3] * sqrt(marg.var[1]) * sqrt(marg.var[3]) * aii^vii * + ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(d13, smoothness = vij, alpha = aij) Sigma31 <- t(Sigma13) vii <- marg.smoothness[2] vjj <- marg.smoothness[3] -vij <- (vii+vjj)/2 -aii <- 1/ranges[2] -ajj <- 1/ranges[3] -aij <- sqrt((aii^2+ajj^2)/2) -Sigma23 <- rho[2,3] * sqrt(marg.var[2]) * sqrt(marg.var[3]) * aii^vii * - ajj^vjj * gamma(vij) / (aij^(2*vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(d23, smoothness=vij, alpha=aij) +vij <- (vii + vjj) / 2 +aii <- 1 / ranges[2] +ajj <- 1 / ranges[3] +aij <- sqrt((aii^2 + ajj^2) / 2) +Sigma23 <- rho[2, 3] * sqrt(marg.var[2]) * sqrt(marg.var[3]) * aii^vii * + ajj^vjj * gamma(vij) / (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(d23, smoothness = vij, alpha = aij) Sigma32 <- t(Sigma23) -#Combine into the super Multivariate covariance matrix -Sigma.All <- rbind( cbind(Sigma11,Sigma12,Sigma13), - cbind(Sigma21,Sigma22,Sigma23), - cbind(Sigma31,Sigma32,Sigma33) - ) +# Combine into the super Multivariate covariance matrix +Sigma.All <- rbind( + cbind(Sigma11, Sigma12, Sigma13), + cbind(Sigma21, Sigma22, Sigma23), + cbind(Sigma31, Sigma32, Sigma33) +) # Cholesky decomposition L.C <- chol(Sigma.All) diff --git a/tests/testthat/sim_vecchia.R b/tests/testthat/sim_vecchia.R index 992e3f0..63c77cf 100755 --- a/tests/testthat/sim_vecchia.R +++ b/tests/testthat/sim_vecchia.R @@ -7,10 +7,10 @@ n.spatial.xy <- 50 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,18 +22,22 @@ ranges <- runif(1, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) -Sigma.All <- marg.var*Matern(rdist(locs), range=ranges, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locs), + range = ranges, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^2) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^2) +nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance) + x.variance +) diff --git a/tests/testthat/sim_vecchia_st.R b/tests/testthat/sim_vecchia_st.R index 9438ee3..da65aa4 100755 --- a/tests/testthat/sim_vecchia_st.R +++ b/tests/testthat/sim_vecchia_st.R @@ -7,10 +7,10 @@ n.spatial.xy <- 15 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^3,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^3, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,23 +22,27 @@ ranges <- runif(2, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc3 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc3 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2, loc3)) locss <- locs -locss[,1:2] <- locss[,1:2] / ranges[1] -locss[,3] <- locs[,3] / ranges[2] +locss[, 1:2] <- locss[, 1:2] / ranges[1] +locss[, 3] <- locs[, 3] / ranges[2] -Sigma.All <- marg.var*Matern(rdist(locss), range=1, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locss), + range = 1, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^3) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^3) +nug.error <- nuggets * rnorm(n.spatial.xy^3) y <- mean.trend.st + st.error + nug.error -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, loc3, L.C, +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, loc3, L.C, st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance, locss) + x.variance, locss +) diff --git a/tests/testthat/test-Log_Likelihood.R b/tests/testthat/test-Log_Likelihood.R index 2933641..db5ec55 100755 --- a/tests/testthat/test-Log_Likelihood.R +++ b/tests/testthat/test-Log_Likelihood.R @@ -55,20 +55,26 @@ test_that("negloglik.full", { d <- rdist(locs) pseq <- create.param.sequence(1) - res.optim.NM <- optim(par=params.init, fn=negloglik.full, d=d, y=y, - param.seq=pseq, control=list(maxit=5000)) + res.optim.NM <- optim( + par = params.init, fn = negloglik.full, d = d, y = y, + param.seq = pseq, control = list(maxit = 5000) + ) LL.full <- negloglik.full(res.optim.NM$par, d, y, pseq) - params.final <- c(exp(res.optim.NM$par[1:2]), - gtools::inv.logit(res.optim.NM$par[3], 0, 2.5), - exp(res.optim.NM$par[4])) + params.final <- c( + exp(res.optim.NM$par[1:2]), + gtools::inv.logit(res.optim.NM$par[3], 0, 2.5), + exp(res.optim.NM$par[4]) + ) - pgp.params <- create.initial.values.flex(params.final[1], - params.final[2], - params.final[3], - params.final[4], - 1, 1) + pgp.params <- create.initial.values.flex( + params.final[1], + params.final[2], + params.final[3], + params.final[4], + 1, 1 + ) LL.full.pgp <- mvnegloglik.full(pgp.params, list(locs), y, pseq) @@ -110,9 +116,11 @@ test_that("mvnegloglik", { ) pseq <- create.param.sequence(P) vec.approx <- vecchia_Mspecify(locs.list, 25) - neg_likelihood <- mvnegloglik(logparams, vec.approx, - unlist(y.list), pseq, P) - expect_equal(34384.23, neg_likelihood, tolerance=1e-2) + neg_likelihood <- mvnegloglik( + logparams, vec.approx, + unlist(y.list), pseq, P + ) + expect_equal(34384.23, neg_likelihood, tolerance = 1e-2) }) test_that("mvnegloglik_ST", { diff --git a/tests/testthat/test-PrestoGP_CreateU_Multivariate.R b/tests/testthat/test-PrestoGP_CreateU_Multivariate.R index cd0fb65..c50fc24 100755 --- a/tests/testthat/test-PrestoGP_CreateU_Multivariate.R +++ b/tests/testthat/test-PrestoGP_CreateU_Multivariate.R @@ -95,29 +95,29 @@ test_that("createUMultivariate", { locs <- as.matrix(expand.grid(dim1, dim2)) - params <- c(3,1.5,0.6,2) + params <- c(3, 1.5, 0.6, 2) - vec.approx <- vecchia_specify(locs, m=5) - U.obj <- createU(vec.approx, covparms=params[1:3], nuggets=params[4]) + vec.approx <- vecchia_specify(locs, m = 5) + U.obj <- createU(vec.approx, covparms = params[1:3], nuggets = params[4]) vec.mapprox <- vecchia_Mspecify(list(locs), 5) - U.mobj <- createUMultivariate(vec.mapprox, c(params,1)) -# Should produce the same ordered locs as GPvecchia in the univariate case + U.mobj <- createUMultivariate(vec.mapprox, c(params, 1)) + # Should produce the same ordered locs as GPvecchia in the univariate case expect_equal(vec.approx$locsord, vec.mapprox$locsord) expect_equal(vec.approx$ord, vec.mapprox$ord) -# Should produce the same U matrix as GPvecchia in the univariate case - expect_equal(sum(abs(U.obj$U-U.mobj$U)), 0, tolerance=1e-3) + # Should produce the same U matrix as GPvecchia in the univariate case + expect_equal(sum(abs(U.obj$U - U.mobj$U)), 0, tolerance = 1e-3) -# Should identify the same nearest neighbors and latents as GPvecchia in the -# univariate case - NNarray <- vec.approx$U.prep$revNNarray[,-6] - cond <- vec.approx$U.prep$revCond[,-6] + # Should identify the same nearest neighbors and latents as GPvecchia in the + # univariate case + NNarray <- vec.approx$U.prep$revNNarray[, -6] + cond <- vec.approx$U.prep$revCond[, -6] for (i in 6:100) { - cur.y <- NNarray[i,cond[i,]] + cur.y <- NNarray[i, cond[i, ]] expect_equal(sort(cur.y), sort(vec.mapprox$q.list$q.y[[i]])) - if (sum(!cond[i,])>0) { - cur.z <- NNarray[i,!cond[i,]] + if (sum(!cond[i, ]) > 0) { + cur.z <- NNarray[i, !cond[i, ]] expect_equal(sort(cur.z), sort(vec.mapprox$q.list$q.z[[i]])) } } @@ -125,20 +125,24 @@ test_that("createUMultivariate", { source("sim_multivariate.R") vec.mapprox <- vecchia_Mspecify(locs.list, 25) - U.mobj <- createUMultivariate(vec.mapprox, c(marg.var, ranges, - marg.smoothness, - nuggets, rho.vec)) + U.mobj <- createUMultivariate(vec.mapprox, c( + marg.var, ranges, + marg.smoothness, + nuggets, rho.vec + )) Umat.tst <- readRDS("Umat.rds") - expect_equal(sum(abs(U.mobj$U-Umat.tst)), 0, tolerance=1e-4) + expect_equal(sum(abs(U.mobj$U - Umat.tst)), 0, tolerance = 1e-4) vec.mapprox.full <- vecchia_Mspecify(locs.list, 374) - U.mobj.full <- createUMultivariate(vec.mapprox.full, c(marg.var, ranges, - marg.smoothness, - nuggets, rho.vec)) + U.mobj.full <- createUMultivariate(vec.mapprox.full, c( + marg.var, ranges, + marg.smoothness, + nuggets, rho.vec + )) Sigma.hat <- solve(U.mobj.full$U %*% t(U.mobj.full$U)) - Sigma.hat <- Sigma.hat[U.mobj.full$latent,U.mobj.full$latent] - Sigma.hat <- Sigma.hat[order(U.mobj.full$ord),order(U.mobj.full$ord)] -# Vecchia approximation of covariance should equal true covariance when m=n-1 - expect_equal(sum(abs(Sigma.All-Sigma.hat)), 0, tolerance=1e-4) + Sigma.hat <- Sigma.hat[U.mobj.full$latent, U.mobj.full$latent] + Sigma.hat <- Sigma.hat[order(U.mobj.full$ord), order(U.mobj.full$ord)] + # Vecchia approximation of covariance should equal true covariance when m=n-1 + expect_equal(sum(abs(Sigma.All - Sigma.hat)), 0, tolerance = 1e-4) }) diff --git a/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R b/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R index 5f1a64e..c4cd698 100755 --- a/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R +++ b/tests/testthat/test-PrestoGP_Multivariate_Vecchia_Model.R @@ -1,163 +1,215 @@ test_that("Invalid locs input", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(as.matrix(1:3)), list(as.matrix(1:3)), - as.matrix(1:3)), - "locs must be a list for multivariate models") + expect_error( + prestogp_fit( + model, list(as.matrix(1:3)), list(as.matrix(1:3)), + as.matrix(1:3) + ), + "locs must be a list for multivariate models" + ) }) test_that("Invalid Y input", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:3), list(as.matrix(1:3)), - list(as.matrix(1:3))), - "Y must be a list for multivariate models") + expect_error( + prestogp_fit( + model, as.matrix(1:3), list(as.matrix(1:3)), + list(as.matrix(1:3)) + ), + "Y must be a list for multivariate models" + ) }) test_that("Invalid X input", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(as.matrix(1:3)), as.matrix(1:3), - list(as.matrix(1:3))), - "X must be a list for multivariate models") + expect_error( + prestogp_fit( + model, list(as.matrix(1:3)), as.matrix(1:3), + list(as.matrix(1:3)) + ), + "X must be a list for multivariate models" + ) }) test_that("locs/Y length mismatch", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:3, 2:4), list(as.matrix(1:3)), - list(as.matrix(1:3))), - "locs and Y must have the same length") + expect_error( + prestogp_fit( + model, list(1:3, 2:4), list(as.matrix(1:3)), + list(as.matrix(1:3)) + ), + "locs and Y must have the same length" + ) }) test_that("locs/X length mismatch", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:3), list(as.matrix(1:3), - as.matrix(2:4)), - list(as.matrix(1:3))), - "locs and X must have the same length") + expect_error( + prestogp_fit( + model, list(1:3), list( + as.matrix(1:3), + as.matrix(2:4) + ), + list(as.matrix(1:3)) + ), + "locs and X must have the same length" + ) }) test_that("locs not a matrix", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:3), list(as.matrix(1:3)), - list(1:3)), - "Each locs must be a matrix") + expect_error( + prestogp_fit( + model, list(1:3), list(as.matrix(1:3)), + list(1:3) + ), + "Each locs must be a matrix" + ) }) test_that("X not a matrix", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:3), list(1:3), - list(as.matrix(1:3))), - "Each X must be a matrix") + expect_error( + prestogp_fit( + model, list(1:3), list(1:3), + list(as.matrix(1:3)) + ), + "Each X must be a matrix" + ) }) test_that("Y not a numeric matrix/vector", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list("foo"), list(as.matrix(1:3)), - list(as.matrix(1:3))), - "Each Y must be a numeric vector or matrix") + expect_error( + prestogp_fit( + model, list("foo"), list(as.matrix(1:3)), + list(as.matrix(1:3)) + ), + "Each Y must be a numeric vector or matrix" + ) }) test_that("locs with differing numbers of columns", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(1:4, 1:4), - list(as.matrix(1:4), as.matrix(1:4)), - list(as.matrix(1:4), matrix(1:4, nrow=2))), - "All locs must have the same number of columns") + expect_error( + prestogp_fit( + model, list(1:4, 1:4), + list(as.matrix(1:4), as.matrix(1:4)), + list(as.matrix(1:4), matrix(1:4, nrow = 2)) + ), + "All locs must have the same number of columns" + ) }) test_that("Y has multiple columns", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(matrix(1:4, nrow=2)), - list(as.matrix(1:4)), list(as.matrix(1:4))), - "Each Y must have only 1 column") + expect_error( + prestogp_fit( + model, list(matrix(1:4, nrow = 2)), + list(as.matrix(1:4)), list(as.matrix(1:4)) + ), + "Each Y must have only 1 column" + ) }) test_that("nrow(Y) != nrow(locs)", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(as.matrix(1:4)), - list(as.matrix(1:4)), list(as.matrix(1:3))), - "Each Y must have the same number of rows as locs") + expect_error( + prestogp_fit( + model, list(as.matrix(1:4)), + list(as.matrix(1:4)), list(as.matrix(1:3)) + ), + "Each Y must have the same number of rows as locs" + ) }) test_that("nrow(Y) != nrow(X)", { model <- new("MultivariateVecchiaModel") - expect_error(prestogp_fit(model, list(as.matrix(1:4)), - list(as.matrix(1:3)), list(as.matrix(1:4))), - "Each Y must have the same number of rows as X") + expect_error( + prestogp_fit( + model, list(as.matrix(1:4)), + list(as.matrix(1:3)), list(as.matrix(1:4)) + ), + "Each Y must have the same number of rows as X" + ) }) test_that("Simulated dataset multivariate spatial", { - load("sim_multivariate_big.RData") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors=25) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list, X.st, locs.list, - scaling = c(1, 1), apanasovich = TRUE, verbose = FALSE, - optim.control = list( - trace = 0, maxit = 5000, - reltol = 1e-3 - ) + load("sim_multivariate_big.RData") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list, X.st, locs.list, + scaling = c(1, 1), apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 ) - beta.out <- as.vector(pgp.mmodel1@beta) - params.out <- pgp.mmodel1@covparams - - expect_length(beta.out, 31) - expect_length(params.out, 15) - expect_equal(beta.out, c(0, 0.96, 0.93, 0.92, 0.89, rep(0, 2), 0.03, - rep(0,3), 0.57, 0.72, 1.11, 1, 0, 0.06, rep(0,4), - 0.93, 0.87, 1.03, 0.92, 0.05, rep(0, 2), 0.01, - 0.05, 0), tolerance=0.05) - expect_equal(params.out[1], 1.7, tolerance=2.5) - expect_equal(params.out[2], 2.9, tolerance=3.3) - expect_equal(params.out[3], 3.2, tolerance=5) - expect_equal(params.out[4], 0.71, tolerance=10) - expect_equal(params.out[5], 0.56, tolerance=2.8) - expect_equal(params.out[6], 0.4, tolerance=1.7) - expect_equal(params.out[7], 0.63, tolerance=1.3) - expect_equal(params.out[8], 0.47, tolerance=1.2) - expect_equal(params.out[9], 0.74, tolerance=1.1) - expect_equal(params.out[10], 1.8, tolerance=1.4) - expect_equal(params.out[11], 2.6, tolerance=2.4) - expect_equal(params.out[12], 1.4, tolerance=0.8) - expect_equal(params.out[13], 0.15, tolerance=0.6) - expect_equal(params.out[14], 0.32, tolerance=0.4) - expect_equal(params.out[15], 0.17, tolerance=0.3) + ) + beta.out <- as.vector(pgp.mmodel1@beta) + params.out <- pgp.mmodel1@covparams + + expect_length(beta.out, 31) + expect_length(params.out, 15) + expect_equal(beta.out, c( + 0, 0.96, 0.93, 0.92, 0.89, rep(0, 2), 0.03, + rep(0, 3), 0.57, 0.72, 1.11, 1, 0, 0.06, rep(0, 4), + 0.93, 0.87, 1.03, 0.92, 0.05, rep(0, 2), 0.01, + 0.05, 0 + ), tolerance = 0.05) + expect_equal(params.out[1], 1.7, tolerance = 2.5) + expect_equal(params.out[2], 2.9, tolerance = 3.3) + expect_equal(params.out[3], 3.2, tolerance = 5) + expect_equal(params.out[4], 0.71, tolerance = 10) + expect_equal(params.out[5], 0.56, tolerance = 2.8) + expect_equal(params.out[6], 0.4, tolerance = 1.7) + expect_equal(params.out[7], 0.63, tolerance = 1.3) + expect_equal(params.out[8], 0.47, tolerance = 1.2) + expect_equal(params.out[9], 0.74, tolerance = 1.1) + expect_equal(params.out[10], 1.8, tolerance = 1.4) + expect_equal(params.out[11], 2.6, tolerance = 2.4) + expect_equal(params.out[12], 1.4, tolerance = 0.8) + expect_equal(params.out[13], 0.15, tolerance = 0.6) + expect_equal(params.out[14], 0.32, tolerance = 0.4) + expect_equal(params.out[15], 0.17, tolerance = 0.3) }) test_that("Simulated dataset multivariate spatiotemporal", { - source("sim_multivariate_big_st.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list, X.st, locs.list, - scaling = c(1, 1, 2), verbose = FALSE, - optim.control = list( - trace = 0, maxit = 5000, - reltol = 1e-3 - ) + source("sim_multivariate_big_st.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list, X.st, locs.list, + scaling = c(1, 1, 2), verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 ) - beta.out <- as.vector(pgp.mmodel1@beta) - params.out <- pgp.mmodel1@covparams - - expect_length(beta.out, 31) - expect_length(params.out, 18) - expect_equal(beta.out, c( - 0, 0.91, 0.86, 0.82, 0.97, rep(0, 6), 0.95, - 0.97, 0.92, 0.78, rep(0, 6), 0.8, 0.97, - 1.04, 0.81, rep(0, 6) - ), tolerance = 1.1) + ) + beta.out <- as.vector(pgp.mmodel1@beta) + params.out <- pgp.mmodel1@covparams + + expect_length(beta.out, 31) + expect_length(params.out, 18) + expect_equal(beta.out, c( + 0, 0.91, 0.86, 0.82, 0.97, rep(0, 6), 0.95, + 0.97, 0.92, 0.78, rep(0, 6), 0.8, 0.97, + 1.04, 0.81, rep(0, 6) + ), tolerance = 1.1) }) test_that("Simulated dataset multivariate spatial prediction", { - source("sim_multivariate_big_pred.R") - pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) - pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, - locs.list.otr, - scaling = c(1, 1), - apanasovich = TRUE, verbose = FALSE, - optim.control = list( - trace = 0, maxit = 5000, - reltol = 1e-3 - ) + source("sim_multivariate_big_pred.R") + pgp.mmodel1 <- new("MultivariateVecchiaModel", n_neighbors = 25) + pgp.mmodel1 <- prestogp_fit(pgp.mmodel1, y.list.otr, X.st.otr, + locs.list.otr, + scaling = c(1, 1), + apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 ) + ) - pgp.mmodel1.pred <- prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst) + pgp.mmodel1.pred <- prestogp_predict(pgp.mmodel1, X.st.otst, locs.list.otst) - mse <- mean((pgp.mmodel1.pred$means - unlist(y.list.otst))^2) + mse <- mean((pgp.mmodel1.pred$means - unlist(y.list.otst))^2) - expect_equal(mse, 1.99, tolerance = 0.1) + expect_equal(mse, 1.99, tolerance = 0.1) }) diff --git a/tests/testthat/test-PrestoGP_Univariate_Model.R b/tests/testthat/test-PrestoGP_Univariate_Model.R index 53a56b1..0eca2c5 100755 --- a/tests/testthat/test-PrestoGP_Univariate_Model.R +++ b/tests/testthat/test-PrestoGP_Univariate_Model.R @@ -1,134 +1,176 @@ test_that("Invalid locs input", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:3), as.matrix(1:3), - 1:3), - "locs must be a matrix") + expect_error( + prestogp_fit( + model, as.matrix(1:3), as.matrix(1:3), + 1:3 + ), + "locs must be a matrix" + ) }) test_that("Invalid Y input", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, "y", as.matrix(1:3), - as.matrix(1:3)), - "Y must be a numeric vector or matrix") + expect_error( + prestogp_fit( + model, "y", as.matrix(1:3), + as.matrix(1:3) + ), + "Y must be a numeric vector or matrix" + ) }) test_that("Invalid X input", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:3), 1:3, - as.matrix(1:3)), - "X must be a matrix") + expect_error( + prestogp_fit( + model, as.matrix(1:3), 1:3, + as.matrix(1:3) + ), + "X must be a matrix" + ) }) test_that("Y has multiple columns", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, matrix(1:4, nrow=2), as.matrix(1:4), - as.matrix(1:4)), - "Y must have only 1 column") + expect_error( + prestogp_fit( + model, matrix(1:4, nrow = 2), as.matrix(1:4), + as.matrix(1:4) + ), + "Y must have only 1 column" + ) }) test_that("nrow(Y) != nrow(locs)", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:4), as.matrix(1:4), - as.matrix(1:3)), - "Y must have the same number of rows as locs") + expect_error( + prestogp_fit( + model, as.matrix(1:4), as.matrix(1:4), + as.matrix(1:3) + ), + "Y must have the same number of rows as locs" + ) }) test_that("nrow(Y) != nrow(X)", { model <- new("VecchiaModel") - expect_error(prestogp_fit(model, as.matrix(1:4), as.matrix(1:3), - as.matrix(1:4)), - "Y must have the same number of rows as X") + expect_error( + prestogp_fit( + model, as.matrix(1:4), as.matrix(1:3), + as.matrix(1:4) + ), + "Y must have the same number of rows as X" + ) }) test_that("Simulated dataset spatial", { - source("sim_vecchia.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=25) - pgp.model1 <- prestogp_fit(pgp.model1, y, X, locs, - scaling=c(1,1), apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - beta.out <- as.vector(pgp.model1@beta) - params.out <- pgp.model1@covparams + source("sim_vecchia.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 25) + pgp.model1 <- prestogp_fit(pgp.model1, y, X, locs, + scaling = c(1, 1), apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + beta.out <- as.vector(pgp.model1@beta) + params.out <- pgp.model1@covparams - expect_length(beta.out, 11) - expect_length(params.out, 5) - expect_equal(beta.out, c(0.01, 0.86, 0.98, 0.94, 0.9, rep(0, 6)), - tolerance=0.03) - expect_equal(params.out[1], 1.6, tolerance=0.5) - expect_equal(params.out[2], 0.4, tolerance=0.2) - expect_equal(params.out[3], 0.59, tolerance=0.2) - expect_equal(params.out[4], 2.0, tolerance=0.15) + expect_length(beta.out, 11) + expect_length(params.out, 5) + expect_equal(beta.out, c(0.01, 0.86, 0.98, 0.94, 0.9, rep(0, 6)), + tolerance = 0.03 + ) + expect_equal(params.out[1], 1.6, tolerance = 0.5) + expect_equal(params.out[2], 0.4, tolerance = 0.2) + expect_equal(params.out[3], 0.59, tolerance = 0.2) + expect_equal(params.out[4], 2.0, tolerance = 0.15) - pgp.model2 <- new("FullModel") - pgp.model2 <- prestogp_fit(pgp.model2, y, X, locs, - scaling=c(1,1), apanasovich=TRUE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - beta.out2 <- as.vector(pgp.model2@beta) - params.out2 <- pgp.model2@covparams + pgp.model2 <- new("FullModel") + pgp.model2 <- prestogp_fit(pgp.model2, y, X, locs, + scaling = c(1, 1), apanasovich = TRUE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + beta.out2 <- as.vector(pgp.model2@beta) + params.out2 <- pgp.model2@covparams - expect_length(beta.out2, 11) - expect_length(params.out2, 5) - expect_equal(beta.out2, c(0.01, 0.86, 0.98, 0.95, 0.9, rep(0, 6)), - tolerance=0.03) - expect_equal(params.out2[1], 1.5, tolerance=0.6) - expect_equal(params.out2[2], 0.4, tolerance=0.15) - expect_equal(params.out2[3], 0.62, tolerance=0.2) - expect_equal(params.out2[4], 2.0, tolerance=0.15) + expect_length(beta.out2, 11) + expect_length(params.out2, 5) + expect_equal(beta.out2, c(0.01, 0.86, 0.98, 0.95, 0.9, rep(0, 6)), + tolerance = 0.03 + ) + expect_equal(params.out2[1], 1.5, tolerance = 0.6) + expect_equal(params.out2[2], 0.4, tolerance = 0.15) + expect_equal(params.out2[3], 0.62, tolerance = 0.2) + expect_equal(params.out2[4], 2.0, tolerance = 0.15) - # Vecchia and full models should be approximately equal - expect_equal(beta.out[1], beta.out2[1], tolerance=0.07) - expect_equal(beta.out[-1], beta.out2[-1], tolerance=0.04) - expect_equal(params.out[1], params.out2[1], tolerance=1) - expect_equal(params.out[2]-params.out2[2], 0, tolerance=0.2) - expect_equal(params.out[3], params.out2[3], tolerance=0.3) - expect_equal(params.out[4], params.out2[4], tolerance=0.2) + # Vecchia and full models should be approximately equal + expect_equal(beta.out[1], beta.out2[1], tolerance = 0.07) + expect_equal(beta.out[-1], beta.out2[-1], tolerance = 0.04) + expect_equal(params.out[1], params.out2[1], tolerance = 1) + expect_equal(params.out[2] - params.out2[2], 0, tolerance = 0.2) + expect_equal(params.out[3], params.out2[3], tolerance = 0.3) + expect_equal(params.out[4], params.out2[4], tolerance = 0.2) }) test_that("Simulated dataset spatiotemporal", { - source("sim_vecchia_st.R") - pgp.model1 <- new("VecchiaModel", n_neighbors=25) - pgp.model1 <- prestogp_fit(pgp.model1, y, X, locs, scaling=c(1,1,2), - apanasovich=FALSE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - beta.out <- as.vector(pgp.model1@beta) - params.out <- pgp.model1@covparams + source("sim_vecchia_st.R") + pgp.model1 <- new("VecchiaModel", n_neighbors = 25) + pgp.model1 <- prestogp_fit(pgp.model1, y, X, locs, + scaling = c(1, 1, 2), + apanasovich = FALSE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + beta.out <- as.vector(pgp.model1@beta) + params.out <- pgp.model1@covparams - expect_length(beta.out, 11) - expect_length(params.out, 6) - expect_equal(beta.out, c(0.01, 0.93, 1.01, 0.91, 0.99, rep(0, 6)), - tolerance=0.015) - expect_equal(params.out[1], 1.7, tolerance=0.55) - expect_equal(params.out[2]-0.19, 0, tolerance=0.05) - expect_equal(params.out[3], 0.22, tolerance=0.05) - expect_equal(params.out[4], 1.12, tolerance=0.3) - expect_equal(params.out[5], 0.62, tolerance=0.05) + expect_length(beta.out, 11) + expect_length(params.out, 6) + expect_equal(beta.out, c(0.01, 0.93, 1.01, 0.91, 0.99, rep(0, 6)), + tolerance = 0.015 + ) + expect_equal(params.out[1], 1.7, tolerance = 0.55) + expect_equal(params.out[2] - 0.19, 0, tolerance = 0.05) + expect_equal(params.out[3], 0.22, tolerance = 0.05) + expect_equal(params.out[4], 1.12, tolerance = 0.3) + expect_equal(params.out[5], 0.62, tolerance = 0.05) - pgp.model2 <- new("FullModel", n_neighbors=25) - pgp.model2 <- prestogp_fit(pgp.model2, y, X, locs, scaling=c(1,1,2), - apanasovich=FALSE, verbose=FALSE, - optim.control=list(trace=0,maxit=5000, - reltol=1e-3)) - beta.out2 <- as.vector(pgp.model2@beta) - params.out2 <- pgp.model2@covparams + pgp.model2 <- new("FullModel", n_neighbors = 25) + pgp.model2 <- prestogp_fit(pgp.model2, y, X, locs, + scaling = c(1, 1, 2), + apanasovich = FALSE, verbose = FALSE, + optim.control = list( + trace = 0, maxit = 5000, + reltol = 1e-3 + ) + ) + beta.out2 <- as.vector(pgp.model2@beta) + params.out2 <- pgp.model2@covparams - expect_length(beta.out2, 11) - expect_length(params.out2, 6) - expect_equal(beta.out2, c(-0.03, 0.93, 1, 0.91, 0.98, rep(0, 6)), - tolerance=0.02) - expect_equal(params.out2[1], 1.6, tolerance=0.5) - expect_equal(params.out2[2]-0.19, 0, tolerance=0.05) - expect_equal(params.out2[3]-0.22, 0, tolerance=0.05) - expect_equal(params.out2[4], 1.18, tolerance=0.15) - expect_equal(params.out2[5], 0.64, tolerance=0.05) + expect_length(beta.out2, 11) + expect_length(params.out2, 6) + expect_equal(beta.out2, c(-0.03, 0.93, 1, 0.91, 0.98, rep(0, 6)), + tolerance = 0.02 + ) + expect_equal(params.out2[1], 1.6, tolerance = 0.5) + expect_equal(params.out2[2] - 0.19, 0, tolerance = 0.05) + expect_equal(params.out2[3] - 0.22, 0, tolerance = 0.05) + expect_equal(params.out2[4], 1.18, tolerance = 0.15) + expect_equal(params.out2[5], 0.64, tolerance = 0.05) - # Vecchia and full models should be approximately equal - expect_equal(beta.out[1], beta.out2[1], tolerance=0.08) - expect_equal(beta.out[-1], beta.out2[-1], tolerance=0.03) - expect_equal(params.out[1], params.out2[1], tolerance=0.6) - expect_equal(params.out[2], params.out2[2], tolerance=0.06) - expect_equal(params.out[3]-params.out2[3], 0, tolerance=0.06) - expect_equal(params.out[4], params.out2[4], tolerance=0.3) - expect_equal(params.out[5], params.out2[5], tolerance=0.1) + # Vecchia and full models should be approximately equal + expect_equal(beta.out[1], beta.out2[1], tolerance = 0.08) + expect_equal(beta.out[-1], beta.out2[-1], tolerance = 0.03) + expect_equal(params.out[1], params.out2[1], tolerance = 0.6) + expect_equal(params.out[2], params.out2[2], tolerance = 0.06) + expect_equal(params.out[3] - params.out2[3], 0, tolerance = 0.06) + expect_equal(params.out[4], params.out2[4], tolerance = 0.3) + expect_equal(params.out[5], params.out2[5], tolerance = 0.1) }) From c37ad44393a2539005c65dde3403d49ba31ce4ea Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Wed, 3 Jan 2024 14:36:09 -0500 Subject: [PATCH 3/9] workaround for lintr bug: https://github.com/REditorSupport/languageserver/issues/89 --- .lintr | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/.lintr b/.lintr index 16c0bd2..46b8c0c 100644 --- a/.lintr +++ b/.lintr @@ -1,5 +1,4 @@ linters: linters_with_defaults( line_length_linter = NULL, - commented_code_linter = NULL -) + commented_code_linter = NULL) encoding: "UTF-8" From 4a8a7a7b5fa01e62dcc0115ced9fd4e117274218 Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Wed, 3 Jan 2024 14:49:23 -0500 Subject: [PATCH 4/9] Remove line length linter to see the remaining errors/warnings --- .github/workflows/lint.yaml | 2 +- .lintr | 4 +++- inst/doc/test.R | 1 - 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml index cf3d710..bc1531a 100644 --- a/.github/workflows/lint.yaml +++ b/.github/workflows/lint.yaml @@ -27,7 +27,7 @@ jobs: - name: Lint run: lintr::lint_package(linters = lintr::linters_with_defaults(object_name_linter = NULL, - commented_code_linter = NULL, cyclocomp_linter = NULL)) + commented_code_linter = NULL, cyclocomp_linter = NULL, line_length_linter = NULL)) shell: Rscript {0} env: LINTR_ERROR_ON_LINT: true diff --git a/.lintr b/.lintr index 46b8c0c..323f6fb 100644 --- a/.lintr +++ b/.lintr @@ -1,4 +1,6 @@ linters: linters_with_defaults( line_length_linter = NULL, - commented_code_linter = NULL) + commented_code_linter = NULL, + object_name_linter = NULL, + cyclocomp_linter = NULL) encoding: "UTF-8" diff --git a/inst/doc/test.R b/inst/doc/test.R index 4c0587d..5992b7d 100644 --- a/inst/doc/test.R +++ b/inst/doc/test.R @@ -1,3 +1,2 @@ ## ----------------------------------------------------------------------------- cat("test", "\n") - From 766f17d0e899c3f3800c5abcc46058b9b6fb8339 Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Wed, 3 Jan 2024 15:30:19 -0500 Subject: [PATCH 5/9] remove linter options from lint action as we have added .lintr project file. Fix all vector_logic_linter warnings --- .github/workflows/lint.yaml | 3 +-- R/Log_Likelihood.R | 2 +- R/PrestoGP_CreateU_Multivariate.R | 6 +++--- R/PrestoGP_Util_Functions.R | 4 ++-- 4 files changed, 7 insertions(+), 8 deletions(-) diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml index bc1531a..f71964d 100644 --- a/.github/workflows/lint.yaml +++ b/.github/workflows/lint.yaml @@ -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, line_length_linter = NULL)) + run: lintr::lint_package(linters = lintr::linters_with_defaults()) shell: Rscript {0} env: LINTR_ERROR_ON_LINT: true diff --git a/R/Log_Likelihood.R b/R/Log_Likelihood.R index 7529cd9..4e628e6 100755 --- a/R/Log_Likelihood.R +++ b/R/Log_Likelihood.R @@ -154,7 +154,7 @@ mvnegloglik_ST <- function(logparams, vecchia.approx, y, param.seq, P, scaling, for (j in 1:nscale) { locs.scaled[vecchia.approx$ondx == i, scaling == j] <- locs.scaled[vecchia.approx$ondx == i, scaling == j] / - params[param.seq[2, 1] + nscale * (i - 1) + j - 1] + params[param.seq[2, 1] + nscale * (i - 1) + j - 1] } } vecchia.approx$locsord <- locs.scaled diff --git a/R/PrestoGP_CreateU_Multivariate.R b/R/PrestoGP_CreateU_Multivariate.R index fdbbe24..78b698f 100755 --- a/R/PrestoGP_CreateU_Multivariate.R +++ b/R/PrestoGP_CreateU_Multivariate.R @@ -160,7 +160,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 } @@ -235,7 +235,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 @@ -264,7 +264,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( diff --git a/R/PrestoGP_Util_Functions.R b/R/PrestoGP_Util_Functions.R index 0a7db69..497a8db 100644 --- a/R/PrestoGP_Util_Functions.R +++ b/R/PrestoGP_Util_Functions.R @@ -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) } From 40b0305a66b443b20b9f39d4cab639ca8c1b3d3f Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Wed, 3 Jan 2024 15:35:05 -0500 Subject: [PATCH 6/9] Add ignore rules for .lintr and .github for R build --- .Rbuildignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.Rbuildignore b/.Rbuildignore index 7fec1b8..c134606 100755 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -5,3 +5,5 @@ ^\.Rproj\.user$ ^doc$ ^Meta$ +^\.lintr +^\.github \ No newline at end of file From 86d7d512c05a2b5ffb5f89ddc2139250252c8f02 Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Wed, 3 Jan 2024 15:40:37 -0500 Subject: [PATCH 7/9] WIP - linting --- .github/workflows/lint.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/lint.yaml b/.github/workflows/lint.yaml index f71964d..d7d07a4 100644 --- a/.github/workflows/lint.yaml +++ b/.github/workflows/lint.yaml @@ -26,7 +26,7 @@ jobs: needs: lint - name: Lint - run: lintr::lint_package(linters = lintr::linters_with_defaults()) + run: lintr::lint_package() shell: Rscript {0} env: LINTR_ERROR_ON_LINT: true From 209a21b17da49e5c5e7929333eb2c64987469972 Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Thu, 11 Jan 2024 13:19:13 -0500 Subject: [PATCH 8/9] Rerun auto-formatter --- R/Log_Likelihood.R | 42 ++++--- R/PrestoGP_CreateU_Multivariate.R | 33 +++-- R/PrestoGP_Multivariate_Vecchia.R | 9 +- R/PrestoGP_Vecchia.R | 28 ++--- R/RcppExports.R | 5 +- src/RcppExports.cpp | 83 ++++++++----- src/createUMC.cpp | 105 +++++++--------- tests/testthat/sim_multivariate_small_pred.R | 73 +++++------ tests/testthat/sim_vecchia_pred.R | 34 ++--- tests/testthat/sim_vecchia_small.R | 24 ++-- tests/testthat/sim_vecchia_small_pred.R | 34 ++--- tests/testthat/test-Log_Likelihood.R | 124 +++++++++++++++---- 12 files changed, 353 insertions(+), 241 deletions(-) diff --git a/R/Log_Likelihood.R b/R/Log_Likelihood.R index 06fabec..4899262 100755 --- a/R/Log_Likelihood.R +++ b/R/Log_Likelihood.R @@ -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 @@ -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 @@ -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)) } @@ -160,7 +168,7 @@ mvnegloglik_ST <- function(logparams, vecchia.approx, y, param.seq, P, scaling, for (j in 1:nscale) { locs.scaled[vecchia.approx$ondx == i, scaling == j] <- locs.scaled[vecchia.approx$ondx == i, scaling == j] / - params[param.seq[2, 1] + nscale * (i - 1) + j - 1] + params[param.seq[2, 1] + nscale * (i - 1) + j - 1] } } vecchia.approx$locsord <- locs.scaled @@ -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 } @@ -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] + ) } diff --git a/R/PrestoGP_CreateU_Multivariate.R b/R/PrestoGP_CreateU_Multivariate.R index c8ba99c..f4b9c6f 100755 --- a/R/PrestoGP_CreateU_Multivariate.R +++ b/R/PrestoGP_CreateU_Multivariate.R @@ -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] } diff --git a/R/PrestoGP_Multivariate_Vecchia.R b/R/PrestoGP_Multivariate_Vecchia.R index 08dab5b..e0d86b2 100755 --- a/R/PrestoGP_Multivariate_Vecchia.R +++ b/R/PrestoGP_Multivariate_Vecchia.R @@ -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) }) diff --git a/R/PrestoGP_Vecchia.R b/R/PrestoGP_Vecchia.R index ebd28ec..0af98b4 100755 --- a/R/PrestoGP_Vecchia.R +++ b/R/PrestoGP_Vecchia.R @@ -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) @@ -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) }) diff --git a/R/RcppExports.R b/R/RcppExports.R index dd16aaf..0e90b6b 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -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) } - diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 3f445f7..da9c5b3 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -1,55 +1,70 @@ // Generated by using Rcpp::compileAttributes() -> do not edit by hand // Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 -#include #include +#include using namespace Rcpp; #ifdef RCPP_USE_GLOBAL_ROSTREAM -Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); -Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); +Rcpp::Rostream &Rcpp::Rcout = Rcpp::Rcpp_cout_get(); +Rcpp::Rostream &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::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::type olocs(olocsSEXP); + Rcpp::traits::input_parameter::type ondx(ondxSEXP); + Rcpp::traits::input_parameter::type curqys(curqysSEXP); + Rcpp::traits::input_parameter::type curqzs(curqzsSEXP); + Rcpp::traits::input_parameter::type vijs(vijsSEXP); + Rcpp::traits::input_parameter::type aijs(aijsSEXP); + Rcpp::traits::input_parameter::type full_const( + full_constSEXP); + Rcpp::traits::input_parameter::type nugget(nuggetSEXP); + Rcpp::traits::input_parameter::type sig2(sig2SEXP); + Rcpp::traits::input_parameter::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); } diff --git a/src/createUMC.cpp b/src/createUMC.cpp index 24b3a7a..213581b 100644 --- a/src/createUMC.cpp +++ b/src/createUMC.cpp @@ -4,59 +4,47 @@ // [[Rcpp::plugins(cpp17)]] #include -double matern1(const double &x, const double &smooth, const double &alpha) -{ - if (smooth == 0.5) - { +double matern1(const double &x, const double &smooth, const double &alpha) { + if (smooth == 0.5) { return exp(-x * alpha); - } - else if (smooth == 1.5) - { + } else if (smooth == 1.5) { double normcon = (pow(alpha, 3) / sqrt(M_PI)) * (tgamma(2.0) / tgamma(1.5)); - return normcon * 0.5 * M_PI * pow(alpha, -3) * exp(-x * alpha) * (1 + x * alpha); - } - else if (smooth == 2.5) - { + return normcon * 0.5 * M_PI * pow(alpha, -3) * exp(-x * alpha) * + (1 + x * alpha); + } else if (smooth == 2.5) { double normcon = (pow(alpha, 5) / sqrt(M_PI)) * (tgamma(3.0) / tgamma(2.5)); double scaled_x = x * alpha; - return normcon * 0.125 * M_PI * pow(alpha, -5) * exp(-scaled_x) * (3 + 3 * scaled_x + scaled_x * scaled_x); - } - else if (x == 0.0) - { + return normcon * 0.125 * M_PI * pow(alpha, -5) * exp(-scaled_x) * + (3 + 3 * scaled_x + scaled_x * scaled_x); + } else if (x == 0.0) { return 1.0; - } - else - { + } else { double normcon = pow(2.0, 1.0 - smooth) / tgamma(smooth); double scaled_x = x * alpha; return normcon * pow(scaled_x, smooth) * R::bessel_k(scaled_x, smooth, 1); } } -// computes Euclidean distance between 2 vectors ([slightly slower] helper function for rdist in C++) -double dist1(const arma::rowvec &a, const arma::rowvec &b) -{ +// computes Euclidean distance between 2 vectors ([slightly slower] helper +// function for rdist in C++) +double dist1(const arma::rowvec &a, const arma::rowvec &b) { double temp = 0.0; // loop over each element and square the distance, add to the placeholder - for (arma::uword i = 0; i < a.n_elem; i++) - { + for (arma::uword i = 0; i < a.n_elem; i++) { temp += (a(i) - b(i)) * (a(i) - b(i)); } return sqrt(temp); } // [[Rcpp::export]] -arma::vec na_omit_c(arma::vec x) -{ +arma::vec na_omit_c(arma::vec x) { // make placeholder vector r arma::vec r(x.size()); // make counter for number of non-NAs int k = 0; - for (unsigned j = 0; j < x.size(); j++) - { + for (unsigned j = 0; j < x.size(); j++) { // if not NA, add 1 to counter and set r(j) = x(j) - if (x(j) == x(j)) - { + if (x(j) == x(j)) { r(j) = x(j); k++; } @@ -70,11 +58,12 @@ arma::vec na_omit_c(arma::vec x) // [[Rcpp::export]] 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) -{ + 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) { int n = arma::as_scalar(ondx.n_elem); int m = arma::as_scalar(curqys.n_rows); int n_inds = 2 * n * (m + 3); @@ -85,16 +74,14 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, // vals.rows(0, 6) = U_beginning; // feeder.cols(0, 6) = {{1, 1, 2, 1, 3, 3, 4}, // {1, 2, 2, 3, 3, 4, 4}}; - ndx.cols(0, 6) = {{0, 0, 1, 0, 2, 2, 3}, - {0, 1, 1, 2, 2, 3, 3}}; + ndx.cols(0, 6) = {{0, 0, 1, 0, 2, 2, 3}, {0, 1, 1, 2, 2, 3, 3}}; vals(arma::span(0, 6)) = U_beginning; // feeder.cols(0,6) = {{0, 0, 1, 0, 2, 2, 3}, // {0, 1, 1, 2, 2, 3, 3}, {1,1,1,1,1,1,1}}; // feeder(2, arma::span(0,6)) = U_beginning.t(); int ind = 7; // int ind = 0; - for (arma::uword i = 3; i < (ondx.n_elem + 1); i++) - { + for (arma::uword i = 3; i < (ondx.n_elem + 1); i++) { double temppi = ondx(i - 1); arma::vec cy = na_omit_c(curqys.col(i - 1)); arma::vec cz = (curqzs.col(i - 1)); @@ -104,22 +91,21 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, arma::mat k2(nq, nq); // Rcpp::Rcout << i << " "; // Rcpp::Rcout << cq; - for (arma::uword j = 0; j < nq; j++) - { + for (arma::uword j = 0; j < nq; j++) { double temppj = ondx(cq(j) - 1); - k1(j) = full_const(temppi - 1, temppj - 1) * - matern1(dist1(olocs.row(i - 1), olocs.row(cq(j) - 1)), - vijs(temppi - 1, temppj - 1), aijs(temppi - 1, temppj - 1)); - for (arma::uword k = j; k < nq; k++) - { + k1(j) = + full_const(temppi - 1, temppj - 1) * + matern1(dist1(olocs.row(i - 1), olocs.row(cq(j) - 1)), + vijs(temppi - 1, temppj - 1), aijs(temppi - 1, temppj - 1)); + for (arma::uword k = j; k < nq; k++) { double temppk = ondx(cq(k) - 1); - k2(j, k) = full_const(temppj - 1, temppk - 1) * - matern1(dist1(olocs.row(cq(j) - 1), olocs.row(cq(k) - 1)), - vijs(temppj - 1, temppk - 1), aijs(temppj - 1, temppk - 1)); + k2(j, k) = + full_const(temppj - 1, temppk - 1) * + matern1(dist1(olocs.row(cq(j) - 1), olocs.row(cq(k) - 1)), + vijs(temppj - 1, temppk - 1), aijs(temppj - 1, temppk - 1)); } // // add nugget if neighbor is not latent - if (j >= cy.n_elem) - { + if (j >= cy.n_elem) { k2(j, j) += nugget(temppj - 1); } } @@ -131,16 +117,14 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, arma::vec bi(nq); bool status = arma::solve(bi, k2, k1); double epsilon = 1e-9; - while (!status) - { + while (!status) { epsilon *= 2; k2.diag() += epsilon; status = arma::solve(bi, k2, k1); } double ritemp = arma::as_scalar(sig2(temppi - 1) - bi.t() * k1); epsilon = 1e-9; - while (ritemp <= 0) - { + while (ritemp <= 0) { epsilon *= 2; k2.diag() += epsilon; arma::solve(bi, k2, k1); @@ -154,11 +138,9 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, // uhat(2*i - 2, 2*i - 1) = -1.0 * uhat(2*i - 1, 2*i - 1); // uhat(2*i - 2, 2*i - 2) = 1.0 / ri; - for (arma::uword j = 0; j < cq.n_elem; j++) - { + for (arma::uword j = 0; j < cq.n_elem; j++) { arma::uword xind = 2 * cq(j) - 2; - if (j >= cy.n_elem) - { + if (j >= cy.n_elem) { xind += 1; } // feeder.col(ind) = {(double)xind,(double) 2*i - 2, bi(j)}; @@ -169,7 +151,9 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, // uhat(xind, 2.0*i - 2.0) = bi(j); } // vals[ind] = 1.0 / ri; - double feed_temp = pow(nugget(temppi - 1), -0.5); // used twice so temp storage to avoid recomputation + double feed_temp = + pow(nugget(temppi - 1), + -0.5); // used twice so temp storage to avoid recomputation // feeder.col(ind) = {(double)2*i - 2,(double)2*i - 2, 1.0 / ri}; // feeder.col(ind + 2) = {(double)2*i - 1,(double) 2*i - 1, feed_temp}; ndx.col(ind) = {2 * i - 2, 2 * i - 2}; @@ -178,7 +162,8 @@ arma::sp_mat createU_helper_mat(const arma::mat &olocs, const arma::vec &ondx, vals(ind + 2) = feed_temp; // vals[ind + 2] = pow(nugget(temppi - 1), -0.5); // vals[ind + 1] = -1.0 * vals[ind + 2]; - // feeder.col(ind + 1) = {(double)2*i - 2,(double) 2*i - 1, -1.0 * feed_temp}; + // feeder.col(ind + 1) = {(double)2*i - 2,(double) 2*i - 1, -1.0 * + // feed_temp}; ndx.col(ind + 1) = {2 * i - 2, 2 * i - 1}; vals(ind + 1) = -1.0 * feed_temp; ind += 3; diff --git a/tests/testthat/sim_multivariate_small_pred.R b/tests/testthat/sim_multivariate_small_pred.R index 0b5101e..0fda2db 100755 --- a/tests/testthat/sim_multivariate_small_pred.R +++ b/tests/testthat/sim_multivariate_small_pred.R @@ -9,20 +9,20 @@ library(MASS) library(fields) library(psych) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) beta.all <- rep(beta1, ny) X.st <- list() for (i in 1:ny) { - Sigma.X <- exp(-rdist(sample(1:p))/3) - X.st[[i]] <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) + Sigma.X <- exp(-rdist(sample(1:p)) / 3) + X.st[[i]] <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) } X.all <- superMatrix(X.st) mean.trend.st <- X.all %*% beta.all n.rho <- choose(ny, 2) rho.vec <- runif(n.rho, 0.2, 0.8) -rho <- matrix(0, nrow=ny, ncol=ny) +rho <- matrix(0, nrow = ny, ncol = ny) rho[upper.tri(rho)] <- rho.vec rho <- rho + t(rho) + diag(1, ny) @@ -36,36 +36,37 @@ params.all <- c(x.variance, ranges, marg.smoothness, nuggets, rho.vec) locs.list <- list() for (i in 1:ny) { - loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) - loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) + loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) + loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs.list[[i]] <- as.matrix(expand.grid(loc1, loc2)) } -Sigma.All <- matrix(nrow=ny*n.spatial.xy^2, ncol=ny*n.spatial.xy^2) +Sigma.All <- matrix(nrow = ny * n.spatial.xy^2, ncol = ny * n.spatial.xy^2) for (i in 1:ny) { for (j in i:ny) { - if (i==j) { - ndx1 <- ((i-1)*n.spatial.xy^2+1):(n.spatial.xy^2*i) + if (i == j) { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) dij <- fields::rdist(locs.list[[i]]) - Sigma.All[ndx1,ndx1] <- marg.var[i] * - fields::Matern(dij, range = ranges[i], - smoothness = marg.smoothness[i]) - } - else { - ndx1 <- ((i-1)*n.spatial.xy^2+1):(n.spatial.xy^2*i) - ndx2 <- ((j-1)*n.spatial.xy^2+1):(n.spatial.xy^2*j) + Sigma.All[ndx1, ndx1] <- marg.var[i] * + fields::Matern(dij, + range = ranges[i], + smoothness = marg.smoothness[i] + ) + } else { + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) + ndx2 <- ((j - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * j) dij <- fields::rdist(locs.list[[i]], locs.list[[j]]) vii <- marg.smoothness[i] vjj <- marg.smoothness[j] - vij <- (vii+vjj)/2 - aii <- 1/ranges[i] - ajj <- 1/ranges[j] - aij <- sqrt((aii^2+ajj^2)/2) - Sigma.All[ndx1,ndx2] <- rho[i,j] * sqrt(marg.var[i]) * + vij <- (vii + vjj) / 2 + aii <- 1 / ranges[i] + ajj <- 1 / ranges[j] + aij <- sqrt((aii^2 + ajj^2) / 2) + Sigma.All[ndx1, ndx2] <- rho[i, j] * sqrt(marg.var[i]) * sqrt(marg.var[j]) * aii^vii * ajj^vjj * gamma(vij) / - (aij^(2*vij) * sqrt(gamma(vii) * gamma(vjj))) * - fields::Matern(dij, smoothness=vij, alpha=aij) - Sigma.All[ndx2,ndx1] <- t(Sigma.All[ndx1,ndx2]) + (aij^(2 * vij) * sqrt(gamma(vii) * gamma(vjj))) * + fields::Matern(dij, smoothness = vij, alpha = aij) + Sigma.All[ndx2, ndx1] <- t(Sigma.All[ndx1, ndx2]) } } } @@ -76,7 +77,7 @@ st.error <- rnorm(n.spatial.xy^2 * ny) %*% L.C nug.error <- NULL for (i in 1:ny) { - nug.error <- c(nug.error, nuggets[i]*rnorm(n.spatial.xy^2)) + nug.error <- c(nug.error, nuggets[i] * rnorm(n.spatial.xy^2)) } y.list <- list() @@ -87,20 +88,22 @@ y.list.otst <- list() X.st.otr <- list() X.st.otst <- list() for (i in 1:ny) { - ndx1 <- ((i-1)*n.spatial.xy^2+1):(n.spatial.xy^2*i) + ndx1 <- ((i - 1) * n.spatial.xy^2 + 1):(n.spatial.xy^2 * i) y.list[[i]] <- mean.trend.st[ndx1] + st.error[ndx1] + nug.error[ndx1] nn <- n.spatial.xy^2 otr <- rep(FALSE, nn) - otr[sample(1:nn, size=floor(nn/2))] <- TRUE - locs.list.otr[[i]] <- locs.list[[i]][otr,] - locs.list.otst[[i]] <- locs.list[[i]][!otr,] + otr[sample(1:nn, size = floor(nn / 2))] <- TRUE + locs.list.otr[[i]] <- locs.list[[i]][otr, ] + locs.list.otst[[i]] <- locs.list[[i]][!otr, ] y.list.otr[[i]] <- y.list[[i]][otr] y.list.otst[[i]] <- y.list[[i]][!otr] - X.st.otr[[i]] <- X.st[[i]][otr,] - X.st.otst[[i]] <- X.st[[i]][!otr,] + X.st.otr[[i]] <- X.st[[i]][otr, ] + X.st.otst[[i]] <- X.st[[i]][!otr, ] } -rm(ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, - loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, - st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, - marg.smoothness, marg.var, x.variance, nn, otr) +rm( + ny, p, p.nz, n.spatial.xy, beta1, i, j, Sigma.X, mean.trend.st, n.rho, + loc1, loc2, ndx1, ndx2, dij, vii, vjj, vij, aii, ajj, aij, L.C, + st.error, nug.error, X.all, rho, rho.vec, ranges, Sigma.All, nuggets, + marg.smoothness, marg.var, x.variance, nn, otr +) diff --git a/tests/testthat/sim_vecchia_pred.R b/tests/testthat/sim_vecchia_pred.R index ffcb97e..dbda626 100755 --- a/tests/testthat/sim_vecchia_pred.R +++ b/tests/testthat/sim_vecchia_pred.R @@ -7,10 +7,10 @@ n.spatial.xy <- 50 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,29 +22,33 @@ ranges <- runif(1, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) -Sigma.All <- marg.var*Matern(rdist(locs), range=ranges, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locs), + range = ranges, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^2) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^2) +nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error nn <- n.spatial.xy^2 otr <- rep(FALSE, nn) -otr[sample(1:nn, size=floor(nn/2))] <- TRUE +otr[sample(1:nn, size = floor(nn / 2))] <- TRUE -locs.otr <- locs[otr,] -locs.otst <- locs[!otr,] +locs.otr <- locs[otr, ] +locs.otst <- locs[!otr, ] y.otr <- y[otr] y.otst <- y[!otr] -X.otr <- X[otr,] -X.otst <- X[!otr,] +X.otr <- X[otr, ] +X.otst <- X[!otr, ] -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance, nn, otr) + x.variance, nn, otr +) diff --git a/tests/testthat/sim_vecchia_small.R b/tests/testthat/sim_vecchia_small.R index 9c884eb..39905dc 100755 --- a/tests/testthat/sim_vecchia_small.R +++ b/tests/testthat/sim_vecchia_small.R @@ -7,10 +7,10 @@ n.spatial.xy <- 10 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,18 +22,22 @@ ranges <- runif(1, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) -Sigma.All <- marg.var*Matern(rdist(locs), range=ranges, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locs), + range = ranges, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^2) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^2) +nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance) + x.variance +) diff --git a/tests/testthat/sim_vecchia_small_pred.R b/tests/testthat/sim_vecchia_small_pred.R index 7f25507..36ba071 100755 --- a/tests/testthat/sim_vecchia_small_pred.R +++ b/tests/testthat/sim_vecchia_small_pred.R @@ -7,10 +7,10 @@ n.spatial.xy <- 20 # number of spatial coordinates per dimension library(MASS) library(fields) -beta1 <- c(rep(1,p.nz), rep(0,p-p.nz)) +beta1 <- c(rep(1, p.nz), rep(0, p - p.nz)) -Sigma.X <- exp(-rdist(sample(1:p))/3) -X <- mvrnorm(n.spatial.xy^2,rep(0,p),Sigma.X) +Sigma.X <- exp(-rdist(sample(1:p)) / 3) +X <- mvrnorm(n.spatial.xy^2, rep(0, p), Sigma.X) mean.trend.st <- as.vector(X %*% beta1) marg.smoothness <- 0.5 + rnorm(1, 0.1, 0.05) @@ -22,29 +22,33 @@ ranges <- runif(1, 0.5, 1.2) params.all <- c(x.variance, ranges, marg.smoothness, nuggets) -loc1 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) -loc2 <- seq(0,1,length.out = n.spatial.xy)+rnorm(n.spatial.xy,0,0.001) +loc1 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) +loc2 <- seq(0, 1, length.out = n.spatial.xy) + rnorm(n.spatial.xy, 0, 0.001) locs <- as.matrix(expand.grid(loc1, loc2)) -Sigma.All <- marg.var*Matern(rdist(locs), range=ranges, - smoothness=marg.smoothness) +Sigma.All <- marg.var * Matern(rdist(locs), + range = ranges, + smoothness = marg.smoothness +) L.C <- chol(Sigma.All) st.error <- as.vector(rnorm(n.spatial.xy^2) %*% L.C) -nug.error <- nuggets*rnorm(n.spatial.xy^2) +nug.error <- nuggets * rnorm(n.spatial.xy^2) y <- mean.trend.st + st.error + nug.error nn <- n.spatial.xy^2 otr <- rep(FALSE, nn) -otr[sample(1:nn, size=floor(nn/2))] <- TRUE +otr[sample(1:nn, size = floor(nn / 2))] <- TRUE -locs.otr <- locs[otr,] -locs.otst <- locs[!otr,] +locs.otr <- locs[otr, ] +locs.otst <- locs[!otr, ] y.otr <- y[otr] y.otst <- y[!otr] -X.otr <- X[otr,] -X.otst <- X[!otr,] +X.otr <- X[otr, ] +X.otst <- X[!otr, ] -rm(p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, +rm( + p, p.nz, n.spatial.xy, beta1, Sigma.X, mean.trend.st, loc1, loc2, L.C, st.error, nug.error, ranges, Sigma.All, nuggets, marg.smoothness, marg.var, - x.variance, nn, otr) + x.variance, nn, otr +) diff --git a/tests/testthat/test-Log_Likelihood.R b/tests/testthat/test-Log_Likelihood.R index db5ec55..943d9ac 100755 --- a/tests/testthat/test-Log_Likelihood.R +++ b/tests/testthat/test-Log_Likelihood.R @@ -1,28 +1,74 @@ +test_that("negloglik_vecchia", { + set.seed(1234) + + y.orig <- rnorm(100) + + locs.dim1 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) + locs.dim2 <- seq(1, 10, length = 10) + rnorm(10, 0, 0.1) + locs <- as.matrix(expand.grid(locs.dim1, locs.dim2)) + + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(100) + + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) + + logparams <- create.initial.values.flex(0.9, 0.5, 0.5, 0.1, 1, 1) + pseq <- create.param.sequence(1) + vec.approx <- vecchia_specify(locs, 5) + LL.vecchia <- negloglik_vecchia(logparams, y, vec.approx, pseq) + expect_equal(194.75, LL.vecchia, tolerance = 1e-2) + + vec.mapprox <- vecchia_Mspecify(list(locs), 5) + LL.mvecchia <- mvnegloglik(logparams, vec.mapprox, y, pseq, 1) + # Univariate likelihood should equal multivariate likelihood + expect_equal(LL.vecchia, LL.mvecchia, tolerance = 1e-1) + + logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) + vec.approx2 <- vec.approx + vec.approx2$locsord <- vec.approx$locsord / 0.5 + LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq) + LL.vecchia.st <- negloglik_vecchia_ST(logparams, y, vec.approx, pseq, 1, 1) + vec.mapprox2 <- vec.mapprox + vec.mapprox2$locsord <- vec.mapprox$locsord / 0.5 + LL.mvecchia2 <- mvnegloglik_ST(logparams2, vec.mapprox2, y, pseq, 1, 1, 1) + + # Likelihood should equal spatiotemporal likelihood after scaling locs + expect_equal(LL.vecchia, LL.vecchia2, tolerance = 1e-3) + expect_equal(LL.vecchia, LL.vecchia.st, tolerance = 1e-3) + expect_equal(LL.vecchia, LL.mvecchia2, tolerance = 1e-1) +}) + test_that("negloglik_vecchia_ST", { - set.seed(7919) - load("small_sim.Rdata") - return(1) # this test isn't useful at the moment - params <- c(0.5, 0.5, 0.5, 0.5) - locs <- locs_train - X <- X_train - Y <- Y_train - locs.scaled <- locs / c(params[2], params[2], params[3]) - vecchia_approx <- vecchia_specify(locs.scaled, 25) - beta.hat <- rep(0, ncol(X)) - Y.hat <- as.matrix(X) %*% beta.hat - res <- as.double(Y - Y.hat) - ord_locs <- locs[vecchia_approx$ord, ] - locs_mat <- cbind(ord_locs[, 1], ord_locs[, 2], ord_locs[, 3]) - result <- negloglik_vecchia_ST(log(params), locs_mat, res, vecchia_approx) - expect_equal(1597.9808688, result, tolerance = 10e-5) - vecchia.result <- optim( - par = log(params), fn = negloglik_vecchia_ST, - locs = locs_mat, res = res, vecchia.approx = vecchia_approx, method = "Nelder-Mead", - control = list(trace = 0) + set.seed(1234) + + y.orig <- rnorm(125) + + locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) + + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) + + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) + + logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) + pseq <- create.param.sequence(1, 2) + vec.approx <- vecchia_specify(locs, 5) + LL.vecchia <- negloglik_vecchia_ST( + logparams, y, vec.approx, pseq, + c(1, 1, 2), 2 ) - params_optim <- exp(vecchia.result$par) - result_optim <- negloglik_vecchia_ST(log(params_optim), locs_mat, res, vecchia_approx) - expect_equal(-6861.65048095, result_optim, tolerance = 1e-5) + expect_equal(284.73, LL.vecchia, tolerance = 1e-2) + + logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) + pseq2 <- create.param.sequence(1) + vec.approx2 <- vec.approx + vec.approx2$locsord[, 1:2] <- vec.approx$locsord[, 1:2] / 2 + vec.approx2$locsord[, 3] <- vec.approx$locsord[, 3] / 3 + LL.vecchia2 <- negloglik_vecchia(logparams2, y, vec.approx2, pseq2) + expect_equal(LL.vecchia, LL.vecchia2, tolerance = 1e-3) }) test_that("negloglik.full", { @@ -99,6 +145,38 @@ test_that("negloglik.full", { expect_equal(LL.full, LL.pgp, tolerance = 1e-3) }) +test_that("negloglik_full_ST", { + set.seed(1234) + + y.orig <- rnorm(125) + + locs.dim1 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim2 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs.dim3 <- seq(1, 5, length = 5) + rnorm(5, 0, 0.1) + locs <- as.matrix(expand.grid(locs.dim1, locs.dim2, locs.dim3)) + + covmat.true <- Matern(rdist(locs), range = 1, smoothness = 0.5) + diag(125) + + y <- y.orig %*% chol(covmat.true) + y <- as.vector(y) + + logparams <- create.initial.values.flex(0.9, c(2, 3), 0.5, 0.1, 1, 1) + pseq <- create.param.sequence(1, 2) + LL.full <- negloglik_full_ST( + logparams, locs, y, pseq, + c(1, 1, 2), 2 + ) + expect_equal(286.84, LL.full, tolerance = 1e-2) + + logparams2 <- create.initial.values.flex(0.9, 1, 0.5, 0.1, 1, 1) + pseq2 <- create.param.sequence(1) + locs2 <- locs + locs2[, 1:2] <- locs[, 1:2] / 2 + locs2[, 3] <- locs[, 3] / 3 + LL.full2 <- negloglik.full(logparams2, rdist(locs2), y, pseq2) + expect_equal(LL.full, LL.full2, tolerance = 1e-3) +}) + test_that("mvnegloglik", { load("sim_multivariate_big.RData") set.seed(1234) From 43ae2729489f613a8f559179e26a40377459d5c2 Mon Sep 17 00:00:00 2001 From: sciome-bot Date: Thu, 11 Jan 2024 14:14:09 -0500 Subject: [PATCH 9/9] Add missing comma in imports section of DESCRIPTION --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 9705ed4..a2d1b69 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -45,7 +45,7 @@ Imports: rlang, mvtnorm, spam, - psych + psych, doParallel License: GPL-3 Encoding: UTF-8