#' @include RcppExports.R
NULL
#' Affine transformation
#'
#' @noRd
affineTrans <- function(OldValue, OldMax, OldMin, NewMax, NewMin) {
OldRange <- (OldMax - OldMin)
NewRange <- (NewMax - NewMin)
NewValue <- (((OldValue - OldMin) * NewRange) / OldRange) + NewMin
NewValue
}
#' Normal niche simulation function
#'
#' @noRd
normal_niche <- function(x, y) {
40 * mvtnorm::dmvnorm(
x = matrix(c(x, y), ncol = 2),
mean = c(0, 0),
sigma = matrix(c(8.5, 0, 0, 8.5), ncol = 2)
)
}
#' Uniform niche simulation function
#'
#' @noRd
uniform_niche <- function(x, y) {
rep(0.5, length(x))
}
#' Bimodal niche simulation function
#'
#' @noRd
bimodal_niche <- function(x, y) {
apply(
matrix(
c(
mvtnorm::dmvnorm(
x = matrix(c(x, y), ncol = 2),
mean = c(-2, -2),
sigma = matrix(c(5, 0, 0, 5), ncol = 2)
) * 30,
mvtnorm::dmvnorm(
x = matrix(c(x, y), ncol = 2),
mean = c(3, 3),
sigma = matrix(c(3, 0, 0, 3), ncol = 2)
) * 12
),
ncol = 2),
1,
max
)
}
#' Function to hash function call
#'
#' @noRd
hashCall <- function(expr, skipargs = c(), env = parent.frame()) {
expr <- expr[c((skipargs * -1L) - 1L)]
expr <- expr[which(names(expr) != "force.reset")]
for (i in seq_along(names(expr)))
if (inherits(expr[[i]], c("name")))
expr[[i]] <- eval(expr[[i]], envir = env)
paste(deparse(expr), collapse = ";")
}
#' Make pretty geoplot
#' @noRd
prettyGeoplot <- function(polygons, col, basemap, main, fun, beside = TRUE,
border = NULL, lwd = 1) {
# make layout
defpar <- graphics::par(no.readonly = TRUE)
graphics::par(mar = c(1, 1, 1, 1), oma = c(0, 0, 0, 0))
if (beside) {
graphics::layout(
matrix(c(1, 1, 3, 2), ncol = 2, byrow = TRUE),
widths = c(0.8, 0.2),
heights = c(0.1, 0.9)
)
} else {
graphics::layout(
matrix(c(1, 3, 2), ncol = 1, byrow = TRUE),
widths = c(1),
heights = c(0.1, 0.8, 0.1)
)
}
# convert to list if not
if (!inherits(polygons, "list")) {
polygons <- list(polygons)
col <- list(col)
border <- list(border)
lwd <- list(lwd)
}
# make title
graphics::plot(
1, 1, type = "n", xlim = c(-1, 1), ylim = c(-1, 1),
axes = FALSE, xlab = "", ylab = ""
)
graphics::mtext(side = 1, line = -0.5, main, cex = 1.5)
# make legend
graphics::plot(
1, 1, type = "n", xlim = c(-1, 1), ylim = c(-1, 1),
axes = FALSE, xlab = "", ylab = ""
)
fun()
# make geoplot
if (is.list(basemap)) {
assertthat::assert_that(
requireNamespace("RgoogleMaps", quietly = TRUE),
msg = "please install the \"RgoogleMaps\" package"
)
RgoogleMaps::PlotOnStaticMap(basemap)
for (i in seq_along(polygons)) {
suppressWarnings(RgoogleMaps::PlotPolysOnStaticMap(
basemap, polygons[[i]], col = col[[i]], border = border[[i]],
add = TRUE, lwd = lwd[[i]]))
}
} else {
allpolygons <- do.call(rbind, polygons)
PBSmapping::plotPolys(
polygons[[1]], col = col[[1]], axes = FALSE,
xlab = "", ylab = "", border = border[[1]],
lwd = lwd[[1]], xlim = range(allpolygons$X),
ylim = range(allpolygons$Y)
)
for (i in seq_along(polygons)[-1]) {
if (nrow(polygons[[i]]) > 0) {
suppressWarnings(
PBSmapping::addPolys(
polygons[[i]], col = col[[i]], xlab = "", ylab = "",
border = border[[i]], lwd = lwd[[i]]
)
)
}
}
}
graphics::par(defpar)
invisible()
}
#' Extract colors from RcolorBrewer palette functions
#'
#' @noRd
brewerCols <- function(values, pal, alpha = 1, n = NULL) {
if (is.null(n) & length(pal) == 1) {
idx <- which(rownames(RColorBrewer::brewer.pal.info) == pal)
n <- RColorBrewer::brewer.pal.info$maxcolors[idx]
} else {
n <- length(values)
}
if (length(pal) == 1) {
suppressWarnings({
r <- grDevices::colorRamp(RColorBrewer::brewer.pal(n, pal))(values)
})
} else{
suppressWarnings({
r <- grDevices::colorRamp(pal)(values)
})
}
grDevices::rgb(
r, maxColorValue = 255,
alpha = scales::rescale(alpha, from = c(0, 1), to = c(0, 255))
)
}
#' Add continuous legend to plot
#'
#' @noRd
continuousLegend <- function(values, pal, posx, posy, center = FALSE,
endlabs = NULL) {
function() {
if (center) {
vabs <- max(abs(range(values)))
values <- seq(-vabs, vabs, length.out = 100)
}
xdiff <- diff(graphics::par()$usr[1:2])
ydiff <- diff(graphics::par()$usr[3:4])
zvals <- pretty(values)
zvals <- zvals[which(zvals > min(values) & zvals < max(values))]
if (max(zvals) < 1) {
digit <- 2
} else {
digit <- 1
}
shape::colorlegend(
zlim = range(values), digit = digit,
col = brewerCols(seq(0, 1, length.out = 100), pal),
zval = zvals, posx = posx, posy = posy, xpd = TRUE
)
if (!is.null(endlabs)) {
xcoord <-
graphics::par()$usr[1] +
mean(graphics::par()$usr[2:1] * posx * 2.2)
ycoords <-
(graphics::par()$usr[3] +
diff(graphics::par()$usr[3:4]) * posy) +
(diff(graphics::par()$usr[3:4]) * c(-0.02, 0.02))
graphics::text(
x = rep(xcoord, 2), y = ycoords, rev(endlabs), cex = 1.2,
ad = c(0.5, 0.5)
)
}
}
}
#' Add categorical legend to plot
#'
#' @noRd
categoricalLegend <- function(col, labels, ncol = 1) {
function() {
if (ncol == 1) {
graphics::legend(
"top", bg = "white", legend = labels, fill = col,
bty = "n", horiz = TRUE, cex = 1.5
)
} else {
graphics::legend(
y = graphics::par()$usr[3] + (diff(graphics::par()$usr[3:4]) * 0.6),
x = graphics::par()$usr[1] + (diff(graphics::par()$usr[1:2]) * 0.5),
bg = "white", legend = labels, fill = col, bty = "n",
ncol = ncol, cex = 1.5, xjust = 0.5, yjust = 0.5,
xpd = TRUE
)
}
}
}
#' Create 1-dimensional demand points
#'
#' @noRd
demand.points.density1d <- function(pts, n, quantile = 0.95, ...) {
# transform pts
curr.mean <- mean(pts[, 1])
curr.sd <- stats::sd(pts[, 1])
pts[, 1] <- (pts[, 1] - curr.mean) / curr.sd
# generate points
quants <- stats::quantile(
pts[, 1], c( (1 - quantile) / 2, quantile + (1 - quantile) / 2)
)
dp <- stats::runif(n, quants[[1]], quants[[2]])
# density kernel
est <- ks::kde(pts[, 1], eval.points = dp, ...)
# back-transform demand point coordinates
dp.pts <- (matrix(est$eval.points, ncol = 1) * curr.sd) + curr.mean
# return object
list(coords = dp.pts, weights = est$estimate)
}
#' Create 2-dimensional demand points
#'
#' @noRd
demand.points.density2d <- function(pts, n, quantile = 0.95, ...) {
# transform pts
curr.mean <- apply(pts, 2, mean)
curr.sd <- apply(pts, 2, stats::sd)
pts <- sweep(pts, MARGIN = 2, FUN = "-", curr.mean)
pts <- sweep(pts, MARGIN = 2, FUN = "/", curr.sd)
# generate points
dp <- sp::spsample(
adehabitatHR::mcp(sp::SpatialPoints(coords = pts),
percent = quantile * 100, unin = c("m"), unout = c("m2")),
n * 1.1, type = "random")@coords[seq_len(n), ]
# fit density kernel
est <- ks::kde(pts, eval.points = dp, ...)
# back-transform dps
dp <- sweep(dp, MARGIN = 2, FUN = "*", curr.sd)
dp <- sweep(dp, MARGIN = 2, FUN = "+", curr.mean)
# prepare data to return
list(coords = dp, weights = est$estimate)
}
#' Create n-dimensional demand points
#'
#' @noRd
demand.points.hypervolume <- function(pts, n, quantile = 0.95, ...) {
# transform pts
curr.mean <- apply(pts, 2, mean)
curr.sd <- apply(pts, 2, stats::sd)
pts <- sweep(pts, MARGIN = 2, FUN = "-", curr.mean)
pts <- sweep(pts, MARGIN = 2, FUN = "/", curr.sd)
# fit density kernel
args <- list(...)
if (!"samples.per.point" %in% names(args))
args$samples.per.point <- 500 * ncol(pts)
if (args$samples.per.point * nrow(pts) < n) {
stop(paste0(
"argument to n.demand.points is too high. Set a higher value ",
"in the argument to samples.per.point (defaults to ",
"500 * dimensions in attribute space)."
)
)
}
# estimate bandwidth for kernel
if (!"kde.bandwidth" %in% names(args)) {
args$kde.bandwidth <- hypervolume::estimate_bandwidth(pts)
}
# fit kernel
hv <- do.call(
hypervolume::hypervolume,
append(list(data = pts, quantile.requested = quantile), args)
)
# extract random points
rndpos <- sample.int(nrow(hv@RandomPoints), n)
# extract coordinates and back-transform
dp <- hv@RandomPoints[rndpos,, drop = FALSE]
dp <- sweep(dp, MARGIN = 2, FUN = "*", curr.sd)
dp <- sweep(dp, MARGIN = 2, FUN = "+", curr.mean)
# return object
list(coords = dp, weights = hv@ValueAtRandomPoints[rndpos])
}
#' Calculate zonal-means
#'
#' @noRd
zonalMean <- function(x, y, ids = names(y)) {
assertthat::assert_that(
inherits(x, "SpatRaster"),
inherits(y, "SpatRaster"),
is.character(ids) || is.numeric(ids),
assertthat::noNA(ids)
)
# rename layers for consistency
names(y) <- ids
names(x) <- "pu"
# compute zonal stats
d <- terra::zonal(y, x, fun = "mean")
# convert wide data frame to long format and exclude zeros
d <- do.call(
rbind,
lapply(ids, function(x) {
d$species <- x
r <- d[, c("species", "pu", x), drop = FALSE]
names(r)[3] <- "value"
# remove zeros
r[r$value > 0, , drop = FALSE]
})
)
# coerce integers
if (is.numeric(ids)) {
d$species <- as.integer(d$species)
}
d$pu <- as.integer(d$pu)
# exclude missing values
d <- d[is.finite(d$value), , drop = FALSE]
# return result
d
}
#' Merge list of RapResults into a single object
#'
#' @noRd
mergeRapResults <- function(x) {
x <- RapResults(
summary = do.call(rbind, lapply(x, methods::slot, name = "summary")),
selections = do.call(rbind, lapply(x, methods::slot, name = "selections")),
amount.held = do.call(
rbind,
lapply(x, methods::slot, name = "amount.held")
),
space.held = do.call(
rbind,
lapply(x, methods::slot, name = "space.held")
),
logging.file = sapply(x, methods::slot, name = "logging.file"))
x@summary$Run_Number <- seq_len(nrow(x@summary))
x
}
#' Read RAP results
#'
#' This function reads files output from Gurobi and returns a `RapResults`
#' object.
#'
#' @param opts `RapReliableOpts` or `RapUnreliableOpts` object
#'
#' @param data `RapData` object
#'
#' @param model.list `list` object containing Gurobi model data
#'
#' @param logging.file `character` Gurobi log files.
#'
#' @param solution.list `list` object containing Gurobi solution data.
#'
#' @param verbose `logical` print progress messages? Defaults to
#' `FALSE`.
#'
#' @keywords internal
#'
#' @return `RapResults` object
#'
#' @seealso [RapReliableOpts()], [RapUnreliableOpts()],
#' [RapData()], [RapResults()].
read.RapResults <- function(opts, data, model.list, logging.file,
solution.list, verbose = FALSE) {
x <- rcpp_extract_model_object(
opts, inherits(opts, "RapUnreliableOpts"), data, model.list, logging.file,
solution.list, verbose
)
x@.cache <- new.env()
x
}
#' Compare Rap objects
#'
#' This function checks objects to see if they share the same input data.
#'
#' @param x `RapData`, `RapUnsolved`, or `RapSolved` object.
#'
#' @param y `RapData`, `RapUnsolved`, or `RapSolved` object.
#'
#' @return `logical` are the objects based on the same data?
#'
#' @keywords internal
#'
#' @seealso [RapData-class], [RapUnsolved-class],
#' [RapSolved-class].
#'
#' @name is.comparable
methods::setGeneric("is.comparable",
function(x, y) methods::standardGeneric("is.comparable"))
#' Basemap
#'
#' This function retrieves google map data for planning units. The google map
#' data is cached to provide fast plotting capabilities.
#'
#' @param x `RapData`, `RapUnsolved`, `RapSolved` object.
#'
#' @param basemap `character` type of base map to display. Valid names are
#' `"roadmap"`, `"mobile"`, `"satellite"`, `"terrain"`,
#' `"hybrid"`, `"mapmaker-roadmap"`, `"mapmaker-hybrid"`.
#'
#' @param grayscale `logical` should base map be gray scale?
#'
#' @param force.reset `logical` ignore data in cache? Setting this as
#' ignore will make function slower but may avoid bugs in cache system.
#'
#' @return `list` with google map data.
#'
#' @keywords internal
#'
#' @seealso `RgoogleMaps::GetMap.bbox()`, [plot()].
basemap <- function(x, basemap = "hybrid", grayscale = FALSE,
force.reset = FALSE) UseMethod("basemap")
#' Test if hash is cached in a Rap object
#'
#' Tests if hash is cached in Rap object.
#'
#' @param x `RapData` or `RapResults` object
#'
#' @param name `character` hash.
#'
#' @note caches are implemented using environments, the hash is used as the
#' name of the object in the environment.
#' @return `logical` Is it cached?
#'
#' @keywords internal
#'
#' @name is.cached
#'
#' @aliases is.cached,RapData,character-method is.cached,RapResults,character-method
methods::setGeneric("is.cached",
function(x, name) methods::standardGeneric("is.cached"))
#' Get and set cache methods
#'
#' Getter and setter methods for caches in RapData and RapResults object.
#'
#' @param x `RapData` or `RapResults` object
#'
#' @param name `character` hash.
#'
#' @param y if `ANY` this object gets cached with name, else if
#' `missing` the object hashed at name gets returned.
#'
#' @note caches are implemented using environments, the hash is used as the
#' name of the object in the environment.
#'
#' @return `ANY` or `NULL` depends on `y` argument.
#'
#' @keywords internal
#'
#' @name cache
#'
#' @aliases cache,RapData,character,ANY-method cache,RapData,character,missing-method cache,RapResults,character,ANY-method cache,RapResults,character,missing-method
methods::setGeneric("cache",
function(x, name, y) methods::standardGeneric("cache"))
#' Plot a 1-dimensional attribute space
#'
#' @noRd
spacePlot.1d <- function(pu, dp, pu.color.palette, main) {
# create X2 vars
pu$X2 <- 0
dp$X2 <- 0
# create colors
if (length(pu.color.palette) == 1)
pu.color.palette <- brewerCols(seq(0, 1, 0.25), pu.color.palette, 4)
# make plot
ggplot2::ggplot() +
ggplot2::geom_point(
ggplot2::aes(
x = .data[["X1"]], y = .data[["X2"]], alpha = .data[["weights"]]
),
data = dp, color = "darkblue", size = 5,
position = ggplot2::position_jitter(width = 0, height = 5)
) +
ggplot2::scale_alpha_continuous(name = "Demand point weight") +
ggplot2::geom_point(
ggplot2::aes(
x = .data[["X1"]], y = .data[["X2"]],
color = .data[["status"]], size = .data[["status"]]
),
data = pu,
position = ggplot2::position_jitter(width = 0, height = 5)) +
ggplot2::scale_color_manual(
name = "Planning unit status",
values = c(
"Locked Out" = pu.color.palette[4],
"Not Selected" = pu.color.palette[1],
"Selected" = pu.color.palette[2],
"Locked In" = pu.color.palette[3]
)
) +
ggplot2::scale_size_manual(
values = c(
"Locked Out" = 2, "Not Selected" = 2,
"Selected" = 4.5, "Locked In" = 4.5
),
guide = "none"
) +
ggplot2::guides(
alpha = ggplot2::guide_legend(order = 1),
color = ggplot2::guide_legend(order = 2)
) +
ggplot2::theme_classic() +
ggplot2::coord_equal() +
ggplot2::theme(
legend.position = "right",
axis.title.y = ggplot2::element_blank(),
axis.ticks.y = ggplot2::element_blank(),
axis.text.y = ggplot2::element_blank(),
axis.line.y = ggplot2::element_line(),
axis.line.x = ggplot2::element_line()
) +
ggplot2::ggtitle(main) +
ggplot2::xlab("Dimension 1") +
ggplot2::ylab("")
}
#' Plot a 2-d attribute space
#'
#' @noRd
spacePlot.2d <- function(pu, dp, pu.color.palette, main) {
# create colors
if (length(pu.color.palette) == 1)
pu.color.palette <- brewerCols(seq(0, 1, 0.25), pu.color.palette, 4)
# make plot
ggplot2::ggplot() +
ggplot2::geom_point(
ggplot2::aes(
x = .data[["X1"]], y = .data[["X2"]], alpha = .data[["weights"]]
),
data = dp, color = "darkblue", size = 5
) +
ggplot2::scale_alpha_continuous(name = "Demand point weight") +
ggplot2::geom_point(
ggplot2::aes(
x = .data[["X1"]], y = .data[["X2"]],
color = .data[["status"]], size = .data[["status"]]
),
data = pu
) +
ggplot2::scale_color_manual(
name = "Planning unit status",
values = c(
"Locked Out" = pu.color.palette[4],
"Not Selected" = pu.color.palette[1],
"Selected" = pu.color.palette[2],
"Locked In" = pu.color.palette[3]
)
) +
ggplot2::scale_size_manual(
values = c(
"Locked Out" = 2, "Not Selected" = 2,
"Selected" = 4.5, "Locked In" = 4.5
),
guide = "none"
) +
ggplot2::guides(
alpha = ggplot2::guide_legend(order = 1),
color = ggplot2::guide_legend(order = 2)
) +
ggplot2::theme_classic() +
ggplot2::coord_equal() +
ggplot2::theme(
legend.position = "right",
axis.line.y = ggplot2::element_line(),
axis.line.x = ggplot2::element_line()
) +
ggplot2::ggtitle(main) +
ggplot2::xlab("Dimension 1") +
ggplot2::ylab("Dimension 2")
}
#' Plot a 3-d attribute space
#'
#' @noRd
spacePlot.3d <- function(pu, dp, pu.color.palette, main) {
# check if rgl is installed
assertthat::assert_that(
requireNamespace("rgl", quietly = TRUE),
msg = "please install the \"rgl\" package"
)
# create frame
rgl::open3d()
# add pu points
if (length(pu.color.palette) == 1)
pu.color.palette <- brewerCols(seq(0, 1, 0.25), pu.color.palette, 4)
pu.cols <- character(nrow(pu))
pu.cols[which(pu$status == "Not Selected")] <- pu.color.palette[1]
pu.cols[which(pu$status == "Selected")] <- pu.color.palette[2]
pu.cols[which(pu$status == "Locked In")] <- pu.color.palette[3]
pu.cols[which(pu$status == "Locked Out")] <- pu.color.palette[4]
rgl::points3d(as.matrix(pu[, seq_len(3)]), col = pu.cols)
# add dp points
dp.cols <- ggplot2::alpha(
rep("darkblue", nrow(dp)),
affineTrans(dp$weights, min(dp$weights), max(dp$weights), 0.1, 1)
)
rgl::points3d(as.matrix(dp[, seq_len(3)]), col = dp.cols)
rgl::title3d(main)
}
#' General argument parsing function
#'
#' @noRd
parseArgs <- function(fn, object = NULL, skip = -1, ...) {
if (!is.null(object))
fn <- paste0(fn, ".", class(object))
ellipses.args <- list(...)
fn.args <- names(formals(fn))
if (!is.null(skip))
fn.args <- fn.args[skip]
ellipses.args[intersect(names(ellipses.args), fn.args)]
}
#' Alternative argument parsing function
#'
#' @noRd
parseArgs2 <- function(args, ...) {
ellipses.args <- list(...)
ellipses.args[intersect(names(ellipses.args), args)]
}
#' Create an empty PolySet object
#' @noRd
emptyPolySet <- function() {
structure(
list(
PID = integer(0),
SID = integer(0),
POS = integer(0),
X = numeric(0),
Y = numeric(0)
),
.Names = c("PID", "SID", "POS", "X", "Y"),
row.names = integer(0),
class = c("PolySet", "data.frame")
)
}
#' Calculate distances between points using URAP
#'
#' @noRd
urap.squared.distance <- function(x, y, y.weights = rep(1, nrow(y))) {
assertthat::assert_that(
inherits(x, "matrix"),
inherits(y, "matrix"),
is.numeric(y.weights),
nrow(y) == length(y.weights),
all(is.finite(c(x))),
all(is.finite(c(y))),
all(is.finite(c(y.weights))), all(y.weights >= 0)
)
rcpp_squared_distance(x, y, y.weights)
}
#' Calculate distances between points using RRAP
#'
#' @noRd
rrap.squared.distance <- function(pu.coordinates, pu.probabilities,
dp.coordinates, dp.weights, failure.distance,
maximum.r.level =
as.integer(length(pu.probabilities))) {
# data integreity checks
assertthat::assert_that(
inherits(pu.coordinates, "matrix"),
inherits(dp.coordinates, "matrix"),
inherits(pu.probabilities, "numeric"),
is.numeric(dp.weights),
assertthat::is.scalar(failure.distance),
assertthat::is.count(maximum.r.level),
nrow(pu.coordinates) == length(pu.probabilities),
nrow(dp.coordinates) == length(dp.weights),
ncol(dp.coordinates) == ncol(pu.coordinates),
all(is.finite(c(dp.weights))),
all(is.finite(c(pu.probabilities))),
all(is.finite(c(pu.coordinates))),
all(is.finite(c(dp.coordinates))),
all(is.finite(c(failure.distance))),
all(is.finite(c(maximum.r.level))),
maximum.r.level <= nrow(pu.coordinates),
failure.distance >= 0,
nrow(pu.coordinates) >= 1,
nrow(dp.coordinates) >= 1
)
# main processing
rcpp_rrap_squared_distance(
pu.coordinates,
pu.probabilities,
dp.coordinates,
dp.weights,
failure.distance,
maximum.r.level
)
}
#' Dump object from model cache
#'
#' @noRd
dump_object <- function(x, mode = c("numeric", "integer", "character")) {
assertthat::assert_that(inherits(x, "externalptr"))
mode <- match.arg(mode)
if (mode == "numeric") {
return(rcpp_dump_numeric_object(x))
}
if (mode == "integer") {
return(rcpp_dump_integer_object(x))
}
if (mode == "character") {
return(rcpp_dump_character_object(x))
}
}
#' Simulate a Gaussian random field
#'
#' @param n `integer` Number of simulations to generate.
#'
#' @param coords `matrix` Matrix containing coordinates for simulating data.
#'
#' @param mu `numeric` Parameter for Gaussian simulations.
#'
#' @param scale `numeric` Parameter for Gaussian simulations.
#'
#' @details
#' This function is largely inspired by: https://rpubs.com/jguelat/autocorr.
#'
#' @return `matrix` object with simulations.
#'
#' @noRd
simulate_gaussian_random_field <- function(n, coords, mu, scale) {
# assert valid arguments
assertthat::assert_that(
assertthat::is.count(n),
assertthat::noNA(n),
is.matrix(coords),
assertthat::noNA(c(coords)),
nrow(coords) >= 1,
ncol(coords) == 2,
assertthat::is.number(mu),
assertthat::noNA(mu),
assertthat::is.number(scale),
assertthat::noNA(scale)
)
# main processing
mu <- rep(0, nrow(coords))
p <- nrow(coords)
chol_d <- chol(exp(-scale * as.matrix(stats::dist(coords))))
out <- t(
matrix(stats::rnorm(n * p), ncol = p) %*%
chol_d + rep(mu, rep(n, p))
)
# ensure matrix output
if (!is.matrix(out)) {
out <- matrix(out, ncol = 1)
}
# return result
out
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.