diff --git a/.gitignore b/.gitignore index 87b9636a..b7327483 100644 --- a/.gitignore +++ b/.gitignore @@ -25,6 +25,11 @@ vignettes/*.html vignettes/*.pdf vignettes/*_files/ +# Rendered output from the data-raw notebooks (sources are tracked, output +# is regenerated by knitting them) +data-raw/*.html +data-raw/*_files/ + # OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3 .httr-oauth diff --git a/DESCRIPTION b/DESCRIPTION index a38914a5..1b77684f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BirdFlowR Title: Predict and Visualize Bird Movement -Version: 0.1.0.9083 +Version: 0.1.0.9084 Authors@R: c(person("Ethan", "Plunkett", email = "plunkett@umass.edu", role = c("aut", "cre"), comment = c(ORCID = "0000-0003-4405-2251")), diff --git a/NAMESPACE b/NAMESPACE index a6383343..330fe9d2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -122,6 +122,7 @@ export(validate_BirdFlow) export(validate_BirdFlowIntervals) export(validate_BirdFlowRoutes) export(validate_Routes) +export(visualize_distance_weights) export(write_intervals) export(write_routes) export(x_to_col) diff --git a/NEWS.md b/NEWS.md index 73a29002..1d1f05ca 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,68 @@ +# BirdFlowR 0.1.0.9084 +2026-07-27 + +## Continuous BMTR + +The focus of this update is calculating Bird Flow Migration Traffic (BMTR) +with a detection rate that varies continuously, as opposed to the default +method which has a binary detection rate. + +Previously `calc_bmtr()` did have both methods but the continuous method +was prohibitively computationally expensive. This update makes it usable +with non-trivial models. + +**Breaking change** This update changes the arguments used to call the +continuous version of the function dropping `weighted` and `euclidean` +arguments. I believe no users are affected by this change, but +I apologize if you were. + +## Specific changes +* `calc_bmtr()`: replaced the `weighted`/`euclidean` arguments with a single + `method` argument (`"binary"`, `"continuous"`, or `"continuous-spherical"`) + that selects the detection model used to assess whether a transition + passes by a point. `"continuous"` uses fast planar (Euclidean) geometry + and is the recommended weighted option; `"continuous-spherical"` is the + slower great-circle version, kept to assess the impact of switching + geometries. +* `calc_bmtr()`: new `...` argument, forwarded to `calc_dist_weights()` for + the two continuous methods, to tune the spread kernel (`kernel`, `gamma`, + `kl`, `s1`) used to model uncertainty in a bird's path. +* `calc_dist_weights()`: renamed the `method` argument to `kernel` to avoid + colliding with `calc_bmtr()`'s new `method` argument. +* Fixed a bug where `calc_bmtr()`'s weighted detection methods ignored a + user-supplied `points` argument due to a positional argument-binding bug. +* Fixed a bug where `calc_bmtr()`'s weighted detection methods ignored a + user-supplied `radius` argument, always using `mean(res(bf))` instead. +* Replaced `weight_between()` and the internal `calc_detection_rate()` + dispatcher with `calc_spherical_detection_rate()` and + `calc_euclidean_detection_rate()`. +* `calc_dist_weights()`: `kernel` now supports 3 additional spread kernels + (`"m1"`, `"m5"`, `"sq"`) alongside the existing `"m3"` and `"bb"`. +* New `visualize_distance_weights()`: plots the shape of `calc_dist_weights()`'s + spread kernel, as an SD envelope or a weight raster, to help choose + `kernel` and its hyperparameters before running `calc_bmtr()`. +* Fixed a units bug in `calc_dist_weights()`'s `"m1"`, `"m3"`, `"m5"`, and + `"sq"` kernels: `dist_along_line`/`line_lengths`/`res_m` were converted to + km while `gamma`/`kl` were left in meters, making the kernel's contribution + to the spread negligible next to the fixed nugget at any realistic line + length. Everything is now kept in meters throughout, matching `gamma`/`kl`. + This changes the numeric output of `calc_bmtr()`'s `"continuous"`/ + `"continuous-spherical"` methods and of `calc_euclidean_detection_rate()`/ + `calc_spherical_detection_rate()` for all kernels, including `"bb"` (whose + `t`/`len` inputs are now meters rather than km, so its `s1` hyperparameter + means something different than before). +* `calc_dist_weights()`/`visualize_distance_weights()`: `gamma`/`kl`/`s1` + now default to `3e10`/`9e5`/`200` instead of the old fixed literals + (`40000`/`2000`/`10`) tuned for a ~1-2 km demo scale. The three + hyperparameters have different units (`gamma` is a variance, `kl` a + length, `s1` a sqrt-length). The new values were chosen by eye against + `visualize_distance_weights()`'s envelope plots. + See `data-raw/kernel_parameter_tuning.Rmd`. + This changes the default numeric output of `calc_bmtr()`'s `"continuous"`/ + `"continuous-spherical"` and of + `calc_euclidean_detection_rate()`/`calc_spherical_detection_rate()`. + # BirdFlowR 0.1.0.9083 2026-06-15 diff --git a/R/calc_bmtr.R b/R/calc_bmtr.R index eb1b33ea..a47608ce 100644 --- a/R/calc_bmtr.R +++ b/R/calc_bmtr.R @@ -18,10 +18,14 @@ #' #' @section Limitations: #' -#' `calc_bmtr()` makes the incorrect simplifying assumption -#' that birds follow the shortest (great circle) path -#' between the center of the the source and destination raster cells. Caution -#' should be used when interpreting the results especially around +#' `calc_bmtr()` makes an incorrect simplifying assumption about the path +#' birds take: either that they follow the shortest path between the centers +#' of the source and destination raster cells (`method = "binary"`), or that +#' they are distributed symmetrically in a Gaussian distribution around that +#' path (the two continuous methods). That path is a great circle in all +#' cases except `method = "continuous"`, which uses a straight line in the +#' model's projected CRS. +#' Caution should be used when interpreting the results especially around #' major geographic features such as coasts, large lakes, mountain ranges, and #' ecological system boundaries that might result in non-linear migration paths. #' @@ -41,16 +45,19 @@ #' standardizing the units based on the entire cell area that that point #' represents. #' -#' #' @param bf A BirdFlow model #' @param points A set of points to calculate movement through. If `points` is #' `NULL` they will default to the BirdFlow model cells that are either active #' or fall between two active cells. Otherwise a data frame with `x` and `y` #' columns containing point coordinates in [crs(bf)][terra::crs()]. -#' @param radius The radius in meters around the points used to assess whether -#' a movement line passes by (or through) the point. If a point is farther than -#' `radius` from a great circle line between two cells centers then it is not -#' between them. +#' @param radius The radius in meters around the points used to assess the +#' detection rate for a movement at the point. +#' With `method = "binary"`, if a point is within `radius` of the great +#' circle line between two cell centers then the movement is detected at +#' that point. For the two continuous detection methods there is a +#' probability distribution for the location of the bird as it passes by the +#' point, and the radius defines the band over which that probability is +#' integrated, giving the detection rate. #' @param n_directions The number of directional bins to use for recording #' movement direction. Must be either `1` indicating no direction information #' or an even number. This is a placeholder, currently only `1` is supported. @@ -71,32 +78,71 @@ #' transition.} #'} #' @inheritParams is_between -#' @param weighted If `FALSE` use the original and quicker version of bmtr -#' that sums all the marginal probability for transitions that pass within a -#' fixed distance of the point. If `TRUE` assign a weight to the point and -#' transition combo that then is multiplied by the marginal probability before -#' summing. This argument is experimental but the default value is identical -#' to the old version. The argument name and behavior when set to `TRUE` may -#' change. +#' @param method The detection model used to determine how much of a +#' transition's movement counts towards a point's BMTR: +#' \describe{ +#' \item{`"binary"`}{(default) Fast and deterministic. A movement line either +#' does or does not pass within `radius` of the point, per [is_between()].} +#' \item{`"continuous"`}{Assigns a continuous weight (0 to 1) based on the +#' probability that a bird's actual path, modeled as spreading away from the +#' straight line between two cells, passes within `radius` of the point. Uses +#' planar (Euclidean) geometry in the model's native CRS, and is the +#' recommended detection model when continuous weighting is desired.} +#' \item{`"continuous-spherical"`}{The same continuous weighting as +#' `"continuous"`, but computed with great-circle (spherical) geometry +#' instead. Much slower, and not recommended for routine use; kept to allow +#' assessing the impact of switching from spherical to Euclidean geometry.} +#' } +#' @param ... For `method = "continuous"` or `"continuous-spherical"`, +#' additional arguments forwarded to [calc_dist_weights()] to control the +#' spread kernel used to model uncertainty in a bird's path: `kernel`, +#' `gamma`, `kl`, and `s1`. See [calc_dist_weights()] for the full list of +#' supported kernels and their hyperparameters, and +#' [visualize_distance_weights()] to explore how they shape the spread +#' before running this (potentially expensive) function. Not applicable, +#' and an error, for `method = "binary"`. #' #' @return See `format` argument. +#' @seealso [visualize_distance_weights()] to explore the spread kernel +#' hyperparameters accepted via `...`. #' @export #' #' @examples #' #' \dontrun{ #' bf <- BirdFlowModels::amewoo -#' bmtr <- calc_bmtr(bf) #' +#' # Binary detection along shortest path +#' bmtr <- calc_bmtr(bf) #' plot_bmtr(bmtr, bf) #' #' animate_bmtr(bmtr, bf) +#' +#' # Continuous detection with a wider spread kernel +#' # (gamma defaults to 3e10) +#' bmtr2 <- calc_bmtr(bf, method = "continuous", gamma = 6e10) +#' animate_bmtr(bmtr2, bf) +#' +#' # Visualize spread for the continuous detection with wider spread +#' visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf), +#' gamma = 6e10, res_m = xres(bf)) +#' +#' # Visualize spread for the continuous detection with default values +#' visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf), +#' res_m = xres(bf)) +#' #' } #' +#' +#' +#' calc_bmtr <- function(bf, points = NULL, radius = NULL, n_directions = 1, format = NULL, batch_size = 5e5, check_radius = TRUE, - weighted = FALSE) { + method = c("binary", "continuous", + "continuous-spherical"), + ...) { + method <- match.arg(method) if (!requireNamespace("SparseArray", quietly = TRUE)) { stop("The SparseArray package is required to use calc_bmtr(). ", @@ -117,15 +163,23 @@ calc_bmtr <- function(bf, points = NULL, radius = NULL, n_directions = 1, format <- tolower(format) stopifnot(format %in% c("points", "spatraster", "dataframe")) - # The only difference between is_between() and weight_between() return formats - # is in the "between" component. In the first it's logical (TRUE is between) - # in the second it is a weight between 0 and 1 (non zero indicates some - # level of betweeness". - if (weighted) { - result <- weight_between(bf, points, radius, n_directions) - } else { - result <- is_between(bf, points, radius, n_directions) - } + # The only difference between is_between() and the continuous detection + # functions' return formats is in the "between" component. In the first + # it's logical (TRUE is between) in the second it is a weight between 0 + # and 1 (non zero indicates some level of betweenness). + result <- switch( + method, + binary = is_between(bf, points = points, radius = radius, + n_directions = n_directions, + batch_size = batch_size, check_radius = check_radius, + ...), + continuous = calc_euclidean_detection_rate( + bf, points = points, radius = radius, n_directions = n_directions, + batch_size = batch_size, check_radius = check_radius, ...), + "continuous-spherical" = calc_spherical_detection_rate( + bf, points = points, radius = radius, n_directions = n_directions, + batch_size = batch_size, check_radius = check_radius, ...) + ) between <- result$between points <- result$points diff --git a/R/calc_distance_weights.R b/R/calc_distance_weights.R index f7242657..fc666bbf 100644 --- a/R/calc_distance_weights.R +++ b/R/calc_distance_weights.R @@ -24,62 +24,92 @@ #' determine the nugget added to the variance to represent the uncertainty in #' the starting and ending location of the transition. #' -#' @param method The method used for calculating the standard deviation -#' in the probability distribution. Currently `"m3"`, Martern 3/2; and -#' `"bb"`, Brownian bridge are supported. +#' @param kernel The kernel used for calculating the standard deviation +#' in the probability distribution: +#' \describe{ +#' \item{`"m3"`}{(default) Matern 3/2. Driven by `gamma` and `kl`.} +#' \item{`"m1"`}{Matern 1/2. Driven by `gamma` and `kl`.} +#' \item{`"m5"`}{Matern 5/2. Driven by `gamma` and `kl`.} +#' \item{`"sq"`}{Squared-exponential (Gaussian). Driven by `gamma` and `kl`.} +#' \item{`"bb"`}{Brownian bridge. Driven by `s1`.} +#' } +#' See [visualize_distance_weights()] to explore how these kernels and their +#' hyperparameters shape the spread. +#' @param gamma Spread magnitude hyperparameter (m^2, a variance) for the +#' Matern-family and squared-exponential kernels (`"m1"`, `"m3"`, `"m5"`, +#' `"sq"`). Ignored for `kernel = "bb"`. The default was tuned by eye with +#' `visualize_distance_weights()`. +#' @param kl Lengthscale hyperparameter (m) for the Matern-family and +#' squared-exponential kernels (`"m1"`, `"m3"`, `"m5"`, `"sq"`). Ignored for +#' `kernel = "bb"`. The default was tuned by eye using +#' `visualize_distance_weights()`. +#' @param s1 Spread magnitude hyperparameter (units of sqrt(m), not m) for +#' the `"bb"` kernel. Ignored for all other kernels. The default was tuned by +#' eye against `visualize_distance_weights()`. #' #' @return A vector of weights of the same length as the first three arguments. +#' @seealso [visualize_distance_weights()] #' @keywords internal calc_dist_weights <- function(dist_to_line, dist_along_line, line_lengths, - radius_m, res_m, method = "m3") { + radius_m, res_m, + kernel = c("m3", "bb", "m1", "m5", "sq"), + gamma = 3e10, kl = 9e5, s1 = 200) { + kernel <- match.arg(kernel) + sd <- calc_dist_weights_sd(dist_along_line, line_lengths, res_m, + kernel = kernel, gamma = gamma, kl = kl, s1 = s1) - # Spatial information - len <- line_lengths / 1000 # length of great circle in km - t <- dist_along_line / 1000 # Distance from start of great circle to where - # the point projects onto the great circle in km - d <- dist_to_line / 1000 # Distance from reference point to the great circle - - res_km <- res_m / 1000 - - r <- radius_m / 1000 # Radius of transact in KM. - - # s2 is the standard deviation of the nugget - s2 <- res_km / 4 # This initial setting of 1/4 the cell width converted - # to KM means the cell boundary (orthagonally) is at 2 - # standard deviations from the cell center, the corners - # will be at 2.8 SD. - - valid_methods <- c("bb", # Brownian Bridge - "m3") # Martern 3/2 + # Calculate weight + weight <- rep(0, length(sd)) + in_range <- (dist_to_line - radius_m) < 1.96 * sd + low_cum_prob <- stats::pnorm(dist_to_line[in_range] - radius_m, + sd = sd[in_range]) + high_cum_prob <- stats::pnorm(dist_to_line[in_range] + radius_m, + sd = sd[in_range]) + weight[in_range] <- high_cum_prob - low_cum_prob - if (!method %in% valid_methods) - stop("Method shold be one of ", paste(valid_methods, collapse = ", ")) + return(weight) +} - if (method == "bb") { +#' Standard deviation of the spread kernel used by [calc_dist_weights()] +#' +#' Computes the standard deviation (in m) of the probability distribution +#' representing uncertainty in a bird's location at a given point along a +#' transition line, for one of the kernels supported by +#' [calc_dist_weights()]. Factored out of [calc_dist_weights()] so +#' [visualize_distance_weights()] can plot the same spread it uses without +#' duplicating the kernel-selection logic. +#' +#' @param dist_along_line How far along the line the point of interest is (m) +#' @param line_lengths How long the line is (m) +#' @param res_m The resolution of the associated bird flow model, used to +#' determine the nugget added to the variance to represent the uncertainty in +#' the starting and ending location of the transition. +#' @inheritParams calc_dist_weights +#' @return A vector of standard deviations (m), the same length as +#' `dist_along_line`. +#' @keywords internal +calc_dist_weights_sd <- function(dist_along_line, line_lengths, res_m, + kernel, gamma, kl, s1) { - # Hyperparameters - s1 <- 10 # Tune this visually looking at a few resulting distributions? + len <- line_lengths # length of great circle (m) + t <- dist_along_line # Distance from start of great circle to where the + # point projects onto the great circle (m) + # s2 is the standard deviation of the nugget + s2 <- res_m / 4 # This initial setting of 1/4 the cell width means the + # cell boundary (orthagonally) is at 2 standard deviations from the cell + # center, the corners will be at 2.8 SD. + if (kernel == "bb") { # Standard deviation - sd <- sqrt(s1^2 * t * (len - t) / len + s2^2) # SD in KM for this point. - - } - - if (method == "m3") { - sd <- sqrt(calc_martern_variance(t, len, k_m3, gamma = 40000, kl = 2000) + + sd <- sqrt(s1^2 * t * (len - t) / len + s2^2) # SD (m) for this point. + } else { + k_fun <- get(paste0("k_", kernel)) + sd <- sqrt(calc_martern_variance(t, len, k_fun, gamma = gamma, kl = kl) + s2^2) } - - # Calculate weight - weight <- rep(0, length(sd)) - in_range <- (d - r) < 1.96 * sd - low_cum_prob <- stats::pnorm(d[in_range] - r, sd = sd[in_range]) - high_cum_prob <- stats::pnorm(d[in_range] + r, sd = sd[in_range]) - weight[in_range] <- high_cum_prob - low_cum_prob - - return(weight) + sd } diff --git a/R/calc_euclidean_detection_rate.R b/R/calc_euclidean_detection_rate.R new file mode 100644 index 00000000..98df8d26 --- /dev/null +++ b/R/calc_euclidean_detection_rate.R @@ -0,0 +1,306 @@ + +#' Assign weights to points based on how close they are to the planar +#' (Euclidean) line associated with each transition +#' +#' This internal function is used to create a betweenness weights array +#' with dimensions [n_active(bf)](n_active()), +#' [n_active(bf)](n_active()), and `length(points)`. +#' The first two dimensions represent from and to cells of possible connections +#' between pairs of locations within the BirdFlow model and both have +#' `n_active()` elements. +#' The third dimension represents reference points that might be between each +#' connection. +#' Cell values are the weight to use when adding the transition probabilities +#' to the reference point BMTR +#' +#' If `points` are `NULL` they default to the center of all cells within +#' the BirdFlow extent that fall between any active cells. +#' This includes all cell centers within a convex hull +#' (in the model's native, projected CRS) around the active cells in `bf` +#' and thus is almost always more than just the active cells. +#' +#' Unlike [calc_spherical_detection_rate()], this projects points onto the +#' line segment connecting each pair of active cells directly in the model's +#' native (planar) CRS, using vectorized matrix algebra +#' ([project_points_simple()]) rather than great-circle math. This is +#' substantially faster and is the recommended continuous detection method. +#' +#' @param bf A BirdFlow model +#' @param points The points to evaluate betweenness on. If NULL the cell +#' centers of all the raster cells within the BirdFlow model that are between +#' active cells in the model will be used. This is calculated by comparing the +#' cell centers to a buffered convex hull around the active cell centers. +#' @param radius The probability density along a transect perpendicular to +#' the line and intersecting the point is summed over the interval within +#' `radius` meters of the point. `radius` defaults to half the cell +#' size (`mean(res(bf))/2`). +#' @param n_directions The number of (equally spaced) directional bins to +#' classify bearings into. Currently only `1` is supported. +#' @param skip_unconnected If `TRUE` then only connections that exist in `bf` +#' will be evaluated, and between matrix will erroneously indicate that +#' the weights associated with transitions that aren't used is always 0. +#' The resulting array can still be used with the model it was built for because +#' those missing connections would always have zero probability. +#' @param batch_size controls the number of movement lines that are processed +#' at a time. A smaller +#' `batch_size` will conserve memory at a slight performance cost. The number +#' of batches will be less than or equal to `n_active(bf)^2 / batch_size`. +#' @param check_radius If `TRUE` an error will be thrown if the radius +#' is not between the resolution and 1/4 the resolution of `bf`. Outside of +#' that range the algorithm is likely to yield distorted results. +#' `0.5 * mean(res(bf))` is the default, and recommended radius. +#' @param ... Additional arguments forwarded to [calc_dist_weights()], +#' allowing the spread kernel (`kernel`) and its hyperparameters +#' (`gamma`, `kl`, `s1`) to be tuned. +#' @return A list with: +#' \item{between}{An array with dimensions representing the +#' "from" location, the "to" location, and the `points`. Cells are weights and +#' will be non-zero if the radius around the point intersects 1.96 standard +#' deviations of the normal distribution of probabilities around the line.} +#' \item{points}{A data,frame of points that define the third dimension +#' in `between`. It is identical to the input `points` if they are not `NULL`. +#' Otherwise it will be a data frame with columns `x`, `y`, and `i` +#' corresponding to the third dimension in `between`. +#' `i` will be `NA` for points that are not within the mask but +#' fall between active cells.} +#' \item{radius}{The radius of the circle in meters.} +#' @seealso [is_between()], [calc_spherical_detection_rate()], and +#' [calc_bmtr()] +#' @keywords internal +calc_euclidean_detection_rate <- function(bf, points = NULL, radius = NULL, + n_directions = 1, + skip_unconnected = TRUE, + batch_size = 1e5, + check_radius = TRUE, ...) { + bf_msg("Generating between weights.\n") + + if (!requireNamespace("SparseArray", quietly = TRUE)) { + stop("The SparseArray package is required to use is_between(). ", + "Please install it prior to calling this function.") + } + + if (!n_directions == 1) { + stop("Currently only one direction is supported.") + } + # Dimensions + # 1 from location n_active() + # 2 to location n_active() + # 3 points (if NULL use cell centers) + # 4 (pending) n_direction, directional bins. + + if (is.null(radius)) { + radius <- mean(res(bf)) / 2 + } else if (check_radius) { + # Analysis of the effect of changing radius is in test-calc_bmtr.R + radius_cells <- radius / mean(res(bf)) # radius converted to cells + if (radius_cells <= 0.25 || radius_cells >= 1) { + stop("radius should be less than the resolution and more than 1/4 the ", + "resolution or BMTR is likely to be biased. ", + "Set check_radius to FALSE to ignore this advice.") + } + } + + # Hull buffer is only used for selecting default points, not for the + # detection radius used in the weighting below. + hull_buffer <- 2 * mean(res(bf)) + + if (is.null(points)) { + bf_msg(" Creating points\n") + # Create points at all cell centers in the rectangular raster + points <- rasterize_distr(get_distr(bf, 1), bf, format = "dataframe") + points <- points[, c("x", "y", "i")] + active <- points[!is.na(points$i), , drop = FALSE] + + bf_msg(" Preparing sf objects\n") + # Convert active and points to sf objects + active_sf <- active |> + sf::st_as_sf(coords = c("x", "y"), crs = crs(bf)) + points_sf <- points |> + sf::st_as_sf(coords = c("x", "y"), crs = crs(bf)) + + bf_msg(" Creating buffered convex hull\n") + # Make a buffered convex hull around the active cells + hull <- sf::st_union(active_sf) |> + sf::st_convex_hull() |> + sf::st_buffer(dist = units::set_units(hull_buffer, "m")) + + bf_msg(" Selecting points inside hull\n") + # Keep the points inside the hull + sv <- points_sf |> + sf::st_intersects(y = hull, sparse = FALSE) |> + as.vector() + + if (FALSE) { + # visualize + plot(hull) + plot(points, col = "grey", add = TRUE) + plot(points[sv, ], col = "black", add = TRUE) + plot(active, col = "red", add = TRUE) + } + + # Subset to the cells that are active or between other active cells + # Note this is in back in original projection + points <- points[sv, , drop = FALSE] + } + + bf_msg(" Initializing arrays. \n") + + # Initialize sparse array (All FALSE) to hold betweenness + # This is a hack because there's no creation method that allows setting a + # dimension, but I can make an empty "random" array. + between <- SparseArray::randomSparseArray( + dim = c(n_active(bf), n_active(bf), nrow(points)), + density = 0) != 0 + dimnames(between) <- list(from = paste0("F_", seq_len(n_active(bf))), + to = paste0("T_", seq_len(n_active(bf))), + loc = paste0("L_", seq_len(nrow(points)))) + + + # Generate table of active cell x, y, and i (in birdflow CRS) + active <- i_to_xy(seq_len(n_active(bf)), bf) + active$i <- seq_len(n_active(bf)) + + # All possible pairs of active points + all_pairs <- expand.grid(from = seq_len(n_active(bf)), + to = seq_len(n_active(bf))) + all_pairs$id <- NA_character_ + sv <- all_pairs$from < all_pairs$to + all_pairs$id[sv] <- paste(all_pairs$from[sv], "-", all_pairs$to[sv]) + all_pairs$id[!sv] <- paste(all_pairs$to[!sv], "-", all_pairs$from[!sv]) + + # Unique pairs of points disregarding order + pairs <- all_pairs[!duplicated(all_pairs$id), , drop = FALSE] + + # Drop self, self pairs. + # They represent stop overs or seasonal residence, not migratory movement + pairs <- pairs[pairs$from != pairs$to, , drop = FALSE] + + # Drop pairs that aren't ever connected + if (skip_unconnected) { + # Create matrix indicating which active cells are connected to + # each other via a non-zero marginal at any timestep + # With sparse models this will eliminate a lot of connections + # With non-sparse models it will still eliminate some + # due to dynamic masking. + ever_connected <- matrix(FALSE, n_active(bf), n_active(bf)) + dm <- get_dynamic_mask(bf) + mi <- bf$marginals$index + mi <- mi[mi$direction == "forward", ] + + for (i in seq_len(nrow(mi))) { + from_dm <- dm[, mi$from[i]] + to_dm <- dm[, mi$to[i]] + marg <- get_marginal(bf, mi$marginal[i]) + ever_connected[from_dm, to_dm] <- as.matrix(marg != 0) + } + ever_connected <- ever_connected | t(ever_connected) # backwards counts + + connected_pairs <- data.frame(from = row(ever_connected)[ever_connected], + to = col(ever_connected)[ever_connected]) + connected_pairs$id <- with(connected_pairs, paste(from, "-", to)) + + pairs <- pairs[pairs$id %in% connected_pairs$id, , drop = FALSE] + + } + + # Project points onto line segments in batches to bound memory use of the + # (n points) x (m lines) projection matrices. + p_mat <- as.matrix(points[, c("x", "y")]) + m <- nrow(p_mat) # number of points + n_total <- nrow(pairs) + batch_ids <- split(seq_len(n_total), + ceiling(seq_len(n_total) / batch_size)) + + for (b in seq_along(batch_ids)) { + idx <- batch_ids[[b]] + pairs_batch <- pairs[idx, , drop = FALSE] + + s_mat <- as.matrix(active[pairs_batch$from, c("x", "y")]) + e_mat <- as.matrix(active[pairs_batch$to, c("x", "y")]) + + n <- nrow(s_mat) + + v_mat <- e_mat - s_mat + vv <- rowSums(v_mat^2) + + projection <- project_points_simple(s = s_mat, e = e_mat, p = p_mat, + m = m, n = n, v = v_mat, vv = vv) + dist_to_line <- as.vector(projection$dist_to_line) + dist_along_line <- as.vector(projection$dist_along) + line_lengths <- rep(sqrt(vv), each = m) + + line_index <- rep(seq_len(n), each = m) + point_index <- rep(seq_len(m), times = n) + + valid <- !is.na(dist_along_line) & dist_along_line > 0 & + dist_along_line < line_lengths + + weights <- calc_dist_weights(dist_to_line[valid], + dist_along_line[valid], + line_lengths[valid], + res_m = mean(res(bf)), + radius_m = radius, + ...) + + add <- cbind(pairs_batch$from[line_index[valid]], + pairs_batch$to[line_index[valid]], + point_index[valid], + weights) + add <- add[add[, 4] != 0, , drop = FALSE] + + between[add[, c(1, 2, 3)]] <- add[, 4] + between[add[, c(2, 1, 3)]] <- add[, 4] + + bf_msg(" ", round(b / length(batch_ids) * 100, 2), "%\n") + } + + return(list(between = between, points = points, radius = radius)) +} + +#' Project points onto line segments (Euclidean/planar geometry) +#' +#' Vectorized computation of the distance from each point in `P` to each +#' line segment defined by `S` (starts) and `E` (ends), and the distance +#' along the segment to the projection of the point. +#' +#' @param s An n x 2 matrix of segment start coordinates. +#' @param e An n x 2 matrix of segment end coordinates. +#' @param p An m x 2 matrix of point coordinates. +#' @param m Number of points (`nrow(p)`). +#' @param n Number of line segments (`nrow(s)`). +#' @param v `e - s`, precomputed segment vectors. +#' @param vv `rowSums(v^2)`, precomputed squared segment lengths. +#' @param clamp If `TRUE` constrain the projected fraction to `[0, 1]`. If +#' `FALSE` (default) projections falling outside the segment are set to `NA`. +#' @param tol Segments with squared length at or below this threshold are +#' treated as having zero length (projected fraction set to 0). +#' @return A list with `dist_to_line` and `dist_along`, each a vector of +#' length `m * n` (point index varying fastest). +#' @keywords internal +project_points_simple <- function(s, e, p, m, n, v, vv, + clamp = FALSE, tol = 0) { + + # 1. Vectorized dot products (m x n matrix) giving the fraction of each + # line's length at which each point projects onto that line. + t_frac <- (p %*% t(v) - rep(rowSums(s * v), each = m)) / rep(vv, each = m) + + # Zero-length lines project everything to their start. + t_frac[, vv <= tol] <- 0 + + if (clamp) { + t_frac <- pmax(0, pmin(1, t_frac)) + } else { + t_frac[t_frac < 0 | t_frac > 1] <- NA + } + + # Coordinates of the projected point for every (point, line) combination. + projx <- rep(s[, 1], each = m) + t_frac * rep(v[, 1], each = m) + projy <- rep(s[, 2], each = m) + t_frac * rep(v[, 2], each = m) + + dist_to_line <- sqrt((rep(p[, 1], times = n) - projx)^2 + + (rep(p[, 2], times = n) - projy)^2) + + dist_along <- t_frac * rep(sqrt(vv), each = m) + + list(dist_to_line = dist_to_line, dist_along = dist_along) +} diff --git a/R/weight_between.R b/R/calc_spherical_detection_rate.R similarity index 85% rename from R/weight_between.R rename to R/calc_spherical_detection_rate.R index ab314b97..de0cb943 100644 --- a/R/weight_between.R +++ b/R/calc_spherical_detection_rate.R @@ -26,13 +26,17 @@ if (FALSE) { #' (in spherical coordinates) around the active cells in `bf` and thus is #' almost always more than just the active cells. #' -#' `weight_between()` and `is_between()` should only differ slightly in -#' their results when calculated on the same model. +#' `calc_spherical_detection_rate()` and `is_between()` should only differ +#' slightly in their results when calculated on the same model. #' #' 1, The number of included reference points may differ. #' 2. The betweenness array will be real weights from 0 to 1 with -#' `weight_between()` and logical with `is_between()` +#' `calc_spherical_detection_rate()` and logical with `is_between()` #' +#' This is the great-circle (spherical) counterpart to +#' [calc_euclidean_detection_rate()]. It's much slower and not recommended +#' for production use, but is kept to allow assessing the impact of switching +#' from spherical to Euclidean geometry. #' #' @param bf A BirdFlow model #' @param points The points to evaluate betweenness on. If NULL the cell @@ -54,6 +58,13 @@ if (FALSE) { #' at a time. A smaller #' `batch_size` will conserve memory at a slight performance cost. The number #' of batches will be less than or equal to `n_active(bf)^2 / batch_size`. +#' @param check_radius If `TRUE` an error will be thrown if the radius +#' is not between the resolution and 1/4 the resolution of `bf`. Outside of +#' that range the algorithm is likely to yield distorted results. +#' `0.5 * mean(res(bf))` is the default, and recommended radius. +#' @param ... Additional arguments forwarded to [calc_dist_weights()], +#' allowing the spread kernel (`kernel`) and its hyperparameters +#' (`gamma`, `kl`, `s1`) to be tuned. #' @return A list with: #' \item{between}{An array with dimensions representing the #' "from" location, the "to" location, and the `points`. Cells are weights and @@ -67,11 +78,14 @@ if (FALSE) { #' `i` will be `NA` for points that are not within the mask but #' fall between active cells.} #' \item{radius}{The radius of the circle in meters.} -#' @seealso [is_between()] and [calc_bmtr()] +#' @seealso [is_between()], [calc_euclidean_detection_rate()], and +#' [calc_bmtr()] #' @keywords internal -weight_between <- function(bf, weight_fun = NULL, points = NULL, radius = NULL, - n_directions = 1, skip_unconnected = TRUE, - batch_size = 1e5, ...) { +calc_spherical_detection_rate <- function(bf, points = NULL, radius = NULL, + n_directions = 1, + skip_unconnected = TRUE, + batch_size = 1e5, + check_radius = TRUE, ...) { bf_msg("Generating between weights.\n") @@ -94,11 +108,21 @@ weight_between <- function(bf, weight_fun = NULL, points = NULL, radius = NULL, # 3 points (if NULL use cell centers) # 4 (pending) n_direction, directional bins. - # Radius is just used for buffering the convex hull - # I'm not sure what the buffer should be so am just setting - # it to the cell dimension. I may have to loop back later and figure - # something else out. - radius <- mean(res(bf)) + if (is.null(radius)) { + radius <- mean(res(bf)) / 2 + } else if (check_radius) { + # Analysis of the effect of changing radius is in test-calc_bmtr.R + radius_cells <- radius / mean(res(bf)) # radius converted to cells + if (radius_cells <= 0.25 || radius_cells >= 1) { + stop("radius should be less than the resolution and more than 1/4 the ", + "resolution or BMTR is likely to be biased. ", + "Set check_radius to FALSE to ignore this advice.") + } + } + + # Hull buffer is only used for selecting default points, not for the + # detection radius used in the weighting below. + hull_buffer <- mean(res(bf)) if (is.null(points)) { bf_msg(" Creating points\n") @@ -121,7 +145,7 @@ weight_between <- function(bf, weight_fun = NULL, points = NULL, radius = NULL, # Make a buffered convex hull around the active cells in lat lon hull <- sf::st_union(active_sph) |> sf::st_convex_hull() |> - sf::st_buffer(hull, dist = units::set_units(radius, "m")) + sf::st_buffer(dist = units::set_units(hull_buffer, "m")) # Selection vector for points that are active or between active cells sv <- points_sph |> @@ -271,14 +295,13 @@ weight_between <- function(bf, weight_fun = NULL, points = NULL, radius = NULL, valid_line_lengths <- rep(line_lengths, each = n_points)[not_end] valid_line_index <- rep(seq_len(n_lines), each = n_points)[not_end] valid_point_index <- rep(seq_len(n_points), times = n_lines)[not_end] - # Placeholder weight fun weights <- calc_dist_weights(valid_dist_to_line, valid_dist_along_line, valid_line_lengths, res_m = mean(res(bf)), radius_m = radius, - method = "m3") + ...) # Add is a matrix defining the new cells with weight to @@ -322,11 +345,6 @@ weight_between <- function(bf, weight_fun = NULL, points = NULL, radius = NULL, # extract weights associated with all points relative to this line points_sph$weight <- between[pairs$from[i], pairs$to[i], ] - # PLot touched points in green - #plot(points_sph[ points_sph$weight > 0, ], cex = 2, - # col = "green", add = TRUE) - - # Plot touched points colored by weight pal <- ebirdst::ebirdst_palettes(n = 256, "weekly") col_index <- points_sph$weight / max(points_sph$weight) * 255 + 1 @@ -337,7 +355,7 @@ weight_between <- function(bf, weight_fun = NULL, points = NULL, radius = NULL, plot(points_sph, col = cols, add = TRUE, pch = 19) n_legend <- 10 - legend_index <- seq(from = 1, to = 256, length.out = nlegend) |> round() + legend_index <- seq(from = 1, to = 256, length.out = n_legend) |> round() legend_colors <- pal[legend_index] legend_values <- (legend_index - 1) / 255 * max(points_sph$weight) legend_values <- signif(legend_values, 2) diff --git a/R/visualize_distance_weights.R b/R/visualize_distance_weights.R new file mode 100644 index 00000000..2d5ec971 --- /dev/null +++ b/R/visualize_distance_weights.R @@ -0,0 +1,147 @@ +#' Visualize the spread kernel used by [calc_dist_weights()] +#' +#' `visualize_distance_weights()` plots the shape of the spread kernel used +#' by [calc_dist_weights()] (and, through it, the `"continuous"` and +#' `"continuous-spherical"` methods of [calc_bmtr()]) for one or more +#' transition-line lengths. It's meant as a tool for choosing `kernel` and +#' its hyperparameters (`gamma`, `kl`, `s1`) before running an expensive +#' `calc_bmtr()` call. +#' +#' @param line_lengths A vector of transition-line lengths (m) to visualize +#' the spread kernel for. +#' @param radius_m The detection radius (m). Defaults to `res_m / 2` (the +#' same default used by [calc_bmtr()]'s continuous methods). For +#' `type = "envelope"` this is only shown as a reference line, and only if +#' explicitly supplied. For `type = "raster"` it's used, as in +#' [calc_dist_weights()], to determine the band of probability density +#' that forms the plotted weight. +#' @param res_m The resolution of the associated bird flow model (m), used +#' to determine the nugget added to the variance. See [calc_dist_weights()]. +#' @inheritParams calc_dist_weights +#' @param type The style of plot: +#' \describe{ +#' \item{`"envelope"`}{(default) For each line length, plots the region +#' within 1.96 standard deviations of the line (in the kernel's probability +#' distribution) as a band along the line. This represents 95 percent of +#' the density.} +#' \item{`"raster"`}{For each line length, plots the actual weight returned +#' by [calc_dist_weights()] over a grid of points around the line, as a +#' heatmap. Slower, but reflects `radius_m`'s effect on the weight as well +#' as the kernel's spread.} +#' } +#' @param n The number of grid points used along (and, for `type = "raster"`, +#' across) each line. +#' +#' @return A `ggplot2` object. It can be displayed with `print()`. +#' @seealso [calc_dist_weights()], [calc_bmtr()] +#' @export +#' +#' @examples +#' visualize_distance_weights() +#' visualize_distance_weights(type = "raster", kernel = "bb", s1 = 20) +visualize_distance_weights <- function(line_lengths = c(2000, 1000, 500, 200), + radius_m = NULL, + res_m = 1000, + kernel = c("m3", "bb", "m1", "m5", "sq"), + gamma = 3e10, kl = 9e5, s1 = 200, + type = c("envelope", "raster"), + n = 100) { + kernel <- match.arg(kernel) + type <- match.arg(type) + + radius_supplied <- !is.null(radius_m) + if (!radius_supplied) { + radius_m <- res_m / 2 + } + + kernel_names <- c(m1 = "Matern 1/2", m3 = "Matern 3/2", m5 = "Matern 5/2", + sq = "Squared exponential", bb = "Brownian bridge") + hp_text <- if (kernel == "bb") { + paste0("s1 = ", s1) + } else { + paste0("gamma = ", gamma, ", kl = ", kl) + } + plot_title <- paste0("Spread kernel: ", kernel_names[[kernel]]) + plot_subtitle <- paste0("kernel = \"", kernel, "\" (", hp_text, ")") + + line_lengths <- sort(line_lengths, decreasing = TRUE) + length_levels <- paste(line_lengths, "m") + length_label <- factor(length_levels, levels = length_levels) + + if (type == "envelope") { + bands <- lapply(seq_along(line_lengths), function(i) { + len <- line_lengths[i] + dist_along <- seq(0, len, length.out = n) + sd_km <- calc_dist_weights_sd(dist_along, rep(len, n), res_m, + kernel = kernel, gamma = gamma, kl = kl, + s1 = s1) / 1000 + data.frame(dist_along_km = dist_along / 1000, + ymin = -1.96 * sd_km, ymax = 1.96 * sd_km, + length_label = length_label[i]) + }) + df <- do.call(rbind, bands) + + p <- ggplot2::ggplot(df, ggplot2::aes(x = .data$dist_along_km, + ymin = .data$ymin, + ymax = .data$ymax, + fill = .data$length_label)) + + ggplot2::geom_ribbon(alpha = 0.5) + + ggplot2::scale_fill_viridis_d(name = "Line length") + + if (radius_supplied) { + p <- p + ggplot2::geom_hline(yintercept = c(-radius_m, radius_m) / 1000, + linetype = "dashed") + } + } else { + # Raster + max_sd_m <- vapply(line_lengths, function(len) { + dist_along <- seq(0, len, length.out = n) + max(calc_dist_weights_sd(dist_along, rep(len, n), res_m, + kernel = kernel, gamma = gamma, kl = kl, + s1 = s1)) + }, numeric(1)) + y_max_m <- 3 * max(max_sd_m) + dist_to_seq <- seq(-y_max_m, y_max_m, length.out = n) + + dist_to_width_km <- diff(dist_to_seq)[1] / 1000 + + grids <- lapply(seq_along(line_lengths), function(i) { + len <- line_lengths[i] + dist_along <- seq(0, len, length.out = n) + grid <- expand.grid(dist_along_line = dist_along, + dist_to_line = dist_to_seq) + grid$weight <- calc_dist_weights( + dist_to_line = abs(grid$dist_to_line), + dist_along_line = grid$dist_along_line, + line_lengths = len, radius_m = radius_m, res_m = res_m, + kernel = kernel, gamma = gamma, kl = kl, s1 = s1) + grid$length_label <- length_label[i] + grid$dist_along_width_km <- (len / (n - 1)) / 1000 + grid + }) + df <- do.call(rbind, grids) + df$dist_along_km <- df$dist_along_line / 1000 + df$dist_to_km <- df$dist_to_line / 1000 + + # geom_tile() (rather than geom_raster()) since the along-line grid + # spacing differs between facets (each spans a different line length). + p <- ggplot2::ggplot(df, ggplot2::aes(x = .data$dist_along_km, + y = .data$dist_to_km, + width = .data$dist_along_width_km, + height = dist_to_width_km, + fill = .data$weight)) + + ggplot2::geom_tile() + + ggplot2::facet_wrap(~length_label) + + ggplot2::scale_fill_viridis_c(name = "Weight") + } + + p + + ggplot2::labs(x = "Distance along line (km)", + y = "Perpendicular distance from line (km)", + title = plot_title, subtitle = plot_subtitle) + + ggplot2::coord_fixed() + + ggplot2::scale_y_continuous( + breaks = scales::breaks_width(res_m / 1000 * 2)) + + ggplot2::scale_x_continuous( + breaks = scales::breaks_width(res_m / 1000 * 2)) +} diff --git a/data-raw/kernel_parameter_tuning.Rmd b/data-raw/kernel_parameter_tuning.Rmd new file mode 100644 index 00000000..7f5012c5 --- /dev/null +++ b/data-raw/kernel_parameter_tuning.Rmd @@ -0,0 +1,215 @@ +--- +title: "BMTR spread kernel: choosing gamma/kl/s1" +date: "2026-07-28" +output: html_document +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE, + fig.width = 7, fig.height = 5) +devtools::load_all("..", quiet = TRUE) +library(ggplot2) +``` + +## What this is + +A record of how the spread kernel hyperparameter defaults (`gamma`, `kl`, +`s1`) used by `calc_dist_weights()` were chosen, and the visualizations they +were chosen against. It is the counterpart to +`data-raw/martern_parameters.R`, which served the same purpose for the +earlier defaults. + +This is a development record, not package documentation — `data-raw/` is +`.Rbuildignore`d, so nothing here ships with the package or is built by +`R CMD check`. It calls package internals directly via `devtools::load_all()`. + +## The defaults + +```r +gamma = 3e10 # "m3", "m1", "m5", "sq" kernels +kl = 9e5 # "m3", "m1", "m5", "sq" kernels +s1 = 200 # "bb" kernel +``` + +Note `gamma` and `kl` where only tuned for the default `m3` kernel. +Likely anytime another of the kernels are used they will need to +be specified as well. + +These are the values set in `R/calc_distance_weights.R` and +`R/visualize_distance_weights.R`. They were chosen by eye against the +envelope plots below, nudging each until the tapering shape looked right. + + +**These hyperparameters are not resolution dependent.** They describe our +uncertainty about a biological process — how far a bird's actual path is +likely to stray from the straight line between two cells — and that +uncertainty is a property of the bird, not of the grid we happen to model it +on. So they are fixed and don't depend on `res_m`. + +Model resolution enters the kernel in exactly one place, separately from +these three: the nugget (`res_m / 4`), which represents uncertainty about +where within a cell the transition actually started and ended. That part +*is* resolution dependent, and it scales itself. + +`res_m = 100000` is used throughout this document only as a reference scale +for reading the plots — a familiar size to think about distances against, +since BirdFlow models typically use 100 or 150 km cells. + + +## Checking the shape at a reference scale + +```{r} +res_m <- 100000 # reference scale only: typical BirdFlow cell size +cells <- c(1, 2, 5, 10, 15, 30) +lens <- cells * res_m + +gamma <- 3e10 +kl <- 9e5 +s1 <- 200 +nugget <- res_m / 4 + +sd_table <- do.call(rbind, lapply(lens, function(len) { + data.frame( + cells = len / res_m, + len_km = len / 1000, + m3_sd_end = calc_dist_weights_sd(0, len, res_m, kernel = "m3", + gamma = gamma, kl = kl, s1 = NA), + m3_sd_mid = calc_dist_weights_sd(len / 2, len, res_m, kernel = "m3", + gamma = gamma, kl = kl, s1 = NA), + bb_sd_mid = calc_dist_weights_sd(len / 2, len, res_m, kernel = "bb", + gamma = NA, kl = NA, s1 = s1) + ) +})) +sd_table +``` + +`m3_sd_end` sits exactly at the nugget (`res_m / 4`) for every length — the +bridge conditioning pins the variance to the nugget at the endpoints, so +this is the one column that would move if the model resolution changed. The +`m3_sd_mid` column grows with line length and asymptotes once the line is +long relative to `kl` — toward `sqrt(gamma + nugget^2)` = +`r round(sqrt(gamma + nugget^2) / 1000)` km — but at these `gamma`/`kl` +values that approach is gradual: it's still climbing at a 1500 km line +(~120 km) and doesn't flatten out until well beyond the range shown here. +`bb_sd_mid` keeps growing with `sqrt(length)` (no plateau; that's inherent +to a Brownian bridge with no lengthscale), reaching ~125 km by 1500 km. +Both are biologically plausible: a bird's path straying on the order of +100 km from the straight line over a multi-week migration leg, tapering to +the endpoints and widest in between. + +## Envelope plots + +The envelope plots show the region within 1.96 standard deviations of the +line — 95% of the kernel's density. These are the plots the defaults +were actually tuned against. + +```{r, fig.height=6, fig.width=10} +n_cells <- c(1, 2, 5, 10, 15, 30) +line_lengths <- n_cells * res_m + +visualize_distance_weights(res_m = res_m, line_lengths = line_lengths, + type = "envelope", kernel = "m3", + gamma = gamma, kl = kl) +``` + +```{r, fig.height=4, fig.width=10} +visualize_distance_weights(res_m = res_m, line_lengths = line_lengths, + type = "envelope", kernel = "bb", s1 = s1) +``` + +Both taper from zero at the endpoints and widen toward the middle. + +## Weight rasters + +The envelope traces the kernel's spread but says nothing about the weight +`calc_dist_weights()` actually returns, which also depends on the detection +radius: the weight is the probability mass falling within `radius_m` of the +point, so it peaks on the line and falls off with the density. The raster +plots below show that weight directly, with `radius_m` at its default of +`res_m / 2`. Fewer and shorter line lengths than above, since each facet is +a 100x100 grid of `calc_dist_weights()` calls. + +```{r, fig.height=4, fig.width=10} +raster_lengths <- c(2, 5, 10) * res_m + +visualize_distance_weights(res_m = res_m, line_lengths = raster_lengths, + radius_m = res_m / 2, type = "raster", + kernel = "m3", gamma = gamma, kl = kl) +``` + +```{r, fig.height=4, fig.width=10} +visualize_distance_weights(res_m = res_m, line_lengths = raster_lengths, + radius_m = res_m / 2, type = "raster", + kernel = "bb", s1 = s1) +``` + +Weight is highest near the start and end of each line, where the spread is +pinched to the nugget and so nearly all the density falls inside `radius_m`, +and lowest mid-line where the same density is spread over a much wider band. +That is the intended behavior: a transition is most informative about +locations near its endpoints. + +## Sweeping `gamma` at a fixed `kl` + +A closer look at how sensitive the `"m3"` envelope is to `gamma`, holding +`kl` at the chosen default (`9e5`): three line lengths as column facets, +`gamma` from `1e10` to `5e10` in `1e10` steps as row facets. + +```{r, fig.height=9, fig.width=8} +n_cells_sweep <- c(2, 5, 15) +line_lengths_sweep <- n_cells_sweep * res_m +gamma_values <- seq(1e10, 5e10, by = 1e10) +kl_fixed <- 9e5 + +length_levels <- paste(line_lengths_sweep / 1000, "km") +gamma_levels <- paste0("gamma = ", format(gamma_values, scientific = TRUE)) + +sweep <- do.call(rbind, lapply(gamma_values, function(g) { + do.call(rbind, lapply(line_lengths_sweep, function(len) { + dist_along <- seq(0, len, length.out = 100) + sd_m <- calc_dist_weights_sd(dist_along, rep(len, 100), res_m, + kernel = "m3", gamma = g, kl = kl_fixed, + s1 = NA) + data.frame( + dist_along_km = dist_along / 1000, + ymin = -1.96 * sd_m / 1000, + ymax = 1.96 * sd_m / 1000, + length_label = factor(paste(len / 1000, "km"), levels = length_levels), + gamma_label = factor(paste0("gamma = ", format(g, scientific = TRUE)), + levels = gamma_levels) + ) + })) +})) + +ggplot(sweep, aes(x = dist_along_km, ymin = ymin, ymax = ymax)) + + geom_ribbon(fill = "steelblue", alpha = 0.6) + + facet_grid(gamma_label ~ length_label) + + coord_fixed() + + labs(x = "Distance along line (km)", + y = "Perpendicular distance from line (km)", + title = "\"m3\" kernel envelope: gamma sweep at kl = 9e5", + subtitle = paste0("res_m = ", res_m)) + + scale_y_continuous(breaks = scales::breaks_width(res_m / 1000 * 2)) + + scale_x_continuous(breaks = scales::breaks_width(res_m / 1000 * 2)) +``` + +Across this range, `gamma` mainly sets how wide the envelope gets once the +line is long enough for the kernel to fully decorrelate (the 1500 km +column), while the 200 and 500 km columns are still partway through the +taper and change more subtly. `kl = 9e5` is what produces that gradual +widen-then-plateau shape across all three lengths, rather than saturating +immediately (small `kl`) or barely moving at all (`kl` much larger than the +longest line). + +## Takeaway + +- Final defaults, chosen by eye against the envelope plots above: + `gamma = 3e10`, `kl = 9e5` for the `"m3"`/`"m1"`/`"m5"`/`"sq"` kernels, + and `s1 = 200` for `"bb"`. +- Applied as fixed literal defaults in `R/calc_distance_weights.R` and + `R/visualize_distance_weights.R`. Fixed rather than `res_m`-derived + because they describe biological movement uncertainty, which is + independent of the model grid; only the nugget (`res_m / 4`) is + resolution dependent. +- Both kernels give spread on the order of a grid cell, tapering to the + nugget at the endpoints, with weight concentrated near those endpoints. diff --git a/data-raw/martern_parameters.R b/data-raw/martern_parameters.R index 102a46f4..77517a3f 100644 --- a/data-raw/martern_parameters.R +++ b/data-raw/martern_parameters.R @@ -1,5 +1,9 @@ # nolint start +# This document has been superceded by data-raw/kernel_parameter_tuning.Rmd please use that one for understanding +# current parameter values. This one may not run properly or the interpretation of units (km vs m) may +# have changed for some of the parameters. + # Note kl was originally l but that lead to warnings about partial argument # matching ambiguity so I changed it to kl for kernel length hyper parameter @@ -10,7 +14,7 @@ plot_95 <- function(T, k, gamma, kl, add = TRUE) { # kl <- 1/2 * T # gamma <- 4 * T xs <- seq(0, T, length.out = 50) - ys <- sqrt(calc_variance(xs, T, k = k, gamma = gamma, kl = kl)) * 1.96 + ys <- sqrt(calc_martern_variance(xs, T, k = k, gamma = gamma, kl = kl)) * 1.96 # xs <- xs + 0.5 * T if (!add) { diff --git a/inst/WORDLIST b/inst/WORDLIST index bc51d42e..ec78882b 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -171,3 +171,8 @@ vectorized lossy reconstructible unclipped +Matern +lengthscale +Lengthscale +decorrelated +decorrelates diff --git a/man/animate_bmtr.Rd b/man/animate_bmtr.Rd index 0f5450c9..72d93baa 100644 --- a/man/animate_bmtr.Rd +++ b/man/animate_bmtr.Rd @@ -57,13 +57,31 @@ Animate migration traffic rates produced by \code{\link[=calc_bmtr]{calc_bmtr()} \dontrun{ bf <- BirdFlowModels::amewoo -bmtr <- calc_bmtr(bf) +# Binary detection along shortest path +bmtr <- calc_bmtr(bf) plot_bmtr(bmtr, bf) animate_bmtr(bmtr, bf) + +# Continuous detection with a wider spread kernel +# (gamma defaults to 3e10) +bmtr2 <- calc_bmtr(bf, method = "continuous", gamma = 6e10) +animate_bmtr(bmtr2, bf) + +# Visualize spread for the continuous detection with wider spread +visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf), + gamma = 6e10, res_m = xres(bf)) + +# Visualize spread for the continuous detection with default values +visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf), + res_m = xres(bf)) + } + + + } \seealso{ \code{\link[=calc_bmtr]{calc_bmtr()}},\code{\link[=plot_bmtr]{plot_bmtr()}} diff --git a/man/calc_bmtr.Rd b/man/calc_bmtr.Rd index a54968c1..30806f74 100644 --- a/man/calc_bmtr.Rd +++ b/man/calc_bmtr.Rd @@ -12,7 +12,8 @@ calc_bmtr( format = NULL, batch_size = 5e+05, check_radius = TRUE, - weighted = FALSE + method = c("binary", "continuous", "continuous-spherical"), + ... ) } \arguments{ @@ -23,10 +24,14 @@ calc_bmtr( or fall between two active cells. Otherwise a data frame with \code{x} and \code{y} columns containing point coordinates in \link[terra:crs]{crs(bf)}.} -\item{radius}{The radius in meters around the points used to assess whether -a movement line passes by (or through) the point. If a point is farther than -\code{radius} from a great circle line between two cells centers then it is not -between them.} +\item{radius}{The radius in meters around the points used to assess the +detection rate for a movement at the point. +With \code{method = "binary"}, if a point is within \code{radius} of the great +circle line between two cell centers then the movement is detected at +that point. For the two continuous detection methods there is a +probability distribution for the location of the bird as it passes by the +point, and the radius defines the band over which that probability is +integrated, giving the detection rate.} \item{n_directions}{The number of directional bins to use for recording movement direction. Must be either \code{1} indicating no direction information @@ -61,13 +66,30 @@ is not between the resolution and 1/4 the resolution of \code{bf}. Outside of that range the algorithm is likely to yield distorted results. \code{0.5 * mean(res(bf))} is the default, and recommended radius.} -\item{weighted}{If \code{FALSE} use the original and quicker version of bmtr -that sums all the marginal probability for transitions that pass within a -fixed distance of the point. If \code{TRUE} assign a weight to the point and -transition combo that then is multiplied by the marginal probability before -summing. This argument is experimental but the default value is identical -to the old version. The argument name and behavior when set to \code{TRUE} may -change.} +\item{method}{The detection model used to determine how much of a +transition's movement counts towards a point's BMTR: +\describe{ +\item{\code{"binary"}}{(default) Fast and deterministic. A movement line either +does or does not pass within \code{radius} of the point, per \code{\link[=is_between]{is_between()}}.} +\item{\code{"continuous"}}{Assigns a continuous weight (0 to 1) based on the +probability that a bird's actual path, modeled as spreading away from the +straight line between two cells, passes within \code{radius} of the point. Uses +planar (Euclidean) geometry in the model's native CRS, and is the +recommended detection model when continuous weighting is desired.} +\item{\code{"continuous-spherical"}}{The same continuous weighting as +\code{"continuous"}, but computed with great-circle (spherical) geometry +instead. Much slower, and not recommended for routine use; kept to allow +assessing the impact of switching from spherical to Euclidean geometry.} +}} + +\item{...}{For \code{method = "continuous"} or \code{"continuous-spherical"}, +additional arguments forwarded to \code{\link[=calc_dist_weights]{calc_dist_weights()}} to control the +spread kernel used to model uncertainty in a bird's path: \code{kernel}, +\code{gamma}, \code{kl}, and \code{s1}. See \code{\link[=calc_dist_weights]{calc_dist_weights()}} for the full list of +supported kernels and their hyperparameters, and +\code{\link[=visualize_distance_weights]{visualize_distance_weights()}} to explore how they shape the spread +before running this (potentially expensive) function. Not applicable, +and an error, for \code{method = "binary"}.} } \value{ See \code{format} argument. @@ -94,10 +116,14 @@ Multiplying the result by the total population would yield: \section{Limitations}{ -\code{calc_bmtr()} makes the incorrect simplifying assumption -that birds follow the shortest (great circle) path -between the center of the the source and destination raster cells. Caution -should be used when interpreting the results especially around +\code{calc_bmtr()} makes an incorrect simplifying assumption about the path +birds take: either that they follow the shortest path between the centers +of the source and destination raster cells (\code{method = "binary"}), or that +they are distributed symmetrically in a Gaussian distribution around that +path (the two continuous methods). That path is a great circle in all +cases except \code{method = "continuous"}, which uses a straight line in the +model's projected CRS. +Caution should be used when interpreting the results especially around major geographic features such as coasts, large lakes, mountain ranges, and ecological system boundaries that might result in non-linear migration paths. @@ -122,11 +148,33 @@ represents. \dontrun{ bf <- BirdFlowModels::amewoo -bmtr <- calc_bmtr(bf) +# Binary detection along shortest path +bmtr <- calc_bmtr(bf) plot_bmtr(bmtr, bf) animate_bmtr(bmtr, bf) + +# Continuous detection with a wider spread kernel +# (gamma defaults to 3e10) +bmtr2 <- calc_bmtr(bf, method = "continuous", gamma = 6e10) +animate_bmtr(bmtr2, bf) + +# Visualize spread for the continuous detection with wider spread +visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf), + gamma = 6e10, res_m = xres(bf)) + +# Visualize spread for the continuous detection with default values +visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf), + res_m = xres(bf)) + } + + + +} +\seealso{ +\code{\link[=visualize_distance_weights]{visualize_distance_weights()}} to explore the spread kernel +hyperparameters accepted via \code{...}. } diff --git a/man/calc_dist_weights.Rd b/man/calc_dist_weights.Rd index 19b72fd0..f1b74563 100644 --- a/man/calc_dist_weights.Rd +++ b/man/calc_dist_weights.Rd @@ -10,7 +10,10 @@ calc_dist_weights( line_lengths, radius_m, res_m, - method = "m3" + kernel = c("m3", "bb", "m1", "m5", "sq"), + gamma = 3e+10, + kl = 9e+05, + s1 = 200 ) } \arguments{ @@ -29,9 +32,31 @@ the weight.} determine the nugget added to the variance to represent the uncertainty in the starting and ending location of the transition.} -\item{method}{The method used for calculating the standard deviation -in the probability distribution. Currently \code{"m3"}, Martern 3/2; and -\code{"bb"}, Brownian bridge are supported.} +\item{kernel}{The kernel used for calculating the standard deviation +in the probability distribution: +\describe{ +\item{\code{"m3"}}{(default) Matern 3/2. Driven by \code{gamma} and \code{kl}.} +\item{\code{"m1"}}{Matern 1/2. Driven by \code{gamma} and \code{kl}.} +\item{\code{"m5"}}{Matern 5/2. Driven by \code{gamma} and \code{kl}.} +\item{\code{"sq"}}{Squared-exponential (Gaussian). Driven by \code{gamma} and \code{kl}.} +\item{\code{"bb"}}{Brownian bridge. Driven by \code{s1}.} +} +See \code{\link[=visualize_distance_weights]{visualize_distance_weights()}} to explore how these kernels and their +hyperparameters shape the spread.} + +\item{gamma}{Spread magnitude hyperparameter (m^2, a variance) for the +Matern-family and squared-exponential kernels (\code{"m1"}, \code{"m3"}, \code{"m5"}, +\code{"sq"}). Ignored for \code{kernel = "bb"}. The default was tuned by eye with +\code{visualize_distance_weights()}.} + +\item{kl}{Lengthscale hyperparameter (m) for the Matern-family and +squared-exponential kernels (\code{"m1"}, \code{"m3"}, \code{"m5"}, \code{"sq"}). Ignored for +\code{kernel = "bb"}. The default was tuned by eye using +\code{visualize_distance_weights()}.} + +\item{s1}{Spread magnitude hyperparameter (units of sqrt(m), not m) for +the \code{"bb"} kernel. Ignored for all other kernels. The default was tuned by +eye against \code{visualize_distance_weights()}.} } \value{ A vector of weights of the same length as the first three arguments. @@ -47,4 +72,7 @@ will be vectorized over the corresponding elements. This is a preliminary version of the function and will likely change. } +\seealso{ +\code{\link[=visualize_distance_weights]{visualize_distance_weights()}} +} \keyword{internal} diff --git a/man/calc_dist_weights_sd.Rd b/man/calc_dist_weights_sd.Rd new file mode 100644 index 00000000..ca69896d --- /dev/null +++ b/man/calc_dist_weights_sd.Rd @@ -0,0 +1,64 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calc_distance_weights.R +\name{calc_dist_weights_sd} +\alias{calc_dist_weights_sd} +\title{Standard deviation of the spread kernel used by \code{\link[=calc_dist_weights]{calc_dist_weights()}}} +\usage{ +calc_dist_weights_sd( + dist_along_line, + line_lengths, + res_m, + kernel, + gamma, + kl, + s1 +) +} +\arguments{ +\item{dist_along_line}{How far along the line the point of interest is (m)} + +\item{line_lengths}{How long the line is (m)} + +\item{res_m}{The resolution of the associated bird flow model, used to +determine the nugget added to the variance to represent the uncertainty in +the starting and ending location of the transition.} + +\item{kernel}{The kernel used for calculating the standard deviation +in the probability distribution: +\describe{ +\item{\code{"m3"}}{(default) Matern 3/2. Driven by \code{gamma} and \code{kl}.} +\item{\code{"m1"}}{Matern 1/2. Driven by \code{gamma} and \code{kl}.} +\item{\code{"m5"}}{Matern 5/2. Driven by \code{gamma} and \code{kl}.} +\item{\code{"sq"}}{Squared-exponential (Gaussian). Driven by \code{gamma} and \code{kl}.} +\item{\code{"bb"}}{Brownian bridge. Driven by \code{s1}.} +} +See \code{\link[=visualize_distance_weights]{visualize_distance_weights()}} to explore how these kernels and their +hyperparameters shape the spread.} + +\item{gamma}{Spread magnitude hyperparameter (m^2, a variance) for the +Matern-family and squared-exponential kernels (\code{"m1"}, \code{"m3"}, \code{"m5"}, +\code{"sq"}). Ignored for \code{kernel = "bb"}. The default was tuned by eye with +\code{visualize_distance_weights()}.} + +\item{kl}{Lengthscale hyperparameter (m) for the Matern-family and +squared-exponential kernels (\code{"m1"}, \code{"m3"}, \code{"m5"}, \code{"sq"}). Ignored for +\code{kernel = "bb"}. The default was tuned by eye using +\code{visualize_distance_weights()}.} + +\item{s1}{Spread magnitude hyperparameter (units of sqrt(m), not m) for +the \code{"bb"} kernel. Ignored for all other kernels. The default was tuned by +eye against \code{visualize_distance_weights()}.} +} +\value{ +A vector of standard deviations (m), the same length as +\code{dist_along_line}. +} +\description{ +Computes the standard deviation (in m) of the probability distribution +representing uncertainty in a bird's location at a given point along a +transition line, for one of the kernels supported by +\code{\link[=calc_dist_weights]{calc_dist_weights()}}. Factored out of \code{\link[=calc_dist_weights]{calc_dist_weights()}} so +\code{\link[=visualize_distance_weights]{visualize_distance_weights()}} can plot the same spread it uses without +duplicating the kernel-selection logic. +} +\keyword{internal} diff --git a/man/calc_euclidean_detection_rate.Rd b/man/calc_euclidean_detection_rate.Rd new file mode 100644 index 00000000..de990fdc --- /dev/null +++ b/man/calc_euclidean_detection_rate.Rd @@ -0,0 +1,98 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calc_euclidean_detection_rate.R +\name{calc_euclidean_detection_rate} +\alias{calc_euclidean_detection_rate} +\title{Assign weights to points based on how close they are to the planar +(Euclidean) line associated with each transition} +\usage{ +calc_euclidean_detection_rate( + bf, + points = NULL, + radius = NULL, + n_directions = 1, + skip_unconnected = TRUE, + batch_size = 1e+05, + check_radius = TRUE, + ... +) +} +\arguments{ +\item{bf}{A BirdFlow model} + +\item{points}{The points to evaluate betweenness on. If NULL the cell +centers of all the raster cells within the BirdFlow model that are between +active cells in the model will be used. This is calculated by comparing the +cell centers to a buffered convex hull around the active cell centers.} + +\item{radius}{The probability density along a transect perpendicular to +the line and intersecting the point is summed over the interval within +\code{radius} meters of the point. \code{radius} defaults to half the cell +size (\code{mean(res(bf))/2}).} + +\item{n_directions}{The number of (equally spaced) directional bins to +classify bearings into. Currently only \code{1} is supported.} + +\item{skip_unconnected}{If \code{TRUE} then only connections that exist in \code{bf} +will be evaluated, and between matrix will erroneously indicate that +the weights associated with transitions that aren't used is always 0. +The resulting array can still be used with the model it was built for because +those missing connections would always have zero probability.} + +\item{batch_size}{controls the number of movement lines that are processed +at a time. A smaller +\code{batch_size} will conserve memory at a slight performance cost. The number +of batches will be less than or equal to \code{n_active(bf)^2 / batch_size}.} + +\item{check_radius}{If \code{TRUE} an error will be thrown if the radius +is not between the resolution and 1/4 the resolution of \code{bf}. Outside of +that range the algorithm is likely to yield distorted results. +\code{0.5 * mean(res(bf))} is the default, and recommended radius.} + +\item{...}{Additional arguments forwarded to \code{\link[=calc_dist_weights]{calc_dist_weights()}}, +allowing the spread kernel (\code{kernel}) and its hyperparameters +(\code{gamma}, \code{kl}, \code{s1}) to be tuned.} +} +\value{ +A list with: +\item{between}{An array with dimensions representing the +"from" location, the "to" location, and the \code{points}. Cells are weights and +will be non-zero if the radius around the point intersects 1.96 standard +deviations of the normal distribution of probabilities around the line.} +\item{points}{A data,frame of points that define the third dimension +in \code{between}. It is identical to the input \code{points} if they are not \code{NULL}. +Otherwise it will be a data frame with columns \code{x}, \code{y}, and \code{i} +corresponding to the third dimension in \code{between}. +\code{i} will be \code{NA} for points that are not within the mask but +fall between active cells.} +\item{radius}{The radius of the circle in meters.} +} +\description{ +This internal function is used to create a betweenness weights array +with dimensions \href{n_active()}{n_active(bf)}, +\href{n_active()}{n_active(bf)}, and \code{length(points)}. +The first two dimensions represent from and to cells of possible connections +between pairs of locations within the BirdFlow model and both have +\code{n_active()} elements. +The third dimension represents reference points that might be between each +connection. +Cell values are the weight to use when adding the transition probabilities +to the reference point BMTR +} +\details{ +If \code{points} are \code{NULL} they default to the center of all cells within +the BirdFlow extent that fall between any active cells. +This includes all cell centers within a convex hull +(in the model's native, projected CRS) around the active cells in \code{bf} +and thus is almost always more than just the active cells. + +Unlike \code{\link[=calc_spherical_detection_rate]{calc_spherical_detection_rate()}}, this projects points onto the +line segment connecting each pair of active cells directly in the model's +native (planar) CRS, using vectorized matrix algebra +(\code{\link[=project_points_simple]{project_points_simple()}}) rather than great-circle math. This is +substantially faster and is the recommended continuous detection method. +} +\seealso{ +\code{\link[=is_between]{is_between()}}, \code{\link[=calc_spherical_detection_rate]{calc_spherical_detection_rate()}}, and +\code{\link[=calc_bmtr]{calc_bmtr()}} +} +\keyword{internal} diff --git a/man/calc_flux.Rd b/man/calc_flux.Rd index 22210081..96f1182b 100644 --- a/man/calc_flux.Rd +++ b/man/calc_flux.Rd @@ -15,10 +15,14 @@ calc_flux(...) \code{NULL} they will default to the BirdFlow model cells that are either active or fall between two active cells. Otherwise a data frame with \code{x} and \code{y} columns containing point coordinates in \link[terra:crs]{crs(bf)}.} - \item{\code{radius}}{The radius in meters around the points used to assess whether -a movement line passes by (or through) the point. If a point is farther than -\code{radius} from a great circle line between two cells centers then it is not -between them.} + \item{\code{radius}}{The radius in meters around the points used to assess the +detection rate for a movement at the point. +With \code{method = "binary"}, if a point is within \code{radius} of the great +circle line between two cell centers then the movement is detected at +that point. For the two continuous detection methods there is a +probability distribution for the location of the bird as it passes by the +point, and the radius defines the band over which that probability is +integrated, giving the detection rate.} \item{\code{n_directions}}{The number of directional bins to use for recording movement direction. Must be either \code{1} indicating no direction information or an even number. This is a placeholder, currently only \code{1} is supported.} @@ -40,13 +44,21 @@ between timesteps. \item{\code{"SpatRaster"}}{Returns a \code{terra::SpatRaster} with layers for each transition.} }} - \item{\code{weighted}}{If \code{FALSE} use the original and quicker version of bmtr -that sums all the marginal probability for transitions that pass within a -fixed distance of the point. If \code{TRUE} assign a weight to the point and -transition combo that then is multiplied by the marginal probability before -summing. This argument is experimental but the default value is identical -to the old version. The argument name and behavior when set to \code{TRUE} may -change.} + \item{\code{method}}{The detection model used to determine how much of a +transition's movement counts towards a point's BMTR: +\describe{ +\item{\code{"binary"}}{(default) Fast and deterministic. A movement line either +does or does not pass within \code{radius} of the point, per \code{\link[=is_between]{is_between()}}.} +\item{\code{"continuous"}}{Assigns a continuous weight (0 to 1) based on the +probability that a bird's actual path, modeled as spreading away from the +straight line between two cells, passes within \code{radius} of the point. Uses +planar (Euclidean) geometry in the model's native CRS, and is the +recommended detection model when continuous weighting is desired.} +\item{\code{"continuous-spherical"}}{The same continuous weighting as +\code{"continuous"}, but computed with great-circle (spherical) geometry +instead. Much slower, and not recommended for routine use; kept to allow +assessing the impact of switching from spherical to Euclidean geometry.} +}} \item{\code{batch_size}}{controls the number of movement lines that are processed at a time. A smaller \code{batch_size} will conserve memory at a slight performance cost. The number diff --git a/man/weight_between.Rd b/man/calc_spherical_detection_rate.Rd similarity index 71% rename from man/weight_between.Rd rename to man/calc_spherical_detection_rate.Rd index 94d11da6..187d85b9 100644 --- a/man/weight_between.Rd +++ b/man/calc_spherical_detection_rate.Rd @@ -1,18 +1,18 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/weight_between.R -\name{weight_between} -\alias{weight_between} +% Please edit documentation in R/calc_spherical_detection_rate.R +\name{calc_spherical_detection_rate} +\alias{calc_spherical_detection_rate} \title{Assign weights to points based on how close they are to the great circle line associated with each transition} \usage{ -weight_between( +calc_spherical_detection_rate( bf, - weight_fun = NULL, points = NULL, radius = NULL, n_directions = 1, skip_unconnected = TRUE, batch_size = 1e+05, + check_radius = TRUE, ... ) } @@ -42,6 +42,15 @@ those missing connections would always have zero probability.} at a time. A smaller \code{batch_size} will conserve memory at a slight performance cost. The number of batches will be less than or equal to \code{n_active(bf)^2 / batch_size}.} + +\item{check_radius}{If \code{TRUE} an error will be thrown if the radius +is not between the resolution and 1/4 the resolution of \code{bf}. Outside of +that range the algorithm is likely to yield distorted results. +\code{0.5 * mean(res(bf))} is the default, and recommended radius.} + +\item{...}{Additional arguments forwarded to \code{\link[=calc_dist_weights]{calc_dist_weights()}}, +allowing the spread kernel (\code{kernel}) and its hyperparameters +(\code{gamma}, \code{kl}, \code{s1}) to be tuned.} } \value{ A list with: @@ -77,14 +86,20 @@ This includes all cell centers within a convex hull (in spherical coordinates) around the active cells in \code{bf} and thus is almost always more than just the active cells. -\code{weight_between()} and \code{is_between()} should only differ slightly in -their results when calculated on the same model. +\code{calc_spherical_detection_rate()} and \code{is_between()} should only differ +slightly in their results when calculated on the same model. 1, The number of included reference points may differ. 2. The betweenness array will be real weights from 0 to 1 with -\code{weight_between()} and logical with \code{is_between()} +\code{calc_spherical_detection_rate()} and logical with \code{is_between()} + +This is the great-circle (spherical) counterpart to +\code{\link[=calc_euclidean_detection_rate]{calc_euclidean_detection_rate()}}. It's much slower and not recommended +for production use, but is kept to allow assessing the impact of switching +from spherical to Euclidean geometry. } \seealso{ -\code{\link[=is_between]{is_between()}} and \code{\link[=calc_bmtr]{calc_bmtr()}} +\code{\link[=is_between]{is_between()}}, \code{\link[=calc_euclidean_detection_rate]{calc_euclidean_detection_rate()}}, and +\code{\link[=calc_bmtr]{calc_bmtr()}} } \keyword{internal} diff --git a/man/plot_bmtr.Rd b/man/plot_bmtr.Rd index bcdfb91b..bf13200f 100644 --- a/man/plot_bmtr.Rd +++ b/man/plot_bmtr.Rd @@ -75,11 +75,29 @@ Plot BirdFlow Migration Traffic Rate (BMTR) \dontrun{ bf <- BirdFlowModels::amewoo -bmtr <- calc_bmtr(bf) +# Binary detection along shortest path +bmtr <- calc_bmtr(bf) plot_bmtr(bmtr, bf) animate_bmtr(bmtr, bf) + +# Continuous detection with a wider spread kernel +# (gamma defaults to 3e10) +bmtr2 <- calc_bmtr(bf, method = "continuous", gamma = 6e10) +animate_bmtr(bmtr2, bf) + +# Visualize spread for the continuous detection with wider spread +visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf), + gamma = 6e10, res_m = xres(bf)) + +# Visualize spread for the continuous detection with default values +visualize_distance_weights(line_lengths = c(2, 5, 10) * xres(bf), + res_m = xres(bf)) + } + + + } diff --git a/man/project_points_simple.Rd b/man/project_points_simple.Rd new file mode 100644 index 00000000..1c090c2d --- /dev/null +++ b/man/project_points_simple.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calc_euclidean_detection_rate.R +\name{project_points_simple} +\alias{project_points_simple} +\title{Project points onto line segments (Euclidean/planar geometry)} +\usage{ +project_points_simple(s, e, p, m, n, v, vv, clamp = FALSE, tol = 0) +} +\arguments{ +\item{s}{An n x 2 matrix of segment start coordinates.} + +\item{e}{An n x 2 matrix of segment end coordinates.} + +\item{p}{An m x 2 matrix of point coordinates.} + +\item{m}{Number of points (\code{nrow(p)}).} + +\item{n}{Number of line segments (\code{nrow(s)}).} + +\item{v}{\code{e - s}, precomputed segment vectors.} + +\item{vv}{\code{rowSums(v^2)}, precomputed squared segment lengths.} + +\item{clamp}{If \code{TRUE} constrain the projected fraction to \verb{[0, 1]}. If +\code{FALSE} (default) projections falling outside the segment are set to \code{NA}.} + +\item{tol}{Segments with squared length at or below this threshold are +treated as having zero length (projected fraction set to 0).} +} +\value{ +A list with \code{dist_to_line} and \code{dist_along}, each a vector of +length \code{m * n} (point index varying fastest). +} +\description{ +Vectorized computation of the distance from each point in \code{P} to each +line segment defined by \code{S} (starts) and \code{E} (ends), and the distance +along the segment to the projection of the point. +} +\keyword{internal} diff --git a/man/visualize_distance_weights.Rd b/man/visualize_distance_weights.Rd new file mode 100644 index 00000000..ab06db98 --- /dev/null +++ b/man/visualize_distance_weights.Rd @@ -0,0 +1,91 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/visualize_distance_weights.R +\name{visualize_distance_weights} +\alias{visualize_distance_weights} +\title{Visualize the spread kernel used by \code{\link[=calc_dist_weights]{calc_dist_weights()}}} +\usage{ +visualize_distance_weights( + line_lengths = c(2000, 1000, 500, 200), + radius_m = NULL, + res_m = 1000, + kernel = c("m3", "bb", "m1", "m5", "sq"), + gamma = 3e+10, + kl = 9e+05, + s1 = 200, + type = c("envelope", "raster"), + n = 100 +) +} +\arguments{ +\item{line_lengths}{A vector of transition-line lengths (m) to visualize +the spread kernel for.} + +\item{radius_m}{The detection radius (m). Defaults to \code{res_m / 2} (the +same default used by \code{\link[=calc_bmtr]{calc_bmtr()}}'s continuous methods). For +\code{type = "envelope"} this is only shown as a reference line, and only if +explicitly supplied. For \code{type = "raster"} it's used, as in +\code{\link[=calc_dist_weights]{calc_dist_weights()}}, to determine the band of probability density +that forms the plotted weight.} + +\item{res_m}{The resolution of the associated bird flow model (m), used +to determine the nugget added to the variance. See \code{\link[=calc_dist_weights]{calc_dist_weights()}}.} + +\item{kernel}{The kernel used for calculating the standard deviation +in the probability distribution: +\describe{ +\item{\code{"m3"}}{(default) Matern 3/2. Driven by \code{gamma} and \code{kl}.} +\item{\code{"m1"}}{Matern 1/2. Driven by \code{gamma} and \code{kl}.} +\item{\code{"m5"}}{Matern 5/2. Driven by \code{gamma} and \code{kl}.} +\item{\code{"sq"}}{Squared-exponential (Gaussian). Driven by \code{gamma} and \code{kl}.} +\item{\code{"bb"}}{Brownian bridge. Driven by \code{s1}.} +} +See \code{\link[=visualize_distance_weights]{visualize_distance_weights()}} to explore how these kernels and their +hyperparameters shape the spread.} + +\item{gamma}{Spread magnitude hyperparameter (m^2, a variance) for the +Matern-family and squared-exponential kernels (\code{"m1"}, \code{"m3"}, \code{"m5"}, +\code{"sq"}). Ignored for \code{kernel = "bb"}. The default was tuned by eye with +\code{visualize_distance_weights()}.} + +\item{kl}{Lengthscale hyperparameter (m) for the Matern-family and +squared-exponential kernels (\code{"m1"}, \code{"m3"}, \code{"m5"}, \code{"sq"}). Ignored for +\code{kernel = "bb"}. The default was tuned by eye using +\code{visualize_distance_weights()}.} + +\item{s1}{Spread magnitude hyperparameter (units of sqrt(m), not m) for +the \code{"bb"} kernel. Ignored for all other kernels. The default was tuned by +eye against \code{visualize_distance_weights()}.} + +\item{type}{The style of plot: +\describe{ +\item{\code{"envelope"}}{(default) For each line length, plots the region +within 1.96 standard deviations of the line (in the kernel's probability +distribution) as a band along the line. This represents 95 percent of +the density.} +\item{\code{"raster"}}{For each line length, plots the actual weight returned +by \code{\link[=calc_dist_weights]{calc_dist_weights()}} over a grid of points around the line, as a +heatmap. Slower, but reflects \code{radius_m}'s effect on the weight as well +as the kernel's spread.} +}} + +\item{n}{The number of grid points used along (and, for \code{type = "raster"}, +across) each line.} +} +\value{ +A \code{ggplot2} object. It can be displayed with \code{print()}. +} +\description{ +\code{visualize_distance_weights()} plots the shape of the spread kernel used +by \code{\link[=calc_dist_weights]{calc_dist_weights()}} (and, through it, the \code{"continuous"} and +\code{"continuous-spherical"} methods of \code{\link[=calc_bmtr]{calc_bmtr()}}) for one or more +transition-line lengths. It's meant as a tool for choosing \code{kernel} and +its hyperparameters (\code{gamma}, \code{kl}, \code{s1}) before running an expensive +\code{calc_bmtr()} call. +} +\examples{ +visualize_distance_weights() +visualize_distance_weights(type = "raster", kernel = "bb", s1 = 20) +} +\seealso{ +\code{\link[=calc_dist_weights]{calc_dist_weights()}}, \code{\link[=calc_bmtr]{calc_bmtr()}} +} diff --git a/tests/testthat/_snaps/calc_bmtr.md b/tests/testthat/_snaps/calc_bmtr.md index 89a9d9fc..64e48e99 100644 --- a/tests/testthat/_snaps/calc_bmtr.md +++ b/tests/testthat/_snaps/calc_bmtr.md @@ -11,3 +11,42 @@ 114 -1875000 -375000 T_02-03 1.173e-10 2021-01-14 149 -1725000 375000 T_01-02 2.287e-10 2021-01-07 +# calc_bmtr() works with method = 'binary' + + Code + top + Output + x y transition bmtr date + 65 -2025000 -225000 T_01-02 3.623e-09 2021-01-07 + 66 -2025000 -225000 T_02-03 1.897e-08 2021-01-14 + 109 -1875000 -225000 T_01-02 1.061e-08 2021-01-07 + 110 -1875000 -225000 T_02-03 8.883e-08 2021-01-14 + 114 -1875000 -375000 T_02-03 1.173e-10 2021-01-14 + 149 -1725000 375000 T_01-02 2.287e-10 2021-01-07 + +# calc_bmtr() works with method = 'continuous' + + Code + top + Output + x y transition bmtr date + 133 -2025000 -375000 T_01-02 5.203e-11 2021-01-07 + 134 -2025000 -375000 T_02-03 5.381e-10 2021-01-14 + 185 -1875000 -75000 T_01-02 1.924e-10 2021-01-07 + 186 -1875000 -75000 T_02-03 5.228e-10 2021-01-14 + 189 -1875000 -225000 T_01-02 2.775e-09 2021-01-07 + 190 -1875000 -225000 T_02-03 1.494e-08 2021-01-14 + +# calc_bmtr() works with method = 'continuous-spherical' + + Code + top + Output + x y transition bmtr date + 97 -2025000 -375000 T_01-02 9.882e-11 2021-01-07 + 98 -2025000 -375000 T_02-03 5.768e-10 2021-01-14 + 141 -1875000 -75000 T_01-02 1.969e-10 2021-01-07 + 142 -1875000 -75000 T_02-03 5.583e-10 2021-01-14 + 145 -1875000 -225000 T_01-02 2.756e-09 2021-01-07 + 146 -1875000 -225000 T_02-03 1.488e-08 2021-01-14 + diff --git a/tests/testthat/_snaps/calc_euclidean_detection_rate.md b/tests/testthat/_snaps/calc_euclidean_detection_rate.md new file mode 100644 index 00000000..70641b1b --- /dev/null +++ b/tests/testthat/_snaps/calc_euclidean_detection_rate.md @@ -0,0 +1,14 @@ +# calc_euclidean_detection_rate() works + + Code + sum(between[[1]]) + Output + [1] 8683.589 + +# calc_euclidean_detection_rate() works with kernel = 'bb' + + Code + sum(between[[1]]) + Output + [1] 8925.382 + diff --git a/tests/testthat/_snaps/calc_spherical_detection_rate.md b/tests/testthat/_snaps/calc_spherical_detection_rate.md new file mode 100644 index 00000000..a966d9da --- /dev/null +++ b/tests/testthat/_snaps/calc_spherical_detection_rate.md @@ -0,0 +1,14 @@ +# calc_spherical_detection_rate() works + + Code + sum(between[[1]]) + Output + [1] 8684.696 + +# calc_spherical_detection_rate() works with kernel = 'bb' + + Code + sum(between[[1]]) + Output + [1] 8931.652 + diff --git a/tests/testthat/helper-small_bf.R b/tests/testthat/helper-small_bf.R new file mode 100644 index 00000000..79ecdaf6 --- /dev/null +++ b/tests/testthat/helper-small_bf.R @@ -0,0 +1,13 @@ + +#' A small, fast BirdFlow model for use in tests +#' +#' Truncates and sparsifies `BirdFlowModels::amewoo` down to a fixture +#' that's fast enough to use repeatedly across the test suite. +#' +#' @return A small `BirdFlow` model. +small_test_bf <- function() { + bf <- BirdFlowModels::amewoo + bf <- truncate_birdflow(bf, start = 1, end = 5) + bf <- sparsify(bf, "conditional", .9, p_protected = 0.05) + bf +} diff --git a/tests/testthat/test-calc_bmtr.R b/tests/testthat/test-calc_bmtr.R index f2f03a31..88006099 100644 --- a/tests/testthat/test-calc_bmtr.R +++ b/tests/testthat/test-calc_bmtr.R @@ -1,9 +1,6 @@ test_that("calc_bmtr() works without directionality", { local_quiet() - # Sparsify and truncate to speed things up - bf <- BirdFlowModels::amewoo - bf <- truncate_birdflow(bf, start = 1, end = 5) - bf <- sparsify(bf, "conditional", .9, p_protected = 0.05) + bf <- small_test_bf() expect_no_error(f <- calc_bmtr(bf)) @@ -18,32 +15,126 @@ test_that("calc_bmtr() works without directionality", { }) -test_that("calc_bmtr() works with weights", { +test_that("calc_bmtr() works with method = 'binary'", { local_quiet() - # Sparsify and truncate to speed things up - bf <- BirdFlowModels::amewoo - bf <- truncate_birdflow(bf, start = 1, end = 5) - bf <- sparsify(bf, "conditional", .9, p_protected = 0.05) + bf <- small_test_bf() - expect_no_error(f <- calc_bmtr(bf, weighted = TRUE)) + expect_no_error(f <- calc_bmtr(bf, method = "binary")) - # Snapshot of first 6 non-zero movements top <- head(f[!f$bmtr == 0, ], 6) top$bmtr <- signif(top$bmtr, 4) - # --- expect_snapshot(top) --- wait for final parameters to save snapshot + expect_snapshot(top) + + expect_no_error(plot_bmtr(f, bf)) + expect_no_error(animate_bmtr(f, bf)) +}) + +test_that("calc_bmtr() works with method = 'continuous'", { + local_quiet() + bf <- small_test_bf() + + expect_no_error(f <- calc_bmtr(bf, method = "continuous")) + + top <- head(f[!f$bmtr == 0, ], 6) + top$bmtr <- signif(top$bmtr, 4) + expect_snapshot(top) - # Visualizations expect_no_error(plot_bmtr(f, bf)) expect_no_error(animate_bmtr(f, bf)) +}) + +test_that("calc_bmtr() works with method = 'continuous-spherical'", { + local_quiet() + bf <- small_test_bf() + + expect_no_error(f <- calc_bmtr(bf, method = "continuous-spherical")) + + top <- head(f[!f$bmtr == 0, ], 6) + top$bmtr <- signif(top$bmtr, 4) + expect_snapshot(top) + + expect_no_error(plot_bmtr(f, bf)) + expect_no_error(animate_bmtr(f, bf)) +}) +test_that("calc_bmtr() errors on an invalid method", { + local_quiet() + bf <- small_test_bf() + + expect_error(calc_bmtr(bf, method = "bogus")) +}) + +test_that("calc_bmtr() respects a user-supplied points argument (Bug A)", { + local_quiet() + bf <- small_test_bf() + + pts <- i_to_xy(seq_len(min(5, n_active(bf))), bf)[, c("x", "y")] + + for (method in c("continuous", "continuous-spherical")) { + f <- calc_bmtr(bf, points = pts, method = method, format = "points") + expect_equal(nrow(f$point), nrow(pts)) + expect_equal(sort(f$point$x), sort(pts$x)) + expect_equal(sort(f$point$y), sort(pts$y)) + } +}) + +test_that("calc_bmtr() respects a user-supplied radius (Bug B)", { + local_quiet() + bf <- small_test_bf() + + r1 <- mean(res(bf)) * 0.3 + r2 <- mean(res(bf)) * 0.9 + + for (method in c("continuous", "continuous-spherical")) { + f1 <- calc_bmtr(bf, method = method, radius = r1) + f2 <- calc_bmtr(bf, method = method, radius = r2) + expect_false(isTRUE(all.equal(f1$bmtr, f2$bmtr))) + } +}) + +test_that("calc_bmtr() enforces check_radius bounds for continuous methods", { + local_quiet() + bf <- small_test_bf() + + bad_radius <- mean(res(bf)) * 2 + + for (method in c("continuous", "continuous-spherical")) { + expect_error(calc_bmtr(bf, method = method, radius = bad_radius)) + expect_no_error(calc_bmtr(bf, method = method, radius = bad_radius, + check_radius = FALSE)) + } +}) + +test_that("'continuous' and 'continuous-spherical' agree closely", { + local_quiet() + bf <- small_test_bf() + + euc <- calc_bmtr(bf, method = "continuous") + sph <- calc_bmtr(bf, method = "continuous-spherical") + + euc$id <- paste(euc$x, euc$y, euc$transition) + sph$id <- paste(sph$x, sph$y, sph$transition) + m <- merge(euc, sph, by = "id", suffixes = c("_euc", "_sph")) + + expect_true(nrow(m) > 0) + expect_equal(m$bmtr_euc, m$bmtr_sph, tolerance = 0.05) +}) + +test_that("calc_bmtr() forwards ... to the spread kernel", { + local_quiet() + bf <- small_test_bf() + + default <- calc_bmtr(bf, method = "continuous") + wider <- calc_bmtr(bf, method = "continuous", gamma = 100000) + expect_false(isTRUE(all.equal(default$bmtr, wider$bmtr))) + + expect_error(calc_bmtr(bf, method = "binary", gamma = 100000)) }) test_that("plot_bmtr() subset drops excluded transitions (no empty facets)", { local_quiet() - bf <- BirdFlowModels::amewoo - bf <- truncate_birdflow(bf, start = 1, end = 5) - bf <- sparsify(bf, "conditional", .9, p_protected = 0.05) + bf <- small_test_bf() f <- calc_bmtr(bf) all_trans <- sort(unique(f$transition)) @@ -177,7 +268,7 @@ test_that("Test sensativity of bmtr to radius", { # and the relationship between radius and abundance that intersects the # circle approaches the square of the radius. # 3. (2) begins to break down as the diameter approaches the extent width, - # at that point you've already captured most of the movent at every point + # at that point you've already captured most of the movement at every point # and a larger radius just means you divide that total by a bigger number. # # Two conclusions: diff --git a/tests/testthat/test-calc_distance_weights.R b/tests/testthat/test-calc_distance_weights.R new file mode 100644 index 00000000..50dd59e2 --- /dev/null +++ b/tests/testthat/test-calc_distance_weights.R @@ -0,0 +1,92 @@ +test_that("calc_dist_weights() works with the default 'm3' kernel", { + w <- calc_dist_weights(dist_to_line = c(0, 500, 50000), + dist_along_line = c(500, 500, 500), + line_lengths = c(1000, 1000, 1000), + radius_m = 200, res_m = 1000) + expect_length(w, 3) + expect_true(all(w >= 0 & w <= 1)) + # Closer points get more weight + expect_true(w[1] > w[2]) + expect_equal(w[3], 0) +}) + +test_that("calc_dist_weights() works with the 'bb' kernel", { + w <- calc_dist_weights(dist_to_line = c(0, 500, 50000), + dist_along_line = c(500, 500, 500), + line_lengths = c(1000, 1000, 1000), + radius_m = 200, res_m = 1000, kernel = "bb") + expect_length(w, 3) + expect_true(all(w >= 0 & w <= 1)) + expect_true(w[1] > w[2]) + expect_equal(w[3], 0) +}) + +test_that("calc_dist_weights() works with the 'm1', 'm5', and 'sq' kernels", { + # Regression test: these previously errored with + # "Error in get(kernel) : object 'm1' not found" because the kernel + # functions are named k_m1/k_m5/k_sq, not m1/m5/sq. + for (k in c("m1", "m5", "sq")) { + w <- calc_dist_weights(dist_to_line = c(0, 500, 50000), + dist_along_line = c(500, 500, 500), + line_lengths = c(1000, 1000, 1000), + radius_m = 200, res_m = 1000, kernel = k) + expect_length(w, 3) + expect_true(all(w >= 0 & w <= 1)) + expect_true(w[1] > w[2]) + expect_equal(w[3], 0) + } +}) + +test_that("calc_dist_weights() errors on an invalid kernel", { + expect_error(calc_dist_weights(dist_to_line = 0, dist_along_line = 500, + line_lengths = 1000, radius_m = 200, + res_m = 1000, kernel = "bogus")) +}) + +test_that("calc_dist_weights() gamma/kl affect the 'm3' kernel's spread", { + args <- list(dist_to_line = 300, dist_along_line = 500, + line_lengths = 1000, radius_m = 200, res_m = 1000) + + default <- do.call(calc_dist_weights, args) + wider_gamma <- do.call(calc_dist_weights, c(args, list(gamma = 400000))) + wider_kl <- do.call(calc_dist_weights, c(args, list(kl = 50))) + + expect_false(isTRUE(all.equal(default, wider_gamma))) + expect_false(isTRUE(all.equal(default, wider_kl))) +}) + +test_that("calc_dist_weights_sd() keeps t/len and gamma/kl in consistent units", { + # Regression test for a units bug where dist_along_line/line_lengths were + # converted to km while gamma/kl were left unconverted (meters), making the + # kernel decay ~1000x slower than intended and its contribution to sd + # negligible next to the nugget. With consistent units, once the line is + # long relative to kl the bridge fully decorrelates mid-line and sd should + # approach sqrt(gamma + nugget^2), not just the nugget. + res_m <- 1000 + gamma <- 40000 + kl <- 2000 + nugget <- res_m / 4 + len <- 20 * kl # long enough for the kernel to fully decay mid-line + + sd_end <- calc_dist_weights_sd(0, len, res_m, kernel = "m3", + gamma = gamma, kl = kl, s1 = NA) + sd_mid <- calc_dist_weights_sd(len / 2, len, res_m, kernel = "m3", + gamma = gamma, kl = kl, s1 = NA) + + expect_equal(sd_end, nugget) + expect_equal(sd_mid, sqrt(gamma + nugget^2), tolerance = 1e-3) +}) + +test_that("calc_dist_weights() s1 widens the 'bb' kernel's spread", { + # At a distance well beyond the default kernel's reach, a much larger s1 + # (wider spread) should pick up nonzero weight where the default has none. + args <- list(dist_to_line = 20000, dist_along_line = 500, + line_lengths = 1000, radius_m = 200, res_m = 1000, + kernel = "bb") + + default <- do.call(calc_dist_weights, args) + wider <- do.call(calc_dist_weights, c(args, list(s1 = 1000))) + + expect_equal(default, 0) + expect_true(wider > default) +}) diff --git a/tests/testthat/test-calc_euclidean_detection_rate.R b/tests/testthat/test-calc_euclidean_detection_rate.R new file mode 100644 index 00000000..8f0a9edc --- /dev/null +++ b/tests/testthat/test-calc_euclidean_detection_rate.R @@ -0,0 +1,24 @@ +test_that("calc_euclidean_detection_rate() works", { + local_quiet() + bf <- small_test_bf() + + expect_no_error(between <- calc_euclidean_detection_rate(bf)) + expect_snapshot(sum(between[[1]])) +}) + +test_that("calc_euclidean_detection_rate() works with kernel = 'bb'", { + local_quiet() + bf <- small_test_bf() + + expect_no_error(between <- calc_euclidean_detection_rate(bf, kernel = "bb")) + expect_snapshot(sum(between[[1]])) +}) + +test_that("calc_euclidean_detection_rate() forwards spread overrides", { + local_quiet() + bf <- small_test_bf() + + default <- calc_euclidean_detection_rate(bf) + wider <- calc_euclidean_detection_rate(bf, gamma = 100000, kl = 10000) + expect_false(isTRUE(all.equal(sum(default[[1]]), sum(wider[[1]])))) +}) diff --git a/tests/testthat/test-calc_spherical_detection_rate.R b/tests/testthat/test-calc_spherical_detection_rate.R new file mode 100644 index 00000000..a1988848 --- /dev/null +++ b/tests/testthat/test-calc_spherical_detection_rate.R @@ -0,0 +1,24 @@ +test_that("calc_spherical_detection_rate() works", { + local_quiet() + bf <- small_test_bf() + + expect_no_error(between <- calc_spherical_detection_rate(bf)) + expect_snapshot(sum(between[[1]])) +}) + +test_that("calc_spherical_detection_rate() works with kernel = 'bb'", { + local_quiet() + bf <- small_test_bf() + + expect_no_error(between <- calc_spherical_detection_rate(bf, kernel = "bb")) + expect_snapshot(sum(between[[1]])) +}) + +test_that("calc_spherical_detection_rate() forwards spread overrides", { + local_quiet() + bf <- small_test_bf() + + default <- calc_spherical_detection_rate(bf) + wider <- calc_spherical_detection_rate(bf, gamma = 100000, kl = 10000) + expect_false(isTRUE(all.equal(sum(default[[1]]), sum(wider[[1]])))) +}) diff --git a/tests/testthat/test-is_between.R b/tests/testthat/test-is_between.R index fe60cfd2..19270fc3 100644 --- a/tests/testthat/test-is_between.R +++ b/tests/testthat/test-is_between.R @@ -2,10 +2,7 @@ test_that("is_between() works", { local_quiet() - # Sparsifying and truncating to speed things up - bf <- BirdFlowModels::amewoo - bf <- truncate_birdflow(bf, start = 1, end = 5) - bf <- sparsify(bf, "conditional", .9, p_protected = 0.05) + bf <- small_test_bf() between <- is_between(bf) diff --git a/tests/testthat/test-rasterize_bmtr.R b/tests/testthat/test-rasterize_bmtr.R index a999e7e3..b523e830 100644 --- a/tests/testthat/test-rasterize_bmtr.R +++ b/tests/testthat/test-rasterize_bmtr.R @@ -5,7 +5,7 @@ test_that("rasterize_bmtr() returns a valid, plottable raster with the same exte bf <- BirdFlowModels::amewoo # Unweighted BMTR object - bmtr_uw <- calc_bmtr(bf, weighted = FALSE) + bmtr_uw <- calc_bmtr(bf, method = "binary") raster_uw <- rasterize_bmtr(bmtr_uw, bf) @@ -25,7 +25,7 @@ test_that("rasterization works for weighted BMTR calculation", { bf <- BirdFlowModels::amewoo # Weighted, non-Euclidean BMTR object - bmtr_sph <- calc_bmtr(bf, weighted = TRUE) + bmtr_sph <- calc_bmtr(bf, method = "continuous-spherical") expect_no_error(raster_sph <- rasterize_bmtr(bmtr_sph, bf)) expect_no_error(terra::plot(raster_sph)) diff --git a/tests/testthat/test-visualize_distance_weights.R b/tests/testthat/test-visualize_distance_weights.R new file mode 100644 index 00000000..b2558883 --- /dev/null +++ b/tests/testthat/test-visualize_distance_weights.R @@ -0,0 +1,28 @@ +test_that("visualize_distance_weights() works with type = 'envelope'", { + p <- NULL + expect_no_error(p <- visualize_distance_weights(type = "envelope")) + expect_s3_class(p, "ggplot") + + p_bb <- NULL + expect_no_error( + p_bb <- visualize_distance_weights(type = "envelope", kernel = "bb", + s1 = 20)) + expect_s3_class(p_bb, "ggplot") +}) + +test_that("visualize_distance_weights() works with type = 'raster'", { + p <- NULL + expect_no_error(p <- visualize_distance_weights(type = "raster", n = 10)) + expect_s3_class(p, "ggplot") + + p_bb <- NULL + expect_no_error( + p_bb <- visualize_distance_weights(type = "raster", kernel = "bb", + s1 = 20, n = 10)) + expect_s3_class(p_bb, "ggplot") +}) + +test_that("visualize_distance_weights() errors on invalid kernel or type", { + expect_error(visualize_distance_weights(kernel = "bogus")) + expect_error(visualize_distance_weights(type = "bogus")) +}) diff --git a/tests/testthat/test-weight_between.R b/tests/testthat/test-weight_between.R deleted file mode 100644 index 6b211915..00000000 --- a/tests/testthat/test-weight_between.R +++ /dev/null @@ -1,13 +0,0 @@ -test_that("weight_between() works", { - - local_quiet() - - # Sparsifying and truncating to speed things up - bf <- BirdFlowModels::amewoo - bf <- truncate_birdflow(bf, start = 1, end = 5) - bf <- sparsify(bf, "conditional", .9, p_protected = 0.05) - - expect_no_error(bw <- weight_between(bf)) - - -})