diff --git a/.github/workflows/Render-README.yaml b/.github/workflows/Render-README.yaml index 29015cfc..c32ca370 100644 --- a/.github/workflows/Render-README.yaml +++ b/.github/workflows/Render-README.yaml @@ -38,4 +38,4 @@ jobs: git config --local user.email "$GITHUB_ACTOR@users.noreply.github.com" git add README.md man/figures/README-* git commit -m "Re-build README.md" || echo "No changes to commit" - git push origin || echo "No changes to commit" + git push -f origin || echo "No changes to commit" diff --git a/DESCRIPTION b/DESCRIPTION index 1acbc6c2..08f0aa92 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Type: Package Package: shar Title: Species-Habitat Associations -Version: 2.0.4 +Version: 2.1 Authors@R: c(person("Maximilian H.K.", "Hesselbarth", email = "mhk.hesselbarth@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-1125-9918")), person("Marco", "Sciaini", email = "marco.sciaini@posteo.net", @@ -30,8 +30,8 @@ Imports: grDevices, methods, spatstat.explore, - spatstat.model, spatstat.geom, + spatstat.model, spatstat.random, stats, terra, diff --git a/NAMESPACE b/NAMESPACE index 1f720795..2fbc51de 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,7 +8,6 @@ S3method(print,rd_pat) S3method(print,rd_ras) export(calculate_energy) export(classify_habitats) -export(estimate_pcf_fast) export(fit_point_process) export(list_to_randomized) export(pack_randomized) diff --git a/NEWS.md b/NEWS.md index 5426d263..8e35a6b4 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,12 @@ +# shar 2.1 +* Improvements + * Remove `comp_fast` argument + * Speed improvements of computation + * General updates to code structure +* Bugfixes + * Removed `n_points` and `window` argument from reconstruction due to methodological issues + * Bugfix related to wrap/unwrap raster and printing + # shar 2.0.4 * Improvements * Remove `Sys.sleep` for verbose reconstruction diff --git a/R/calc_gest.R b/R/calc_gest.R new file mode 100644 index 00000000..065fa852 --- /dev/null +++ b/R/calc_gest.R @@ -0,0 +1,32 @@ +#' calc_gest +#' +#' @description Calculate Gest +#' +#' @param dist matrix with distance pairs. +#' @param r vector with distances r. +#' @param n_points numeric with number of points +#' +#' @details +#' Calculates Gest based on distances created with \code{get_dist_pairs}. +#' +#' @seealso +#' \code{\link{get_dist_pairs}} +#' +#' @return data.frame +#' +#' @aliases calc_gest +#' @rdname calc_gest +#' +#' @keywords internal +calc_gest <- function(dist, r, n_points){ + + mat <- matrix(nrow = n_points, ncol = n_points, data = Inf) + mat[dist[, 1:2]] <- dist[, 3] + + distances_min <- apply(X = mat, MARGIN = 2, FUN = min, na.rm = TRUE) + + hist_min <- graphics::hist(distances_min, breaks = r, plot = FALSE) + + data.frame(r = hist_min$mids, edf = cumsum(hist_min$counts) / n_points) + +} diff --git a/R/calculate_energy.R b/R/calculate_energy.R index bc313fe6..20b18d50 100644 --- a/R/calculate_energy.R +++ b/R/calculate_energy.R @@ -6,18 +6,13 @@ #' @param weights Vector with weights used to calculate energy. #' The first number refers to Gest(r), the second number to pcf(r). #' @param return_mean Logical if the mean energy is returned. -#' @param comp_fast Integer with threshold at which summary functions are estimated -#' in a computational fast way. #' @param verbose Logical if progress report is printed. #' #' @details #' The function calculates the mean energy (or deviation) between the observed #' pattern and all reconstructed patterns (for more information see Tscheschel & #' Stoyan (2006) or Wiegand & Moloney (2014)). The pair correlation function and the -#' nearest neighbour distance function are used to describe the patterns. For large -#' patterns \code{comp_fast = TRUE} decreases the computational demand, because no edge -#' correction is used and the pair correlation function is estimated based on Ripley's -#' K-function. For more information see \code{\link{estimate_pcf_fast}}. +#' nearest neighbour distance function are used to describe the patterns. #' #' @seealso #' \code{\link{plot_energy}} \cr @@ -54,9 +49,8 @@ #' #' @export calculate_energy <- function(pattern, - weights = c(0.5, 0.5), + weights = c(1, 1), return_mean = FALSE, - comp_fast = 1000, verbose = TRUE){ # check if class is correct @@ -83,7 +77,7 @@ calculate_energy <- function(pattern, # calculate r sequence r <- seq(from = 0, to = spatstat.explore::rmax.rule(W = pattern_observed$window, - lambda = spatstat.geom::intensity.ppp(pattern_observed)), + lambda = spatstat.geom::intensity.ppp(pattern_observed)), length.out = 250) if (inherits(x = pattern, what = "rd_pat")) { @@ -96,71 +90,22 @@ calculate_energy <- function(pattern, } else { - # check if weights make sense - if (sum(weights) > 1 || sum(weights) == 0) { - - stop("The sum of 'weights' must be 0 < sum(weights) <= 1.", call. = FALSE) - - } - - # check if number of points exceed comp_fast limit - if (pattern_observed$n > comp_fast) { - - comp_fast <- TRUE - - } else { - - comp_fast <- FALSE - - } - # calculate summary functions for observed pattern - if (comp_fast) { - - gest_observed <- spatstat.explore::Gest(X = pattern_observed, correction = "none", - r = r) - - pcf_observed <- estimate_pcf_fast(pattern = pattern_observed, - correction = "none", method = "c", - spar = 0.5, r = r) + gest_observed <- spatstat.explore::Gest(X = pattern_observed, + correction = "none", r = r) - } else { - - gest_observed <- spatstat.explore::Gest(X = pattern_observed, - correction = "han", r = r) - - pcf_observed <- spatstat.explore::pcf(X = pattern_observed, - correction = "best", divisor = "d", r = r) - - } + pcf_observed <- spatstat.explore::pcf(X = pattern_observed, + correction = "none", divisor = "d", r = r) # loop through all reconstructed patterns result <- vapply(seq_along(pattern_randomized), function(x) { - # fast computation of summary stats - if (comp_fast) { - - gest_reconstruction <- spatstat.explore::Gest(X = pattern_randomized[[x]], - correction = "none", - r = r) - - pcf_reconstruction <- estimate_pcf_fast(pattern = pattern_randomized[[x]], - correction = "none", method = "c", - spar = 0.5, r = r) - - # normal computation of summary stats - } else { - - gest_reconstruction <- spatstat.explore::Gest(X = pattern_randomized[[x]], - correction = "han", - r = r) + gest_reconstruction <- spatstat.explore::Gest(X = pattern_randomized[[x]], + correction = "none", r = r) pcf_reconstruction <- spatstat.explore::pcf(X = pattern_randomized[[x]], - correction = "best", - divisor = "d", - r = r) - - } + correction = "none", divisor = "d", + r = r) # difference between observed and reconstructed pattern energy <- (mean(abs(gest_observed[[3]] - gest_reconstruction[[3]]), na.rm = TRUE) * weights[[1]]) + diff --git a/R/data.R b/R/data.R index 3f6ee34d..18470c86 100644 --- a/R/data.R +++ b/R/data.R @@ -34,31 +34,3 @@ #' #' @format A spatstat ppp object. "species_b" - -#' Gamma test -#' -#' Randomized data for species b using the gamma test. -#' -#' @format rd_pat object. -"gamma_test" - -#' Reconstruction -#' -#' Randomized data for species b using pattern reconstruction. -#' -#' @format rd_pat object. -"reconstruction" - -#' Torus trans -#' -#' Torus translation of the classified \code{landscape}. -#' -#' @format rd_ras object. -"torus_trans" - -#' Random walk -#' -#' Randomization of the \code{landscape} using the habitat randomization algorithm. -#' -#' @format rd_ras object. -"random_walk" diff --git a/R/estimate_pcf_fast.R b/R/estimate_pcf_fast.R deleted file mode 100644 index 7b3750a2..00000000 --- a/R/estimate_pcf_fast.R +++ /dev/null @@ -1,46 +0,0 @@ -#' estimate_pcf_fast -#' -#' @description Fast estimation of the pair correlation function -#' -#' @param pattern ppp object with point pattern. -#' @param ... Arguments passed down to \code{link{Kest}} or \code{\link{pcf.fv}}. -#' -#' @details -#' The functions estimates the pair correlation functions based on an estimation -#' of Ripley's K-function. This makes it computationally faster than estimating the -#' pair correlation function directly. It is a wrapper around \code{\link{Kest}} and -#' \code{\link{pcf.fv}}. -#' -#' @seealso -#' \code{\link{Kest}} \cr -#' \code{\link{pcf.fv}} -#' -#' @return fv.object -#' -#' @examples -#' pcf_species_b <- estimate_pcf_fast(species_a) -#' -#' @aliases estimate_pcf_fast -#' @rdname estimate_pcf_fast -#' -#' @references -#' Chiu, S.N., Stoyan, D., Kendall, W.S., Mecke, J., 2013. Stochastic geometry and -#' its applications, 3rd ed, Wiley Series in Probability and Statistics. -#' John Wiley & Sons Inc, Chichester, UK. ISBN 978-0-470-66481-0 -#' -#' Ripley, B.D., 1977. Modelling spatial patterns. Journal of the Royal Statistical -#' Society. Series B (Methodological) 39, 172–192. -#' -#' -#' Stoyan, D., Stoyan, H., 1994. Fractals, random shapes and point fields. -#' John Wiley & Sons, Chichester. ISBN 978-0-471-93757-9 -#' -#' @export -estimate_pcf_fast <- function(pattern, ...){ - - k_fun <- suppressMessages(spatstat.explore::Kest(X = pattern, ...)) # estimate K-fct - - result <- spatstat.explore::pcf.fv(X = k_fun, ...) # estimate pcf from K-fct - - return(result) -} diff --git a/R/get_dist_pairs.R b/R/get_dist_pairs.R new file mode 100644 index 00000000..a6266772 --- /dev/null +++ b/R/get_dist_pairs.R @@ -0,0 +1,26 @@ +#' get_dist_pairs +#' +#' @description Distance between points +#' +#' @param X ppp object +#' @param rmax Numeric with maximum distance +#' +#' @details +#' Returns matrix with point pairs and distances between them. +#' +#' @seealso +#' \code{\link{pcf.ppp}} +#' +#' @return matrix +#' +#' @aliases get_dist_pairs +#' @rdname get_dist_pairs +#' +#' @keywords internal +get_dist_pairs <- function(X, rmax){ + + dist_observed <- spatstat.geom::closepairs(X = X, rmax = rmax, what = "ijd", twice = TRUE) + + cbind(dist_observed$i, dist_observed$j, dist_observed$d) + +} diff --git a/R/pack_randomized.R b/R/pack_randomized.R index 618a4dde..0b3daed0 100644 --- a/R/pack_randomized.R +++ b/R/pack_randomized.R @@ -27,8 +27,9 @@ #' @export pack_randomized <- function(raster) { - # wrap observerd raster - raster$observed <- terra::wrap(raster$observed) + # check if observed is present + # wrap observed raster + if (inherits(x = raster$observed, what = "SpatRaster")) raster$observed <- terra::wrap(raster$observed) # wrap all randomized raster raster$randomized <- lapply(X = raster$randomized, FUN = terra::wrap) diff --git a/R/plot.rd_mar.R b/R/plot.rd_mar.R index 7e1605b0..b4a258f9 100644 --- a/R/plot.rd_mar.R +++ b/R/plot.rd_mar.R @@ -8,8 +8,6 @@ #' @param n Integer with number or vector of ids of randomized pattern to plot. #' See Details section for more information. #' @param probs Vector with quantiles of randomized data used for envelope construction. -#' @param comp_fast Integer with threshold at which summary functions are estimated -#' in a computational fast way. #' @param ask Logical if the user is asked to press before second summary function #' is plotted (only used if \code{what = "sf"}). #' @param verbose Logical if progress report is printed. @@ -18,9 +16,6 @@ #' @details #' The function plots the pair correlation function and the nearest neighbour function of #' the observed pattern and the reconstructed patterns (as "simulation envelopes"). -#' For large patterns \code{comp_fast = TRUE} decreases the computational demand because no edge -#' correction is used and the pair correlation function is estimated based on Ripley's -#' K-function. For more information see \code{\link{estimate_pcf_fast}}. #' #' It is also possible to plot n randomized patterns and the observed pattern #' using \code{what = "pp"}. If \code{n} is a single number, \code{n} randomized @@ -48,7 +43,7 @@ #' @rdname plot.rd_mar #' #' @export -plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_fast = 1000, +plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), ask = TRUE, verbose = TRUE, ...) { # check if class is correct @@ -70,17 +65,6 @@ plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_ if (what == "sf") { - # check if number of points exceed comp_fast limit - if (x$observed$n > comp_fast) { - - comp_fast <- TRUE - - } else { - - comp_fast <- FALSE - - } - name_unit <- spatstat.geom::unitname(x$observed)[[1]] # unit name for labels # calculate r @@ -135,8 +119,7 @@ plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_ # specify quantums g(r) col_kmmr <- ifelse(test = result_observed[, 3] < result_randomized[, 2] | result_observed[, 3] > result_randomized[, 3], - yes = "#1f78b4", - no = "#b2df8a") + yes = "#1f78b4", no = "#b2df8a") # plot results graphics::plot(NULL, xlim = range(r), ylim = yrange, @@ -195,8 +178,11 @@ plot.rd_mar <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_ # convert to dataframe current_pattern <- as.data.frame(subset_pattern[[i]]) + current_pattern$marks <- ((current_pattern$marks - min(current_pattern$marks)) / + (max(current_pattern$marks) - min(current_pattern$marks)) * 1) + 0.25 + # plot points - graphics::plot(x = current_pattern$x, y = current_pattern$y, + graphics::plot(x = current_pattern$x, y = current_pattern$y, cex = current_pattern$marks, type = "p", asp = 1, xlim = x_range, ylim = y_range, axes = FALSE, main = names_pattern[[i]], xlab = "", ylab = "") diff --git a/R/plot.rd_pat.R b/R/plot.rd_pat.R index a0a8ba09..9edba871 100644 --- a/R/plot.rd_pat.R +++ b/R/plot.rd_pat.R @@ -8,8 +8,6 @@ #' @param n Integer with number or vector of ids of randomized pattern to plot. #' See Details section for more information. #' @param probs Vector with quantiles of randomized data used for envelope construction. -#' @param comp_fast Integer with threshold at which summary functions are estimated -#' in a computational fast way. #' @param ask Logical if the user is asked to press before second summary function #' is plotted (only used if \code{what = "sf"}). #' @param verbose Logical if progress report is printed. @@ -18,9 +16,6 @@ #' @details #' The function plots the pair correlation function and the nearest neighbour function of #' the observed pattern and the reconstructed patterns (as "simulation envelopes"). -#' For large patterns \code{comp_fast = TRUE} decreases the computational demand because no edge -#' correction is used and the pair correlation function is estimated based on Ripley's -#' K-function. For more information see \code{\link{estimate_pcf_fast}}. #' #' It is also possible to plot n randomized patterns and the observed pattern #' using \code{what = "pp"}. If \code{n} is a single number, \code{n} randomized @@ -47,7 +42,7 @@ #' @rdname plot.rd_pat #' #' @export -plot.rd_pat <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_fast = 1000, +plot.rd_pat <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), ask = TRUE, verbose = TRUE, ...) { # check if class is correct @@ -69,17 +64,6 @@ plot.rd_pat <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_ if (what == "sf") { - # check if number of points exceed comp_fast limit - if (x$observed$n > comp_fast) { - - comp_fast <- TRUE - - } else { - - comp_fast <- FALSE - - } - name_unit <- spatstat.geom::unitname(x$observed)[[1]] # unit name for labels # calculate r @@ -96,21 +80,10 @@ plot.rd_pat <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), comp_ result <- lapply(seq_along(pattern), function(x) { # calculate summary functions - if (comp_fast) { - - gest_result <- spatstat.explore::Gest(pattern[[x]], correction = "none", r = r) + gest_result <- spatstat.explore::Gest(pattern[[x]], correction = "none", r = r) - pcf_result <- estimate_pcf_fast(pattern[[x]], correction = "none", - method = "c", spar = 0.5, r = r) - - } else { - - gest_result <- spatstat.explore::Gest(pattern[[x]], correction = "han", r = r) - - pcf_result <- spatstat.explore::pcf(pattern[[x]], divisor = "d", - correction = "best", r = r) - - } + pcf_result <- spatstat.explore::pcf(pattern[[x]], divisor = "d", + correction = "none", r = r) gest_df <- as.data.frame(gest_result) # conver to df diff --git a/R/plot.rd_ras.R b/R/plot.rd_ras.R index 0f02e6b7..63a0f7d7 100644 --- a/R/plot.rd_ras.R +++ b/R/plot.rd_ras.R @@ -42,7 +42,7 @@ plot.rd_ras <- function(x, n = NULL, col, verbose = TRUE, nrow, ncol, ...) { } # check if observed is present - if (!methods::is(x$observed, "SpatRaster")) { + if (!inherits(x = x$observed, what = "SpatRaster")) { stop("Input must include 'observed' raster.", call. = FALSE) diff --git a/R/print.rd_ras.R b/R/print.rd_ras.R index ef68054c..147cf442 100644 --- a/R/print.rd_ras.R +++ b/R/print.rd_ras.R @@ -29,8 +29,14 @@ #' @export print.rd_ras <- function(x, ...) { + if (inherits(x = x$randomized[[1]], what = "PackedSpatRaster")) { + + cat("Please use 'unpack_randomized()' to unwrap object.") + + } else { + # check if observed raster is included - if (!methods::is(x$observed, "SpatRaster")) { + if (!inherits(x = x$observed, what = "SpatRaster")) { number_raster_obs <- 0 @@ -58,4 +64,5 @@ print.rd_ras <- function(x, ...) { "Method: ", x$method, "\n", "Observed pattern: ", includes_observed, "\n", "Extent: ", extent_window, " (xmin, xmax, ymin, ymax) \n"), ...) + } } diff --git a/R/reconstruct_algorithm.R b/R/reconstruct_algorithm.R new file mode 100644 index 00000000..d17a37c1 --- /dev/null +++ b/R/reconstruct_algorithm.R @@ -0,0 +1,323 @@ +#' reconstruct_algorithm +#' +#' @description Pattern reconstruction (internal) +#' +#' @param pattern ppp object with pattern. +#' @param n_random Integer with number of randomizations. +#' @param e_threshold Double with minimum energy to stop reconstruction. +#' @param max_runs Integer with maximum number of iterations if \code{e_threshold} +#' is not reached. +#' @param no_change Integer with number of iterations at which the reconstruction will +#' stop if the energy does not decrease. +#' @param annealing Double with probability to keep relocated point even if energy +#' did not decrease. +#' @param weights Vector with weights used to calculate energy. +#' The first number refers to Gest(r), the second number to pcf(r). +#' @param r_length Integer with number of intervals from \code{r = 0} to \code{r = rmax} for which +#' the summary functions are evaluated. +#' @param r_max Double with maximum distance used during calculation of summary functions. If \code{NULL}, +#' will be estimated from data. +#' @param stoyan Coefficient for Stoyan's bandwidth selection rule. +#' @param verbose Logical if progress report is printed. +#' @param plot Logical if pcf(r) function is plotted and updated during optimization. +#' +#' @return list +#' +#' @aliases reconstruct_algorithm +#' @rdname reconstruct_algorithm +#' +#' @keywords internal +reconstruct_algorithm <- function(pattern, + method, + n_random, + e_threshold, + max_runs, + no_change, + annealing, + weights, + r_length, + r_max, + stoyan, + verbose, + plot){ + + # check if n_random is >= 1 + if (n_random < 1) { + + stop("n_random must be >= 1.", call. = FALSE) + + } + + # unmark pattern + if (spatstat.geom::is.marked(pattern)) { + + pattern <- spatstat.geom::unmark(pattern) + + if (verbose) { + warning("Unmarked provided input pattern. For marked pattern, see reconstruct_pattern_marks().", + call. = FALSE) + + } + } + + # grab window and number of points + n_points <- pattern$n + window <- pattern$window + + # calculate intensity + lambda <- n_points / spatstat.geom::area(window) + lambda2area <- (n_points * (n_points - 1)) / spatstat.geom::area(window) + + # calculate bandwidth using + h <- stoyan / sqrt(lambda) + bw <- h / sqrt(5) + + # calculate r from data + if (is.null(r_max)) { + + r <- seq(from = 0, to = spatstat.explore::rmax.rule(W = window, lambda = lambda), + length.out = r_length) + + r_max <- max(r) + + # use provided r_max + } else { + + r <- seq(from = 0, to = r_max, length.out = r_length) + + } + + # create list with arguments fpr pcf + denargs <- list(kernel = "epanechnikov", bw = bw, n = r_length, from = 0, to = r_max) + + # set names of randomization randomized_1 ... randomized_n + names_randomization <- paste0("randomized_", seq_len(n_random)) + + # create empty lists for results + energy_list <- vector("list", length = n_random) + iterations_list <- vector("list", length = n_random) + stop_criterion_list <- as.list(rep("max_runs", times = n_random)) + result_list <- vector("list", length = n_random) + + # set names + names(energy_list) <- names_randomization + names(iterations_list) <- names_randomization + names(stop_criterion_list) <- names_randomization + names(result_list) <- names_randomization + + # calculate summary functions + dist_observed <- get_dist_pairs(X = pattern, rmax = r_max) + + gest_observed <- calc_gest(dist = dist_observed, r = r, n_points = pattern$n) + + pcf_observed <- spatstat.explore::sewpcf(d = dist_observed[, 3], w = 1, denargs = denargs, + lambda2area = lambda2area, divisor = "d") + + # create n_random recondstructed patterns + for (i in seq_len(n_random)) { + + # simulated starting pattern for reconstruction + if (method == "homo") { + + # create Poisson simulation data + simulated <- spatstat.random::runifpoint(n = n_points, nsim = 1, drop = TRUE, + win = window, warn = FALSE) + + } else if (method == "hetero") { + + # estimate lambda(x,y) + lambda_xy <- spatstat.explore::density.ppp(pattern) + + # create starting pattern + simulated <- spatstat.random::rpoint(n = n_points, f = lambda_xy, + win = pattern$window) + + + } else if (method == "cluster") { + + # start with fitted Thomas cluster model + cluster_ppm <- spatstat.model::kppm.ppp(pattern, cluster = "Thomas", statistic = "pcf", + statargs = list(divisor = "d", correction = "best"), + improve.type = "none") + + # simulate clustered pattern + simulated <- spatstat.model::simulate.kppm(cluster_ppm, nsim = 1, drop = TRUE, + window = window, verbose = FALSE) + + # remove points because more points in simulated + if (n_points <= simulated$n) { + + # difference between patterns + difference <- simulated$n - n_points + # id of points to remove + remove_points <- sample(x = seq_len(simulated$n), size = difference) + # remove points + simulated <- simulated[-remove_points] + + # add points because less points in simulated + } else { + + # difference between patterns + difference <- n_points - simulated$n + # create missing points + missing_points <- spatstat.random::runifpoint(n = difference, nsim = 1, drop = TRUE, + win = pattern$window, warn = FALSE) + # add missing points to simulated + simulated <- spatstat.geom::superimpose.ppp(simulated, missing_points, + W = pattern$window, check = FALSE) + + } + } + + # energy before reconstruction + energy <- Inf + + # counter if energy changed + energy_counter <- 0 + + # df for energy + energy_df <- data.frame(i = seq(from = 1, to = max_runs, by = 1), energy = NA) + + # create random number for annealing prob + if (annealing != 0) { + + random_annealing <- stats::runif(n = max_runs, min = 0, max = 1) + + } else { + + random_annealing <- rep(0, max_runs) + + } + + # random ids of pattern + rp_id <- sample(x = seq_len(simulated$n), size = max_runs, replace = TRUE) + + # create random new points + # MH: This could use same method as simulated? + rp_coords <- spatstat.random::runifpoint(n = max_runs, nsim = 1, drop = TRUE, + win = simulated$window, warn = FALSE) + + # # MH: Use this distance to check what changed + # dist_simulated <- get_dist_pairs(X = simulated, rmax = r_max) + + # pattern reconstruction algorithm (optimization of energy) - not longer than max_runs + for (j in seq_len(max_runs)) { + + # data for relocation + relocated <- simulated + + # get current point id + id_current <- rp_id[[j]] + + # relocate point + relocated$x[[id_current]] <- rp_coords$x[[j]] + + relocated$y[[id_current]] <- rp_coords$y[[j]] + + # calculate summary functions relocation + + dist_relocated <- get_dist_pairs(X = relocated, rmax = r_max) + + gest_relocated <- calc_gest(dist = dist_relocated, r = r, n_points = n_points) + + pcf_relocated <- spatstat.explore::sewpcf(d = dist_relocated[, 3], w = 1, denargs = denargs, + lambda2area = lambda2area, divisor = "d") + + # energy after relocation + energy_relocated <- (mean(abs(gest_observed[, 2] - gest_relocated[, 2]), na.rm = TRUE) * weights[[1]]) + + (mean(abs(pcf_observed[, 2] - pcf_relocated[, 2]), na.rm = TRUE) * weights[[2]]) + + # lower energy after relocation + if (energy_relocated < energy || random_annealing[j] < annealing) { + + # keep relocated pattern + simulated <- relocated + + # keep energy_relocated as energy + energy <- energy_relocated + + # set counter since last change back to 0 + energy_counter <- 0 + + # plot observed vs reconstructed + if (plot) { + + # https://support.rstudio.com/hc/en-us/community/posts/200661917-Graph-does-not-update-until-loop-completion + Sys.sleep(0.01) + + graphics::plot(x = pcf_observed$r, y = pcf_observed$g, + type = "l", col = "black", xlab = "r", ylab = "g(r)") + + graphics::abline(h = 1, lty = 2, col = "grey") + + graphics::lines(x = pcf_relocated$r, y = pcf_relocated$g, col = "red") + + graphics::legend("topright", legend = c("observed", "reconstructed"), + col = c("black", "red"), lty = 1, inset = 0.025) + + } + + # increase counter no change + } else { + + energy_counter <- energy_counter + 1 + + } + + # save energy in data frame + energy_df[j, 2] <- energy + + # print progress + if (verbose) { + + message("\r> Progress: n_random: ", i, "/", n_random, + " || max_runs: ", floor(j / max_runs * 100), "%", + " || energy = ", round(energy, 5), "\t\t", + appendLF = FALSE) + + } + + # exit loop if e threshold or no_change counter max is reached + if (energy <= e_threshold || energy_counter > no_change) { + + # set stop criterion due to energy + stop_criterion_list[[i]] <- ifelse(test = energy <= e_threshold, + yes = "e_threshold", no = "no_change") + + break + + } + } + + # close plotting device + if (plot) { + + grDevices::dev.off() + + } + + # remove NAs if stopped due to energy + if (stop_criterion_list[[i]] %in% c("e_threshold", "no_change")) { + + energy_df <- energy_df[1:j, ] + + } + + # save results in lists + energy_list[[i]] <- energy_df + iterations_list[[i]] <- j + result_list[[i]] <- simulated + + } + + # write result in new line if progress was printed + if (verbose) message("\r") + + # combine to one list + reconstruction <- list(randomized = result_list, observed = pattern, + method = method, energy_df = energy_list, + stop_criterion = stop_criterion_list, + iterations = iterations_list) + + return(reconstruction) +} diff --git a/R/reconstruct_pattern.R b/R/reconstruct_pattern.R index 2c17791d..f475e2bc 100644 --- a/R/reconstruct_pattern.R +++ b/R/reconstruct_pattern.R @@ -13,16 +13,13 @@ #' stop if the energy does not decrease. #' @param annealing Double with probability to keep relocated point even if energy #' did not decrease. -#' @param comp_fast Integer with threshold at which summary functions are estimated -#' in a computational fast way. -#' @param n_points Integer with number of points to be simulated. -#' @param window owin object with window of simulated pattern. #' @param weights Vector with weights used to calculate energy. #' The first number refers to Gest(r), the second number to pcf(r). #' @param r_length Integer with number of intervals from \code{r=0} to \code{r=rmax} for which #' the summary functions are evaluated. #' @param r_max Double with maximum distance used during calculation of summary functions. If \code{NULL}, #' will be estimated from data. +#' @param stoyan Coefficient for Stoyan's bandwidth selection rule. #' @param return_input Logical if the original input data is returned. #' @param simplify Logical if only pattern will be returned if \code{n_random=1} #' and \code{return_input=FALSE}. @@ -37,16 +34,12 @@ #' The pair correlation function and the nearest neighbour distance function are #' used to describe the patterns. #' -#' For large patterns (\code{n > comp_fast}) the pair correlation function can be estimated -#' from Ripley's K-function without edge correction. This decreases the computational -#' time. For more information see \code{\link{estimate_pcf_fast}}. -#' #' The reconstruction can be stopped automatically if for n steps the energy does not #' decrease. The number of steps can be controlled by \code{no_change} and is set to #' \code{no_change = Inf} as default to never stop automatically. #' #' The weights must be 0 < sum(weights) <= 1. To weight both summary functions identical, -#' use \code{weights = c(0.5, 0.5)}. +#' use \code{weights = c(1, 1)}. #' #' \code{spatstat} sets \code{r_length} to 513 by default. However, a lower value decreases #' the computational time, while increasing the "bumpiness" of the summary function. @@ -94,71 +87,65 @@ reconstruct_pattern <- function(pattern, method = "homo", n_random = 1, e_threshold = 0.01, - max_runs = 1000, + max_runs, no_change = Inf, annealing = 0.01, - comp_fast = 1000, - n_points = NULL, - window = NULL, - weights = c(0.5, 0.5), - r_length = 250, + weights = c(1, 1), + r_length = 255, r_max = NULL, + stoyan = 0.15, return_input = TRUE, simplify = FALSE, verbose = TRUE, plot = FALSE) { - if (method == "homo") { - - reconstruction <- reconstruct_pattern_homo(pattern = pattern, n_random = n_random, - e_threshold = e_threshold, max_runs = max_runs, - no_change = no_change, annealing = annealing, - comp_fast = comp_fast, n_points = n_points, - window = window, weights = weights, - r_length = r_length, r_max = r_max, return_input = return_input, - simplify = simplify, verbose = verbose, - plot = plot) + # check if correct method is selected + if (!method %in% c("homo", "hetero", "cluster")) stop("Method must be one of the following: 'homo', 'hetero', or 'cluster'.", + call. = FALSE) + reconstruction <- reconstruct_algorithm(pattern = pattern, method = method, n_random = n_random, + e_threshold = e_threshold, max_runs = max_runs, + no_change = no_change, annealing = annealing, + weights = weights, r_length = r_length, r_max = r_max, + stoyan = stoyan, verbose = verbose, plot = plot) + # set class of result + class(reconstruction) <- "rd_pat" - } else if (method == "cluster") { + # remove input if return_input = FALSE + if (!return_input) { - if (verbose && (!is.null(n_points) || !is.null(window) || !is.null(r_max))) { + # set observed to NA + reconstruction$observed <- "NA" - warning("'n_points', 'window', or 'r_max' are not used for method='cluster'.", call. = FALSE) - - } + # check if output should be simplified + if (simplify) { - reconstruction <- reconstruct_pattern_cluster(pattern = pattern, n_random = n_random, - e_threshold = e_threshold, max_runs = max_runs, - no_change = no_change, annealing = annealing, - comp_fast = comp_fast, weights = weights, - r_length = r_length, r_max = r_max, return_input = return_input, - simplify = simplify, verbose = verbose, - plot = plot) + # not possible if more than one pattern is present + if (n_random > 1 && verbose) { - } else if (method == "hetero") { + warning("'simplify = TRUE' not possible for 'n_random > 1'.", + call. = FALSE) - if (verbose && (!is.null(n_points) || !is.null(window) || !is.null(r_max))) { + # only one random pattern is present that should be returend + } else if (n_random == 1) { - warning("'n_points', 'window', or 'r_max' are not used for method='hetero'.", call. = FALSE) + reconstruction <- reconstruction$randomized[[1]] + } } - reconstruction <- reconstruct_pattern_hetero(pattern = pattern, n_random = n_random, - e_threshold = e_threshold, max_runs = max_runs, - no_change = no_change, annealing = annealing, - comp_fast = comp_fast, weights = weights, - r_length = r_length, r_max = r_max, return_input = return_input, - simplify = simplify, verbose = verbose, - plot = plot) - + # return input if return_input = TRUE } else { - stop("Method must be one of the following: 'homo', 'cluster', 'hetero', or 'marks'.", - call. = FALSE) + # return warning if simply = TRUE because not possible if return_input = TRUE (only verbose = TRUE) + if (simplify && verbose) { + + warning("'simplify = TRUE' not possible for 'return_input = TRUE'.", call. = FALSE) + } } return(reconstruction) + } diff --git a/R/reconstruct_pattern_cluster.R b/R/reconstruct_pattern_cluster.R deleted file mode 100644 index 162c98dc..00000000 --- a/R/reconstruct_pattern_cluster.R +++ /dev/null @@ -1,415 +0,0 @@ -#' reconstruct_pattern_cluster -#' -#' @description Pattern reconstruction for clustered patterns -#' -#' @param pattern ppp object with pattern. -#' @param n_random Integer with number of randomizations. -#' @param e_threshold Double with minimum energy to stop reconstruction. -#' @param max_runs Integer with maximum number of iterations if \code{e_threshold} -#' is not reached. -#' @param no_change Integer with number of iterations at which the reconstruction will -#' stop if the energy does not decrease. -#' @param annealing Double with probability to keep relocated point even if energy -#' did not decrease. -#' @param comp_fast Integer with threshold at which summary functions are estimated -#' in a computational fast way. -#' @param weights Vector with weights used to calculate energy. -#' The first number refers to Gest(r), the second number to pcf(r). -#' @param r_length Integer with number of intervals from \code{r = 0} to \code{r = rmax} for which -#' the summary functions are evaluated. -#' @param r_max Double with maximum distance used during calculation of summary functions. If \code{NULL}, -#' will be estimated from data. -#' @param return_input Logical if the original input data is returned. -#' @param simplify Logical if only pattern will be returned if \code{n_random = 1} -#' and \code{return_input = FALSE}. -#' @param verbose Logical if progress report is printed. -#' @param plot Logical if pcf(r) function is plotted and updated during optimization. -#' -#' @return rd_pat -#' -#' @examples -#' \dontrun{ -#' pattern_recon <- reconstruct_pattern_cluster(species_b, n_random = 19, max_runs = 1000) -#' } -#' -#' @aliases reconstruct_pattern_cluster -#' @rdname reconstruct_pattern_cluster -#' -#' @references -#' Kirkpatrick, S., Gelatt, C.D.Jr., Vecchi, M.P., 1983. Optimization by simulated -#' annealing. Science 220, 671–680. -#' -#' Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random point -#' patterns. Computational Statistics and Data Analysis 51, 859–871. -#' -#' -#' Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern analysis in -#' ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN 978-1-4200-8254-8 -#' -#' @keywords internal -reconstruct_pattern_cluster <- function(pattern, - n_random = 1, - e_threshold = 0.01, - max_runs = 1000, - no_change = Inf, - annealing = 0.01, - comp_fast = 1000, - weights = c(0.5, 0.5), - r_length = 250, - r_max = NULL, - return_input = TRUE, - simplify = FALSE, - verbose = TRUE, - plot = FALSE){ - - # check if n_random is >= 1 - if (n_random < 1) { - - stop("n_random must be >= 1.", call. = FALSE) - - } - - # check if number of points exceed comp_fast limit - if (pattern$n > comp_fast) { - - # Print message that summary functions will be computed fast - if (verbose) { - - message("> Using fast compuation of summary functions.") - - } - - comp_fast <- TRUE - - } else { - - comp_fast <- FALSE - - } - - # set names of randomization randomized_1 ... randomized_n - names_randomization <- paste0("randomized_", seq_len(n_random)) - - # create empty lists for results - energy_list <- vector("list", length = n_random) - iterations_list <- vector("list", length = n_random) - stop_criterion_list <- as.list(rep("max_runs", times = n_random)) - result_list <- vector("list", length = n_random) - - # set names - names(energy_list) <- names_randomization - names(iterations_list) <- names_randomization - names(stop_criterion_list) <- names_randomization - names(result_list) <- names_randomization - - # check if weights make sense - if (sum(weights) > 1 || sum(weights) == 0) { - - stop("The sum of 'weights' must be 0 < sum(weights) <= 1.", call. = FALSE) - - } - - # unmark pattern - if (spatstat.geom::is.marked(pattern)) { - - pattern <- spatstat.geom::unmark(pattern) - - if (verbose) { - - warning("Unmarked provided input pattern. For marked pattern, see reconstruct_pattern_marks().", - call. = FALSE) - - } - } - - # calculate r from data - if (is.null(r_max)) { - - r <- seq(from = 0, to = spatstat.explore::rmax.rule(W = pattern$window, lambda = spatstat.geom::intensity.ppp(pattern)), - length.out = r_length) - - # use provided r_max - } else { - - r <- seq(from = 0, to = r_max, length.out = r_length) - - } - - # start with fitted pattern - # fit Thomas process - fitted_process <- spatstat.model::kppm.ppp(pattern, cluster = "Thomas", - statistic = "pcf", - statargs = list(divisor = "d", - correction = "best"), - improve.type = "none") - - # simulte clustered pattern - simulated <- spatstat.model::simulate.kppm(fitted_process, nsim = 1, drop = TRUE, - window = pattern$window, verbose = FALSE) - - # remove points because more points in simulated - if (pattern$n < simulated$n) { - - # difference between patterns - difference <- simulated$n - pattern$n - # id of points to remove - remove_points <- sample(x = seq_len(simulated$n), size = difference) - # remove points - simulated <- simulated[-remove_points] - - } - - # add points because less points in simulated - if (pattern$n > simulated$n) { - - # difference between patterns - difference <- pattern$n - simulated$n - # create missing points - missing_points <- spatstat.random::runifpoint(n = difference, nsim = 1, drop = TRUE, - win = pattern$window, - warn = FALSE) - # add missing points to simulated - simulated <- spatstat.geom::superimpose.ppp(simulated, missing_points, - W = pattern$window, check = FALSE) - - } - - # fast computation of summary functions - if (comp_fast) { - - gest_observed <- spatstat.explore::Gest(pattern, correction = "none", r = r) - - gest_simulated <- spatstat.explore::Gest(simulated, correction = "none", r = r) - - pcf_observed <- estimate_pcf_fast(pattern, correction = "none", - method = "c", spar = 0.5, r = r) - - pcf_simulated <- estimate_pcf_fast(simulated, correction = "none", - method = "c", spar = 0.5, r = r) - - # normal computation of summary functions - } else { - - gest_observed <- spatstat.explore::Gest(X = pattern, correction = "han", r = r) - - gest_simulated <- spatstat.explore::Gest(X = simulated, correction = "han", r = r) - - pcf_observed <- spatstat.explore::pcf.ppp(X = pattern, correction = "best", - divisor = "d", r = r) - - pcf_simulated <- spatstat.explore::pcf.ppp(X = simulated, correction = "best", - divisor = "d", r = r) - - } - - # energy before reconstruction - energy <- (mean(abs(gest_observed[[3]] - gest_simulated[[3]]), na.rm = TRUE) * weights[[1]]) + - (mean(abs(pcf_observed[[3]] - pcf_simulated[[3]]), na.rm = TRUE) * weights[[2]]) - - # create n_random recondstructed patterns - for (current_pattern in seq_len(n_random)) { - - # current simulated - simulated_current <- simulated - energy_current <- energy - - # counter of iterations - iterations <- 0 - - # counter if energy changed - energy_counter <- 0 - - # df for energy - energy_df <- data.frame(i = seq(from = 1, to = max_runs, by = 1), energy = NA) - - # random ids of pattern - rp_id <- sample(x = seq_len(simulated_current$n), size = max_runs, replace = TRUE) - - # create random new points - rp_coords <- spatstat.random::runifpoint(n = max_runs, nsim = 1, drop = TRUE, - win = simulated_current$window, warn = FALSE) - - # create random number for annealing prob - if (annealing != 0) { - - random_annealing <- stats::runif(n = max_runs, min = 0, max = 1) - - } else { - - random_annealing <- rep(0, max_runs) - - } - - # pattern reconstruction algorithm (optimaztion of energy) - not longer than max_runs - for (i in seq_len(max_runs)) { - - # data for relocation - relocated <- simulated_current - - # get current point id - rp_id_current <- rp_id[[i]] - - # relocate point - relocated$x[[rp_id_current]] <- rp_coords$x[[i]] - - relocated$y[[rp_id_current]] <- rp_coords$y[[i]] - - # calculate summary functions after relocation - if (comp_fast) { - - gest_relocated <- spatstat.explore::Gest(relocated, correction = "none", r = r) - - pcf_relocated <- estimate_pcf_fast(relocated, correction = "none", - method = "c", spar = 0.5, r = r) - - } else { - - gest_relocated <- spatstat.explore::Gest(X = relocated, correction = "han", r = r) - - pcf_relocated <- spatstat.explore::pcf.ppp(X = relocated, correction = "best", - divisor = "d", r = r) - - } - - # energy after relocation - energy_relocated <- (mean(abs(gest_observed[[3]] - gest_relocated[[3]]), na.rm = TRUE) * weights[[1]]) + - (mean(abs(pcf_observed[[3]] - pcf_relocated[[3]]), na.rm = TRUE) * weights[[2]]) - - # lower energy after relocation - if (energy_relocated < energy_current || random_annealing[i] < annealing) { - - # keep relocated pattern - simulated_current <- relocated - - # keep energy_relocated as energy - energy_current <- energy_relocated - - # set counter since last change back to 0 - energy_counter <- 0 - - # plot observed vs reconstructed - if (plot) { - - # https://support.rstudio.com/hc/en-us/community/posts/200661917-Graph-does-not-update-until-loop-completion - Sys.sleep(0.01) - - graphics::plot(x = pcf_observed[[1]], y = pcf_observed[[3]], - type = "l", col = "black", xlab = "r", ylab = "g(r)") - - graphics::abline(h = 1, lty = 2, col = "grey") - - graphics::lines(x = pcf_relocated[[1]], y = pcf_relocated[[3]], col = "red") - - graphics::legend("topright", legend = c("observed", "reconstructed"), - col = c("black", "red"), lty = 1, inset = 0.025) - - } - - # increase counter no change - } else { - - energy_counter <- energy_counter + 1 - - } - - # count iterations - iterations <- iterations + 1 - - # save energy in data frame - energy_df[iterations, 2] <- energy_current - - # print progress - if (verbose) { - - message("\r> Progress: n_random: ", current_pattern, "/", n_random, - " || max_runs: ", floor(i / max_runs * 100), "%", - " || energy = ", round(energy_current, 5), "\t\t", - appendLF = FALSE) - - } - - # exit loop if e threshold or no_change counter max is reached - if (energy_current <= e_threshold || energy_counter > no_change) { - - # set stop criterion due to energy - stop_criterion_list[[current_pattern]] <- "e_threshold/no_change" - - break - - } - } - - # close plotting device - if (plot) { - - grDevices::dev.off() - - } - - # remove NAs if stopped due to energy - if (stop_criterion_list[[current_pattern]] == "e_threshold/no_change") { - - energy_df <- energy_df[1:iterations, ] - - } - - # save results in lists - energy_list[[current_pattern]] <- energy_df - iterations_list[[current_pattern]] <- iterations - result_list[[current_pattern]] <- simulated_current - } - - # combine to one list - reconstruction <- list(randomized = result_list, observed = pattern, - method = "reconstruct_pattern_cluster()", - energy_df = energy_list, stop_criterion = stop_criterion_list, - iterations = iterations_list) - - # set class of result - class(reconstruction) <- "rd_pat" - - # remove input if return_input = FALSE - if (!return_input) { - - # set observed to NA - reconstruction$observed <- "NA" - - # check if output should be simplified - if (simplify) { - - # not possible if more than one pattern is present - if (n_random > 1 && verbose) { - - warning("'simplify = TRUE' not possible for 'n_random > 1'.", - call. = FALSE) - - } - - # only one random pattern is present that should be returend - else if (n_random == 1) { - - reconstruction <- reconstruction$randomized[[1]] - - } - } - - # return input if return_input = TRUE - } else { - - # return warning if simply = TRUE because not possible if return_input = TRUE (only verbose = TRUE) - if (simplify && verbose) { - - warning("'simplify = TRUE' not possible for 'return_input = TRUE'.", call. = FALSE) - - } - } - - # write result in new line if progress was printed - if (verbose) { - - message("\r") - - } - - return(reconstruction) -} diff --git a/R/reconstruct_pattern_hetero.R b/R/reconstruct_pattern_hetero.R deleted file mode 100644 index d1a738e9..00000000 --- a/R/reconstruct_pattern_hetero.R +++ /dev/null @@ -1,383 +0,0 @@ -#' reconstruct_pattern_hetero -#' -#' @description Pattern reconstruction for heterogeneous patterns -#' -#' @param pattern ppp object with pattern. -#' @param n_random Integer with number of randomizations. -#' @param e_threshold Double with minimum energy to stop reconstruction. -#' @param max_runs Integer with maximum number of iterations if \code{e_threshold} -#' is not reached. -#' @param no_change Integer with number of iterations at which the reconstruction will -#' stop if the energy does not decrease. -#' @param annealing Double with probability to keep relocated point even if energy -#' did not decrease. -#' @param comp_fast Integer with threshold at which summary functions are estimated -#' in a computational fast way. -#' @param weights Vector with weights used to calculate energy. -#' The first number refers to Gest(r), the second number to pcf(r). -#' @param r_length Integer with number of intervals from \code{r = 0} to \code{r = rmax} for which -#' the summary functions are evaluated. -#' @param r_max Double with maximum distance used during calculation of summary functions. If \code{NULL}, -#' will be estimated from data. -#' @param return_input Logical if the original input data is returned. -#' @param simplify Logical if only pattern will be returned if \code{n_random = 1} -#' and \code{return_input = FALSE}. -#' @param verbose Logical if progress report is printed. -#' @param plot Logical if pcf(r) function is plotted and updated during optimization. -#' -#' @return rd_pat -#' -#' @examples -#' \dontrun{ -#' input_pattern <- spatstat.random::rpoispp(lambda = function(x, y) {100 * exp(-3 * x)}, -#' nsim = 1) -#' pattern_recon <- reconstruct_pattern_hetero(input_pattern, n_random = 19, max_runs = 1000) -#' } -#' -#' @aliases reconstruct_pattern_hetero -#' @rdname reconstruct_pattern_hetero -#' -#' @references -#' Kirkpatrick, S., Gelatt, C.D.Jr., Vecchi, M.P., 1983. Optimization by simulated -#' annealing. Science 220, 671–680. -#' -#' Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random point -#' patterns. Computational Statistics and Data Analysis 51, 859–871. -#' -#' -#' Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern analysis in -#' ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN 978-1-4200-8254-8 -#' -#' @keywords internal -reconstruct_pattern_hetero <- function(pattern, - n_random = 1, - e_threshold = 0.01, - max_runs = 1000, - no_change = Inf, - annealing = 0.01, - comp_fast = 1000, - weights = c(0.5, 0.5), - r_length = 250, - r_max = NULL, - return_input = TRUE, - simplify = FALSE, - verbose = TRUE, - plot = FALSE){ - - # check if n_random is >= 1 - if (n_random < 1) { - - stop("n_random must be >= 1.", call. = FALSE) - - } - - # check if number of points exceed comp_fast limit - if (pattern$n > comp_fast) { - - # Print message that summary functions will be computed fast - if (verbose) { - - message("> Using fast compuation of summary functions.") - - } - - comp_fast <- TRUE - - } else { - - comp_fast <- FALSE - - } - - # set names of randomization randomized_1 ... randomized_n - names_randomization <- paste0("randomized_", seq_len(n_random)) - - # create empty lists for results - energy_list <- vector("list", length = n_random) - iterations_list <- vector("list", length = n_random) - stop_criterion_list <- as.list(rep("max_runs", times = n_random)) - result_list <- vector("list", length = n_random) - - # set names - names(energy_list) <- names_randomization - names(iterations_list) <- names_randomization - names(stop_criterion_list) <- names_randomization - names(result_list) <- names_randomization - - # check if weights make sense - if (sum(weights) > 1 || sum(weights) == 0) { - - stop("The sum of 'weights' must be 0 < sum(weights) <= 1.", call. = FALSE) - - } - - # unmark pattern - if (spatstat.geom::is.marked(pattern)) { - - pattern <- spatstat.geom::unmark(pattern) - - if (verbose) { - - warning("Unmarked provided input pattern. For marked pattern, see reconstruct_pattern_marks().", - call. = FALSE) - - } - } - - # calculate r from data - if (is.null(r_max)) { - - r <- seq(from = 0, to = spatstat.explore::rmax.rule(W = pattern$window, lambda = spatstat.geom::intensity.ppp(pattern)), - length.out = r_length) - - # use provided r_max - } else { - - r <- seq(from = 0, to = r_max, length.out = r_length) - - } - - # estimate lambda(x,y) - lambda <- spatstat.explore::density.ppp(pattern) - - # create starting pattern - simulated <- spatstat.random::rpoint(n = pattern$n, f = lambda, win = pattern$window) - - - # fast computation of summary functions - if (comp_fast) { - - gest_observed <- spatstat.explore::Gest(pattern, correction = "none", r = r) - - gest_simulated <- spatstat.explore::Gest(simulated, correction = "none", r = r) - - pcf_observed <- estimate_pcf_fast(pattern, correction = "none", - method = "c", spar = 0.5, r = r) - - pcf_simulated <- estimate_pcf_fast(simulated, correction = "none", - method = "c", spar = 0.5, r = r) - - # normal computation of summary functions - } else { - - gest_observed <- spatstat.explore::Gest(X = pattern, correction = "han", r = r) - - gest_simulated <- spatstat.explore::Gest(X = simulated, correction = "han", r = r) - - pcf_observed <- spatstat.explore::pcf.ppp(X = pattern, correction = "best", - divisor = "d", r = r) - - pcf_simulated <- spatstat.explore::pcf.ppp(X = simulated, correction = "best", - divisor = "d", r = r) - - } - - # energy before reconstruction - energy <- (mean(abs(gest_observed[[3]] - gest_simulated[[3]]), na.rm = TRUE) * weights[[1]]) + - (mean(abs(pcf_observed[[3]] - pcf_simulated[[3]]), na.rm = TRUE) * weights[[2]]) - - # create n_random recondstructed patterns - for (current_pattern in seq_len(n_random)) { - - # current simulated - simulated_current <- simulated - energy_current <- energy - - # counter of iterations - iterations <- 0 - - # counter if energy changed - energy_counter <- 0 - - # df for energy - energy_df <- data.frame(i = seq(from = 1, to = max_runs, by = 1), energy = NA) - - # random ids of pattern - rp_id <- sample(x = seq_len(simulated_current$n), size = max_runs, replace = TRUE) - - # create random new points - rp_coords <- spatstat.random::rpoint(n = max_runs, f = lambda, win = pattern$window) - - # create random number for annealing prob - if (annealing != 0) { - - random_annealing <- stats::runif(n = max_runs, min = 0, max = 1) - - } else { - - random_annealing <- rep(0, max_runs) - - } - - # pattern reconstruction algorithm (optimaztion of energy) - not longer than max_runs - for (i in seq_len(max_runs)) { - - # data for relocation - relocated <- simulated_current - - # get current point id - rp_id_current <- rp_id[[i]] - - # relocate point - relocated$x[[rp_id_current]] <- rp_coords$x[[i]] - - relocated$y[[rp_id_current]] <- rp_coords$y[[i]] - - # calculate summary functions after relocation - if (comp_fast) { - - gest_relocated <- spatstat.explore::Gest(relocated, correction = "none", r = r) - - pcf_relocated <- estimate_pcf_fast(relocated, correction = "none", - method = "c", spar = 0.5, r = r) - - } else { - - gest_relocated <- spatstat.explore::Gest(X = relocated, correction = "han", r = r) - - pcf_relocated <- spatstat.explore::pcf.ppp(X = relocated, correction = "best", - divisor = "d", r = r) - - } - - # energy after relocation - energy_relocated <- (mean(abs(gest_observed[[3]] - gest_relocated[[3]]), na.rm = TRUE) * weights[[1]]) + - (mean(abs(pcf_observed[[3]] - pcf_relocated[[3]]), na.rm = TRUE) * weights[[2]]) - - # lower energy after relocation - if (energy_relocated < energy_current || random_annealing[i] < annealing) { - - # keep relocated pattern - simulated_current <- relocated - - # keep energy_relocated as energy - energy_current <- energy_relocated - - # set counter since last change back to 0 - energy_counter <- 0 - - # plot observed vs reconstructed - if (plot) { - - # https://support.rstudio.com/hc/en-us/community/posts/200661917-Graph-does-not-update-until-loop-completion - Sys.sleep(0.01) - - graphics::plot(x = pcf_observed[[1]], y = pcf_observed[[3]], - type = "l", col = "black", xlab = "r", ylab = "g(r)") - - graphics::abline(h = 1, lty = 2, col = "grey") - - graphics::lines(x = pcf_relocated[[1]], y = pcf_relocated[[3]], col = "red") - - graphics::legend("topright", legend = c("observed", "reconstructed"), - col = c("black", "red"), lty = 1, inset = 0.025) - - } - - # increase counter no change - } else { - - energy_counter <- energy_counter + 1 - - } - - # count iterations - iterations <- iterations + 1 - - # save energy in data frame - energy_df[iterations, 2] <- energy_current - - # print progress - if (verbose) { - - message("\r> Progress: n_random: ", current_pattern, "/", n_random, - " || max_runs: ", floor(i / max_runs * 100), "%", - " || energy = ", round(energy_current, 5), "\t\t", - appendLF = FALSE) - - } - - # exit loop if e threshold or no_change counter max is reached - if (energy_current <= e_threshold || energy_counter > no_change) { - - # set stop criterion due to energy - stop_criterion_list[[current_pattern]] <- "e_threshold/no_change" - - break - - } - } - - # close plotting device - if (plot) { - - grDevices::dev.off() - - } - - # remove NAs if stopped due to energy - if (stop_criterion_list[[current_pattern]] == "e_threshold/no_change") { - - energy_df <- energy_df[1:iterations, ] - - } - - # save results in lists - energy_list[[current_pattern]] <- energy_df - iterations_list[[current_pattern]] <- iterations - result_list[[current_pattern]] <- simulated_current - - } - - # combine to one list - reconstruction <- list(randomized = result_list, observed = pattern, - method = "reconstruct_pattern_hetero()", - energy_df = energy_list, stop_criterion = stop_criterion_list, - iterations = iterations_list) - - # set class of result - class(reconstruction) <- "rd_pat" - - # remove input if return_input = FALSE - if (!return_input) { - - # set observed to NA - reconstruction$observed <- "NA" - - # check if output should be simplified - if (simplify) { - - # not possible if more than one pattern is present - if (n_random > 1 && verbose) { - - warning("'simplify = TRUE' not possible for 'n_random > 1'.", - call. = FALSE) - - # only one random pattern is present that should be returend - } else if (n_random == 1) { - - reconstruction <- reconstruction$randomized[[1]] - - } - } - - # return input if return_input = TRUE - } else { - - # return warning if simply = TRUE because not possible if return_input = TRUE (only verbose = TRUE) - if (simplify && verbose) { - - warning("'simplify = TRUE' not possible for 'return_input = TRUE'.", call. = FALSE) - - } - } - - # write result in new line if progress was printed - if (verbose) { - - message("\r") - - } - - return(reconstruction) -} diff --git a/R/reconstruct_pattern_homo.R b/R/reconstruct_pattern_homo.R deleted file mode 100644 index d7a96656..00000000 --- a/R/reconstruct_pattern_homo.R +++ /dev/null @@ -1,407 +0,0 @@ -#' reconstruct_pattern_homo -#' -#' @description Pattern reconstruction for homogeneous pattern -#' -#' @param pattern ppp object with pattern. -#' @param n_random Integer with number of randomizations. -#' @param e_threshold Double with minimum energy to stop reconstruction. -#' @param max_runs Integer with maximum number of iterations if \code{e_threshold} -#' is not reached. -#' @param no_change Integer with number of iterations at which the reconstruction will -#' stop if the energy does not decrease. -#' @param annealing Double with probability to keep relocated point even if energy -#' did not decrease. -#' @param n_points Integer with number of points to be simulated. -#' @param window owin object with window of simulated pattern. -#' @param comp_fast Integer with threshold at which summary functions are estimated -#' in a computational fast way. -#' @param weights Vector with weights used to calculate energy. -#' The first number refers to Gest(r), the second number to pcf(r). -#' @param r_length Integer with number of intervals from \code{r = 0} to \code{r = rmax} for which -#' the summary functions are evaluated. -#' @param r_max Double with maximum distance used during calculation of summary functions. If \code{NULL}, -#' will be estimated from data. -#' @param return_input Logical if the original input data is returned. -#' @param simplify Logical if only pattern will be returned if \code{n_random = 1} -#' and \code{return_input = FALSE}. -#' @param verbose Logical if progress report is printed. -#' @param plot Logical if pcf(r) function is plotted and updated during optimization. -#' -#' @return rd_pat -#' -#' @examples -#' \dontrun{ -#' pattern_recon_a <- reconstruct_pattern_homo(species_a, n_random = 19, -#' max_runs = 1000) -#' -#' pattern_recon_b <- reconstruct_pattern_homo(species_a, n_points = 70, -#' n_random = 19, max_runs = 1000) -#' } -#' -#' @aliases reconstruct_pattern_homo -#' @rdname reconstruct_pattern_homo -#' -#' @references -#' Kirkpatrick, S., Gelatt, C.D.Jr., Vecchi, M.P., 1983. Optimization by simulated -#' annealing. Science 220, 671–680. -#' -#' Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random point -#' patterns. Computational Statistics and Data Analysis 51, 859–871. -#' -#' -#' Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern analysis in -#' ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN 978-1-4200-8254-8 -#' -#' @keywords internal -reconstruct_pattern_homo <- function(pattern, - n_random = 1, - e_threshold = 0.01, - max_runs = 1000, - no_change = Inf, - annealing = 0.01, - n_points = NULL, - window = NULL, - comp_fast = 1000, - weights = c(0.5, 0.5), - r_length = 250, - r_max = NULL, - return_input = TRUE, - simplify = FALSE, - verbose = TRUE, - plot = FALSE){ - - # check if n_random is >= 1 - if (n_random < 1) { - - stop("n_random must be >= 1.", call. = FALSE) - - } - - # use number of points of pattern if not provided - if (is.null(n_points)) { - - message("> Using number of points 'pattern'.") - - n_points <- pattern$n - - } - - # use window of pattern if not provided - if (is.null(window)) { - - message("> Using window of 'pattern'.") - - window <- pattern$window - - } - - # calculate intensity - intensity <- n_points / spatstat.geom::area(window) - - # check if number of points exceed comp_fast limit - if (n_points > comp_fast) { - - # Print message that summary functions will be computed fast - if (verbose) { - - message("> Using fast compuation of summary functions.") - - } - - comp_fast <- TRUE - - } else { - - comp_fast <- FALSE - - } - - # set names of randomization randomized_1 ... randomized_n - names_randomization <- paste0("randomized_", seq_len(n_random)) - - # create empty lists for results - energy_list <- vector("list", length = n_random) - iterations_list <- vector("list", length = n_random) - stop_criterion_list <- as.list(rep("max_runs", times = n_random)) - result_list <- vector("list", length = n_random) - - # set names - names(energy_list) <- names_randomization - names(iterations_list) <- names_randomization - names(stop_criterion_list) <- names_randomization - names(result_list) <- names_randomization - - # check if weights make sense - if (sum(weights) > 1 || sum(weights) == 0) { - - stop("The sum of 'weights' must be 0 < sum(weights) <= 1.", call. = FALSE) - - } - - # unmark pattern - if (spatstat.geom::is.marked(pattern)) { - - pattern <- spatstat.geom::unmark(pattern) - - if (verbose) { - warning("Unmarked provided input pattern. For marked pattern, see reconstruct_pattern_marks().", - call. = FALSE) - - } - } - - # calculate r from data - if (is.null(r_max)) { - - r <- seq(from = 0, to = spatstat.explore::rmax.rule(W = window, lambda = intensity), - length.out = r_length) - - # use provided r_max - } else { - - r <- seq(from = 0, to = r_max, length.out = r_length) - - } - - # create Poisson simulation data - simulated <- spatstat.random::runifpoint(n = n_points, nsim = 1, drop = TRUE, - win = window, warn = FALSE) - - # fast computation of summary functions - if (comp_fast) { - - gest_observed <- spatstat.explore::Gest(pattern, correction = "none", r = r) - - gest_simulated <- spatstat.explore::Gest(simulated, correction = "none", r = r) - - pcf_observed <- estimate_pcf_fast(pattern, correction = "none", - method = "c", spar = 0.5, r = r) - - pcf_simulated <- estimate_pcf_fast(simulated, correction = "none", - method = "c", spar = 0.5, r = r) - - # normal computation of summary functions - } else { - - gest_observed <- spatstat.explore::Gest(X = pattern, correction = "han", r = r) - - gest_simulated <- spatstat.explore::Gest(X = simulated, correction = "han", r = r) - - pcf_observed <- spatstat.explore::pcf.ppp(X = pattern, correction = "best", - divisor = "d", r = r) - - pcf_simulated <- spatstat.explore::pcf.ppp(X = simulated, correction = "best", - divisor = "d", r = r) - - } - - # energy before reconstruction - energy <- (mean(abs(gest_observed[[3]] - gest_simulated[[3]]), na.rm = TRUE) * weights[[1]]) + - (mean(abs(pcf_observed[[3]] - pcf_simulated[[3]]), na.rm = TRUE) * weights[[2]]) - - # create n_random recondstructed patterns - for (current_pattern in seq_len(n_random)) { - - # current simulated - simulated_current <- simulated - energy_current <- energy - - # counter of iterations - iterations <- 0 - - # counter if energy changed - energy_counter <- 0 - - # df for energy - energy_df <- data.frame(i = seq(from = 1, to = max_runs, by = 1), energy = NA) - - # random ids of pattern - rp_id <- sample(x = seq_len(simulated_current$n), size = max_runs, replace = TRUE) - - # create random new points - rp_coords <- spatstat.random::runifpoint(n = max_runs, nsim = 1, drop = TRUE, - win = simulated_current$window, - warn = FALSE) - - # create random number for annealing prob - if (annealing != 0) { - - random_annealing <- stats::runif(n = max_runs, min = 0, max = 1) - - } else { - - random_annealing <- rep(0, max_runs) - - } - - # pattern reconstruction algorithm (optimaztion of energy) - not longer than max_runs - for (i in seq_len(max_runs)) { - - # data for relocation - relocated <- simulated_current - - # get current point id - rp_id_current <- rp_id[[i]] - - # relocate point - relocated$x[[rp_id_current]] <- rp_coords$x[[i]] - - relocated$y[[rp_id_current]] <- rp_coords$y[[i]] - - # calculate summary functions after relocation - if (comp_fast) { - - gest_relocated <- spatstat.explore::Gest(relocated, correction = "none", r = r) - - pcf_relocated <- estimate_pcf_fast(relocated, correction = "none", - method = "c", spar = 0.5, r = r) - } else { - - gest_relocated <- spatstat.explore::Gest(X = relocated, correction = "han", r = r) - - pcf_relocated <- spatstat.explore::pcf.ppp(X = relocated, correction = "best", - divisor = "d", r = r) - - } - - # energy after relocation - energy_relocated <- (mean(abs(gest_observed[[3]] - gest_relocated[[3]]), na.rm = TRUE) * weights[[1]]) + - (mean(abs(pcf_observed[[3]] - pcf_relocated[[3]]), na.rm = TRUE) * weights[[2]]) - - # lower energy after relocation - if (energy_relocated < energy_current || random_annealing[i] < annealing) { - - # keep relocated pattern - simulated_current <- relocated - - # keep energy_relocated as energy - energy_current <- energy_relocated - - # set counter since last change back to 0 - energy_counter <- 0 - - # plot observed vs reconstructed - if (plot) { - - # https://support.rstudio.com/hc/en-us/community/posts/200661917-Graph-does-not-update-until-loop-completion - Sys.sleep(0.01) - - graphics::plot(x = pcf_observed[[1]], y = pcf_observed[[3]], - type = "l", col = "black", xlab = "r", ylab = "g(r)") - - graphics::abline(h = 1, lty = 2, col = "grey") - - graphics::lines(x = pcf_relocated[[1]], y = pcf_relocated[[3]], col = "red") - - graphics::legend("topright", legend = c("observed", "reconstructed"), - col = c("black", "red"), lty = 1, inset = 0.025) - - } - - # increase counter no change - } else { - - energy_counter <- energy_counter + 1 - - } - - # count iterations - iterations <- iterations + 1 - - # save energy in data frame - energy_df[iterations, 2] <- energy_current - - # print progress - if (verbose) { - - message("\r> Progress: n_random: ", current_pattern, "/", n_random, - " || max_runs: ", floor(i / max_runs * 100), "%", - " || energy = ", round(energy_current, 5), "\t\t", - appendLF = FALSE) - - } - - # exit loop if e threshold or no_change counter max is reached - if (energy_current <= e_threshold || energy_counter > no_change) { - - # set stop criterion due to energy - stop_criterion_list[[current_pattern]] <- "e_threshold/no_change" - - break - - } - } - - # close plotting device - if (plot) { - - grDevices::dev.off() - - } - - # remove NAs if stopped due to energy - if (stop_criterion_list[[current_pattern]] == "e_threshold/no_change") { - - energy_df <- energy_df[1:iterations, ] - - } - - # save results in lists - energy_list[[current_pattern]] <- energy_df - iterations_list[[current_pattern]] <- iterations - result_list[[current_pattern]] <- simulated_current - - } - - # combine to one list - reconstruction <- list(randomized = result_list, observed = pattern, - method = "reconstruct_pattern_homo()", - energy_df = energy_list, stop_criterion = stop_criterion_list, - iterations = iterations_list) - - # set class of result - class(reconstruction) <- "rd_pat" - - # remove input if return_input = FALSE - if (!return_input) { - - # set observed to NA - reconstruction$observed <- "NA" - - # check if output should be simplified - if (simplify) { - - # not possible if more than one pattern is present - if (n_random > 1 && verbose) { - - warning("'simplify = TRUE' not possible for 'n_random > 1'.", - call. = FALSE) - - # only one random pattern is present that should be returend - } else if (n_random == 1) { - - reconstruction <- reconstruction$randomized[[1]] - - } - } - - # return input if return_input = TRUE - } else { - - # return warning if simply = TRUE because not possible if return_input = TRUE (only verbose = TRUE) - if (simplify && verbose) { - - warning("'simplify = TRUE' not possible for 'return_input = TRUE'.", call. = FALSE) - - } - } - - # write result in new line if progress was printed - if (verbose) { - - message("\r") - - } - - return(reconstruction) -} diff --git a/R/reconstruct_pattern_marks.R b/R/reconstruct_pattern_marks.R index cfd7a484..2bc55026 100644 --- a/R/reconstruct_pattern_marks.R +++ b/R/reconstruct_pattern_marks.R @@ -99,60 +99,48 @@ reconstruct_pattern_marks <- function(pattern, } - # set names of randomization randomized_1 ... randomized_n - names_randomization <- paste0("randomized_", seq_len(n_random)) - - # create empty lists for results - energy_list <- vector("list", length = n_random) - iterations_list <- vector("list", length = n_random) - stop_criterion <- as.list(rep("max_runs", times = n_random)) - result_list <- vector("list", length = n_random) - - # set names - names(energy_list) <- names_randomization - names(result_list) <- names_randomization - names(iterations_list) <- names_randomization - names(stop_criterion) <- names_randomization - # calculate r from data if (is.null(r_max)) { r <- seq(from = 0, to = spatstat.explore::rmax.rule(W = pattern$window, lambda = spatstat.geom::intensity.ppp(pattern)), length.out = r_length) - # use provided r_max + # use provided r_max } else { r <- seq(from = 0, to = r_max, length.out = r_length) } - # create random pattern - simulated <- pattern - - # assign shuffled marks to pattern - spatstat.geom::marks(simulated) <- sample(x = marked_pattern$marks, size = simulated$n, - replace = TRUE) + # set names of randomization randomized_1 ... randomized_n + names_randomization <- paste0("randomized_", seq_len(n_random)) - # calculate summary functions - kmmr_observed <- spatstat.explore::markcorr(marked_pattern, correction = "Ripley", - r = r) + # create empty lists for results + energy_list <- vector("list", length = n_random) + iterations_list <- vector("list", length = n_random) + stop_criterion_list <- as.list(rep("max_runs", times = n_random)) + result_list <- vector("list", length = n_random) - kmmr_simulated <- spatstat.explore::markcorr(simulated, correction = "Ripley", - r = r) + # set names + names(energy_list) <- names_randomization + names(iterations_list) <- names_randomization + names(stop_criterion_list) <- names_randomization + names(result_list) <- names_randomization - # energy before reconstruction - energy <- mean(abs(kmmr_observed[[3]] - kmmr_simulated[[3]]), na.rm = TRUE) + # calculate summary functions + kmmr_observed <- spatstat.explore::markcorr(marked_pattern, correction = "none", r = r) # create n_random recondstructed patterns - for (current_pattern in seq_len(n_random)) { + for (i in seq_len(n_random)) { + + # create random pattern + simulated <- pattern - # current simulated - simulated_current <- simulated - energy_current <- energy + # assign shuffled marks to pattern + spatstat.geom::marks(simulated) <- sample(x = marked_pattern$marks, size = simulated$n, + replace = TRUE) - # counter of iterations - iterations <- 0 + energy <- Inf # counter if energy changed energy_counter <- 0 @@ -161,11 +149,6 @@ reconstruct_pattern_marks <- function(pattern, energy_df <- data.frame(i = seq(from = 1, to = max_runs, by = 1), energy = NA) - # get two random points to switch marks - rp_a <- sample(x = seq_len(simulated_current$n), size = max_runs, replace = TRUE) - - rp_b <- sample(x = seq_len(simulated_current$n), size = max_runs, replace = TRUE) - # create random number for annealing prob if (annealing != 0) { @@ -177,41 +160,46 @@ reconstruct_pattern_marks <- function(pattern, } + # get two random points to switch marks + rp_a <- sample(x = seq_len(simulated$n), size = max_runs, replace = TRUE) + + rp_b <- sample(x = seq_len(simulated$n), size = max_runs, replace = TRUE) + # pattern reconstruction algorithm (optimaztion of energy) - not longer than max_runs - for (i in seq_len(max_runs)) { + for (j in seq_len(max_runs)) { - relocated <- simulated_current # data for relocation + relocated <- simulated # data for relocation # current random points - rp_a_current <- rp_a[[i]] + a_current <- rp_a[[j]] - rp_b_current <- rp_b[[i]] + b_current <- rp_b[[j]] # get marks of the two random points - mark_a <- relocated$marks[[rp_a_current]] + mark_a <- relocated$marks[[a_current]] - mark_b <- relocated$marks[[rp_b_current]] + mark_b <- relocated$marks[[b_current]] # switch the marks of the two points - relocated$marks[[rp_a_current]] <- mark_b + relocated$marks[[a_current]] <- mark_b - relocated$marks[[rp_b_current]] <- mark_a + relocated$marks[[b_current]] <- mark_a # calculate summary functions after relocation - kmmr_relocated <- spatstat.explore::markcorr(relocated, correction = "Ripley", - r = r) + kmmr_relocated <- spatstat.explore::markcorr(relocated, correction = "none", + r = r) # energy after relocation energy_relocated <- mean(abs(kmmr_observed[[3]] - kmmr_relocated[[3]]), na.rm = TRUE) # lower energy after relocation - if (energy_relocated < energy_current || random_annealing[i] < annealing) { + if (energy_relocated < energy || random_annealing[i] < annealing) { # keep relocated pattern - simulated_current <- relocated + simulated <- relocated # keep energy_relocated as energy - energy_current <- energy_relocated + energy <- energy_relocated # plot observed vs reconstructed if (plot) { @@ -238,27 +226,25 @@ reconstruct_pattern_marks <- function(pattern, } - # count iterations - iterations <- iterations + 1 - # save energy in data frame - energy_df[iterations, 2] <- energy_current + energy_df[j, 2] <- energy # print progress if (verbose) { - message("\r> Progress: n_random: ", current_pattern, "/", n_random, - " || max_runs: ", floor(i / max_runs * 100), "%", - " || energy = ", round(energy_current, 5), "\t\t", + message("\r> Progress: n_random: ", i, "/", n_random, + " || max_runs: ", floor(j / max_runs * 100), "%", + " || energy = ", round(energy, 5), "\t\t", appendLF = FALSE) } # exit loop if e threshold or no_change counter max is reached - if (energy_current <= e_threshold || energy_counter > no_change) { + if (energy <= e_threshold || energy_counter > no_change) { # set stop criterion due to energy - stop_criterion[[current_pattern]] <- "e_threshold/no_change" + stop_criterion_list[[i]] <- ifelse(test = energy <= e_threshold, + yes = "e_threshold", no = "no_change") break @@ -272,23 +258,30 @@ reconstruct_pattern_marks <- function(pattern, } # remove NAs if stopped due to energy - if (stop_criterion[[current_pattern]] == "e_threshold/no_change") { + if (stop_criterion_list[[i]] %in% c("e_threshold", "no_change")) { - energy_df <- energy_df[1:iterations, ] + energy_df <- energy_df[1:j, ] } # save results in lists - energy_list[[current_pattern]] <- energy_df - iterations_list[[current_pattern]] <- iterations - result_list[[current_pattern]] <- simulated_current + energy_list[[i]] <- energy_df + iterations_list[[i]] <- j + result_list[[i]] <- simulated + + } + + # write result in new line if progress was printed + if (verbose) { + + message("\r") } # combine to one list reconstruction <- list(randomized = result_list, observed = marked_pattern, - method = "reconstruct_pattern_marks()", - energy_df = energy_list, stop_criterion = stop_criterion, + method = "marks", + energy_df = energy_list, stop_criterion = stop_criterion_list, iterations = iterations_list) # set class of returning object @@ -330,12 +323,5 @@ reconstruct_pattern_marks <- function(pattern, } } - # write result in new line if progress was printed - if (verbose) { - - message("\r") - - } - return(reconstruction) } diff --git a/R/results_habitat_association.R b/R/results_habitat_association.R index 682a9d17..bd9b155b 100644 --- a/R/results_habitat_association.R +++ b/R/results_habitat_association.R @@ -81,7 +81,7 @@ results_habitat_association <- function(pattern, raster, significance_level = 0. if (inherits(x = raster, what = "rd_ras")) { # check if randomized and observed is present - if (!methods::is(raster$observed, "SpatRaster")) { + if (!inherits(x = raster$observed, what = "SpatRaster")) { stop("The observed raster needs to be included in the input 'raster'.", call. = FALSE) diff --git a/R/sysdata.rda b/R/sysdata.rda index bfa6c94a..ae051883 100644 Binary files a/R/sysdata.rda and b/R/sysdata.rda differ diff --git a/R/unpack_randomized.R b/R/unpack_randomized.R index dd634ce9..61338875 100644 --- a/R/unpack_randomized.R +++ b/R/unpack_randomized.R @@ -29,11 +29,12 @@ #' @export unpack_randomized <- function(raster) { - # wrap observerd raster - raster$observed <- terra::rast(raster$observed) + # check if observed is present + # unwrap observed raster + if (inherits(x = raster$observed, what = "PackedSpatRaster")) raster$observed <- terra::unwrap(raster$observed) # wrap all randomized raster - raster$randomized <- lapply(X = raster$randomized, FUN = terra::rast) + raster$randomized <- lapply(X = raster$randomized, FUN = terra::unwrap) return(raster) diff --git a/README.Rmd b/README.Rmd index 6604274d..cd305831 100644 --- a/README.Rmd +++ b/README.Rmd @@ -26,15 +26,9 @@ README Last updated: `r Sys.Date()` -**S**pecies-**h**abitat **a**ssociations in **R** is a `R` package to analyze species-habitat associations. Therefore, information about the location of the species is needed (as a point pattern) and about the environmental conditions (as a raster map). In order to analyse the data for significant habitat associations either the location data or the environmental data is randomized n-times. Then, counts within the habitats are compared between the randomized data and the observed data. Positive or negative associations are present if the observed counts is higher or lower than the randomized counts (using quantile thresholds). Methods are described in Plotkin et al. (2000), Harms et al. (2001) and Wiegand & Moloney (2014). **shar** is mainly based on the [`spatstat`](http://spatstat.org) (Baddeley et al. 2015) and [`terra`](https://rspatial.org/terra/) (Hijmans 2022) package. +**S**pecies-**h**abitat **a**ssociations in **R** provides a toolset of functions in the `R` programming language to analyze species-habitat associations. Therefore, information about the location of the species (as a point pattern) and the environmental conditions (as a raster) is needed. In order to analyse the data for significant habitat associations either the location data or the environmental data is randomized *n*-times. Then, counts within the habitats are compared between the observed and the randomized data. Positive or negative associations are present if the observed counts are higher or lower than the randomized counts (using quantile thresholds). Methods are described in Plotkin et al. (2000), Harms et al. (2001) and Wiegand & Moloney (2014). **shar** is mainly based on the [`spatstat`](http://spatstat.org) (Baddeley et al. 2015) and [`terra`](https://rspatial.org/terra/) (Hijmans 2022) package. -#### Citation - -The **shar** package is part of our academic work. To cite the package or acknowledge its use in publications, please cite the following paper. - -> Hesselbarth, M.H.K., (2021). shar: A R package to analyze species-habitat associations using point pattern analysis. Journal of Open Source Software, 6(67), 3811. https://doi.org/10.21105/joss.03811 - -The get a BibTex entry, please use `citation("shar")`. +You can find more information and help using the corresponding [homepage](https://r-spatialecology.github.io/shar/). ## Installation @@ -51,79 +45,19 @@ And the development version from [GitHub](https://github.com/r-spatialecology/sh remotes::install_github("r-spatialecology/shar") ``` -This also automatically installs all non-base `R` package dependencies, namely the following packages: `classInt`, `raster`, `spatstat.explore`, `spatstat.model`, `spatstat.geom`. +This also automatically installs all non-base `R` package dependencies, namely: `classInt`, `spatstat.explore`, `spatstat.geom`, `spatstat.model`, `spatstat.random`, and `terra`. -## How to use shar - -```{r import-libs, message = FALSE, warning = FALSE} -library(shar) -library(spatstat) -library(terra) - -set.seed(42) -``` - -**shar** comes with build-in example data sets. `species_a` and `species_b` are exemplary location of species, e.g. trees, as `ppp`-objects from the `spatstat` package. `landscape` contains exemplary continuous environmental data. However, all methods depend on discrete data. Therefore we need to classify the data first. -However, all methods require "fully mapped data" in a sense that NA cells of the environmental data are allowed only if simultaneously these areas cannot accommodate any locations of the point pattern (e.g., a water body within a forest area). This needs to be reflected in the observation window of the point pattern. For the torus translation method, no NA values are allowed at all. - -```{r environmental-data} -landscape_classified <- classify_habitats(raster = terra::rast(landscape), n = 5, style = "fisher") -``` +## How to use **shar** -There are two possibilities to randomize the environmental data, both described in Harms et al. (2001). The first shifts the habitat map in all 4 cardinal directions around a torus. The second one assigns the habitat values to an empty map using a random walk algorithm. Both functions return a list with randomized rasters and the observed one. For more information on the methods, please click [here](https://r-spatialecology.github.io/shar/articles/articles/background.html). - -```{r habitat_random, eval = FALSE} -torus_trans <- translate_raster(raster = landscape_classified) - -random_walk <- randomize_raster(raster = landscape_classified, n_random = 39) -``` - -To plot the randomized raster, you can use the plot function and specify the number of raster as as well as the color palette used for the discrete environmental data. - -```{r plot_habitat-random, eval = FALSE} -col_palette <- c("#440154FF", "#3B528BFF", "#21908CFF", "#5DC863FF", "#FDE725FF") - -plot(torus_trans, n = 3, col = col_palette) -``` - -To randomize the point pattern, either use the Gamma test described by Plotkin et al. (2000) or pattern reconstruction (Kirkpatrick et al. 1983; Tscheschel & Stoyan 2006). - -```{r pattern-random, eval = FALSE} -gamma_test <- fit_point_process(pattern = species_b, process = "cluster", n_random = 39) - -# (this can takes some time) -reconstruction <- reconstruct_pattern(pattern = species_b, n_random = 39, e_threshold = 0.05) -``` +Please refer to `vignette("Get started")` and the [homepage](https://r-spatialecology.github.io/shar/) to get an introduction to **shar**. -Of course, there are several utility functions. For example, you can plot the summary function of the observed pattern and the simulation envelopes of randomized patterns (`what = "sf"`) or some randomized and the observed pattern (`what = "pp"`) using the plot function. +### Citation -```{r plot-random_pattern, fig.align = "center", out.height = "100%", out.width = "100%", message = FALSE} -plot(reconstruction, what = "pp") -``` - -Another utility function allows to calculate the differences between the observed pattern and the randomized patterns (also called energy using summary functions). - -```{r calculate-energy, message = FALSE} -calculate_energy(reconstruction, return_mean = TRUE) -``` - -The data was created that `species_a` has a negative association to habitat 4 and `species_b` has a positive association to habitat 5, which is reflected in the results. - -Given the characteristics of the method, a positive association to one habitat inevitably leads to a negative association to at least one of the other habitats (and vice versa; Yamada et al. 2006). For example, a high amount of individual points in the positively associated habitat simultaneously mean that less individual points can be present in the other habitats. - -Furthermore, please be aware that due to the randomization of the null model data, results might slightly differ between different randomization approaches (e.g., `fit_point_process()` vs. `translate_raster()`) and even for repetitions of the same approach. Thus, the exact `lo` and `hi` thresholds might be slightly different when re-running the examples. However, the counts of the observed data should be identical, and general results and trends should be similar. - -```{r result-unpack, echo = FALSE} -torus_trans <- unpack_randomized(torus_trans) -``` - -```{r results} -significance_level <- 0.01 +The **shar** package is part of our academic work. To cite the package or acknowledge its use in publications, please cite the following paper. -results_habitat_association(pattern = species_a, raster = torus_trans, significance_level = significance_level) +> Hesselbarth, M.H.K., (2021). shar: A R package to analyze species-habitat associations using point pattern analysis. Journal of Open Source Software, 6(67), 3811. https://doi.org/10.21105/joss.03811 -results_habitat_association(pattern = reconstruction, raster = landscape_classified, significance_level = significance_level) -``` +The get a BibTex entry, please use `citation("shar")`. ## Contributing and Code of Conduct @@ -133,20 +67,3 @@ Please note that the **shar** package is released with a [Contributor Code of Co To contribute to this project, please see the [Contributing guidelines](CONTRIBUTING.md). -### References - -Baddeley, A., Rubak, E., Turner, R., 2015. Spatial point patterns: Methodology and applications with R. Chapman and Hall/CRC Press, London. ISBN 978-1-4822-1020-0 - -Harms, K.E., Condit, R., Hubbell, S.P., Foster, R.B., 2001. Habitat associations of trees and shrubs in a 50-ha neotropical forest plot. Journal of Ecology 89, 947–959. - -Hijmans, R.J., 2022 terra: Spatial Data Analysis. R package version 1.6-7. . - -Kirkpatrick, S., Gelatt, C.D.Jr., Vecchi, M.P., 1983. Optimization by simulated annealing. Science 220, 671–680. - -Plotkin, J.B., Potts, M.D., Leslie, N., Manokaran, N., LaFrankie, J.V., Ashton, P.S., 2000. Species-area curves, spatial aggregation, and habitat specialization in tropical forests. Journal of Theoretical Biology 207, 81–99. - -Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random point patterns. Computational Statistics and Data Analysis 51, 859–871. - -Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern analysis in ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN 978-1-4200-8254-8 - -Yamada, T., Tomita, A., Itoh, A., Yamakura, T., Ohkubo, T., Kanzaki, M., Tan, S., Ashton, P.S., 2006. Habitat associations of Sterculiaceae trees in a Bornean rain forest plot. Journal of Vegetation Science 17, 559–566. diff --git a/README.md b/README.md index c533402d..f8370248 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ -README Last updated: 2023-03-08 +README Last updated: 2023-08-31 | CI | Development | CRAN | License | |--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------| @@ -14,31 +14,23 @@ README Last updated: 2023-03-08 -**S**pecies-**h**abitat **a**ssociations in **R** is a `R` package to -analyze species-habitat associations. Therefore, information about the -location of the species is needed (as a point pattern) and about the -environmental conditions (as a raster map). In order to analyse the data -for significant habitat associations either the location data or the -environmental data is randomized n-times. Then, counts within the -habitats are compared between the randomized data and the observed data. -Positive or negative associations are present if the observed counts is -higher or lower than the randomized counts (using quantile thresholds). -Methods are described in Plotkin et al. (2000), Harms et al. (2001) and -Wiegand & Moloney (2014). **shar** is mainly based on the +**S**pecies-**h**abitat **a**ssociations in **R** provides a toolset of +functions in the `R` programming language to analyze species-habitat +associations. Therefore, information about the location of the species +(as a point pattern) and the environmental conditions (as a raster) is +needed. In order to analyse the data for significant habitat +associations either the location data or the environmental data is +randomized *n*-times. Then, counts within the habitats are compared +between the observed and the randomized data. Positive or negative +associations are present if the observed counts are higher or lower than +the randomized counts (using quantile thresholds). Methods are described +in Plotkin et al. (2000), Harms et al. (2001) and Wiegand & Moloney +(2014). **shar** is mainly based on the [`spatstat`](http://spatstat.org) (Baddeley et al. 2015) and [`terra`](https://rspatial.org/terra/) (Hijmans 2022) package. -#### Citation - -The **shar** package is part of our academic work. To cite the package -or acknowledge its use in publications, please cite the following paper. - -> Hesselbarth, M.H.K., (2021). shar: A R package to analyze -> species-habitat associations using point pattern analysis. Journal of -> Open Source Software, 6(67), 3811. -> - -The get a BibTex entry, please use `citation("shar")`. +You can find more information and help using the corresponding +[homepage](https://r-spatialecology.github.io/shar/). ## Installation @@ -58,132 +50,26 @@ remotes::install_github("r-spatialecology/shar") ``` This also automatically installs all non-base `R` package dependencies, -namely the following packages: `classInt`, `raster`, `spatstat.explore`, -`spatstat.model`, `spatstat.geom`. +namely: `classInt`, `spatstat.explore`, `spatstat.geom`, +`spatstat.model`, `spatstat.random`, and `terra`. -## How to use shar +## How to use **shar** -``` r -library(shar) -library(spatstat) -library(terra) +Please refer to `vignette("Get started")` and the +[homepage](https://r-spatialecology.github.io/shar/) to get an +introduction to **shar**. -set.seed(42) -``` - -**shar** comes with build-in example data sets. `species_a` and -`species_b` are exemplary location of species, e.g. trees, as -`ppp`-objects from the `spatstat` package. `landscape` contains -exemplary continuous environmental data. However, all methods depend on -discrete data. Therefore we need to classify the data first. However, -all methods require “fully mapped data” in a sense that NA cells of the -environmental data are allowed only if simultaneously these areas cannot -accommodate any locations of the point pattern (e.g., a water body -within a forest area). This needs to be reflected in the observation -window of the point pattern. For the torus translation method, no NA -values are allowed at all. - -``` r -landscape_classified <- classify_habitats(raster = terra::rast(landscape), n = 5, style = "fisher") -``` - -There are two possibilities to randomize the environmental data, both -described in Harms et al. (2001). The first shifts the habitat map in -all 4 cardinal directions around a torus. The second one assigns the -habitat values to an empty map using a random walk algorithm. Both -functions return a list with randomized rasters and the observed one. -For more information on the methods, please click -[here](https://r-spatialecology.github.io/shar/articles/articles/background.html). - -``` r -torus_trans <- translate_raster(raster = landscape_classified) - -random_walk <- randomize_raster(raster = landscape_classified, n_random = 39) -``` - -To plot the randomized raster, you can use the plot function and specify -the number of raster as as well as the color palette used for the -discrete environmental data. - -``` r -col_palette <- c("#440154FF", "#3B528BFF", "#21908CFF", "#5DC863FF", "#FDE725FF") - -plot(torus_trans, n = 3, col = col_palette) -``` - -To randomize the point pattern, either use the Gamma test described by -Plotkin et al. (2000) or pattern reconstruction (Kirkpatrick et -al. 1983; Tscheschel & Stoyan 2006). - -``` r -gamma_test <- fit_point_process(pattern = species_b, process = "cluster", n_random = 39) +### Citation -# (this can takes some time) -reconstruction <- reconstruct_pattern(pattern = species_b, n_random = 39, e_threshold = 0.05) -``` - -Of course, there are several utility functions. For example, you can -plot the summary function of the observed pattern and the simulation -envelopes of randomized patterns (`what = "sf"`) or some randomized and -the observed pattern (`what = "pp"`) using the plot function. - -``` r -plot(reconstruction, what = "pp") -``` - - - -Another utility function allows to calculate the differences between the -observed pattern and the randomized patterns (also called energy using -summary functions). - -``` r -calculate_energy(reconstruction, return_mean = TRUE) -## [1] 0.04907375 -``` +The **shar** package is part of our academic work. To cite the package +or acknowledge its use in publications, please cite the following paper. -The data was created that `species_a` has a negative association to -habitat 4 and `species_b` has a positive association to habitat 5, which -is reflected in the results. - -Given the characteristics of the method, a positive association to one -habitat inevitably leads to a negative association to at least one of -the other habitats (and vice versa; Yamada et al. 2006). For example, a -high amount of individual points in the positively associated habitat -simultaneously mean that less individual points can be present in the -other habitats. - -Furthermore, please be aware that due to the randomization of the null -model data, results might slightly differ between different -randomization approaches (e.g., `fit_point_process()` -vs. `translate_raster()`) and even for repetitions of the same approach. -Thus, the exact `lo` and `hi` thresholds might be slightly different -when re-running the examples. However, the counts of the observed data -should be identical, and general results and trends should be similar. +> Hesselbarth, M.H.K., (2021). shar: A R package to analyze +> species-habitat associations using point pattern analysis. Journal of +> Open Source Software, 6(67), 3811. +> -``` r -significance_level <- 0.01 - -results_habitat_association(pattern = species_a, raster = torus_trans, significance_level = significance_level) -## > Input: randomized raster -## > Quantile thresholds: negative < 0.005 || positive > 0.995 -## habitat breaks count lo hi significance -## 1 1 NA 35 10 35 n.s. -## 2 2 NA 44 19 53 n.s. -## 3 3 NA 36 15 49 n.s. -## 4 4 NA 4 15 58 negative -## 5 5 NA 73 48 90 n.s. - -results_habitat_association(pattern = reconstruction, raster = landscape_classified, significance_level = significance_level) -## > Input: randomized pattern -## > Quantile thresholds: negative < 0.005 || positive > 0.995 -## habitat breaks count lo hi significance -## 1 1 NA 6 2.19 13.43 n.s. -## 2 2 NA 18 17.57 55.86 n.s. -## 3 3 NA 18 33.19 58.00 negative -## 4 4 NA 21 35.38 54.81 negative -## 5 5 NA 129 37.38 84.10 positive -``` +The get a BibTex entry, please use `citation("shar")`. ## Contributing and Code of Conduct @@ -197,39 +83,3 @@ you agree to abide by its terms. To contribute to this project, please see the [Contributing guidelines](CONTRIBUTING.md). - -### References - -Baddeley, A., Rubak, E., Turner, R., 2015. Spatial point patterns: -Methodology and applications with R. Chapman and Hall/CRC Press, London. -ISBN 978-1-4822-1020-0 - -Harms, K.E., Condit, R., Hubbell, S.P., Foster, R.B., 2001. Habitat -associations of trees and shrubs in a 50-ha neotropical forest plot. -Journal of Ecology 89, 947–959. - - -Hijmans, R.J., 2022 terra: Spatial Data Analysis. R package version -1.6-7. . - -Kirkpatrick, S., Gelatt, C.D.Jr., Vecchi, M.P., 1983. Optimization by -simulated annealing. Science 220, 671–680. - - -Plotkin, J.B., Potts, M.D., Leslie, N., Manokaran, N., LaFrankie, J.V., -Ashton, P.S., 2000. Species-area curves, spatial aggregation, and -habitat specialization in tropical forests. Journal of Theoretical -Biology 207, 81–99. - -Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random -point patterns. Computational Statistics and Data Analysis 51, 859–871. - - -Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern -analysis in ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN -978-1-4200-8254-8 - -Yamada, T., Tomita, A., Itoh, A., Yamakura, T., Ohkubo, T., Kanzaki, M., -Tan, S., Ashton, P.S., 2006. Habitat associations of Sterculiaceae trees -in a Bornean rain forest plot. Journal of Vegetation Science 17, -559–566. diff --git a/_pkgdown.yml b/_pkgdown.yml index db289d56..5a8c514a 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -10,10 +10,6 @@ reference: - species_a - species_b - landscape - - gamma_test - - reconstruction - - random_walk - - torus_trans - title: Randomize data contents: - fit_point_process @@ -25,16 +21,15 @@ reference: contents: - calculate_energy - plot_energy - - plot.rd_mar - - plot.rd_pat - - plot.rd_ras - results_habitat_association - title: Various contents: - classify_habitats - - estimate_pcf_fast - list_to_randomized - pack_randomized + - plot.rd_mar + - plot.rd_pat + - plot.rd_ras - print.rd_mar - print.rd_pat - print.rd_ras @@ -87,7 +82,7 @@ navbar: contact: icon: fa-user - hreft: https://mhesselbarth.rbind.io + hreft: https://www.maxhesselbarth.com home: sidebar: diff --git a/codemeta.json b/codemeta.json index 90099f22..09066327 100644 --- a/codemeta.json +++ b/codemeta.json @@ -7,7 +7,7 @@ "codeRepository": "https://r-spatialecology.github.io/shar/", "issueTracker": "https://github.com/r-spatialecology/shar/issues/", "license": "https://spdx.org/licenses/GPL-3.0", - "version": "2.0.4", + "version": "2.1", "programmingLanguage": { "@type": "ComputerLanguage", "name": "R", @@ -250,7 +250,7 @@ }, "SystemRequirements": null }, - "fileSize": "1007.274KB", + "fileSize": "723.822KB", "citation": [ { "@type": "ScholarlyArticle", diff --git a/cran-comments.md b/cran-comments.md index 6af778ff..17ecf01b 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,7 +1,10 @@ For details changes, please see NEWS.md +## shar 2.1 +Speed improvements + ## shar 2.0.4 -Small speed improvments +Small speed improvements ## shar 2.0.3 This is a re-submission. The below error has now been fixed diff --git a/data-raw/example_data.R b/data-raw/example_data.R index 84b73756..bb1bf53d 100644 --- a/data-raw/example_data.R +++ b/data-raw/example_data.R @@ -1,6 +1,7 @@ library(dplyr) library(NLMR) library(usethis) +library(rgbif) library(sf) library(spatstat) library(terra) @@ -75,9 +76,35 @@ gamma_test <- fit_point_process(pattern = species_b, n_random = n_random, proces reconstruction <- reconstruct_pattern(pattern = species_b, n_random = n_random, e_threshold = 0.05) + +#### Vignette Domestica + +#### gbif #### + +# Retrieve key for Cormus domestica +key <- rgbif::name_backbone(name = 'Cormus domestica', kingdom = 'plants') + +# Establish region of interest +roi <- c(xmin = -20, xmax = 45, ymin = 30, ymax = 73) +roi_bbox <- sf::st_bbox(roi, crs = sf::st_crs("EPSG:4326")) +roi_sfc <- sf::st_sfc(sf::st_point(c(roi[["xmin"]], roi[["ymin"]])), + sf::st_point(c(roi[["xmax"]], roi[["ymax"]])), + crs = "EPSG:4326") + +# Retrieve occurrences for the region of interest +# 99,999; 10,000 +res <- rgbif::occ_search(taxonKey = as.numeric(key$usageKey), limit = 99999, + geometry = roi_bbox) + +data_simp_precomp <- data.frame(id = res$data$key, lat = res$data$decimalLatitude, + lon = res$data$decimalLongitude) %>% + dplyr::filter(!is.na(lat) | !is.na(lon)) + +nrow(data_simp_precomp) + #### Save data #### -overwrite <- TRUE +overwrite <- FALSE # save landscape landscape <- terra::wrap(landscape) @@ -88,14 +115,10 @@ usethis::use_data(species_a, overwrite = overwrite) usethis::use_data(species_b, overwrite = overwrite) -# save random landscape data -torus_trans <- pack_randomized(raster = torus_trans) -usethis::use_data(torus_trans, overwrite = overwrite) - -random_walk <- pack_randomized(raster = random_walk) -usethis::use_data(random_walk, overwrite = overwrite) - -# save random point data -usethis::use_data(gamma_test, overwrite = overwrite) +# save interal data +data_internal <- list(torus_trans = pack_randomized(raster = torus_trans), + random_walk = pack_randomized(raster = random_walk), + gamma_test = gamma_test, reconstruction = reconstruction, + data_gbif = data_simp_precomp) -usethis::use_data(reconstruction, overwrite = overwrite) +usethis::use_data(data_internal, overwrite = overwrite, internal = TRUE) diff --git a/data-raw/vignette-domestica.R b/data-raw/vignette-domestica.R deleted file mode 100644 index 1171fcef..00000000 --- a/data-raw/vignette-domestica.R +++ /dev/null @@ -1,28 +0,0 @@ -library(dplyr) # For data wrangling -library(rgbif) # For retrieving species occurrence data -library(sf) # For spatial data operations - -#### gbif #### - -# Retrieve key for Cormus domestica -key <- rgbif::name_backbone(name = 'Cormus domestica', kingdom = 'plants') - -# Establish region of interest -roi <- c(xmin = -20, xmax = 45, ymin = 30, ymax = 73) -roi_bbox <- sf::st_bbox(roi, crs = sf::st_crs("EPSG:4326")) -roi_sfc <- sf::st_sfc(sf::st_point(c(roi[["xmin"]], roi[["ymin"]])), - sf::st_point(c(roi[["xmax"]], roi[["ymax"]])), - crs = "EPSG:4326") - -# Retrieve occurrences for the region of interest -# 99,999; 10,000 -res <- rgbif::occ_search(taxonKey = as.numeric(key$usageKey), limit = 99999, - geometry = roi_bbox) - -data_simp_precomp <- data.frame(id = res$data$key, lat = res$data$decimalLatitude, - lon = res$data$decimalLongitude) %>% - dplyr::filter(!is.na(lat) | !is.na(lon)) - -nrow(data_simp_precomp) - -usethis::use_data(data_simp_precomp, overwrite = TRUE, internal = TRUE) diff --git a/data/gamma_test.rda b/data/gamma_test.rda deleted file mode 100644 index 7ed24fc4..00000000 Binary files a/data/gamma_test.rda and /dev/null differ diff --git a/data/landscape.rda b/data/landscape.rda index 337c92b6..5eab6153 100644 Binary files a/data/landscape.rda and b/data/landscape.rda differ diff --git a/data/random_walk.rda b/data/random_walk.rda deleted file mode 100644 index 3f062f62..00000000 Binary files a/data/random_walk.rda and /dev/null differ diff --git a/data/reconstruction.rda b/data/reconstruction.rda deleted file mode 100644 index dfb6bc17..00000000 Binary files a/data/reconstruction.rda and /dev/null differ diff --git a/data/species_a.rda b/data/species_a.rda index d5b1db5e..f9f099d0 100644 Binary files a/data/species_a.rda and b/data/species_a.rda differ diff --git a/data/species_b.rda b/data/species_b.rda index 745e978d..97a8c14a 100644 Binary files a/data/species_b.rda and b/data/species_b.rda differ diff --git a/data/torus_trans.rda b/data/torus_trans.rda deleted file mode 100644 index c705f4ba..00000000 Binary files a/data/torus_trans.rda and /dev/null differ diff --git a/man/calc_gest.Rd b/man/calc_gest.Rd new file mode 100644 index 00000000..31ac2d72 --- /dev/null +++ b/man/calc_gest.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calc_gest.R +\name{calc_gest} +\alias{calc_gest} +\title{calc_gest} +\usage{ +calc_gest(dist, r, n_points) +} +\arguments{ +\item{dist}{matrix with distance pairs.} + +\item{r}{vector with distances r.} + +\item{n_points}{numeric with number of points} +} +\value{ +data.frame +} +\description{ +Calculate Gest +} +\details{ +Calculates Gest based on distances created with \code{get_dist_pairs}. +} +\seealso{ +\code{\link{get_dist_pairs}} +} +\keyword{internal} diff --git a/man/calculate_energy.Rd b/man/calculate_energy.Rd index 7471de57..b788f268 100644 --- a/man/calculate_energy.Rd +++ b/man/calculate_energy.Rd @@ -6,9 +6,8 @@ \usage{ calculate_energy( pattern, - weights = c(0.5, 0.5), + weights = c(1, 1), return_mean = FALSE, - comp_fast = 1000, verbose = TRUE ) } @@ -20,9 +19,6 @@ The first number refers to Gest(r), the second number to pcf(r).} \item{return_mean}{Logical if the mean energy is returned.} -\item{comp_fast}{Integer with threshold at which summary functions are estimated -in a computational fast way.} - \item{verbose}{Logical if progress report is printed.} } \value{ @@ -35,10 +31,7 @@ Calculate mean energy The function calculates the mean energy (or deviation) between the observed pattern and all reconstructed patterns (for more information see Tscheschel & Stoyan (2006) or Wiegand & Moloney (2014)). The pair correlation function and the -nearest neighbour distance function are used to describe the patterns. For large -patterns \code{comp_fast = TRUE} decreases the computational demand, because no edge -correction is used and the pair correlation function is estimated based on Ripley's -K-function. For more information see \code{\link{estimate_pcf_fast}}. +nearest neighbour distance function are used to describe the patterns. } \examples{ pattern_random <- fit_point_process(species_a, n_random = 19) diff --git a/man/estimate_pcf_fast.Rd b/man/estimate_pcf_fast.Rd deleted file mode 100644 index 767cf46d..00000000 --- a/man/estimate_pcf_fast.Rd +++ /dev/null @@ -1,45 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/estimate_pcf_fast.R -\name{estimate_pcf_fast} -\alias{estimate_pcf_fast} -\title{estimate_pcf_fast} -\usage{ -estimate_pcf_fast(pattern, ...) -} -\arguments{ -\item{pattern}{ppp object with point pattern.} - -\item{...}{Arguments passed down to \code{link{Kest}} or \code{\link{pcf.fv}}.} -} -\value{ -fv.object -} -\description{ -Fast estimation of the pair correlation function -} -\details{ -The functions estimates the pair correlation functions based on an estimation -of Ripley's K-function. This makes it computationally faster than estimating the -pair correlation function directly. It is a wrapper around \code{\link{Kest}} and -\code{\link{pcf.fv}}. -} -\examples{ -pcf_species_b <- estimate_pcf_fast(species_a) - -} -\references{ -Chiu, S.N., Stoyan, D., Kendall, W.S., Mecke, J., 2013. Stochastic geometry and -its applications, 3rd ed, Wiley Series in Probability and Statistics. -John Wiley & Sons Inc, Chichester, UK. ISBN 978-0-470-66481-0 - -Ripley, B.D., 1977. Modelling spatial patterns. Journal of the Royal Statistical -Society. Series B (Methodological) 39, 172–192. - - -Stoyan, D., Stoyan, H., 1994. Fractals, random shapes and point fields. -John Wiley & Sons, Chichester. ISBN 978-0-471-93757-9 -} -\seealso{ -\code{\link{Kest}} \cr -\code{\link{pcf.fv}} -} diff --git a/man/gamma_test.Rd b/man/gamma_test.Rd deleted file mode 100644 index e7f3dc18..00000000 --- a/man/gamma_test.Rd +++ /dev/null @@ -1,16 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data.R -\docType{data} -\name{gamma_test} -\alias{gamma_test} -\title{Gamma test} -\format{ -rd_pat object. -} -\usage{ -gamma_test -} -\description{ -Randomized data for species b using the gamma test. -} -\keyword{datasets} diff --git a/man/get_dist_pairs.Rd b/man/get_dist_pairs.Rd new file mode 100644 index 00000000..6c774fb9 --- /dev/null +++ b/man/get_dist_pairs.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_dist_pairs.R +\name{get_dist_pairs} +\alias{get_dist_pairs} +\title{get_dist_pairs} +\usage{ +get_dist_pairs(X, rmax) +} +\arguments{ +\item{X}{ppp object} + +\item{rmax}{Numeric with maximum distance} +} +\value{ +matrix +} +\description{ +Distance between points +} +\details{ +Returns matrix with point pairs and distances between them. +} +\seealso{ +\code{\link{pcf.ppp}} +} +\keyword{internal} diff --git a/man/plot.rd_mar.Rd b/man/plot.rd_mar.Rd index 126ac3d8..3b8dd405 100644 --- a/man/plot.rd_mar.Rd +++ b/man/plot.rd_mar.Rd @@ -9,7 +9,6 @@ what = "sf", n = NULL, probs = c(0.025, 0.975), - comp_fast = 1000, ask = TRUE, verbose = TRUE, ... @@ -26,9 +25,6 @@ See Details section for more information.} \item{probs}{Vector with quantiles of randomized data used for envelope construction.} -\item{comp_fast}{Integer with threshold at which summary functions are estimated -in a computational fast way.} - \item{ask}{Logical if the user is asked to press before second summary function is plotted (only used if \code{what = "sf"}).} @@ -45,9 +41,6 @@ Plot method for rd_pat object \details{ The function plots the pair correlation function and the nearest neighbour function of the observed pattern and the reconstructed patterns (as "simulation envelopes"). -For large patterns \code{comp_fast = TRUE} decreases the computational demand because no edge -correction is used and the pair correlation function is estimated based on Ripley's -K-function. For more information see \code{\link{estimate_pcf_fast}}. It is also possible to plot n randomized patterns and the observed pattern using \code{what = "pp"}. If \code{n} is a single number, \code{n} randomized diff --git a/man/plot.rd_pat.Rd b/man/plot.rd_pat.Rd index 2949e052..a6528638 100644 --- a/man/plot.rd_pat.Rd +++ b/man/plot.rd_pat.Rd @@ -9,7 +9,6 @@ what = "sf", n = NULL, probs = c(0.025, 0.975), - comp_fast = 1000, ask = TRUE, verbose = TRUE, ... @@ -26,9 +25,6 @@ See Details section for more information.} \item{probs}{Vector with quantiles of randomized data used for envelope construction.} -\item{comp_fast}{Integer with threshold at which summary functions are estimated -in a computational fast way.} - \item{ask}{Logical if the user is asked to press before second summary function is plotted (only used if \code{what = "sf"}).} @@ -45,9 +41,6 @@ Plot method for rd_pat object \details{ The function plots the pair correlation function and the nearest neighbour function of the observed pattern and the reconstructed patterns (as "simulation envelopes"). -For large patterns \code{comp_fast = TRUE} decreases the computational demand because no edge -correction is used and the pair correlation function is estimated based on Ripley's -K-function. For more information see \code{\link{estimate_pcf_fast}}. It is also possible to plot n randomized patterns and the observed pattern using \code{what = "pp"}. If \code{n} is a single number, \code{n} randomized diff --git a/man/random_walk.Rd b/man/random_walk.Rd deleted file mode 100644 index 34c66142..00000000 --- a/man/random_walk.Rd +++ /dev/null @@ -1,16 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data.R -\docType{data} -\name{random_walk} -\alias{random_walk} -\title{Random walk} -\format{ -rd_ras object. -} -\usage{ -random_walk -} -\description{ -Randomization of the \code{landscape} using the habitat randomization algorithm. -} -\keyword{datasets} diff --git a/man/reconstruct_algorithm.Rd b/man/reconstruct_algorithm.Rd new file mode 100644 index 00000000..6f2dc95b --- /dev/null +++ b/man/reconstruct_algorithm.Rd @@ -0,0 +1,60 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/reconstruct_algorithm.R +\name{reconstruct_algorithm} +\alias{reconstruct_algorithm} +\title{reconstruct_algorithm} +\usage{ +reconstruct_algorithm( + pattern, + method, + n_random, + e_threshold, + max_runs, + no_change, + annealing, + weights, + r_length, + r_max, + stoyan, + verbose, + plot +) +} +\arguments{ +\item{pattern}{ppp object with pattern.} + +\item{n_random}{Integer with number of randomizations.} + +\item{e_threshold}{Double with minimum energy to stop reconstruction.} + +\item{max_runs}{Integer with maximum number of iterations if \code{e_threshold} +is not reached.} + +\item{no_change}{Integer with number of iterations at which the reconstruction will +stop if the energy does not decrease.} + +\item{annealing}{Double with probability to keep relocated point even if energy +did not decrease.} + +\item{weights}{Vector with weights used to calculate energy. +The first number refers to Gest(r), the second number to pcf(r).} + +\item{r_length}{Integer with number of intervals from \code{r = 0} to \code{r = rmax} for which +the summary functions are evaluated.} + +\item{r_max}{Double with maximum distance used during calculation of summary functions. If \code{NULL}, +will be estimated from data.} + +\item{stoyan}{Coefficient for Stoyan's bandwidth selection rule.} + +\item{verbose}{Logical if progress report is printed.} + +\item{plot}{Logical if pcf(r) function is plotted and updated during optimization.} +} +\value{ +list +} +\description{ +Pattern reconstruction (internal) +} +\keyword{internal} diff --git a/man/reconstruct_pattern.Rd b/man/reconstruct_pattern.Rd index 5a449733..8237f862 100644 --- a/man/reconstruct_pattern.Rd +++ b/man/reconstruct_pattern.Rd @@ -9,15 +9,13 @@ reconstruct_pattern( method = "homo", n_random = 1, e_threshold = 0.01, - max_runs = 1000, + max_runs, no_change = Inf, annealing = 0.01, - comp_fast = 1000, - n_points = NULL, - window = NULL, - weights = c(0.5, 0.5), - r_length = 250, + weights = c(1, 1), + r_length = 255, r_max = NULL, + stoyan = 0.15, return_input = TRUE, simplify = FALSE, verbose = TRUE, @@ -43,13 +41,6 @@ stop if the energy does not decrease.} \item{annealing}{Double with probability to keep relocated point even if energy did not decrease.} -\item{comp_fast}{Integer with threshold at which summary functions are estimated -in a computational fast way.} - -\item{n_points}{Integer with number of points to be simulated.} - -\item{window}{owin object with window of simulated pattern.} - \item{weights}{Vector with weights used to calculate energy. The first number refers to Gest(r), the second number to pcf(r).} @@ -59,6 +50,8 @@ the summary functions are evaluated.} \item{r_max}{Double with maximum distance used during calculation of summary functions. If \code{NULL}, will be estimated from data.} +\item{stoyan}{Coefficient for Stoyan's bandwidth selection rule.} + \item{return_input}{Logical if the original input data is returned.} \item{simplify}{Logical if only pattern will be returned if \code{n_random=1} @@ -82,16 +75,12 @@ deviation between the observed and the reconstructed pattern decreases. The pair correlation function and the nearest neighbour distance function are used to describe the patterns. -For large patterns (\code{n > comp_fast}) the pair correlation function can be estimated -from Ripley's K-function without edge correction. This decreases the computational -time. For more information see \code{\link{estimate_pcf_fast}}. - The reconstruction can be stopped automatically if for n steps the energy does not decrease. The number of steps can be controlled by \code{no_change} and is set to \code{no_change = Inf} as default to never stop automatically. The weights must be 0 < sum(weights) <= 1. To weight both summary functions identical, -use \code{weights = c(0.5, 0.5)}. +use \code{weights = c(1, 1)}. \code{spatstat} sets \code{r_length} to 513 by default. However, a lower value decreases the computational time, while increasing the "bumpiness" of the summary function. diff --git a/man/reconstruct_pattern_cluster.Rd b/man/reconstruct_pattern_cluster.Rd deleted file mode 100644 index eff3f2e4..00000000 --- a/man/reconstruct_pattern_cluster.Rd +++ /dev/null @@ -1,84 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/reconstruct_pattern_cluster.R -\name{reconstruct_pattern_cluster} -\alias{reconstruct_pattern_cluster} -\title{reconstruct_pattern_cluster} -\usage{ -reconstruct_pattern_cluster( - pattern, - n_random = 1, - e_threshold = 0.01, - max_runs = 1000, - no_change = Inf, - annealing = 0.01, - comp_fast = 1000, - weights = c(0.5, 0.5), - r_length = 250, - r_max = NULL, - return_input = TRUE, - simplify = FALSE, - verbose = TRUE, - plot = FALSE -) -} -\arguments{ -\item{pattern}{ppp object with pattern.} - -\item{n_random}{Integer with number of randomizations.} - -\item{e_threshold}{Double with minimum energy to stop reconstruction.} - -\item{max_runs}{Integer with maximum number of iterations if \code{e_threshold} -is not reached.} - -\item{no_change}{Integer with number of iterations at which the reconstruction will -stop if the energy does not decrease.} - -\item{annealing}{Double with probability to keep relocated point even if energy -did not decrease.} - -\item{comp_fast}{Integer with threshold at which summary functions are estimated -in a computational fast way.} - -\item{weights}{Vector with weights used to calculate energy. -The first number refers to Gest(r), the second number to pcf(r).} - -\item{r_length}{Integer with number of intervals from \code{r = 0} to \code{r = rmax} for which -the summary functions are evaluated.} - -\item{r_max}{Double with maximum distance used during calculation of summary functions. If \code{NULL}, -will be estimated from data.} - -\item{return_input}{Logical if the original input data is returned.} - -\item{simplify}{Logical if only pattern will be returned if \code{n_random = 1} -and \code{return_input = FALSE}.} - -\item{verbose}{Logical if progress report is printed.} - -\item{plot}{Logical if pcf(r) function is plotted and updated during optimization.} -} -\value{ -rd_pat -} -\description{ -Pattern reconstruction for clustered patterns -} -\examples{ -\dontrun{ -pattern_recon <- reconstruct_pattern_cluster(species_b, n_random = 19, max_runs = 1000) -} - -} -\references{ -Kirkpatrick, S., Gelatt, C.D.Jr., Vecchi, M.P., 1983. Optimization by simulated -annealing. Science 220, 671–680. - -Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random point -patterns. Computational Statistics and Data Analysis 51, 859–871. - - -Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern analysis in -ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN 978-1-4200-8254-8 -} -\keyword{internal} diff --git a/man/reconstruct_pattern_hetero.Rd b/man/reconstruct_pattern_hetero.Rd deleted file mode 100644 index 1071ed2d..00000000 --- a/man/reconstruct_pattern_hetero.Rd +++ /dev/null @@ -1,86 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/reconstruct_pattern_hetero.R -\name{reconstruct_pattern_hetero} -\alias{reconstruct_pattern_hetero} -\title{reconstruct_pattern_hetero} -\usage{ -reconstruct_pattern_hetero( - pattern, - n_random = 1, - e_threshold = 0.01, - max_runs = 1000, - no_change = Inf, - annealing = 0.01, - comp_fast = 1000, - weights = c(0.5, 0.5), - r_length = 250, - r_max = NULL, - return_input = TRUE, - simplify = FALSE, - verbose = TRUE, - plot = FALSE -) -} -\arguments{ -\item{pattern}{ppp object with pattern.} - -\item{n_random}{Integer with number of randomizations.} - -\item{e_threshold}{Double with minimum energy to stop reconstruction.} - -\item{max_runs}{Integer with maximum number of iterations if \code{e_threshold} -is not reached.} - -\item{no_change}{Integer with number of iterations at which the reconstruction will -stop if the energy does not decrease.} - -\item{annealing}{Double with probability to keep relocated point even if energy -did not decrease.} - -\item{comp_fast}{Integer with threshold at which summary functions are estimated -in a computational fast way.} - -\item{weights}{Vector with weights used to calculate energy. -The first number refers to Gest(r), the second number to pcf(r).} - -\item{r_length}{Integer with number of intervals from \code{r = 0} to \code{r = rmax} for which -the summary functions are evaluated.} - -\item{r_max}{Double with maximum distance used during calculation of summary functions. If \code{NULL}, -will be estimated from data.} - -\item{return_input}{Logical if the original input data is returned.} - -\item{simplify}{Logical if only pattern will be returned if \code{n_random = 1} -and \code{return_input = FALSE}.} - -\item{verbose}{Logical if progress report is printed.} - -\item{plot}{Logical if pcf(r) function is plotted and updated during optimization.} -} -\value{ -rd_pat -} -\description{ -Pattern reconstruction for heterogeneous patterns -} -\examples{ -\dontrun{ -input_pattern <- spatstat.random::rpoispp(lambda = function(x, y) {100 * exp(-3 * x)}, -nsim = 1) -pattern_recon <- reconstruct_pattern_hetero(input_pattern, n_random = 19, max_runs = 1000) -} - -} -\references{ -Kirkpatrick, S., Gelatt, C.D.Jr., Vecchi, M.P., 1983. Optimization by simulated -annealing. Science 220, 671–680. - -Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random point -patterns. Computational Statistics and Data Analysis 51, 859–871. - - -Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern analysis in -ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN 978-1-4200-8254-8 -} -\keyword{internal} diff --git a/man/reconstruct_pattern_homo.Rd b/man/reconstruct_pattern_homo.Rd deleted file mode 100644 index ad7f5a7c..00000000 --- a/man/reconstruct_pattern_homo.Rd +++ /dev/null @@ -1,94 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/reconstruct_pattern_homo.R -\name{reconstruct_pattern_homo} -\alias{reconstruct_pattern_homo} -\title{reconstruct_pattern_homo} -\usage{ -reconstruct_pattern_homo( - pattern, - n_random = 1, - e_threshold = 0.01, - max_runs = 1000, - no_change = Inf, - annealing = 0.01, - n_points = NULL, - window = NULL, - comp_fast = 1000, - weights = c(0.5, 0.5), - r_length = 250, - r_max = NULL, - return_input = TRUE, - simplify = FALSE, - verbose = TRUE, - plot = FALSE -) -} -\arguments{ -\item{pattern}{ppp object with pattern.} - -\item{n_random}{Integer with number of randomizations.} - -\item{e_threshold}{Double with minimum energy to stop reconstruction.} - -\item{max_runs}{Integer with maximum number of iterations if \code{e_threshold} -is not reached.} - -\item{no_change}{Integer with number of iterations at which the reconstruction will -stop if the energy does not decrease.} - -\item{annealing}{Double with probability to keep relocated point even if energy -did not decrease.} - -\item{n_points}{Integer with number of points to be simulated.} - -\item{window}{owin object with window of simulated pattern.} - -\item{comp_fast}{Integer with threshold at which summary functions are estimated -in a computational fast way.} - -\item{weights}{Vector with weights used to calculate energy. -The first number refers to Gest(r), the second number to pcf(r).} - -\item{r_length}{Integer with number of intervals from \code{r = 0} to \code{r = rmax} for which -the summary functions are evaluated.} - -\item{r_max}{Double with maximum distance used during calculation of summary functions. If \code{NULL}, -will be estimated from data.} - -\item{return_input}{Logical if the original input data is returned.} - -\item{simplify}{Logical if only pattern will be returned if \code{n_random = 1} -and \code{return_input = FALSE}.} - -\item{verbose}{Logical if progress report is printed.} - -\item{plot}{Logical if pcf(r) function is plotted and updated during optimization.} -} -\value{ -rd_pat -} -\description{ -Pattern reconstruction for homogeneous pattern -} -\examples{ -\dontrun{ -pattern_recon_a <- reconstruct_pattern_homo(species_a, n_random = 19, -max_runs = 1000) - -pattern_recon_b <- reconstruct_pattern_homo(species_a, n_points = 70, -n_random = 19, max_runs = 1000) -} - -} -\references{ -Kirkpatrick, S., Gelatt, C.D.Jr., Vecchi, M.P., 1983. Optimization by simulated -annealing. Science 220, 671–680. - -Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random point -patterns. Computational Statistics and Data Analysis 51, 859–871. - - -Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern analysis in -ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN 978-1-4200-8254-8 -} -\keyword{internal} diff --git a/man/reconstruction.Rd b/man/reconstruction.Rd deleted file mode 100644 index c03e9f27..00000000 --- a/man/reconstruction.Rd +++ /dev/null @@ -1,16 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data.R -\docType{data} -\name{reconstruction} -\alias{reconstruction} -\title{Reconstruction} -\format{ -rd_pat object. -} -\usage{ -reconstruction -} -\description{ -Randomized data for species b using pattern reconstruction. -} -\keyword{datasets} diff --git a/man/torus_trans.Rd b/man/torus_trans.Rd deleted file mode 100644 index 636d3817..00000000 --- a/man/torus_trans.Rd +++ /dev/null @@ -1,16 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data.R -\docType{data} -\name{torus_trans} -\alias{torus_trans} -\title{Torus trans} -\format{ -rd_ras object. -} -\usage{ -torus_trans -} -\description{ -Torus translation of the classified \code{landscape}. -} -\keyword{datasets} diff --git a/tests/testthat/test-calculate_energy.R b/tests/testthat/test-calculate_energy.R index ef8d2020..e0333657 100644 --- a/tests/testthat/test-calculate_energy.R +++ b/tests/testthat/test-calculate_energy.R @@ -10,9 +10,8 @@ marks_sub <- spatstat.geom::subset.ppp(species_a, select = dbh) marks_recon <- reconstruct_pattern_marks(pattern_random_a$randomized[[1]], marks_sub, n_random = 3, max_runs = 10, verbose = FALSE) -marks_recon_NA <- marks_recon - -marks_recon_NA$energy_df <- "NA" +marks_recon_na <- marks_recon +marks_recon_na$energy_df <- "NA" ################################################################################ @@ -31,7 +30,6 @@ testthat::test_that("calculate_energy uses weights", { testthat::expect_false(unweighted == weighted) }) - testthat::test_that("calculate_energy returns mean ", { mean_energy <- mean(calculate_energy(pattern_random_a, verbose = FALSE)) @@ -41,13 +39,6 @@ testthat::test_that("calculate_energy returns mean ", { expected = mean_energy) }) -testthat::test_that("calculate_energy can use comp_fast ", { - - testthat::expect_length(calculate_energy(pattern_random_a, comp_fast = 50, - verbose = FALSE), - n = 3) -}) - testthat::test_that("calculate_energy returns works for reconstructed marks", { testthat::expect_length(calculate_energy(marks_recon, verbose = FALSE), n = 3) @@ -57,8 +48,8 @@ testthat::test_that("calculate_energy returns works for reconstructed marks", { testthat::expect_equal(calculate_energy(marks_recon, return_mean = TRUE, verbose = FALSE), expected = mean_energy) - testthat::expect_equal(calculate_energy(marks_recon_NA, return_mean = TRUE, - verbose = FALSE), expected = mean_energy) + testthat::expect_length(calculate_energy(marks_recon_na, verbose = FALSE), n = 3) + }) testthat::test_that("calculate_energy returns error if observed not included", { @@ -74,17 +65,3 @@ testthat::test_that("calculate_energy returns error if wrong class ", { regexp = "Class of 'pattern' must be 'rd_pat' or 'rd_mar'.", fixed = TRUE) }) - -testthat::test_that("calculate_energy returns error if weights are wrong ", { - - testthat::expect_error(calculate_energy(pattern_random_a, weights = c(0, 0), - return_mean = TRUE, verbose = FALSE), - regexp = "The sum of 'weights' must be 0 < sum(weights) <= 1.", - fixed = TRUE) - - testthat::expect_error(calculate_energy(pattern_random_a, weights = c(1, 1), - return_mean = TRUE, verbose = FALSE), - regexp = "The sum of 'weights' must be 0 < sum(weights) <= 1.", - fixed = TRUE) -}) - diff --git a/tests/testthat/test-estimate_pcf_fast.R b/tests/testthat/test-estimate_pcf_fast.R deleted file mode 100644 index f6f61b74..00000000 --- a/tests/testthat/test-estimate_pcf_fast.R +++ /dev/null @@ -1,10 +0,0 @@ -testthat::context("test-estimate_pcf_fast") - -################################################################################ - -testthat::test_that("estimate_pcf returns spatstat.fv object", { - - pcf_est <- estimate_pcf_fast(pattern = species_b) - - testthat::expect_is(pcf_est, "fv") -}) diff --git a/tests/testthat/test-pack-unpack.R b/tests/testthat/test-pack-unpack.R new file mode 100644 index 00000000..86b605bf --- /dev/null +++ b/tests/testthat/test-pack-unpack.R @@ -0,0 +1,37 @@ +testthat::context("test-pack_randomized") + +landscape_classified <- classify_habitats(terra::rast(landscape), n = 5, style = "fisher") +landscape_classified[terra::values(landscape_classified) != 1] <- 2 + +landscape_random <- randomize_raster(landscape_classified, n_random = 2) +landscape_ni <- randomize_raster(landscape_classified, n_random = 2, return_input = FALSE) + +x <- pack_randomized(raster = landscape_random) +x_ni <- pack_randomized(raster = landscape_random) + + +################################################################################ + +testthat::test_that("pack_randomized wraps raster", { + + testthat::expect_s4_class(object = x$observed, class = "PackedSpatRaster") + testthat::expect_true(all(sapply(x$randomized, inherits, what = "PackedSpatRaster"))) + testthat::expect_true(all(sapply(x_ni$randomized, inherits, what = "PackedSpatRaster"))) + +}) + +testthat::context("test-pack_randomized") + +y <- unpack_randomized(raster = x) +y_ni <- unpack_randomized(raster = x_ni) + + +testthat::test_that("unpack_randomized unwraps raster", { + + testthat::expect_s4_class(object = y$observed, class = "SpatRaster") + testthat::expect_true(all(sapply(y$randomized, inherits, what = "SpatRaster"))) + testthat::expect_true(all(sapply(y_ni$randomized, inherits, what = "SpatRaster"))) + +}) + + diff --git a/tests/testthat/test-plot_rd_pat.R b/tests/testthat/test-plot_rd_pat.R index 64f22a9d..e9fa8d8b 100644 --- a/tests/testthat/test-plot_rd_pat.R +++ b/tests/testthat/test-plot_rd_pat.R @@ -30,14 +30,10 @@ testthat::test_that("plot returns error if observed is missing", { fixed = TRUE) }) -testthat::test_that("plot uses comp_fast", { - - testthat::expect_null(plot(pattern_random, comp_fast = 50, verbose = FALSE, ask = FALSE)) -}) - testthat::test_that("plot works for reconstructed marks", { testthat::expect_null(plot(marks_recon, verbose = FALSE)) + testthat::expect_null(plot(marks_recon, what = "pp", verbose = FALSE)) }) testthat::test_that("plot returns error if what is wrong", { diff --git a/tests/testthat/test-reconstruct_pattern.R b/tests/testthat/test-reconstruct_pattern.R index ca32dd1c..d4302c24 100644 --- a/tests/testthat/test-reconstruct_pattern.R +++ b/tests/testthat/test-reconstruct_pattern.R @@ -1,7 +1,7 @@ testthat::context("test-reconstruct_pattern") # normal reconstruction -pattern_recon_homo <- reconstruct_pattern(pattern = species_a, n_random = 3, method = "homo", +pattern_recon_homo <- reconstruct_pattern(pattern = species_a, n_random = 3, max_runs = 1, verbose = FALSE) # cluster reconstruction @@ -10,7 +10,19 @@ pattern_recon_cluster <- reconstruct_pattern(pattern = species_a, n_random = 3, # cluster reconstruction pattern_recon_hetero <- reconstruct_pattern(pattern = species_b, n_random = 3, method = "hetero", - max_runs = 1, verbose = FALSE) + max_runs = 1, verbose = FALSE) + +pattern_recon_ni <- reconstruct_pattern(pattern = species_a, n_random = 3, + max_runs = 1, return_input = FALSE, + verbose = FALSE) + +pattern_recon_energy <- reconstruct_pattern(pattern = species_a, max_runs = 1000, + e_threshold = 0.1, n_random = 3, + verbose = FALSE) + +pattern_recon_simple <- reconstruct_pattern(pattern = species_a, n_random = 1, + max_runs = 1, simplify = TRUE, return_input = FALSE, + verbose = FALSE) ################################################################################ @@ -24,17 +36,66 @@ testthat::test_that("reconstruct_pattern returns correct class", { }) -testthat::test_that("reconstruct_patternreturns warnings", { +testthat::test_that("Output is a long as n_random for reconstruct_pattern", { + + testthat::expect_type(pattern_recon_homo$randomized, type = "list") + + testthat::expect_length(pattern_recon_homo$randomized, n = 3) +}) + +testthat::test_that("Output includes randomizations and original pattern for reconstruct_pattern", { + + testthat::expect_named(pattern_recon_homo$randomized, + expected = paste0("randomized_", c(1:3))) + + testthat::expect_equal(pattern_recon_homo$observed, + expected = spatstat.geom::unmark(species_a)) +}) + +testthat::test_that("Reconstructed patterns have same number of points", { + + testthat::expect_true(all(vapply(pattern_recon_homo$randomized, + FUN.VALUE = logical(1), + function(x) x$n == species_a$n))) +}) + +testthat::test_that("Input pattern can not be returned for reconstruct_pattern", { + + testthat::expect_equal(object = pattern_recon_ni$observed, + expected = "NA") +}) + +testthat::test_that("Reconstruction stops if e_threshold is reached", { + + energy <- calculate_energy(pattern_recon_energy, verbose = FALSE) + + testthat::expect_true(object = all(energy < 0.1)) + + testthat::expect_true(all(pattern_recon_energy$stop_criterion == "e_threshold")) + +}) + +testthat::test_that("simplify works for reconstruct_pattern", { - testthat::expect_warning(reconstruct_pattern(pattern = species_a, n_random = 3, - method = "cluster", n_points = 5, - max_runs = 1), + testthat::expect_is(pattern_recon_simple, "ppp") +}) + +testthat::test_that("reconstruct_pattern returns error if n_random < 1", { + + testthat::expect_error(reconstruct_pattern(pattern = species_a, n_random = -5, + verbose = FALSE), + regexp = "n_random must be >= 1.") +}) - regexp = "'n_points', 'window', or 'r_max' are not used for method='cluster'.") +testthat::test_that("reconstruct_pattern returns warnings", { - testthat::expect_warning(reconstruct_pattern(pattern = species_a, n_random = 3, - method = "hetero", n_points = 5, - max_runs = 1), + testthat::expect_warning(reconstruct_pattern(pattern = species_a, + n_random = 2, max_runs = 1, + return_input = FALSE, simplify = TRUE), + regexp = "'simplify = TRUE' not possible for 'n_random > 1'.") - regexp = "'n_points', 'window', or 'r_max' are not used for method='hetero'.") + testthat::expect_warning(reconstruct_pattern(pattern = species_a, + n_random = 1, max_runs = 1, + simplify = TRUE), + regexp = "'simplify = TRUE' not possible for 'return_input = TRUE'.") }) diff --git a/tests/testthat/test-reconstruct_pattern_cluster.R b/tests/testthat/test-reconstruct_pattern_cluster.R deleted file mode 100644 index 90b1a2b4..00000000 --- a/tests/testthat/test-reconstruct_pattern_cluster.R +++ /dev/null @@ -1,109 +0,0 @@ -testthat::context("test-reconstruct_pattern_cluster") - -# normal reconstruction -pattern_recon <- reconstruct_pattern_cluster(pattern = species_b, n_random = 3, - max_runs = 1, verbose = FALSE) - -pattern_recon_ni <- reconstruct_pattern_cluster(pattern = species_b, n_random = 2, - max_runs = 1, return_input = FALSE, - verbose = FALSE) - -pattern_recon_comp_fast <- reconstruct_pattern_cluster(pattern = species_b, - n_random = 1, max_runs = 1, - comp_fast = 0, verbose = FALSE) - -pattern_recon_energy <- reconstruct_pattern_cluster(pattern = species_b, - e_threshold = 0.1, n_random = 3, - verbose = FALSE) - -pattern_recon_simple <- reconstruct_pattern_cluster(pattern = species_b, - n_random = 1, max_runs = 1, - return_input = FALSE, simplify = TRUE, - verbose = FALSE) - -################################################################################ - -testthat::test_that("Output is a long as n_random for reconstruct_pattern_cluster", { - - testthat::expect_is(pattern_recon, class = "rd_pat") - - testthat::expect_type(pattern_recon$randomized, type = "list") - - testthat::expect_length(pattern_recon$randomized, n = 3) -}) - -testthat::test_that("Output includes randomizations and original pattern for reconstruct_pattern_cluster", { - - testthat::expect_named(pattern_recon$randomized, - expected = paste0("randomized_", c(1:3))) - - testthat::expect_equal(pattern_recon$observed, - expected = spatstat.geom::unmark(species_b)) -}) - -testthat::test_that("Reconstructed patterns have same number of points", { - - testthat::expect_true(all(vapply(pattern_recon$randomized, - FUN.VALUE = logical(1), - function(x) x$n == species_b$n))) -}) - -testthat::test_that("Input pattern can not be returned for reconstruct_pattern_cluster", { - - testthat::expect_equal(object = pattern_recon_ni$observed, - expected = "NA") -}) - -testthat::test_that("Argument comp_fast = TRUE is working", { - - testthat::expect_is(pattern_recon_comp_fast, "rd_pat") -}) - -testthat::test_that("Reconstruction stops if e_threshold is reached", { - - energy <- calculate_energy(pattern_recon_energy, verbose = FALSE) - - testthat::expect_true(object = all(energy < 0.1)) -}) - - -testthat::test_that("simplify works for reconstruct_pattern_cluster", { - - testthat::expect_is(pattern_recon_simple, "ppp") -}) - -testthat::test_that("reconstruct_pattern_cluster returns error if n_random < 1", { - - testthat::expect_error(reconstruct_pattern_cluster(pattern = species_b, - n_random = -5, verbose = FALSE), - regexp = "n_random must be >= 1.", - fixed = TRUE) -}) - -testthat::test_that("reconstruct_pattern_cluster returns error if weights are wrong ", { - - testthat::expect_error(reconstruct_pattern_cluster(pattern = species_b, - weights = c(0, 0), verbose = FALSE), - regexp = "The sum of 'weights' must be 0 < sum(weights) <= 1.", - fixed = TRUE) - - testthat::expect_error(reconstruct_pattern_cluster(pattern = species_b, - weights = c(1, 1), verbose = FALSE), - regexp = "The sum of 'weights' must be 0 < sum(weights) <= 1.", - fixed = TRUE) -}) - -testthat::test_that("reconstruct_pattern_cluster returns warnings", { - - testthat::expect_warning(reconstruct_pattern_cluster(pattern = species_b, - n_random = 2, max_runs = 1, - return_input = FALSE, simplify = TRUE), - regexp = "'simplify = TRUE' not possible for 'n_random > 1'.", - fixed = TRUE) - - testthat::expect_warning(reconstruct_pattern_cluster(pattern = species_b, - n_random = 1, max_runs = 1, - simplify = TRUE), - regexp = "'simplify = TRUE' not possible for 'return_input = TRUE'.", - fixed = TRUE) -}) diff --git a/tests/testthat/test-reconstruct_pattern_hetero.R b/tests/testthat/test-reconstruct_pattern_hetero.R deleted file mode 100644 index 05352aff..00000000 --- a/tests/testthat/test-reconstruct_pattern_hetero.R +++ /dev/null @@ -1,112 +0,0 @@ -testthat::context("test-reconstruct_pattern_hetero") - -input_pattern <- spatstat.random::rpoispp(lambda = function(x, y) {500 * exp(-3 * x)}, nsim = 1) - -# normal reconstruction -pattern_recon <- reconstruct_pattern_hetero(pattern = input_pattern, - n_random = 3, max_runs = 1, - verbose = FALSE) - -pattern_recon_ni <- reconstruct_pattern_hetero(pattern = input_pattern, - n_random = 2, max_runs = 1, - return_input = FALSE, verbose = FALSE) - -pattern_recon_comp_fast <- reconstruct_pattern_hetero(pattern = input_pattern, - n_random = 1, max_runs = 1, - comp_fast = 0, verbose = FALSE) - -pattern_recon_energy <- reconstruct_pattern_hetero(pattern = input_pattern, - e_threshold = 0.1, n_random = 3, - verbose = FALSE) - -pattern_recon_simple <- reconstruct_pattern_hetero(pattern = input_pattern, - n_random = 1, max_runs = 1, - return_input = FALSE, simplify = TRUE, - verbose = FALSE) - -################################################################################ - -testthat::test_that("Output is a long as n_random for reconstruct_pattern_hetero", { - - testthat::expect_is(pattern_recon, class = "rd_pat") - - testthat::expect_type(pattern_recon$randomized, type = "list") - - testthat::expect_length(pattern_recon$randomized, n = 3) -}) - -testthat::test_that("Output includes randomizations and original pattern for reconstruct_pattern_hetero", { - - testthat::expect_named(pattern_recon$randomized, - expected = paste0("randomized_", c(1:3))) - - testthat::expect_equal(pattern_recon$observed, - expected = spatstat.geom::unmark(input_pattern)) -}) - -testthat::test_that("Reconstructed patterns have same number of points", { - - testthat::expect_true(all(vapply(pattern_recon$randomized, - FUN.VALUE = logical(1), - function(x) x$n == input_pattern$n))) -}) - -testthat::test_that("Input pattern can not be returned for reconstruct_pattern_hetero", { - - testthat::expect_equal(object = pattern_recon_ni$observed, - expected = "NA") -}) - -testthat::test_that("Argument comp_fast = TRUE is working", { - - testthat::expect_is(pattern_recon_comp_fast, "rd_pat") -}) - -testthat::test_that("Reconstruction stops if e_threshold is reached", { - - energy <- calculate_energy(pattern_recon_energy, verbose = FALSE) - - testthat::expect_true(object = all(energy < 0.1)) -}) - - -testthat::test_that("simplify works for reconstruct_pattern_hetero", { - - testthat::expect_is(pattern_recon_simple, "ppp") -}) - -testthat::test_that("reconstruct_pattern_hetero returns error if n_random < 1", { - - testthat::expect_error(reconstruct_pattern_hetero(pattern = input_pattern, - n_random = -5, verbose = FALSE), - regexp = "n_random must be >= 1.", - fixed = TRUE) -}) - -testthat::test_that("reconstruct_pattern_hetero returns error if weights are wrong ", { - - testthat::expect_error(reconstruct_pattern_hetero(pattern = input_pattern, - weights = c(0, 0), verbose = FALSE), - regexp = "The sum of 'weights' must be 0 < sum(weights) <= 1.", - fixed = TRUE) - - testthat::expect_error(reconstruct_pattern_hetero(pattern = input_pattern, - weights = c(1, 1), verbose = FALSE), - regexp = "The sum of 'weights' must be 0 < sum(weights) <= 1.", - fixed = TRUE) -}) - -testthat::test_that("reconstruct_pattern_hetero returns warnings", { - - testthat::expect_warning(reconstruct_pattern_hetero(pattern = input_pattern, - n_random = 2, max_runs = 1, - return_input = FALSE, simplify = TRUE), - regexp = "'simplify = TRUE' not possible for 'n_random > 1'.", - fixed = TRUE) - - testthat::expect_warning(reconstruct_pattern_hetero(pattern = input_pattern, - n_random = 1, max_runs = 1, - simplify = TRUE), - regexp = "'simplify = TRUE' not possible for 'return_input = TRUE'.", - fixed = TRUE) -}) diff --git a/tests/testthat/test-reconstruct_pattern_homo.R b/tests/testthat/test-reconstruct_pattern_homo.R deleted file mode 100644 index 387d34d6..00000000 --- a/tests/testthat/test-reconstruct_pattern_homo.R +++ /dev/null @@ -1,132 +0,0 @@ -testthat::context("test-reconstruct_pattern_homo") - -# normal reconstruction -pattern_recon <- reconstruct_pattern_homo(pattern = species_a, n_random = 3, - max_runs = 1, verbose = FALSE) - -pattern_recon_ni <- reconstruct_pattern_homo(pattern = species_a, n_random = 2, - max_runs = 1, return_input = FALSE, - verbose = FALSE) - -pattern_recon_comp_fast <- reconstruct_pattern_homo(pattern = species_a, - n_random = 1, max_runs = 1, - comp_fast = 0, verbose = FALSE) - -pattern_recon_energy <- reconstruct_pattern_homo(pattern = species_a, - e_threshold = 0.1, n_random = 3, - verbose = FALSE) - -pattern_recon_simple <- reconstruct_pattern_homo(pattern = species_a, - n_random = 1, max_runs = 1, - return_input = FALSE, simplify = TRUE, - verbose = FALSE) - -n_points_diff <- 100 -window_diff <- spatstat.geom::owin(xrange = c(0, 1100), yrange = c(0, 1250)) - -pattern_recon_diff <- reconstruct_pattern_homo(pattern = species_a, - n_points = n_points_diff, - window = window_diff, r_max = 350, - n_random = 3, max_runs = 1, - verbose = FALSE) - -################################################################################ - -testthat::test_that("Output is a long as n_random for reconstruct_pattern_homo", { - - testthat::expect_is(pattern_recon, class = "rd_pat") - - testthat::expect_type(pattern_recon$randomized, type = "list") - - testthat::expect_length(pattern_recon$randomized, n = 3) -}) - -testthat::test_that("Output includes randomizations and original pattern for reconstruct_pattern_homo", { - - testthat::expect_named(pattern_recon$randomized, - expected = paste0("randomized_", c(1:3))) - - testthat::expect_equal(pattern_recon$observed, - expected = spatstat.geom::unmark(species_a)) -}) - -testthat::test_that("Reconstructed patterns have same number of points", { - - testthat::expect_true(all(vapply(pattern_recon$randomized, - FUN.VALUE = logical(1), - function(x) x$n == species_a$n))) -}) - -testthat::test_that("Reconstructed patterns have specified number of points and win", { - - testthat::expect_true(all(vapply(pattern_recon_diff$randomized, - FUN.VALUE = logical(1), - function(x) x$n == n_points_diff))) - - testthat::expect_true(all(vapply(pattern_recon_diff$randomized, - FUN.VALUE = logical(1), - function(x) all(x$window$xrange == window_diff$xrange)))) - - testthat::expect_true(all(vapply(pattern_recon_diff$randomized, - FUN.VALUE = logical(1), - function(x) all(x$window$yrange == window_diff$yrange)))) -}) - -testthat::test_that("Input pattern can not be returned for reconstruct_pattern_homo", { - - testthat::expect_equal(object = pattern_recon_ni$observed, - expected = "NA") -}) - -testthat::test_that("Argument comp_fast = TRUE is working", { - - testthat::expect_is(pattern_recon_comp_fast, "rd_pat") -}) - -testthat::test_that("Reconstruction stops if e_threshold is reached", { - - energy <- calculate_energy(pattern_recon_energy, verbose = FALSE) - - testthat::expect_true(object = all(energy < 0.1)) -}) - -testthat::test_that("simplify works for reconstruct_pattern_homo", { - - testthat::expect_is(pattern_recon_simple, "ppp") -}) - -testthat::test_that("reconstruct_pattern_homo returns error if n_random < 1", { - - testthat::expect_error(reconstruct_pattern_homo(pattern = species_a, n_random = -5, - verbose = FALSE), - regexp = "n_random must be >= 1.", - fixed = TRUE) -}) - -testthat::test_that("reconstruct_pattern_homo returns error if weights are wrong ", { - - testthat::expect_error(reconstruct_pattern_homo(pattern = species_a, weights = c(0, 0), - verbose = FALSE), - regexp = "The sum of 'weights' must be 0 < sum(weights) <= 1.", - fixed = TRUE) - - testthat::expect_error(reconstruct_pattern_homo(pattern = species_a, weights = c(1, 1), - verbose = FALSE), - regexp = "The sum of 'weights' must be 0 < sum(weights) <= 1.", - fixed = TRUE) -}) - -testthat::test_that("reconstruct_pattern_homo returns warnings", { - - testthat::expect_warning(reconstruct_pattern_homo(pattern = species_a, - n_random = 2, max_runs = 1, - return_input = FALSE, simplify = TRUE), - regexp = "'simplify = TRUE' not possible for 'n_random > 1'.", - fixed = TRUE) - - testthat::expect_warning(reconstruct_pattern_homo(pattern = species_a, - n_random = 1, max_runs = 1, - simplify = TRUE), - regexp = "'simplify = TRUE' not possible for 'return_input = TRUE'.", - fixed = TRUE) -}) diff --git a/tests/testthat/test-reconstruct_pattern_marks.R b/tests/testthat/test-reconstruct_pattern_marks.R index 6c060e15..b1e926f6 100644 --- a/tests/testthat/test-reconstruct_pattern_marks.R +++ b/tests/testthat/test-reconstruct_pattern_marks.R @@ -1,8 +1,8 @@ testthat::context("test-reconstruct_pattern_marks") -pattern_recon <- reconstruct_pattern_homo(species_a, n_random = 1, return_input = FALSE, - simplify = TRUE, max_runs = 1, - verbose = FALSE) +pattern_recon <- reconstruct_pattern(species_a, n_random = 1, return_input = FALSE, + simplify = TRUE, max_runs = 1, + verbose = FALSE) marks_sub <- spatstat.geom::subset.ppp(species_a, select = dbh) @@ -25,20 +25,6 @@ marks_recon_energy <- reconstruct_pattern_marks(pattern = pattern_recon, marked_ n_random = 3, e_threshold = 0.1, verbose = FALSE) - -n_points_diff <- 100 -window_diff <- spatstat.geom::owin(xrange = c(0, 900), yrange = c(0, 1250)) - -pattern_recon_diff <- reconstruct_pattern_homo(species_a, n_points = n_points_diff, - window = window_diff, n_random = 1, - return_input = FALSE, simplify = TRUE, - max_runs = 1, verbose = FALSE) - -marks_recon_diff <- reconstruct_pattern_marks(pattern = pattern_recon_diff, - marked_pattern = marks_sub, - n_random = 3, max_runs = 1, - verbose = FALSE) - ################################################################################ testthat::test_that("Output is a long as n_random for reconstruct_pattern_marks", { @@ -50,21 +36,6 @@ testthat::test_that("Output is a long as n_random for reconstruct_pattern_marks" testthat::expect_length(marks_recon$randomized, n = 3) }) - -testthat::test_that("reconstruct_pattern_marks works for pattern of difference length", { - - testthat::expect_is(marks_recon_diff, class = "rd_mar") - - testthat::expect_type(marks_recon_diff$randomized, type = "list") - - testthat::expect_length(marks_recon_diff$randomized, n = 3) - - testthat::expect_named(marks_recon_diff$randomized, - expected = paste0("randomized_", c(1:3))) - - testthat::expect_equal(marks_recon_diff$observed, expected = marks_sub) -}) - testthat::test_that("Output includes randomizations and original pattern for reconstruct_pattern_marks", { testthat::expect_named(marks_recon$randomized, @@ -90,7 +61,7 @@ testthat::test_that("Reconstruction stops if e_threshold is reached", { testthat::expect_true(all(energy < 0.1 & energy > 0.01)) - testthat::expect_true(all(marks_recon_energy$stop_criterion == "e_threshold/no_change")) + testthat::expect_true(all(marks_recon_energy$stop_criterion %in% c("e_threshold", "no_change"))) }) testthat::test_that("All errors are returned for reconstruct_pattern_marks", { @@ -116,14 +87,6 @@ testthat::test_that("All errors are returned for reconstruct_pattern_marks", { regexp = "'pattern' must be unmarked and 'marked_pattern' marked", fixed = TRUE) - # testthat::expect_error(reconstruct_pattern_marks(pattern = spatstat.geom::unmark(species_b), - # marked_pattern = marks_sub, - # n_random = 3, - # max_runs = 1, - # verbose = FALSE), - # regexp = "'pattern' and 'pattern' must have same window and number of points", - # fixed = TRUE) - testthat::expect_error(reconstruct_pattern_marks(pattern = pattern_recon, marked_pattern = spatstat.geom::subset.ppp(species_a, select = status), diff --git a/vignettes/articles/cormus_domestica_tmp.Rmd b/vignettes/articles/cormus_domestica_tmp.Rmd index 4e4d21b2..7b33b364 100644 --- a/vignettes/articles/cormus_domestica_tmp.Rmd +++ b/vignettes/articles/cormus_domestica_tmp.Rmd @@ -62,7 +62,7 @@ To retrieve species occurrence data the `R` package `rgbif` (Chamberlain & Boett ```{r gbif, echo=FALSE, eval=TRUE, message=FALSE, warning=FALSE} # Load bundled GBIF occurrence data -data_simp <- shar:::data_simp_precomp +data_simp <- shar:::data_internal$data_gbif ``` ```{r roi, echo=FALSE, eval=TRUE, message=FALSE, warning=FALSE} diff --git a/vignettes/get_started.Rmd b/vignettes/get_started.Rmd index 9ecb6007..962a4520 100644 --- a/vignettes/get_started.Rmd +++ b/vignettes/get_started.Rmd @@ -28,11 +28,15 @@ set.seed(42) ## Design -The core of **shar** are functions to to simulate null model data by randomizing either the environmental data (i.e. raster data) or the locations of species (i.e. point pattern data). The null data is then used to analyse if significant species-habitat associations are present. Additionally, functions to visualize and analyse the results are available as well as some utility functions. The methods are mainly described in Harms et al. (2001), Plotkin et al. (2000) and Wiegand & Moloney (2014). The methods are not necessary complementary, but are rather different approaches for the same result. +The core of **shar** are functions to simulate null model data by randomizing either the environmental data (i.e., raster) or the locations of species (i.e., point pattern). The null model data is then used to analyse if significant species-habitat associations are present. Additionally, functions to visualize and analyse the results are available as well as some utility functions. The methods are mainly described in Harms et al. (2001), Plotkin et al. (2000), and Wiegand & Moloney (2014). The methods are not necessary complementary, but are rather different approaches to achieve the same result. + +### Data + +**shar** comes with build-in example data sets. `species_a` and `species_b` are exemplary locations of individuals, e.g., tree locations wthin a forest study plot (as `ppp`-objects); `landscape` contains exemplary continuous environmental data (as `SpatRast`-object). ## Preprocessing of input data -All functions are designed for discrete habitat classes. Following, in case only continuous data is available, this has to be classified to discrete classes. `classify_habitats` provides several ways to classify the data such as the Fisher-Jenks algorithm (Fisher 1958, Jenks & Caspall 1971) +All functions are designed for discrete habitat classes. Thus, if continuous data is available, this has to be classified to discrete classes first. `classify_habitats` provides several ways to classify the data, such as the Fisher-Jenks algorithm (Fisher 1958, Jenks & Caspall 1971). ```{r classify_habitats} landscape_discrete <- classify_habitats(raster = terra::rast(landscape), n = 5, style = "fisher") @@ -40,9 +44,9 @@ landscape_discrete <- classify_habitats(raster = terra::rast(landscape), n = 5, ## Randomize environmental data -There are two functions to randomize the environmental data: `translate_raster()` and `randomize_raster()`. The first function is a torus translation of the raster, shifting the habitat map in all four cardinal directions. This is only possible for rectangular observation areas and results in `n_random <- (terra::nrow(terra::rast(landscape)) + 1) * (terra::ncol(terra::rast(landscape)) + 1) - 4` randomized rasters. The other function randomizes the environmental data using a random-walk algorithm. Here, the number of randomized rasters can be specified using the `n_random` argument. +There are two functions to randomize the environmental data: `translate_raster()` and `randomize_raster()`. The first function is a torus translation of the raster, shifting the habitat map in all four cardinal directions. This is only possible for rectangular observation areas and results in `R (terra::nrow(terra::rast(landscape)) + 1) * (terra::ncol(terra::rast(landscape)) + 1) - 4` randomized raster (based on number of rows and cols). The other function randomizes the environmental data using a random-walk algorithm. Here, the number of randomized raster can be specified using the `n_random` argument. -However, all methods require "fully mapped data" in a sense that NA cells of the environmental data are allowed only if simultaneously these areas cannot accommodate any locations of the point pattern (e.g., a water body within a forest area). This needs to be reflected in the observation window of the point pattern. For the torus translation method, no NA values are allowed at all. +All methods require "fully mapped data" in a sense that `NA` cells of the environmental data are allowed only if these cells cannot accommodate any locations of individuals (e.g., a water body within a forest area). This needs to be reflected in the observation window of the point pattern. For the torus translation method, no `NA` values are allowed at all. ```{r randomize_raster, eval = FALSE} torus_trans <- translate_raster(raster = landscape_discrete) @@ -50,9 +54,14 @@ torus_trans <- translate_raster(raster = landscape_discrete) random_walk <- randomize_raster(raster = landscape_discrete, n_random = 39) ``` +```{r randomize_raster_import, echo = FALSE, eval = TRUE} +torus_trans <- unpack_randomized(raster = shar:::data_internal$torus_trans) +random_walk <- unpack_randomized(raster = shar:::data_internal$random_walk) +``` + ## Randomize location data -To randomize the location data (i.e. the point pattern) either `fit_point_process()` or `reconstruct_pattern()` are available. The first fits either a Poisson process or a cluster process to the data. The difference to solutions from the `spatstat` package is that the number of points is always identical. The second functions reconstructs the spatial characteristics of the data using pattern reconstruction (Tscheschel & Stoyan 2006). This is advantageous for point patterns not describable by simple point process models. For both function, the number of patterns can be specified by the `n_random` argument. +To randomize the location data (i.e. the point pattern) either `fit_point_process()` or `reconstruct_pattern()` are available. The first fits a Poisson process or a cluster process to the data. The second functions reconstructs the spatial characteristics of the data using pattern reconstruction (Kirkpatrick et al. 1983; Tscheschel & Stoyan 2006). This is advantageous for point patterns not describable by simple point process models. For both function, the number of patterns can be specified by the `n_random` argument. ```{r randomize_pp, eval = FALSE} gamma_test <- fit_point_process(pattern = species_b, process = "cluster", n_random = 39) @@ -62,15 +71,16 @@ reconstruction <- reconstruct_pattern(pattern = species_b, n_random = 39, e_threshold = 0.05, method = "cluster") ``` -## Analyse results +```{r randomize_import, echo = FALSE, eval = TRUE} +gamma_test <- shar:::data_internal$gamma_test +reconstruction <- shar:::data_internal$reconstruction +``` -The most important function to analyse results is `results_habitat_association()`. This function compares the observed data to the null model data and by that is able to show significant species-habitat associations. The functions work for both, randomized environmental data or randomized location data. +## Analyse results -Please be aware that due to the randomization of the null model data, results might slightly differ between different randomization approaches (e.g., `fit_point_process()` vs. `translate_raster()`) and even for repetitions of the same approach. However, the counts of the observed data should be identical, and general results and trends should be similar. +The most important function to analyse results is `results_habitat_association()`. This function compares the observed data to the null model data and by that is able to show significant species-habitat associations. The function work for both, randomized environmental data or randomized location data. -```{r result_import, echo = FALSE, eval = TRUE} -random_walk <- unpack_randomized(raster = random_walk) -``` +Please be aware that due to the randomization of the null model data, results might slightly differ between different randomization approaches (e.g., `fit_point_process()` vs. `translate_raster()`), and even for repetitions of the same approach. Thus, the exact `lo` and `hi` thresholds might be slightly different when re-running the examples. However, the counts of the observed data should be identical, and general results and trends should be similar. ```{r results} results_habitat_association(pattern = species_a, raster = random_walk) @@ -78,22 +88,34 @@ results_habitat_association(pattern = species_a, raster = random_walk) results_habitat_association(pattern = reconstruction, raster = landscape_discrete) ``` -There is also the possibility to visualize the randomized data using the `plot()` function. +The data was created that `species_a` has a negative association to habitat 4 and `species_b` has a positive association to habitat 5, which is reflected in the results. + +Given the characteristics of the method, a positive association to one habitat inevitably leads to a negative association to at least one of the other habitats (and vice versa; Yamada et al. 2006). For example, a high amount of individual points in the positively associated habitat simultaneously mean that less individual points can be present in the other habitats. -```{r plotting, fig.align = "center", out.height = "100%", out.width = "100%", message = FALSE} -plot(random_walk) +## Utility functions +There is also the possibility to visualize the randomized data using the `plot()` function. + +```{r plot_recon, eval = FALSE} plot(reconstruction, ask = FALSE) ``` -For the randomized point pattern data, it is also possible to show the "difference" in terms of the energy (Tscheschel & Stoyan 2006) between the patterns. - -```{r energy, message = FALSE} -calculate_energy(pattern = gamma_test, return_mean = TRUE) +```{r plot_tors, fig.align = "center", out.height = "65%", out.width = "65%", message = FALSE} +col_palette <- c("#440154FF", "#3B528BFF", "#21908CFF", "#5DC863FF", "#FDE725FF") -calculate_energy(pattern = reconstruction, return_mean = TRUE) +plot(torus_trans, col = col_palette) ``` +There are many more functions, which can be found [here](https://r-spatialecology.github.io/shar/reference/index.html). + +### Citation + +The **shar** package is part of our academic work. To cite the package or acknowledge its use in publications, please cite the following paper. + +> Hesselbarth, M.H.K., (2021). shar: A R package to analyze species-habitat associations using point pattern analysis. Journal of Open Source Software, 6(67), 3811. https://doi.org/10.21105/joss.03811 + +The get a BibTex entry, please use `citation("shar")`. + ### References Fisher, W.D., 1958. On grouping for maximum homogeneity. Journal of the American Statistical Association 53, 789–798. @@ -102,8 +124,12 @@ Harms, K.E., Condit, R., Hubbell, S.P., Foster, R.B., 2001. Habitat associations Jenks, G.F., Caspall, F.C., 1971. Error in choroplethic maps: Definition, measurement, reduction. Annals of the Association of American Geographers 61, 217–244. +Kirkpatrick, S., Gelatt, C.D.Jr., Vecchi, M.P., 1983. Optimization by simulated annealing. Science 220, 671–680. + Plotkin, J.B., Potts, M.D., Leslie, N., Manokaran, N., LaFrankie, J.V., Ashton, P.S., 2000. Species-area curves, spatial aggregation, and habitat specialization in tropical forests. Journal of Theoretical Biology 207, 81–99. Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random point patterns. Computational Statistics and Data Analysis 51, 859–871. Wiegand, T., Moloney, K.A., 2014. Handbook of spatial point-pattern analysis in ecology. Chapman and Hall/CRC Press, Boca Raton. ISBN 978-1-4200-8254-8 + +Yamada, T., Tomita, A., Itoh, A., Yamakura, T., Ohkubo, T., Kanzaki, M., Tan, S., Ashton, P.S., 2006. Habitat associations of Sterculiaceae trees in a Bornean rain forest plot. Journal of Vegetation Science 17, 559–566.