diff --git a/.Rbuildignore b/.Rbuildignore index cacd3cb..b9c03b9 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -12,3 +12,4 @@ ^README\.Rmd$ ^GeoTwilightR\.Rproj$ ^codemeta\.json$ +^LICENSE$ diff --git a/CITATION.cff b/CITATION.cff index ff32d98..ca94d0f 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -8,7 +8,7 @@ message: 'To cite package "GeoPathSampleR" in publications use:' type: software license: GPL-3.0-or-later title: 'GeoPathSampleR: Posterior Sampling of Geolocator Paths' -version: 0.1.0.9000 +version: 0.2.0 doi: 10.5281/zenodo.21944194 abstract: Reconstructs bird migration trajectories from stationary-period likelihood maps built with GeoPressureR. A Gibbs sampler combines residence-dependent departure, @@ -89,8 +89,8 @@ references: given-names: Mathieu email: mathieu@mgravey.com orcid: https://orcid.org/0000-0002-0871-1507 - year: '2026' version: '>= 3.6.0' + year: '2026' - type: software title: ggplot2 abstract: 'ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics' diff --git a/DESCRIPTION b/DESCRIPTION index e476435..f1687e7 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: GeoPathSampleR Title: Posterior Sampling of Geolocator Paths -Version: 0.1.0.9000 +Version: 0.2.0 Authors@R: person("Raphaël", "Nussbaumer", , "rafnuss@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-8185-1020")) diff --git a/NEWS.md b/NEWS.md index 2b9f3d8..d142af5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,12 +2,20 @@ ## Main +- [Add boundary sweeps to resample movement timing while retaining residence-segment locations](https://github.com/GeoPressure/GeoPathSampleR/commit/f15e91e), enabled every iteration by default with `boundary_interval = 1`. Use `0` or `Inf` to disable them. + - [Update the default ground-speed and residence-dependent departure parameters](https://github.com/GeoPressure/GeoPathSampleR/commit/d40999d) from the pooled multi-sensor calibration. Results sampled with the previous defaults should be recomputed before comparison. - [Replace the route-detour model with a Gamma prior on route excess](https://github.com/GeoPressure/GeoPathSampleR/commit/7d4fb3f) whose scale depends on interval duration and direct endpoint distance ([3798978](https://github.com/GeoPressure/GeoPathSampleR/commit/3798978)). The prior applies to intervals with at least two movement days and endpoints at least 300 km apart, with predictors limited to the empirical calibration support. - Export `sampling_path_default_movement()` and document the movement model, `route_detour` and the modular cut of `long_period_light_only` in `sampling_path()`. ## Minor +- Exclude the redundant repository GPL text from the R source package to resolve its package-check note. + +- [Fix known locations on water cells, stabilise the prior for direct routes, and reject excluded stationary periods](https://github.com/GeoPressure/GeoPathSampleR/commit/49e7735). +- [Speed up route-prior evaluation, path summaries and diagnostics](https://github.com/GeoPressure/GeoPathSampleR/commit/6da2d15), and [avoid replicating scalar movement endpoints](https://github.com/GeoPressure/GeoPathSampleR/commit/0b989d1). +- Document boundary sweeps, the complete sampler update schedule, route-prior support and computational controls in the function reference and getting-started vignette. + - [Fix the diagnostic report rendering and browser opening](https://github.com/GeoPressure/GeoPathSampleR/commit/fb0d61f), and [declare its optional report dependencies](https://github.com/GeoPressure/GeoPathSampleR/commit/3ff7a89). - [Align the README and vignette with the movement model and sampler controls](https://github.com/GeoPressure/GeoPathSampleR/commit/ae3f0de), and add a provisional citation for the submitted methods manuscript alongside the software citation. - [Format the elapsed time compactly in the sampling message](https://github.com/GeoPressure/GeoPathSampleR/commit/38fe161). diff --git a/R/path_collapse.R b/R/path_collapse.R index fc00ae6..a3ef9b2 100644 --- a/R/path_collapse.R +++ b/R/path_collapse.R @@ -11,7 +11,7 @@ #' @param stap A stationary-period data frame containing `stap_id`, `start`, #' and `end`. If `NULL`, these columns are taken from `paths`. #' @param label_stap_id logical; if `TRUE`, replace `stap_id` in the output -#' with a character label of the form `"{stap_id_start}->{stap_id_end}`. +#' with a character label of the form `"{stap_id_start}->{stap_id_end}"`. #' The numeric `stap_id_start` and `stap_id_end` columns are always kept. #' #' @return A data frame with one row per collapsed stay and columns diff --git a/R/path_summary.R b/R/path_summary.R index 9d075b6..e87353f 100644 --- a/R/path_summary.R +++ b/R/path_summary.R @@ -80,23 +80,31 @@ path_summary <- function( } if (by == "stap") { + row_by_stap <- split( + seq_len(nrow(paths)), + factor(paths$stap_id, levels = stap_id) + ) out <- do.call( rbind, - lapply(stap_id, function(id) { - df <- paths[paths$stap_id == id, , drop = FALSE] - pos <- position_summary(df) - data.frame( - stap_id = id, - stap_id_start = id, - stap_id_end = id, - ind = as.integer(pos[["ind"]]), - lat = pos[["lat"]], - lon = pos[["lon"]], - start = stap$start[match(id, stap$stap_id)], - end = stap$end[match(id, stap$stap_id)], - n_path = length(unique(path_id[paths$stap_id == id])) - ) - }) + Map( + function(id, row_i) { + df <- paths[row_i, , drop = FALSE] + pos <- position_summary(df) + data.frame( + stap_id = id, + stap_id_start = id, + stap_id_end = id, + ind = as.integer(pos[["ind"]]), + lat = pos[["lat"]], + lon = pos[["lon"]], + start = stap$start[match(id, stap$stap_id)], + end = stap$end[match(id, stap$stap_id)], + n_path = length(unique(path_id[row_i])) + ) + }, + stap_id, + row_by_stap + ) ) row.names(out) <- NULL return(out) @@ -114,11 +122,16 @@ path_summary <- function( group_start <- match(seq_len(max(group)), group) group_end <- c(group_start[-1L] - 1L, n_stap) + row_by_group <- split( + seq_len(nrow(paths)), + group[match(paths$stap_id, stap_id)] + ) + out <- do.call( rbind, lapply(seq_along(group_start), function(i) { ids <- stap_id[group_start[i]:group_end[i]] - df <- paths[paths$stap_id %in% ids, , drop = FALSE] + df <- paths[row_by_group[[i]], , drop = FALSE] pos <- position_summary(df) data.frame( stap_id = ids[1L], diff --git a/R/sampling_path.R b/R/sampling_path.R index ebd3872..6df5e1b 100644 --- a/R/sampling_path.R +++ b/R/sampling_path.R @@ -75,6 +75,7 @@ sampling_path_with_route_model <- function( warmup = floor(iter / 4), thin = 1, block_interval = 2, + boundary_interval = 1, thr_likelihood = 0.99, thr_gs = 2000 / 24, refresh = 10, @@ -118,6 +119,11 @@ sampling_path_with_route_model <- function( block_interval >= 2, is.infinite(block_interval) || block_interval == floor(block_interval), + is.numeric(boundary_interval), + length(boundary_interval) == 1, + boundary_interval >= 0, + is.infinite(boundary_interval) || + boundary_interval == floor(boundary_interval), iter > warmup, is.null(seed) || (is.numeric(seed) && length(seed) == 1 && seed == floor(seed)), @@ -141,8 +147,21 @@ sampling_path_with_route_model <- function( } else { block_interval <- as.integer(block_interval) } + boundary_interval <- if (is.infinite(boundary_interval)) { + 0L + } else { + as.integer(boundary_interval) + } workers <- min(as.integer(workers), chains) + # Each transition is one daily step, so excluded periods cannot be bridged. + if (!all(tag$stap$include %||% TRUE)) { + cli::cli_abort(c( + "Excluded stationary periods are not supported: {.field stap$include} is {.val FALSE} for stap_id {.val {tag$stap$stap_id[!tag$stap$include]}}.", + "i" = "Re-run {.fun GeoPressureR::tag_set_map} without {.arg include_stap_id} or {.arg include_min_duration}." + )) + } + prep <- sampling_path_prepare( tag = tag, likelihood = likelihood, @@ -195,7 +214,6 @@ sampling_path_with_route_model <- function( use_parallel <- workers > 1L && chains > 1L && .Platform$OS.type != "windows" progress_bar <- NULL progress_closed <- FALSE - progress_pos <- 0L progress_fn <- NULL if (!quiet) { @@ -221,8 +239,10 @@ sampling_path_with_route_model <- function( add = TRUE ) progress_fn <- function(amount) { - progress_pos <<- progress_pos + amount - utils::setTxtProgressBar(progress_bar, progress_pos) + utils::setTxtProgressBar( + progress_bar, + utils::getTxtProgressBar(progress_bar) + amount + ) } } @@ -240,6 +260,7 @@ sampling_path_with_route_model <- function( warmup = warmup, thin = thin, block_interval = block_interval, + boundary_interval = boundary_interval, long_period_i = long_period_i, component_weights = component_weights, route_prior = prepared_route_prior, @@ -259,6 +280,7 @@ sampling_path_with_route_model <- function( warmup = warmup, thin = thin, block_interval = block_interval, + boundary_interval = boundary_interval, long_period_i = long_period_i, component_weights = component_weights, route_prior = prepared_route_prior, @@ -369,6 +391,7 @@ sampling_path_with_route_model <- function( thin = thin, refresh = refresh, block_interval = block_interval, + boundary_interval = boundary_interval, long_period_light_only = long_period_light_only, n_light_only_long_periods = length(long_period_i), n_repair_iterations = n_repair_iterations, @@ -423,7 +446,8 @@ sampling_path_with_route_model <- function( #' #' @param tag GeoPressureR tag object containing likelihood maps and #' `tag_set_map` parameters. Long periods are identified by `tag$stap$stap0`, -#' as created by [GeoPressureR::tag_stap_daily()]. +#' as created by [GeoPressureR::tag_stap_daily()]. All stationary periods must +#' be included; excluded periods cannot be bridged by daily transitions. #' @param iter number of Gibbs iterations per chain, including `warmup`. #' @param likelihood tag field containing the light likelihood map. #' @param movement daily movement model: a ground-speed kernel passed to @@ -452,6 +476,12 @@ sampling_path_with_route_model <- function( #' current residence segment as a single block. All other iterations are #' site iterations that update each eligible period once. Use `0` or `Inf` to #' disable block iterations; values must otherwise be at least `2`. +#' @param boundary_interval interval between boundary sweeps, which resample +#' the date of every boundary between consecutive residence segments while +#' keeping both locations. A sweep runs after the site or block update of +#' every `boundary_interval`-th iteration. The default `1` sweeps every +#' iteration; use `0` or `Inf` to disable them. Otherwise use a positive +#' integer. See Details for the conditional update. #' @param thr_likelihood retained light-likelihood mass in each stationary #' period. Lower values keep fewer cells and speed up the sampler, but exclude #' more of the light likelihood surface entirely. The default `0.99` is a @@ -486,16 +516,18 @@ sampling_path_with_route_model <- function( #' where \eqn{T_\ell} is the interval duration in days and \eqn{D_\ell} is in #' km. Both covariates are clamped to the range of the calibration data. The #' term is omitted when \eqn{D_\ell < 300} km or when the interval contains at -#' most one daily displacement. +#' most one non-zero daily displacement. #' #' Because `route_detour` multiplies the Gamma scale, the prior mean and every #' prior quantile of the route excess \eqn{R_\ell - 1} scale linearly with it. #' For example, an interval of about 31 days with endpoints about 3,500 km #' apart has \eqn{\theta_\ell \approx 0.15}, i.e. a prior mean route ratio of #' about 1.19; `route_detour = 2` raises it to about 1.37 and `0.5` lowers it -#' to about 1.09. `route_detour = 0` collapses the scale, so any detour is -#' strongly penalised. The prior remains soft: the realised detour of the -#' posterior also depends on the light and movement components. +#' to about 1.09. `route_detour = 0` uses machine epsilon as the scale, so any +#' detour is strongly penalised. Route excess is also floored at machine +#' epsilon when evaluating the log density. The prior remains soft: the +#' realised detour of the posterior also depends on the light and movement +#' components. #' #' **Long periods and modular cut.** Twilights accurately locate long periods #' but leave daily-period locations ambiguous. With @@ -510,6 +542,27 @@ sampling_path_with_route_model <- function( #' the movement model cannot pull long-period locations away from their light #' likelihood. Known-location periods are always fixed. #' +#' **Site, block and boundary updates.** By default, odd iterations update +#' each eligible period separately and even iterations update each current +#' residence segment at one shared location (`block_interval = 2`). Site +#' sweeps alternate forwards and backwards. Block updates use the intersection +#' of the periods' retained light supports and cannot span a long period when +#' the route prior is active or include a long period drawn from light alone. +#' +#' After either update, `boundary_interval = 1` adds a forward sweep over the +#' boundaries between consecutive residence segments. Each boundary is drawn +#' over all positions within its two segments that keep both non-empty, +#' retain their two locations, respect each period's light support, and keep +#' fixed periods at their current location. This changes departure and arrival +#' timing without changing the number or sequence of residence segments. +#' Candidate positions are weighted by the light likelihood and +#' residence-dependent departure terms, with `component_weights` applied. +#' The transition between the two locations is unchanged, and fixed +#' long-period endpoints keep cumulative distance within each route interval +#' unchanged, so the movement kernel and route prior cancel in this update. +#' Boundary sweeps run during warmup as well as sampling, before draws are +#' saved. Their frequency changes the update schedule, not the target model. +#' #' **Computational approximations.** Impossible geography, such as water after #' masking, is removed from the state space. The arguments `thr_likelihood`, #' `thr_gs` and `movement$approx_eps` control how aggressively the sampler @@ -578,6 +631,7 @@ sampling_path <- function( warmup = floor(iter / 4), thin = 1, block_interval = 2, + boundary_interval = 1, thr_likelihood = 0.99, thr_gs = 2000 / 24, refresh = 10, @@ -597,6 +651,7 @@ sampling_path <- function( warmup = warmup, thin = thin, block_interval = block_interval, + boundary_interval = boundary_interval, thr_likelihood = thr_likelihood, thr_gs = thr_gs, refresh = refresh, @@ -645,6 +700,7 @@ sampling_path_run_chain <- function( warmup, thin, block_interval, + boundary_interval, long_period_i, component_weights, route_prior, @@ -675,6 +731,8 @@ sampling_path_run_chain <- function( )) } + boundary_fixed <- seq_len(nstap) %in% block_excluded_i + # Initialize the path so every adjacent stap is movement-feasible before the # stochastic iterations start. path_idx <- sampling_path_initialize( @@ -819,6 +877,30 @@ sampling_path_run_chain <- function( } } + # Boundary sweeps resample when segments change, not where they are. + if ( + boundary_interval > 0L && + iter_i %% boundary_interval == 0L && + nstap > 1L + ) { + path_idx <- sampling_path_update_boundary( + path_idx = path_idx, + fixed = boundary_fixed, + lk = lk, + kt = kt, + component_weights = component_weights + ) + if (!is.null(route_edge_distance)) { + route_distance <- sampling_path_route_distance_state( + path_idx, + kt, + route_prior + ) + route_edge_distance <- route_distance$edge + route_interval_distance <- route_distance$interval + } + } + if (iter_i > warmup && ((iter_i - warmup) %% thin == 0L)) { save_idx <- save_idx + 1L samples[save_idx, ] <- path_idx @@ -932,9 +1014,10 @@ sampling_path_update_block <- function( ) } - state_next$log_prob <- state_next$log_prob + + feasible <- which(is.finite(state_next$log_prob)) + state_next$log_prob[feasible] <- state_next$log_prob[feasible] + sampling_path_route_log_prior( - candidate_idx = state_next$idx, + candidate_idx = state_next$idx[feasible], path_idx = path_idx, replace_i = block_start:block_end, kt = kt, @@ -962,6 +1045,88 @@ sampling_path_update_block <- function( sampled_idx } +#' Update Residence-Segment Boundaries +#' +#' @description +#' Internal Gibbs sweep used by `sampling_path_run_chain()`. Each boundary +#' between two consecutive residence segments is resampled over every position +#' that keeps both segments non-empty and fixed periods on their current cell. +#' Both cells, and hence the number of segments, are unchanged. Only the light +#' likelihood of the reassigned periods and the residence-dependent departure +#' terms of the two segments depend on the position: the movement kernel +#' between the two cells is shared, and fixed long periods keep the move inside +#' one route interval, so the route prior cancels. +#' +#' @param fixed logical vector flagging periods that keep their cell. +#' @return Integer vector of path indices. +#' @noRd +sampling_path_update_boundary <- function( + path_idx, + fixed, + lk, + kt, + component_weights +) { + nstap <- length(path_idx) + # log_stay[L]: staying through a segment of L periods; log_move[L]: leaving it. + log_stay <- c(0, cumsum(log1p(-kt$move_probability))) + log_move <- log(kt$move_probability) + n_boundary <- sum(path_idx[-1L] != path_idx[-nstap]) + + for (boundary_i in seq_len(n_boundary)) { + t <- which(path_idx[-1L] != path_idx[-nstap])[boundary_i] + s <- t + while (s > 1L && path_idx[s - 1L] == path_idx[t]) { + s <- s - 1L + } + u <- t + 1L + while (u < nstap && path_idx[u + 1L] == path_idx[t + 1L]) { + u <- u + 1L + } + window <- s:u + cell <- path_idx[c(t, t + 1L)] + + # Light of every window period at either cell, cumulated from each side. + # Supports are sorted, so findInterval() avoids hashing large supports. + light <- vapply( + window, + function(i) { + pos <- pmax( + findInterval(cell, lk[[i]]$idx, checkSorted = FALSE, checkNA = FALSE), + 1L + ) + ifelse(lk[[i]]$idx[pos] == cell, lk[[i]]$log_prob[pos], -Inf) + }, + numeric(2) + ) + end_1 <- s:(u - 1L) + len_1 <- end_1 - s + 1L + len_2 <- u - end_1 + light_prob <- cumsum(light[1L, ])[len_1] + + rev(cumsum(rev(light[2L, ])))[len_1 + 1L] + move_prob <- log_stay[len_1] + + log_move[len_1] + + log_stay[len_2] + + if (u < nstap) log_move[len_2] else 0 + + log_prob <- sampling_path_weight_log_prob( + light_prob, + component_weights[["light"]] + ) + + sampling_path_weight_log_prob(move_prob, component_weights[["movement"]]) + fixed_i <- window[fixed[window]] + log_prob[ + end_1 < max(s, fixed_i[path_idx[fixed_i] == cell[1L]]) | + end_1 >= min(u, fixed_i[path_idx[fixed_i] == cell[2L]]) + ] <- -Inf + + new_end <- end_1[sampling_path_sample_prob(exp(log_prob - max(log_prob)))] + path_idx[window] <- rep(cell, c(new_end - s + 1L, u - new_end)) + } + + path_idx +} + #' Update One Site #' #' @description @@ -1006,9 +1171,11 @@ sampling_path_update_site <- function( ), component_weights[["movement"]] ) - state$log_prob <- state$log_prob + + # The route prior cannot revive candidates excluded by hard support. + feasible <- which(is.finite(state$log_prob)) + state$log_prob[feasible] <- state$log_prob[feasible] + sampling_path_route_log_prior( - candidate_idx = state$idx, + candidate_idx = state$idx[feasible], path_idx = path_idx, replace_i = stap_i, kt = kt, @@ -1559,7 +1726,8 @@ sampling_path_route_log_prior <- function( ) for (interval_i in affected_interval) { - interval <- route_prior$intervals[interval_i, ] + interval_start_i <- route_prior$intervals$start_i[interval_i] + interval_end_i <- route_prior$intervals$end_i[interval_i] edge_i <- route_prior$interval_edges[[interval_i]] current_edge_distance <- if (is.null(route_edge_distance)) { sampling_path_index_distance_km( @@ -1609,15 +1777,15 @@ sampling_path_route_log_prior <- function( n_move_days <- n_move_days + (candidate_edge_distance > 0) } - start_idx <- if (interval$start_i %in% replace_i) { + start_idx <- if (interval_start_i %in% replace_i) { candidate_idx } else { - path_idx[interval$start_i] + path_idx[interval_start_i] } - end_idx <- if (interval$end_i %in% replace_i) { + end_idx <- if (interval_end_i %in% replace_i) { candidate_idx } else { - path_idx[interval$end_i] + path_idx[interval_end_i] } direct_distance <- sampling_path_index_distance_km( start_idx, @@ -1629,9 +1797,14 @@ sampling_path_route_log_prior <- function( is.finite(direct_distance) & direct_distance >= route_prior$min_direct_distance_km & n_move_days > 1L - route_ratio <- pmax(route_distance / direct_distance, 1) + # Floor the excess so exactly straight grid routes (e.g. along one column) + # keep the continuous limit of the density instead of a zero-density cliff. + route_excess <- pmax( + route_distance / direct_distance - 1, + .Machine$double.eps + ) support_covariates <- sampling_path_route_support_projection( - interval$duration_days, + route_prior$intervals$duration_days[interval_i], direct_distance, route_prior$support %||% NULL ) @@ -1641,16 +1814,18 @@ sampling_path_route_log_prior <- function( (support_covariates$duration_log - route_prior$duration_center) + route_prior$distance_coefficient * (support_covariates$distance_log - route_prior$distance_center) - log_density <- stats::dgamma( - route_ratio - 1, - shape = route_prior$gamma_shape, - scale = pmax( - (route_prior$detour %||% 1) * exp(linear_predictor), - .Machine$double.eps - ), - log = TRUE + # Closed-form Gamma log density: much cheaper than stats::dgamma(). + scale <- pmax( + (route_prior$detour %||% 1) * exp(linear_predictor), + .Machine$double.eps ) - log_density <- rep_len(log_density, length(route_ratio)) + shape <- route_prior$gamma_shape + log_density <- (shape - 1) * + log(route_excess) - + route_excess / scale - + lgamma(shape) - + shape * log(scale) + log_density <- rep_len(log_density, length(route_excess)) log_density[!is.finite(log_density)] <- log(.Machine$double.xmin) out[apply_prior] <- out[apply_prior] + log_density[apply_prior] } @@ -1701,17 +1876,12 @@ sampling_path_route_support_projection <- function( #' @return Numeric distance in kilometres. #' @noRd sampling_path_index_distance_km <- function(from_idx, to_idx, kt) { - n <- max(length(from_idx), length(to_idx)) - from_idx <- rep_len(from_idx, n) - to_idx <- rep_len(to_idx, n) + # A scalar endpoint recycles in the arithmetic instead of being replicated. from_row <- (from_idx - 1L) %% kt$nrow + 1L to_row <- (to_idx - 1L) %% kt$nrow + 1L - from_lat <- kt$lat_rad[from_row] - to_lat <- kt$lat_rad[to_row] - from_lon <- kt$lon_rad[(from_idx - 1L) %/% kt$nrow + 1L] - to_lon <- kt$lon_rad[(to_idx - 1L) %/% kt$nrow + 1L] - dlat <- to_lat - from_lat - dlon <- to_lon - from_lon + dlat <- kt$lat_rad[to_row] - kt$lat_rad[from_row] + dlon <- kt$lon_rad[(to_idx - 1L) %/% kt$nrow + 1L] - + kt$lon_rad[(from_idx - 1L) %/% kt$nrow + 1L] a <- sin(dlat / 2)^2 + kt$cos_lat[from_row] * kt$cos_lat[to_row] * sin(dlon / 2)^2 2 * 6378.137 * atan2(sqrt(pmin(a, 1)), sqrt(pmax(1 - a, 0))) @@ -1833,27 +2003,26 @@ sampling_path_transition_log_prob <- function( kt ) { kernel <- kt$movement_kernel - n <- max(length(source_idx), length(destination_idx)) - source_idx <- rep_len(source_idx, n) - destination_idx <- rep_len(destination_idx, n) + # A scalar endpoint recycles in the arithmetic instead of being replicated. source_row <- kt$cell_row[source_idx] drow <- kt$cell_row[destination_idx] - source_row dcol <- kt$cell_col[destination_idx] - kt$cell_col[source_idx] - out <- rep(-Inf, n) move_pos <- source_idx != destination_idx & abs(drow) <= kernel$max_drow & abs(dcol) <= kernel$max_dcol & kt$struct_mask[destination_idx] + out <- rep(-Inf, length(move_pos)) move_i <- which(move_pos) - lookup_idx <- source_row[move_i] + - (drow[move_i] + kernel$max_drow) * kt$nrow + - (dcol[move_i] + kernel$max_dcol) * kt$nrow * kernel$lookup_ndrow + lookup_idx <- (source_row + + (drow + kernel$max_drow) * kt$nrow + + (dcol + kernel$max_dcol) * kt$nrow * kernel$lookup_ndrow)[move_i] move_log_prob <- kernel$move_log_lookup[lookup_idx] finite_move <- is.finite(move_log_prob) move_i <- move_i[finite_move] + move_log_norm <- kernel$move_log_norm[source_idx] out[move_i] <- move_log_prob[finite_move] - - kernel$move_log_norm[source_idx[move_i]] + if (length(move_log_norm) == 1L) move_log_norm else move_log_norm[move_i] out } diff --git a/R/sampling_path_diagnostic.R b/R/sampling_path_diagnostic.R index 3253cc9..20be2f5 100644 --- a/R/sampling_path_diagnostic.R +++ b/R/sampling_path_diagnostic.R @@ -131,16 +131,21 @@ sampling_path_diagnostic <- function( out$iteration <- do.call( rbind, lapply(seq_along(checkpoint_n), function(i) { - checkpoint <- sampling_path_diagnostic_compute( - paths[paths$j <= saved_j[checkpoint_n[i]], , drop = FALSE], - stap, - spatial = FALSE, - sample_summary = out$sample_summary[ - out$sample_summary$j <= saved_j[checkpoint_n[i]], - , - drop = FALSE - ] - ) + # The full checkpoint equals the main computation. + checkpoint <- if (checkpoint_n[i] == length(saved_j)) { + out + } else { + sampling_path_diagnostic_compute( + paths[paths$j <= saved_j[checkpoint_n[i]], , drop = FALSE], + stap, + spatial = FALSE, + sample_summary = out$sample_summary[ + out$sample_summary$j <= saved_j[checkpoint_n[i]], + , + drop = FALSE + ] + ) + } sampling_path_diagnostic_checkpoint( checkpoint, checkpoint_label[i], diff --git a/R/sampling_path_prepare.R b/R/sampling_path_prepare.R index 54d1dcf..9dcea56 100644 --- a/R/sampling_path_prepare.R +++ b/R/sampling_path_prepare.R @@ -45,12 +45,39 @@ sampling_path_prepare <- function( map$stap$stap_id[map$stap$include] ) + # Pin known periods to the cell of their known location. A coastal site can + # fall on a water cell, which GeoPressureR replaces with a uniform land map; + # keep that cell in the structural domain so the path stays anchored. + included_stap <- which(map$stap$include) + known_i <- which( + !is.na(map$stap$known_lat[included_stap]) & + !is.na(map$stap$known_lon[included_stap]) + ) + known_idx <- vapply( + included_stap[known_i], + function(i) { + which.min(abs(g$lat - map$stap$known_lat[i])) + + (which.min(abs(g$lon - map$stap$known_lon[i])) - 1L) * + as.integer(g$dim[1]) + }, + integer(1) + ) + struct_mask[known_idx] <- TRUE + lk[known_i] <- Map( + function(state, idx) { + state$idx <- idx + state$log_prob <- 0 + state + }, + lk[known_i], + known_idx + ) + # A likelihood concentrated on the outer grid edge can indicate that the map # extent truncated the light surface before sampler support was constructed. map_boundary <- matrix(FALSE, nrow = g$dim[1], ncol = g$dim[2]) map_boundary[c(1L, g$dim[1]), ] <- TRUE map_boundary[, c(1L, g$dim[2])] <- TRUE - included_stap <- which(map$stap$include) unknown_stap <- is.na(map$stap$known_lat[included_stap]) | is.na(map$stap$known_lon[included_stap]) boundary_mass <- vapply( diff --git a/README.Rmd b/README.Rmd index 33516e4..e33a7c5 100644 --- a/README.Rmd +++ b/README.Rmd @@ -20,11 +20,7 @@ knitr::opts_chunk$set(collapse = TRUE, comment = "#>") GeoPathSampleR reconstructs bird migration trajectories from stationary-period -likelihood maps built with GeoPressureR. Its Gibbs sampler combines -residence-dependent departure, a ground-speed distribution conditional on -movement, and a prior on cumulative route distance between long periods. -Unknown long-period locations can be sampled from light evidence alone. The -package also provides posterior stay summaries and sampling diagnostics. +likelihood maps built with GeoPressureR. Its Gibbs sampler combines residence-dependent departure, a ground-speed distribution conditional on movement, and a prior on cumulative route distance between long periods. Unknown long-period locations can be sampled from light evidence alone. Site and block updates sample locations, while boundary sweeps resample movement timing with segment locations held fixed. The package also provides posterior stay summaries and sampling diagnostics. ## 📦 Installation diff --git a/README.md b/README.md index a9ca657..631b160 100644 --- a/README.md +++ b/README.md @@ -18,8 +18,10 @@ stationary-period likelihood maps built with GeoPressureR. Its Gibbs sampler combines residence-dependent departure, a ground-speed distribution conditional on movement, and a prior on cumulative route distance between long periods. Unknown long-period locations can be -sampled from light evidence alone. The package also provides posterior -stay summaries and sampling diagnostics. +sampled from light evidence alone. Site and block updates sample +locations, while boundary sweeps resample movement timing with segment +locations held fixed. The package also provides posterior stay summaries +and sampling diagnostics. ## 📦 Installation diff --git a/codemeta.json b/codemeta.json index f0ce151..5944535 100644 --- a/codemeta.json +++ b/codemeta.json @@ -8,7 +8,7 @@ "codeRepository": "https://github.com/GeoPressure/GeoPathSampleR", "issueTracker": "https://github.com/GeoPressure/GeoPathSampleR/issues", "license": "https://spdx.org/licenses/GPL-3.0", - "version": "0.1.0.9000", + "version": "0.2.0", "programmingLanguage": { "@type": "ComputerLanguage", "name": "R", @@ -315,7 +315,7 @@ }, "SystemRequirements": null }, - "fileSize": "52570.547KB", + "fileSize": "52621.971KB", "citation": [ { "@type": "SoftwareSourceCode", diff --git a/man/path_collapse.Rd b/man/path_collapse.Rd index 6bf9bf6..c13cc48 100644 --- a/man/path_collapse.Rd +++ b/man/path_collapse.Rd @@ -14,7 +14,7 @@ path_collapse(paths, stap = NULL, label_stap_id = FALSE) and \code{end}. If \code{NULL}, these columns are taken from \code{paths}.} \item{label_stap_id}{logical; if \code{TRUE}, replace \code{stap_id} in the output -with a character label of the form \verb{"\{stap_id_start\}->\{stap_id_end\}}. +with a character label of the form \code{"{stap_id_start}->{stap_id_end}"}. The numeric \code{stap_id_start} and \code{stap_id_end} columns are always kept.} } \value{ diff --git a/man/sampling_path.Rd b/man/sampling_path.Rd index 6a5dd22..6b1acf2 100644 --- a/man/sampling_path.Rd +++ b/man/sampling_path.Rd @@ -16,6 +16,7 @@ sampling_path( warmup = floor(iter/4), thin = 1, block_interval = 2, + boundary_interval = 1, thr_likelihood = 0.99, thr_gs = 2000/24, refresh = 10, @@ -27,7 +28,8 @@ sampling_path( \arguments{ \item{tag}{GeoPressureR tag object containing likelihood maps and \code{tag_set_map} parameters. Long periods are identified by \code{tag$stap$stap0}, -as created by \code{\link[GeoPressureR:tag_stap_daily]{GeoPressureR::tag_stap_daily()}}.} +as created by \code{\link[GeoPressureR:tag_stap_daily]{GeoPressureR::tag_stap_daily()}}. All stationary periods must +be included; excluded periods cannot be bridged by daily transitions.} \item{iter}{number of Gibbs iterations per chain, including \code{warmup}.} @@ -67,6 +69,13 @@ current residence segment as a single block. All other iterations are site iterations that update each eligible period once. Use \code{0} or \code{Inf} to disable block iterations; values must otherwise be at least \code{2}.} +\item{boundary_interval}{interval between boundary sweeps, which resample +the date of every boundary between consecutive residence segments while +keeping both locations. A sweep runs after the site or block update of +every \code{boundary_interval}-th iteration. The default \code{1} sweeps every +iteration; use \code{0} or \code{Inf} to disable them. Otherwise use a positive +integer. See Details for the conditional update.} + \item{thr_likelihood}{retained light-likelihood mass in each stationary period. Lower values keep fewer cells and speed up the sampler, but exclude more of the light likelihood surface entirely. The default \code{0.99} is a @@ -115,16 +124,18 @@ with shape \eqn{\alpha_R = 1.22} and where \eqn{T_\ell} is the interval duration in days and \eqn{D_\ell} is in km. Both covariates are clamped to the range of the calibration data. The term is omitted when \eqn{D_\ell < 300} km or when the interval contains at -most one daily displacement. +most one non-zero daily displacement. Because \code{route_detour} multiplies the Gamma scale, the prior mean and every prior quantile of the route excess \eqn{R_\ell - 1} scale linearly with it. For example, an interval of about 31 days with endpoints about 3,500 km apart has \eqn{\theta_\ell \approx 0.15}, i.e. a prior mean route ratio of about 1.19; \code{route_detour = 2} raises it to about 1.37 and \code{0.5} lowers it -to about 1.09. \code{route_detour = 0} collapses the scale, so any detour is -strongly penalised. The prior remains soft: the realised detour of the -posterior also depends on the light and movement components. +to about 1.09. \code{route_detour = 0} uses machine epsilon as the scale, so any +detour is strongly penalised. Route excess is also floored at machine +epsilon when evaluating the log density. The prior remains soft: the +realised detour of the posterior also depends on the light and movement +components. \strong{Long periods and modular cut.} Twilights accurately locate long periods but leave daily-period locations ambiguous. With @@ -139,6 +150,27 @@ information flows from the long-period likelihoods to the daily path, but the movement model cannot pull long-period locations away from their light likelihood. Known-location periods are always fixed. +\strong{Site, block and boundary updates.} By default, odd iterations update +each eligible period separately and even iterations update each current +residence segment at one shared location (\code{block_interval = 2}). Site +sweeps alternate forwards and backwards. Block updates use the intersection +of the periods' retained light supports and cannot span a long period when +the route prior is active or include a long period drawn from light alone. + +After either update, \code{boundary_interval = 1} adds a forward sweep over the +boundaries between consecutive residence segments. Each boundary is drawn +over all positions within its two segments that keep both non-empty, +retain their two locations, respect each period's light support, and keep +fixed periods at their current location. This changes departure and arrival +timing without changing the number or sequence of residence segments. +Candidate positions are weighted by the light likelihood and +residence-dependent departure terms, with \code{component_weights} applied. +The transition between the two locations is unchanged, and fixed +long-period endpoints keep cumulative distance within each route interval +unchanged, so the movement kernel and route prior cancel in this update. +Boundary sweeps run during warmup as well as sampling, before draws are +saved. Their frequency changes the update schedule, not the target model. + \strong{Computational approximations.} Impossible geography, such as water after masking, is removed from the state space. The arguments \code{thr_likelihood}, \code{thr_gs} and \code{movement$approx_eps} control how aggressively the sampler diff --git a/vignettes/getting-started.qmd b/vignettes/getting-started.qmd index 0d368b9..81f1fd6 100644 --- a/vignettes/getting-started.qmd +++ b/vignettes/getting-started.qmd @@ -32,7 +32,7 @@ knitr::opts_knit$set( ### Set-up your project folder -We recommend using GeoPressureTemplate as the starting point for your analysis. It provides a folder structure where all files for the analysis have their standardized location. Check the instruction on the [dedicated GeoPressureManual chapter](https://geopressure.org/GeoPressureManual/geopressuretemplate-intro.html). +We recommend using GeoPressureTemplate as the starting point for your analysis. It provides a folder structure where all files for the analysis have their standardised location. Check the instructions in the [dedicated GeoPressureManual chapter](https://geopressure.org/GeoPressureManual/geopressuretemplate-intro.html). Here, we use the `14OI` example contained in the GeoPathSampleR package. The vignette working directory is set to the package's `extdata` folder, so GeoPressureR can use its standard relative paths. @@ -162,7 +162,7 @@ plot_twl_calib(tag) ``` ### Add the water mask -For landbird, a water mask removes cells from the set of locations that can be sampled. +For land birds, a water mask removes cells from the set of locations that can be sampled. ```{r} tag$map_light <- map_add_mask_water(tag$map_light) @@ -177,7 +177,7 @@ The movement model implemented in `sampling_path()` combines three sources of in - The light likelihood describes where each stationary period could be located. - The movement model describes plausible speeds and the probability of departing after a given residence time. -- The route-length prior provides information about the expected detour of long intervals between stationary periods. +- The prior on the route-distance ratio describes cumulative travel relative to direct endpoint distance between consecutive long periods. ### Defining parameters @@ -188,7 +188,7 @@ The movement model defines the daily transition between consecutive stationary p - A ground-speed kernel (`method`, `shape`, `scale`, `low_speed_fix`, `zero_speed_ratio`), evaluated with `GeoPressureR::speed2prob()`. Speed is the displacement between two consecutive periods divided by 24 h. The default is a Gamma distribution with shape `1.26` and scale `10.3` km/h. - A residence-dependent departure probability `move_stay(t)`, where `t` is the number of consecutive days already spent at the current location. The default declines exponentially from `p_0 = 0.607` just after arrival to a plateau `p_inf = 0.137` with time constant `tau = 0.932` days. `move_stay_parameters` only stores these values for `plot_movement()`; the sampler uses `move_stay`. -The default `movement` argument of `sampling_path()` is `sampling_path_default_movement()`, calibrated on multi-sensor tracks of land birds. We write it out below to show its structure; `movement <- sampling_path_default_movement()` gives the same list. See `?sampling_path_default_movement` for details. +The default `movement` argument of `sampling_path()` is `sampling_path_default_movement()`, calibrated on multi-sensor reference trajectories of land birds. We write it out below to show its structure; `movement <- sampling_path_default_movement()` gives the same list. See `?sampling_path_default_movement` for details. ```{r} movement <- list( @@ -212,8 +212,8 @@ plot_movement(movement) The other arguments of `sampling_path()` that change the scientific model are: - `component_weights = c(light = 1, movement = 1, route = 1)` uses each component without rescaling. Each weight scales that component's log probability, so changing even all three weights by the same factor changes the sampled distribution. Setting a weight to zero removes that component's soft contribution while retaining hard support constraints; use this only for an explicit sensitivity analysis. -- `route_detour = 1` keeps the calibrated route prior. The value multiplies the expected excess route distance (route length over direct distance, minus one) between consecutive long periods. For instance, for a prior mean route is 19% longer than direct; `2` doubles this excess, `0.5` halves it and `0` favours direct routes. -- `long_period_light_only = TRUE` applies a modular cut: unknown long-period locations are drawn from their light likelihood alone, and movement only informs the daily periods between them. `FALSE` lets movement also update long periods. +- `route_detour = 1` keeps the calibrated route prior. The value multiplies the expected excess route distance (route length over direct distance, minus one) between consecutive long periods. For an interval of about 31 days with endpoints about 3,500 km apart, the prior mean route is about 19% longer than direct; `2` doubles this excess, `0.5` halves it and `0` strongly penalises detours. The prior applies only when endpoints are at least 300 km apart and the trajectory contains at least two non-zero daily displacements. +- `long_period_light_only = TRUE` applies a modular cut: unknown long-period locations are drawn from their light likelihood alone, and movement only informs the daily periods between them. `FALSE` lets the movement and route components also update unknown long periods. Known locations remain fixed. #### Computational controls @@ -225,15 +225,18 @@ The following arguments control computation and diagnostics rather than the scie | `warmup` | `floor(iter / 4)` | Initial iterations discarded while the chain reaches its typical posterior region. | | `thin` | `1` | Saving interval. The number saved is `floor((iter - warmup) / thin)`. Keep `1` unless storage is limiting; thinning does not improve mixing. | | `chains` | `1` | Number of independent chains. Use at least four for final inference. | -| `block_interval` | `2` | Frequency of block updates, which move a whole stay at once. Keep the optimized default unless explicitly comparing sampler efficiency. | +| `block_interval` | `2` | Every second iteration updates whole residence segments at shared locations; the others update individual periods. Use `0` or `Inf` to disable block iterations, or an integer of at least `2`. | +| `boundary_interval` | `1` | Sweep over residence-segment boundaries after every site or block iteration to resample movement timing while retaining segment locations. Use `0` or `Inf` to disable sweeps, or a positive integer to run them less often. | | `thr_likelihood` | `0.99` | Retained light-likelihood mass used to reduce the state space. The conservative default is appropriate for final inference. | | `thr_gs` | `2000 / 24` | Hard truncation of the speed kernel, in km/h (2000 km per day). Speeds are already weighted by the movement model, so set it above the speeds the kernel supports rather than at typical speeds; the default excludes less than 0.1% of the default kernel. Larger values cost computation; smaller values forbid transitions outright. | | `workers` | `1` | Number of parallel workers for independent chains. Use `1` while debugging and increase only when appropriate for the system. | | `seed` | `NULL` | Random-number seed for reproducible chains. | +All stationary periods must remain included in `tag_set_map()`: the sampler treats each transition as one daily step and cannot bridge excluded periods. + ### First run -The sampler uses Gibbs updates: it repeatedly updates one stationary-period location, or a block of locations, conditional on the current values of the other periods. The resulting draws are correlated, so the run must be long enough to explore the posterior and must be checked with multiple independent chains. +The sampler uses Gibbs updates: it alternates individual-period and residence-segment location updates, then resamples the boundaries between consecutive segments. Each boundary update keeps both segments non-empty and their locations unchanged, and draws a new boundary from the light likelihood and residence-dependent departure probabilities. Known locations and long periods drawn from light alone remain fixed. The route-distance prior cancels because boundary shifts leave cumulative distance within each long-period interval unchanged. See `?sampling_path` for the full update schedule and constraints. The resulting draws are correlated, so the run must be long enough to explore the posterior and must be checked with multiple independent chains. Begin with a short run to verify that the likelihood, movement support, and route prior are compatible. This run is for checking the workflow, not for inference.