From 6da2d1588fb34d83d3aa29265d56bb2bccdfa92f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Sun, 25 Oct 2026 09:46:24 +0100 Subject: [PATCH 1/9] Speed up route prior, path summaries and diagnostics Evaluate the route prior only on movement-feasible candidates with a closed-form Gamma log density, split rows once in path_summary(), and reuse the full diagnostic for the final checkpoint. Outputs are unchanged. Co-Authored-By: Claude Opus 5.5 --- R/path_summary.R | 45 +++++++++++++++++++++++------------- R/sampling_path.R | 42 ++++++++++++++++++--------------- R/sampling_path_diagnostic.R | 25 ++++++++++++-------- 3 files changed, 68 insertions(+), 44 deletions(-) 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..16e8853 100644 --- a/R/sampling_path.R +++ b/R/sampling_path.R @@ -932,9 +932,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, @@ -1006,9 +1007,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 +1562,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 +1613,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, @@ -1631,7 +1635,7 @@ sampling_path_route_log_prior <- function( n_move_days > 1L route_ratio <- pmax(route_distance / direct_distance, 1) support_covariates <- sampling_path_route_support_projection( - interval$duration_days, + route_prior$intervals$duration_days[interval_i], direct_distance, route_prior$support %||% NULL ) @@ -1641,15 +1645,17 @@ 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 ) + shape <- route_prior$gamma_shape + log_density <- (shape - 1) * + log(route_ratio - 1) - + (route_ratio - 1) / scale - + lgamma(shape) - + shape * log(scale) log_density <- rep_len(log_density, length(route_ratio)) log_density[!is.finite(log_density)] <- log(.Machine$double.xmin) out[apply_prior] <- out[apply_prior] + log_density[apply_prior] 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], From 49e7735ed548979987b060af50c0fab223951a95 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Sun, 25 Oct 2026 09:49:33 +0100 Subject: [PATCH 2/9] Fix known water sites, straight-route prior and excluded periods - Pin known periods to the cell of their known location and keep it in the structural domain, so a coastal site falling on a water cell no longer becomes a uniform land map. - Floor the route excess at machine epsilon so routes along one grid column are not assigned a near-zero prior density. - Abort when stationary periods are excluded, since each transition is modelled as one daily step. Co-Authored-By: Claude Opus 5.5 --- R/sampling_path.R | 21 +++++++++++++++++---- R/sampling_path_prepare.R | 29 ++++++++++++++++++++++++++++- 2 files changed, 45 insertions(+), 5 deletions(-) diff --git a/R/sampling_path.R b/R/sampling_path.R index 16e8853..492c3b0 100644 --- a/R/sampling_path.R +++ b/R/sampling_path.R @@ -143,6 +143,14 @@ sampling_path_with_route_model <- function( } 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, @@ -1633,7 +1641,12 @@ 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( route_prior$intervals$duration_days[interval_i], direct_distance, @@ -1652,11 +1665,11 @@ sampling_path_route_log_prior <- function( ) shape <- route_prior$gamma_shape log_density <- (shape - 1) * - log(route_ratio - 1) - - (route_ratio - 1) / scale - + log(route_excess) - + route_excess / scale - lgamma(shape) - shape * log(scale) - log_density <- rep_len(log_density, length(route_ratio)) + 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] } 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( From 88fa26788f379f57eeca7fc161de726dc4f00bba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Sun, 25 Oct 2026 12:13:31 +0100 Subject: [PATCH 3/9] Track sampling progress without <<- Read the position back from the progress bar so jarl's undesirable operator rule passes. Co-Authored-By: Claude Opus 5.5 --- R/sampling_path.R | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/R/sampling_path.R b/R/sampling_path.R index 492c3b0..9cff627 100644 --- a/R/sampling_path.R +++ b/R/sampling_path.R @@ -203,7 +203,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) { @@ -229,8 +228,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 + ) } } From 0b989d112e77cf6d8ec1bb4e4cbc34bea3fa9751 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Sun, 25 Oct 2026 13:59:02 +0100 Subject: [PATCH 4/9] Avoid replicating scalar endpoints in movement and distance lookups A single neighbouring cell now recycles in the arithmetic instead of being copied to the candidate length and re-gathered. Results are identical. Co-Authored-By: Claude Opus 5.5 --- R/sampling_path.R | 28 +++++++++++----------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/R/sampling_path.R b/R/sampling_path.R index 9cff627..65ac70a 100644 --- a/R/sampling_path.R +++ b/R/sampling_path.R @@ -1721,17 +1721,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))) @@ -1853,27 +1848,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 } From f15e91eae33f09211222533f9b97c3ec3cd7fd24 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Sun, 25 Oct 2026 14:12:01 +0100 Subject: [PATCH 5/9] Add boundary sweeps to the path sampler A boundary sweep resamples, by Gibbs, the date of every boundary between consecutive residence segments while keeping both locations. It complements block updates, which move a segment without changing its dates. Controlled by boundary_interval (0 keeps the previous sampler). Co-Authored-By: Claude Opus 5.5 --- R/sampling_path.R | 123 +++++++++++++++++++++++++++++++++++++++++++ man/sampling_path.Rd | 6 +++ 2 files changed, 129 insertions(+) diff --git a/R/sampling_path.R b/R/sampling_path.R index 65ac70a..3a170e6 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 = 0, 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,6 +147,11 @@ 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. @@ -249,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, @@ -268,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, @@ -378,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, @@ -461,6 +475,10 @@ 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. Use `0` or `Inf` to disable them. #' @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 @@ -587,6 +605,7 @@ sampling_path <- function( warmup = floor(iter / 4), thin = 1, block_interval = 2, + boundary_interval = 0, thr_likelihood = 0.99, thr_gs = 2000 / 24, refresh = 10, @@ -606,6 +625,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, @@ -654,6 +674,7 @@ sampling_path_run_chain <- function( warmup, thin, block_interval, + boundary_interval, long_period_i, component_weights, route_prior, @@ -684,6 +705,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( @@ -828,6 +851,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 @@ -972,6 +1019,82 @@ 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. + light <- vapply( + window, + function(i) lk[[i]]$log_prob[match(cell, lk[[i]]$idx)], + numeric(2) + ) + light[is.na(light)] <- -Inf + 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 diff --git a/man/sampling_path.Rd b/man/sampling_path.Rd index 6a5dd22..716fa07 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 = 0, thr_likelihood = 0.99, thr_gs = 2000/24, refresh = 10, @@ -67,6 +68,11 @@ 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. Use \code{0} or \code{Inf} to disable them.} + \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 From a0799ab64e495405acdcb91c237fde6cce99a8e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Sun, 25 Oct 2026 18:04:41 +0100 Subject: [PATCH 6/9] Look up boundary-sweep light values by binary search Supports are sorted, so findInterval() replaces match(), which hashed the whole support for every period. Results are identical. Co-Authored-By: Claude Opus 5.5 --- R/sampling_path.R | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/R/sampling_path.R b/R/sampling_path.R index 3a170e6..e0b0742 100644 --- a/R/sampling_path.R +++ b/R/sampling_path.R @@ -1061,12 +1061,18 @@ sampling_path_update_boundary <- function( 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) lk[[i]]$log_prob[match(cell, lk[[i]]$idx)], + 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) ) - light[is.na(light)] <- -Inf end_1 <- s:(u - 1L) len_1 <- end_1 - s + 1L len_2 <- u - end_1 From 9f8cf753042a082211ce38af1667db95e2914f09 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Sun, 25 Oct 2026 20:17:42 +0100 Subject: [PATCH 7/9] Run boundary sweeps every iteration by default Across the 21 validation tags, boundary sweeps let every tag meet the convergence criterion (19 without) at about 4% of a site sweep's cost. Co-Authored-By: Claude Opus 5.5 --- R/sampling_path.R | 7 ++++--- man/sampling_path.Rd | 5 +++-- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/R/sampling_path.R b/R/sampling_path.R index e0b0742..f1c7425 100644 --- a/R/sampling_path.R +++ b/R/sampling_path.R @@ -75,7 +75,7 @@ sampling_path_with_route_model <- function( warmup = floor(iter / 4), thin = 1, block_interval = 2, - boundary_interval = 0, + boundary_interval = 1, thr_likelihood = 0.99, thr_gs = 2000 / 24, refresh = 10, @@ -478,7 +478,8 @@ sampling_path_with_route_model <- function( #' @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. Use `0` or `Inf` to disable them. +#' every `boundary_interval`-th iteration. The default `1` sweeps every +#' iteration at a small cost; use `0` or `Inf` to disable them. #' @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 @@ -605,7 +606,7 @@ sampling_path <- function( warmup = floor(iter / 4), thin = 1, block_interval = 2, - boundary_interval = 0, + boundary_interval = 1, thr_likelihood = 0.99, thr_gs = 2000 / 24, refresh = 10, diff --git a/man/sampling_path.Rd b/man/sampling_path.Rd index 716fa07..c5fcdd2 100644 --- a/man/sampling_path.Rd +++ b/man/sampling_path.Rd @@ -16,7 +16,7 @@ sampling_path( warmup = floor(iter/4), thin = 1, block_interval = 2, - boundary_interval = 0, + boundary_interval = 1, thr_likelihood = 0.99, thr_gs = 2000/24, refresh = 10, @@ -71,7 +71,8 @@ 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. Use \code{0} or \code{Inf} to disable them.} +every \code{boundary_interval}-th iteration. The default \code{1} sweeps every +iteration at a small cost; use \code{0} or \code{Inf} to disable them.} \item{thr_likelihood}{retained light-likelihood mass in each stationary period. Lower values keep fewer cells and speed up the sampler, but exclude From 5aedeb71671f9f3863636636544065fd99fe3e2d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Mon, 26 Oct 2026 06:57:50 +0100 Subject: [PATCH 8/9] Document boundary sweeps and prepare v0.2.0 --- CITATION.cff | 4 ++-- DESCRIPTION | 2 +- NEWS.md | 6 ++++++ R/path_collapse.R | 2 +- R/sampling_path.R | 37 +++++++++++++++++++++++++++++------ README.Rmd | 6 +----- README.md | 6 ++++-- codemeta.json | 4 ++-- man/path_collapse.Rd | 2 +- man/sampling_path.Rd | 37 +++++++++++++++++++++++++++++------ vignettes/getting-started.qmd | 19 ++++++++++-------- 11 files changed, 91 insertions(+), 34 deletions(-) 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..b98a979 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,12 +2,18 @@ ## 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 +- [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/sampling_path.R b/R/sampling_path.R index f1c7425..6df5e1b 100644 --- a/R/sampling_path.R +++ b/R/sampling_path.R @@ -446,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 @@ -479,7 +480,8 @@ sampling_path_with_route_model <- function( #' 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 at a small cost; use `0` or `Inf` to disable them. +#' 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 @@ -514,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 @@ -538,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 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 c5fcdd2..6b1acf2 100644 --- a/man/sampling_path.Rd +++ b/man/sampling_path.Rd @@ -28,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}.} @@ -72,7 +73,8 @@ disable block iterations; values must otherwise be at least \code{2}.} 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 at a small cost; use \code{0} or \code{Inf} to disable them.} +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 @@ -122,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 @@ -146,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. From f0053b814ce0b7d470a94c7ab4d5ce36129dabed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Rapha=C3=ABl=20Nussbaumer?= Date: Mon, 26 Oct 2026 06:59:48 +0100 Subject: [PATCH 9/9] Exclude redundant GPL text from package builds --- .Rbuildignore | 1 + NEWS.md | 2 ++ 2 files changed, 3 insertions(+) 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/NEWS.md b/NEWS.md index b98a979..d142af5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,6 +10,8 @@ ## 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.