From a136f582cd3816c2eee6ec7ee0e166c42f4899e3 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Tue, 31 Oct 2023 11:30:34 +0100 Subject: [PATCH 01/39] Make some tests less verbose --- tests/testthat/test-list_to_randomized.R | 9 ++++++--- tests/testthat/test-pack-unpack.R | 6 +++--- tests/testthat/test-plot_rd_ras.R | 2 +- tests/testthat/test-reconstruct_pattern.R | 1 - 4 files changed, 10 insertions(+), 8 deletions(-) diff --git a/tests/testthat/test-list_to_randomized.R b/tests/testthat/test-list_to_randomized.R index f14ac20a..1c97a819 100644 --- a/tests/testthat/test-list_to_randomized.R +++ b/tests/testthat/test-list_to_randomized.R @@ -7,7 +7,7 @@ pattern_random <- lapply(X = 1:3, function(i) { }) pattern_conv <- list_to_randomized(list = pattern_random, - observed = species_b) + observed = species_b) landscape_classified <- classify_habitats(raster = terra::rast(landscape), n = 3, style = "fisher") @@ -67,12 +67,15 @@ testthat::test_that("list_to_randomized returns errors", { testthat::test_that("list_to_randomized works with results_habitat_associations", { - res_a <- results_habitat_association(pattern = pattern_conv, raster = landscape_classified) + res_a <- results_habitat_association(pattern = pattern_conv, raster = landscape_classified, + verbose = FALSE) - res_b <- results_habitat_association(pattern = species_b, raster = raster_conv) + res_b <- results_habitat_association(pattern = species_b, raster = raster_conv, + verbose = FALSE) testthat::expect_s3_class(object = res_a, class = "data.frame") testthat::expect_s3_class(object = res_b, class = "data.frame") }) + diff --git a/tests/testthat/test-pack-unpack.R b/tests/testthat/test-pack-unpack.R index 0379b2fa..e0ba7667 100644 --- a/tests/testthat/test-pack-unpack.R +++ b/tests/testthat/test-pack-unpack.R @@ -3,13 +3,13 @@ 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) +landscape_random <- randomize_raster(landscape_classified, n_random = 2, verbose = FALSE) +landscape_ni <- randomize_raster(landscape_classified, n_random = 2, + return_input = FALSE, verbose = FALSE) x <- pack_randomized(raster = landscape_random) x_ni <- pack_randomized(raster = landscape_random) - ################################################################################ testthat::test_that("pack_randomized wraps raster", { diff --git a/tests/testthat/test-plot_rd_ras.R b/tests/testthat/test-plot_rd_ras.R index 5145eb3a..eb65151b 100644 --- a/tests/testthat/test-plot_rd_ras.R +++ b/tests/testthat/test-plot_rd_ras.R @@ -48,6 +48,6 @@ testthat::test_that("plot returns error if wrong id are selected ", { testthat::test_that("plot returns warning if more than 10 classes are present", { - testthat::expect_warning(plot(raster_random_cont), + testthat::expect_warning(plot(raster_random_cont, n = 1), regexp = "The raster has more than 10 classes. Please make sure discrete classes are provided.") }) diff --git a/tests/testthat/test-reconstruct_pattern.R b/tests/testthat/test-reconstruct_pattern.R index 6dd1bfea..29bf38e3 100644 --- a/tests/testthat/test-reconstruct_pattern.R +++ b/tests/testthat/test-reconstruct_pattern.R @@ -105,4 +105,3 @@ testthat::test_that("reconstruct_pattern returns warnings", { verbose = FALSE), regexp = "'simplify = TRUE' not possible for 'return_input = TRUE'.") }) - From 25a8d684bf73c6e01d840c7506cdf844c1d1cc37 Mon Sep 17 00:00:00 2001 From: ChrisWudel Date: Thu, 16 Nov 2023 14:04:41 +0100 Subject: [PATCH 02/39] Add function for PPR as presented in Wudel et al. 2023 --- R/reconstruct_pattern_multi_trait_marks.R | 687 ++++++++++++++++++++++ 1 file changed, 687 insertions(+) create mode 100644 R/reconstruct_pattern_multi_trait_marks.R diff --git a/R/reconstruct_pattern_multi_trait_marks.R b/R/reconstruct_pattern_multi_trait_marks.R new file mode 100644 index 00000000..44555340 --- /dev/null +++ b/R/reconstruct_pattern_multi_trait_marks.R @@ -0,0 +1,687 @@ +#' reconstruct_pattern_multi_trait_marks +#' +#' @description Pattern reconstruction of a pattern marked by multiple traits. +#' +#' @param marked_pattern ppp object with marked pattern. See Details section +#' for more information. +#' @param n_repetitions Integer representing the number of simulations to be +#' performed. +#' @param max_steps Maximum number simulation steps. +#' @param no_changes Integer representing the number of iterations +#' (per 1000 simulation steps) +#' after which the reconstruction is terminated if the energy does not decrease. +#' @param rcount Integer representing the number of intervals for which the +#' summary statistics are evaluated. +#' @param rmax Maximum distance [m] at which the summary statistics are +#' evaluated. +#' @param issue Integer that determines after how many simulations steps an +#' output occurs. +#' @param divisor Divisor in the smoothing kernel, d or r. +#' @param kernel_arg The kernel used to calculate the energy, possible kernels +#' can be: Gaussian, Epanechnikov, Rectangular, Cumulative. +#' @param timing Logical value: The computation timeis measured if this is TRUE +#' @param energy_evaluation Logical value: If this is TRUE, the procedure stores +#' the energy shares of the total energy per simulation step. +#' @param show_graphic Logical value: If this is TRUE, the procedure records the +#' point pattern during optimisation and updated +#' @param Lp Distance measure for the calculation of the energy function +#' (Lp distance, 1 ≤ p < ∞). +#' @param bw Bandwidth [m] with which the kernels are scaled, so that this is +#' the standard deviation of the smoothing kernel. +#' @param sd This is the standard deviation [m] used in the move_coordinate +#' action. +#' @param steps_tol After the value steps_tol it is checked whether the energy +#' change is smaller than tol. +#' @param tol Tolerance: stops the procedure of energy more than 1 − tol +#' times no_changes. +#' @param w_markcorr Vector of possible weightings of individual mcf's. +#' \code{c(1, 1, 1, 1, 1, 1, 1, 1)} +#' @param prob_of_actions Vector of probabilities for the actions performed. +#' \code{c(move_coordinate = 0.4, switch_coords = 0.1, exchange_mark_one = 0.1, +#' exchange_mark_two = 0.1, pick_mark_one = 0.2, pick_mark_two = 0.1)} +#' @param k Vector of values k; used only if Dk is included in w_statistics +#' elow. +#' @param w_statistics vector of named weights for optional spatial statistics +#' from the ‘spatstat’package to be included in the energy calculation.This may +#' include Dk, K, Hs, pcf. +#' +#' @details +#' A novel approach carries out a pattern reconstruction of marked dot patterns +#' as described by Tscheschel and Stoyan (2006) and Wiegand and Moloney (2014). +#' +#' One particular feature is the simultaneous consideration of both marks, +#' accounting for their correlation during reconstruction. +#' +#' The marked point pattern (PPP object) must is currently structured as follows: +#' X-coordinate, Y-coordinate, metric mark (e.g. diameter at breast height), +#' and nominal mark (e.g. tree species).It is calculated in the unit metre [m]. +#' +#' A combination of the mark correlation function and pair correlation function +#' is used for pattern description. Additional summary statistics may be +#' considered.Two randomly selected marks are chosen in each iteration, and one +#' of various actions is performed. Changes will only be retained if the +#' difference between the observed and reconstructed pattern decreases +#' (minimizing energy). +#' +#' This method is currently only suitable for homogeneous point patterns. +#' +#' A comprehensive description of the method can be +#' found in Wudel et al. (2023). Future developments will involve considering +#' more than two markers, expanding and reducing the observation window under +#' using an Edge Correction. Additionally, data input will be automatically +#' detected. +#' +#' @seealso +#' \code{\link{fit_point_process}} \cr +#' \code{\link{reconstruct_pattern}} \cr +#' \code{\link{reconstruct_pattern_marks}} +#' +#' @return rd_mar +#' +#' @examples +#' \dontrun{ +#' +#' # Random example data set +#' xr <- 500 +#' yr <- 1000 +#' N <- 400 +#' y <- runif(N, min = 0, max = yr) +#' x <- runif(N, min = 0, max = xr) +#' species <- as.factor(c("A","B")) +#' species <- sample(species, N, replace = TRUE) +#' diameter <- runif(N, 1, 40) +#' plot(x,y, col = species, pch = 19, cex = diameter*0.02) +#' random <- data.frame(x, y, diameter/ 0.1,species) +#' colnames(random) <- c("x", "y", "dbh [mm]", "Tree species") +#' +#' # Conversion to a ppp object and conversion of the metric mark to metres. +#' w <- owin(c(0, 500), c(0, 1000)) +#' marked_pattern <- as.ppp(data.frame(random), W = w) +#' marked_pattern$marks$dbh..mm.<-marked_pattern$marks$dbh..mm.*0.001 +#' +#' # Reconstruction function +#' reconstruction <- Multi_trait_point_pattern_reconstruction( +#' marked_pattern, n_repetitions = 1, max_steps = 100000, no_changes = 5, +#' rcount = 250, rmax = 25, issue = 1000, divisor = "r", +#' kernel_arg = "epanechnikov", timing = TRUE, energy_evaluation = TRUE, +#' show_graphic = TRUE, Lp = 1, bw = 0.5, sd = "step", steps_tol = 1000, +#' tol = 1e-4, w_markcorr = c(d_d=1, all=1, d_all=1, all_all=1, d_d0=1, +#' all0=1, d_all0=1, all_all0=1), prob_of_actions = c(move_coordinate = 0.4, +#' switch_coords = 0.1, exchange_mark_one = 0.1, exchange_mark_two = 0.1, +#' pick_mark_one = 0.2, pick_mark_two = 0.1), +#' k = 1, w_statistics = c() +#' ) +#' } +#' +#' @aliases reconstruct_pattern_multi_trait_marks +#' @rdname reconstruct_pattern_multi_trait_marks +#' +#' @references +#' Tscheschel, A., Stoyan, D., 2006. Statistical reconstruction of random point +#' patterns. Computational Statistics and Data Analysis 51, 859–871. +#' +#' +#' Wiegand, T., He, F., & Hubbell, S. P. (2013). A systematic comparison of +#' summary characteristics for quantifying point patterns in ecology. Ecography, 36, 92–103. +#' https://doi.org/10.1111/j.1600- 0587.2012.07361.x +#' +#' 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 +#' +#' Wudel, C., Schlicht, R., & Berger, U. (2023). Multi-trait point pattern +#' reconstruction of plant ecosystems. Methods in Ecology and Evolution, 14, 2668–2679. +#' https://doi.org/10.1111/2041-210X.14206 +#' +#' @export +reconstruct_pattern_multi_trait_marks <- function(marked_pattern, + xr = marked_pattern$window$xrange, + yr = marked_pattern$window$yrange, + n_repetitions = 1, + max_steps = 10000, + no_changes = 5, + rcount = 250, + rmax = 25, + issue = 1000 , + divisor = "r", + kernel_arg = "epanechnikov", + timing = FALSE, + energy_evaluation = FALSE, + show_graphic = FALSE, + Lp = 1, + bw = if (divisor %in% c("r","d")) 0.5 else 5, + sd = "step", + steps_tol = 1000, + tol = 1e-4, + w_markcorr = c(d_d=1, all=1, d_all=1, all_all=1, d_d0=1, all0=1, d_all0=1, all_all0=1), + prob_of_actions = c(move_coordinate = 0.4, switch_coords = 0.1, exchange_mark_one = 0.1, exchange_mark_two = 0.1, pick_mark_one = 0.2, pick_mark_two = 0.1), + k = 1, + w_statistics = c() + ) +{ + +# If several reconstructions are to be carried out, a list is created here in which the results are then saved continuously. + if(n_repetitions > 1) { + names_reconstruction <- paste0("reconstruction_", seq_len(n_repetitions)) + reconstruction_list <- vector("list", length = n_repetitions) + names(reconstruction_list) <- names_reconstruction + } +# Loop to perform multiple reconstructions. + for (t in seq_len(n_repetitions)) { + +# Definition of the product-moment function for calculating the contribution of a point at the coordinates x, y with marking. + calc_moments <- function(fn, p, exclude=NULL, x, y, mark, kernel, rmax_bw, r) { + d2 <- (p$x-x)^2 + (p$y-y)^2 + use <- d2 <= rmax_bw^2 + use[exclude] <- FALSE + z <- crossprod(p$mark[use, , drop = FALSE], + outer(sqrt(d2[use]), r, function(d, r) kernel(r, d))) + z[fn$i, , drop = FALSE] * mark[fn$j] + z[fn$j, , drop = FALSE] * mark[fn$i] + } + +# Definition of the calc_moments_full function to calculate the calc_moments function for the whole pattern. + calc_moments_full <- function(fn, p, kernel, rmax_bw, r) { + f <- 0 + for (i in seq_len(nrow(p))) { + f <- f + calc_moments(fn, p, i:nrow(p), p$x[i], p$y[i], p$mark[i, ], + kernel, rmax_bw, r) + } + rownames(f) <- paste(names(fn$i), names(fn$j), sep = ":") + f + } + +# Function for the transformation of variables to dummy variables and back + to.dummy <- function(f) { + x <- matrix(0, length(f), nlevels(f), dimnames=list(names(f), levels(f))) + x[cbind(seq_along(f), as.integer(f))] <- 1 + x + } + from.dummy <- function(x, levels=colnames(x)) { + f <- as.integer(x %*% seq_along(levels)) + levels(f) <- levels + class(f) <- "factor" + f + } + +# Compute optional spatial statistics using the spatstat package. + compute_statistics <- function(x, y, k, xr, yr) { + stat <- names(w_statistics) + names(stat) <- stat + lapply(stat, function(name) switch(name, + +# Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. + Dk = { + nnd_ <- as.matrix(nndist(x, y, k=k)) + apply(nnd_, 2, function(z) cumsum(hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) + }, +# Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. + K = { + kest<-Kest(ppp(x,y,window=owin(xr,yr)), rmax=rmax, correction="none")# , breaks = c(-Inf, r) + kest$un + }, +# Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. + pcf = { + pcfest<-pcf(ppp(x,y,window=owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") + pcfest$un + }, +# Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. + Hs = { + hest<-Hest(ppp(x,y,window=owin(xr,yr)), correction="none") #, breaks = c(-Inf, r) + hest$raw + }, + stop("unknown statistic") + )) + } + +# Defining the Energy_fun function to calculate the "energy" of the pattern (where a lower energy indicates a better match). + Energy_fun <- function(f, f0, statistics, f_, f0_, statistics_) { + result <- c( + f = sum(fn$w * rowMeans(abs( + f / nrow(p) - + f_ / nrow(p_) + )^Lp)), + f0 = sum(fn$w0 * abs( + f0 / nrow(p) - + f0_ / nrow(p_) + )^Lp), + if (length(w_statistics)) + sapply(seq_along(w_statistics), function(i) w_statistics[i] * + mean(abs(statistics[[i]] - statistics_[[i]])^Lp, na.rm = TRUE), + USE.NAMES=FALSE + ) + ) + c(energy = sum(result), result) + } + +# Load the reference point pattern as a data frame with the components x, y, mark, where x, y are the coordinates of the point and mark is a matrix representing the marks or their dummy values. + p_ <- data.frame( + x = marked_pattern$x, + y = marked_pattern$y) + + p_$mark <- cbind(`1`= 1, + diameter = marked_pattern$marks[[1]], + to.dummy(marked_pattern$marks[[2]])) + + marknames <- colnames(p_$mark) + diameter <- 2 + species <- seq_along(marknames)[-(1:2)] + +# Check whether certain requirements are met; if not, the reconstruction is aborted and an error message is displayed. + if (is.null(marked_pattern[["marks"]])) { + stop("'marked_pattern' must be marked", call. = FALSE) + } + + if (class(marked_pattern[["marks"]][[1]]) != "numeric") { + stop("mark one must be 'numeric', an example would be the DBH + (Diameter at breast height).", call. = FALSE) + } + + if (class(marked_pattern[["marks"]][[2]]) != "factor") { + stop("mark two must be a 'factor', an example would be the tree species.", + call. = FALSE) + } + + if (n_repetitions < 1) { + stop("n_repetitions must be at least 1 for the function to be executed.", + call. = FALSE) + } + +# Definition of parameters for the estimation of correlation functions. + rmin <- rmax / rcount + r <- seq(rmin, rmax, , rcount) + +# Calculation of the kernels. + kernel <- switch(kernel_arg, + epanechnikov = { + a <- bw * sqrt(5) + rmax_bw <- rmax + a + switch(divisor, + { + rmax_bw <- sqrt(rmax^2 + a/pi) + function(r, d) pmax.int(0, 1 - ((r^2-d^2)*pi/a)^2) * 0.75/a + }, + none = function(r, d) pmax.int(0, 1 - ((r-d)/a)^2) * 0.75/a, + r = function(r, d) pmax.int(0, 1 - ((r-d)/a)^2) * 0.75/(a*2*pi*r), + d = function(r, d) pmax.int(0, 1 - ((r-d)/a)^2) * 0.75/(a*2*pi*d) + ) + }, + rectangular =, box = { + a <- bw * sqrt(3) + rmax_bw <- rmax + a + switch(divisor, + { + rmax_bw <- sqrt(rmax^2 + a/pi) + function(r, d) dunif((r^2-d^2)*pi,-a,+a) + }, + none = function(r, d) dunif(r,d-a,d+a), + r = function(r, d) dunif(r,d-a,d+a)/(2*pi*r), + d = function(r, d) dunif(r,d-a,d+a)/(2*pi*d) + ) + }, + gaussian = { + rmax_bw <- Inf + switch(divisor, + function(r, d) dnorm((r^2-d^2)*pi,0,sd=bw), + none = function(r, d) dnorm(r,d,sd = bw), + r = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*r), + d = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*d) + ) + }, + + cumulative = { + rmax_bw <- rmax + switch(divisor, + function(r, d) as.numeric(d <= r), + none = function(r, d) as.numeric(d <= r), + r = function(r, d) (d <= r) / (2*pi*r), + d = function(r, d) (d <= r) / (2*pi*d) + ) + } + ) + +# If separate summary statistics (Dk, K, Hs, pcf) are to be included in the energy calculation, the "spatstat" package is loaded here and an error message is displayed if it is not installed. + if (length(w_statistics)) { + if("spatstat" %in% rownames(installed.packages()) == FALSE) { + stop("the package 'spatstat' must be installed on your computer for the + application if the w_statistics option is used.", call. = FALSE) + } + require (spatstat) + } + +# Determination of the weightings of the mark correlation functions. + fn <- list() + for (i in seq_along(marknames)) for (j in seq_along(marknames)) if (i <= j) { + fn$i <- c(fn$i,i) + fn$j <- c(fn$j,j) + fn$w <- c(fn$w, + if (i == diameter && j == diameter) w_markcorr["d_d"] + else if(i == 1 || j ==1) w_markcorr["all"] + else if(i == diameter || j == diameter) w_markcorr["d_all"] + else w_markcorr["all_all"]) + fn$w0<-c(fn$w0, + if (i == diameter && j == diameter) w_markcorr["d_d0"] + else if(i == 1 || j == 1) w_markcorr["all0"] + else if(i == diameter || j == diameter) w_markcorr["d_all0"] + else w_markcorr["all_all0"]) + } + names(fn$i) <- marknames[fn$i] + names(fn$j) <- marknames[fn$j] + +# Defines the initial state of the new ponit pattern. + n <- nrow(p_) + xwr <- obs_window$xrange + ywr <- obs_window$yrange + p <- p_[sample.int(nrow(p_),ceiling(nrow(p_)*((diff(xwr)*diff(ywr))/(diff(xr)*diff(yr)))),TRUE), ]# rpois + p$x <- runif(nrow(p), xwr[1], xwr[2]) + p$y <- runif(nrow(p), ywr[1], ywr[2]) + p$mark[, diameter] <- quantile(p_$mark[, diameter], + probs = runif(nrow(p), 0, 1), type = 4) + p$mark[, species] <- p_$mark[, species, drop = FALSE][ + sample.int(nrow(p_), nrow(p), replace = TRUE),, drop = FALSE] + +# Calculates the functions for the reference and the new dot pattern as well as calculating the "energy" that measures their distance. + f_ <- calc_moments_full(fn, p_, kernel, rmax_bw, r) + f0_ <- colSums(p_$mark[, fn$i] * p_$mark[, fn$j]) + names(f0_) <- rownames(f_) + statistics_<- compute_statistics(p_$x, p_$y, k, xr, yr) + f <- calc_moments_full(fn, p, kernel, rmax_bw, r) + f0 <- colSums(p$mark[, fn$i] * p$mark[, fn$j]) + names(f0) <- rownames(f) + statistics <- compute_statistics(p$x, p$y, k, xwr, ywr) + +# Prepare the graphical output. + if(show_graphic == TRUE) { + par(mfrow = 1:2) + num_species <- from.dummy (p_$mark[, species, drop = FALSE]) + plot(y~x, p_, pch=19, col= 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter], xlim = xr, + ylim = yr, xaxs ="i", yaxs ="i", main ="Reference", xlab ="x [m]", + ylab ="y [m]") + text(p_$x, p_$y, as.integer(num_species), cex=0.7) + + plot(y~x, p, type = "n", + xlim = xwr, ylim = ywr, xaxs = "i", yaxs = "i", main = "Reconstructed", + xlab = "x [m]", ylab = "y [m]") + clip(xwr[1], xwr[2], ywr[1], ywr[2]) + } + +# Show warning if certain distances between pairs of trees are not present. + energy <- Energy_fun(f, f0, statistics, f_, f0_, statistics_)["energy"] + +# Prepares variables for the storage of progress. + energy_launch <- as.vector(energy) + energy_course <- data.frame(i = seq(from = 1, to = max_steps,by = 1), + energy = NA) + energy_improvement <- 0L + number_of_actions <- integer(0) + number_of_actions_with_energy_improvement <- integer(0) + no_changes_energy <- Inf + no_changes_counter <- 0L + step_list <- integer(0) + action_list <- character(0) + Energy_list <- numeric(0) + +# loop to improve the newly created dot pattern by reducing its energy. + step <- 0L + system.time(repeat { + energy_course[step, 2] <- energy + +# Updating the graphical output of all "issue" steps. + if (step %% issue == 0) { + if(show_graphic == TRUE) { + rect(xwr[1], ywr[1], xwr[2], ywr[2], col="white") + num_species <- from.dummy (p$mark[, species, drop = FALSE]) + + points(y~x, p, pch = 19, col = 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter]) + text(p$x, p$y, as.integer(num_species), cex = 0.7) + } + +# Generates text output with the current calculated values (for example the energy), this is updated every "issue" simulation step. + message("> Progress:", if(n_repetitions > 1) names_reconstruction[[t]], " || iterations: ", step, + " || Simulation progress: ", floor(step/max_steps * 100), "%", + " || energy = ", round(energy, 5), " || energy improvement = ", + energy_improvement, "\t\t\r", appendLF = FALSE) + + Sys.sleep(0) + flush.console() + } +# the next code block aborts the reconstruction if the energy decreases by less than tol in "no_changes" intervals of steps_tol simulation steps. + if (step %% steps_tol == 0) { + if(energy > no_changes_energy * (1-tol)) { + no_changes_counter <- no_changes_counter + 1L + if(no_changes_counter == no_changes) { + message("the simulation was terminated, because the energy did not + decrease in ", no_changes * issue, " simulation steps.") + stop_criterion <- "no_changes" + break + } + } else { + no_changes_counter<-0L + stop_criterion<-"max_steps" + } + no_changes_energy <- energy + } + + if (step < max_steps) step <- step + 1 else break + action <- sample(c("move_coordinate", "switch_coords", "pick_mark_one", + "pick_mark_two", "exchange_mark_one", "exchange_mark_two"), + 1,, prob_of_actions) + number_of_actions[action] <- (if (is.na(number_of_actions[action])) + 0L else number_of_actions[action]) + 1L + + statistics.new <- statistics + +# Switch selection for the possible changes to the reconstructed point pattern for energy minimisation (probabilities of how often an action is taken: 40%, 10%, 20%, 10%, 10%). + switch(action, + # Displacement of coordinates of a point in the new point pattern, is applied in xx% of the cases. + move_coordinate = { + i <- sample.int(nrow(p), 1, replace = TRUE) + s <- if(sd=="step") nrow(p) * 1 / step else sd + x <- rnorm(1, p$x[i], diff(xwr) * s) %% xwr[2] + y <- rnorm(1, p$y[i], diff(ywr) * s) %% ywr[2] + mdiff <- p$mark[i, ] + f.new <- f - calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, + kernel, rmax_bw, r) + + calc_moments(fn, p, i, x, y, mdiff, kernel, rmax_bw, r) + f0.new<- f0 + statistics.new <- compute_statistics(replace(p$x, i, x), replace(p$y, i, y), k, xwr, ywr) + }, + + # Swaps the coordinates of two randomly drawn points from the new point pattern, applied in xx% of the trap. + switch_coords = { + i <- sample.int(nrow(p), 2, replace = FALSE) + mdiff <- p$mark[i[1], ] - p$mark[i[2], ] + f.new <- f - calc_moments(fn, p, i, p$x[i[1]], p$y[i[1]], + mdiff, + kernel, rmax_bw, r) + + calc_moments(fn, p, i, p$x[i[2]], + p$y[i[2]], mdiff, kernel, rmax_bw, r) + f0.new<- f0 + }, + + # Displacement of coordinates of a point in the new point pattern, applied in xx% of the cases. + exchange_mark_one = { + i <- sample.int(nrow(p), 2, replace = FALSE) + m <- p$mark[i, ] + m[, diameter] <- m[2:1, diameter] + mdiff <- m[1, ] - p$mark[i[1], ] + q <- p[i, ] + q$mark[1, ] <- m[1, ] + f.new <- f + calc_moments(fn, p, i[1], p$x[i[1]], + p$y[i[1]], mdiff, kernel, + rmax_bw, r) - + calc_moments(fn, p, i, p$x[i[2]], p$y[i[2]], + mdiff, kernel, rmax_bw, r) - + calc_moments(fn, q, 2, q$x[2], q$y[2], mdiff, + kernel, rmax_bw, r) + f0.new <- f0 + m[1,fn$i] * m[1, fn$j] - + p$mark[i[1], fn$i] * p$mark[i[1], fn$j] + m[2, fn$i] * + m[2, fn$j] - p$mark[i[2], fn$i] * + p$mark[i[2], fn$j] + }, + + # Swaps the type assignment of two randomly drawn points from the new point pattern, applied in 10% of the trap. + exchange_mark_two = { + i <- sample.int(nrow(p), 2, replace = FALSE) + m <- p$mark[i, ] + m[, species] <- m[2:1, species] + mdiff <- m[1, ] - p$mark[i[1], ] + q <- p[i, ] + q$mark <- m + + f.new <- f + calc_moments(fn, p, i[1], p$x[i[1]], p$y[i[1]], + mdiff, kernel, rmax_bw, r) - + calc_moments(fn, p, i, p$x[i[2]], p$y[i[2]], + mdiff, kernel, rmax_bw, r) - + calc_moments(fn, q, 2, q$x[2], q$y[2], mdiff, + kernel, rmax_bw, r) + f0.new <- f0 + m[1, fn$i] * m[1, fn$j] - p$mark[i[1], fn$i] * + p$mark[i[1], fn$j] + m[2, fn$i] * m[2, fn$j] - + p$mark[i[2], fn$i] * p$mark[i[2], fn$j] + }, + + # If the distribution (continuous function) of the diameter of the reference pattern generates a randomly drawn value for a randomly selected point in the new point pattern, the trap is applied in 20%. + pick_mark_one = { + i <- sample.int(nrow(p), 1, replace = TRUE) + m <- p$mark[i, ] + m[diameter] <-quantile(p_$mark[,diameter],probs = runif(1,0,1), + type = 4) + mdiff <- m - p$mark[i, ] + f.new <- f + calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, + kernel, rmax_bw, r) + f0.new<- f0 + m[fn$i] * m[fn$j] - p$mark[i, fn$i] * + p$mark[i, fn$j] + }, + + # Draws a random value for a point from the new point pattern from the type distribution (discrete function) of the reference pattern, is applied in 10% of the traps. + pick_mark_two = { + i <- sample.int(nrow(p), 1, replace = TRUE) + j <- sample.int(nrow(p_), 1, replace = TRUE) + m <- p$mark[i, ] + m[species] <- p_$mark[j, species] + mdiff <- m - p$mark[i, ] + f.new <- f + calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, + kernel, rmax_bw, r) + f0.new <- f0 + m[fn$i] * m[fn$j] - p$mark[i, fn$i] * + p$mark[i, fn$j] + }, + stop("undefined case") + ) +# calculates the energy + Energy <- Energy_fun(f.new, f0.new, statistics.new, f_, f0_, statistics_) + energy.new<-Energy[["energy"]] + +# Sets the currently calculated energy as the new reference value if it is less than the previous energy. + if(energy.new >= energy) next + f <- f.new + f0 <- f0.new + statistics <- statistics.new + energy <- energy.new + +# Executes the previously defined actions. + switch(action, + move_coordinate = { + p$x[i] <- x + p$y[i] <- y + }, + switch_coords = { + p$x[i] <- p$x[rev(i)] + p$y[i] <- p$y[rev(i)] + }, + pick_mark_one =, + pick_mark =, + pick_mark_two =, + exchange_mark_one =, + exchange_mark_two = + { + p$mark[i, ] <- m + }, + stop("undefined case") + ) + +# Saves the intermediate results and increases running numbers. + if (energy_evaluation == TRUE) { + step_list <-c(step_list,step) + action_list <-c(action_list, action) + Energy_list <-rbind(Energy_list, Energy) + number_of_actions_with_energy_improvement[action] <- + (if (is.na(number_of_actions_with_energy_improvement[action])) + 0L else number_of_actions_with_energy_improvement[action]) + 1L + } + energy_improvement <- energy_improvement + 1L + +# End of reconstruction loop. + +}) -> process.time + message("\n") + + # Saves all results Transfers them to the "reconstruction" list. + p_$species <- from.dummy(p_$mark[, species, drop = FALSE]) + p$species <- from.dummy(p$mark[, species, drop = FALSE]) + method <- "Reconstruction of a homogeneous point pattern" + Parameter_setting <- list(n_repetitions=n_repetitions, max_steps=max_steps, + no_changes=no_changes, rcount=rcount, rmax=rmax, + issue=issue, divisor=divisor, kernel_arg=kernel_arg, + timing=timing, energy_evaluation=energy_evaluation, + show_graphic=show_graphic, Lp=Lp, k=k, bw=bw, sd=sd, + prob_of_actions=prob_of_actions, + w_markcorr=w_markcorr, w_statistics=w_statistics) + iterations <- step + energy_current <- energy_course[step, 2] + adapted_p_<- data.frame(p_$x, p_$y, p_$mark[,2], p_$species) + colnames(adapted_p_)<-c("x", "y", "diameter", "species") + adapted_p<- data.frame(p$x, p$y, p$mark[,2], p$species) + colnames(adapted_p)<-c("x", "y", "diameter", "species") + win_change<- if (xr[1] != xwr[1] | xr[2] != xwr[2] | yr[1] != ywr[1] | yr[2] != ywr[2]) {TRUE}else{FALSE} + +# Saves the results. + reconstruction <- + list( reference = adapted_p_, + reconstructed = adapted_p, + window = c(xr, yr), + obs_window = + if (win_change == TRUE) { + c(xwr, ywr) + }, + r = r, + f_reference = f_, + f_reconstructed = f, + Parameter_setting = Parameter_setting, + method = method, + stop_criterion = stop_criterion, + iterations = step, + simulation_time = + if (timing == TRUE) { + paste(round(process.time[3], 2), "s") + }, + energy_launch = energy_launch, + energy_course = energy_course, + energy_current = energy_current, + energy_improvement = energy_improvement, + number_of_actions = + if(energy_evaluation == TRUE) { + number_of_actions + }, + number_of_actions_with_energy_improvement = + if(energy_evaluation == TRUE) { + number_of_actions_with_energy_improvement + }, + energy_details = + if(energy_evaluation == TRUE) { + data.frame(step_list, action_list, Energy_list) + }) + + if (!(timing && energy_evaluation && win_change)) { + reconstruction <- reconstruction[-which(sapply(reconstruction, is.null))] + } + +# Adds the results of further reconstructions to the "reconstruction" list if several are performed. + if (n_repetitions > 1) { + reconstruction_list[[t]] <- reconstruction + } + } + + if(n_repetitions > 1) { + reconstruction_list + } else { + reconstruction + } + } From 46dd486f25a0762617a440f360b8792c77fa8aa6 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 17 Nov 2023 15:49:03 +0100 Subject: [PATCH 03/39] Update DESCRIPTION, CITATION etc. --- DESCRIPTION | 4 +- NEWS.md | 4 + SHAR.Rproj | 44 +++++------ cran-comments.md | 185 ++++++++++++++++++++++++----------------------- inst/CITATION | 38 ++++++---- 5 files changed, 148 insertions(+), 127 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index f3aed639..3a714fc5 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,11 +1,13 @@ Type: Package Package: shar Title: Species-Habitat Associations -Version: 2.1.1 +Version: 2.2 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", role = "aut", comment = c(ORCID = "0000-0002-3042-5435")), + person("Chris", "Wudel", email = "chris.wudel@tu-dresden.de", + role = "aut", comment = c(ORCID = "0000-0003-0446-4665")), person("Zeke", "Marshall", email = "ee18zm@leeds.ac.uk", role = "ctb", comment = c(ORCID = "0000-0001-9260-7827")), person("Thomas", "Etherington", email = "teth001@aucklanduni.ac.nz", diff --git a/NEWS.md b/NEWS.md index aa20ed0f..2bd619cc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# shar 2.2 +* Improvements + * Added a new function `reconstruct_pattern_multi` + # shar 2.1.1 * Bugfixes (thanks to @baddstats) diff --git a/SHAR.Rproj b/SHAR.Rproj index 8f290f92..3509a420 100644 --- a/SHAR.Rproj +++ b/SHAR.Rproj @@ -1,22 +1,22 @@ -Version: 1.0 - -RestoreWorkspace: Default -SaveWorkspace: Default -AlwaysSaveHistory: Default - -EnableCodeIndexing: Yes -UseSpacesForTab: Yes -NumSpacesForTab: 2 -Encoding: UTF-8 - -RnwWeave: Sweave -LaTeX: pdfLaTeX - -AutoAppendNewline: Yes -StripTrailingWhitespace: Yes - -BuildType: Package -PackageUseDevtools: Yes -PackageInstallArgs: --no-multiarch --with-keep.source -PackageCheckArgs: --no-multiarch --as-cran -PackageRoxygenize: rd,collate,namespace,vignette +Version: 1.0 + +RestoreWorkspace: Default +SaveWorkspace: Default +AlwaysSaveHistory: Default + +EnableCodeIndexing: Yes +UseSpacesForTab: Yes +NumSpacesForTab: 2 +Encoding: UTF-8 + +RnwWeave: Sweave +LaTeX: pdfLaTeX + +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes + +BuildType: Package +PackageUseDevtools: Yes +PackageInstallArgs: --no-multiarch --with-keep.source +PackageCheckArgs: --no-multiarch --as-cran +PackageRoxygenize: rd,collate,namespace,vignette diff --git a/cran-comments.md b/cran-comments.md index 26ae2cd2..2c3e03e0 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -1,91 +1,94 @@ -For details changes, please see NEWS.md - -## shar 2.1.1 -Minor bug fixes - -## shar 2.1 -Speed improvements - -## shar 2.0.4 -Small speed improvements - -## shar 2.0.3 -This is a re-submission. The below error has now been fixed - -``` -Package CITATION file contains call(s) to old-style citEntry(). Please use bibentry() instead. -``` - -Improvements of existing functions - -## shar 2.0.2 -Update `spatstat` dependency - -## shar 2.0.1 -Fixing external repo, typo maintainer name - -## shar 2.0.0 -* Replace `raster` with `terra` and several code changes related to this -* Fixed issue of earlier submission "Unknown, possibly misspelled, fields in DESCRIPTION: 'Remotes'" - -## shar 1.3.2 -Improvements of existing functions and `spatstat` update - -## shar 1.3.1 -Bug fixes - -## shar 1.3 -Improvements of general package structure - -## shar 1.2.1 -Improvement of existing functions - -## shar 1.2 -Update `spatstat` dependencies - -## shar 1.1.1 -Minor improvements and new license - -## shar 1.1 -Improvements of existing functions - -## shar 1.0.1 -Some minor improvements to existing functions - -## shar 1.0 -New functionality and renaming/splitting of some functions - -## shar 0.5 -Improvements of existing functions - -## shar 0.4 -Improvements of general package structure - -## shar 0.3.1 -Some minor bugfixes - -## shar 0.3 -Some bugfixes and improvements of existing functions as well as new functions - -## shar 0.2 -Some bugfixes and improvements of existing functions - -## Review CRAN submission -1. Thanks. Please omit the redundant "in R". - -* Done as suggested - -2. Is there some reference about the method you can add in the Description field in the form Authors (year) ? - -* Added "Methods are mainly based on Plotkin et al. (2000) and Harms et al. (2001) ." to the description field - -3. Thanks, we see: - running tests for arch 'i386' ... [320s] OK - running tests for arch 'x64' ... [356s] OK -which is together alreadyv more than 10 min which is the CRAN threshold for a package check. Can you pls simplify the test cases or run less important tests only conditionally if some env var is set that you only define on your machine? - -* The test run faster now (checked on win-builder.r-project) - running tests for arch 'i386' ... [152s] OK - running tests for arch 'x64' ... [164s] OK - -* Renamed package from `SHAR` to `shar` +For details changes, please see NEWS.md + +# shar 2.2 +Added new functions + +## shar 2.1.1 +Minor bug fixes + +## shar 2.1 +Speed improvements + +## shar 2.0.4 +Small speed improvements + +## shar 2.0.3 +This is a re-submission. The below error has now been fixed + +``` +Package CITATION file contains call(s) to old-style citEntry(). Please use bibentry() instead. +``` + +Improvements of existing functions + +## shar 2.0.2 +Update `spatstat` dependency + +## shar 2.0.1 +Fixing external repo, typo maintainer name + +## shar 2.0.0 +* Replace `raster` with `terra` and several code changes related to this +* Fixed issue of earlier submission "Unknown, possibly misspelled, fields in DESCRIPTION: 'Remotes'" + +## shar 1.3.2 +Improvements of existing functions and `spatstat` update + +## shar 1.3.1 +Bug fixes + +## shar 1.3 +Improvements of general package structure + +## shar 1.2.1 +Improvement of existing functions + +## shar 1.2 +Update `spatstat` dependencies + +## shar 1.1.1 +Minor improvements and new license + +## shar 1.1 +Improvements of existing functions + +## shar 1.0.1 +Some minor improvements to existing functions + +## shar 1.0 +New functionality and renaming/splitting of some functions + +## shar 0.5 +Improvements of existing functions + +## shar 0.4 +Improvements of general package structure + +## shar 0.3.1 +Some minor bugfixes + +## shar 0.3 +Some bugfixes and improvements of existing functions as well as new functions + +## shar 0.2 +Some bugfixes and improvements of existing functions + +## Review CRAN submission +1. Thanks. Please omit the redundant "in R". + +* Done as suggested + +2. Is there some reference about the method you can add in the Description field in the form Authors (year) ? + +* Added "Methods are mainly based on Plotkin et al. (2000) and Harms et al. (2001) ." to the description field + +3. Thanks, we see: + running tests for arch 'i386' ... [320s] OK + running tests for arch 'x64' ... [356s] OK +which is together alreadyv more than 10 min which is the CRAN threshold for a package check. Can you pls simplify the test cases or run less important tests only conditionally if some env var is set that you only define on your machine? + +* The test run faster now (checked on win-builder.r-project) + running tests for arch 'i386' ... [152s] OK + running tests for arch 'x64' ... [164s] OK + +* Renamed package from `SHAR` to `shar` diff --git a/inst/CITATION b/inst/CITATION index 9272504b..a17e8c59 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -1,13 +1,25 @@ -citHeader("To cite 'shar' in publications and acknowledge its use, please cite the following paper:") - -bibentry( - bibtype = "Article", - title = "shar: A R package to analyze species-habitat associations using point pattern analysis", - author = "Maximilian H.K. Hesselbarth", - journal = "Journal of Open Source Software", - year = "2021", - volume = "6", - number = "67", - pages = "3811", - doi = "10.21105/joss.03811", - textVersion = "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") +c( + bibentry( + bibtype = "Article", + title = "shar: A R package to analyze species-habitat associations using point pattern analysis", + author = "Maximilian H.K. Hesselbarth", + journal = "Journal of Open Source Software", + year = "2021", + volume = "6", + number = "67", + pages = "3811", + doi = "10.21105/joss.03811", + header = "To cite 'shar' in publications and acknowledge its use, please cite the following paper:" + ), + bibentry( + bibtype = "Article", + title = "Multi-trait point pattern reconstruction of plant ecosystems.", + author = c(person("Chris", "Wudel"), person("Robert", "Schlicht"), person("Uta", "Berger")), + journal = "Methods in Ecology and Evolution", + year = "2023", + volume = "14", + pages = "2668–2679", + doi = "https://doi.org/10.1111/2041-210X.14206", + header = "If you use the 'reconstruct_pattern_multi' function, please also cite:" + ) +) From ef4fb9de50bcbb7ce9599d6a680ac812149bd8a0 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 17 Nov 2023 17:16:27 +0100 Subject: [PATCH 04/39] Rename and import updates to reconstruct_pattern_multi --- NAMESPACE | 1 + ...it_marks.R => reconstruct_pattern_multi.R} | 282 +++++++++--------- man/reconstruct_pattern_multi.Rd | 188 ++++++++++++ man/shar.Rd | 1 + 4 files changed, 330 insertions(+), 142 deletions(-) rename R/{reconstruct_pattern_multi_trait_marks.R => reconstruct_pattern_multi.R} (74%) create mode 100644 man/reconstruct_pattern_multi.Rd diff --git a/NAMESPACE b/NAMESPACE index 2fbc51de..bd1b9e47 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -15,6 +15,7 @@ export(plot_energy) export(randomize_raster) export(reconstruct_pattern) export(reconstruct_pattern_marks) +export(reconstruct_pattern_multi) export(results_habitat_association) export(translate_raster) export(unpack_randomized) diff --git a/R/reconstruct_pattern_multi_trait_marks.R b/R/reconstruct_pattern_multi.R similarity index 74% rename from R/reconstruct_pattern_multi_trait_marks.R rename to R/reconstruct_pattern_multi.R index 44555340..11acc36b 100644 --- a/R/reconstruct_pattern_multi_trait_marks.R +++ b/R/reconstruct_pattern_multi.R @@ -1,9 +1,10 @@ -#' reconstruct_pattern_multi_trait_marks +#' reconstruct_pattern_multi #' #' @description Pattern reconstruction of a pattern marked by multiple traits. #' #' @param marked_pattern ppp object with marked pattern. See Details section #' for more information. +#' @param xr,yr Add description #' @param n_repetitions Integer representing the number of simulations to be #' performed. #' @param max_steps Maximum number simulation steps. @@ -100,7 +101,7 @@ #' marked_pattern$marks$dbh..mm.<-marked_pattern$marks$dbh..mm.*0.001 #' #' # Reconstruction function -#' reconstruction <- Multi_trait_point_pattern_reconstruction( +#' reconstruction <- reconstruct_pattern_multi( #' marked_pattern, n_repetitions = 1, max_steps = 100000, no_changes = 5, #' rcount = 250, rmax = 25, issue = 1000, divisor = "r", #' kernel_arg = "epanechnikov", timing = TRUE, energy_evaluation = TRUE, @@ -113,18 +114,17 @@ #' ) #' } #' -#' @aliases reconstruct_pattern_multi_trait_marks -#' @rdname reconstruct_pattern_multi_trait_marks +#' @aliases reconstruct_pattern_multi +#' @rdname reconstruct_pattern_multi #' #' @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., He, F., & Hubbell, S. P. (2013). A systematic comparison of -#' summary characteristics for quantifying point patterns in ecology. Ecography, 36, 92–103. -#' https://doi.org/10.1111/j.1600- 0587.2012.07361.x -#' #' 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 #' @@ -133,31 +133,29 @@ #' https://doi.org/10.1111/2041-210X.14206 #' #' @export -reconstruct_pattern_multi_trait_marks <- function(marked_pattern, - xr = marked_pattern$window$xrange, - yr = marked_pattern$window$yrange, - n_repetitions = 1, - max_steps = 10000, - no_changes = 5, - rcount = 250, - rmax = 25, - issue = 1000 , - divisor = "r", - kernel_arg = "epanechnikov", - timing = FALSE, - energy_evaluation = FALSE, - show_graphic = FALSE, - Lp = 1, - bw = if (divisor %in% c("r","d")) 0.5 else 5, - sd = "step", - steps_tol = 1000, - tol = 1e-4, - w_markcorr = c(d_d=1, all=1, d_all=1, all_all=1, d_d0=1, all0=1, d_all0=1, all_all0=1), - prob_of_actions = c(move_coordinate = 0.4, switch_coords = 0.1, exchange_mark_one = 0.1, exchange_mark_two = 0.1, pick_mark_one = 0.2, pick_mark_two = 0.1), - k = 1, - w_statistics = c() - ) -{ +reconstruct_pattern_multi <- function(marked_pattern, + xr = marked_pattern$window$xrange, + yr = marked_pattern$window$yrange, + n_repetitions = 1, + max_steps = 10000, + no_changes = 5, + rcount = 250, + rmax = 25, + issue = 1000 , + divisor = "r", + kernel_arg = "epanechnikov", + timing = FALSE, + energy_evaluation = FALSE, + show_graphic = FALSE, + Lp = 1, + bw = if (divisor %in% c("r","d")) 0.5 else 5, + sd = "step", + steps_tol = 1000, + tol = 1e-4, + w_markcorr = c(d_d=1, all=1, d_all=1, all_all=1, d_d0=1, all0=1, d_all0=1, all_all0=1), + prob_of_actions = c(move_coordinate = 0.4, switch_coords = 0.1, exchange_mark_one = 0.1, exchange_mark_two = 0.1, pick_mark_one = 0.2, pick_mark_two = 0.1), + k = 1, + w_statistics= c()) { # If several reconstructions are to be carried out, a list is created here in which the results are then saved continuously. if(n_repetitions > 1) { @@ -168,89 +166,89 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, # Loop to perform multiple reconstructions. for (t in seq_len(n_repetitions)) { -# Definition of the product-moment function for calculating the contribution of a point at the coordinates x, y with marking. - calc_moments <- function(fn, p, exclude=NULL, x, y, mark, kernel, rmax_bw, r) { - d2 <- (p$x-x)^2 + (p$y-y)^2 - use <- d2 <= rmax_bw^2 - use[exclude] <- FALSE - z <- crossprod(p$mark[use, , drop = FALSE], - outer(sqrt(d2[use]), r, function(d, r) kernel(r, d))) - z[fn$i, , drop = FALSE] * mark[fn$j] + z[fn$j, , drop = FALSE] * mark[fn$i] - } + # Definition of the product-moment function for calculating the contribution of a point at the coordinates x, y with marking. + calc_moments <- function(fn, p, exclude=NULL, x, y, mark, kernel, rmax_bw, r) { + d2 <- (p$x-x)^2 + (p$y-y)^2 + use <- d2 <= rmax_bw^2 + use[exclude] <- FALSE + z <- crossprod(p$mark[use, , drop = FALSE], + outer(sqrt(d2[use]), r, function(d, r) kernel(r, d))) + z[fn$i, , drop = FALSE] * mark[fn$j] + z[fn$j, , drop = FALSE] * mark[fn$i] + } -# Definition of the calc_moments_full function to calculate the calc_moments function for the whole pattern. - calc_moments_full <- function(fn, p, kernel, rmax_bw, r) { - f <- 0 - for (i in seq_len(nrow(p))) { - f <- f + calc_moments(fn, p, i:nrow(p), p$x[i], p$y[i], p$mark[i, ], - kernel, rmax_bw, r) + # Definition of the calc_moments_full function to calculate the calc_moments function for the whole pattern. + calc_moments_full <- function(fn, p, kernel, rmax_bw, r) { + f <- 0 + for (i in seq_len(nrow(p))) { + f <- f + calc_moments(fn, p, i:nrow(p), p$x[i], p$y[i], p$mark[i, ], + kernel, rmax_bw, r) + } + rownames(f) <- paste(names(fn$i), names(fn$j), sep = ":") + f } - rownames(f) <- paste(names(fn$i), names(fn$j), sep = ":") - f - } -# Function for the transformation of variables to dummy variables and back - to.dummy <- function(f) { - x <- matrix(0, length(f), nlevels(f), dimnames=list(names(f), levels(f))) - x[cbind(seq_along(f), as.integer(f))] <- 1 - x + # Function for the transformation of variables to dummy variables and back + to_dummy <- function(f) { + x <- matrix(0, length(f), nlevels(f), dimnames=list(names(f), levels(f))) + x[cbind(seq_along(f), as.integer(f))] <- 1 + x } - from.dummy <- function(x, levels=colnames(x)) { - f <- as.integer(x %*% seq_along(levels)) - levels(f) <- levels - class(f) <- "factor" - f + from_dummy <- function(x, levels=colnames(x)) { + f <- as.integer(x %*% seq_along(levels)) + levels(f) <- levels + class(f) <- "factor" + f } -# Compute optional spatial statistics using the spatstat package. - compute_statistics <- function(x, y, k, xr, yr) { - stat <- names(w_statistics) - names(stat) <- stat - lapply(stat, function(name) switch(name, - -# Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. - Dk = { - nnd_ <- as.matrix(nndist(x, y, k=k)) - apply(nnd_, 2, function(z) cumsum(hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) - }, -# Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. - K = { - kest<-Kest(ppp(x,y,window=owin(xr,yr)), rmax=rmax, correction="none")# , breaks = c(-Inf, r) - kest$un - }, -# Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. - pcf = { - pcfest<-pcf(ppp(x,y,window=owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") - pcfest$un - }, -# Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. - Hs = { - hest<-Hest(ppp(x,y,window=owin(xr,yr)), correction="none") #, breaks = c(-Inf, r) - hest$raw - }, - stop("unknown statistic") - )) - } + # Compute optional spatial statistics using the spatstat package. + compute_statistics <- function(x, y, k, xr, yr) { + stat <- names(w_statistics) + names(stat) <- stat + lapply(stat, function(name) switch(name, + + # Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. + Dk = { + nnd_ <- as.matrix(nndist(x, y, k=k)) + apply(nnd_, 2, function(z) cumsum(hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) + }, + # Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. + K = { + kest<-Kest(ppp(x,y,window=owin(xr,yr)), rmax=rmax, correction="none")# , breaks = c(-Inf, r) + kest$un + }, + # Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. + pcf = { + pcfest<-pcf(ppp(x,y,window=owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") + pcfest$un + }, + # Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. + Hs = { + hest<-Hest(ppp(x,y,window=owin(xr,yr)), correction="none") #, breaks = c(-Inf, r) + hest$raw + }, + stop("unknown statistic") + )) + } -# Defining the Energy_fun function to calculate the "energy" of the pattern (where a lower energy indicates a better match). - Energy_fun <- function(f, f0, statistics, f_, f0_, statistics_) { - result <- c( - f = sum(fn$w * rowMeans(abs( - f / nrow(p) - - f_ / nrow(p_) - )^Lp)), - f0 = sum(fn$w0 * abs( - f0 / nrow(p) - - f0_ / nrow(p_) - )^Lp), - if (length(w_statistics)) - sapply(seq_along(w_statistics), function(i) w_statistics[i] * - mean(abs(statistics[[i]] - statistics_[[i]])^Lp, na.rm = TRUE), - USE.NAMES=FALSE - ) - ) - c(energy = sum(result), result) - } + # Defining the Energy_fun function to calculate the "energy" of the pattern (where a lower energy indicates a better match). + energy_fun <- function(f, f0, statistics, f_, f0_, statistics_) { + result <- c( + f = sum(fn$w * rowMeans(abs( + f / nrow(p) - + f_ / nrow(p_) + )^Lp)), + f0 = sum(fn$w0 * abs( + f0 / nrow(p) - + f0_ / nrow(p_) + )^Lp), + if (length(w_statistics)) + sapply(seq_along(w_statistics), function(i) w_statistics[i] * + mean(abs(statistics[[i]] - statistics_[[i]])^Lp, na.rm = TRUE), + USE.NAMES=FALSE + ) + ) + c(energy = sum(result), result) + } # Load the reference point pattern as a data frame with the components x, y, mark, where x, y are the coordinates of the point and mark is a matrix representing the marks or their dummy values. p_ <- data.frame( @@ -259,7 +257,7 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, p_$mark <- cbind(`1`= 1, diameter = marked_pattern$marks[[1]], - to.dummy(marked_pattern$marks[[2]])) + to_dummy(marked_pattern$marks[[2]])) marknames <- colnames(p_$mark) diameter <- 2 @@ -270,12 +268,12 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, stop("'marked_pattern' must be marked", call. = FALSE) } - if (class(marked_pattern[["marks"]][[1]]) != "numeric") { + if (!inherits(marked_pattern[["marks"]][[1]], "numeric")) { stop("mark one must be 'numeric', an example would be the DBH (Diameter at breast height).", call. = FALSE) } - if (class(marked_pattern[["marks"]][[2]]) != "factor") { + if (!inherits(marked_pattern[["marks"]][[2]], "factor")) { stop("mark two must be a 'factor', an example would be the tree species.", call. = FALSE) } @@ -310,20 +308,20 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, switch(divisor, { rmax_bw <- sqrt(rmax^2 + a/pi) - function(r, d) dunif((r^2-d^2)*pi,-a,+a) + function(r, d) stats::dunif((r^2-d^2)*pi,-a,+a) }, - none = function(r, d) dunif(r,d-a,d+a), - r = function(r, d) dunif(r,d-a,d+a)/(2*pi*r), - d = function(r, d) dunif(r,d-a,d+a)/(2*pi*d) + none = function(r, d) stats::dunif(r,d-a,d+a), + r = function(r, d) stats::dunif(r,d-a,d+a)/(2*pi*r), + d = function(r, d) stats::dunif(r,d-a,d+a)/(2*pi*d) ) }, gaussian = { rmax_bw <- Inf switch(divisor, - function(r, d) dnorm((r^2-d^2)*pi,0,sd=bw), - none = function(r, d) dnorm(r,d,sd = bw), - r = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*r), - d = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*d) + function(r, d) stats::dnorm((r^2-d^2)*pi,0,sd=bw), + none = function(r, d) stats::dnorm(r,d,sd = bw), + r = function(r, d) stats::dnorm(r,d,sd = bw)/ (2*pi*r), + d = function(r, d) stats::dnorm(r,d,sd = bw)/ (2*pi*d) ) }, @@ -340,7 +338,7 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, # If separate summary statistics (Dk, K, Hs, pcf) are to be included in the energy calculation, the "spatstat" package is loaded here and an error message is displayed if it is not installed. if (length(w_statistics)) { - if("spatstat" %in% rownames(installed.packages()) == FALSE) { + if("spatstat" %in% rownames(utils::installed.packages()) == FALSE) { stop("the package 'spatstat' must be installed on your computer for the application if the w_statistics option is used.", call. = FALSE) } @@ -371,10 +369,10 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, xwr <- obs_window$xrange ywr <- obs_window$yrange p <- p_[sample.int(nrow(p_),ceiling(nrow(p_)*((diff(xwr)*diff(ywr))/(diff(xr)*diff(yr)))),TRUE), ]# rpois - p$x <- runif(nrow(p), xwr[1], xwr[2]) - p$y <- runif(nrow(p), ywr[1], ywr[2]) - p$mark[, diameter] <- quantile(p_$mark[, diameter], - probs = runif(nrow(p), 0, 1), type = 4) + p$x <- stats::runif(nrow(p), xwr[1], xwr[2]) + p$y <- stats::runif(nrow(p), ywr[1], ywr[2]) + p$mark[, diameter] <- stats::quantile(p_$mark[, diameter], + probs = stats::runif(nrow(p), 0, 1), type = 4) p$mark[, species] <- p_$mark[, species, drop = FALSE][ sample.int(nrow(p_), nrow(p), replace = TRUE),, drop = FALSE] @@ -390,21 +388,21 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, # Prepare the graphical output. if(show_graphic == TRUE) { - par(mfrow = 1:2) - num_species <- from.dummy (p_$mark[, species, drop = FALSE]) + graphics::par(mfrow = 1:2) + num_species <- from_dummy (p_$mark[, species, drop = FALSE]) plot(y~x, p_, pch=19, col= 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter], xlim = xr, ylim = yr, xaxs ="i", yaxs ="i", main ="Reference", xlab ="x [m]", ylab ="y [m]") - text(p_$x, p_$y, as.integer(num_species), cex=0.7) + graphics::text(p_$x, p_$y, as.integer(num_species), cex=0.7) - plot(y~x, p, type = "n", + graphics::plot(y~x, p, type = "n", xlim = xwr, ylim = ywr, xaxs = "i", yaxs = "i", main = "Reconstructed", xlab = "x [m]", ylab = "y [m]") - clip(xwr[1], xwr[2], ywr[1], ywr[2]) + graphics::clip(xwr[1], xwr[2], ywr[1], ywr[2]) } # Show warning if certain distances between pairs of trees are not present. - energy <- Energy_fun(f, f0, statistics, f_, f0_, statistics_)["energy"] + energy <- energy_fun(f, f0, statistics, f_, f0_, statistics_)["energy"] # Prepares variables for the storage of progress. energy_launch <- as.vector(energy) @@ -427,11 +425,11 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, # Updating the graphical output of all "issue" steps. if (step %% issue == 0) { if(show_graphic == TRUE) { - rect(xwr[1], ywr[1], xwr[2], ywr[2], col="white") - num_species <- from.dummy (p$mark[, species, drop = FALSE]) + graphics::rect(xwr[1], ywr[1], xwr[2], ywr[2], col="white") + num_species <- from_dummy (p$mark[, species, drop = FALSE]) - points(y~x, p, pch = 19, col = 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter]) - text(p$x, p$y, as.integer(num_species), cex = 0.7) + graphics::points(y~x, p, pch = 19, col = 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter]) + graphics::text(p$x, p$y, as.integer(num_species), cex = 0.7) } # Generates text output with the current calculated values (for example the energy), this is updated every "issue" simulation step. @@ -441,7 +439,7 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, energy_improvement, "\t\t\r", appendLF = FALSE) Sys.sleep(0) - flush.console() + utils::flush.console() } # the next code block aborts the reconstruction if the energy decreases by less than tol in "no_changes" intervals of steps_tol simulation steps. if (step %% steps_tol == 0) { @@ -475,8 +473,8 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, move_coordinate = { i <- sample.int(nrow(p), 1, replace = TRUE) s <- if(sd=="step") nrow(p) * 1 / step else sd - x <- rnorm(1, p$x[i], diff(xwr) * s) %% xwr[2] - y <- rnorm(1, p$y[i], diff(ywr) * s) %% ywr[2] + x <- stats::rnorm(1, p$x[i], diff(xwr) * s) %% xwr[2] + y <- stats::rnorm(1, p$y[i], diff(ywr) * s) %% ywr[2] mdiff <- p$mark[i, ] f.new <- f - calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, kernel, rmax_bw, r) + @@ -542,7 +540,7 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, pick_mark_one = { i <- sample.int(nrow(p), 1, replace = TRUE) m <- p$mark[i, ] - m[diameter] <-quantile(p_$mark[,diameter],probs = runif(1,0,1), + m[diameter] <- stats::quantile(p_$mark[,diameter],probs = stats::runif(1,0,1), type = 4) mdiff <- m - p$mark[i, ] f.new <- f + calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, @@ -566,7 +564,7 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, stop("undefined case") ) # calculates the energy - Energy <- Energy_fun(f.new, f0.new, statistics.new, f_, f0_, statistics_) + Energy <- energy_fun(f.new, f0.new, statistics.new, f_, f0_, statistics_) energy.new<-Energy[["energy"]] # Sets the currently calculated energy as the new reference value if it is less than the previous energy. @@ -614,8 +612,8 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, message("\n") # Saves all results Transfers them to the "reconstruction" list. - p_$species <- from.dummy(p_$mark[, species, drop = FALSE]) - p$species <- from.dummy(p$mark[, species, drop = FALSE]) + p_$species <- from_dummy(p_$mark[, species, drop = FALSE]) + p$species <- from_dummy(p$mark[, species, drop = FALSE]) method <- "Reconstruction of a homogeneous point pattern" Parameter_setting <- list(n_repetitions=n_repetitions, max_steps=max_steps, no_changes=no_changes, rcount=rcount, rmax=rmax, diff --git a/man/reconstruct_pattern_multi.Rd b/man/reconstruct_pattern_multi.Rd new file mode 100644 index 00000000..73014754 --- /dev/null +++ b/man/reconstruct_pattern_multi.Rd @@ -0,0 +1,188 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/reconstruct_pattern_multi.R +\name{reconstruct_pattern_multi} +\alias{reconstruct_pattern_multi} +\title{reconstruct_pattern_multi} +\usage{ +reconstruct_pattern_multi( + marked_pattern, + xr = marked_pattern$window$xrange, + yr = marked_pattern$window$yrange, + n_repetitions = 1, + max_steps = 10000, + no_changes = 5, + rcount = 250, + rmax = 25, + issue = 1000, + divisor = "r", + kernel_arg = "epanechnikov", + timing = FALSE, + energy_evaluation = FALSE, + show_graphic = FALSE, + Lp = 1, + bw = if (divisor \%in\% c("r", "d")) 0.5 else 5, + sd = "step", + steps_tol = 1000, + tol = 1e-04, + w_markcorr = c(d_d = 1, all = 1, d_all = 1, all_all = 1, d_d0 = 1, all0 = 1, d_all0 = + 1, all_all0 = 1), + prob_of_actions = c(move_coordinate = 0.4, switch_coords = 0.1, exchange_mark_one = + 0.1, exchange_mark_two = 0.1, pick_mark_one = 0.2, pick_mark_two = 0.1), + k = 1, + w_statistics = c() +) +} +\arguments{ +\item{marked_pattern}{ppp object with marked pattern. See Details section +for more information.} + +\item{xr, yr}{Add description} + +\item{n_repetitions}{Integer representing the number of simulations to be +performed.} + +\item{max_steps}{Maximum number simulation steps.} + +\item{no_changes}{Integer representing the number of iterations +(per 1000 simulation steps) +after which the reconstruction is terminated if the energy does not decrease.} + +\item{rcount}{Integer representing the number of intervals for which the +summary statistics are evaluated.} + +\item{rmax}{Maximum distance [m] at which the summary statistics are +evaluated.} + +\item{issue}{Integer that determines after how many simulations steps an +output occurs.} + +\item{divisor}{Divisor in the smoothing kernel, d or r.} + +\item{kernel_arg}{The kernel used to calculate the energy, possible kernels +can be: Gaussian, Epanechnikov, Rectangular, Cumulative.} + +\item{timing}{Logical value: The computation timeis measured if this is TRUE} + +\item{energy_evaluation}{Logical value: If this is TRUE, the procedure stores +the energy shares of the total energy per simulation step.} + +\item{show_graphic}{Logical value: If this is TRUE, the procedure records the +point pattern during optimisation and updated} + +\item{Lp}{Distance measure for the calculation of the energy function +(Lp distance, 1 ≤ p < ∞).} + +\item{bw}{Bandwidth [m] with which the kernels are scaled, so that this is +the standard deviation of the smoothing kernel.} + +\item{sd}{This is the standard deviation [m] used in the move_coordinate +action.} + +\item{steps_tol}{After the value steps_tol it is checked whether the energy +change is smaller than tol.} + +\item{tol}{Tolerance: stops the procedure of energy more than 1 − tol +times no_changes.} + +\item{w_markcorr}{Vector of possible weightings of individual mcf's. +\code{c(1, 1, 1, 1, 1, 1, 1, 1)}} + +\item{prob_of_actions}{Vector of probabilities for the actions performed. +\code{c(move_coordinate = 0.4, switch_coords = 0.1, exchange_mark_one = 0.1, +exchange_mark_two = 0.1, pick_mark_one = 0.2, pick_mark_two = 0.1)}} + +\item{k}{Vector of values k; used only if Dk is included in w_statistics +elow.} + +\item{w_statistics}{vector of named weights for optional spatial statistics +from the ‘spatstat’package to be included in the energy calculation.This may +include Dk, K, Hs, pcf.} +} +\value{ +rd_mar +} +\description{ +Pattern reconstruction of a pattern marked by multiple traits. +} +\details{ +A novel approach carries out a pattern reconstruction of marked dot patterns +as described by Tscheschel and Stoyan (2006) and Wiegand and Moloney (2014). + +One particular feature is the simultaneous consideration of both marks, +accounting for their correlation during reconstruction. + +The marked point pattern (PPP object) must is currently structured as follows: +X-coordinate, Y-coordinate, metric mark (e.g. diameter at breast height), +and nominal mark (e.g. tree species).It is calculated in the unit metre [m]. + +A combination of the mark correlation function and pair correlation function +is used for pattern description. Additional summary statistics may be +considered.Two randomly selected marks are chosen in each iteration, and one +of various actions is performed. Changes will only be retained if the +difference between the observed and reconstructed pattern decreases +(minimizing energy). + +This method is currently only suitable for homogeneous point patterns. + +A comprehensive description of the method can be +found in Wudel et al. (2023). Future developments will involve considering +more than two markers, expanding and reducing the observation window under +using an Edge Correction. Additionally, data input will be automatically +detected. +} +\examples{ +\dontrun{ + +# Random example data set + xr <- 500 + yr <- 1000 + N <- 400 + y <- runif(N, min = 0, max = yr) + x <- runif(N, min = 0, max = xr) + species <- as.factor(c("A","B")) + species <- sample(species, N, replace = TRUE) + diameter <- runif(N, 1, 40) + plot(x,y, col = species, pch = 19, cex = diameter*0.02) + random <- data.frame(x, y, diameter/ 0.1,species) + colnames(random) <- c("x", "y", "dbh [mm]", "Tree species") + +# Conversion to a ppp object and conversion of the metric mark to metres. + w <- owin(c(0, 500), c(0, 1000)) + marked_pattern <- as.ppp(data.frame(random), W = w) + marked_pattern$marks$dbh..mm.<-marked_pattern$marks$dbh..mm.*0.001 + +# Reconstruction function +reconstruction <- reconstruct_pattern_multi( + marked_pattern, n_repetitions = 1, max_steps = 100000, no_changes = 5, + rcount = 250, rmax = 25, issue = 1000, divisor = "r", + kernel_arg = "epanechnikov", timing = TRUE, energy_evaluation = TRUE, + show_graphic = TRUE, Lp = 1, bw = 0.5, sd = "step", steps_tol = 1000, + tol = 1e-4, w_markcorr = c(d_d=1, all=1, d_all=1, all_all=1, d_d0=1, + all0=1, d_all0=1, all_all0=1), prob_of_actions = c(move_coordinate = 0.4, + switch_coords = 0.1, exchange_mark_one = 0.1, exchange_mark_two = 0.1, + pick_mark_one = 0.2, pick_mark_two = 0.1), + k = 1, w_statistics = c() +) +} + +} +\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 + +Wudel, C., Schlicht, R., & Berger, U. (2023). Multi-trait point pattern +reconstruction of plant ecosystems. Methods in Ecology and Evolution, 14, 2668–2679. +https://doi.org/10.1111/2041-210X.14206 +} +\seealso{ +\code{\link{fit_point_process}} \cr +\code{\link{reconstruct_pattern}} \cr +\code{\link{reconstruct_pattern_marks}} +} diff --git a/man/shar.Rd b/man/shar.Rd index 5d64a00e..668db0da 100644 --- a/man/shar.Rd +++ b/man/shar.Rd @@ -26,6 +26,7 @@ Useful links: Authors: \itemize{ \item Marco Sciaini \email{marco.sciaini@posteo.net} (\href{https://orcid.org/0000-0002-3042-5435}{ORCID}) + \item Chris Wudel \email{chris.wudel@tu-dresden.de} (\href{https://orcid.org/0000-0003-0446-4665}{ORCID}) } Other contributors: From a9f1b85b50ef856c51ad9dfa7cfb26924ff3b46f Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 24 Nov 2023 09:07:36 +0100 Subject: [PATCH 05/39] :bug: Bug-fix in reconstruction --- R/reconstruct_algorithm.R | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/R/reconstruct_algorithm.R b/R/reconstruct_algorithm.R index 20727bf6..31416fa4 100644 --- a/R/reconstruct_algorithm.R +++ b/R/reconstruct_algorithm.R @@ -60,9 +60,7 @@ reconstruct_algorithm <- function(pattern, window <- pattern$window # check if pattern is emtpy - if (n_points == 0){ - stop("The observed pattern contains no points.", call. = FALSE) - } + if (n_points == 0) stop("The observed pattern contains no points.", call. = FALSE) # calculate intensity lambda <- n_points / spatstat.geom::area(window) @@ -145,23 +143,27 @@ reconstruct_algorithm <- function(pattern, window = window, verbose = FALSE) # remove points because more points in simulated - if (n_points <= simulated$n) { + 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 { + # add points because less points in simulated + } else if (n_points > simulated$n) { # 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) @@ -172,6 +174,9 @@ reconstruct_algorithm <- function(pattern, # check if simulated is empty if (simulated$n == 0) stop("The simulated pattern contains no points.", call. = FALSE) + # check if simulated has same points as observed + if (n_points != simulated$n) stop("The simulated pattern contains the same amount as observed pattern.", call. = FALSE) + # energy before reconstruction energy <- Inf From 2ce7d18773d7d9df3389d160b522b70735d38838 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 24 Nov 2023 10:34:56 +0100 Subject: [PATCH 06/39] Fix R-CMD-Checks --- R/reconstruct_pattern_multi.R | 19 +++++-------------- R/shar-package.R | 2 ++ 2 files changed, 7 insertions(+), 14 deletions(-) diff --git a/R/reconstruct_pattern_multi.R b/R/reconstruct_pattern_multi.R index 11acc36b..6f969044 100644 --- a/R/reconstruct_pattern_multi.R +++ b/R/reconstruct_pattern_multi.R @@ -208,22 +208,22 @@ reconstruct_pattern_multi <- function(marked_pattern, # Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. Dk = { - nnd_ <- as.matrix(nndist(x, y, k=k)) - apply(nnd_, 2, function(z) cumsum(hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) + nnd_ <- as.matrix(spatstat.geom::nndist(x, y, k=k)) + apply(nnd_, 2, function(z) cumsum(graphics::hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) }, # Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. K = { - kest<-Kest(ppp(x,y,window=owin(xr,yr)), rmax=rmax, correction="none")# , breaks = c(-Inf, r) + kest<-spatstat.explore::Kest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), rmax=rmax, correction="none")# , breaks = c(-Inf, r) kest$un }, # Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. pcf = { - pcfest<-pcf(ppp(x,y,window=owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") + pcfest<-spatstat.explore::pcf(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") pcfest$un }, # Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. Hs = { - hest<-Hest(ppp(x,y,window=owin(xr,yr)), correction="none") #, breaks = c(-Inf, r) + hest<-spatstat.explore::Hest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), correction="none") #, breaks = c(-Inf, r) hest$raw }, stop("unknown statistic") @@ -336,15 +336,6 @@ reconstruct_pattern_multi <- function(marked_pattern, } ) -# If separate summary statistics (Dk, K, Hs, pcf) are to be included in the energy calculation, the "spatstat" package is loaded here and an error message is displayed if it is not installed. - if (length(w_statistics)) { - if("spatstat" %in% rownames(utils::installed.packages()) == FALSE) { - stop("the package 'spatstat' must be installed on your computer for the - application if the w_statistics option is used.", call. = FALSE) - } - require (spatstat) - } - # Determination of the weightings of the mark correlation functions. fn <- list() for (i in seq_along(marknames)) for (j in seq_along(marknames)) if (i <= j) { diff --git a/R/shar-package.R b/R/shar-package.R index fe37eaa9..f490c5d4 100644 --- a/R/shar-package.R +++ b/R/shar-package.R @@ -18,6 +18,8 @@ globalVariables(c("count", "hi", "iso", "lo", + "mark", + "obs_window", "r", "summary_function", "theo", From 7bcc6727794207fba5b7572df7db5f9bca0f05e9 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 1 Dec 2023 15:14:19 +0100 Subject: [PATCH 07/39] Update citation everywhere --- R/data.R | 2 +- README.Rmd | 6 +++++- inst/CITATION | 2 +- man/landscape.Rd | 2 +- vignettes/articles/publication_record.Rmd | 8 +++++--- vignettes/get_started.Rmd | 6 +++++- 6 files changed, 18 insertions(+), 8 deletions(-) diff --git a/R/data.R b/R/data.R index 18470c86..0a27607e 100644 --- a/R/data.R +++ b/R/data.R @@ -4,7 +4,7 @@ #' generated with the \code{NLMR::nlm_fbm()} algorithm. #' #' @format A SpatRaster object. -#' @source Simulated neutral landscape model with R. https://github.com/ropensci/NLMR/ +#' @source Simulated neutral landscape model with R. "landscape" #' Species a diff --git a/README.Rmd b/README.Rmd index cd305831..be19212b 100644 --- a/README.Rmd +++ b/README.Rmd @@ -55,7 +55,11 @@ Please refer to `vignette("Get started")` and the [homepage](https://r-spatialec 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 +> 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. . + +If you use the `reconstruct_pattern_multi` function, please also cite. + +> Wudel C., Schlicht R., Berger U. (2023). Multi-trait point pattern reconstruction of plant ecosystems. Methods in Ecology and Evolution, 14, 2668–2679. . The get a BibTex entry, please use `citation("shar")`. diff --git a/inst/CITATION b/inst/CITATION index a17e8c59..fda66fbe 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -2,7 +2,7 @@ c( bibentry( bibtype = "Article", title = "shar: A R package to analyze species-habitat associations using point pattern analysis", - author = "Maximilian H.K. Hesselbarth", + author = person("Maximilian H.K.", "Hesselbarth"), journal = "Journal of Open Source Software", year = "2021", volume = "6", diff --git a/man/landscape.Rd b/man/landscape.Rd index fdcce8f9..037aa47c 100644 --- a/man/landscape.Rd +++ b/man/landscape.Rd @@ -8,7 +8,7 @@ A SpatRaster object. } \source{ -Simulated neutral landscape model with R. https://github.com/ropensci/NLMR/ +Simulated neutral landscape model with R. } \usage{ landscape diff --git a/vignettes/articles/publication_record.Rmd b/vignettes/articles/publication_record.Rmd index 8a3c6dc5..062eea5c 100644 --- a/vignettes/articles/publication_record.Rmd +++ b/vignettes/articles/publication_record.Rmd @@ -8,14 +8,16 @@ vignette: > %\VignetteEncoding{UTF-8} --- -We try to maintain a list of publications that have used the **shar** package. If you have used **shar** in your own work, we would love to hear from you. Please [open an issue on GitHub](https://github.com/r-spatialecology/shar/issues) or [drop us an e-mail](mailto:mhk.hesselbarth@gmail.com) so we can add your work to this list. - 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 +If you use the `reconstruct_pattern_multi` function, please also cite. + +> Wudel C., Schlicht R., Berger U. (2023). Multi-trait point pattern reconstruction of plant ecosystems. Methods in Ecology and Evolution, 14, 2668–2679. . + The get a BibTex entry, please use `citation("shar")`. ## A list of publications -...Work-in-Progress... +To see all papers that already used **landscapemetrics**, please [click here](https://scholar.google.com/scholar?cites=1225822175035945633&as_sdt=2005&sciodt=0,5&hl=en) diff --git a/vignettes/get_started.Rmd b/vignettes/get_started.Rmd index 962a4520..db04bb6f 100644 --- a/vignettes/get_started.Rmd +++ b/vignettes/get_started.Rmd @@ -112,7 +112,11 @@ There are many more functions, which can be found [here](https://r-spatialecolog 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 +> 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. . + +If you use the `reconstruct_pattern_multi` function, please also cite. + +> Wudel C., Schlicht R., Berger U. (2023). Multi-trait point pattern reconstruction of plant ecosystems. Methods in Ecology and Evolution, 14, 2668–2679. . The get a BibTex entry, please use `citation("shar")`. From 8d80e72f6e52f3a9b04755da9ba5046b0a38495b Mon Sep 17 00:00:00 2001 From: ChrisWudel Date: Tue, 12 Dec 2023 13:54:49 +0100 Subject: [PATCH 08/39] outsourcing of auxiliary functions --- R/calc_moments.R | 64 ++++++++++ R/compute_statistics.R | 54 ++++++++ R/dummy_transf.R | 31 +++++ R/energy_fun.R | 59 +++++++++ R/reconstruct_pattern_multi_trait_marks.R | 145 ++-------------------- R/select_kernel.R | 76 ++++++++++++ 6 files changed, 293 insertions(+), 136 deletions(-) create mode 100644 R/calc_moments.R create mode 100644 R/compute_statistics.R create mode 100644 R/dummy_transf.R create mode 100644 R/energy_fun.R create mode 100644 R/select_kernel.R diff --git a/R/calc_moments.R b/R/calc_moments.R new file mode 100644 index 00000000..36c6d282 --- /dev/null +++ b/R/calc_moments.R @@ -0,0 +1,64 @@ +#' calc_moments +#' +#' @description calc_moments (internal) +#' +#' @param fn Determination of the weightings of the mark correlation functions. +#' @param p Defines the initial state of the new ponit pattern. +#' @param x x Coordinate of the points from the reference point pattern. +#' @param y y Coordinate of the points from the reference point pattern. +#' @param mark Marks the currently viewed point pattern. +#' @param kernel Result of the kernel calculation, calculated with the +#' calc_kernels function. +#' @param rmax_bw Maximum distance at which the summary statistics are +#' evaluated + Bandwidth with which the kernels are scaled, so that this is the +#' standard deviation of the smoothing kernel. +#' @param r is a sequence from rmin to rmax in rcount steps. +#' @param exclude +#' +#' @details +#' Definition of the product-moment function for calculating the contribution +#' of a point at the coordinates x, y with marking. +#' +#' @return matrix +#' +#' @aliases calc_moments +#' @rdname calc_moments +#' +#' @aliases reconstruct_algorithm +#' @rdname reconstruct_algorithm +#' +#' @keywords internal +#' +#' +calc_moments <- function(fn, + p, + exclude=NULL, + x, + y, + mark, + kernel, + rmax_bw, + r) + { + d2 <- (p$x-x)^2 + (p$y-y)^2 + use <- d2 <= rmax_bw^2 + use[exclude] <- FALSE + z <- crossprod(p$mark[use, , drop = FALSE], + outer(sqrt(d2[use]), r, function(d, r) kernel(r, d))) + z[fn$i, , drop = FALSE] * mark[fn$j] + z[fn$j, , drop = FALSE] * mark[fn$i] + } + +calc_moments_full <- function(fn, + p, + kernel, + rmax_bw, + r) + { + f <- 0 + for (i in seq_len(nrow(p))) { + f <- f + calc_moments(fn, p, i:nrow(p), p$x[i], p$y[i], p$mark[i, ], + kernel, rmax_bw, r) + } + rownames(f) <- paste(names(fn$i), names(fn$j), sep = ":") + f + } diff --git a/R/compute_statistics.R b/R/compute_statistics.R new file mode 100644 index 00000000..5a08ee63 --- /dev/null +++ b/R/compute_statistics.R @@ -0,0 +1,54 @@ +#' compute_statistics +#' +#' @description compute_statistics (internal) +#' +#' @param x x Coordinate of the points from the reference point pattern. +#' @param y y Coordinate of the points from the reference point pattern. +#' @param k Vector of values k; used only if Dk is included in w_statistics +#' elow. +#' @param xr x extension of the observation window (start, end) +#' @param yr y extension of the observation window (start, end) +#' @param w_statistics +#' +#' @details +#' Compute optional spatial statistics using the spatstat package. +#' +#' @return list +#' +#' @aliases compute_statistics +#' @rdname compute_statistics +#' +#' @keywords internal +#' + +compute_statistics <- function(x, y, k, xr, yr, w_statistics) { + stat <- names(w_statistics) + names(stat) <- stat + lapply(stat, function(name) switch(name, + # Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. + + Dk = { + nnd_ <- as.matrix(nndist(x, y, k=k)) + apply(nnd_, 2, function(z) cumsum(hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) + }, + # Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. + + K = { + kest<-Kest(ppp(x,y,window=owin(xr,yr)), rmax=rmax, correction="none") + kest$un + }, + # Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. + + pcf = { + pcfest<-pcf(ppp(x,y,window=owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") + pcfest$un + }, + # Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. + + Hs = { + hest<-Hest(ppp(x,y,window=owin(xr,yr)), correction="none") + hest$raw + }, + stop("unknown statistic") + )) + } diff --git a/R/dummy_transf.R b/R/dummy_transf.R new file mode 100644 index 00000000..c85cddd0 --- /dev/null +++ b/R/dummy_transf.R @@ -0,0 +1,31 @@ +#' dummy_transf +#' +#' @description dummy_transf (internal) +#' +#' @param f Result of the calc_moments_full function which represents +#' product-moment contribution of a point at coordinates x, y with marks, +#' for the whole point pattern. +#' +#' @details +#' Function for the transformation of variables to dummy variables and back +#' +#' @return matrix +#' +#' @aliases dummy_transf +#' @rdname dummy_transf +#' +#' @keywords internal +#' +to.dummy <- function(f) { + x <- matrix(0, length(f), nlevels(f), dimnames=list(names(f), levels(f))) + x[cbind(seq_along(f), as.integer(f))] <- 1 + x + } + +from.dummy <- function(x, levels=colnames(x)) { + f <- as.integer(x %*% seq_along(levels)) + levels(f) <- levels + class(f) <- "factor" + f + } + diff --git a/R/energy_fun.R b/R/energy_fun.R new file mode 100644 index 00000000..19b23716 --- /dev/null +++ b/R/energy_fun.R @@ -0,0 +1,59 @@ +#' energy_fun +#' +#' @description energy_fun (internal) +#' +#' @param f_ Result of the calc_moments_full function which represents +#' product-moment contribution of a point at coordinates x, y with marks, +#' for the whole reference point pattern. +#' @param f0_ Column sums of the weights of the brand correlation functions of +#' the reference point pattern. +#' @param statistics_ Results of the compute_statistics function for the +#' reference point pattern (calculation of optional spatial statistics). +#' @details +#' Defining the Energy_fun function to calculate the "energy" of the pattern +#' (where a lower energy indicates a better match). +#' @param f Result of the calc_moments_full function which represents +#' product-moment contribution of a point at coordinates x, y with marks, +#' for the whole new ponit pattern. +#' @param f0 Column sums of the weights of the brand correlation functions of +#' the new point pattern. +#' @param statistics Results of the compute_statistics function for the +#' new point pattern (calculation of optional spatial statistics). +#' @param fn Determination of the weightings of the mark correlation functions. +#' @param p Defines the initial state of the new ponit pattern. +#' @param p_ Reference point pattern. +#' @param Lp Distance measure for the calculation of the energy function +#' (Lp distance, 1 ≤ p < ∞). +#' @param w_statistics ector of named weights for optional spatial statistics +#' from the ‘spatstat’package to be included in the energy calculation.This may +#' include Dk, K, Hs, pcf. +#' +#' @details +#' Defining the Energy_fun function to calculate the "energy" of the pattern +#' (where a lower energy indicates a better match). +#' +#' @return vector +#' +#' @aliases energy_fun +#' @rdname energy_fun +#' +#' @keywords internal +#' +Energy_fun <- function(f, f0, statistics, f_, f0_, statistics_, fn, p, p_, Lp, w_statistics) { + result <- c( + f = sum(fn$w * rowMeans(abs( + f / nrow(p) - + f_ / nrow(p_) + )^Lp)), + f0 = sum(fn$w0 * abs( + f0 / nrow(p) - + f0_ / nrow(p_) + )^Lp), + if (length(w_statistics)) + sapply(seq_along(w_statistics), function(i) w_statistics[i] * + mean(abs(statistics[[i]] - statistics_[[i]])^Lp, na.rm = TRUE), + USE.NAMES=FALSE + ) + ) + c(energy = sum(result), result) +} diff --git a/R/reconstruct_pattern_multi_trait_marks.R b/R/reconstruct_pattern_multi_trait_marks.R index 44555340..a5ce1006 100644 --- a/R/reconstruct_pattern_multi_trait_marks.R +++ b/R/reconstruct_pattern_multi_trait_marks.R @@ -136,6 +136,7 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, xr = marked_pattern$window$xrange, yr = marked_pattern$window$yrange, + obs_window = owin(c(xr),c(yr)), n_repetitions = 1, max_steps = 10000, no_changes = 5, @@ -168,90 +169,6 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, # Loop to perform multiple reconstructions. for (t in seq_len(n_repetitions)) { -# Definition of the product-moment function for calculating the contribution of a point at the coordinates x, y with marking. - calc_moments <- function(fn, p, exclude=NULL, x, y, mark, kernel, rmax_bw, r) { - d2 <- (p$x-x)^2 + (p$y-y)^2 - use <- d2 <= rmax_bw^2 - use[exclude] <- FALSE - z <- crossprod(p$mark[use, , drop = FALSE], - outer(sqrt(d2[use]), r, function(d, r) kernel(r, d))) - z[fn$i, , drop = FALSE] * mark[fn$j] + z[fn$j, , drop = FALSE] * mark[fn$i] - } - -# Definition of the calc_moments_full function to calculate the calc_moments function for the whole pattern. - calc_moments_full <- function(fn, p, kernel, rmax_bw, r) { - f <- 0 - for (i in seq_len(nrow(p))) { - f <- f + calc_moments(fn, p, i:nrow(p), p$x[i], p$y[i], p$mark[i, ], - kernel, rmax_bw, r) - } - rownames(f) <- paste(names(fn$i), names(fn$j), sep = ":") - f - } - -# Function for the transformation of variables to dummy variables and back - to.dummy <- function(f) { - x <- matrix(0, length(f), nlevels(f), dimnames=list(names(f), levels(f))) - x[cbind(seq_along(f), as.integer(f))] <- 1 - x - } - from.dummy <- function(x, levels=colnames(x)) { - f <- as.integer(x %*% seq_along(levels)) - levels(f) <- levels - class(f) <- "factor" - f - } - -# Compute optional spatial statistics using the spatstat package. - compute_statistics <- function(x, y, k, xr, yr) { - stat <- names(w_statistics) - names(stat) <- stat - lapply(stat, function(name) switch(name, - -# Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. - Dk = { - nnd_ <- as.matrix(nndist(x, y, k=k)) - apply(nnd_, 2, function(z) cumsum(hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) - }, -# Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. - K = { - kest<-Kest(ppp(x,y,window=owin(xr,yr)), rmax=rmax, correction="none")# , breaks = c(-Inf, r) - kest$un - }, -# Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. - pcf = { - pcfest<-pcf(ppp(x,y,window=owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") - pcfest$un - }, -# Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. - Hs = { - hest<-Hest(ppp(x,y,window=owin(xr,yr)), correction="none") #, breaks = c(-Inf, r) - hest$raw - }, - stop("unknown statistic") - )) - } - -# Defining the Energy_fun function to calculate the "energy" of the pattern (where a lower energy indicates a better match). - Energy_fun <- function(f, f0, statistics, f_, f0_, statistics_) { - result <- c( - f = sum(fn$w * rowMeans(abs( - f / nrow(p) - - f_ / nrow(p_) - )^Lp)), - f0 = sum(fn$w0 * abs( - f0 / nrow(p) - - f0_ / nrow(p_) - )^Lp), - if (length(w_statistics)) - sapply(seq_along(w_statistics), function(i) w_statistics[i] * - mean(abs(statistics[[i]] - statistics_[[i]])^Lp, na.rm = TRUE), - USE.NAMES=FALSE - ) - ) - c(energy = sum(result), result) - } - # Load the reference point pattern as a data frame with the components x, y, mark, where x, y are the coordinates of the point and mark is a matrix representing the marks or their dummy values. p_ <- data.frame( x = marked_pattern$x, @@ -290,53 +207,9 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, r <- seq(rmin, rmax, , rcount) # Calculation of the kernels. - kernel <- switch(kernel_arg, - epanechnikov = { - a <- bw * sqrt(5) - rmax_bw <- rmax + a - switch(divisor, - { - rmax_bw <- sqrt(rmax^2 + a/pi) - function(r, d) pmax.int(0, 1 - ((r^2-d^2)*pi/a)^2) * 0.75/a - }, - none = function(r, d) pmax.int(0, 1 - ((r-d)/a)^2) * 0.75/a, - r = function(r, d) pmax.int(0, 1 - ((r-d)/a)^2) * 0.75/(a*2*pi*r), - d = function(r, d) pmax.int(0, 1 - ((r-d)/a)^2) * 0.75/(a*2*pi*d) - ) - }, - rectangular =, box = { - a <- bw * sqrt(3) - rmax_bw <- rmax + a - switch(divisor, - { - rmax_bw <- sqrt(rmax^2 + a/pi) - function(r, d) dunif((r^2-d^2)*pi,-a,+a) - }, - none = function(r, d) dunif(r,d-a,d+a), - r = function(r, d) dunif(r,d-a,d+a)/(2*pi*r), - d = function(r, d) dunif(r,d-a,d+a)/(2*pi*d) - ) - }, - gaussian = { - rmax_bw <- Inf - switch(divisor, - function(r, d) dnorm((r^2-d^2)*pi,0,sd=bw), - none = function(r, d) dnorm(r,d,sd = bw), - r = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*r), - d = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*d) - ) - }, - - cumulative = { - rmax_bw <- rmax - switch(divisor, - function(r, d) as.numeric(d <= r), - none = function(r, d) as.numeric(d <= r), - r = function(r, d) (d <= r) / (2*pi*r), - d = function(r, d) (d <= r) / (2*pi*d) - ) - } - ) + sel_kernel <- select_kernel(kernel_arg, bw, rmax, divisor) + kernel <- as.function(sel_kernel[[1]]) + rmax_bw <- as.numeric(sel_kernel[[2]]) # If separate summary statistics (Dk, K, Hs, pcf) are to be included in the energy calculation, the "spatstat" package is loaded here and an error message is displayed if it is not installed. if (length(w_statistics)) { @@ -382,11 +255,11 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, f_ <- calc_moments_full(fn, p_, kernel, rmax_bw, r) f0_ <- colSums(p_$mark[, fn$i] * p_$mark[, fn$j]) names(f0_) <- rownames(f_) - statistics_<- compute_statistics(p_$x, p_$y, k, xr, yr) + statistics_<- compute_statistics(p_$x, p_$y, k, xr, yr, w_statistics) f <- calc_moments_full(fn, p, kernel, rmax_bw, r) f0 <- colSums(p$mark[, fn$i] * p$mark[, fn$j]) names(f0) <- rownames(f) - statistics <- compute_statistics(p$x, p$y, k, xwr, ywr) + statistics <- compute_statistics(p$x, p$y, k, xwr, ywr, w_statistics) # Prepare the graphical output. if(show_graphic == TRUE) { @@ -404,7 +277,7 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, } # Show warning if certain distances between pairs of trees are not present. - energy <- Energy_fun(f, f0, statistics, f_, f0_, statistics_)["energy"] + energy <- Energy_fun(f, f0, statistics, f_, f0_, statistics_, fn, p, p_, Lp, w_statistics)["energy"] # Prepares variables for the storage of progress. energy_launch <- as.vector(energy) @@ -482,7 +355,7 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, kernel, rmax_bw, r) + calc_moments(fn, p, i, x, y, mdiff, kernel, rmax_bw, r) f0.new<- f0 - statistics.new <- compute_statistics(replace(p$x, i, x), replace(p$y, i, y), k, xwr, ywr) + statistics.new <- compute_statistics(replace(p$x, i, x), replace(p$y, i, y), k, xwr, ywr, w_statistics) }, # Swaps the coordinates of two randomly drawn points from the new point pattern, applied in xx% of the trap. @@ -566,7 +439,7 @@ reconstruct_pattern_multi_trait_marks <- function(marked_pattern, stop("undefined case") ) # calculates the energy - Energy <- Energy_fun(f.new, f0.new, statistics.new, f_, f0_, statistics_) + Energy <- Energy_fun(f.new, f0.new, statistics.new, f_, f0_, statistics_, fn, p, p_, Lp, w_statistics) energy.new<-Energy[["energy"]] # Sets the currently calculated energy as the new reference value if it is less than the previous energy. diff --git a/R/select_kernel.R b/R/select_kernel.R new file mode 100644 index 00000000..d14cea14 --- /dev/null +++ b/R/select_kernel.R @@ -0,0 +1,76 @@ +#' select_kernel +#' +#' @description select_kernel (internal) +#' +#' @param kernel_arg Parameter of the function +#' reconstruct_pattern_multi_trait_marks, specifies the kernel to be used to +#' calculate the energy, possible kernels can be: Gaussian, Epanechnikov, +#' Rectangular, Cumulative. +#' @param bw Bandwidth with which the kernels are scaled, so that this is the +#' standard deviation of the smoothing kernel. +#' @param rmax Maximum distance at which the summary statistics are evaluated. +#' @param divisor Divisor in the smoothing kernel, d or r. +#' +#' @details +#' Returns the function of the selected kernel, which is then used to +#' calculate the kernel. +#' +#' @return function +#' +#' @aliases select_kernel +#' @rdname select_kernel +#' +#' @keywords internal +#' +select_kernel <- function(kernel_arg, bw, rmax, divisor) { + kernel <- switch(kernel_arg, + epanechnikov = { + a <- bw * sqrt(5) + rmax_bw <- rmax + a + switch(divisor, + { + rmax_bw <- sqrt(rmax^2 + a/pi) + function(r, d) pmax.int(0, 1 - ((r^2-d^2)*pi/a)^2) * 0.75/a + }, + none = function(r, d) pmax.int(0, 1 - ((r-d)/a)^2) * 0.75/a, + r = function(r, d) pmax.int(0, 1 - ((r-d)/a)^2) * 0.75/(a*2*pi*r), + d = function(r, d) pmax.int(0, 1 - ((r-d)/a)^2) * 0.75/(a*2*pi*d) + ) + }, + rectangular =, box = { + a <- bw * sqrt(3) + rmax_bw <- rmax + a + switch(divisor, + { + rmax_bw <- sqrt(rmax^2 + a/pi) + function(r, d) dunif((r^2-d^2)*pi,-a,+a) + }, + none = function(r, d) dunif(r,d-a,d+a), + r = function(r, d) dunif(r,d-a,d+a)/(2*pi*r), + d = function(r, d) dunif(r,d-a,d+a)/(2*pi*d) + ) + }, + gaussian = { + rmax_bw <- Inf + switch(divisor, + function(r, d) dnorm((r^2-d^2)*pi,0,sd=bw), + none = function(r, d) dnorm(r,d,sd = bw), + r = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*r), + d = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*d) + ) + }, + + cumulative = { + rmax_bw <- rmax + switch(divisor, + function(r, d) as.numeric(d <= r), + none = function(r, d) as.numeric(d <= r), + r = function(r, d) (d <= r) / (2*pi*r), + d = function(r, d) (d <= r) / (2*pi*d) + ) + } + ) + sel_kernel <- list(kernel, rmax_bw) +} + + From 76e1cd5268046586e2904642cfd211f1de8a4c7b Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Thu, 14 Dec 2023 11:13:49 +0100 Subject: [PATCH 09/39] Fix couple of RCMD checks --- R/calc_moments.R | 5 --- R/compute_statistics.R | 12 +++--- R/dummy_transf.R | 4 +- R/reconstruct_pattern_multi.R | 4 +- R/select_kernel.R | 16 ++++---- man/calc_moments.Rd | 41 +++++++++++++++++++++ man/compute_statistics.Rd | 32 ++++++++++++++++ man/dummy_transf.Rd | 24 ++++++++++++ man/energy_fun.Rd | 69 +++++++++++++++++++++++++++++++++++ man/select_kernel.Rd | 32 ++++++++++++++++ 10 files changed, 215 insertions(+), 24 deletions(-) create mode 100644 man/calc_moments.Rd create mode 100644 man/compute_statistics.Rd create mode 100644 man/dummy_transf.Rd create mode 100644 man/energy_fun.Rd create mode 100644 man/select_kernel.Rd diff --git a/R/calc_moments.R b/R/calc_moments.R index 36c6d282..1d5ac81b 100644 --- a/R/calc_moments.R +++ b/R/calc_moments.R @@ -24,12 +24,7 @@ #' @aliases calc_moments #' @rdname calc_moments #' -#' @aliases reconstruct_algorithm -#' @rdname reconstruct_algorithm -#' #' @keywords internal -#' -#' calc_moments <- function(fn, p, exclude=NULL, diff --git a/R/compute_statistics.R b/R/compute_statistics.R index 5a08ee63..7657e17f 100644 --- a/R/compute_statistics.R +++ b/R/compute_statistics.R @@ -19,8 +19,6 @@ #' @rdname compute_statistics #' #' @keywords internal -#' - compute_statistics <- function(x, y, k, xr, yr, w_statistics) { stat <- names(w_statistics) names(stat) <- stat @@ -28,25 +26,25 @@ compute_statistics <- function(x, y, k, xr, yr, w_statistics) { # Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. Dk = { - nnd_ <- as.matrix(nndist(x, y, k=k)) - apply(nnd_, 2, function(z) cumsum(hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) + nnd_ <- as.matrix(spatstat.geom::nndist(x, y, k=k)) + apply(nnd_, 2, function(z) cumsum(graphics::hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) }, # Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. K = { - kest<-Kest(ppp(x,y,window=owin(xr,yr)), rmax=rmax, correction="none") + kest<-spatstat.explore::Kest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), rmax=rmax, correction="none") kest$un }, # Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. pcf = { - pcfest<-pcf(ppp(x,y,window=owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") + pcfest<-spatstat.explore::pcf(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") pcfest$un }, # Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. Hs = { - hest<-Hest(ppp(x,y,window=owin(xr,yr)), correction="none") + hest<-spatstat.explore::Hest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), correction="none") hest$raw }, stop("unknown statistic") diff --git a/R/dummy_transf.R b/R/dummy_transf.R index c85cddd0..ce823236 100644 --- a/R/dummy_transf.R +++ b/R/dummy_transf.R @@ -16,13 +16,13 @@ #' #' @keywords internal #' -to.dummy <- function(f) { +to_dummy <- function(f) { x <- matrix(0, length(f), nlevels(f), dimnames=list(names(f), levels(f))) x[cbind(seq_along(f), as.integer(f))] <- 1 x } -from.dummy <- function(x, levels=colnames(x)) { +from_dummy <- function(x, levels=colnames(x)) { f <- as.integer(x %*% seq_along(levels)) levels(f) <- levels class(f) <- "factor" diff --git a/R/reconstruct_pattern_multi.R b/R/reconstruct_pattern_multi.R index a8edf612..7554b1b7 100644 --- a/R/reconstruct_pattern_multi.R +++ b/R/reconstruct_pattern_multi.R @@ -252,7 +252,7 @@ reconstruct_pattern_multi <- function(marked_pattern, # Prepare the graphical output. if(show_graphic == TRUE) { graphics::par(mfrow = 1:2) - num_species <- from_dummy (p_$mark[, species, drop = FALSE]) + num_species <- from_dummy(p_$mark[, species, drop = FALSE]) plot(y~x, p_, pch=19, col= 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter], xlim = xr, ylim = yr, xaxs ="i", yaxs ="i", main ="Reference", xlab ="x [m]", ylab ="y [m]") @@ -289,7 +289,7 @@ reconstruct_pattern_multi <- function(marked_pattern, if (step %% issue == 0) { if(show_graphic == TRUE) { graphics::rect(xwr[1], ywr[1], xwr[2], ywr[2], col="white") - num_species <- from_dummy (p$mark[, species, drop = FALSE]) + num_species <- from_dummy(p$mark[, species, drop = FALSE]) graphics::points(y~x, p, pch = 19, col = 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter]) graphics::text(p$x, p$y, as.integer(num_species), cex = 0.7) diff --git a/R/select_kernel.R b/R/select_kernel.R index d14cea14..052acb29 100644 --- a/R/select_kernel.R +++ b/R/select_kernel.R @@ -43,20 +43,20 @@ select_kernel <- function(kernel_arg, bw, rmax, divisor) { switch(divisor, { rmax_bw <- sqrt(rmax^2 + a/pi) - function(r, d) dunif((r^2-d^2)*pi,-a,+a) + function(r, d) stats::dunif((r^2-d^2)*pi,-a,+a) }, - none = function(r, d) dunif(r,d-a,d+a), - r = function(r, d) dunif(r,d-a,d+a)/(2*pi*r), - d = function(r, d) dunif(r,d-a,d+a)/(2*pi*d) + none = function(r, d) stats::dunif(r,d-a,d+a), + r = function(r, d) stats::dunif(r,d-a,d+a)/(2*pi*r), + d = function(r, d) stats::dunif(r,d-a,d+a)/(2*pi*d) ) }, gaussian = { rmax_bw <- Inf switch(divisor, - function(r, d) dnorm((r^2-d^2)*pi,0,sd=bw), - none = function(r, d) dnorm(r,d,sd = bw), - r = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*r), - d = function(r, d) dnorm(r,d,sd = bw)/ (2*pi*d) + function(r, d) stats::dnorm((r^2-d^2)*pi,0,sd=bw), + none = function(r, d) stats::dnorm(r,d,sd = bw), + r = function(r, d) stats::dnorm(r,d,sd = bw)/ (2*pi*r), + d = function(r, d) stats::dnorm(r,d,sd = bw)/ (2*pi*d) ) }, diff --git a/man/calc_moments.Rd b/man/calc_moments.Rd new file mode 100644 index 00000000..cd4ed616 --- /dev/null +++ b/man/calc_moments.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calc_moments.R +\name{calc_moments} +\alias{calc_moments} +\title{calc_moments} +\usage{ +calc_moments(fn, p, exclude = NULL, x, y, mark, kernel, rmax_bw, r) +} +\arguments{ +\item{fn}{Determination of the weightings of the mark correlation functions.} + +\item{p}{Defines the initial state of the new ponit pattern.} + +\item{exclude}{} + +\item{x}{x Coordinate of the points from the reference point pattern.} + +\item{y}{y Coordinate of the points from the reference point pattern.} + +\item{mark}{Marks the currently viewed point pattern.} + +\item{kernel}{Result of the kernel calculation, calculated with the +calc_kernels function.} + +\item{rmax_bw}{Maximum distance at which the summary statistics are +evaluated + Bandwidth with which the kernels are scaled, so that this is the + standard deviation of the smoothing kernel.} + +\item{r}{is a sequence from rmin to rmax in rcount steps.} +} +\value{ +matrix +} +\description{ +calc_moments (internal) +} +\details{ +Definition of the product-moment function for calculating the contribution +of a point at the coordinates x, y with marking. +} +\keyword{internal} diff --git a/man/compute_statistics.Rd b/man/compute_statistics.Rd new file mode 100644 index 00000000..6416b35b --- /dev/null +++ b/man/compute_statistics.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/compute_statistics.R +\name{compute_statistics} +\alias{compute_statistics} +\title{compute_statistics} +\usage{ +compute_statistics(x, y, k, xr, yr, w_statistics) +} +\arguments{ +\item{x}{x Coordinate of the points from the reference point pattern.} + +\item{y}{y Coordinate of the points from the reference point pattern.} + +\item{k}{Vector of values k; used only if Dk is included in w_statistics +elow.} + +\item{xr}{x extension of the observation window (start, end)} + +\item{yr}{y extension of the observation window (start, end)} + +\item{w_statistics}{} +} +\value{ +list +} +\description{ +compute_statistics (internal) +} +\details{ +Compute optional spatial statistics using the spatstat package. +} +\keyword{internal} diff --git a/man/dummy_transf.Rd b/man/dummy_transf.Rd new file mode 100644 index 00000000..bfbe52c9 --- /dev/null +++ b/man/dummy_transf.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/dummy_transf.R +\name{to_dummy} +\alias{to_dummy} +\alias{dummy_transf} +\title{dummy_transf} +\usage{ +to_dummy(f) +} +\arguments{ +\item{f}{Result of the calc_moments_full function which represents +product-moment contribution of a point at coordinates x, y with marks, +for the whole point pattern.} +} +\value{ +matrix +} +\description{ +dummy_transf (internal) +} +\details{ +Function for the transformation of variables to dummy variables and back +} +\keyword{internal} diff --git a/man/energy_fun.Rd b/man/energy_fun.Rd new file mode 100644 index 00000000..4f38ad92 --- /dev/null +++ b/man/energy_fun.Rd @@ -0,0 +1,69 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/energy_fun.R +\name{Energy_fun} +\alias{Energy_fun} +\alias{energy_fun} +\title{energy_fun} +\usage{ +Energy_fun( + f, + f0, + statistics, + f_, + f0_, + statistics_, + fn, + p, + p_, + Lp, + w_statistics +) +} +\arguments{ +\item{f}{Result of the calc_moments_full function which represents +product-moment contribution of a point at coordinates x, y with marks, +for the whole new ponit pattern.} + +\item{f0}{Column sums of the weights of the brand correlation functions of +the new point pattern.} + +\item{statistics}{Results of the compute_statistics function for the +new point pattern (calculation of optional spatial statistics).} + +\item{f_}{Result of the calc_moments_full function which represents +product-moment contribution of a point at coordinates x, y with marks, +for the whole reference point pattern.} + +\item{f0_}{Column sums of the weights of the brand correlation functions of +the reference point pattern.} + +\item{statistics_}{Results of the compute_statistics function for the +reference point pattern (calculation of optional spatial statistics).} + +\item{fn}{Determination of the weightings of the mark correlation functions.} + +\item{p}{Defines the initial state of the new ponit pattern.} + +\item{p_}{Reference point pattern.} + +\item{Lp}{Distance measure for the calculation of the energy function +(Lp distance, 1 ≤ p < ∞).} + +\item{w_statistics}{ector of named weights for optional spatial statistics +from the ‘spatstat’package to be included in the energy calculation.This may +include Dk, K, Hs, pcf.} +} +\value{ +vector +} +\description{ +energy_fun (internal) +} +\details{ +Defining the Energy_fun function to calculate the "energy" of the pattern +(where a lower energy indicates a better match). + +Defining the Energy_fun function to calculate the "energy" of the pattern +(where a lower energy indicates a better match). +} +\keyword{internal} diff --git a/man/select_kernel.Rd b/man/select_kernel.Rd new file mode 100644 index 00000000..2d02f608 --- /dev/null +++ b/man/select_kernel.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/select_kernel.R +\name{select_kernel} +\alias{select_kernel} +\title{select_kernel} +\usage{ +select_kernel(kernel_arg, bw, rmax, divisor) +} +\arguments{ +\item{kernel_arg}{Parameter of the function +reconstruct_pattern_multi_trait_marks, specifies the kernel to be used to +calculate the energy, possible kernels can be: Gaussian, Epanechnikov, +Rectangular, Cumulative.} + +\item{bw}{Bandwidth with which the kernels are scaled, so that this is the +standard deviation of the smoothing kernel.} + +\item{rmax}{Maximum distance at which the summary statistics are evaluated.} + +\item{divisor}{Divisor in the smoothing kernel, d or r.} +} +\value{ +function +} +\description{ +select_kernel (internal) +} +\details{ +Returns the function of the selected kernel, which is then used to +calculate the kernel. +} +\keyword{internal} From a95db7defce2d3f4602700628614af19e6792fb1 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Thu, 14 Dec 2023 10:16:55 +0000 Subject: [PATCH 10/39] Re-build README.md --- README.md | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index f8370248..4b62f918 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ -README Last updated: 2023-08-31 +README Last updated: 2023-12-14 | CI | Development | CRAN | License | |--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------| @@ -67,7 +67,13 @@ 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. -> +> . + +If you use the `reconstruct_pattern_multi` function, please also cite. + +> Wudel C., Schlicht R., Berger U. (2023). Multi-trait point pattern +> reconstruction of plant ecosystems. Methods in Ecology and Evolution, +> 14, 2668–2679. . The get a BibTex entry, please use `citation("shar")`. From 2d7786f04a4aa649098fd933a0a61000df5ca652 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Thu, 14 Dec 2023 11:45:14 +0100 Subject: [PATCH 11/39] Update gh action test workflow --- .github/workflows/R-CMD-check.yaml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 8c5e168f..415a110b 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -17,7 +17,7 @@ jobs: name: ${{ matrix.config.os }} (${{ matrix.config.r }}) strategy: - fail-fast: false + fail-fast: TRUE matrix: config: - {os: macOS-latest, r: 'release'} @@ -46,4 +46,5 @@ jobs: - uses: r-lib/actions/check-r-package@v2 with: - upload-snapshots: true + args: 'c("--as-cran")' + error-on: '"note"' From 7920cac537ef2f7667692ad606b6774f8cdf9b9a Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Thu, 14 Dec 2023 15:10:21 +0100 Subject: [PATCH 12/39] Some cosmetic changes and fix RCMD --- R/calc_moments.R | 30 +- R/compute_statistics.R | 69 +-- R/dummy_transf.R | 8 +- R/energy_fun.R | 18 +- R/reconstruct_pattern_multi.R | 773 +++++++++++++++---------------- R/select_kernel.R | 7 +- man/calc_moments.Rd | 15 +- man/compute_statistics.Rd | 30 +- man/dummy_transf.Rd | 2 +- man/energy_fun.Rd | 21 +- man/reconstruct_pattern_multi.Rd | 56 +-- man/select_kernel.Rd | 4 +- 12 files changed, 498 insertions(+), 535 deletions(-) diff --git a/R/calc_moments.R b/R/calc_moments.R index 1d5ac81b..f46856c2 100644 --- a/R/calc_moments.R +++ b/R/calc_moments.R @@ -1,19 +1,17 @@ #' calc_moments #' -#' @description calc_moments (internal) +#' @description Calculate moments #' #' @param fn Determination of the weightings of the mark correlation functions. #' @param p Defines the initial state of the new ponit pattern. -#' @param x x Coordinate of the points from the reference point pattern. -#' @param y y Coordinate of the points from the reference point pattern. +#' @param x,y x and y coordinates of the points from the reference point pattern. #' @param mark Marks the currently viewed point pattern. -#' @param kernel Result of the kernel calculation, calculated with the -#' calc_kernels function. +#' @param kernel Result of the kernel calculation, calculated with the calc_kernels function. #' @param rmax_bw Maximum distance at which the summary statistics are #' evaluated + Bandwidth with which the kernels are scaled, so that this is the -#' standard deviation of the smoothing kernel. -#' @param r is a sequence from rmin to rmax in rcount steps. -#' @param exclude +#' standard deviation of the smoothing kernel. +#' @param r Sequence from rmin to rmax in rcount steps. +#' @param exclude Vector indicating which values not to use. #' #' @details #' Definition of the product-moment function for calculating the contribution @@ -27,33 +25,33 @@ #' @keywords internal calc_moments <- function(fn, p, - exclude=NULL, + exclude = NULL, x, y, mark, kernel, rmax_bw, - r) - { + r) { + d2 <- (p$x-x)^2 + (p$y-y)^2 use <- d2 <= rmax_bw^2 use[exclude] <- FALSE z <- crossprod(p$mark[use, , drop = FALSE], outer(sqrt(d2[use]), r, function(d, r) kernel(r, d))) z[fn$i, , drop = FALSE] * mark[fn$j] + z[fn$j, , drop = FALSE] * mark[fn$i] - } +} calc_moments_full <- function(fn, p, kernel, rmax_bw, - r) - { + r) { + f <- 0 for (i in seq_len(nrow(p))) { f <- f + calc_moments(fn, p, i:nrow(p), p$x[i], p$y[i], p$mark[i, ], kernel, rmax_bw, r) - } + } rownames(f) <- paste(names(fn$i), names(fn$j), sep = ":") f - } +} diff --git a/R/compute_statistics.R b/R/compute_statistics.R index 7657e17f..e3bcf20a 100644 --- a/R/compute_statistics.R +++ b/R/compute_statistics.R @@ -1,14 +1,14 @@ #' compute_statistics #' -#' @description compute_statistics (internal) +#' @description Compute summary statistics #' -#' @param x x Coordinate of the points from the reference point pattern. -#' @param y y Coordinate of the points from the reference point pattern. -#' @param k Vector of values k; used only if Dk is included in w_statistics -#' elow. -#' @param xr x extension of the observation window (start, end) -#' @param yr y extension of the observation window (start, end) -#' @param w_statistics +#' @param x,y x and y coordinates of the points from the reference point pattern. +#' @param k Vector of values k; used only if Dk is included in w_statistics below. +#' @param xr,yr x and y extension of the observation window (start, end). +#' @param w_statistics vector of named weights for optional spatial statistics +#' from the \code{spatstat} package to be included in the energy calculation. This may +#' include Dk, K, Hs, pcf. +#' @param bw,divisor,kernel_arg,rmax Several parameters related to summary function. #' #' @details #' Compute optional spatial statistics using the spatstat package. @@ -19,34 +19,35 @@ #' @rdname compute_statistics #' #' @keywords internal -compute_statistics <- function(x, y, k, xr, yr, w_statistics) { +compute_statistics <- function(x, y, k, xr, yr, w_statistics, bw, divisor, kernel_arg, rmax) { + stat <- names(w_statistics) names(stat) <- stat lapply(stat, function(name) switch(name, - # Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. - - Dk = { - nnd_ <- as.matrix(spatstat.geom::nndist(x, y, k=k)) - apply(nnd_, 2, function(z) cumsum(graphics::hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) - }, - # Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. - - K = { - kest<-spatstat.explore::Kest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), rmax=rmax, correction="none") - kest$un - }, - # Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. + # Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. + Dk = { + nnd_ <- as.matrix(spatstat.geom::nndist(x, y, k=k)) + apply(nnd_, 2, function(z) cumsum(graphics::hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) + }, - pcf = { - pcfest<-spatstat.explore::pcf(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") - pcfest$un - }, - # Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. + # Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. + K = { + kest<-spatstat.explore::Kest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), rmax=rmax, correction="none") + kest$un + }, - Hs = { - hest<-spatstat.explore::Hest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), correction="none") - hest$raw - }, - stop("unknown statistic") - )) - } + # Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. + pcf = { + pcfest<-spatstat.explore::pcf(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") + pcfest$un + }, + # Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. + Hs = { + hest<-spatstat.explore::Hest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), correction="none") + hest$raw + }, + # wrong selection + stop("unknown statistic") + ) + ) +} diff --git a/R/dummy_transf.R b/R/dummy_transf.R index ce823236..9a3db7e9 100644 --- a/R/dummy_transf.R +++ b/R/dummy_transf.R @@ -1,6 +1,6 @@ #' dummy_transf #' -#' @description dummy_transf (internal) +#' @description Tranfsorm to dummy variables #' #' @param f Result of the calc_moments_full function which represents #' product-moment contribution of a point at coordinates x, y with marks, @@ -15,17 +15,15 @@ #' @rdname dummy_transf #' #' @keywords internal -#' to_dummy <- function(f) { x <- matrix(0, length(f), nlevels(f), dimnames=list(names(f), levels(f))) x[cbind(seq_along(f), as.integer(f))] <- 1 x - } +} from_dummy <- function(x, levels=colnames(x)) { f <- as.integer(x %*% seq_along(levels)) levels(f) <- levels class(f) <- "factor" f - } - +} diff --git a/R/energy_fun.R b/R/energy_fun.R index 19b23716..bff33994 100644 --- a/R/energy_fun.R +++ b/R/energy_fun.R @@ -1,17 +1,7 @@ #' energy_fun #' -#' @description energy_fun (internal) +#' @description Energy function #' -#' @param f_ Result of the calc_moments_full function which represents -#' product-moment contribution of a point at coordinates x, y with marks, -#' for the whole reference point pattern. -#' @param f0_ Column sums of the weights of the brand correlation functions of -#' the reference point pattern. -#' @param statistics_ Results of the compute_statistics function for the -#' reference point pattern (calculation of optional spatial statistics). -#' @details -#' Defining the Energy_fun function to calculate the "energy" of the pattern -#' (where a lower energy indicates a better match). #' @param f Result of the calc_moments_full function which represents #' product-moment contribution of a point at coordinates x, y with marks, #' for the whole new ponit pattern. @@ -23,9 +13,9 @@ #' @param p Defines the initial state of the new ponit pattern. #' @param p_ Reference point pattern. #' @param Lp Distance measure for the calculation of the energy function -#' (Lp distance, 1 ≤ p < ∞). -#' @param w_statistics ector of named weights for optional spatial statistics -#' from the ‘spatstat’package to be included in the energy calculation.This may +#' (Lp distance, 1 <= p 1) { names_reconstruction <- paste0("reconstruction_", seq_len(n_repetitions)) reconstruction_list <- vector("list", length = n_repetitions) names(reconstruction_list) <- names_reconstruction } -# Loop to perform multiple reconstructions. + + # Loop to perform multiple reconstructions. for (t in seq_len(n_repetitions)) { - # Load the reference point pattern as a data frame with the components x, y, mark, where x, y are the coordinates of the point and mark is a matrix representing the marks or their dummy values. - p_ <- data.frame( - x = marked_pattern$x, - y = marked_pattern$y) + # Load the reference point pattern as a data frame with the components x, y, + # mark, where x, y are the coordinates of the point and mark is a matrix representing the marks or their dummy values. + p_ <- data.frame(x = marked_pattern$x, y = marked_pattern$y) - p_$mark <- cbind(`1`= 1, - diameter = marked_pattern$marks[[1]], - to_dummy(marked_pattern$marks[[2]])) + p_$mark <- cbind(`1`= 1, diameter = marked_pattern$marks[[1]], + to_dummy(marked_pattern$marks[[2]])) - marknames <- colnames(p_$mark) - diameter <- 2 - species <- seq_along(marknames)[-(1:2)] + marknames <- colnames(p_$mark) + diameter <- 2 + species <- seq_along(marknames)[-(1:2)] -# Check whether certain requirements are met; if not, the reconstruction is aborted and an error message is displayed. - if (is.null(marked_pattern[["marks"]])) { - stop("'marked_pattern' must be marked", call. = FALSE) - } + # Check whether certain requirements are met; if not, the reconstruction is + # aborted and an error message is displayed. + if (is.null(marked_pattern[["marks"]])) { + stop("'marked_pattern' must be marked", call. = FALSE) + } - if (!inherits(marked_pattern[["marks"]][[1]], "numeric")) { - stop("mark one must be 'numeric', an example would be the DBH - (Diameter at breast height).", call. = FALSE) - } + if (!inherits(marked_pattern[["marks"]][[1]], "numeric")) { + stop("mark one must be 'numeric', an example would be the DBH (Diameter at breast height).", + call. = FALSE) + } - if (!inherits(marked_pattern[["marks"]][[2]], "factor")) { - stop("mark two must be a 'factor', an example would be the tree species.", - call. = FALSE) - } + if (!inherits(marked_pattern[["marks"]][[2]], "factor")) { + stop("mark two must be a 'factor', an example would be the tree species.", + call. = FALSE) + } - if (n_repetitions < 1) { - stop("n_repetitions must be at least 1 for the function to be executed.", - call. = FALSE) - } + if (n_repetitions < 1) { + stop("n_repetitions must be at least 1 for the function to be executed.", + call. = FALSE) + } -# Definition of parameters for the estimation of correlation functions. - rmin <- rmax / rcount - r <- seq(rmin, rmax, , rcount) # MH: Is rcount by or length.out argument? - -# Calculation of the kernels. - sel_kernel <- select_kernel(kernel_arg, bw, rmax, divisor) - kernel <- as.function(sel_kernel[[1]]) - rmax_bw <- as.numeric(sel_kernel[[2]]) - -# Determination of the weightings of the mark correlation functions. - fn <- list() - for (i in seq_along(marknames)) for (j in seq_along(marknames)) if (i <= j) { - fn$i <- c(fn$i,i) - fn$j <- c(fn$j,j) - fn$w <- c(fn$w, - if (i == diameter && j == diameter) w_markcorr["d_d"] - else if(i == 1 || j ==1) w_markcorr["all"] - else if(i == diameter || j == diameter) w_markcorr["d_all"] - else w_markcorr["all_all"]) - fn$w0<-c(fn$w0, - if (i == diameter && j == diameter) w_markcorr["d_d0"] - else if(i == 1 || j == 1) w_markcorr["all0"] - else if(i == diameter || j == diameter) w_markcorr["d_all0"] - else w_markcorr["all_all0"]) - } - names(fn$i) <- marknames[fn$i] - names(fn$j) <- marknames[fn$j] - - # Defines the initial state of the new ponit pattern. - n <- nrow(p_) - xwr <- obs_window$xrange - ywr <- obs_window$yrange - p <- p_[sample.int(nrow(p_),ceiling(nrow(p_)*((diff(xwr)*diff(ywr))/(diff(xr)*diff(yr)))),TRUE), ]# rpois - p$x <- stats::runif(nrow(p), xwr[1], xwr[2]) - p$y <- stats::runif(nrow(p), ywr[1], ywr[2]) - p$mark[, diameter] <- stats::quantile(p_$mark[, diameter], - probs = stats::runif(nrow(p), 0, 1), type = 4) - p$mark[, species] <- p_$mark[, species, drop = FALSE][ + # Definition of parameters for the estimation of correlation functions. + rmin <- rmax / rcount + r <- seq(rmin, rmax, length.out = rcount) + + # Calculation of the kernels. + sel_kernel <- select_kernel(kernel_arg, bw, rmax, divisor) + kernel <- as.function(sel_kernel[[1]]) + rmax_bw <- as.numeric(sel_kernel[[2]]) + + # get observation window of data + obs_window <- marked_pattern$window + + # Determination of the weightings of the mark correlation functions. + fn <- list() + for (i in seq_along(marknames)) for (j in seq_along(marknames)) if (i <= j) { + fn$i <- c(fn$i,i) + fn$j <- c(fn$j,j) + fn$w <- c(fn$w, + if (i == diameter && j == diameter) w_markcorr["d_d"] + else if(i == 1 || j ==1) w_markcorr["all"] + else if(i == diameter || j == diameter) w_markcorr["d_all"] + else w_markcorr["all_all"]) + fn$w0<-c(fn$w0, + if (i == diameter && j == diameter) w_markcorr["d_d0"] + else if(i == 1 || j == 1) w_markcorr["all0"] + else if(i == diameter || j == diameter) w_markcorr["d_all0"] + else w_markcorr["all_all0"]) + } + names(fn$i) <- marknames[fn$i] + names(fn$j) <- marknames[fn$j] + + # Defines the initial state of the new ponit pattern. + xwr <- obs_window$xrange + ywr <- obs_window$yrange + p <- p_[sample.int(nrow(p_),ceiling(nrow(p_)*((diff(xwr)*diff(ywr))/(diff(xr)*diff(yr)))),TRUE), ]# rpois + p$x <- stats::runif(nrow(p), xwr[1], xwr[2]) + p$y <- stats::runif(nrow(p), ywr[1], ywr[2]) + p$mark[, diameter] <- stats::quantile(p_$mark[, diameter], + probs = stats::runif(nrow(p), 0, 1), type = 4) + p$mark[, species] <- p_$mark[, species, drop = FALSE][ sample.int(nrow(p_), nrow(p), replace = TRUE),, drop = FALSE] -# Calculates the functions for the reference and the new dot pattern as well as calculating the "energy" that measures their distance. - f_ <- calc_moments_full(fn, p_, kernel, rmax_bw, r) - f0_ <- colSums(p_$mark[, fn$i] * p_$mark[, fn$j]) - names(f0_) <- rownames(f_) - statistics_<- compute_statistics(p_$x, p_$y, k, xr, yr, w_statistics) - f <- calc_moments_full(fn, p, kernel, rmax_bw, r) - f0 <- colSums(p$mark[, fn$i] * p$mark[, fn$j]) - names(f0) <- rownames(f) - statistics <- compute_statistics(p$x, p$y, k, xwr, ywr, w_statistics) - -# Prepare the graphical output. - if(show_graphic == TRUE) { - graphics::par(mfrow = 1:2) - num_species <- from_dummy(p_$mark[, species, drop = FALSE]) - plot(y~x, p_, pch=19, col= 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter], xlim = xr, - ylim = yr, xaxs ="i", yaxs ="i", main ="Reference", xlab ="x [m]", - ylab ="y [m]") - graphics::text(p_$x, p_$y, as.integer(num_species), cex=0.7) - - graphics::plot(y~x, p, type = "n", - xlim = xwr, ylim = ywr, xaxs = "i", yaxs = "i", main = "Reconstructed", - xlab = "x [m]", ylab = "y [m]") - graphics::clip(xwr[1], xwr[2], ywr[1], ywr[2]) - } - -# Show warning if certain distances between pairs of trees are not present. - energy <- Energy_fun(f, f0, statistics, f_, f0_, statistics_, fn, p, p_, Lp, w_statistics)["energy"] - -# Prepares variables for the storage of progress. - energy_launch <- as.vector(energy) - energy_course <- data.frame(i = seq(from = 1, to = max_steps,by = 1), - energy = NA) - energy_improvement <- 0L - number_of_actions <- integer(0) - number_of_actions_with_energy_improvement <- integer(0) - no_changes_energy <- Inf - no_changes_counter <- 0L - step_list <- integer(0) - action_list <- character(0) - Energy_list <- numeric(0) - -# loop to improve the newly created dot pattern by reducing its energy. - step <- 0L - system.time(repeat { - energy_course[step, 2] <- energy - -# Updating the graphical output of all "issue" steps. - if (step %% issue == 0) { - if(show_graphic == TRUE) { - graphics::rect(xwr[1], ywr[1], xwr[2], ywr[2], col="white") - num_species <- from_dummy(p$mark[, species, drop = FALSE]) - - graphics::points(y~x, p, pch = 19, col = 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter]) - graphics::text(p$x, p$y, as.integer(num_species), cex = 0.7) - } - -# Generates text output with the current calculated values (for example the energy), this is updated every "issue" simulation step. - message("> Progress:", if(n_repetitions > 1) names_reconstruction[[t]], " || iterations: ", step, - " || Simulation progress: ", floor(step/max_steps * 100), "%", - " || energy = ", round(energy, 5), " || energy improvement = ", - energy_improvement, "\t\t\r", appendLF = FALSE) - - Sys.sleep(0) - utils::flush.console() + # Calculates the functions for the reference and the new dot pattern as well as calculating the "energy" that measures their distance. + f_ <- calc_moments_full(fn, p_, kernel, rmax_bw, r) + f0_ <- colSums(p_$mark[, fn$i] * p_$mark[, fn$j]) + names(f0_) <- rownames(f_) + statistics_<- compute_statistics(p_$x, p_$y, k, xr, yr, w_statistics, bw, divisor, kernel_arg, rmax) + f <- calc_moments_full(fn, p, kernel, rmax_bw, r) + f0 <- colSums(p$mark[, fn$i] * p$mark[, fn$j]) + names(f0) <- rownames(f) + statistics <- compute_statistics(p$x, p$y, k, xwr, ywr, w_statistics) + + # Prepare the graphical output. + if(plot == TRUE) { + graphics::par(mfrow = 1:2) + num_species <- from_dummy(p_$mark[, species, drop = FALSE]) + graphics::plot(y~x, p_, pch=19, col= 2L + as.integer(num_species), + cex = 1.3 + 4 * mark[, diameter], xlim = xr, ylim = yr, + xaxs ="i", yaxs ="i", main ="Reference", xlab ="x [m]", ylab ="y [m]") + graphics::text(p_$x, p_$y, as.integer(num_species), cex=0.7) + + graphics::plot(y~x, p, type = "n", xlim = xwr, ylim = ywr, xaxs = "i", + yaxs = "i", main = "Reconstructed", xlab = "x [m]", ylab = "y [m]") + graphics::clip(xwr[1], xwr[2], ywr[1], ywr[2]) } -# the next code block aborts the reconstruction if the energy decreases by less than tol in "no_changes" intervals of steps_tol simulation steps. - if (step %% steps_tol == 0) { - if(energy > no_changes_energy * (1-tol)) { - no_changes_counter <- no_changes_counter + 1L - if(no_changes_counter == no_changes) { - message("the simulation was terminated, because the energy did not - decrease in ", no_changes * issue, " simulation steps.") - stop_criterion <- "no_changes" - break + + # Show warning if certain distances between pairs of trees are not present. + energy <- Energy_fun(f, f0, statistics, f_, f0_, statistics_, fn, p, p_, Lp, w_statistics)["energy"] + + # Prepares variables for the storage of progress. + energy_launch <- as.vector(energy) + energy_course <- data.frame(i = seq(from = 1, to = max_steps,by = 1), + energy = NA) + energy_improvement <- 0L + number_of_actions <- integer(0) + no_change_energy <- Inf + no_change_counter <- 0L + step_list <- integer(0) + action_list <- character(0) + Energy_list <- numeric(0) + number_of_actions_with_energy_improvement <- integer(0) + + # loop to improve the newly created dot pattern by reducing its energy. + step <- 0L + system.time(repeat { + + energy_course[step, 2] <- energy + + # Updating the graphical output of all "issue" steps. + if (step %% issue == 0) { + if(plot == TRUE) { + graphics::rect(xwr[1], ywr[1], xwr[2], ywr[2], col="white") + num_species <- from_dummy(p$mark[, species, drop = FALSE]) + + graphics::points(y~x, p, pch = 19, col = 2L + as.integer(num_species), cex = 1.3 + 4 * mark[, diameter]) + graphics::text(p$x, p$y, as.integer(num_species), cex = 0.7) } - } else { - no_changes_counter<-0L - stop_criterion<-"max_steps" - } - no_changes_energy <- energy - } - if (step < max_steps) step <- step + 1 else break - action <- sample(c("move_coordinate", "switch_coords", "pick_mark_one", - "pick_mark_two", "exchange_mark_one", "exchange_mark_two"), - 1,, prob_of_actions) - number_of_actions[action] <- (if (is.na(number_of_actions[action])) - 0L else number_of_actions[action]) + 1L - - statistics.new <- statistics - -# Switch selection for the possible changes to the reconstructed point pattern for energy minimisation (probabilities of how often an action is taken: 40%, 10%, 20%, 10%, 10%). - switch(action, - # Displacement of coordinates of a point in the new point pattern, is applied in xx% of the cases. - move_coordinate = { - i <- sample.int(nrow(p), 1, replace = TRUE) - s <- if(sd=="step") nrow(p) * 1 / step else sd - x <- stats::rnorm(1, p$x[i], diff(xwr) * s) %% xwr[2] - y <- stats::rnorm(1, p$y[i], diff(ywr) * s) %% ywr[2] - mdiff <- p$mark[i, ] - f.new <- f - calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, - kernel, rmax_bw, r) + - calc_moments(fn, p, i, x, y, mdiff, kernel, rmax_bw, r) - f0.new<- f0 - statistics.new <- compute_statistics(replace(p$x, i, x), replace(p$y, i, y), k, xwr, ywr, w_statistics) - }, + # Generates text output with the current calculated values (for example the energy), this is updated every "issue" simulation step. + message("> Progress:", if(n_repetitions > 1) names_reconstruction[[t]], " || iterations: ", step, + " || Simulation progress: ", floor(step/max_steps * 100), "%", + " || energy = ", round(energy, 5), " || energy improvement = ", + energy_improvement, "\t\t\r", appendLF = FALSE) - # Swaps the coordinates of two randomly drawn points from the new point pattern, applied in xx% of the trap. - switch_coords = { - i <- sample.int(nrow(p), 2, replace = FALSE) - mdiff <- p$mark[i[1], ] - p$mark[i[2], ] - f.new <- f - calc_moments(fn, p, i, p$x[i[1]], p$y[i[1]], - mdiff, - kernel, rmax_bw, r) + - calc_moments(fn, p, i, p$x[i[2]], - p$y[i[2]], mdiff, kernel, rmax_bw, r) - f0.new<- f0 - }, + Sys.sleep(0) + utils::flush.console() + } - # Displacement of coordinates of a point in the new point pattern, applied in xx% of the cases. - exchange_mark_one = { - i <- sample.int(nrow(p), 2, replace = FALSE) - m <- p$mark[i, ] - m[, diameter] <- m[2:1, diameter] - mdiff <- m[1, ] - p$mark[i[1], ] - q <- p[i, ] - q$mark[1, ] <- m[1, ] - f.new <- f + calc_moments(fn, p, i[1], p$x[i[1]], - p$y[i[1]], mdiff, kernel, - rmax_bw, r) - - calc_moments(fn, p, i, p$x[i[2]], p$y[i[2]], - mdiff, kernel, rmax_bw, r) - - calc_moments(fn, q, 2, q$x[2], q$y[2], mdiff, - kernel, rmax_bw, r) - f0.new <- f0 + m[1,fn$i] * m[1, fn$j] - - p$mark[i[1], fn$i] * p$mark[i[1], fn$j] + m[2, fn$i] * - m[2, fn$j] - p$mark[i[2], fn$i] * - p$mark[i[2], fn$j] - }, + # the next code block aborts the reconstruction if the energy decreases by less than tol in "no_change" intervals of steps_tol simulation steps. + if (step %% steps_tol == 0) { + if(energy > no_change_energy * (1-tol)) { + no_change_counter <- no_change_counter + 1L + if(no_change_counter == no_change) { + message("the simulation was terminated, because the energy did not + decrease in ", no_change * issue, " simulation steps.") + stop_criterion <- "no_change" + break + } + } else { + no_change_counter<-0L + stop_criterion<-"max_steps" + } + no_change_energy <- energy + } - # Swaps the type assignment of two randomly drawn points from the new point pattern, applied in 10% of the trap. - exchange_mark_two = { - i <- sample.int(nrow(p), 2, replace = FALSE) - m <- p$mark[i, ] - m[, species] <- m[2:1, species] - mdiff <- m[1, ] - p$mark[i[1], ] - q <- p[i, ] - q$mark <- m - - f.new <- f + calc_moments(fn, p, i[1], p$x[i[1]], p$y[i[1]], - mdiff, kernel, rmax_bw, r) - - calc_moments(fn, p, i, p$x[i[2]], p$y[i[2]], - mdiff, kernel, rmax_bw, r) - - calc_moments(fn, q, 2, q$x[2], q$y[2], mdiff, - kernel, rmax_bw, r) - f0.new <- f0 + m[1, fn$i] * m[1, fn$j] - p$mark[i[1], fn$i] * - p$mark[i[1], fn$j] + m[2, fn$i] * m[2, fn$j] - - p$mark[i[2], fn$i] * p$mark[i[2], fn$j] - }, - - # If the distribution (continuous function) of the diameter of the reference pattern generates a randomly drawn value for a randomly selected point in the new point pattern, the trap is applied in 20%. - pick_mark_one = { - i <- sample.int(nrow(p), 1, replace = TRUE) - m <- p$mark[i, ] - m[diameter] <- stats::quantile(p_$mark[,diameter],probs = stats::runif(1,0,1), - type = 4) - mdiff <- m - p$mark[i, ] - f.new <- f + calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, - kernel, rmax_bw, r) - f0.new<- f0 + m[fn$i] * m[fn$j] - p$mark[i, fn$i] * - p$mark[i, fn$j] + if (step < max_steps) step <- step + 1 else break + action <- sample(c("move_coordinate", "switch_coords", "pick_mark_one", + "pick_mark_two", "exchange_mark_one", "exchange_mark_two"), + 1, prob = prob_of_actions) + number_of_actions[action] <- (if (is.na(number_of_actions[action])) 0L else number_of_actions[action]) + 1L + + statistics.new <- statistics + + # Switch selection for the possible changes to the reconstructed point pattern for energy minimisation (probabilities of how often an action is taken: 40%, 10%, 20%, 10%, 10%). + switch(action, + # Displacement of coordinates of a point in the new point pattern, is applied in xx% of the cases. + move_coordinate = { + i <- sample.int(nrow(p), 1, replace = TRUE) + s <- if(sd=="step") nrow(p) * 1 / step else sd + x <- stats::rnorm(1, p$x[i], diff(xwr) * s) %% xwr[2] + y <- stats::rnorm(1, p$y[i], diff(ywr) * s) %% ywr[2] + mdiff <- p$mark[i, ] + f.new <- f - calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, + kernel, rmax_bw, r) + + calc_moments(fn, p, i, x, y, mdiff, kernel, rmax_bw, r) + f0.new <- f0 + statistics.new <- compute_statistics(replace(p$x, i, x), replace(p$y, i, y), k, xwr, ywr, w_statistics) + }, + # Swaps the coordinates of two randomly drawn points from the new point pattern, applied in xx% of the trap. + switch_coords = { + i <- sample.int(nrow(p), 2, replace = FALSE) + mdiff <- p$mark[i[1], ] - p$mark[i[2], ] + f.new <- f - calc_moments(fn, p, i, p$x[i[1]], p$y[i[1]], + mdiff, + kernel, rmax_bw, r) + + calc_moments(fn, p, i, p$x[i[2]], + p$y[i[2]], mdiff, kernel, rmax_bw, r) + f0.new<- f0 + }, + # Displacement of coordinates of a point in the new point pattern, applied in xx% of the cases. + exchange_mark_one = { + i <- sample.int(nrow(p), 2, replace = FALSE) + m <- p$mark[i, ] + m[, diameter] <- m[2:1, diameter] + mdiff <- m[1, ] - p$mark[i[1], ] + q <- p[i, ] + q$mark[1, ] <- m[1, ] + f.new <- f + calc_moments(fn, p, i[1], p$x[i[1]], + p$y[i[1]], mdiff, kernel, + rmax_bw, r) - + calc_moments(fn, p, i, p$x[i[2]], p$y[i[2]], + mdiff, kernel, rmax_bw, r) - + calc_moments(fn, q, 2, q$x[2], q$y[2], mdiff, + kernel, rmax_bw, r) + f0.new <- f0 + m[1,fn$i] * m[1, fn$j] - + p$mark[i[1], fn$i] * p$mark[i[1], fn$j] + m[2, fn$i] * + m[2, fn$j] - p$mark[i[2], fn$i] * + p$mark[i[2], fn$j] + }, + # Swaps the type assignment of two randomly drawn points from the new point + # pattern, applied in 10% of the trap. + exchange_mark_two = { + i <- sample.int(nrow(p), 2, replace = FALSE) + m <- p$mark[i, ] + m[, species] <- m[2:1, species] + mdiff <- m[1, ] - p$mark[i[1], ] + q <- p[i, ] + q$mark <- m + + f.new <- f + calc_moments(fn, p, i[1], p$x[i[1]], p$y[i[1]], + mdiff, kernel, rmax_bw, r) - + calc_moments(fn, p, i, p$x[i[2]], p$y[i[2]], + mdiff, kernel, rmax_bw, r) - + calc_moments(fn, q, 2, q$x[2], q$y[2], mdiff, + kernel, rmax_bw, r) + f0.new <- f0 + m[1, fn$i] * m[1, fn$j] - p$mark[i[1], fn$i] * + p$mark[i[1], fn$j] + m[2, fn$i] * m[2, fn$j] - + p$mark[i[2], fn$i] * p$mark[i[2], fn$j] + }, + # If the distribution (continuous function) of the diameter of the reference pattern generates a randomly drawn value for a randomly selected point in the new point pattern, the trap is applied in 20%. + pick_mark_one = { + i <- sample.int(nrow(p), 1, replace = TRUE) + m <- p$mark[i, ] + m[diameter] <- stats::quantile(p_$mark[,diameter],probs = stats::runif(1,0,1), + type = 4) + mdiff <- m - p$mark[i, ] + f.new <- f + calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, + kernel, rmax_bw, r) + f0.new<- f0 + m[fn$i] * m[fn$j] - p$mark[i, fn$i] * + p$mark[i, fn$j] + }, + # Draws a random value for a point from the new point pattern from the type distribution (discrete function) of the reference pattern, is applied in 10% of the traps. + pick_mark_two = { + i <- sample.int(nrow(p), 1, replace = TRUE) + j <- sample.int(nrow(p_), 1, replace = TRUE) + m <- p$mark[i, ] + m[species] <- p_$mark[j, species] + mdiff <- m - p$mark[i, ] + f.new <- f + calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, + kernel, rmax_bw, r) + f0.new <- f0 + m[fn$i] * m[fn$j] - p$mark[i, fn$i] * + p$mark[i, fn$j] + }, + # return error + stop("undefined case") + ) + + # calculates the energy + Energy <- Energy_fun(f.new, f0.new, statistics.new, f_, f0_, statistics_, fn, p, p_, Lp, w_statistics) + + energy.new<-Energy[["energy"]] + + # Sets the currently calculated energy as the new reference value if it is less than the previous energy. + if(energy.new >= energy) next + f <- f.new + f0 <- f0.new + statistics <- statistics.new + energy <- energy.new + + # Executes the previously defined actions. + switch(action, + # move + move_coordinate = { + p$x[i] <- x + p$y[i] <- y + }, + # switch + switch_coords = { + p$x[i] <- p$x[rev(i)] + p$y[i] <- p$y[rev(i)] + }, + pick_mark_one =, + pick_mark =, + pick_mark_two =, + exchange_mark_one =, + exchange_mark_two = { + p$mark[i, ] <- m + }, + # no action defined + stop("undefined case") + ) + + # Saves the intermediate results and increases running numbers. + if (energy_evaluation == TRUE) { + step_list <-c(step_list,step) + action_list <-c(action_list, action) + Energy_list <-rbind(Energy_list, Energy) + number_of_actions_with_energy_improvement[action] <- + (if (is.na(number_of_actions_with_energy_improvement[action])) + 0L else number_of_actions_with_energy_improvement[action]) + 1L + } + energy_improvement <- energy_improvement + 1L + + # End of reconstruction loop. + }) -> process.time + + message("\n") + + # Saves all results Transfers them to the "reconstruction" list. + p_$species <- from_dummy(p_$mark[, species, drop = FALSE]) + p$species <- from_dummy(p$mark[, species, drop = FALSE]) + method <- "Reconstruction of a homogeneous point pattern" + Parameter_setting <- list(n_repetitions=n_repetitions, max_steps=max_steps, + no_change=no_change, rcount=rcount, rmax=rmax, + issue=issue, divisor=divisor, kernel_arg=kernel_arg, + timing=timing, energy_evaluation=energy_evaluation, + plot=plot, Lp=Lp, k=k, bw=bw, sd=sd, + prob_of_actions=prob_of_actions, + w_markcorr=w_markcorr, w_statistics=w_statistics) + iterations <- step + energy_current <- energy_course[step, 2] + adapted_p_<- data.frame(p_$x, p_$y, p_$mark[,2], p_$species) + colnames(adapted_p_)<-c("x", "y", "diameter", "species") + adapted_p<- data.frame(p$x, p$y, p$mark[,2], p$species) + colnames(adapted_p)<-c("x", "y", "diameter", "species") + win_change<- if (xr[1] != xwr[1] | xr[2] != xwr[2] | yr[1] != ywr[1] | yr[2] != ywr[2]) {TRUE}else{FALSE} + + # Saves the results. + reconstruction <- + list(reference = adapted_p_, + reconstructed = adapted_p, + window = c(xr, yr), + obs_window = + if (win_change == TRUE) { + c(xwr, ywr) }, - - # Draws a random value for a point from the new point pattern from the type distribution (discrete function) of the reference pattern, is applied in 10% of the traps. - pick_mark_two = { - i <- sample.int(nrow(p), 1, replace = TRUE) - j <- sample.int(nrow(p_), 1, replace = TRUE) - m <- p$mark[i, ] - m[species] <- p_$mark[j, species] - mdiff <- m - p$mark[i, ] - f.new <- f + calc_moments(fn, p, i, p$x[i], p$y[i], mdiff, - kernel, rmax_bw, r) - f0.new <- f0 + m[fn$i] * m[fn$j] - p$mark[i, fn$i] * - p$mark[i, fn$j] + r = r, + f_reference = f_, + f_reconstructed = f, + Parameter_setting = Parameter_setting, + method = method, + stop_criterion = stop_criterion, + iterations = step, + simulation_time = + if (timing == TRUE) { + paste(round(process.time[3], 2), "s") }, - stop("undefined case") - ) -# calculates the energy - Energy <- Energy_fun(f.new, f0.new, statistics.new, f_, f0_, statistics_, fn, p, p_, Lp, w_statistics) - - energy.new<-Energy[["energy"]] - -# Sets the currently calculated energy as the new reference value if it is less than the previous energy. - if(energy.new >= energy) next - f <- f.new - f0 <- f0.new - statistics <- statistics.new - energy <- energy.new - -# Executes the previously defined actions. - switch(action, - move_coordinate = { - p$x[i] <- x - p$y[i] <- y + energy_launch = energy_launch, + energy_course = energy_course, + energy_current = energy_current, + energy_improvement = energy_improvement, + number_of_actions = + if(energy_evaluation == TRUE) { + number_of_actions }, - switch_coords = { - p$x[i] <- p$x[rev(i)] - p$y[i] <- p$y[rev(i)] + number_of_actions_with_energy_improvement = + if(energy_evaluation == TRUE) { + number_of_actions_with_energy_improvement }, - pick_mark_one =, - pick_mark =, - pick_mark_two =, - exchange_mark_one =, - exchange_mark_two = - { - p$mark[i, ] <- m - }, - stop("undefined case") - ) - -# Saves the intermediate results and increases running numbers. - if (energy_evaluation == TRUE) { - step_list <-c(step_list,step) - action_list <-c(action_list, action) - Energy_list <-rbind(Energy_list, Energy) - number_of_actions_with_energy_improvement[action] <- - (if (is.na(number_of_actions_with_energy_improvement[action])) - 0L else number_of_actions_with_energy_improvement[action]) + 1L - } - energy_improvement <- energy_improvement + 1L - -# End of reconstruction loop. - -}) -> process.time - message("\n") - - # Saves all results Transfers them to the "reconstruction" list. - p_$species <- from_dummy(p_$mark[, species, drop = FALSE]) - p$species <- from_dummy(p$mark[, species, drop = FALSE]) - method <- "Reconstruction of a homogeneous point pattern" - Parameter_setting <- list(n_repetitions=n_repetitions, max_steps=max_steps, - no_changes=no_changes, rcount=rcount, rmax=rmax, - issue=issue, divisor=divisor, kernel_arg=kernel_arg, - timing=timing, energy_evaluation=energy_evaluation, - show_graphic=show_graphic, Lp=Lp, k=k, bw=bw, sd=sd, - prob_of_actions=prob_of_actions, - w_markcorr=w_markcorr, w_statistics=w_statistics) - iterations <- step - energy_current <- energy_course[step, 2] - adapted_p_<- data.frame(p_$x, p_$y, p_$mark[,2], p_$species) - colnames(adapted_p_)<-c("x", "y", "diameter", "species") - adapted_p<- data.frame(p$x, p$y, p$mark[,2], p$species) - colnames(adapted_p)<-c("x", "y", "diameter", "species") - win_change<- if (xr[1] != xwr[1] | xr[2] != xwr[2] | yr[1] != ywr[1] | yr[2] != ywr[2]) {TRUE}else{FALSE} - -# Saves the results. - reconstruction <- - list( reference = adapted_p_, - reconstructed = adapted_p, - window = c(xr, yr), - obs_window = - if (win_change == TRUE) { - c(xwr, ywr) - }, - r = r, - f_reference = f_, - f_reconstructed = f, - Parameter_setting = Parameter_setting, - method = method, - stop_criterion = stop_criterion, - iterations = step, - simulation_time = - if (timing == TRUE) { - paste(round(process.time[3], 2), "s") - }, - energy_launch = energy_launch, - energy_course = energy_course, - energy_current = energy_current, - energy_improvement = energy_improvement, - number_of_actions = - if(energy_evaluation == TRUE) { - number_of_actions - }, - number_of_actions_with_energy_improvement = - if(energy_evaluation == TRUE) { - number_of_actions_with_energy_improvement - }, - energy_details = - if(energy_evaluation == TRUE) { + energy_details = + if(energy_evaluation == TRUE) { data.frame(step_list, action_list, Energy_list) - }) + }) - if (!(timing && energy_evaluation && win_change)) { - reconstruction <- reconstruction[-which(sapply(reconstruction, is.null))] - } + if (!(timing && energy_evaluation && win_change)) { + reconstruction <- reconstruction[-which(sapply(reconstruction, is.null))] + } -# Adds the results of further reconstructions to the "reconstruction" list if several are performed. - if (n_repetitions > 1) { - reconstruction_list[[t]] <- reconstruction - } + # Adds the results of further reconstructions to the "reconstruction" list if several are performed. + if (n_repetitions > 1) { + reconstruction_list[[t]] <- reconstruction + } } if(n_repetitions > 1) { - reconstruction_list - } else { - reconstruction - } + return(reconstruction_list) + } else { + return(reconstruction) } +} diff --git a/R/select_kernel.R b/R/select_kernel.R index 052acb29..774beb6c 100644 --- a/R/select_kernel.R +++ b/R/select_kernel.R @@ -1,6 +1,6 @@ #' select_kernel #' -#' @description select_kernel (internal) +#' @description Kernel selection #' #' @param kernel_arg Parameter of the function #' reconstruct_pattern_multi_trait_marks, specifies the kernel to be used to @@ -15,7 +15,7 @@ #' Returns the function of the selected kernel, which is then used to #' calculate the kernel. #' -#' @return function +#' @return list #' #' @aliases select_kernel #' @rdname select_kernel @@ -59,7 +59,6 @@ select_kernel <- function(kernel_arg, bw, rmax, divisor) { d = function(r, d) stats::dnorm(r,d,sd = bw)/ (2*pi*d) ) }, - cumulative = { rmax_bw <- rmax switch(divisor, @@ -70,7 +69,7 @@ select_kernel <- function(kernel_arg, bw, rmax, divisor) { ) } ) - sel_kernel <- list(kernel, rmax_bw) + list(kernel, rmax_bw) } diff --git a/man/calc_moments.Rd b/man/calc_moments.Rd index cd4ed616..40bc84f5 100644 --- a/man/calc_moments.Rd +++ b/man/calc_moments.Rd @@ -11,28 +11,25 @@ calc_moments(fn, p, exclude = NULL, x, y, mark, kernel, rmax_bw, r) \item{p}{Defines the initial state of the new ponit pattern.} -\item{exclude}{} +\item{exclude}{Vector indicating which values not to use.} -\item{x}{x Coordinate of the points from the reference point pattern.} - -\item{y}{y Coordinate of the points from the reference point pattern.} +\item{x, y}{x and y coordinates of the points from the reference point pattern.} \item{mark}{Marks the currently viewed point pattern.} -\item{kernel}{Result of the kernel calculation, calculated with the -calc_kernels function.} +\item{kernel}{Result of the kernel calculation, calculated with the calc_kernels function.} \item{rmax_bw}{Maximum distance at which the summary statistics are evaluated + Bandwidth with which the kernels are scaled, so that this is the - standard deviation of the smoothing kernel.} +standard deviation of the smoothing kernel.} -\item{r}{is a sequence from rmin to rmax in rcount steps.} +\item{r}{Sequence from rmin to rmax in rcount steps.} } \value{ matrix } \description{ -calc_moments (internal) +Calculate moments } \details{ Definition of the product-moment function for calculating the contribution diff --git a/man/compute_statistics.Rd b/man/compute_statistics.Rd index 6416b35b..b2776c7a 100644 --- a/man/compute_statistics.Rd +++ b/man/compute_statistics.Rd @@ -4,27 +4,37 @@ \alias{compute_statistics} \title{compute_statistics} \usage{ -compute_statistics(x, y, k, xr, yr, w_statistics) +compute_statistics( + x, + y, + k, + xr, + yr, + w_statistics, + bw, + divisor, + kernel_arg, + rmax +) } \arguments{ -\item{x}{x Coordinate of the points from the reference point pattern.} +\item{x, y}{x and y coordinates of the points from the reference point pattern.} -\item{y}{y Coordinate of the points from the reference point pattern.} +\item{k}{Vector of values k; used only if Dk is included in w_statistics below.} -\item{k}{Vector of values k; used only if Dk is included in w_statistics -elow.} +\item{xr, yr}{x and y extension of the observation window (start, end).} -\item{xr}{x extension of the observation window (start, end)} +\item{w_statistics}{vector of named weights for optional spatial statistics +from the \code{spatstat} package to be included in the energy calculation. This may +include Dk, K, Hs, pcf.} -\item{yr}{y extension of the observation window (start, end)} - -\item{w_statistics}{} +\item{bw, divisor, kernel_arg, rmax}{Several parameters related to summary function.} } \value{ list } \description{ -compute_statistics (internal) +Compute summary statistics } \details{ Compute optional spatial statistics using the spatstat package. diff --git a/man/dummy_transf.Rd b/man/dummy_transf.Rd index bfbe52c9..e03c2745 100644 --- a/man/dummy_transf.Rd +++ b/man/dummy_transf.Rd @@ -16,7 +16,7 @@ for the whole point pattern.} matrix } \description{ -dummy_transf (internal) +Tranfsorm to dummy variables } \details{ Function for the transformation of variables to dummy variables and back diff --git a/man/energy_fun.Rd b/man/energy_fun.Rd index 4f38ad92..6bcbbf65 100644 --- a/man/energy_fun.Rd +++ b/man/energy_fun.Rd @@ -30,16 +30,6 @@ the new point pattern.} \item{statistics}{Results of the compute_statistics function for the new point pattern (calculation of optional spatial statistics).} -\item{f_}{Result of the calc_moments_full function which represents -product-moment contribution of a point at coordinates x, y with marks, -for the whole reference point pattern.} - -\item{f0_}{Column sums of the weights of the brand correlation functions of -the reference point pattern.} - -\item{statistics_}{Results of the compute_statistics function for the -reference point pattern (calculation of optional spatial statistics).} - \item{fn}{Determination of the weightings of the mark correlation functions.} \item{p}{Defines the initial state of the new ponit pattern.} @@ -47,22 +37,19 @@ reference point pattern (calculation of optional spatial statistics).} \item{p_}{Reference point pattern.} \item{Lp}{Distance measure for the calculation of the energy function -(Lp distance, 1 ≤ p < ∞).} +(Lp distance, 1 <= p Date: Thu, 14 Dec 2023 16:25:51 +0100 Subject: [PATCH 13/39] Adding tinytex to rcmd check --- .github/workflows/R-CMD-check.yaml | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 415a110b..337c9779 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -33,6 +33,8 @@ jobs: - uses: r-lib/actions/setup-pandoc@v2 + - uses: r-lib/actions/setup-tinytex@v2 + - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} From 74e0b77ac58543c782fc9d81b28051d035723417 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Thu, 14 Dec 2023 16:41:37 +0100 Subject: [PATCH 14/39] :clown: gh-actions --- .github/workflows/R-CMD-check.yaml | 9 +++++---- _pkgdown.yml | 1 + 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 337c9779..16399751 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -17,7 +17,7 @@ jobs: name: ${{ matrix.config.os }} (${{ matrix.config.r }}) strategy: - fail-fast: TRUE + fail-fast: FALSE matrix: config: - {os: macOS-latest, r: 'release'} @@ -33,7 +33,7 @@ jobs: - uses: r-lib/actions/setup-pandoc@v2 - - uses: r-lib/actions/setup-tinytex@v2 + # - uses: r-lib/actions/setup-tinytex@v2 - uses: r-lib/actions/setup-r@v2 with: @@ -48,5 +48,6 @@ jobs: - uses: r-lib/actions/check-r-package@v2 with: - args: 'c("--as-cran")' - error-on: '"note"' + args: 'c("--no-manual", "--as-cran")' + error_on: '"warning"' + upload-snapshots: true diff --git a/_pkgdown.yml b/_pkgdown.yml index 5a8c514a..46034b36 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -16,6 +16,7 @@ reference: - randomize_raster - reconstruct_pattern - reconstruct_pattern_marks + - reconstruct_pattern_multi - translate_raster - title: Analyse results contents: From 17560d7d7dd06c5a1020186a789229bdce54fe04 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Thu, 14 Dec 2023 16:52:14 +0100 Subject: [PATCH 15/39] Still fun with the gh actions --- .github/workflows/R-CMD-check.yaml | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 16399751..5d9732b0 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -33,8 +33,6 @@ jobs: - uses: r-lib/actions/setup-pandoc@v2 - # - uses: r-lib/actions/setup-tinytex@v2 - - uses: r-lib/actions/setup-r@v2 with: r-version: ${{ matrix.config.r }} @@ -48,6 +46,5 @@ jobs: - uses: r-lib/actions/check-r-package@v2 with: - args: 'c("--no-manual", "--as-cran")' - error_on: '"warning"' + args: 'c("--as-cran")' upload-snapshots: true From bf8f37bc80270e7d17b6bc57e510f9ca34dd4a57 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 15 Dec 2023 08:31:49 +0100 Subject: [PATCH 16/39] Wait? Is it the --as-cran? --- .github/workflows/R-CMD-check.yaml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 5d9732b0..ee53424f 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -46,5 +46,4 @@ jobs: - uses: r-lib/actions/check-r-package@v2 with: - args: 'c("--as-cran")' upload-snapshots: true From e48ae3a9a76795474c89e81e402c6deb1d215cea Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 15 Dec 2023 08:55:48 +0100 Subject: [PATCH 17/39] Make tests less verbose --- R/plot.rd_ras.R | 8 -------- tests/testthat/test-list_to_randomized.R | 1 - tests/testthat/test-plot_rd_ras.R | 9 --------- tests/testthat/test-randomize_raster.R | 9 +++++---- tests/testthat/test-reconstruct_pattern_marks.R | 6 +++--- tests/testthat/test-results_habitat_association.R | 3 +-- tests/testthat/test-translate_raster.R | 7 ++++--- 7 files changed, 13 insertions(+), 30 deletions(-) diff --git a/R/plot.rd_ras.R b/R/plot.rd_ras.R index 7c4e96e5..60cf6b74 100644 --- a/R/plot.rd_ras.R +++ b/R/plot.rd_ras.R @@ -50,14 +50,6 @@ plot.rd_ras <- function(x, n = NULL, col, verbose = TRUE, nrow, ncol, ...) { habitats <- sort(table(terra::values(x$observed), useNA = "no")) # get table of habitats - # print warning if more than 10 classes are present - if (length(habitats) > 10) { - - warning("The raster has more than 10 classes. Please make sure discrete classes are provided.", - call. = FALSE) - - } - # get randomized pattern subset_raster <- sample_randomized(randomized = x$randomized, n = n, verbose = verbose) diff --git a/tests/testthat/test-list_to_randomized.R b/tests/testthat/test-list_to_randomized.R index 1c97a819..5f2cf5d2 100644 --- a/tests/testthat/test-list_to_randomized.R +++ b/tests/testthat/test-list_to_randomized.R @@ -78,4 +78,3 @@ testthat::test_that("list_to_randomized works with results_habitat_associations" testthat::expect_s3_class(object = res_b, class = "data.frame") }) - diff --git a/tests/testthat/test-plot_rd_ras.R b/tests/testthat/test-plot_rd_ras.R index eb65151b..97a6f92e 100644 --- a/tests/testthat/test-plot_rd_ras.R +++ b/tests/testthat/test-plot_rd_ras.R @@ -13,9 +13,6 @@ raster_random_ni <- translate_raster(raster = landscape_classified, steps_x = 1:1, steps_y = 1:1, return_input = FALSE, verbose = FALSE) -raster_random_cont <- translate_raster(raster = terra::rast(landscape), - steps_x = 1:1, steps_y = 1:1, verbose = FALSE) - ################################################################################ testthat::test_that("plot returns plot", { @@ -45,9 +42,3 @@ testthat::test_that("plot returns error if wrong id are selected ", { testthat::expect_error(plot(raster_random, n = c(100, 101, 102), verbose = FALSE), regexp = "Please provide at least on valid ID for n.") }) - -testthat::test_that("plot returns warning if more than 10 classes are present", { - - testthat::expect_warning(plot(raster_random_cont, n = 1), - regexp = "The raster has more than 10 classes. Please make sure discrete classes are provided.") -}) diff --git a/tests/testthat/test-randomize_raster.R b/tests/testthat/test-randomize_raster.R index 520eebcf..a4fd5559 100644 --- a/tests/testthat/test-randomize_raster.R +++ b/tests/testthat/test-randomize_raster.R @@ -58,19 +58,20 @@ testthat::test_that("randomize_raster returns error of n_random < 1", { testthat::test_that("randomize_raster returns all warnings", { testthat::expect_warning(randomize_raster(raster = landscape_classified, n_random = 1, - simplify = TRUE), + simplify = TRUE, verbose = FALSE), regexp = "'simplify = TRUE' not possible for 'return_input = TRUE'.") testthat::expect_warning(randomize_raster(raster = landscape_classified, n_random = 2, - simplify = TRUE, return_input = FALSE), + simplify = TRUE, return_input = FALSE, verbose = FALSE), regexp = "'simplify = TRUE' not possible for 'n_random > 1'.") - testthat::expect_warning(randomize_raster(raster = terra::rast(landscape), n_random = 1), + testthat::expect_warning(randomize_raster(raster = terra::rast(landscape), n_random = 1, + verbose = FALSE), regexp = "The raster has more than 10 classes. Please make sure discrete classes are provided.") }) testthat::test_that("Warning if NA are present", { - testthat::expect_warning(randomize_raster(raster = landscape_wrong, n_random = 1), + testthat::expect_warning(randomize_raster(raster = landscape_wrong, n_random = 1, verbose = FALSE), regexp = "NA values present. Please make sure the observation window of the point pattern reflects this.") }) diff --git a/tests/testthat/test-reconstruct_pattern_marks.R b/tests/testthat/test-reconstruct_pattern_marks.R index 2407d352..7364913c 100644 --- a/tests/testthat/test-reconstruct_pattern_marks.R +++ b/tests/testthat/test-reconstruct_pattern_marks.R @@ -104,13 +104,13 @@ testthat::test_that("All warnings are returned for reconstruct_pattern_marks", { testthat::expect_warning(reconstruct_pattern_marks(pattern = pattern_recon, marked_pattern = marks_sub, n_random = 2, max_runs = 1, - return_input = FALSE, - simplify = TRUE), + return_input = FALSE, simplify = TRUE, + verbose = FALSE), regexp = "'simplify = TRUE' not possible for 'n_random > 1'") testthat::expect_warning(reconstruct_pattern_marks(pattern = pattern_recon, marked_pattern = marks_sub, n_random = 1, max_runs = 1, - simplify = TRUE), + simplify = TRUE, verbose = FALSE), regexp = "'simplify = TRUE' not possible for 'return_input = TRUE'") }) diff --git a/tests/testthat/test-results_habitat_association.R b/tests/testthat/test-results_habitat_association.R index ccec218c..580416cc 100644 --- a/tests/testthat/test-results_habitat_association.R +++ b/tests/testthat/test-results_habitat_association.R @@ -82,8 +82,7 @@ testthat::test_that("results_habitat_association returns warning if significance testthat::expect_warning(results_habitat_association(raster = landscape_classified$raster, pattern = random_a, significance_level = 0.75, verbose = FALSE), - regexp = "Make sure 'signifcance_level' is meaningful (e.g. 'significance_level = 0.05').", - fixed = TRUE) + regexp = "Make sure 'signifcance_level' is meaningful \\(e.g. 'significance_level = 0.05'\\).") }) testthat::test_that("results_habitat_association returns warning if more than 25 classes are present", { diff --git a/tests/testthat/test-translate_raster.R b/tests/testthat/test-translate_raster.R index 6a78c588..1fc23938 100644 --- a/tests/testthat/test-translate_raster.R +++ b/tests/testthat/test-translate_raster.R @@ -68,13 +68,14 @@ testthat::test_that("simplify is working for translate_raster", { testthat::test_that("Warning if more than 10 classes are present for translate_raster", { - testthat::expect_warning(translate_raster(raster = terra::rast(landscape), steps_x = 5, steps_y = 5), + testthat::expect_warning(translate_raster(raster = terra::rast(landscape), steps_x = 5, steps_y = 5, + verbose = FALSE), regexp = "The raster has more than 10 classes. Please make sure discrete classes are provided.") }) testthat::test_that("Stop if NA are present", { testthat::expect_error(translate_raster(raster = landscape_wrong, steps_x = 5, steps_y = 5), - regexp = "NA values are not allowed for 'translate_raster()'.", - fixed = TRUE) + regexp = "NA values are not allowed for 'translate_raster\\()'.") }) + From b0ed3dfe6605c2459578c3d34b6aea69ae9bda4f Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 15 Dec 2023 10:00:02 +0100 Subject: [PATCH 18/39] [skip-ci] Make example multi simpler --- R/reconstruct_pattern_multi.R | 37 ++++++++++++----------------------- 1 file changed, 13 insertions(+), 24 deletions(-) diff --git a/R/reconstruct_pattern_multi.R b/R/reconstruct_pattern_multi.R index fabff8ff..6668a556 100644 --- a/R/reconstruct_pattern_multi.R +++ b/R/reconstruct_pattern_multi.R @@ -75,33 +75,22 @@ #' \dontrun{ #' #' # Random example data set -#' xr <- 500 -#' yr <- 1000 -#' N <- 400 -#' y <- runif(N, min = 0, max = yr) -#' x <- runif(N, min = 0, max = xr) -#' species <- as.factor(c("A","B")) -#' species <- sample(species, N, replace = TRUE) -#' diameter <- runif(N, 1, 40) -#' random <- data.frame(x, y, diameter/ 0.1,species) -#' colnames(random) <- c("x", "y", "dbh [mm]", "Tree species") +#' xr <- 500 +#' yr <- 1000 +#' N <- 400 +#' y <- runif(N, min = 0, max = yr) +#' x <- runif(N, min = 0, max = xr) #' -#' # Conversion to a ppp object and conversion of the metric mark to metres. -#' w <- owin(c(0, 500), c(0, 1000)) -#' marked_pattern <- as.ppp(data.frame(random), W = w) -#' marked_pattern$marks$dbh..mm.<-marked_pattern$marks$dbh..mm.*0.001 +#' species <- sample(c("A","B"), N, replace = TRUE) +#' diameter <- runif(N, 0.1, 0.4) +#' +#' random <- data.frame(x = x, y = y, dbh = diameter, species = factor(species)) +#' +#' marked_pattern <- spatstat.geom::as.ppp(random, W = spatstat.geom::owin(c(0, xr), c(0, yr))) #' #' # Reconstruction function -#' reconstruction <- reconstruct_pattern_multi( -#' marked_pattern, n_repetitions = 1, max_steps = 100000, no_change = 5, -#' rcount = 250, rmax = 25, issue = 1000, divisor = "r", -#' kernel_arg = "epanechnikov", timing = TRUE, energy_evaluation = TRUE, -#' plot = FALSE, Lp = 1, bw = 0.5, sd = "step", steps_tol = 1000, -#' tol = 1e-4, w_markcorr = c(d_d=1, all=1, d_all=1, all_all=1, d_d0=1, -#' all0=1, d_all0=1, all_all0=1), prob_of_actions = c(move_coordinate = 0.4, -#' switch_coords = 0.1, exchange_mark_one = 0.1, exchange_mark_two = 0.1, -#' pick_mark_one = 0.2, pick_mark_two = 0.1), -#' k = 1, w_statistics = c()) +#' reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 19, +#' max_steps = 10000) #' } #' #' @aliases reconstruct_pattern_multi From fc6b9acb9475043c2c7741bbc734a451b9633968 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 15 Dec 2023 10:06:48 +0100 Subject: [PATCH 19/39] [skip-ci] Adding verbose to multi reconstruction --- R/reconstruct_pattern_multi.R | 30 ++++++++++++---------- man/reconstruct_pattern_multi.Rd | 44 +++++++++++++------------------- 2 files changed, 35 insertions(+), 39 deletions(-) diff --git a/R/reconstruct_pattern_multi.R b/R/reconstruct_pattern_multi.R index 6668a556..4b70389b 100644 --- a/R/reconstruct_pattern_multi.R +++ b/R/reconstruct_pattern_multi.R @@ -41,6 +41,7 @@ #' @param w_statistics vector of named weights for optional spatial statistics #' from the \code{spatstat} package to be included in the energy calculation. This may #' include Dk, K, Hs, pcf. +#' @param verbose Logical if progress report is printed. #' #' @details #' A novel approach carries out a pattern reconstruction of marked dot patterns @@ -134,7 +135,8 @@ reconstruct_pattern_multi <- function(marked_pattern, w_markcorr = c(d_d=1, all=1, d_all=1, all_all=1, d_d0=1, all0=1, d_all0=1, all_all0=1), prob_of_actions = c(move_coordinate = 0.4, switch_coords = 0.1, exchange_mark_one = 0.1, exchange_mark_two = 0.1, pick_mark_one = 0.2, pick_mark_two = 0.1), k = 1, - w_statistics= c()) { + w_statistics= c(), + verbose = TRUE) { # If several reconstructions are to be carried out, a list is created here in # which the results are then saved continuously. @@ -232,7 +234,7 @@ reconstruct_pattern_multi <- function(marked_pattern, statistics <- compute_statistics(p$x, p$y, k, xwr, ywr, w_statistics) # Prepare the graphical output. - if(plot == TRUE) { + if(plot) { graphics::par(mfrow = 1:2) num_species <- from_dummy(p_$mark[, species, drop = FALSE]) graphics::plot(y~x, p_, pch=19, col= 2L + as.integer(num_species), @@ -269,7 +271,7 @@ reconstruct_pattern_multi <- function(marked_pattern, # Updating the graphical output of all "issue" steps. if (step %% issue == 0) { - if(plot == TRUE) { + if(plot) { graphics::rect(xwr[1], ywr[1], xwr[2], ywr[2], col="white") num_species <- from_dummy(p$mark[, species, drop = FALSE]) @@ -277,14 +279,15 @@ reconstruct_pattern_multi <- function(marked_pattern, graphics::text(p$x, p$y, as.integer(num_species), cex = 0.7) } - # Generates text output with the current calculated values (for example the energy), this is updated every "issue" simulation step. - message("> Progress:", if(n_repetitions > 1) names_reconstruction[[t]], " || iterations: ", step, - " || Simulation progress: ", floor(step/max_steps * 100), "%", - " || energy = ", round(energy, 5), " || energy improvement = ", - energy_improvement, "\t\t\r", appendLF = FALSE) + if (verbose) { + # Generates text output with the current calculated values (for example the energy), this is updated every "issue" simulation step. + message("\r> Progress: ", if(n_repetitions > 1) names_reconstruction[[t]], " || iterations: ", step, + " || Simulation progress: ", floor(step/max_steps * 100), "%", + " || energy = ", round(energy, 5), " || energy improvement = ", + energy_improvement, "\t\t\r", appendLF = FALSE) - Sys.sleep(0) - utils::flush.console() + Sys.sleep(0) + } } # the next code block aborts the reconstruction if the energy decreases by less than tol in "no_change" intervals of steps_tol simulation steps. @@ -292,8 +295,9 @@ reconstruct_pattern_multi <- function(marked_pattern, if(energy > no_change_energy * (1-tol)) { no_change_counter <- no_change_counter + 1L if(no_change_counter == no_change) { - message("the simulation was terminated, because the energy did not - decrease in ", no_change * issue, " simulation steps.") + if (verbose) { + message("the simulation was terminated, because the energy did not decrease in ", no_change * issue, " simulation steps.") + } stop_criterion <- "no_change" break } @@ -455,7 +459,7 @@ reconstruct_pattern_multi <- function(marked_pattern, # End of reconstruction loop. }) -> process.time - message("\n") + if (verbose) message("\n") # Saves all results Transfers them to the "reconstruction" list. p_$species <- from_dummy(p_$mark[, species, drop = FALSE]) diff --git a/man/reconstruct_pattern_multi.Rd b/man/reconstruct_pattern_multi.Rd index a7e3eb97..9dce31d5 100644 --- a/man/reconstruct_pattern_multi.Rd +++ b/man/reconstruct_pattern_multi.Rd @@ -29,7 +29,8 @@ reconstruct_pattern_multi( prob_of_actions = c(move_coordinate = 0.4, switch_coords = 0.1, exchange_mark_one = 0.1, exchange_mark_two = 0.1, pick_mark_one = 0.2, pick_mark_two = 0.1), k = 1, - w_statistics = c() + w_statistics = c(), + verbose = TRUE ) } \arguments{ @@ -93,6 +94,8 @@ exchange_mark_two = 0.1, pick_mark_one = 0.2, pick_mark_two = 0.1)}.} \item{w_statistics}{vector of named weights for optional spatial statistics from the \code{spatstat} package to be included in the energy calculation. This may include Dk, K, Hs, pcf.} + +\item{verbose}{Logical if progress report is printed.} } \value{ list @@ -126,33 +129,22 @@ A comprehensive description of the method can be found in Wudel et al. (2023). \dontrun{ # Random example data set - xr <- 500 - yr <- 1000 - N <- 400 - y <- runif(N, min = 0, max = yr) - x <- runif(N, min = 0, max = xr) - species <- as.factor(c("A","B")) - species <- sample(species, N, replace = TRUE) - diameter <- runif(N, 1, 40) - random <- data.frame(x, y, diameter/ 0.1,species) - colnames(random) <- c("x", "y", "dbh [mm]", "Tree species") - -# Conversion to a ppp object and conversion of the metric mark to metres. - w <- owin(c(0, 500), c(0, 1000)) - marked_pattern <- as.ppp(data.frame(random), W = w) - marked_pattern$marks$dbh..mm.<-marked_pattern$marks$dbh..mm.*0.001 +xr <- 500 +yr <- 1000 +N <- 400 +y <- runif(N, min = 0, max = yr) +x <- runif(N, min = 0, max = xr) + +species <- sample(c("A","B"), N, replace = TRUE) +diameter <- runif(N, 0.1, 0.4) + +random <- data.frame(x = x, y = y, dbh = diameter, species = factor(species)) + +marked_pattern <- spatstat.geom::as.ppp(random, W = spatstat.geom::owin(c(0, xr), c(0, yr))) # Reconstruction function -reconstruction <- reconstruct_pattern_multi( - marked_pattern, n_repetitions = 1, max_steps = 100000, no_change = 5, - rcount = 250, rmax = 25, issue = 1000, divisor = "r", - kernel_arg = "epanechnikov", timing = TRUE, energy_evaluation = TRUE, - plot = FALSE, Lp = 1, bw = 0.5, sd = "step", steps_tol = 1000, - tol = 1e-4, w_markcorr = c(d_d=1, all=1, d_all=1, all_all=1, d_d0=1, - all0=1, d_all0=1, all_all0=1), prob_of_actions = c(move_coordinate = 0.4, - switch_coords = 0.1, exchange_mark_one = 0.1, exchange_mark_two = 0.1, - pick_mark_one = 0.2, pick_mark_two = 0.1), - k = 1, w_statistics = c()) +reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 19, +max_steps = 10000) } } From 00ab8454a09cdeff9a5e62769892a74f300c5322 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 15 Dec 2023 10:35:10 +0100 Subject: [PATCH 20/39] Adding some simple tests for reconstruct_multi --- .../testthat/test-reconstruct_pattern_multi.R | 53 +++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 tests/testthat/test-reconstruct_pattern_multi.R diff --git a/tests/testthat/test-reconstruct_pattern_multi.R b/tests/testthat/test-reconstruct_pattern_multi.R new file mode 100644 index 00000000..e4ac0874 --- /dev/null +++ b/tests/testthat/test-reconstruct_pattern_multi.R @@ -0,0 +1,53 @@ +# testthat::context("test-reconstruct_pattern_multi") + +# create random data +xr <- 500 +yr <- 1000 + +N <- 400 + +y <- runif(N, min = 0, max = yr) +x <- runif(N, min = 0, max = xr) + +species <- sample(c("A","B"), N, replace = TRUE) +diameter <- runif(N, 0.1, 0.4) + +random <- data.frame(x = x, y = y, dbh = diameter, species = factor(species)) + +# Conversion to a ppp object and conversion of the metric mark to metres. +marked_pattern <- spatstat.geom::as.ppp(random, W = spatstat.geom::owin(c(0, 500), c(0, 1000))) + +# Reconstruction function +multi_recon <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 3, max_steps = 10, + verbose = FALSE) + +multi_recon_simple <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 100, + verbose = FALSE) + +################################################################################ + +testthat::test_that("Output is a long as n_random for reconstruct_pattern_multi", { + + testthat::expect_type(multi_recon, type = "list") + + testthat::expect_length(multi_recon, n = 3) +}) + +testthat::test_that("Output includes randomizations and original pattern for reconstruct_pattern_multi", { + + testthat::expect_true(all(sapply(multi_recon, function(i) i$reference == random))) + + testthat::expect_true(all(sapply(multi_recon, function(i) nrow(i$reconstructed) == N))) +}) + + +testthat::test_that("Only one pattern returned for n = 1", { + + testthat::expect_length(multi_recon_simple, n = 14) +}) + +testthat::test_that("Energy decresead for for reconstruct_pattern_multi", { + + testthat::expect_lt(object = multi_recon_simple$energy_current, + expected = multi_recon_simple$energy_launch) +}) From 9003c640533fa40f0cca0761e79ec352620bc829 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 15 Dec 2023 10:48:59 +0100 Subject: [PATCH 21/39] Adding vignette template --- README.Rmd | 2 +- _pkgdown.yml | 9 ++++++--- inst/CITATION | 2 +- vignettes/articles/publication_record.Rmd | 2 +- .../articles/reconstruct_multi_patterns.Rmd | 17 +++++++++++++++++ vignettes/get_started.Rmd | 2 +- 6 files changed, 27 insertions(+), 7 deletions(-) create mode 100644 vignettes/articles/reconstruct_multi_patterns.Rmd diff --git a/README.Rmd b/README.Rmd index be19212b..e5fdb1dd 100644 --- a/README.Rmd +++ b/README.Rmd @@ -57,7 +57,7 @@ The **shar** package is part of our academic work. To cite the package or acknow > 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. . -If you use the `reconstruct_pattern_multi` function, please also cite. +If you use the `reconstruct_pattern_multi()` function, please also cite. > Wudel C., Schlicht R., Berger U. (2023). Multi-trait point pattern reconstruction of plant ecosystems. Methods in Ecology and Evolution, 14, 2668–2679. . diff --git a/_pkgdown.yml b/_pkgdown.yml index 46034b36..57189a46 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -60,12 +60,15 @@ navbar: icon: fa-graduation-cap text: "More" menu: - - text: "Reconstruct multiple pattern" - icon: fa-book - href: articles/reconstruct_multiple_patterns.html - text: "Parallelization" icon: fa-book href: articles/parallelization.html + - text: "Reconstruct multiple patterns" + icon: fa-book + href: articles/reconstruct_multiple_patterns.html + - text: "Multi-trait pattern reconstruction" + icon: fa-book + href: articles/reconstruct_multiple_patterns.html - text: "Analysing the climatic niche of Cormus domestica" icon: fa-book href: articles/cormus_domestica_tmp.html diff --git a/inst/CITATION b/inst/CITATION index fda66fbe..69780568 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -20,6 +20,6 @@ c( volume = "14", pages = "2668–2679", doi = "https://doi.org/10.1111/2041-210X.14206", - header = "If you use the 'reconstruct_pattern_multi' function, please also cite:" + header = "If you use the 'reconstruct_pattern_multi()' function, please also cite:" ) ) diff --git a/vignettes/articles/publication_record.Rmd b/vignettes/articles/publication_record.Rmd index 062eea5c..d44226a1 100644 --- a/vignettes/articles/publication_record.Rmd +++ b/vignettes/articles/publication_record.Rmd @@ -12,7 +12,7 @@ The **shar** package is part of our academic work. To cite the package or acknow > 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 -If you use the `reconstruct_pattern_multi` function, please also cite. +If you use the `reconstruct_pattern_multi()` function, please also cite. > Wudel C., Schlicht R., Berger U. (2023). Multi-trait point pattern reconstruction of plant ecosystems. Methods in Ecology and Evolution, 14, 2668–2679. . diff --git a/vignettes/articles/reconstruct_multi_patterns.Rmd b/vignettes/articles/reconstruct_multi_patterns.Rmd new file mode 100644 index 00000000..bfd23a39 --- /dev/null +++ b/vignettes/articles/reconstruct_multi_patterns.Rmd @@ -0,0 +1,17 @@ +--- +title: "Multi-trait point pattern reconstruction of plant ecosystems" +author: "" +date: "`r Sys.Date()`" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{How to reconstruct multiple patterns} + %\VignetteEncoding{UTF-8} + %\VignetteEngine{knitr::rmarkdown} +editor_options: + chunk_output_type: console +--- + +...Work-in progress... + +If you want more information about multi-trait point pattern reconstruction, please refer to the [corresponding paper](https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.14206 +). diff --git a/vignettes/get_started.Rmd b/vignettes/get_started.Rmd index db04bb6f..a3c5b310 100644 --- a/vignettes/get_started.Rmd +++ b/vignettes/get_started.Rmd @@ -114,7 +114,7 @@ The **shar** package is part of our academic work. To cite the package or acknow > 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. . -If you use the `reconstruct_pattern_multi` function, please also cite. +If you use the `reconstruct_pattern_multi()` function, please also cite. > Wudel C., Schlicht R., Berger U. (2023). Multi-trait point pattern reconstruction of plant ecosystems. Methods in Ecology and Evolution, 14, 2668–2679. . From 92dd3d8c00c06219c07e8a2f2479ceb9a1ccb58c Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 15 Dec 2023 09:50:47 +0000 Subject: [PATCH 22/39] Re-build README.md --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 4b62f918..a2e866e7 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ -README Last updated: 2023-12-14 +README Last updated: 2023-12-15 | CI | Development | CRAN | License | |--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|----------------------------------------------------------------------------------------------------------------------------|-------------------------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------| @@ -69,7 +69,7 @@ or acknowledge its use in publications, please cite the following paper. > Open Source Software, 6(67), 3811. > . -If you use the `reconstruct_pattern_multi` function, please also cite. +If you use the `reconstruct_pattern_multi()` function, please also cite. > Wudel C., Schlicht R., Berger U. (2023). Multi-trait point pattern > reconstruction of plant ecosystems. Methods in Ecology and Evolution, From 0cd90fc2fd16fa88bfe447bbc0b2011ec4eba85e Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Fri, 15 Dec 2023 11:01:40 +0100 Subject: [PATCH 23/39] Fix homepage links --- _pkgdown.yml | 6 +++--- vignettes/articles/reconstruct_multi_patterns.Rmd | 2 +- ...ltiple_patterns.Rmd => reconstruct_several_patterns.Rmd} | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) rename vignettes/articles/{reconstruct_multiple_patterns.Rmd => reconstruct_several_patterns.Rmd} (96%) diff --git a/_pkgdown.yml b/_pkgdown.yml index 57189a46..80eb5fed 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -63,12 +63,12 @@ navbar: - text: "Parallelization" icon: fa-book href: articles/parallelization.html - - text: "Reconstruct multiple patterns" + - text: "Reconstruct several patterns" icon: fa-book - href: articles/reconstruct_multiple_patterns.html + href: articles/reconstruct_several_patterns.html - text: "Multi-trait pattern reconstruction" icon: fa-book - href: articles/reconstruct_multiple_patterns.html + href: articles/reconstruct_multi_patterns.html - text: "Analysing the climatic niche of Cormus domestica" icon: fa-book href: articles/cormus_domestica_tmp.html diff --git a/vignettes/articles/reconstruct_multi_patterns.Rmd b/vignettes/articles/reconstruct_multi_patterns.Rmd index bfd23a39..fc461e9d 100644 --- a/vignettes/articles/reconstruct_multi_patterns.Rmd +++ b/vignettes/articles/reconstruct_multi_patterns.Rmd @@ -4,7 +4,7 @@ author: "" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{How to reconstruct multiple patterns} + %\VignetteIndexEntry{Multi-trait point pattern reconstruction of plant ecosystems} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: diff --git a/vignettes/articles/reconstruct_multiple_patterns.Rmd b/vignettes/articles/reconstruct_several_patterns.Rmd similarity index 96% rename from vignettes/articles/reconstruct_multiple_patterns.Rmd rename to vignettes/articles/reconstruct_several_patterns.Rmd index 751b1066..a23a978e 100644 --- a/vignettes/articles/reconstruct_multiple_patterns.Rmd +++ b/vignettes/articles/reconstruct_several_patterns.Rmd @@ -1,10 +1,10 @@ --- -title: "How to reconstruct multiple patterns" +title: "How to reconstruct several patterns" author: "Maximilian H.K. Hesselbarth" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > - %\VignetteIndexEntry{How to reconstruct multiple patterns} + %\VignetteIndexEntry{How to reconstruct several patterns} %\VignetteEncoding{UTF-8} %\VignetteEngine{knitr::rmarkdown} editor_options: From d9cf61e982ae986a55b8ca3cb8af4fd953250ff7 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 08:31:57 +0100 Subject: [PATCH 24/39] Update code cov --- .github/workflows/Test-coverage.yaml | 21 ++++++++++++++++++++- 1 file changed, 20 insertions(+), 1 deletion(-) diff --git a/.github/workflows/Test-coverage.yaml b/.github/workflows/Test-coverage.yaml index e100d3d7..7e6f85fd 100644 --- a/.github/workflows/Test-coverage.yaml +++ b/.github/workflows/Test-coverage.yaml @@ -31,5 +31,24 @@ jobs: needs: coverage - name: Test coverage - run: covr::codecov(quiet = FALSE) + run: | + covr::codecov( + quiet = FALSE, + clean = FALSE, + install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package") + ) shell: Rscript {0} + + - name: Show testthat output + if: always() + run: | + ## -------------------------------------------------------------------- + find ${{ runner.temp }}/package -name 'testthat.Rout*' -exec cat '{}' \; || true + shell: bash + + - name: Upload test results + if: failure() + uses: actions/upload-artifact@v3 + with: + name: coverage-test-failures + path: ${{ runner.temp }}/package From 5b30eb26116beee94624fea76f34ca5f6295996a Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 08:52:40 +0100 Subject: [PATCH 25/39] Fix bug related to compute_statistics --- R/compute_statistics.R | 2 +- R/reconstruct_pattern_multi.R | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/compute_statistics.R b/R/compute_statistics.R index e3bcf20a..99674a1f 100644 --- a/R/compute_statistics.R +++ b/R/compute_statistics.R @@ -27,7 +27,7 @@ compute_statistics <- function(x, y, k, xr, yr, w_statistics, bw, divisor, kerne # Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. Dk = { nnd_ <- as.matrix(spatstat.geom::nndist(x, y, k=k)) - apply(nnd_, 2, function(z) cumsum(graphics::hist(z[z <= rmax], breaks = c(-Inf, r), plot = FALSE) $ count) / length(z)) + apply(nnd_, 2, function(z) cumsum(graphics::hist(z[z <= rmax], breaks = c(-Inf, rmax), plot = FALSE) $ count) / length(z)) }, # Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. diff --git a/R/reconstruct_pattern_multi.R b/R/reconstruct_pattern_multi.R index 4b70389b..64780ac5 100644 --- a/R/reconstruct_pattern_multi.R +++ b/R/reconstruct_pattern_multi.R @@ -231,7 +231,7 @@ reconstruct_pattern_multi <- function(marked_pattern, f <- calc_moments_full(fn, p, kernel, rmax_bw, r) f0 <- colSums(p$mark[, fn$i] * p$mark[, fn$j]) names(f0) <- rownames(f) - statistics <- compute_statistics(p$x, p$y, k, xwr, ywr, w_statistics) + statistics <- compute_statistics(p$x, p$y, k, xwr, ywr, w_statistics, bw, divisor, kernel_arg, rmax) # Prepare the graphical output. if(plot) { From 8e6f96a255db337e7df29643034e94a66d8e6059 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 09:13:59 +0100 Subject: [PATCH 26/39] Fix bug in reconstruct_multi and compute_stats --- R/compute_statistics.R | 8 ++++---- R/reconstruct_pattern_multi.R | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/R/compute_statistics.R b/R/compute_statistics.R index 99674a1f..1bc7e0b7 100644 --- a/R/compute_statistics.R +++ b/R/compute_statistics.R @@ -8,7 +8,7 @@ #' @param w_statistics vector of named weights for optional spatial statistics #' from the \code{spatstat} package to be included in the energy calculation. This may #' include Dk, K, Hs, pcf. -#' @param bw,divisor,kernel_arg,rmax Several parameters related to summary function. +#' @param bw,divisor,kernel_arg,r Several parameters related to summary function. #' #' @details #' Compute optional spatial statistics using the spatstat package. @@ -19,7 +19,7 @@ #' @rdname compute_statistics #' #' @keywords internal -compute_statistics <- function(x, y, k, xr, yr, w_statistics, bw, divisor, kernel_arg, rmax) { +compute_statistics <- function(x, y, k, xr, yr, w_statistics, bw, divisor, kernel_arg, r) { stat <- names(w_statistics) names(stat) <- stat @@ -27,12 +27,12 @@ compute_statistics <- function(x, y, k, xr, yr, w_statistics, bw, divisor, kerne # Calculation of the Dk(r)-function, if this is to be taken into account for the energy calculation. Dk = { nnd_ <- as.matrix(spatstat.geom::nndist(x, y, k=k)) - apply(nnd_, 2, function(z) cumsum(graphics::hist(z[z <= rmax], breaks = c(-Inf, rmax), plot = FALSE) $ count) / length(z)) + apply(nnd_, 2, function(z) cumsum(graphics::hist(z[z <= max(r)], breaks = c(-Inf, max(r)), plot = FALSE) $ count) / length(z)) }, # Calculation of the K(r)-function, if this is to be taken into account for the energy calculation. K = { - kest<-spatstat.explore::Kest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), rmax=rmax, correction="none") + kest<-spatstat.explore::Kest(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), rmax=max(r), correction="none") kest$un }, diff --git a/R/reconstruct_pattern_multi.R b/R/reconstruct_pattern_multi.R index 64780ac5..cc0eeb9a 100644 --- a/R/reconstruct_pattern_multi.R +++ b/R/reconstruct_pattern_multi.R @@ -227,11 +227,11 @@ reconstruct_pattern_multi <- function(marked_pattern, f_ <- calc_moments_full(fn, p_, kernel, rmax_bw, r) f0_ <- colSums(p_$mark[, fn$i] * p_$mark[, fn$j]) names(f0_) <- rownames(f_) - statistics_<- compute_statistics(p_$x, p_$y, k, xr, yr, w_statistics, bw, divisor, kernel_arg, rmax) + statistics_<- compute_statistics(p_$x, p_$y, k, xr, yr, w_statistics, bw, divisor, kernel_arg, r) f <- calc_moments_full(fn, p, kernel, rmax_bw, r) f0 <- colSums(p$mark[, fn$i] * p$mark[, fn$j]) names(f0) <- rownames(f) - statistics <- compute_statistics(p$x, p$y, k, xwr, ywr, w_statistics, bw, divisor, kernel_arg, rmax) + statistics <- compute_statistics(p$x, p$y, k, xwr, ywr, w_statistics, bw, divisor, kernel_arg, r) # Prepare the graphical output. if(plot) { @@ -329,7 +329,7 @@ reconstruct_pattern_multi <- function(marked_pattern, kernel, rmax_bw, r) + calc_moments(fn, p, i, x, y, mdiff, kernel, rmax_bw, r) f0.new <- f0 - statistics.new <- compute_statistics(replace(p$x, i, x), replace(p$y, i, y), k, xwr, ywr, w_statistics) + statistics.new <- compute_statistics(replace(p$x, i, x), replace(p$y, i, y), k, xwr, ywr, w_statistics, bw, divisor, kernel_arg, r) }, # Swaps the coordinates of two randomly drawn points from the new point pattern, applied in xx% of the trap. switch_coords = { From a5f778fe572722896580536dcb46c89e5f786199 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 09:27:44 +0100 Subject: [PATCH 27/39] Use pcf.ppp() fun --- R/calculate_energy.R | 9 ++++----- R/compute_statistics.R | 2 +- R/plot.rd_pat.R | 3 +-- 3 files changed, 6 insertions(+), 8 deletions(-) diff --git a/R/calculate_energy.R b/R/calculate_energy.R index 20b18d50..38865f93 100644 --- a/R/calculate_energy.R +++ b/R/calculate_energy.R @@ -94,8 +94,7 @@ calculate_energy <- function(pattern, gest_observed <- spatstat.explore::Gest(X = pattern_observed, correction = "none", r = r) - pcf_observed <- spatstat.explore::pcf(X = pattern_observed, - correction = "none", divisor = "d", r = r) + pcf_observed <- spatstat.explore::pcf.ppp(X = pattern_observed, correction = "none", divisor = "d", r = r) # loop through all reconstructed patterns result <- vapply(seq_along(pattern_randomized), function(x) { @@ -103,9 +102,9 @@ calculate_energy <- function(pattern, gest_reconstruction <- spatstat.explore::Gest(X = pattern_randomized[[x]], correction = "none", r = r) - pcf_reconstruction <- spatstat.explore::pcf(X = pattern_randomized[[x]], - correction = "none", divisor = "d", - r = r) + pcf_reconstruction <- spatstat.explore::pcf.ppp(X = pattern_randomized[[x]], + 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/compute_statistics.R b/R/compute_statistics.R index 1bc7e0b7..0a5bb73d 100644 --- a/R/compute_statistics.R +++ b/R/compute_statistics.R @@ -38,7 +38,7 @@ compute_statistics <- function(x, y, k, xr, yr, w_statistics, bw, divisor, kerne # Calculation of the pcf(r)-function (spherical contact distribution), if this is to be taken into account for the energy calculation. pcf = { - pcfest<-spatstat.explore::pcf(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") + pcfest<-spatstat.explore::pcf.ppp(spatstat.geom::ppp(x,y,window=spatstat.geom::owin(xr,yr)), r=c(0,r), kernel=kernel_arg, divisor=divisor, bw=bw, correction="none") pcfest$un }, # Calculation of the Hs(r)-function (pair correlation function), if this is to be taken into account for the energy calculation. diff --git a/R/plot.rd_pat.R b/R/plot.rd_pat.R index 9edba871..a2f87cc4 100644 --- a/R/plot.rd_pat.R +++ b/R/plot.rd_pat.R @@ -82,8 +82,7 @@ plot.rd_pat <- function(x, what = "sf", n = NULL, probs = c(0.025, 0.975), # calculate summary functions gest_result <- spatstat.explore::Gest(pattern[[x]], correction = "none", r = r) - pcf_result <- spatstat.explore::pcf(pattern[[x]], divisor = "d", - correction = "none", r = r) + pcf_result <- spatstat.explore::pcf.ppp(pattern[[x]], divisor = "d", correction = "none", r = r) gest_df <- as.data.frame(gest_result) # conver to df From 68c97824b08f553ca3221bccd2987eb51f1cbd1c Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 09:28:01 +0100 Subject: [PATCH 28/39] Update tests and docs --- man/compute_statistics.Rd | 15 ++------------- tests/testthat/test-reconstruct_pattern_multi.R | 17 +++++++++++++++++ 2 files changed, 19 insertions(+), 13 deletions(-) diff --git a/man/compute_statistics.Rd b/man/compute_statistics.Rd index b2776c7a..de418d8b 100644 --- a/man/compute_statistics.Rd +++ b/man/compute_statistics.Rd @@ -4,18 +4,7 @@ \alias{compute_statistics} \title{compute_statistics} \usage{ -compute_statistics( - x, - y, - k, - xr, - yr, - w_statistics, - bw, - divisor, - kernel_arg, - rmax -) +compute_statistics(x, y, k, xr, yr, w_statistics, bw, divisor, kernel_arg, r) } \arguments{ \item{x, y}{x and y coordinates of the points from the reference point pattern.} @@ -28,7 +17,7 @@ compute_statistics( from the \code{spatstat} package to be included in the energy calculation. This may include Dk, K, Hs, pcf.} -\item{bw, divisor, kernel_arg, rmax}{Several parameters related to summary function.} +\item{bw, divisor, kernel_arg, r}{Several parameters related to summary function.} } \value{ list diff --git a/tests/testthat/test-reconstruct_pattern_multi.R b/tests/testthat/test-reconstruct_pattern_multi.R index e4ac0874..7ec619bf 100644 --- a/tests/testthat/test-reconstruct_pattern_multi.R +++ b/tests/testthat/test-reconstruct_pattern_multi.R @@ -24,6 +24,10 @@ multi_recon <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 3, max_ multi_recon_simple <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 100, verbose = FALSE) +multi_recon_fun <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 100, + verbose = FALSE, w_statistics = c("Dk" = 1, "K" = 0.5, "Hs" = 0.5, "pcf" = 1), + divisor = "d", kernel = "gaussian") + ################################################################################ testthat::test_that("Output is a long as n_random for reconstruct_pattern_multi", { @@ -51,3 +55,16 @@ testthat::test_that("Energy decresead for for reconstruct_pattern_multi", { testthat::expect_lt(object = multi_recon_simple$energy_current, expected = multi_recon_simple$energy_launch) }) + +testthat::test_that("Test additional arguments of reconstruct_pattern_multi", { + + expect_equal(object = multi_recon_fun$Parameter_setting$w_statistics, + expected = c("Dk" = 1, "K" = 0.5, "Hs" = 0.5, "pcf" = 1)) + + expect_equal(object = multi_recon_fun$Parameter_setting$divisor, + expected = "d") + + expect_equal(object = multi_recon_fun$Parameter_setting$kernel_arg, + expected = "gaussian") + +}) From 2eebfe8c24ee96c629a31d925366fcdb69bf8d22 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 10:46:36 +0100 Subject: [PATCH 29/39] Fix DOI citation --- inst/CITATION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/CITATION b/inst/CITATION index 69780568..53574f04 100644 --- a/inst/CITATION +++ b/inst/CITATION @@ -19,7 +19,7 @@ c( year = "2023", volume = "14", pages = "2668–2679", - doi = "https://doi.org/10.1111/2041-210X.14206", + doi = "10.1111/2041-210X.14206", header = "If you use the 'reconstruct_pattern_multi()' function, please also cite:" ) ) From 7f78fba1d70cc7125901b516397ef59128f02296 Mon Sep 17 00:00:00 2001 From: ChrisWudel Date: Mon, 18 Dec 2023 12:03:42 +0100 Subject: [PATCH 30/39] Vignette and function for plotting summary statistics for reconstruct_pattern_multi. --- R/plot_sum_stat.R | 225 ++++++++++++++++++ R/reconstruct_pattern_multi.R | 2 +- .../articles/reconstruct_multi_patterns.Rmd | 78 ++++++ 3 files changed, 304 insertions(+), 1 deletion(-) create mode 100644 R/plot_sum_stat.R diff --git a/R/plot_sum_stat.R b/R/plot_sum_stat.R new file mode 100644 index 00000000..e3bbb087 --- /dev/null +++ b/R/plot_sum_stat.R @@ -0,0 +1,225 @@ +#' plot_sum_stat +#' +#' @description plot_sum_stat +#' +#' @param reconstruction Result list of the dot pattern recostruction with +#' multiple marks. +#' +#' @details +#' Calculates and visualises various summary statistics for the results of +#' multi-marks point pattern reconstruction. +#' +#' @return ggplot object +#' +#' @aliases plot_sum_stat +#' @rdname plot_sum_stat +#' +#' @keywords internal + +plot_sum_stat <- function(reconstruction) { + +library(ggplot2) +library(spatstat) +library(reshape) + +# Data import from the results of point pattern reconstruction. + +k_func_recon<-list() +pcf_func_recon<-list() +Hs_func_recon<-list() +dbh_markcorr_func_recon<-list() +species_markcorr_func_recon<-list() + +rmin <- 0 +divisor <- "r" +kernel_arg <- "epanechnikov" + +if(nchar(names(reconstruction[1])) == 9){ + + n_repetitions <- 1 + rmax <- as.numeric(reconstruction$Parameter_setting$rmax) + rcount <- as.numeric(reconstruction$Parameter_setting$rcount) + rpresented <- min(c(as.numeric(reconstruction$window[2]),as.numeric(reconstruction$window[4]))) + r <- seq(rmin, if(rpresented >= rmax*2 ){rmax*2}else{rpresented}, length.out = rcount) + bw <- as.numeric(reconstruction$Parameter_setting$bw) + + ppp_reference <- as.ppp(reconstruction$reference, reconstruction$window) + }else { + n_repetitions <-reconstruction[[1]]$Parameter_setting$n_repetitions + rmax <- as.numeric(reconstruction[[n_repetitions]]$Parameter_setting$rmax) + rpresented <- min(c(as.numeric(reconstruction[[n_repetitions]]$window[2]),as.numeric(reconstruction[[n_repetitions]]$window[4]))) + rcount <- as.numeric(reconstruction[[n_repetitions]]$Parameter_setting$rcount) + r <- seq(rmin, if(rpresented <= rmax*2 ){rmax*2}else{rpresented}, length.out = rcount) + bw <- as.numeric(reconstruction[[n_repetitions]]$Parameter_setting$bw) + + ppp_reference <- as.ppp(reconstruction$reconstruction_1$reference, reconstruction$reconstruction_1$window) + } + +for (i in seq_len(n_repetitions)) { +message("Progress in the creation of the figures:", i/n_repetitions*100,"%\t\t\r", appendLF = FALSE) +if (n_repetitions > 1) { + ppp_reconstructed<-as.ppp(reconstruction[[i]]$reconstructed, reconstruction[[i]]$window) + }else { + ppp_reconstructed<-as.ppp(reconstruction$reconstructed, reconstruction$window) + } + pcf_func<-spatstat.explore::pcf.ppp(ppp_reconstructed, bw = bw, + kernel=kernel_arg, correction = "none", divisor = divisor, r = r) + pcf_func_recon[[i]]<- pcf_func$un + + k_func <- Kest(ppp_reconstructed, correction="none", r = r) + k_func_recon[[i]]<- k_func$un + + ppp_reconstructed$marks$mark<-ppp_reconstructed$marks$mark[,2] + markcorr_func <- markcorr(ppp_reconstructed, correction = "none", r = r) + + dbh_markcorr_func_recon[[i]] <- markcorr_func$diameter$un + species_markcorr_func_recon[[i]] <- markcorr_func$species$un +} +# Visualising the k function of the results of the point pattern reconstruction. +k_func_recon<-as.data.frame(k_func_recon) +name <- c(sprintf("k_func%03d", seq_len(n_repetitions))) +colnames(k_func_recon) <- name + +k_func <- spatstat.explore::Kest(ppp_reference, correction="none", r = r) + +if(n_repetitions > 1){ + k_func_recon_mean <- rowMeans(k_func_recon) + k_func_all <- data.frame(cbind(k_func$un,k_func_recon, k_func_recon_mean, k_func$r)) + }else { + k_func_all <- data.frame(cbind(k_func$un,k_func_recon, k_func$r)) + } + +colnames(k_func_all)[1] <- c("Reference") +colnames(k_func_all)[length(k_func_all)] <- c("r") + +df_k_func <- na.omit(data.frame(melt(k_func_all,"r"))) + +ggp_k_func_all <- +ggplot(data = df_k_func, aes(x = r, y = value)) + + geom_line(aes(group = variable), col = "grey") + + geom_line(data = subset(df_k_func, variable == "Reference"), col = "black") + + geom_line(data = subset(df_k_func, variable == "k_func_recon_mean"), col = "black",linetype = "dashed") + + geom_vline(xintercept = rmax, linetype='dashed')+ + theme_bw() + + theme_classic() + + theme(panel.border = element_rect(colour = "black", fill=NA, size=1))+ + ggtitle("K-function") + + labs(y = "K(r)", x = "r [m]") + + theme(legend.title = element_blank()) + +# Visualising the Pair correlation functions (pcf) of the results of the point pattern reconstruction. +pcf_func_recon<-as.data.frame(pcf_func_recon) +name <- c(sprintf("pcf_func_recon%03d", seq_len(n_repetitions))) +colnames(pcf_func_recon) <- name + +pcf_func_recon_min <- apply(pcf_func_recon, 1, FUN = min) +pcf_func_recon_max <- apply(pcf_func_recon, 1, FUN = max) + +pcf_func <- spatstat.explore::pcf.ppp(ppp_reference, bw = bw, kernel=kernel_arg, + correction = "none", r = r, divisor = divisor) +if(n_repetitions > 1){ + pcf_func_recon_mean <- rowMeans(pcf_func_recon) + pcf_func_all <- data.frame(cbind(pcf_func$un,pcf_func_recon, pcf_func_recon_min, pcf_func_recon_max, pcf_func_recon_mean, pcf_func$r)) + + }else { + pcf_func_all <- data.frame(cbind(pcf_func$un,pcf_func_recon, pcf_func_recon_min, pcf_func_recon_max, pcf_func$r)) + + } + +colnames(pcf_func_all)[1] <- c("Reference") +colnames(pcf_func_all)[length(pcf_func_all)] <- c("r") + +df_pcf_func <- na.omit(data.frame(melt(pcf_func_all,"r"))) + +ggp_pcf_func_all <- +ggplot(data = df_pcf_func, aes(x = r, y = value)) + + geom_line(aes(group = variable), col = "grey") + + geom_line(data = subset(df_pcf_func, variable == "Reference"), col = "black") + + geom_line(data = subset(df_pcf_func, variable == "pcf_func_recon_mean"), col = "black",linetype = "dashed") + + geom_vline(xintercept = rmax, linetype='dashed')+ + theme_bw() + + theme_classic() + + theme(panel.border = element_rect(colour = "black", fill=NA, size=1))+ + ggtitle("Pair correlation function") + + labs(y = "g(r)", x = "r [m]") + + theme(legend.title = element_blank()) + +# Visualising the mark correlation functions (mcf`s) of the results of the point pattern reconstruction. +dbh_markcorr_func_recon<-as.data.frame(dbh_markcorr_func_recon) +name <- c(sprintf("dbh_markcorr_func_recon%03d", seq_len(n_repetitions))) +colnames(dbh_markcorr_func_recon) <- name + +species_markcorr_func_recon<-as.data.frame(species_markcorr_func_recon) +name <- c(sprintf("species_markcorr_func_recon%03d", seq_len(n_repetitions))) +colnames(species_markcorr_func_recon) <- name + +markcorr_func <- markcorr(ppp_reference, correction = "none", r = r) +dbh_markcorr_func <- markcorr_func$diameter$un +species_markcorr_func <- markcorr_func$species$un + +if(n_repetitions > 1){ + dbh_markcorr_func_recon_mean <- rowMeans(dbh_markcorr_func_recon) + dbh_markcorr_all <- data.frame(cbind(dbh_markcorr_func, + dbh_markcorr_func_recon, + dbh_markcorr_func_recon_mean, + markcorr_func$diameter$r)) + }else { + dbh_markcorr_all <- data.frame(cbind(dbh_markcorr_func, + dbh_markcorr_func_recon, + markcorr_func$diameter$r)) + } + +if(n_repetitions > 1){ + species_markcorr_func_recon_mean <- rowMeans(species_markcorr_func_recon) + species_markcorr_all <- data.frame(cbind(species_markcorr_func, + species_markcorr_func_recon, + species_markcorr_func_recon_mean, + markcorr_func$diameter$r)) + }else { + species_markcorr_all <- data.frame(cbind(species_markcorr_func, + species_markcorr_func_recon, + markcorr_func$diameter$r)) + } + +colnames(dbh_markcorr_all)[1] <- c("Reference mark dbh") +colnames(dbh_markcorr_all)[length(dbh_markcorr_all)] <- c("r") + +colnames(species_markcorr_all)[1] <- c("Reference mark species") +colnames(species_markcorr_all)[length(species_markcorr_all)] <- c("r") + +df_dbh_markcorr_func_all <- na.omit(data.frame(melt(dbh_markcorr_all,"r"))) +df_species_markcorr_func_all <- na.omit(data.frame(melt(species_markcorr_all,"r"))) + +ggp_dbh_markcorr_func_all <- + ggplot(data = df_dbh_markcorr_func_all, aes(x = r, y = value)) + + geom_line(aes(group = variable ), col = "grey") + + geom_line(data = subset(df_dbh_markcorr_func_all, variable == "Reference mark dbh"), col = "black") + + geom_line(data = subset(df_dbh_markcorr_func_all, variable == "dbh_markcorr_func_recon_mean"), col = "black",linetype = "dashed") + + geom_vline(xintercept = rmax, linetype='dashed')+ + theme_bw() + + theme_classic() + + theme(panel.border = element_rect(colour = "black", fill=NA, size=1))+ + ggtitle("Mark Correlation Function (dbh)") + + labs(y = "kmm(r)",x = "r [m]") + + theme(legend.title = element_blank()) + + ggp_species_markcorr_func_all <- + ggplot(data = df_species_markcorr_func_all, aes(x = r, y = value)) + + geom_line(aes(group = variable), col = "grey") + + geom_line(data = subset(df_species_markcorr_func_all, variable == "Reference mark species"), col = "black") + + geom_line(data = subset(df_species_markcorr_func_all, variable == "species_markcorr_func_recon_mean"), col = "black",linetype = "dashed") + + geom_vline(xintercept = rmax, linetype='dashed')+ + theme_bw() + + theme_classic() + + theme(panel.border = element_rect(colour = "black", fill=NA, size=1))+ + ggtitle("Mark Correlation Function (species)") + + labs(y = "kmm(r)",x = "r [m]") + + theme(legend.title = element_blank()) + +result <- list(ggp_k_func_all, ggp_pcf_func_all, ggp_species_markcorr_func_all,ggp_dbh_markcorr_func_all) + +cat(sep="\n\n") +print("look under Plots to see the result.") + +return(result) +} diff --git a/R/reconstruct_pattern_multi.R b/R/reconstruct_pattern_multi.R index cc0eeb9a..72849ca9 100644 --- a/R/reconstruct_pattern_multi.R +++ b/R/reconstruct_pattern_multi.R @@ -90,7 +90,7 @@ #' marked_pattern <- spatstat.geom::as.ppp(random, W = spatstat.geom::owin(c(0, xr), c(0, yr))) #' #' # Reconstruction function -#' reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 19, +#' reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 2, #' max_steps = 10000) #' } #' diff --git a/vignettes/articles/reconstruct_multi_patterns.Rmd b/vignettes/articles/reconstruct_multi_patterns.Rmd index fc461e9d..d57f20d7 100644 --- a/vignettes/articles/reconstruct_multi_patterns.Rmd +++ b/vignettes/articles/reconstruct_multi_patterns.Rmd @@ -15,3 +15,81 @@ editor_options: If you want more information about multi-trait point pattern reconstruction, please refer to the [corresponding paper](https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.14206 ). + +If you wish to include several marks simultaneously in a reconstruction, you can +use the following code. The libraries used must first be loaded. + +```{r} +library(shar) +library(spatstat) +``` + +The next step is to load the point pattern, here is an example of a random point +pattern with several marks to show the structure of the data used. + +```{r} +xr <- 500 +yr <- 1000 +N <- 400 +y <- runif(N, min = 0, max = yr) +x <- runif(N, min = 0, max = xr) +species <- sample(c("A","B"), N, replace = TRUE) +diameter <- runif(N, 0.1, 0.4) +random <- data.frame(x = x, y = y, dbh = diameter, species = factor(species)) +marked_pattern <- as.ppp(random, W = owin(c(0, xr), c(0, yr))) +``` + +The point pattern must contain the following data An x and y coordinate, a +metric mark (in metres) and a nominal mark defined as a factor. The order must +be respected. Now the reconstruction with several marks can be started with the +following code. Note that the maximum number of iterations has been set to +max_steps = 10000 to keep the computation time for this example to a minimum. +For an application, this value should be increased according to the number of +points in the pattern. + +```{r} +reconstruction <- reconstruct_pattern_multi( +marked_pattern, +n_repetitions = 1, +max_steps = 10000) +``` + +As a result, you will receive a list containing a variety of information, for +example, the reference pattern, the reconstructed pattern, the number of +successful actions, the energy development and much more. If you wish to +perform several reconstructions of the same reference pattern, you must +increase n_repetitions to the desired number. + +```{r} +reconstruction_2 <- reconstruct_pattern_multi( +marked_pattern, +n_repetitions = 2, +max_steps = 10000) +``` + +To activate a visualisation of the reconstruction that shows the changes in the +pattern at the relevant time, you must proceed as follows. + +```{r} +reconstruction_3 <- reconstruct_pattern_multi( +marked_pattern, +n_repetitions = 1, +max_steps = 10000, +plot = TRUE) +``` + +Finally, you can use the following function to view different summary statistics +of the reference pattern (black line) compared to the reconstructed pattern +(grey line). For this, however, the listed libraries must be loaded +additionally. + +```{r} +library(ggplot2) +library(reshape) + +plot_sum_stat(reconstruction) +``` + + + + From ce291cb0cd8284b7b5c8dc4292dc70ab295dd5fe Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 14:51:58 +0100 Subject: [PATCH 31/39] Fix plot_sum_stat (had to remove ggplot2) and upadate vignette --- NAMESPACE | 1 + R/plot_sum_stat.R | 304 +++++++----------- man/plot_sum_stat.Rd | 22 ++ man/reconstruct_pattern_multi.Rd | 2 +- .../articles/reconstruct_multi_patterns.Rmd | 9 +- 5 files changed, 150 insertions(+), 188 deletions(-) create mode 100644 man/plot_sum_stat.Rd diff --git a/NAMESPACE b/NAMESPACE index bd1b9e47..d3e1296f 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -12,6 +12,7 @@ export(fit_point_process) export(list_to_randomized) export(pack_randomized) export(plot_energy) +export(plot_sum_stat) export(randomize_raster) export(reconstruct_pattern) export(reconstruct_pattern_marks) diff --git a/R/plot_sum_stat.R b/R/plot_sum_stat.R index e3bbb087..dda0991e 100644 --- a/R/plot_sum_stat.R +++ b/R/plot_sum_stat.R @@ -2,7 +2,7 @@ #' #' @description plot_sum_stat #' -#' @param reconstruction Result list of the dot pattern recostruction with +#' @param reconstruction Result list of the dot pattern reconstruction with #' multiple marks. #' #' @details @@ -14,212 +14,158 @@ #' @aliases plot_sum_stat #' @rdname plot_sum_stat #' -#' @keywords internal - +#' @export plot_sum_stat <- function(reconstruction) { -library(ggplot2) -library(spatstat) -library(reshape) - -# Data import from the results of point pattern reconstruction. + # Data import from the results of point pattern reconstruction. + k_func_recon <- list() + pcf_func_recon <- list() + dbh_markcorr_func_recon <- list() + species_markcorr_func_recon <- list() -k_func_recon<-list() -pcf_func_recon<-list() -Hs_func_recon<-list() -dbh_markcorr_func_recon<-list() -species_markcorr_func_recon<-list() + rmin <- 0 + divisor <- "r" + kernel_arg <- "epanechnikov" -rmin <- 0 -divisor <- "r" -kernel_arg <- "epanechnikov" + if(nchar(names(reconstruction[1])) == 9){ -if(nchar(names(reconstruction[1])) == 9){ + n_repetitions <- 1 + rmax <- as.numeric(reconstruction$Parameter_setting$rmax) + rcount <- as.numeric(reconstruction$Parameter_setting$rcount) + rpresented <- min(c(as.numeric(reconstruction$window[2]),as.numeric(reconstruction$window[4]))) + r <- seq(rmin, if(rpresented >= rmax*2 ){rmax*2}else{rpresented}, length.out = rcount) + bw <- as.numeric(reconstruction$Parameter_setting$bw) - n_repetitions <- 1 - rmax <- as.numeric(reconstruction$Parameter_setting$rmax) - rcount <- as.numeric(reconstruction$Parameter_setting$rcount) - rpresented <- min(c(as.numeric(reconstruction$window[2]),as.numeric(reconstruction$window[4]))) - r <- seq(rmin, if(rpresented >= rmax*2 ){rmax*2}else{rpresented}, length.out = rcount) - bw <- as.numeric(reconstruction$Parameter_setting$bw) + ppp_reference <- spatstat.geom::as.ppp(reconstruction$reference, reconstruction$window) + reconstruction <- list(reconstruction) + } else { - ppp_reference <- as.ppp(reconstruction$reference, reconstruction$window) - }else { - n_repetitions <-reconstruction[[1]]$Parameter_setting$n_repetitions + n_repetitions <- reconstruction[[1]]$Parameter_setting$n_repetitions rmax <- as.numeric(reconstruction[[n_repetitions]]$Parameter_setting$rmax) rpresented <- min(c(as.numeric(reconstruction[[n_repetitions]]$window[2]),as.numeric(reconstruction[[n_repetitions]]$window[4]))) rcount <- as.numeric(reconstruction[[n_repetitions]]$Parameter_setting$rcount) r <- seq(rmin, if(rpresented <= rmax*2 ){rmax*2}else{rpresented}, length.out = rcount) bw <- as.numeric(reconstruction[[n_repetitions]]$Parameter_setting$bw) - ppp_reference <- as.ppp(reconstruction$reconstruction_1$reference, reconstruction$reconstruction_1$window) - } + ppp_reference <- spatstat.geom::as.ppp(reconstruction$reconstruction_1$reference, reconstruction$reconstruction_1$window) + } -for (i in seq_len(n_repetitions)) { -message("Progress in the creation of the figures:", i/n_repetitions*100,"%\t\t\r", appendLF = FALSE) -if (n_repetitions > 1) { - ppp_reconstructed<-as.ppp(reconstruction[[i]]$reconstructed, reconstruction[[i]]$window) - }else { - ppp_reconstructed<-as.ppp(reconstruction$reconstructed, reconstruction$window) - } - pcf_func<-spatstat.explore::pcf.ppp(ppp_reconstructed, bw = bw, - kernel=kernel_arg, correction = "none", divisor = divisor, r = r) - pcf_func_recon[[i]]<- pcf_func$un + for (i in seq_len(n_repetitions)) { - k_func <- Kest(ppp_reconstructed, correction="none", r = r) - k_func_recon[[i]]<- k_func$un + message("Progress in the creation of the figures: ", i/n_repetitions*100,"%\t\t\r", appendLF = FALSE) - ppp_reconstructed$marks$mark<-ppp_reconstructed$marks$mark[,2] - markcorr_func <- markcorr(ppp_reconstructed, correction = "none", r = r) + ppp_reconstructed <- spatstat.geom::as.ppp(reconstruction[[i]][[2]], reconstruction[[i]][[3]]) - dbh_markcorr_func_recon[[i]] <- markcorr_func$diameter$un - species_markcorr_func_recon[[i]] <- markcorr_func$species$un -} -# Visualising the k function of the results of the point pattern reconstruction. -k_func_recon<-as.data.frame(k_func_recon) -name <- c(sprintf("k_func%03d", seq_len(n_repetitions))) -colnames(k_func_recon) <- name + pcf_func <-spatstat.explore::pcf.ppp(ppp_reconstructed, bw = bw, kernel=kernel_arg, + correction = "none", divisor = divisor, r = r) + pcf_func_recon[[i]] <- pcf_func$un -k_func <- spatstat.explore::Kest(ppp_reference, correction="none", r = r) + k_func <- spatstat.explore::Kest(ppp_reconstructed, correction="none", r = r) + k_func_recon[[i]]<- k_func$un -if(n_repetitions > 1){ - k_func_recon_mean <- rowMeans(k_func_recon) - k_func_all <- data.frame(cbind(k_func$un,k_func_recon, k_func_recon_mean, k_func$r)) - }else { - k_func_all <- data.frame(cbind(k_func$un,k_func_recon, k_func$r)) - } - -colnames(k_func_all)[1] <- c("Reference") -colnames(k_func_all)[length(k_func_all)] <- c("r") - -df_k_func <- na.omit(data.frame(melt(k_func_all,"r"))) - -ggp_k_func_all <- -ggplot(data = df_k_func, aes(x = r, y = value)) + - geom_line(aes(group = variable), col = "grey") + - geom_line(data = subset(df_k_func, variable == "Reference"), col = "black") + - geom_line(data = subset(df_k_func, variable == "k_func_recon_mean"), col = "black",linetype = "dashed") + - geom_vline(xintercept = rmax, linetype='dashed')+ - theme_bw() + - theme_classic() + - theme(panel.border = element_rect(colour = "black", fill=NA, size=1))+ - ggtitle("K-function") + - labs(y = "K(r)", x = "r [m]") + - theme(legend.title = element_blank()) - -# Visualising the Pair correlation functions (pcf) of the results of the point pattern reconstruction. -pcf_func_recon<-as.data.frame(pcf_func_recon) -name <- c(sprintf("pcf_func_recon%03d", seq_len(n_repetitions))) -colnames(pcf_func_recon) <- name - -pcf_func_recon_min <- apply(pcf_func_recon, 1, FUN = min) -pcf_func_recon_max <- apply(pcf_func_recon, 1, FUN = max) - -pcf_func <- spatstat.explore::pcf.ppp(ppp_reference, bw = bw, kernel=kernel_arg, - correction = "none", r = r, divisor = divisor) -if(n_repetitions > 1){ - pcf_func_recon_mean <- rowMeans(pcf_func_recon) - pcf_func_all <- data.frame(cbind(pcf_func$un,pcf_func_recon, pcf_func_recon_min, pcf_func_recon_max, pcf_func_recon_mean, pcf_func$r)) + ppp_reconstructed$marks$mark <- ppp_reconstructed$marks$mark[,2] + markcorr_func <- spatstat.explore::markcorr(ppp_reconstructed, correction = "none", r = r) - }else { - pcf_func_all <- data.frame(cbind(pcf_func$un,pcf_func_recon, pcf_func_recon_min, pcf_func_recon_max, pcf_func$r)) + dbh_markcorr_func_recon[[i]] <- markcorr_func$diameter$un + species_markcorr_func_recon[[i]] <- markcorr_func$species$un + } - } + # Visualising the k function of the results of the point pattern reconstruction. + k_func_recon <- as.data.frame(k_func_recon) + name <- c(sprintf("k_func%03d", seq_len(n_repetitions))) + colnames(k_func_recon) <- name -colnames(pcf_func_all)[1] <- c("Reference") -colnames(pcf_func_all)[length(pcf_func_all)] <- c("r") + k_func <- spatstat.explore::Kest(ppp_reference, correction="none", r = r) + k_func_recon_mean <- rowMeans(k_func_recon) + k_func_all <- data.frame(cbind(k_func$un, k_func_recon, k_func_recon_mean, k_func$r)) -df_pcf_func <- na.omit(data.frame(melt(pcf_func_all,"r"))) + colnames(k_func_all)[1] <- c("Reference") + colnames(k_func_all)[length(k_func_all)] <- c("r") -ggp_pcf_func_all <- -ggplot(data = df_pcf_func, aes(x = r, y = value)) + - geom_line(aes(group = variable), col = "grey") + - geom_line(data = subset(df_pcf_func, variable == "Reference"), col = "black") + - geom_line(data = subset(df_pcf_func, variable == "pcf_func_recon_mean"), col = "black",linetype = "dashed") + - geom_vline(xintercept = rmax, linetype='dashed')+ - theme_bw() + - theme_classic() + - theme(panel.border = element_rect(colour = "black", fill=NA, size=1))+ - ggtitle("Pair correlation function") + - labs(y = "g(r)", x = "r [m]") + - theme(legend.title = element_blank()) + # Visualising the Pair correlation functions (pcf) of the results of the point pattern reconstruction. + pcf_func_recon <- as.data.frame(pcf_func_recon) + name <- c(sprintf("pcf_func_recon%03d", seq_len(n_repetitions))) + colnames(pcf_func_recon) <- name -# Visualising the mark correlation functions (mcf`s) of the results of the point pattern reconstruction. -dbh_markcorr_func_recon<-as.data.frame(dbh_markcorr_func_recon) -name <- c(sprintf("dbh_markcorr_func_recon%03d", seq_len(n_repetitions))) -colnames(dbh_markcorr_func_recon) <- name + pcf_func <- spatstat.explore::pcf.ppp(ppp_reference, bw = bw, kernel=kernel_arg, + correction = "none", r = r, divisor = divisor) + pcf_func_recon_mean <- rowMeans(pcf_func_recon) + pcf_func_all <- data.frame(cbind(pcf_func$un, pcf_func_recon, pcf_func_recon_mean, pcf_func$r)) -species_markcorr_func_recon<-as.data.frame(species_markcorr_func_recon) -name <- c(sprintf("species_markcorr_func_recon%03d", seq_len(n_repetitions))) -colnames(species_markcorr_func_recon) <- name + colnames(pcf_func_all)[1] <- c("Reference") + colnames(pcf_func_all)[length(pcf_func_all)] <- c("r") -markcorr_func <- markcorr(ppp_reference, correction = "none", r = r) -dbh_markcorr_func <- markcorr_func$diameter$un -species_markcorr_func <- markcorr_func$species$un + # Visualising the mark correlation functions (mcf`s) of the results of the point pattern reconstruction. + dbh_markcorr_func_recon<-as.data.frame(dbh_markcorr_func_recon) + name <- c(sprintf("dbh_markcorr_func_recon%03d", seq_len(n_repetitions))) + colnames(dbh_markcorr_func_recon) <- name + + species_markcorr_func_recon <- as.data.frame(species_markcorr_func_recon) + name <- c(sprintf("species_markcorr_func_recon%03d", seq_len(n_repetitions))) + colnames(species_markcorr_func_recon) <- name + + markcorr_func <- spatstat.explore::markcorr(ppp_reference, correction = "none", r = r) + dbh_markcorr_func <- markcorr_func$diameter$un + species_markcorr_func <- markcorr_func$species$un -if(n_repetitions > 1){ dbh_markcorr_func_recon_mean <- rowMeans(dbh_markcorr_func_recon) - dbh_markcorr_all <- data.frame(cbind(dbh_markcorr_func, - dbh_markcorr_func_recon, - dbh_markcorr_func_recon_mean, - markcorr_func$diameter$r)) - }else { - dbh_markcorr_all <- data.frame(cbind(dbh_markcorr_func, - dbh_markcorr_func_recon, - markcorr_func$diameter$r)) - } - -if(n_repetitions > 1){ + dbh_markcorr_all <- data.frame(cbind(dbh_markcorr_func, dbh_markcorr_func_recon, + dbh_markcorr_func_recon_mean, markcorr_func$diameter$r)) + species_markcorr_func_recon_mean <- rowMeans(species_markcorr_func_recon) - species_markcorr_all <- data.frame(cbind(species_markcorr_func, - species_markcorr_func_recon, - species_markcorr_func_recon_mean, - markcorr_func$diameter$r)) - }else { - species_markcorr_all <- data.frame(cbind(species_markcorr_func, - species_markcorr_func_recon, - markcorr_func$diameter$r)) - } - -colnames(dbh_markcorr_all)[1] <- c("Reference mark dbh") -colnames(dbh_markcorr_all)[length(dbh_markcorr_all)] <- c("r") - -colnames(species_markcorr_all)[1] <- c("Reference mark species") -colnames(species_markcorr_all)[length(species_markcorr_all)] <- c("r") - -df_dbh_markcorr_func_all <- na.omit(data.frame(melt(dbh_markcorr_all,"r"))) -df_species_markcorr_func_all <- na.omit(data.frame(melt(species_markcorr_all,"r"))) - -ggp_dbh_markcorr_func_all <- - ggplot(data = df_dbh_markcorr_func_all, aes(x = r, y = value)) + - geom_line(aes(group = variable ), col = "grey") + - geom_line(data = subset(df_dbh_markcorr_func_all, variable == "Reference mark dbh"), col = "black") + - geom_line(data = subset(df_dbh_markcorr_func_all, variable == "dbh_markcorr_func_recon_mean"), col = "black",linetype = "dashed") + - geom_vline(xintercept = rmax, linetype='dashed')+ - theme_bw() + - theme_classic() + - theme(panel.border = element_rect(colour = "black", fill=NA, size=1))+ - ggtitle("Mark Correlation Function (dbh)") + - labs(y = "kmm(r)",x = "r [m]") + - theme(legend.title = element_blank()) - - ggp_species_markcorr_func_all <- - ggplot(data = df_species_markcorr_func_all, aes(x = r, y = value)) + - geom_line(aes(group = variable), col = "grey") + - geom_line(data = subset(df_species_markcorr_func_all, variable == "Reference mark species"), col = "black") + - geom_line(data = subset(df_species_markcorr_func_all, variable == "species_markcorr_func_recon_mean"), col = "black",linetype = "dashed") + - geom_vline(xintercept = rmax, linetype='dashed')+ - theme_bw() + - theme_classic() + - theme(panel.border = element_rect(colour = "black", fill=NA, size=1))+ - ggtitle("Mark Correlation Function (species)") + - labs(y = "kmm(r)",x = "r [m]") + - theme(legend.title = element_blank()) - -result <- list(ggp_k_func_all, ggp_pcf_func_all, ggp_species_markcorr_func_all,ggp_dbh_markcorr_func_all) - -cat(sep="\n\n") -print("look under Plots to see the result.") - -return(result) + species_markcorr_all <- data.frame(cbind(species_markcorr_func, species_markcorr_func_recon, + species_markcorr_func_recon_mean, markcorr_func$diameter$r)) + + colnames(dbh_markcorr_all)[1] <- c("Reference_mark_dbh") + colnames(dbh_markcorr_all)[length(dbh_markcorr_all)] <- c("r") + + colnames(species_markcorr_all)[1] <- c("Reference_mark_species") + colnames(species_markcorr_all)[length(species_markcorr_all)] <- c("r") + + graphics::par(mfrow = c(2, 2)) + + # plot Kest + graphics::plot(NULL, xlim = range(r), ylim = range(c(k_func_all$Reference, + k_func_all$k_func_recon_mean)), + main = "K-function", xlab = "r", ylab = "K(r)") + + graphics::lines(x = k_func_all$r, y = k_func_all$Reference, lty = "dashed") + graphics::lines(x = k_func_all$r, y = k_func_all$k_func_recon_mean) + graphics::abline(v = rmax, lty = "dotted", col = "grey") + + graphics::legend(x = "topleft", legend = c("Reference", "Reconstructed"), + lty = c(2, 1), inset = 0.015) + + # plot pcf + graphics::plot(NULL, xlim = range(r), ylim = range(c(pcf_func_all$Reference, + pcf_func_all$pcf_func_recon_mean), + finite = TRUE), + main = "Pair correlation function", xlab = "r", ylab = "g(r)") + + graphics::lines(x = pcf_func_all$r, y = pcf_func_all$Reference, lty = "dashed") + graphics::lines(x = pcf_func_all$r, y = pcf_func_all$pcf_func_recon_mean) + graphics::abline(v = rmax, lty = "dotted", col = "grey") + + # # ask user to hit enter + # graphics::par(ask = TRUE) + + # plot kmm(r) + graphics::plot(NULL, xlim = range(r), ylim = range(c(dbh_markcorr_all$Reference_mark_dbh, + dbh_markcorr_all$dbh_markcorr_func_recon_mean)), + main = "Mark Correlation Function (dbh)", xlab = "r", ylab = "kmm(r)") + + graphics::lines(x = dbh_markcorr_all$r, y = dbh_markcorr_all$Reference_mark_dbh, lty = "dashed") + graphics::lines(x = dbh_markcorr_all$r, y = dbh_markcorr_all$dbh_markcorr_func_recon_mean) + graphics::abline(v = rmax, lty = "dotted", col = "grey") + + graphics::plot(NULL, xlim = range(r), ylim = range(c(species_markcorr_all$Reference_mark_species, + species_markcorr_all$species_markcorr_func_recon_mean)), + main = "Mark Correlation Function (species)", xlab = "r", ylab = "kmm(r)") + + graphics::lines(x = species_markcorr_all$r, y = species_markcorr_all$Reference_mark_species, lty = "dashed") + graphics::lines(x = species_markcorr_all$r, y = species_markcorr_all$species_markcorr_func_recon_mean) + graphics::abline(v = rmax, lty = "dotted", col = "grey") + + graphics::par(mfrow = c(1, 1)) + } diff --git a/man/plot_sum_stat.Rd b/man/plot_sum_stat.Rd new file mode 100644 index 00000000..44db8913 --- /dev/null +++ b/man/plot_sum_stat.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot_sum_stat.R +\name{plot_sum_stat} +\alias{plot_sum_stat} +\title{plot_sum_stat} +\usage{ +plot_sum_stat(reconstruction) +} +\arguments{ +\item{reconstruction}{Result list of the dot pattern reconstruction with +multiple marks.} +} +\value{ +ggplot object +} +\description{ +plot_sum_stat +} +\details{ +Calculates and visualises various summary statistics for the results of +multi-marks point pattern reconstruction. +} diff --git a/man/reconstruct_pattern_multi.Rd b/man/reconstruct_pattern_multi.Rd index 9dce31d5..dd652c5c 100644 --- a/man/reconstruct_pattern_multi.Rd +++ b/man/reconstruct_pattern_multi.Rd @@ -143,7 +143,7 @@ random <- data.frame(x = x, y = y, dbh = diameter, species = factor(species)) marked_pattern <- spatstat.geom::as.ppp(random, W = spatstat.geom::owin(c(0, xr), c(0, yr))) # Reconstruction function -reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 19, +reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 2, max_steps = 10000) } diff --git a/vignettes/articles/reconstruct_multi_patterns.Rmd b/vignettes/articles/reconstruct_multi_patterns.Rmd index d57f20d7..670c9cf2 100644 --- a/vignettes/articles/reconstruct_multi_patterns.Rmd +++ b/vignettes/articles/reconstruct_multi_patterns.Rmd @@ -19,7 +19,7 @@ If you want more information about multi-trait point pattern reconstruction, ple If you wish to include several marks simultaneously in a reconstruction, you can use the following code. The libraries used must first be loaded. -```{r} +```{r message = FALSE, warning = FALSE} library(shar) library(spatstat) ``` @@ -84,12 +84,5 @@ of the reference pattern (black line) compared to the reconstructed pattern additionally. ```{r} -library(ggplot2) -library(reshape) - plot_sum_stat(reconstruction) ``` - - - - From 2557c5f9cc3607e997bdd963bfa04798fcf7c518 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 14:56:13 +0100 Subject: [PATCH 32/39] [skip-ci] Update NEWS --- NEWS.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 2bd619cc..06bc32f2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,7 @@ # shar 2.2 * Improvements - * Added a new function `reconstruct_pattern_multi` + * Added a new function `reconstruct_pattern_multi()` including several internal functions + * Added new function `plot_sum_stat()` # shar 2.1.1 * Bugfixes (thanks to @baddstats) From 0c78761a5734d41a3643a1166f939b475f6fea73 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 15:16:27 +0100 Subject: [PATCH 33/39] Update plot for multi --- NAMESPACE | 2 +- NEWS.md | 3 +- R/{plot_sum_stat.R => plot.rd_multi.R} | 47 +++++++++++-------- R/reconstruct_pattern_multi.R | 4 +- man/plot.rd_multi.Rd | 29 ++++++++++++ man/plot_sum_stat.Rd | 22 --------- man/reconstruct_pattern_multi.Rd | 2 +- tests/testthat/test-plot_rd_multi.R | 34 ++++++++++++++ .../articles/reconstruct_multi_patterns.Rmd | 2 +- 9 files changed, 97 insertions(+), 48 deletions(-) rename R/{plot_sum_stat.R => plot.rd_multi.R} (81%) create mode 100644 man/plot.rd_multi.Rd delete mode 100644 man/plot_sum_stat.Rd create mode 100644 tests/testthat/test-plot_rd_multi.R diff --git a/NAMESPACE b/NAMESPACE index d3e1296f..65fcb239 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand S3method(plot,rd_mar) +S3method(plot,rd_multi) S3method(plot,rd_pat) S3method(plot,rd_ras) S3method(print,rd_mar) @@ -12,7 +13,6 @@ export(fit_point_process) export(list_to_randomized) export(pack_randomized) export(plot_energy) -export(plot_sum_stat) export(randomize_raster) export(reconstruct_pattern) export(reconstruct_pattern_marks) diff --git a/NEWS.md b/NEWS.md index 06bc32f2..380f48ca 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,7 +1,6 @@ # shar 2.2 * Improvements - * Added a new function `reconstruct_pattern_multi()` including several internal functions - * Added new function `plot_sum_stat()` + * Added a new function `reconstruct_pattern_multi()` including several internal functions and methods # shar 2.1.1 * Bugfixes (thanks to @baddstats) diff --git a/R/plot_sum_stat.R b/R/plot.rd_multi.R similarity index 81% rename from R/plot_sum_stat.R rename to R/plot.rd_multi.R index dda0991e..8fef96a0 100644 --- a/R/plot_sum_stat.R +++ b/R/plot.rd_multi.R @@ -2,20 +2,25 @@ #' #' @description plot_sum_stat #' -#' @param reconstruction Result list of the dot pattern reconstruction with +#' @param x rd_multi Object created with \code{reconstruct_pattern_multi}. #' multiple marks. +#' @param verbose Logical if progress should be printed. +#' @param ... Currently not used #' #' @details #' Calculates and visualises various summary statistics for the results of #' multi-marks point pattern reconstruction. #' -#' @return ggplot object +#' @seealso +#' \code{\link{reconstruct_pattern_multi}} #' -#' @aliases plot_sum_stat -#' @rdname plot_sum_stat +#' @return void +#' +#' @aliases plot.rd_multi +#' @rdname plot.rd_multi #' #' @export -plot_sum_stat <- function(reconstruction) { +plot.rd_multi <- function(x, verbose = TRUE, ...) { # Data import from the results of point pattern reconstruction. k_func_recon <- list() @@ -27,34 +32,34 @@ plot_sum_stat <- function(reconstruction) { divisor <- "r" kernel_arg <- "epanechnikov" - if(nchar(names(reconstruction[1])) == 9){ + if(nchar(names(x[1])) == 9){ n_repetitions <- 1 - rmax <- as.numeric(reconstruction$Parameter_setting$rmax) - rcount <- as.numeric(reconstruction$Parameter_setting$rcount) - rpresented <- min(c(as.numeric(reconstruction$window[2]),as.numeric(reconstruction$window[4]))) + rmax <- as.numeric(x$Parameter_setting$rmax) + rcount <- as.numeric(x$Parameter_setting$rcount) + rpresented <- min(c(as.numeric(x$window[2]),as.numeric(x$window[4]))) r <- seq(rmin, if(rpresented >= rmax*2 ){rmax*2}else{rpresented}, length.out = rcount) - bw <- as.numeric(reconstruction$Parameter_setting$bw) + bw <- as.numeric(x$Parameter_setting$bw) - ppp_reference <- spatstat.geom::as.ppp(reconstruction$reference, reconstruction$window) - reconstruction <- list(reconstruction) + ppp_reference <- spatstat.geom::as.ppp(x$reference, x$window) + x <- list(x) } else { - n_repetitions <- reconstruction[[1]]$Parameter_setting$n_repetitions - rmax <- as.numeric(reconstruction[[n_repetitions]]$Parameter_setting$rmax) - rpresented <- min(c(as.numeric(reconstruction[[n_repetitions]]$window[2]),as.numeric(reconstruction[[n_repetitions]]$window[4]))) - rcount <- as.numeric(reconstruction[[n_repetitions]]$Parameter_setting$rcount) + n_repetitions <- x[[1]]$Parameter_setting$n_repetitions + rmax <- as.numeric(x[[n_repetitions]]$Parameter_setting$rmax) + rpresented <- min(c(as.numeric(x[[n_repetitions]]$window[2]),as.numeric(x[[n_repetitions]]$window[4]))) + rcount <- as.numeric(x[[n_repetitions]]$Parameter_setting$rcount) r <- seq(rmin, if(rpresented <= rmax*2 ){rmax*2}else{rpresented}, length.out = rcount) - bw <- as.numeric(reconstruction[[n_repetitions]]$Parameter_setting$bw) + bw <- as.numeric(x[[n_repetitions]]$Parameter_setting$bw) - ppp_reference <- spatstat.geom::as.ppp(reconstruction$reconstruction_1$reference, reconstruction$reconstruction_1$window) + ppp_reference <- spatstat.geom::as.ppp(x$reconstruction_1$reference, x$reconstruction_1$window) } for (i in seq_len(n_repetitions)) { - message("Progress in the creation of the figures: ", i/n_repetitions*100,"%\t\t\r", appendLF = FALSE) + if (verbose) message("Progress in the creation of the figures: ", round(i/n_repetitions*100),"%\t\t\r", appendLF = FALSE) - ppp_reconstructed <- spatstat.geom::as.ppp(reconstruction[[i]][[2]], reconstruction[[i]][[3]]) + ppp_reconstructed <- spatstat.geom::as.ppp(x[[i]][[2]], x[[i]][[3]]) pcf_func <-spatstat.explore::pcf.ppp(ppp_reconstructed, bw = bw, kernel=kernel_arg, correction = "none", divisor = divisor, r = r) @@ -168,4 +173,6 @@ plot_sum_stat <- function(reconstruction) { graphics::par(mfrow = c(1, 1)) + invisible() + } diff --git a/R/reconstruct_pattern_multi.R b/R/reconstruct_pattern_multi.R index 72849ca9..ba9b9252 100644 --- a/R/reconstruct_pattern_multi.R +++ b/R/reconstruct_pattern_multi.R @@ -70,7 +70,7 @@ #' \code{\link{reconstruct_pattern}} \cr #' \code{\link{reconstruct_pattern_marks}} #' -#' @return list +#' @return rd_multi #' #' @examples #' \dontrun{ @@ -528,8 +528,10 @@ reconstruct_pattern_multi <- function(marked_pattern, } if(n_repetitions > 1) { + class(reconstruction_list) <- "rd_multi" return(reconstruction_list) } else { + class(reconstruction) <- "rd_multi" return(reconstruction) } } diff --git a/man/plot.rd_multi.Rd b/man/plot.rd_multi.Rd new file mode 100644 index 00000000..d6c46f4c --- /dev/null +++ b/man/plot.rd_multi.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/plot.rd_multi.R +\name{plot.rd_multi} +\alias{plot.rd_multi} +\title{plot_sum_stat} +\usage{ +\method{plot}{rd_multi}(x, verbose = TRUE, ...) +} +\arguments{ +\item{x}{rd_multi Object created with \code{reconstruct_pattern_multi}. +multiple marks.} + +\item{verbose}{Logical if progress should be printed.} + +\item{...}{Currently not used} +} +\value{ +void +} +\description{ +plot_sum_stat +} +\details{ +Calculates and visualises various summary statistics for the results of +multi-marks point pattern reconstruction. +} +\seealso{ +\code{\link{reconstruct_pattern_multi}} +} diff --git a/man/plot_sum_stat.Rd b/man/plot_sum_stat.Rd deleted file mode 100644 index 44db8913..00000000 --- a/man/plot_sum_stat.Rd +++ /dev/null @@ -1,22 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/plot_sum_stat.R -\name{plot_sum_stat} -\alias{plot_sum_stat} -\title{plot_sum_stat} -\usage{ -plot_sum_stat(reconstruction) -} -\arguments{ -\item{reconstruction}{Result list of the dot pattern reconstruction with -multiple marks.} -} -\value{ -ggplot object -} -\description{ -plot_sum_stat -} -\details{ -Calculates and visualises various summary statistics for the results of -multi-marks point pattern reconstruction. -} diff --git a/man/reconstruct_pattern_multi.Rd b/man/reconstruct_pattern_multi.Rd index dd652c5c..28941ac2 100644 --- a/man/reconstruct_pattern_multi.Rd +++ b/man/reconstruct_pattern_multi.Rd @@ -98,7 +98,7 @@ include Dk, K, Hs, pcf.} \item{verbose}{Logical if progress report is printed.} } \value{ -list +rd_multi } \description{ Pattern reconstruction of a pattern marked by multiple traits. diff --git a/tests/testthat/test-plot_rd_multi.R b/tests/testthat/test-plot_rd_multi.R new file mode 100644 index 00000000..e724155e --- /dev/null +++ b/tests/testthat/test-plot_rd_multi.R @@ -0,0 +1,34 @@ +# testthat::context("test-plot_rd_multi") + +# create random data +xr <- 500 +yr <- 1000 + +N <- 400 + +y <- runif(N, min = 0, max = yr) +x <- runif(N, min = 0, max = xr) + +species <- sample(c("A","B"), N, replace = TRUE) +diameter <- runif(N, 0.1, 0.4) + +random <- data.frame(x = x, y = y, dbh = diameter, species = factor(species)) + +# Conversion to a ppp object and conversion of the metric mark to metres. +marked_pattern <- spatstat.geom::as.ppp(random, W = spatstat.geom::owin(c(0, 500), c(0, 1000))) + +# Reconstruction function +multi_recon <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 3, max_steps = 10, + verbose = FALSE) + +multi_recon_simple <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10, + verbose = FALSE) + +################################################################################ + +testthat::test_that("plot returns plot", { + + testthat::expect_null(plot(multi_recon, verbose = FALSE)) + testthat::expect_null(plot(multi_recon_simple, verbose = FALSE)) + +}) diff --git a/vignettes/articles/reconstruct_multi_patterns.Rmd b/vignettes/articles/reconstruct_multi_patterns.Rmd index 670c9cf2..42d5c953 100644 --- a/vignettes/articles/reconstruct_multi_patterns.Rmd +++ b/vignettes/articles/reconstruct_multi_patterns.Rmd @@ -84,5 +84,5 @@ of the reference pattern (black line) compared to the reconstructed pattern additionally. ```{r} -plot_sum_stat(reconstruction) +plot(reconstruction) ``` From cf387aacaf0d41c9e74011551c9f62246148a380 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 15:17:33 +0100 Subject: [PATCH 34/39] [skip-ci] Fix homepage --- _pkgdown.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index 80eb5fed..a48a63fe 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -31,6 +31,7 @@ reference: - plot.rd_mar - plot.rd_pat - plot.rd_ras + - plot.rd_multi - print.rd_mar - print.rd_pat - print.rd_ras From 05fbffcbfc4d7032223bfdb112b86487f9456fd0 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Mon, 18 Dec 2023 15:24:21 +0100 Subject: [PATCH 35/39] As always forgot some docs :) --- R/plot.rd_multi.R | 4 ++-- man/plot.rd_multi.Rd | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/R/plot.rd_multi.R b/R/plot.rd_multi.R index 8fef96a0..a4b65a73 100644 --- a/R/plot.rd_multi.R +++ b/R/plot.rd_multi.R @@ -1,6 +1,6 @@ -#' plot_sum_stat +#' plot.rd_multi #' -#' @description plot_sum_stat +#' @description Plot method for rd_multi object #' #' @param x rd_multi Object created with \code{reconstruct_pattern_multi}. #' multiple marks. diff --git a/man/plot.rd_multi.Rd b/man/plot.rd_multi.Rd index d6c46f4c..3c6cbe7a 100644 --- a/man/plot.rd_multi.Rd +++ b/man/plot.rd_multi.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/plot.rd_multi.R \name{plot.rd_multi} \alias{plot.rd_multi} -\title{plot_sum_stat} +\title{plot.rd_multi} \usage{ \method{plot}{rd_multi}(x, verbose = TRUE, ...) } @@ -18,7 +18,7 @@ multiple marks.} void } \description{ -plot_sum_stat +Plot method for rd_multi object } \details{ Calculates and visualises various summary statistics for the results of From 11b3256cbc3eb925c7902f4f54e4df311d7122b2 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Tue, 19 Dec 2023 08:47:28 +0100 Subject: [PATCH 36/39] [skip-ci] Some minor cosmetic changes to vignettes --- .../articles/reconstruct_multi_patterns.Rmd | 49 ++++++------------- .../articles/reconstruct_several_patterns.Rmd | 10 ++-- 2 files changed, 20 insertions(+), 39 deletions(-) diff --git a/vignettes/articles/reconstruct_multi_patterns.Rmd b/vignettes/articles/reconstruct_multi_patterns.Rmd index 42d5c953..73c8e451 100644 --- a/vignettes/articles/reconstruct_multi_patterns.Rmd +++ b/vignettes/articles/reconstruct_multi_patterns.Rmd @@ -16,16 +16,16 @@ editor_options: If you want more information about multi-trait point pattern reconstruction, please refer to the [corresponding paper](https://besjournals.onlinelibrary.wiley.com/doi/10.1111/2041-210X.14206 ). -If you wish to include several marks simultaneously in a reconstruction, you can -use the following code. The libraries used must first be loaded. +If you wish to include several marks simultaneously in a reconstruction, you can use the following code. The libraries used must first be loaded. ```{r message = FALSE, warning = FALSE} library(shar) library(spatstat) ``` -The next step is to load the point pattern, here is an example of a random point -pattern with several marks to show the structure of the data used. +> Please note that the maximum number of iterations has been set to `max_steps = 10000` and `n_repetitions = 1`/`n_repetitions = 3` to keep computational time low for this example. For real-world applications, it is advisable to raise these values. Additionally, we set `verbose = FALSE` in the vignette to minimize printed output. We recommend using the default setting `verbose = TRUE` when executing the code to view progress reports. + +The next step is to load the point pattern, here is an example of a random point pattern with several marks to show the structure of the data used. ```{r} xr <- 500 @@ -39,49 +39,28 @@ random <- data.frame(x = x, y = y, dbh = diameter, species = factor(species)) marked_pattern <- as.ppp(random, W = owin(c(0, xr), c(0, yr))) ``` -The point pattern must contain the following data An x and y coordinate, a -metric mark (in metres) and a nominal mark defined as a factor. The order must -be respected. Now the reconstruction with several marks can be started with the -following code. Note that the maximum number of iterations has been set to -max_steps = 10000 to keep the computation time for this example to a minimum. -For an application, this value should be increased according to the number of -points in the pattern. +The point pattern must contain the following data An x and y coordinate, a metric mark (in metres) and a nominal mark defined as a factor. The order must be respected. Now the reconstruction with several marks can be started with the following code. Note that the maximum number of iterations has been set to max_steps = 10000 to keep the computation time for this example to a minimum. For an application, this value should be increased according to the number of points in the pattern. ```{r} -reconstruction <- reconstruct_pattern_multi( -marked_pattern, -n_repetitions = 1, -max_steps = 10000) +reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000, + verbose = FALSE) ``` -As a result, you will receive a list containing a variety of information, for -example, the reference pattern, the reconstructed pattern, the number of -successful actions, the energy development and much more. If you wish to -perform several reconstructions of the same reference pattern, you must -increase n_repetitions to the desired number. +As a result, you will receive a list containing a variety of information, for example, the reference pattern, the reconstructed pattern, the number of successful actions, the energy development and much more. If you wish to perform several reconstructions of the same reference pattern, you must increase n_repetitions to the desired number. ```{r} -reconstruction_2 <- reconstruct_pattern_multi( -marked_pattern, -n_repetitions = 2, -max_steps = 10000) +reconstruction_2 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 3, max_steps = 10000, + verbose = FALSE) ``` -To activate a visualisation of the reconstruction that shows the changes in the -pattern at the relevant time, you must proceed as follows. +To activate a visualisation of the reconstruction that shows the changes in the pattern at the relevant time, you must proceed as follows. ```{r} -reconstruction_3 <- reconstruct_pattern_multi( -marked_pattern, -n_repetitions = 1, -max_steps = 10000, -plot = TRUE) +reconstruction_3 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000, plot = TRUE, + verbose = FALSE) ``` -Finally, you can use the following function to view different summary statistics -of the reference pattern (black line) compared to the reconstructed pattern -(grey line). For this, however, the listed libraries must be loaded -additionally. +Finally, you can use the following function to view different summary statistics of the reference pattern (black line) compared to the reconstructed pattern (grey line). For this, however, the listed libraries must be loaded additionally. ```{r} plot(reconstruction) diff --git a/vignettes/articles/reconstruct_several_patterns.Rmd b/vignettes/articles/reconstruct_several_patterns.Rmd index a23a978e..92d4e406 100644 --- a/vignettes/articles/reconstruct_several_patterns.Rmd +++ b/vignettes/articles/reconstruct_several_patterns.Rmd @@ -18,7 +18,9 @@ knitr::opts_chunk$set( ) ``` -In case you want to reconstruct several patterns at once (e.g. for different points in time if repeated censuses are available), you can use the following code. Please be aware that the maximum number of iterations was set to `max_runs = 10` to keep the computational time low for this example. For an applied case, this value should be increased. +In case you want to reconstruct several patterns at once (e.g. for different points in time if repeated censuses are available), you can use the following code. + +> Please note that the maximum number of iterations has been set to `max_runs = 1000` and `n_random = 3` to keep computational time low for this example. For real-world applications, it is advisable to raise these values. Additionally, we set `verbose = FALSE` in the vignette to minimize printed output. We recommend using the default setting `verbose = TRUE` when executing the code to view progress reports. ```{r load-packages, message = FALSE, warning = FALSE} library(shar) @@ -34,7 +36,7 @@ list_pattern <- list(species_a, species_b) # reconstruct all patterns in list result <- lapply(list_pattern, function(x) reconstruct_pattern(pattern = x, n_random = 3, - max_runs = 10, verbose = FALSE)) + max_runs = 1000, verbose = FALSE)) ``` @@ -53,7 +55,7 @@ Firstly, reconstruct only the spatial characteristics *n* times. The observed pa ```{r reconstruct-pattern} # reconstruct spatial strucutre reconstructed_pattern <- reconstruct_pattern(species_a, n_random = 3, - max_runs = 10, return_input = FALSE, + max_runs = 1000, return_input = FALSE, verbose = FALSE) ``` @@ -67,7 +69,7 @@ species_a_marks <- subset(species_a, select = dbh) result_marks <- lapply(reconstructed_pattern$randomized, function(x) reconstruct_pattern_marks(pattern = x, marked_pattern = species_a_marks, - max_runs = 10, + max_runs = 1000, n_random = 3, verbose = FALSE)) ``` From fa9fbe9b5b1e663e7d014dd09e3f97f96dc5e0d9 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Tue, 19 Dec 2023 14:01:14 +0100 Subject: [PATCH 37/39] Better legend placement --- R/plot.rd_multi.R | 6 +++--- vignettes/articles/reconstruct_multi_patterns.Rmd | 13 +++++-------- 2 files changed, 8 insertions(+), 11 deletions(-) diff --git a/R/plot.rd_multi.R b/R/plot.rd_multi.R index a4b65a73..d69ae098 100644 --- a/R/plot.rd_multi.R +++ b/R/plot.rd_multi.R @@ -127,7 +127,7 @@ plot.rd_multi <- function(x, verbose = TRUE, ...) { colnames(species_markcorr_all)[1] <- c("Reference_mark_species") colnames(species_markcorr_all)[length(species_markcorr_all)] <- c("r") - graphics::par(mfrow = c(2, 2)) + graphics::par(mfrow = c(2, 2), mar = c(4.0, 4.0, 1.5, 0.5)) # bltr # plot Kest graphics::plot(NULL, xlim = range(r), ylim = range(c(k_func_all$Reference, @@ -139,7 +139,7 @@ plot.rd_multi <- function(x, verbose = TRUE, ...) { graphics::abline(v = rmax, lty = "dotted", col = "grey") graphics::legend(x = "topleft", legend = c("Reference", "Reconstructed"), - lty = c(2, 1), inset = 0.015) + lty = c(2, 1), bty = "n", cex = 0.8, ncol = 2) # plot pcf graphics::plot(NULL, xlim = range(r), ylim = range(c(pcf_func_all$Reference, @@ -171,7 +171,7 @@ plot.rd_multi <- function(x, verbose = TRUE, ...) { graphics::lines(x = species_markcorr_all$r, y = species_markcorr_all$species_markcorr_func_recon_mean) graphics::abline(v = rmax, lty = "dotted", col = "grey") - graphics::par(mfrow = c(1, 1)) + graphics::par(mfrow = c(1, 1), c(5.1, 4.1, 4.1, 2.1)) invisible() diff --git a/vignettes/articles/reconstruct_multi_patterns.Rmd b/vignettes/articles/reconstruct_multi_patterns.Rmd index 73c8e451..f7ce7bef 100644 --- a/vignettes/articles/reconstruct_multi_patterns.Rmd +++ b/vignettes/articles/reconstruct_multi_patterns.Rmd @@ -23,7 +23,7 @@ library(shar) library(spatstat) ``` -> Please note that the maximum number of iterations has been set to `max_steps = 10000` and `n_repetitions = 1`/`n_repetitions = 3` to keep computational time low for this example. For real-world applications, it is advisable to raise these values. Additionally, we set `verbose = FALSE` in the vignette to minimize printed output. We recommend using the default setting `verbose = TRUE` when executing the code to view progress reports. +> Please note that the maximum number of iterations has been set to `max_steps = 10000` and `n_repetitions = 1`/`n_repetitions = 3` to keep computational time low for this example. For real-world applications, it is advisable to raise these values. The next step is to load the point pattern, here is an example of a random point pattern with several marks to show the structure of the data used. @@ -42,25 +42,22 @@ marked_pattern <- as.ppp(random, W = owin(c(0, xr), c(0, yr))) The point pattern must contain the following data An x and y coordinate, a metric mark (in metres) and a nominal mark defined as a factor. The order must be respected. Now the reconstruction with several marks can be started with the following code. Note that the maximum number of iterations has been set to max_steps = 10000 to keep the computation time for this example to a minimum. For an application, this value should be increased according to the number of points in the pattern. ```{r} -reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000, - verbose = FALSE) +reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000) ``` As a result, you will receive a list containing a variety of information, for example, the reference pattern, the reconstructed pattern, the number of successful actions, the energy development and much more. If you wish to perform several reconstructions of the same reference pattern, you must increase n_repetitions to the desired number. ```{r} -reconstruction_2 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 3, max_steps = 10000, - verbose = FALSE) +reconstruction_2 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 3, max_steps = 10000) ``` To activate a visualisation of the reconstruction that shows the changes in the pattern at the relevant time, you must proceed as follows. ```{r} -reconstruction_3 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000, plot = TRUE, - verbose = FALSE) +reconstruction_3 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000, plot = TRUE) ``` -Finally, you can use the following function to view different summary statistics of the reference pattern (black line) compared to the reconstructed pattern (grey line). For this, however, the listed libraries must be loaded additionally. +Finally, you can use the `plot()` function to view different summary statistics of the reference pattern (dashed line) compared to the reconstructed pattern (solid line). ```{r} plot(reconstruction) From 5b03d68601d99163e14dfa9aeb709d54a7d1a00c Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Tue, 19 Dec 2023 14:03:22 +0100 Subject: [PATCH 38/39] Upsidaisy..forgot mar argument --- R/plot.rd_multi.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/plot.rd_multi.R b/R/plot.rd_multi.R index d69ae098..41ef11c2 100644 --- a/R/plot.rd_multi.R +++ b/R/plot.rd_multi.R @@ -171,7 +171,7 @@ plot.rd_multi <- function(x, verbose = TRUE, ...) { graphics::lines(x = species_markcorr_all$r, y = species_markcorr_all$species_markcorr_func_recon_mean) graphics::abline(v = rmax, lty = "dotted", col = "grey") - graphics::par(mfrow = c(1, 1), c(5.1, 4.1, 4.1, 2.1)) + graphics::par(mfrow = c(1, 1), mar = c(5.1, 4.1, 4.1, 2.1)) invisible() From 8b37fdebc6b1f2118e4a9d6ef3bcb46aece34894 Mon Sep 17 00:00:00 2001 From: mhesselbarth Date: Thu, 21 Dec 2023 09:45:12 +0100 Subject: [PATCH 39/39] Last fix vignettes --- vignettes/articles/reconstruct_multi_patterns.Rmd | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/vignettes/articles/reconstruct_multi_patterns.Rmd b/vignettes/articles/reconstruct_multi_patterns.Rmd index f7ce7bef..53a92138 100644 --- a/vignettes/articles/reconstruct_multi_patterns.Rmd +++ b/vignettes/articles/reconstruct_multi_patterns.Rmd @@ -42,23 +42,23 @@ marked_pattern <- as.ppp(random, W = owin(c(0, xr), c(0, yr))) The point pattern must contain the following data An x and y coordinate, a metric mark (in metres) and a nominal mark defined as a factor. The order must be respected. Now the reconstruction with several marks can be started with the following code. Note that the maximum number of iterations has been set to max_steps = 10000 to keep the computation time for this example to a minimum. For an application, this value should be increased according to the number of points in the pattern. ```{r} -reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000) +reconstruction <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000, issue = 5000) ``` As a result, you will receive a list containing a variety of information, for example, the reference pattern, the reconstructed pattern, the number of successful actions, the energy development and much more. If you wish to perform several reconstructions of the same reference pattern, you must increase n_repetitions to the desired number. ```{r} -reconstruction_2 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 3, max_steps = 10000) +reconstruction_2 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 3, max_steps = 10000, issue = 5000) ``` To activate a visualisation of the reconstruction that shows the changes in the pattern at the relevant time, you must proceed as follows. ```{r} -reconstruction_3 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000, plot = TRUE) +reconstruction_3 <- reconstruct_pattern_multi(marked_pattern, n_repetitions = 1, max_steps = 10000, issue = 5000, plot = TRUE) ``` Finally, you can use the `plot()` function to view different summary statistics of the reference pattern (dashed line) compared to the reconstructed pattern (solid line). ```{r} -plot(reconstruction) +plot(reconstruction, verbose = FALSE) ```