From 7a02eea0d167a02b71b6d013cc228f6203a68162 Mon Sep 17 00:00:00 2001 From: soukhova Date: Thu, 23 Oct 2025 17:57:50 -0400 Subject: [PATCH 1/4] feature suggestion: total_constrained() accessibility measure from forthcoming Soukhov et al. 2025. NOTE: attempted to copy the formatting pattern from spatial_availability.R. --- R/total_constrained.R | 156 ++++++++++++++ inst/REFERENCES.bib | 14 ++ man/total_constrained.Rd | 123 +++++++++++ tests/testthat/test-total_constrained.R | 276 ++++++++++++++++++++++++ 4 files changed, 569 insertions(+) create mode 100644 R/total_constrained.R create mode 100644 man/total_constrained.Rd create mode 100644 tests/testthat/test-total_constrained.R diff --git a/R/total_constrained.R b/R/total_constrained.R new file mode 100644 index 0000000..d3f8c95 --- /dev/null +++ b/R/total_constrained.R @@ -0,0 +1,156 @@ +#' Total constrained +#' +#' Calculates total constrained accessibility, part of a family of accessibility +#'measures proposed in +#'\insertCite{soukhovFamilyAccessibilityMeasures2025;textual}{accessibility} that +#' proportionally allocates the total sum of opportunities in the region based +#' on the relative travel impedance between zones. It is linearly proportion to +#' the conventional gravity-based Hansen-type accessibility measure and units of +#' the resulting values can be understood as the number of opportunities in the +#' region that can be potentially interacted with. +#' @template description_generic_cost +#' +#' @template travel_matrix +#' @template land_use_data +#' @template opportunity +#' @template travel_cost +#' @template decay_function +#' @template group_by +#' @template fill_missing_ids_combinations +#' @param detailed_results A `logical`. Whether to return total constrained +#' results aggregated by origin-destination pair (`TRUE`) or by origin +#' (`FALSE`, the default). When `TRUE`, the output also includes the balancing +#' factor `kappa_tot` used in the calculation and the unconstrained +#' accessibility value `weighted_opportunity`. Please note that the argument +#' `fill_missing_ids` does not affect the output when `detailed_results` is +#' `TRUE`. +#' +#' @template return_accessibility +#' +#' @references +#' \insertAllCited{} +#' +#' @examplesIf identical(tolower(Sys.getenv("NOT_CRAN")), "true") +#' # the example below is based on Soukhov et al. (2023) paper +#' +#' travel_matrix <- data.table::data.table( +#' from_id = rep(c("A", "B", "C"), each = 3), +#' to_id = as.character(rep(1:3, 3)), +#' travel_time = c(15, 30, 100, 30, 15, 100, 100, 100, 15) +#' ) +#' land_use_data <- data.table::data.table( +#' id = c("A", "B", "C", "1", "2", "3"), +#' population = c(50000, 150000, 10000, 0, 0, 0), +#' jobs = c(0, 0, 0, 100000, 100000, 10000) +#' ) +#' +#' df <- total_constrained( +#' travel_matrix, +#' land_use_data, +#' opportunity = "jobs", +#' travel_cost = "travel_time", +#' decay_function = decay_exponential(decay_value = 0.1) +#' ) +#' df +#' +#' detailed_df <- total_constrained( +#' travel_matrix, +#' land_use_data, +#' opportunity = "jobs", +#' travel_cost = "travel_time", +#' decay_function = decay_exponential(decay_value = 0.1), +#' detailed_results = TRUE +#' ) +#' detailed_df +#' +#' @export +total_constrained <- function(travel_matrix, + land_use_data, + opportunity, + travel_cost, + decay_function, + group_by = character(0), + fill_missing_ids = TRUE, + detailed_results = FALSE) { + checkmate::assert_string(opportunity) + checkmate::assert_string(travel_cost) + checkmate::assert_logical(detailed_results, len = 1, any.missing = FALSE) + assert_decay_function(decay_function) + assert_group_by(group_by) + assert_detailed_fill_missing_ids(fill_missing_ids, detailed_results) + assert_travel_matrix(travel_matrix, travel_cost, group_by) + assert_land_use_data( + land_use_data, + travel_matrix, + opportunity) + + if (!inherits(travel_matrix, "data.table")) { + original_class <- class(travel_matrix) + data <- data.table::as.data.table(travel_matrix) + } else { + data <- data.table::copy(travel_matrix) + } + + if (!inherits(land_use_data, "data.table")) { + land_use_data <- data.table::as.data.table(land_use_data) + } + + merge_by_reference(data, land_use_data, opportunity, left_df_idcol = "to_id") + + data <- apply_gravity_measure(data, decay_function, travel_cost) + + groups <- c("from_id", group_by) + if ("decay_function_arg" %in% names(data)) { + groups <- c(groups, "decay_function_arg") + group_by <- c(group_by, "decay_function_arg") + } + + warn_extra_cols( + travel_matrix, + travel_cost, + group_id = "from_id", + groups = groups + ) + + # Calculate unconstrained accessibility (weighted_opportunities) per ij: W_j * f(c_ij) + .opportunity_colname <- opportunity + data[, weighted_opportunity := get(.opportunity_colname) * opp_weight] + + # Calculate kappa_tot for each ij: kappa_tot_ij = (W_j * f(c_ij)) / sum_over_all_ij (W_j * f(c_ij)) + total_weighted_system <- data[, sum(weighted_opportunity)] + data[, kappa_tot := weighted_opportunity / total_weighted_system] + + # Total opportunities in the region + total_opportunities_region <- sum(land_use_data[[.opportunity_colname]]) + + # Calculate total constrained opportunity: kappa_tot_ij * total regional opportunities + data[, constrained_opportunity := kappa_tot * total_opportunities_region] + + if (detailed_results) { + # Return detailed results with all intermediate calculations + data[, total_opportunities_region := total_opportunities_region] + + # Drop unnecessary columns and rename + cols_to_drop <- c(travel_cost, opportunity, "opp_weight", total_opportunities_region) + access <- data[, (cols_to_drop) := NULL] + + data.table::setnames(access, "constrained_opportunity", opportunity) + } else { + # Return aggregated results by origin + access <- data[ + , + .(access = sum(constrained_opportunity)), + by = c("from_id", group_by) + ] + + if (fill_missing_ids) { + access <- fill_missing_ids(access, travel_matrix, groups) + } + + data.table::setnames(access, c("from_id", "access"), c("id", opportunity)) + } + + if (exists("original_class")) class(access) <- original_class + + return(access[]) +} diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 183a2dd..2f32593 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -162,3 +162,17 @@ @article{tomasiello2023time keywords = {Cumulative accessibility,Equity,MTUP,Public transport,Time interval,Transport policy,Travel time}, file = {/home/dhersz/Zotero/storage/NNWIAXLP/S0143622823001388.html} } + + +@misc{soukhovFamilyAccessibilityMeasures2025, + title = {A family of accessibility measures derived from spatial interaction principles}, + url = {https://osf.io/a9dxb_v1}, + doi = {10.31219/osf.io/a9dxb_v1}, + abstract = {Transportation planning has long prioritized the efficiency of movement, or mobility. However, the concept of accessibility represents a more comprehensive evolution, shifting focus from mere movement to the potential to reach (i.e., spatially interact) with desired destinations. Despite growing recognition of accessibility-based planning approaches, the concept remains fragmented, with inconsistent definitions and unclear interpretations. This work's aim is to clarify and unify the concept of accessibility by connecting it into spatial interaction modeling. We demonstrate that widely used mobility and accessibility models, such as gravity-based accessibility and spatial interaction models, share common theoretical roots. From this foundation, this paper offers three contributions: (A) we introduce a family of accessibility measures within the principles of spatial interaction, and (B) formally define four members of the family, namely the 'unconstrained' measure (i.e., Hansen-type accessibility), the 'total constrained' measure (i.e., a constrained version of the Hansen-type accessibility), the 'singly constrained' measure (i.e., related to the popular two step floating catchment approach - 2SFCA), and the 'doubly constrained' measure representing realized interactions or 'access', effectively equal to the doubly constrained spatial interaction model; and (C) we demonstrate the interpretability advantages of the family, as these constrained accessibility measures yield values in units of the number of potential "opportunities for spatial interaction" or "population for spatial interaction" for each zone and zonal flow. The family of accessibility measures proposed here clarifies the concept of 'potential' in accessibility, demonstrates theoretical and formulaic linkages across popular accessibility and spatial interaction models, and reintroduces measurement units into accessibility measures. By doing so, we believe this family of measures can unlock a clearer, more interpretable, and cohesive foundation for accessibility analysis.}, + publisher = {{OSF}}, + author = {Soukhov, Anastasia and Pereira, Rafael and Higgins, Christopher and Paez, Antonio}, + urldate = {2025-09-19}, + date = {2025-05-23}, + langid = {english}, + keywords = {accessibility, balancing factors, competitive accessibility, constraints, doubly constrained, singly constrained, spatial interaction modelling, total constrained}, +} diff --git a/man/total_constrained.Rd b/man/total_constrained.Rd new file mode 100644 index 0000000..a033e6f --- /dev/null +++ b/man/total_constrained.Rd @@ -0,0 +1,123 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/total_constrained.R +\name{total_constrained} +\alias{total_constrained} +\title{Total constrained} +\usage{ +total_constrained( + travel_matrix, + land_use_data, + opportunity, + travel_cost, + decay_function, + group_by = character(0), + fill_missing_ids = TRUE, + detailed_results = FALSE +) +} +\arguments{ +\item{travel_matrix}{A data frame. The travel matrix describing the costs +(i.e. travel time, distance, monetary cost, etc.) between the origins and +destinations in the study area. Must contain the columns \code{from_id}, \code{to_id} +and any others specified in \code{travel_cost}.} + +\item{land_use_data}{A data frame. The distribution of opportunities within +the study area cells. Must contain the columns \code{id} and any others +specified in \code{opportunity}.} + +\item{opportunity}{A string. The name of the column in \code{land_use_data} +with the number of opportunities/resources/services to be considered when +calculating accessibility levels.} + +\item{travel_cost}{A string. The name of the column in \code{travel_matrix} +with the travel cost between origins and destinations.} + +\item{decay_function}{A \code{fuction} that converts travel cost into an +impedance factor used to weight opportunities. This function should take a +\code{numeric} vector and also return a \code{numeric} vector as output, with the +same length as the input. For convenience, the package currently includes +the following functions: \code{\link[=decay_binary]{decay_binary()}}, \code{\link[=decay_exponential]{decay_exponential()}}, +\code{\link[=decay_power]{decay_power()}} and \code{\link[=decay_stepped]{decay_stepped()}}. See the documentation of each decay +function for more details.} + +\item{group_by}{A \code{character} vector. When not \code{character(0)} (the default), +indicates the \code{travel_matrix} columns that should be used to group the +accessibility estimates by. For example, if \code{travel_matrix} includes a +departure time column, that specifies the departure time of each entry in +the data frame, passing \code{"departure_time"} to this parameter results in +accessibility estimates grouped by origin and by departure time.} + +\item{fill_missing_ids}{A \code{logical}. When calculating grouped accessibility +estimates (i.e. when \code{by_col} is not \code{NULL}), some combinations of groups +and origins may be missing. For example, if a single trip can depart from +origin \code{A} at 7:15am and reach destination \code{B} within 55 minutes, but no +trips departing from \code{A} at 7:30am can be completed at all, this second +combination will not be included in the output. When \code{TRUE} (the default), +the function identifies which combinations would be left out and fills +their respective accessibility values with 0, which incurs in a +performance penalty.} + +\item{detailed_results}{A \code{logical}. Whether to return total constrained +results aggregated by origin-destination pair (\code{TRUE}) or by origin +(\code{FALSE}, the default). When \code{TRUE}, the output also includes the balancing +factor \code{kappa_tot} used in the calculation and the unconstrained +accessibility value \code{weighted_opportunity}. Please note that the argument +\code{fill_missing_ids} does not affect the output when \code{detailed_results} is +\code{TRUE}.} +} +\value{ +A data frame containing the accessibility estimates for each +origin/destination (depending if \code{active} is \code{TRUE} or \code{FALSE}) in the +travel matrix. +} +\description{ +Calculates total constrained accessibility, part of a family of accessibility +measures proposed in +\insertCite{soukhovFamilyAccessibilityMeasures2025;textual}{accessibility} that +proportionally allocates the total sum of opportunities in the region based +on the relative travel impedance between zones. It is linearly proportion to +the conventional gravity-based Hansen-type accessibility measure and units of +the resulting values can be understood as the number of opportunities in the +region that can be potentially interacted with. + +This function is generic over any kind of numeric travel cost, +such as distance, time and money. +} +\examples{ +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} +# the example below is based on Soukhov et al. (2023) paper + +travel_matrix <- data.table::data.table( + from_id = rep(c("A", "B", "C"), each = 3), + to_id = as.character(rep(1:3, 3)), + travel_time = c(15, 30, 100, 30, 15, 100, 100, 100, 15) +) +land_use_data <- data.table::data.table( + id = c("A", "B", "C", "1", "2", "3"), + population = c(50000, 150000, 10000, 0, 0, 0), + jobs = c(0, 0, 0, 100000, 100000, 10000) +) + +df <- total_constrained( + travel_matrix, + land_use_data, + opportunity = "jobs", + travel_cost = "travel_time", + decay_function = decay_exponential(decay_value = 0.1) +) +df + +detailed_df <- total_constrained( + travel_matrix, + land_use_data, + opportunity = "jobs", + travel_cost = "travel_time", + decay_function = decay_exponential(decay_value = 0.1), + detailed_results = TRUE +) +detailed_df +\dontshow{\}) # examplesIf} +} +\references{ +\insertAllCited{} +} diff --git a/tests/testthat/test-total_constrained.R b/tests/testthat/test-total_constrained.R new file mode 100644 index 0000000..c1082a7 --- /dev/null +++ b/tests/testthat/test-total_constrained.R @@ -0,0 +1,276 @@ +# if running manually, please run the following line first: +# source("tests/testthat/setup.R") + +tester <- function( + travel_matrix = smaller_matrix, + land_use_data = get("land_use_data", envir = parent.frame()), + opportunity = "jobs", + travel_cost = "travel_time", + decay_function = decay_exponential(0.1), + group_by = "mode", + fill_missing_ids = TRUE, + detailed_results = FALSE +) { + total_constrained( + travel_matrix, + land_use_data, + opportunity, + travel_cost, + decay_function, + group_by, + fill_missing_ids, + detailed_results + ) +} + +test_that("raises errors due to incorrect input", { + expect_error(tester(decay_function = "a")) + expect_error(tester(decay_function = mean)) + expect_error(tester(decay_function = get)) + + expect_error(tester(opportunity = 1)) + expect_error(tester(opportunity = c("schools", "jobs"))) + + expect_error(tester(travel_cost = 1)) + expect_error(tester(travel_cost = c("travel_time", "monetary_cost"))) + + expect_error(tester(group_by = 1)) + expect_error(tester(group_by = NA)) + expect_error(tester(group_by = "from_id")) + expect_error(tester(group_by = c("mode", "mode"))) + + expect_error(tester(fill_missing_ids = 1)) + expect_error(tester(fill_missing_ids = c(TRUE, TRUE))) + expect_error(tester(fill_missing_ids = NA)) + + expect_error(tester(detailed_results = 1)) + expect_error(tester(detailed_results = c(TRUE, TRUE))) + expect_error(tester(detailed_results = NA)) + + expect_error(tester(as.list(travel_matrix))) + expect_error(tester(travel_matrix[, .(oi = from_id, to_id, travel_time)])) + expect_error(tester(travel_matrix[, .(from_id, oi = to_id, travel_time)])) + expect_error( + tester( + travel_matrix[, .(from_id, to_id, oi = travel_time)], + travel_cost = "travel_time" + ) + ) + expect_error( + tester( + travel_matrix[, .(from_id, to_id, travel_time, oi = mode)], + group_by = "mode" + ) + ) + + expect_error(tester(as.list(land_use_data))) + expect_error( + tester(land_use_data = land_use_data[, .(oi = id, jobs)]) + ) + expect_error( + tester( + land_use_data = land_use_data[, .(id, oi = jobs)], + opportunity = "jobs" + ) + ) +}) + +test_that("throws warning if travel_matrix has extra col", { + # i.e. col not listed in travel_cost and by_col + expect_warning(tester(group_by = character(0))) +}) + +test_that("returns a dataframe whose class is the same as travel_matrix's", { + result <- tester() + expect_is(result, "data.table") + result <- tester(land_use_data = as.data.frame(land_use_data)) + expect_is(result, "data.table") + + result <- tester(as.data.frame(smaller_matrix)) + expect_false(inherits(result, "data.table")) + expect_is(result, "data.frame") + result <- tester(as.data.frame(smaller_matrix), as.data.frame(land_use_data)) + expect_false(inherits(result, "data.table")) + expect_is(result, "data.frame") +}) + +test_that("result has correct structure", { + result <- tester() + expect_true(ncol(result) == 3) + expect_is(result$id, "character") + expect_is(result$mode, "character") + expect_is(result$jobs, "numeric") + + result <- tester(opportunity = "schools") + expect_true(ncol(result) == 3) + expect_is(result$id, "character") + expect_is(result$mode, "character") + expect_is(result$schools, "numeric") + + suppressWarnings(result <- tester(group_by = character(0))) + expect_true(ncol(result) == 2) + expect_is(result$id, "character") + expect_is(result$jobs, "numeric") + + result <- tester( + data.table::data.table( + mode = character(), + from_id = character(), + to_id = character(), + travel_time = integer() + ) + ) + expect_true(ncol(result) == 3) + expect_true(nrow(result) == 0) + expect_is(result$id, "character") + expect_is(result$mode, "character") + expect_is(result$jobs, "numeric") +}) + +test_that("input data sets remain unchanged", { + original_smaller_matrix <- travel_matrix[1:10] + original_land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) + + result <- tester() + + # Reset all indices for comparison + data.table::setindex(smaller_matrix, NULL) + data.table::setindex(land_use_data, NULL) + data.table::setindex(original_smaller_matrix, NULL) + data.table::setindex(original_land_use_data, NULL) + + expect_equal(original_smaller_matrix, smaller_matrix) + expect_equal(original_land_use_data, land_use_data) +}) + +test_that("fill_missing_ids arg works correctly", { + small_travel_matrix <- travel_matrix[ + from_id %in% c("89a88cdb57bffff", "89a88cdb597ffff") + ] + small_travel_matrix <- small_travel_matrix[ + !(from_id == "89a88cdb57bffff" & mode == "transit2") + ] + + result <- tester(small_travel_matrix, fill_missing_ids = TRUE) + result[, jobs := as.integer(jobs)] + data.table::setkey(result, NULL) + expect_identical( + result, + data.table::data.table( + id = rep(c("89a88cdb57bffff", "89a88cdb597ffff"), each = 2), + mode = rep(c("transit", "transit2"), times = 2), + jobs = c(164981L, 0L, 165553L, 165553L) # Update these expected values + ) + ) + + result <- tester(small_travel_matrix, fill_missing_ids = FALSE) + result[, jobs := as.integer(jobs)] + expect_identical( + result, + data.table::data.table( + id = c("89a88cdb57bffff", "89a88cdb597ffff", "89a88cdb597ffff"), + mode = c("transit", "transit", "transit2"), + jobs = c(164981L, 165553L, 165553L) # Update these expected values + ) + ) +}) + +test_that("accepts custom decay function", { + selected_ids <- c( + "89a88cdb57bffff", + "89a88cdb597ffff", + "89a88cdb5b3ffff", + "89a88cdb5cfffff", + "89a88cd909bffff" + ) + smaller_travel_matrix <- travel_matrix[ + from_id %in% selected_ids & to_id %in% selected_ids + ] + + custom_function <- function(travel_cost) rep(1L, length(travel_cost)) + + result <- tester(smaller_travel_matrix, decay_function = custom_function) + result[, jobs := round(jobs, digits = 1)] #for comparision purposes + expect_identical( + result, + data.table::data.table( + id = rep(selected_ids, 2), + mode = rep(c("transit", "transit2"), each = 5), + jobs = rep(49608.8, 10) # All values are the same + ) + ) +}) + +test_that("calculates total constrained accessibility correctly", { + # Use the same test data from spatial_availability tests + paper_travel_matrix <- data.table::data.table( + from_id = rep(c("A", "B", "C"), each = 3), + to_id = as.character(rep(1:3, 3)), + travel_time = c(15, 30, 100, 30, 15, 100, 100, 100, 15) + ) + paper_land_use_data <- data.table::data.table( + id = c("A", "B", "C", "1", "2", "3"), + jobs = c(0, 0, 0, 100000, 100000, 10000) + ) + + result <- tester( + paper_travel_matrix, + paper_land_use_data, + group_by = character(0) + ) + #result[, jobs := round(jobs, digits = 0)] + + # Test conservation property - sum should equal total opportunities + total_jobs <- sum(paper_land_use_data[id %in% c("1", "2", "3"), jobs]) + expect_equal(sum(result$jobs), total_jobs)#, tolerance = 1) + + # Test that all origins are present + expect_setequal(result$id, c("A", "B", "C")) +}) + +test_that("results are grouped by decay_function_arg when needed", { + small_travel_matrix <- travel_matrix[ + from_id %in% c("89a88cdb57bffff", "89a88cdb597ffff") & + mode != "transit2" + ] + + result <- tester( + small_travel_matrix, + decay_function = decay_exponential(c(0.5, 0.6)) + ) + #result[, jobs := round(jobs, 1)] + + expect_true("decay_function_arg" %in% names(result)) + expect_equal(unique(result$decay_function_arg), c(0.5, 0.6)) +}) + +test_that("throws warning w/ fill_missing_ids = FALSE with detailed_results", { + expect_warning(tester(fill_missing_ids = FALSE, detailed_results = TRUE)) +}) + +test_that("result has correct structure with detailed_results = TRUE", { + result <- tester(detailed_results = TRUE) + expect_true(ncol(result) >= 5) # from_id, to_id, group_cols, kappa_tot, jobs + expect_is(result$mode, "character") + expect_is(result$from_id, "character") + expect_is(result$to_id, "character") + expect_is(result$kappa_tot, "numeric") + expect_is(result$jobs, "numeric") + + result <- tester(detailed_results = TRUE, opportunity = "schools") + expect_true(ncol(result) >= 5) + expect_is(result$mode, "character") + expect_is(result$from_id, "character") + expect_is(result$to_id, "character") + expect_is(result$kappa_tot, "numeric") + expect_is(result$schools, "numeric") + + result <- tester(smaller_matrix[0], detailed_results = TRUE) + expect_true(ncol(result) >= 5) + expect_true(nrow(result) == 0) + expect_is(result$mode, "character") + expect_is(result$from_id, "character") + expect_is(result$to_id, "character") + expect_is(result$kappa_tot, "numeric") + expect_is(result$jobs, "numeric") +}) From 9ab7d14640316def233099fc6028f6c8d94edc83 Mon Sep 17 00:00:00 2001 From: soukhova Date: Thu, 23 Oct 2025 18:01:16 -0400 Subject: [PATCH 2/4] changes to other documents, when testing and checking the package. --- DESCRIPTION | 2 +- NAMESPACE | 1 + man/accessibility.Rd | 1 - man/balancing_cost.Rd | 2 +- man/concentration_index.Rd | 2 +- man/cost_to_closest.Rd | 2 +- man/cumulative_cutoff.Rd | 2 +- man/cumulative_interval.Rd | 2 +- man/decay_binary.Rd | 2 +- man/decay_exponential.Rd | 2 +- man/decay_linear.Rd | 2 +- man/decay_logistic.Rd | 2 +- man/decay_power.Rd | 2 +- man/decay_stepped.Rd | 2 +- man/fgt_poverty.Rd | 2 +- man/floating_catchment_area.Rd | 2 +- man/gini_index.Rd | 2 +- man/gravity.Rd | 2 +- man/palma_ratio.Rd | 2 +- man/spatial_availability.Rd | 2 +- man/theil_t.Rd | 2 +- 21 files changed, 20 insertions(+), 20 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 79a198a..2b71fe8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -48,4 +48,4 @@ RdMacros: Encoding: UTF-8 NeedsCompilation: no Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.3 diff --git a/NAMESPACE b/NAMESPACE index 3d55ae3..b99153c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -18,6 +18,7 @@ export(gravity) export(palma_ratio) export(spatial_availability) export(theil_t) +export(total_constrained) importFrom(Rdpack,reprompt) importFrom(data.table,"%chin%") importFrom(data.table,":=") diff --git a/man/accessibility.Rd b/man/accessibility.Rd index a657e17..b74e655 100644 --- a/man/accessibility.Rd +++ b/man/accessibility.Rd @@ -3,7 +3,6 @@ \docType{package} \name{accessibility} \alias{accessibility} -\alias{_PACKAGE} \alias{accessibility-package} \title{accessibility: Transport accessibility measures} \description{ diff --git a/man/balancing_cost.Rd b/man/balancing_cost.Rd index e1b826e..c1b224f 100644 --- a/man/balancing_cost.Rd +++ b/man/balancing_cost.Rd @@ -85,7 +85,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/concentration_index.Rd b/man/concentration_index.Rd index bbdedda..fc419ee 100644 --- a/man/concentration_index.Rd +++ b/man/concentration_index.Rd @@ -69,7 +69,7 @@ calculating the standard relative CI and the corrected CI, as proposed by \insertCite{erreygers2009correcting;textual}{accessibility}. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/cost_to_closest.Rd b/man/cost_to_closest.Rd index 9447b01..1d55599 100644 --- a/man/cost_to_closest.Rd +++ b/man/cost_to_closest.Rd @@ -70,7 +70,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/cumulative_cutoff.Rd b/man/cumulative_cutoff.Rd index 134ba56..c417de0 100644 --- a/man/cumulative_cutoff.Rd +++ b/man/cumulative_cutoff.Rd @@ -80,7 +80,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/cumulative_interval.Rd b/man/cumulative_interval.Rd index 9b31536..645deb6 100644 --- a/man/cumulative_interval.Rd +++ b/man/cumulative_interval.Rd @@ -76,7 +76,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/decay_binary.Rd b/man/decay_binary.Rd index 3ab3d88..95e91ff 100644 --- a/man/decay_binary.Rd +++ b/man/decay_binary.Rd @@ -24,7 +24,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} weighting_function <- decay_binary(cutoff = 30) weighting_function(c(20, 35)) diff --git a/man/decay_exponential.Rd b/man/decay_exponential.Rd index 3f6a882..efd2f0a 100644 --- a/man/decay_exponential.Rd +++ b/man/decay_exponential.Rd @@ -24,7 +24,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} weighting_function <- decay_exponential(decay_value = 0.1) weighting_function(c(20, 30)) diff --git a/man/decay_linear.Rd b/man/decay_linear.Rd index d02ddb2..f4a083d 100644 --- a/man/decay_linear.Rd +++ b/man/decay_linear.Rd @@ -24,7 +24,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} weighting_function <- decay_linear(cutoff = 30) weighting_function(c(20, 35)) diff --git a/man/decay_logistic.Rd b/man/decay_logistic.Rd index 00e517a..85e21be 100644 --- a/man/decay_logistic.Rd +++ b/man/decay_logistic.Rd @@ -41,7 +41,7 @@ condition f(0) = 1 is met (i.e. the weight of an opportunity whose cost to reach is 0 is 1). } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} weighting_function <- decay_logistic(cutoff = 30, sd = 5) weighting_function(c(0, 30, 45, 60)) diff --git a/man/decay_power.Rd b/man/decay_power.Rd index cc0490a..ba47beb 100644 --- a/man/decay_power.Rd +++ b/man/decay_power.Rd @@ -23,7 +23,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} weighting_function <- decay_power(decay_value = 0.1) weighting_function(c(20, 35)) diff --git a/man/decay_stepped.Rd b/man/decay_stepped.Rd index 13d5cf4..d8d71b4 100644 --- a/man/decay_stepped.Rd +++ b/man/decay_stepped.Rd @@ -44,7 +44,7 @@ the output will be named \code{"s(10,20,30);w(0.66,0.33,0)"}. } } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} weighting_function <- decay_stepped( c(10, 20, 30, 40), weights = c(0.75, 0.5, 0.25, 0) diff --git a/man/fgt_poverty.Rd b/man/fgt_poverty.Rd index 7927eb8..fbb0916 100644 --- a/man/fgt_poverty.Rd +++ b/man/fgt_poverty.Rd @@ -76,7 +76,7 @@ however, FGT0 value is 1 and FGT1 and FGT2 values approach 1. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/floating_catchment_area.Rd b/man/floating_catchment_area.Rd index 0d0aeb5..29b1644 100644 --- a/man/floating_catchment_area.Rd +++ b/man/floating_catchment_area.Rd @@ -96,7 +96,7 @@ present in other FCA measures. It was originally proposed by } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/gini_index.Rd b/man/gini_index.Rd index 3f832d8..0e9618a 100644 --- a/man/gini_index.Rd +++ b/man/gini_index.Rd @@ -44,7 +44,7 @@ A data frame containing the inequality estimates for the study area. Calculates the Gini Index of a given accessibility distribution. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/gravity.Rd b/man/gravity.Rd index 340fd0d..fbd14cd 100644 --- a/man/gravity.Rd +++ b/man/gravity.Rd @@ -74,7 +74,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/palma_ratio.Rd b/man/palma_ratio.Rd index 4ed0e70..9108338 100644 --- a/man/palma_ratio.Rd +++ b/man/palma_ratio.Rd @@ -56,7 +56,7 @@ planning as the average accessibility of the richest 10\% divided by the average accessibility of the poorest 40\%. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) diff --git a/man/spatial_availability.Rd b/man/spatial_availability.Rd index 22059cd..9f744d3 100644 --- a/man/spatial_availability.Rd +++ b/man/spatial_availability.Rd @@ -91,7 +91,7 @@ This function is generic over any kind of numeric travel cost, such as distance, time and money. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} # the example below is based on Soukhov et al. (2023) paper travel_matrix <- data.table::data.table( diff --git a/man/theil_t.Rd b/man/theil_t.Rd index bbb8c5f..b3987a2 100644 --- a/man/theil_t.Rd +++ b/man/theil_t.Rd @@ -63,7 +63,7 @@ decomposed into a between-groups inequaliy component and a within-groups component. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) (if (getRversion() >= "3.4") withAutoprint else force)(\{ # examplesIf} +\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} data_dir <- system.file("extdata", package = "accessibility") travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) From c6d29f3703f981fbc24e20c76911d651bdd89392 Mon Sep 17 00:00:00 2001 From: soukhova Date: Thu, 4 Dec 2025 19:01:33 -0500 Subject: [PATCH 3/4] add-constrained-accessibility-wrapper --- .Rbuildignore | 1 + NAMESPACE | 2 +- R/constrained_accessibility.R | 144 ++++++++ R/doubly_constrained.R | 148 ++++++++ R/singly_constrained.R | 104 ++++++ R/total_constrained.R | 214 +++++------- man/balancing_cost.Rd | 4 +- man/constrained_accessibility.Rd | 140 ++++++++ man/doubly_constrained.Rd | 19 ++ man/floating_catchment_area.Rd | 4 +- man/roxygen/templates/demand.R | 4 +- man/roxygen/templates/supply.R | 3 + man/singly_constrained.Rd | 16 + man/spatial_availability.Rd | 4 +- man/total_constrained.Rd | 118 +------ testing_constrained.Rmd | 337 +++++++++++++++++++ tests/testthat/test-constrained-essentials.R | 151 +++++++++ tests/testthat/test-constrained-wrapper.R | 45 +++ tests/testthat/test-total_constrained.R | 276 --------------- 19 files changed, 1214 insertions(+), 520 deletions(-) create mode 100644 R/constrained_accessibility.R create mode 100644 R/doubly_constrained.R create mode 100644 R/singly_constrained.R create mode 100644 man/constrained_accessibility.Rd create mode 100644 man/doubly_constrained.Rd create mode 100644 man/roxygen/templates/supply.R create mode 100644 man/singly_constrained.Rd create mode 100644 testing_constrained.Rmd create mode 100644 tests/testthat/test-constrained-essentials.R create mode 100644 tests/testthat/test-constrained-wrapper.R delete mode 100644 tests/testthat/test-total_constrained.R diff --git a/.Rbuildignore b/.Rbuildignore index 3a804fd..1d906a7 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -11,3 +11,4 @@ ^codemeta\.json$ ^data-raw$ ^\.pre-commit-config\.yaml$ +^testing_constrained\.Rmd$ diff --git a/NAMESPACE b/NAMESPACE index b99153c..b2f77f4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,6 +2,7 @@ export(balancing_cost) export(concentration_index) +export(constrained_accessibility) export(cost_to_closest) export(cumulative_cutoff) export(cumulative_interval) @@ -18,7 +19,6 @@ export(gravity) export(palma_ratio) export(spatial_availability) export(theil_t) -export(total_constrained) importFrom(Rdpack,reprompt) importFrom(data.table,"%chin%") importFrom(data.table,":=") diff --git a/R/constrained_accessibility.R b/R/constrained_accessibility.R new file mode 100644 index 0000000..d868d62 --- /dev/null +++ b/R/constrained_accessibility.R @@ -0,0 +1,144 @@ +#' Constrained accessibility +#' +#' Calculates accessibility using constrained gravity models: +#' - `"total"`: Allocates total opportunities proportionally based on travel impedance. +#' - `"singly"`: Allocates opportunities proportionally, constrained on one side (demand or supply). +#' - `"doubly"`: Allocates flows so origin totals equal demand and destination totals equal supply. +#' +#' @template description_generic_cost +#' @template travel_matrix +#' @template land_use_data +#' @template travel_cost +#' @template decay_function +#' @template demand +#' @template supply +#' @param constraint A string. One of `"total"`, `"singly"`, or `"doubly"`. +#' @param return_demand_side Logical for `"total"` and `"singly"`, must be `NULL` for `"doubly"`. +#' @param error_threshold Numeric. Convergence criterion for doubly-constrained case. +#' @param improvement_threshold Numeric. Convergence criterion for improvement. +#' @param max_iterations Integer. Maximum iterations for doubly-constrained calibration. +#' @template group_by +#' @template fill_missing_ids_combinations +#' @param detailed_results Logical. Whether to return detailed OD-level results. +#' +#' @details +#' See individual function documentation for mathematical details: +#' [total_constrained()], [singly_constrained()], [doubly_constrained()]. +#' +#' @family Constrained accessibility +#' +#' @examples +#' # Load demo data shipped with the package (used for 'total' and 'singly') +#' data_dir <- system.file("extdata", package = "accessibility") +#' travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) +#' land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) +#' +#' # Total-constrained (supply-side) +#' constrained_accessibility("total", travel_matrix, land_use_data, +#' travel_cost = "travel_time", +#' decay_function = decay_exponential(0.1), +#' demand = NULL, +#' supply = "jobs", +#' return_demand_side = FALSE +#' ) +#' +#' # Singly-constrained (demand-side) +#' constrained_accessibility("singly", travel_matrix, land_use_data, +#' travel_cost = "travel_time", +#' decay_function = decay_exponential(0.1), +#' demand = "population", +#' supply = "jobs", +#' return_demand_side = TRUE +#' ) +#' +#' # Doubly-constrained: use a small toy dataset with matching totals +#' tm_small <- data.table::data.table( +#' expand.grid(from_id = c("1","2","3"), to_id = c("1","2","3")) +#' ) +#' tm_small[, travel_time := c(10, 30, 15, 30, 10, 25, 15, 25, 10)] +#' lu_small <- data.table::data.table( +#' id = c("1","2","3"), +#' population = c(4, 10, 6), # sum = 20 +#' jobs = c(7, 5, 8) # sum = 20 +#' ) +#' +#' constrained_accessibility("doubly", tm_small, lu_small, +#' travel_cost = "travel_time", +#' decay_function = decay_exponential(0.1), +#' demand = "population", +#' supply = "jobs", +#' return_demand_side = NULL +#' ) +#' @export +constrained_accessibility <- function(constraint, + travel_matrix, + land_use_data, + travel_cost, + decay_function, + demand, + supply, + return_demand_side = NULL, + error_threshold = 0.001, + improvement_threshold = 1e-6, + max_iterations = 1000, + group_by = character(0), + fill_missing_ids = TRUE, + detailed_results = FALSE) { + + checkmate::assert_choice(constraint, c("total", "singly", "doubly")) + + + if (constraint == "doubly") { + if (!is.null(return_demand_side)) { + stop("For 'doubly', return_demand_side must be NULL.") + } + return(doubly_constrained( + travel_matrix = travel_matrix, + land_use_data = land_use_data, + travel_cost = travel_cost, + decay_function = decay_function, + demand = demand, + supply = supply, + error_threshold = error_threshold, + improvement_threshold = improvement_threshold, + max_iterations = max_iterations, + group_by = group_by, + fill_missing_ids = fill_missing_ids, + detailed_results = detailed_results + )) + } + + if (is.null(return_demand_side)) { + stop(sprintf("For '%s', return_demand_side must be TRUE or FALSE.", constraint)) + } + + if (constraint == "total") { + return(total_constrained( + travel_matrix = travel_matrix, + land_use_data = land_use_data, + travel_cost = travel_cost, + decay_function = decay_function, + group_by = group_by, + fill_missing_ids = fill_missing_ids, + detailed_results = detailed_results, + return_demand_side = return_demand_side, + demand = if (return_demand_side) demand else NULL, + supply = if (!return_demand_side) supply else NULL + )) + } + + if (constraint == "singly") { + return(singly_constrained( + travel_matrix = travel_matrix, + land_use_data = land_use_data, + travel_cost = travel_cost, + decay_function = decay_function, + demand = demand, + supply = supply, + return_demand_side = return_demand_side, + group_by = group_by, + fill_missing_ids = fill_missing_ids, + detailed_results = detailed_results + )) + } +} diff --git a/R/doubly_constrained.R b/R/doubly_constrained.R new file mode 100644 index 0000000..1fc71ee --- /dev/null +++ b/R/doubly_constrained.R @@ -0,0 +1,148 @@ +#' Doubly constrained accessibility +#' +#' Calculates accessibility using Wilson's doubly-constrained gravity model. +#' This measure allocates flows between origins and destinations such that +#' origin totals equal demand and destination totals equal supply. +#' Internal helper used by [constrained_accessibility()] when `constraint = "doubly"`. +#' +#' @name doubly_constrained +#' @keywords internal +#' @return A `data.table`/`data.frame` with either OD-level flows (`detailed_results = TRUE`) +#' or marginals; see Details. +#' @examples NULL +#' @importFrom utils globalVariables +utils::globalVariables(c( + # common ids / columns + "from_id", "to_id", "id", "opp_weight", "supply", "demand", "weighted_demand", "weighted_supply", + # total-constrained + "kappa_total", "hatkappa_total", "K_total", "hatK_total", "constrained_opportunity", + # singly-constrained + "denom_i", "denom_j", "kappa_singly", "hatkappa_singly", "A_i", "B_j", "singly_access", + # doubly-constrained + "kappa_doubly", "flow", "error" +)) +doubly_constrained <- function(travel_matrix, + land_use_data, + travel_cost, + decay_function, + demand, + supply, + return_demand_side = NULL, + error_threshold = 0.001, + improvement_threshold = 1e-6, + max_iterations = 1000, + group_by = character(0), + fill_missing_ids = TRUE, + detailed_results = FALSE) { + # Validate inputs + checkmate::assert_string(travel_cost) + checkmate::assert_string(demand) + checkmate::assert_string(supply) + checkmate::assert_null(return_demand_side) + checkmate::assert_number(error_threshold, lower = 0) + checkmate::assert_number(improvement_threshold, lower = 0) + checkmate::assert_int(max_iterations, lower = 1) + checkmate::assert_logical(detailed_results, len = 1) + assert_decay_function(decay_function) + assert_group_by(group_by) + assert_detailed_fill_missing_ids(fill_missing_ids, detailed_results) + assert_travel_matrix(travel_matrix, travel_cost, group_by) + assert_land_use_data(land_use_data, travel_matrix, opportunity = supply, demand = demand) + + # Convert to data.table + if (!inherits(travel_matrix, "data.table")) { + original_class <- class(travel_matrix) + data <- data.table::as.data.table(travel_matrix) + } else { + data <- data.table::copy(travel_matrix) + } + if (!inherits(land_use_data, "data.table")) { + land_use_data <- data.table::as.data.table(land_use_data) + } + + # Merge demand and supply + merge_by_reference(data, land_use_data, demand, left_df_idcol = "from_id") + merge_by_reference(data, land_use_data, supply, left_df_idcol = "to_id") + + # Apply decay + data <- apply_gravity_measure(data, decay_function, travel_cost) + + # Prepare vectors and impedance matrix + origins <- unique(data$from_id) + destinations <- unique(data$to_id) + O <- land_use_data[id %in% origins, get(demand)] + D <- land_use_data[id %in% destinations, get(supply)] + + # Validate totals match + if (abs(sum(O) - sum(D)) > 1e-6) { + stop("For doubly-constrained, the sum of origins must equal the sum of destinations.") + } + + # Build impedance matrix + impedance_matrix <- matrix(data$opp_weight, nrow = length(origins), ncol = length(destinations), byrow = FALSE) + + # Iterative proportional fitting + Ai <- rep(1, length(O)) + Bj <- rep(1, length(D)) + previous_error <- Inf + iteration_count <- 0 + stop_reason <- "" + + repeat { + iteration_count <- iteration_count + 1 + + # Update Ai + for (i in seq_along(O)) { + Ai[i] <- 1 / (sum(Bj * D * impedance_matrix[i, ]) + 1e-9) + } + + # Update Bj + Bj_new <- numeric(length(D)) + for (j in seq_along(D)) { + Bj_new[j] <- 1 / (sum(Ai * O * impedance_matrix[, j]) + 1e-9) + } + + # Compute flows + Tij <- outer(Ai * O, Bj_new * D) * impedance_matrix + + # Compute error + error <- (sum(abs(O - rowSums(Tij))) + sum(abs(D - colSums(Tij)))) / sum(O) + error_change <- abs(previous_error - error) + + if (error < error_threshold || error_change < improvement_threshold || iteration_count >= max_iterations) { + Bj <- Bj_new + stop_reason <- if (iteration_count >= max_iterations) "Max iterations reached" else if (error < error_threshold) "Error threshold met" else "Slow improvement" + break + } + + previous_error <- error + Bj <- Bj_new + } + + # Compute kappa_doubly +# kappa_ij^D = A_i * B_j * O_i * f(c_ij) + kappa_matrix <- outer(Ai * O, Bj) * impedance_matrix + + # Prepare output + detailed_dt <- data.table::data.table( + from_id = rep(origins, each = length(destinations)), + to_id = rep(destinations, times = length(origins)), + flow = as.vector(t(Tij)), + A_i = rep(Ai, each = length(destinations)), + B_j = rep(Bj, times = length(origins)), + kappa_doubly = as.vector(t(kappa_matrix)), + error = error + ) + + if (detailed_results) { + return(detailed_dt[]) + } else { + warning("Aggregated results equal marginals (O_i or D_j). Interpret with caution.") + agg <- data.table::rbindlist(list( + data.table::data.table(id = origins, flow = O), + data.table::data.table(id = destinations, flow = D) + )) + return(agg[]) + } + +} diff --git a/R/singly_constrained.R b/R/singly_constrained.R new file mode 100644 index 0000000..1e07b17 --- /dev/null +++ b/R/singly_constrained.R @@ -0,0 +1,104 @@ +#' Singly constrained accessibility +#' +#' Allocates opportunities at each destination proportionally based on travel impedance and population at the origin. Uses the logic of Wilon's single constraint. Returns values as either 'demand' or 'supply'. +#' Internal helper used by [constrained_accessibility()] when `constraint = "singly"`. +#' +#' @name singly_constrained +#' @keywords internal +#' @return A `data.table`/`data.frame` with results (structure mirrors the wrapper). +#' @examples NULL +#' @importFrom utils globalVariables +utils::globalVariables(c( + # common ids / columns + "from_id", "to_id", "id", "opp_weight", "supply", "demand", "weighted_demand", "weighted_supply", + # total-constrained + "kappa_total", "hatkappa_total", "K_total", "hatK_total", "constrained_opportunity", + # singly-constrained + "denom_i", "denom_j", "kappa_singly", "hatkappa_singly", "A_i", "B_j", "singly_access", + # doubly-constrained + "kappa_doubly", "flow", "error" +)) +singly_constrained <- function(travel_matrix, + land_use_data, + travel_cost, + decay_function, + demand, + supply, + return_demand_side = FALSE, + group_by = character(0), + fill_missing_ids = TRUE, + detailed_results = FALSE) { + # Validate + checkmate::assert_string(travel_cost) + checkmate::assert_string(demand) + checkmate::assert_string(supply) + checkmate::assert_logical(return_demand_side, len = 1) + checkmate::assert_logical(detailed_results, len = 1) + assert_decay_function(decay_function) + assert_group_by(group_by) + assert_detailed_fill_missing_ids(fill_missing_ids, detailed_results) + assert_travel_matrix(travel_matrix, travel_cost, group_by) + assert_land_use_data(land_use_data, travel_matrix, opportunity = supply, demand = demand) + + # Convert to data.table + if (!inherits(travel_matrix, "data.table")) { + original_class <- class(travel_matrix) + data <- data.table::as.data.table(travel_matrix) + } else { + data <- data.table::copy(travel_matrix) + } + if (!inherits(land_use_data, "data.table")) { + land_use_data <- data.table::as.data.table(land_use_data) + } + + # Merge demand and supply + merge_by_reference(data, land_use_data, demand, left_df_idcol = "from_id") + merge_by_reference(data, land_use_data, supply, left_df_idcol = "to_id") + + # Apply decay + data <- apply_gravity_measure(data, decay_function, travel_cost) + + if (!return_demand_side) { + # Supply-constrained (returns the number of supply or JOBS) + data[, weighted_demand := get(demand) * opp_weight] + data[, denom_j := sum(weighted_demand), by = c("to_id", group_by)] + data[, kappa_singly := weighted_demand / denom_j] + data[, singly_access := kappa_singly * get(supply)] + + if (detailed_results) { + access <- data[, .( + from_id, + to_id, + kappa_singly, + supply = singly_access, + B_j = 1 / denom_j + )] + } else { + access <- data[, .(supply = sum(singly_access)), by = c("from_id", group_by)] + if (fill_missing_ids) access <- fill_missing_ids(access, travel_matrix, c("from_id", group_by)) + } + + } else { + # Demand-constrained (returns the number of demand or POPULATION) + data[, weighted_supply := get(supply) * opp_weight] + data[, denom_i := sum(weighted_supply), by = c("from_id", group_by)] + data[, hatkappa_singly := weighted_supply / denom_i] + data[, singly_access := hatkappa_singly * get(demand)] + + if (detailed_results) { + access <- data[, .( + from_id, + to_id, + hatkappa_singly, + demand = singly_access, + A_i = 1 / denom_i + )] + } else { + access <- data[, .(demand = sum(singly_access)), by = c("to_id", group_by)] + if (fill_missing_ids) access <- fill_missing_ids(access, travel_matrix, c("to_id", group_by)) + } + } + + if (exists("original_class")) class(access) <- original_class + return(access[]) +} diff --git a/R/total_constrained.R b/R/total_constrained.R index d3f8c95..be424d8 100644 --- a/R/total_constrained.R +++ b/R/total_constrained.R @@ -1,156 +1,124 @@ -#' Total constrained +#' Total constrained accessibility #' -#' Calculates total constrained accessibility, part of a family of accessibility -#'measures proposed in -#'\insertCite{soukhovFamilyAccessibilityMeasures2025;textual}{accessibility} that -#' proportionally allocates the total sum of opportunities in the region based -#' on the relative travel impedance between zones. It is linearly proportion to -#' the conventional gravity-based Hansen-type accessibility measure and units of -#' the resulting values can be understood as the number of opportunities in the -#' region that can be potentially interacted with. -#' @template description_generic_cost +#' Allocates total opportunities in the region proportionally based on travel impedance. Uses the logic of a total (or unconstrained by Wilon's terms) constraint. Returns values as either 'demand' or 'supply' +#' Internal helper used by [constrained_accessibility()] when `constraint = "total"`. +#' @name total_constrained +#' @keywords internal +#' @return A `data.table`/`data.frame` with results (structure mirrors the wrapper). +#' @examples NULL #' -#' @template travel_matrix -#' @template land_use_data -#' @template opportunity -#' @template travel_cost -#' @template decay_function -#' @template group_by -#' @template fill_missing_ids_combinations -#' @param detailed_results A `logical`. Whether to return total constrained -#' results aggregated by origin-destination pair (`TRUE`) or by origin -#' (`FALSE`, the default). When `TRUE`, the output also includes the balancing -#' factor `kappa_tot` used in the calculation and the unconstrained -#' accessibility value `weighted_opportunity`. Please note that the argument -#' `fill_missing_ids` does not affect the output when `detailed_results` is -#' `TRUE`. +#' @importFrom utils globalVariables #' -#' @template return_accessibility -#' -#' @references -#' \insertAllCited{} -#' -#' @examplesIf identical(tolower(Sys.getenv("NOT_CRAN")), "true") -#' # the example below is based on Soukhov et al. (2023) paper -#' -#' travel_matrix <- data.table::data.table( -#' from_id = rep(c("A", "B", "C"), each = 3), -#' to_id = as.character(rep(1:3, 3)), -#' travel_time = c(15, 30, 100, 30, 15, 100, 100, 100, 15) -#' ) -#' land_use_data <- data.table::data.table( -#' id = c("A", "B", "C", "1", "2", "3"), -#' population = c(50000, 150000, 10000, 0, 0, 0), -#' jobs = c(0, 0, 0, 100000, 100000, 10000) -#' ) -#' -#' df <- total_constrained( -#' travel_matrix, -#' land_use_data, -#' opportunity = "jobs", -#' travel_cost = "travel_time", -#' decay_function = decay_exponential(decay_value = 0.1) -#' ) -#' df -#' -#' detailed_df <- total_constrained( -#' travel_matrix, -#' land_use_data, -#' opportunity = "jobs", -#' travel_cost = "travel_time", -#' decay_function = decay_exponential(decay_value = 0.1), -#' detailed_results = TRUE -#' ) -#' detailed_df -#' -#' @export +utils::globalVariables(c( + # common ids / columns + "from_id", "to_id", "id", "opp_weight", "supply", "demand", "weighted_demand", "weighted_supply", + # total-constrained + "kappa_total", "hatkappa_total", "K_total", "hatK_total", "constrained_opportunity", + # singly-constrained + "denom_i", "denom_j", "kappa_singly", "hatkappa_singly", "A_i", "B_j", "singly_access", + # doubly-constrained + "kappa_doubly", "flow", "error" +)) total_constrained <- function(travel_matrix, - land_use_data, - opportunity, - travel_cost, - decay_function, - group_by = character(0), - fill_missing_ids = TRUE, - detailed_results = FALSE) { - checkmate::assert_string(opportunity) + land_use_data, + travel_cost, + decay_function, + group_by = character(0), + fill_missing_ids = TRUE, + detailed_results = FALSE, + return_demand_side = FALSE, + demand = NULL, # population + supply = NULL) { # jobs + + # Validate inputs checkmate::assert_string(travel_cost) checkmate::assert_logical(detailed_results, len = 1, any.missing = FALSE) + checkmate::assert_logical(return_demand_side, len = 1, any.missing = FALSE) + + if (return_demand_side) { + if (is.null(demand) || !is.null(supply)) { + stop("For return_demand_side = TRUE, demand must be specified and supply must be NULL.") + } + merge_id <- "from_id" + group_id <- "to_id" + weighted_col <- "weighted_demand" + kappa_col <- "hatkappa_total" + total_col <- "hatK_total" + result_col <- "demand" + opportunity_col <- demand + } else { + if (is.null(supply) || !is.null(demand)) { + stop("For return_demand_side = FALSE, supply must be specified and demand must be NULL.") + } + merge_id <- "to_id" + group_id <- "from_id" + weighted_col <- "weighted_supply" + kappa_col <- "kappa_total" + total_col <- "K_total" + result_col <- "supply" + opportunity_col <- supply + } + assert_decay_function(decay_function) assert_group_by(group_by) assert_detailed_fill_missing_ids(fill_missing_ids, detailed_results) assert_travel_matrix(travel_matrix, travel_cost, group_by) - assert_land_use_data( - land_use_data, - travel_matrix, - opportunity) + assert_land_use_data(land_use_data, travel_matrix, opportunity_col) + # Convert to data.table if (!inherits(travel_matrix, "data.table")) { original_class <- class(travel_matrix) data <- data.table::as.data.table(travel_matrix) } else { data <- data.table::copy(travel_matrix) } - if (!inherits(land_use_data, "data.table")) { land_use_data <- data.table::as.data.table(land_use_data) } - merge_by_reference(data, land_use_data, opportunity, left_df_idcol = "to_id") + # Merge land use data + merge_by_reference(data, land_use_data, opportunity_col, left_df_idcol = merge_id) + # Apply decay data <- apply_gravity_measure(data, decay_function, travel_cost) - groups <- c("from_id", group_by) - if ("decay_function_arg" %in% names(data)) { - groups <- c(groups, "decay_function_arg") - group_by <- c(group_by, "decay_function_arg") - } - - warn_extra_cols( - travel_matrix, - travel_cost, - group_id = "from_id", - groups = groups - ) - - # Calculate unconstrained accessibility (weighted_opportunities) per ij: W_j * f(c_ij) - .opportunity_colname <- opportunity - data[, weighted_opportunity := get(.opportunity_colname) * opp_weight] - - # Calculate kappa_tot for each ij: kappa_tot_ij = (W_j * f(c_ij)) / sum_over_all_ij (W_j * f(c_ij)) - total_weighted_system <- data[, sum(weighted_opportunity)] - data[, kappa_tot := weighted_opportunity / total_weighted_system] - - # Total opportunities in the region - total_opportunities_region <- sum(land_use_data[[.opportunity_colname]]) - - # Calculate total constrained opportunity: kappa_tot_ij * total regional opportunities - data[, constrained_opportunity := kappa_tot * total_opportunities_region] + # Core calculation + data[, (weighted_col) := get(opportunity_col) * opp_weight] + total_weighted_system <- data[, sum(get(weighted_col))] + data[, (kappa_col) := get(weighted_col) / total_weighted_system] + total_opportunities_region <- sum(land_use_data[[opportunity_col]]) + data[, constrained_opportunity := get(kappa_col) * total_opportunities_region] if (detailed_results) { - # Return detailed results with all intermediate calculations - data[, total_opportunities_region := total_opportunities_region] - - # Drop unnecessary columns and rename - cols_to_drop <- c(travel_cost, opportunity, "opp_weight", total_opportunities_region) - access <- data[, (cols_to_drop) := NULL] - - data.table::setnames(access, "constrained_opportunity", opportunity) + if (return_demand_side) { + access <- data[, .( + from_id, + to_id, + weighted_demand, + hatkappa_total, + demand = constrained_opportunity, + hatK_total = total_opportunities_region / total_weighted_system + )] + } else { + access <- data[, .( + from_id, + to_id, + weighted_supply, + kappa_total, + supply = constrained_opportunity, + K_total = total_opportunities_region / total_weighted_system + )] + } } else { - # Return aggregated results by origin - access <- data[ - , - .(access = sum(constrained_opportunity)), - by = c("from_id", group_by) - ] - - if (fill_missing_ids) { - access <- fill_missing_ids(access, travel_matrix, groups) + if (return_demand_side) { + access <- data[, .(demand = sum(constrained_opportunity)), by = c("to_id", group_by)] + if (fill_missing_ids) access <- fill_missing_ids(access, travel_matrix, c("to_id", group_by)) + } else { + access <- data[, .(supply = sum(constrained_opportunity)), by = c("from_id", group_by)] + if (fill_missing_ids) access <- fill_missing_ids(access, travel_matrix, c("from_id", group_by)) } - - data.table::setnames(access, c("from_id", "access"), c("id", opportunity)) } if (exists("original_class")) class(access) <- original_class - return(access[]) } diff --git a/man/balancing_cost.Rd b/man/balancing_cost.Rd index c1b224f..1a3042d 100644 --- a/man/balancing_cost.Rd +++ b/man/balancing_cost.Rd @@ -33,8 +33,8 @@ calculating accessibility levels.} with the travel cost between origins and destinations.} \item{demand}{A string. The name of the column in \code{land_use_data} with the -number of people in each origin that will be considered potential -competitors.} +number of opportunity-demanders at each origin (e.g., people) that will be +considered.} \item{cost_increment}{A number. The cost increment that should be used when defining the travel cost distribution from which the potential balancing diff --git a/man/constrained_accessibility.Rd b/man/constrained_accessibility.Rd new file mode 100644 index 0000000..c999955 --- /dev/null +++ b/man/constrained_accessibility.Rd @@ -0,0 +1,140 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/constrained_accessibility.R +\name{constrained_accessibility} +\alias{constrained_accessibility} +\title{Constrained accessibility} +\usage{ +constrained_accessibility( + constraint, + travel_matrix, + land_use_data, + travel_cost, + decay_function, + demand, + supply, + return_demand_side = NULL, + error_threshold = 0.001, + improvement_threshold = 1e-06, + max_iterations = 1000, + group_by = character(0), + fill_missing_ids = TRUE, + detailed_results = FALSE +) +} +\arguments{ +\item{constraint}{A string. One of \code{"total"}, \code{"singly"}, or \code{"doubly"}.} + +\item{travel_matrix}{A data frame. The travel matrix describing the costs +(i.e. travel time, distance, monetary cost, etc.) between the origins and +destinations in the study area. Must contain the columns \code{from_id}, \code{to_id} +and any others specified in \code{travel_cost}.} + +\item{land_use_data}{A data frame. The distribution of opportunities within +the study area cells. Must contain the columns \code{id} and any others +specified in \code{opportunity}.} + +\item{travel_cost}{A string. The name of the column in \code{travel_matrix} +with the travel cost between origins and destinations.} + +\item{decay_function}{A \code{fuction} that converts travel cost into an +impedance factor used to weight opportunities. This function should take a +\code{numeric} vector and also return a \code{numeric} vector as output, with the +same length as the input. For convenience, the package currently includes +the following functions: \code{\link[=decay_binary]{decay_binary()}}, \code{\link[=decay_exponential]{decay_exponential()}}, +\code{\link[=decay_power]{decay_power()}} and \code{\link[=decay_stepped]{decay_stepped()}}. See the documentation of each decay +function for more details.} + +\item{demand}{A string. The name of the column in \code{land_use_data} with the +number of opportunity-demanders at each origin (e.g., people) that will be +considered.} + +\item{supply}{A string. The name of the column in \code{land_use_data} with the +number of opportunity supply at each destination (e.g., jobs, school-seats) +that will be considered.} + +\item{return_demand_side}{Logical for \code{"total"} and \code{"singly"}, must be \code{NULL} for \code{"doubly"}.} + +\item{error_threshold}{Numeric. Convergence criterion for doubly-constrained case.} + +\item{improvement_threshold}{Numeric. Convergence criterion for improvement.} + +\item{max_iterations}{Integer. Maximum iterations for doubly-constrained calibration.} + +\item{group_by}{A \code{character} vector. When not \code{character(0)} (the default), +indicates the \code{travel_matrix} columns that should be used to group the +accessibility estimates by. For example, if \code{travel_matrix} includes a +departure time column, that specifies the departure time of each entry in +the data frame, passing \code{"departure_time"} to this parameter results in +accessibility estimates grouped by origin and by departure time.} + +\item{fill_missing_ids}{A \code{logical}. When calculating grouped accessibility +estimates (i.e. when \code{by_col} is not \code{NULL}), some combinations of groups +and origins may be missing. For example, if a single trip can depart from +origin \code{A} at 7:15am and reach destination \code{B} within 55 minutes, but no +trips departing from \code{A} at 7:30am can be completed at all, this second +combination will not be included in the output. When \code{TRUE} (the default), +the function identifies which combinations would be left out and fills +their respective accessibility values with 0, which incurs in a +performance penalty.} + +\item{detailed_results}{Logical. Whether to return detailed OD-level results.} +} +\description{ +Calculates accessibility using constrained gravity models: +\itemize{ +\item \code{"total"}: Allocates total opportunities proportionally based on travel impedance. +\item \code{"singly"}: Allocates opportunities proportionally, constrained on one side (demand or supply). +\item \code{"doubly"}: Allocates flows so origin totals equal demand and destination totals equal supply. +} + +This function is generic over any kind of numeric travel cost, +such as distance, time and money. +} +\details{ +See individual function documentation for mathematical details: +\code{\link[=total_constrained]{total_constrained()}}, \code{\link[=singly_constrained]{singly_constrained()}}, \code{\link[=doubly_constrained]{doubly_constrained()}}. +} +\examples{ +# Load demo data shipped with the package (used for 'total' and 'singly') +data_dir <- system.file("extdata", package = "accessibility") +travel_matrix <- readRDS(file.path(data_dir, "travel_matrix.rds")) +land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) + +# Total-constrained (supply-side) +constrained_accessibility("total", travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_exponential(0.1), + demand = NULL, + supply = "jobs", + return_demand_side = FALSE +) + +# Singly-constrained (demand-side) +constrained_accessibility("singly", travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_exponential(0.1), + demand = "population", + supply = "jobs", + return_demand_side = TRUE +) + +# Doubly-constrained: use a small toy dataset with matching totals +tm_small <- data.table::data.table( + expand.grid(from_id = c("1","2","3"), to_id = c("1","2","3")) +) +tm_small[, travel_time := c(10, 30, 15, 30, 10, 25, 15, 25, 10)] +lu_small <- data.table::data.table( + id = c("1","2","3"), + population = c(4, 10, 6), # sum = 20 + jobs = c(7, 5, 8) # sum = 20 +) + +constrained_accessibility("doubly", tm_small, lu_small, + travel_cost = "travel_time", + decay_function = decay_exponential(0.1), + demand = "population", + supply = "jobs", + return_demand_side = NULL +) +} +\concept{Constrained accessibility} diff --git a/man/doubly_constrained.Rd b/man/doubly_constrained.Rd new file mode 100644 index 0000000..ed8df2d --- /dev/null +++ b/man/doubly_constrained.Rd @@ -0,0 +1,19 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/doubly_constrained.R +\name{doubly_constrained} +\alias{doubly_constrained} +\title{Doubly constrained accessibility} +\value{ +A \code{data.table}/\code{data.frame} with either OD-level flows (\code{detailed_results = TRUE}) +or marginals; see Details. +} +\description{ +Calculates accessibility using Wilson's doubly-constrained gravity model. +This measure allocates flows between origins and destinations such that +origin totals equal demand and destination totals equal supply. +Internal helper used by \code{\link[=constrained_accessibility]{constrained_accessibility()}} when \code{constraint = "doubly"}. +} +\examples{ +NULL +} +\keyword{internal} diff --git a/man/floating_catchment_area.Rd b/man/floating_catchment_area.Rd index 29b1644..718faae 100644 --- a/man/floating_catchment_area.Rd +++ b/man/floating_catchment_area.Rd @@ -34,8 +34,8 @@ calculating accessibility levels.} with the travel cost between origins and destinations.} \item{demand}{A string. The name of the column in \code{land_use_data} with the -number of people in each origin that will be considered potential -competitors.} +number of opportunity-demanders at each origin (e.g., people) that will be +considered.} \item{method}{A string. Which floating catchment area measure to use. Current available options are \code{"2sfca"} and \code{"bfca"}. More info in the diff --git a/man/roxygen/templates/demand.R b/man/roxygen/templates/demand.R index cf22cc3..5f23260 100644 --- a/man/roxygen/templates/demand.R +++ b/man/roxygen/templates/demand.R @@ -1,3 +1,3 @@ #' @param demand A string. The name of the column in `land_use_data` with the -#' number of people in each origin that will be considered potential -#' competitors. +#' number of opportunity-demanders at each origin (e.g., people) that will be +#' considered. diff --git a/man/roxygen/templates/supply.R b/man/roxygen/templates/supply.R new file mode 100644 index 0000000..0dc752b --- /dev/null +++ b/man/roxygen/templates/supply.R @@ -0,0 +1,3 @@ +#' @param supply A string. The name of the column in `land_use_data` with the +#' number of opportunity supply at each destination (e.g., jobs, school-seats) +#' that will be considered. diff --git a/man/singly_constrained.Rd b/man/singly_constrained.Rd new file mode 100644 index 0000000..ad0cebf --- /dev/null +++ b/man/singly_constrained.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/singly_constrained.R +\name{singly_constrained} +\alias{singly_constrained} +\title{Singly constrained accessibility} +\value{ +A \code{data.table}/\code{data.frame} with results (structure mirrors the wrapper). +} +\description{ +Allocates opportunities at each destination proportionally based on travel impedance and population at the origin. Uses the logic of Wilon's single constraint. Returns values as either 'demand' or 'supply'. +Internal helper used by \code{\link[=constrained_accessibility]{constrained_accessibility()}} when \code{constraint = "singly"}. +} +\examples{ +NULL +} +\keyword{internal} diff --git a/man/spatial_availability.Rd b/man/spatial_availability.Rd index 9f744d3..d44d6be 100644 --- a/man/spatial_availability.Rd +++ b/man/spatial_availability.Rd @@ -35,8 +35,8 @@ calculating accessibility levels.} with the travel cost between origins and destinations.} \item{demand}{A string. The name of the column in \code{land_use_data} with the -number of people in each origin that will be considered potential -competitors.} +number of opportunity-demanders at each origin (e.g., people) that will be +considered.} \item{decay_function}{A \code{fuction} that converts travel cost into an impedance factor used to weight opportunities. This function should take a diff --git a/man/total_constrained.Rd b/man/total_constrained.Rd index a033e6f..a793997 100644 --- a/man/total_constrained.Rd +++ b/man/total_constrained.Rd @@ -2,122 +2,16 @@ % Please edit documentation in R/total_constrained.R \name{total_constrained} \alias{total_constrained} -\title{Total constrained} -\usage{ -total_constrained( - travel_matrix, - land_use_data, - opportunity, - travel_cost, - decay_function, - group_by = character(0), - fill_missing_ids = TRUE, - detailed_results = FALSE -) -} -\arguments{ -\item{travel_matrix}{A data frame. The travel matrix describing the costs -(i.e. travel time, distance, monetary cost, etc.) between the origins and -destinations in the study area. Must contain the columns \code{from_id}, \code{to_id} -and any others specified in \code{travel_cost}.} - -\item{land_use_data}{A data frame. The distribution of opportunities within -the study area cells. Must contain the columns \code{id} and any others -specified in \code{opportunity}.} - -\item{opportunity}{A string. The name of the column in \code{land_use_data} -with the number of opportunities/resources/services to be considered when -calculating accessibility levels.} - -\item{travel_cost}{A string. The name of the column in \code{travel_matrix} -with the travel cost between origins and destinations.} - -\item{decay_function}{A \code{fuction} that converts travel cost into an -impedance factor used to weight opportunities. This function should take a -\code{numeric} vector and also return a \code{numeric} vector as output, with the -same length as the input. For convenience, the package currently includes -the following functions: \code{\link[=decay_binary]{decay_binary()}}, \code{\link[=decay_exponential]{decay_exponential()}}, -\code{\link[=decay_power]{decay_power()}} and \code{\link[=decay_stepped]{decay_stepped()}}. See the documentation of each decay -function for more details.} - -\item{group_by}{A \code{character} vector. When not \code{character(0)} (the default), -indicates the \code{travel_matrix} columns that should be used to group the -accessibility estimates by. For example, if \code{travel_matrix} includes a -departure time column, that specifies the departure time of each entry in -the data frame, passing \code{"departure_time"} to this parameter results in -accessibility estimates grouped by origin and by departure time.} - -\item{fill_missing_ids}{A \code{logical}. When calculating grouped accessibility -estimates (i.e. when \code{by_col} is not \code{NULL}), some combinations of groups -and origins may be missing. For example, if a single trip can depart from -origin \code{A} at 7:15am and reach destination \code{B} within 55 minutes, but no -trips departing from \code{A} at 7:30am can be completed at all, this second -combination will not be included in the output. When \code{TRUE} (the default), -the function identifies which combinations would be left out and fills -their respective accessibility values with 0, which incurs in a -performance penalty.} - -\item{detailed_results}{A \code{logical}. Whether to return total constrained -results aggregated by origin-destination pair (\code{TRUE}) or by origin -(\code{FALSE}, the default). When \code{TRUE}, the output also includes the balancing -factor \code{kappa_tot} used in the calculation and the unconstrained -accessibility value \code{weighted_opportunity}. Please note that the argument -\code{fill_missing_ids} does not affect the output when \code{detailed_results} is -\code{TRUE}.} -} +\title{Total constrained accessibility} \value{ -A data frame containing the accessibility estimates for each -origin/destination (depending if \code{active} is \code{TRUE} or \code{FALSE}) in the -travel matrix. +A \code{data.table}/\code{data.frame} with results (structure mirrors the wrapper). } \description{ -Calculates total constrained accessibility, part of a family of accessibility -measures proposed in -\insertCite{soukhovFamilyAccessibilityMeasures2025;textual}{accessibility} that -proportionally allocates the total sum of opportunities in the region based -on the relative travel impedance between zones. It is linearly proportion to -the conventional gravity-based Hansen-type accessibility measure and units of -the resulting values can be understood as the number of opportunities in the -region that can be potentially interacted with. - -This function is generic over any kind of numeric travel cost, -such as distance, time and money. +Allocates total opportunities in the region proportionally based on travel impedance. Uses the logic of a total (or unconstrained by Wilon's terms) constraint. Returns values as either 'demand' or 'supply' +Internal helper used by \code{\link[=constrained_accessibility]{constrained_accessibility()}} when \code{constraint = "total"}. } \examples{ -\dontshow{if (identical(tolower(Sys.getenv("NOT_CRAN")), "true")) withAutoprint(\{ # examplesIf} -# the example below is based on Soukhov et al. (2023) paper +NULL -travel_matrix <- data.table::data.table( - from_id = rep(c("A", "B", "C"), each = 3), - to_id = as.character(rep(1:3, 3)), - travel_time = c(15, 30, 100, 30, 15, 100, 100, 100, 15) -) -land_use_data <- data.table::data.table( - id = c("A", "B", "C", "1", "2", "3"), - population = c(50000, 150000, 10000, 0, 0, 0), - jobs = c(0, 0, 0, 100000, 100000, 10000) -) - -df <- total_constrained( - travel_matrix, - land_use_data, - opportunity = "jobs", - travel_cost = "travel_time", - decay_function = decay_exponential(decay_value = 0.1) -) -df - -detailed_df <- total_constrained( - travel_matrix, - land_use_data, - opportunity = "jobs", - travel_cost = "travel_time", - decay_function = decay_exponential(decay_value = 0.1), - detailed_results = TRUE -) -detailed_df -\dontshow{\}) # examplesIf} -} -\references{ -\insertAllCited{} } +\keyword{internal} diff --git a/testing_constrained.Rmd b/testing_constrained.Rmd new file mode 100644 index 0000000..a66b389 --- /dev/null +++ b/testing_constrained.Rmd @@ -0,0 +1,337 @@ +--- +title: "testing_constrained" +output: html_document +--- + +```{r} +library(data.table) +library(dplyr) +source("R/assert_and_assign.R") +source("R/merge_by_reference.R") +source("R/apply_gravity_measure.R") +source("R/warn_extra_cols.R") +source("R/fill_missing_ids.R") + + +# source("R/total_constrained.R") +# source("R/singly_constrained.R") +#source("R/doubly_constrained.R") +source("R/constrained_accessibility.R") + +source("R/decay_exponential.R") +source("R/decay_power.R") +source("R/spatial_availability.R") +source("R/floating_catchment_area.R") +``` + +# Test 1. Toy data from the spaital availability (2023) paper. +## TOT CONST +```{r} +# Toy travel matrix and land use data +travel_matrix <- expand.grid(from_id = c("1", "2", "3"), + to_id = c("1", "2", "3")) |> + dplyr::mutate(travel_time = c(10, 30, 15, 30, 10, 25, 15, 25, 10)) |> + data.table::data.table() + +land_use_data <- data.table::data.table( + id = c("1", "2", "3"), + population = c(4, 10, 6), # supply + jobs = c(160, 150, 180) # demand +) + +# Decay function +decay <- decay_exponential(decay_value = 0.1) +``` + + +```{r} +agg_false <- constrained_accessibility("total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay, + demand = NULL, + supply = "jobs", + detailed_results = FALSE, + return_demand_side = FALSE) +agg_false + +agg_true <- constrained_accessibility("total", travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay, + demand = "population", + supply = NULL, + detailed_results = FALSE, + return_demand_side = TRUE) +agg_true + +det_false <- constrained_accessibility("total", travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay, + demand = NULL, + supply = "jobs", + detailed_results = TRUE, + return_demand_side = FALSE) +det_false +det_true <- constrained_accessibility("total", travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay, + demand = "population", + supply = NULL, + detailed_results = TRUE, + return_demand_side = TRUE) +det_true +``` + +```{r} +det_true$demand |> sum() +det_false$supply |> sum() +agg_true$demand |> sum() +agg_false$supply |> sum() + +det_false$kappa_total |> sum() +det_true$hatkappa_total |> sum() +``` + +## SINGLE CONST +```{r} +agg_false <- constrained_accessibility("singly", travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay, + demand = "population", + supply = "jobs", + detailed_results = FALSE, + return_demand_side = FALSE) +agg_false + +agg_true <- constrained_accessibility("singly", travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay, + demand = "population", + supply = "jobs", + detailed_results = FALSE, + return_demand_side = TRUE) +agg_true + +det_false <- constrained_accessibility("singly", travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay, + demand = "population", + supply = "jobs", + detailed_results = TRUE, + return_demand_side = FALSE) +det_false + + +det_true <- constrained_accessibility("singly", travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay, + demand = "population", + supply = "jobs", + detailed_results = TRUE, + return_demand_side = TRUE) +det_true +``` + +Checks! +```{r} +det_true$demand |> sum() #accessible population (490) +det_false$supply |> sum() #accessible opps (20) +agg_true$demand |> sum() #accessible population (490) +agg_false$supply |> sum() #accessible opps (20) +``` + + +```{r} +det_false$kappa_singly |> sum() #should be the number of destinations +det_false[, sum(kappa_singly), by = to_id] #should all be 1 + +det_true$hatkappa_singly |> sum() #should be the number of origins +det_true[, sum(hatkappa_singly), by = from_id] #should all be 1 +``` + +# Test 2. Testing results with total/singly constrained toy example in 2025 paper +Data: +```{r test-data} +travel_matrix <- expand.grid(from_id = c("1", "2", "3"), + to_id = c("1", "2", "3")) |> + mutate(travel_time = c(10, 30, 15, 30, 10, 25, 15, 25, 10)) |> data.table() + +land_use_data <- data.table( + id = c("1", "2", "3"), + population = c(4, 10, 6), + jobs = c(160, 150, 180) +) +``` + +11/12 --> total_constrained works! I tested it for all decay scenarios (3, 2, 0.1). +```{r} +constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 3), + demand = NULL, + supply = "jobs", + detailed_results = TRUE, + return_demand_side = FALSE) + +constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 3), + demand = NULL, + supply = "jobs", + detailed_results = FALSE, + return_demand_side = FALSE) + + +constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 3), + demand = "population", + supply = NULL, + detailed_results = TRUE, + return_demand_side = TRUE) + +constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 3), + demand = "population", + supply = NULL, + detailed_results = FALSE, + return_demand_side = TRUE) +``` + +11/12 --> Singly-constrained works! I tested it for all decay scenarios (3, 2, 0.1)! +```{r} +constrained_accessibility("singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 3), + demand = "population", + supply = "jobs", + detailed_results = TRUE, + return_demand_side = FALSE) + +constrained_accessibility("singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 3), + demand = "population", + supply = "jobs", + detailed_results = FALSE, + return_demand_side = FALSE) +# +``` + +# Test 3. Testing results with doubly constrained toy example in 2025 paper + +```{r} +# Toy dataset where totals match (sum demand = sum supply = 20) +travel_matrix <- expand.grid(from_id = c("1", "2", "3"), + to_id = c("1", "2", "3")) |> + dplyr::mutate(travel_time = c(10, 30, 15, 30, 10, 25, 15, 25, 10)) |> + data.table() + +land_use_data <- data.table( + id = c("1", "2", "3"), + population = c(4, 10, 6), + jobs = c(7, 5, 8) +) +``` +#-- +11/12 --> tested this for decay_value =3 (table 10), and check zone 2 for all scenarios (3, 2, 0.1) in Table 11 and they match! Though.. the paper accidentally flips the columns. +```{r} +det <- + constrained_accessibility("doubly", + #doubly_constrained( + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = "population", + supply = "jobs", + detailed_results = TRUE, + return_demand_side = NULL) + +det + +agg <- + constrained_accessibility("doubly", + #doubly_constrained( + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = "population", + supply = "jobs", + detailed_results = FALSE, + return_demand_side = NULL) +``` + + +```{r} +agg$flow |> sum() #should be double, +det$flow |> sum() + +det$kappa_doubly |> sum() #should be the number of destinations or origins +``` + +# Test 4. Testing singly constrained's equality to spatial_availbility() and 2SFCA(). +```{r} +travel_matrix <- expand.grid(from_id = c("1", "2", "3"), + to_id = c("1", "2", "3")) |> + mutate(travel_time = c(10, 30, 15, 30, 10, 25, 15, 25, 10)) |> data.table() + +land_use_data <- data.table( + id = c("1", "2", "3"), + population = c(4, 10, 6), + jobs = c(160, 150, 180) +) +``` + +Spatial_availability() equality: +```{r} +constrained_accessibility("singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = "population", + supply = "jobs", + detailed_results = FALSE, + return_demand_side = FALSE) + +spatial_availability(travel_matrix, + land_use_data, + opportunity = "jobs", + travel_cost = "travel_time", + demand = "population", + decay_function = decay_power(decay_value = 2), + alpha = 1 + ) +``` + +floating_catchment_area(method="2sfca") equality: +```{r} +result_singly <- constrained_accessibility("singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = "population", + supply = "jobs", + detailed_results = FALSE, + return_demand_side = FALSE) + +round(result_singly$supply / land_use_data$population, 6) + +floating_catchment_area( + travel_matrix, + land_use_data, + opportunity = "jobs", + travel_cost = "travel_time", + demand = "population", + method = "2sfca", + decay_function = decay_power(decay_value = 2) + ) + +``` +yay! diff --git a/tests/testthat/test-constrained-essentials.R b/tests/testthat/test-constrained-essentials.R new file mode 100644 index 0000000..405c382 --- /dev/null +++ b/tests/testthat/test-constrained-essentials.R @@ -0,0 +1,151 @@ +# Small internal helpers +make_tm_3x3 <- function() { + tm <- data.table::data.table( + expand.grid(from_id = c("1","2","3"), to_id = c("1","2","3")) + ) + tm[, travel_time := c(10, 30, 15, 30, 10, 25, 15, 25, 10)] + tm +} +make_lu_large <- function() { + data.table::data.table( + id = c("1","2","3"), + population = c(4, 10, 6), + jobs = c(160, 150, 180) + ) +} +make_lu_match20 <- function() { + data.table::data.table( + id = c("1","2","3"), + population = c(4, 10, 6), # sum = 20 + jobs = c(7, 5, 8) # sum = 20 + ) +} + +exp_decay_01 <- decay_exponential(decay_value = 0.1) +pow_decay_3 <- decay_power(decay_value = 3) + +# ------------------------------- TOTAL ------------------------------- +test_that("total-constrained: preserves totals on supply/demand sides", { + tm <- make_tm_3x3() + lu <- make_lu_large() + + # Supply-side output (returns accessible opportunities) + det_supply <- constrained_accessibility("total", tm, lu, + travel_cost = "travel_time", decay_function = exp_decay_01, + demand = NULL, supply = "jobs", + detailed_results = TRUE, return_demand_side = FALSE + ) + agg_supply <- constrained_accessibility("total", tm, lu, + travel_cost = "travel_time", decay_function = exp_decay_01, + demand = NULL, supply = "jobs", + detailed_results = FALSE, return_demand_side = FALSE + ) + + # Demand-side output (returns accessible population) + det_demand <- constrained_accessibility("total", tm, lu, + travel_cost = "travel_time", decay_function = exp_decay_01, + demand = "population", supply = NULL, + detailed_results = TRUE, return_demand_side = TRUE + ) + agg_demand <- constrained_accessibility("total", tm, lu, + travel_cost = "travel_time", decay_function = exp_decay_01, + demand = "population", supply = NULL, + detailed_results = FALSE, return_demand_side = TRUE + ) + + expect_equal(sum(det_supply$supply), sum(lu$jobs)) + expect_equal(sum(agg_supply$supply), sum(lu$jobs)) + expect_equal(sum(det_demand$demand), sum(lu$population)) + expect_equal(sum(agg_demand$demand), sum(lu$population)) +}) + +# ------------------------------ SINGLY ------------------------------- +test_that("singly-constrained: preserves totals and kappa properties", { + tm <- make_tm_3x3() + lu <- make_lu_large() + + # Supply-constrained (returns accessible opportunities) + det_supply <- constrained_accessibility("singly", tm, lu, + travel_cost = "travel_time", decay_function = exp_decay_01, + demand = "population", supply = "jobs", + detailed_results = TRUE, return_demand_side = FALSE + ) + agg_supply <- constrained_accessibility("singly", tm, lu, + travel_cost = "travel_time", decay_function = exp_decay_01, + demand = "population", supply = "jobs", + detailed_results = FALSE, return_demand_side = FALSE + ) + + # Demand-constrained (returns accessible population) + det_demand <- constrained_accessibility("singly", tm, lu, + travel_cost = "travel_time", decay_function = exp_decay_01, + demand = "population", supply = "jobs", + detailed_results = TRUE, return_demand_side = TRUE + ) + agg_demand <- constrained_accessibility("singly", tm, lu, + travel_cost = "travel_time", decay_function = exp_decay_01, + demand = "population", supply = "jobs", + detailed_results = FALSE, return_demand_side = TRUE + ) + + # Totals + expect_equal(sum(det_supply$supply), sum(lu$jobs)) + expect_equal(sum(agg_supply$supply), sum(lu$jobs)) + expect_equal(sum(det_demand$demand), sum(lu$population)) + expect_equal(sum(agg_demand$demand), sum(lu$population)) + + # κ normalization properties (these DO hold for singly) + by_to <- det_supply[, .(sum_k = sum(kappa_singly)), by = to_id] + by_from <- det_demand[, .(sum_k = sum(hatkappa_singly)), by = from_id] + expect_true(all(abs(by_to$sum_k - 1) < 1e-8)) # per-destination sums == 1 + expect_true(all(abs(by_from$sum_k - 1) < 1e-8)) # per-origin sums == 1 + expect_equal(sum(det_supply$kappa_singly), length(unique(det_supply$to_id))) + expect_equal(sum(det_demand$hatkappa_singly), length(unique(det_demand$from_id))) +}) + +# ------------------------------ DOUBLY ------------------------------- +test_that("doubly-constrained: errors when totals mismatch", { + tm <- make_tm_3x3() + lu_bad <- data.table::copy(make_lu_match20()) + lu_bad[id == "3", jobs := jobs + 1] # break equality + + expect_error( + constrained_accessibility("doubly", tm, lu_bad, + travel_cost = "travel_time", decay_function = exp_decay_01, + demand = "population", supply = "jobs", + detailed_results = TRUE, return_demand_side = NULL + ), + "sum of origins must equal the sum of destinations" + ) +}) + + + +test_that("doubly-constrained: OD flows match marginals", { + tm <- make_tm_3x3() + lu <- make_lu_match20() + decay <- pow_decay_3 + + det <- constrained_accessibility("doubly", tm, lu, + travel_cost = "travel_time", decay_function = decay, + demand = "population", supply = "jobs", + detailed_results = TRUE, return_demand_side = NULL + ) + + # Robustly pick the numeric flow column + candidate_cols <- c("flow", "supply", "demand") + value_col <- intersect(candidate_cols, names(det))[1] + expect_true(!is.null(value_col), info = "No flow/supply/demand column in detailed output") + expect_true(is.numeric(det[[value_col]]), info = "Flow column must be numeric") + + # Sum flows by origin and destination + by_origin <- det[, .(sum_flow = sum(get(value_col))), by = from_id] + by_dest <- det[, .(sum_flow = sum(get(value_col))), by = to_id] + + # Build expected marginals from land use + O <- lu[match(by_origin$from_id, lu$id), population] + D <- lu[match(by_dest$to_id, lu$id), jobs] + + expect_equal(by_origin$sum_flow, O, tolerance = 1e-3) + expect_equal(by_dest$sum_flow, D, tolerance = 1e-3) +}) diff --git a/tests/testthat/test-constrained-wrapper.R b/tests/testthat/test-constrained-wrapper.R new file mode 100644 index 0000000..f47105b --- /dev/null +++ b/tests/testthat/test-constrained-wrapper.R @@ -0,0 +1,45 @@ + +testthat::local_edition(3) + +test_that("constrained_accessibility: wrapper argument rules are enforced", { + tm <- data.table::data.table( + expand.grid(from_id = c("1","2"), to_id = c("1","2")) + ) + tm[, travel_time := c(10, 20, 15, 10)] + lu <- data.table::data.table( + id = c("1","2"), + population = c(5, 7), + jobs = c(100, 120) + ) + decay <- decay_exponential(decay_value = 0.1) + + # invalid constraint -> any error is fine + expect_error( + constrained_accessibility("invalid", tm, lu, + travel_cost = "travel_time", decay_function = decay, + demand = "population", supply = "jobs", return_demand_side = TRUE + ) + ) + + # total/singly require TRUE or FALSE (not NULL) + expect_error( + constrained_accessibility("total", tm, lu, + travel_cost = "travel_time", decay_function = decay, + demand = "population", supply = "jobs", return_demand_side = NULL + ) + ) + expect_error( + constrained_accessibility("singly", tm, lu, + travel_cost = "travel_time", decay_function = decay, + demand = "population", supply = "jobs", return_demand_side = NULL + ) + ) + + # doubly requires NULL (error if TRUE) + expect_error( + constrained_accessibility("doubly", tm, lu, + travel_cost = "travel_time", decay_function = decay, + demand = "population", supply = "jobs", return_demand_side = TRUE + ) + ) +}) diff --git a/tests/testthat/test-total_constrained.R b/tests/testthat/test-total_constrained.R deleted file mode 100644 index c1082a7..0000000 --- a/tests/testthat/test-total_constrained.R +++ /dev/null @@ -1,276 +0,0 @@ -# if running manually, please run the following line first: -# source("tests/testthat/setup.R") - -tester <- function( - travel_matrix = smaller_matrix, - land_use_data = get("land_use_data", envir = parent.frame()), - opportunity = "jobs", - travel_cost = "travel_time", - decay_function = decay_exponential(0.1), - group_by = "mode", - fill_missing_ids = TRUE, - detailed_results = FALSE -) { - total_constrained( - travel_matrix, - land_use_data, - opportunity, - travel_cost, - decay_function, - group_by, - fill_missing_ids, - detailed_results - ) -} - -test_that("raises errors due to incorrect input", { - expect_error(tester(decay_function = "a")) - expect_error(tester(decay_function = mean)) - expect_error(tester(decay_function = get)) - - expect_error(tester(opportunity = 1)) - expect_error(tester(opportunity = c("schools", "jobs"))) - - expect_error(tester(travel_cost = 1)) - expect_error(tester(travel_cost = c("travel_time", "monetary_cost"))) - - expect_error(tester(group_by = 1)) - expect_error(tester(group_by = NA)) - expect_error(tester(group_by = "from_id")) - expect_error(tester(group_by = c("mode", "mode"))) - - expect_error(tester(fill_missing_ids = 1)) - expect_error(tester(fill_missing_ids = c(TRUE, TRUE))) - expect_error(tester(fill_missing_ids = NA)) - - expect_error(tester(detailed_results = 1)) - expect_error(tester(detailed_results = c(TRUE, TRUE))) - expect_error(tester(detailed_results = NA)) - - expect_error(tester(as.list(travel_matrix))) - expect_error(tester(travel_matrix[, .(oi = from_id, to_id, travel_time)])) - expect_error(tester(travel_matrix[, .(from_id, oi = to_id, travel_time)])) - expect_error( - tester( - travel_matrix[, .(from_id, to_id, oi = travel_time)], - travel_cost = "travel_time" - ) - ) - expect_error( - tester( - travel_matrix[, .(from_id, to_id, travel_time, oi = mode)], - group_by = "mode" - ) - ) - - expect_error(tester(as.list(land_use_data))) - expect_error( - tester(land_use_data = land_use_data[, .(oi = id, jobs)]) - ) - expect_error( - tester( - land_use_data = land_use_data[, .(id, oi = jobs)], - opportunity = "jobs" - ) - ) -}) - -test_that("throws warning if travel_matrix has extra col", { - # i.e. col not listed in travel_cost and by_col - expect_warning(tester(group_by = character(0))) -}) - -test_that("returns a dataframe whose class is the same as travel_matrix's", { - result <- tester() - expect_is(result, "data.table") - result <- tester(land_use_data = as.data.frame(land_use_data)) - expect_is(result, "data.table") - - result <- tester(as.data.frame(smaller_matrix)) - expect_false(inherits(result, "data.table")) - expect_is(result, "data.frame") - result <- tester(as.data.frame(smaller_matrix), as.data.frame(land_use_data)) - expect_false(inherits(result, "data.table")) - expect_is(result, "data.frame") -}) - -test_that("result has correct structure", { - result <- tester() - expect_true(ncol(result) == 3) - expect_is(result$id, "character") - expect_is(result$mode, "character") - expect_is(result$jobs, "numeric") - - result <- tester(opportunity = "schools") - expect_true(ncol(result) == 3) - expect_is(result$id, "character") - expect_is(result$mode, "character") - expect_is(result$schools, "numeric") - - suppressWarnings(result <- tester(group_by = character(0))) - expect_true(ncol(result) == 2) - expect_is(result$id, "character") - expect_is(result$jobs, "numeric") - - result <- tester( - data.table::data.table( - mode = character(), - from_id = character(), - to_id = character(), - travel_time = integer() - ) - ) - expect_true(ncol(result) == 3) - expect_true(nrow(result) == 0) - expect_is(result$id, "character") - expect_is(result$mode, "character") - expect_is(result$jobs, "numeric") -}) - -test_that("input data sets remain unchanged", { - original_smaller_matrix <- travel_matrix[1:10] - original_land_use_data <- readRDS(file.path(data_dir, "land_use_data.rds")) - - result <- tester() - - # Reset all indices for comparison - data.table::setindex(smaller_matrix, NULL) - data.table::setindex(land_use_data, NULL) - data.table::setindex(original_smaller_matrix, NULL) - data.table::setindex(original_land_use_data, NULL) - - expect_equal(original_smaller_matrix, smaller_matrix) - expect_equal(original_land_use_data, land_use_data) -}) - -test_that("fill_missing_ids arg works correctly", { - small_travel_matrix <- travel_matrix[ - from_id %in% c("89a88cdb57bffff", "89a88cdb597ffff") - ] - small_travel_matrix <- small_travel_matrix[ - !(from_id == "89a88cdb57bffff" & mode == "transit2") - ] - - result <- tester(small_travel_matrix, fill_missing_ids = TRUE) - result[, jobs := as.integer(jobs)] - data.table::setkey(result, NULL) - expect_identical( - result, - data.table::data.table( - id = rep(c("89a88cdb57bffff", "89a88cdb597ffff"), each = 2), - mode = rep(c("transit", "transit2"), times = 2), - jobs = c(164981L, 0L, 165553L, 165553L) # Update these expected values - ) - ) - - result <- tester(small_travel_matrix, fill_missing_ids = FALSE) - result[, jobs := as.integer(jobs)] - expect_identical( - result, - data.table::data.table( - id = c("89a88cdb57bffff", "89a88cdb597ffff", "89a88cdb597ffff"), - mode = c("transit", "transit", "transit2"), - jobs = c(164981L, 165553L, 165553L) # Update these expected values - ) - ) -}) - -test_that("accepts custom decay function", { - selected_ids <- c( - "89a88cdb57bffff", - "89a88cdb597ffff", - "89a88cdb5b3ffff", - "89a88cdb5cfffff", - "89a88cd909bffff" - ) - smaller_travel_matrix <- travel_matrix[ - from_id %in% selected_ids & to_id %in% selected_ids - ] - - custom_function <- function(travel_cost) rep(1L, length(travel_cost)) - - result <- tester(smaller_travel_matrix, decay_function = custom_function) - result[, jobs := round(jobs, digits = 1)] #for comparision purposes - expect_identical( - result, - data.table::data.table( - id = rep(selected_ids, 2), - mode = rep(c("transit", "transit2"), each = 5), - jobs = rep(49608.8, 10) # All values are the same - ) - ) -}) - -test_that("calculates total constrained accessibility correctly", { - # Use the same test data from spatial_availability tests - paper_travel_matrix <- data.table::data.table( - from_id = rep(c("A", "B", "C"), each = 3), - to_id = as.character(rep(1:3, 3)), - travel_time = c(15, 30, 100, 30, 15, 100, 100, 100, 15) - ) - paper_land_use_data <- data.table::data.table( - id = c("A", "B", "C", "1", "2", "3"), - jobs = c(0, 0, 0, 100000, 100000, 10000) - ) - - result <- tester( - paper_travel_matrix, - paper_land_use_data, - group_by = character(0) - ) - #result[, jobs := round(jobs, digits = 0)] - - # Test conservation property - sum should equal total opportunities - total_jobs <- sum(paper_land_use_data[id %in% c("1", "2", "3"), jobs]) - expect_equal(sum(result$jobs), total_jobs)#, tolerance = 1) - - # Test that all origins are present - expect_setequal(result$id, c("A", "B", "C")) -}) - -test_that("results are grouped by decay_function_arg when needed", { - small_travel_matrix <- travel_matrix[ - from_id %in% c("89a88cdb57bffff", "89a88cdb597ffff") & - mode != "transit2" - ] - - result <- tester( - small_travel_matrix, - decay_function = decay_exponential(c(0.5, 0.6)) - ) - #result[, jobs := round(jobs, 1)] - - expect_true("decay_function_arg" %in% names(result)) - expect_equal(unique(result$decay_function_arg), c(0.5, 0.6)) -}) - -test_that("throws warning w/ fill_missing_ids = FALSE with detailed_results", { - expect_warning(tester(fill_missing_ids = FALSE, detailed_results = TRUE)) -}) - -test_that("result has correct structure with detailed_results = TRUE", { - result <- tester(detailed_results = TRUE) - expect_true(ncol(result) >= 5) # from_id, to_id, group_cols, kappa_tot, jobs - expect_is(result$mode, "character") - expect_is(result$from_id, "character") - expect_is(result$to_id, "character") - expect_is(result$kappa_tot, "numeric") - expect_is(result$jobs, "numeric") - - result <- tester(detailed_results = TRUE, opportunity = "schools") - expect_true(ncol(result) >= 5) - expect_is(result$mode, "character") - expect_is(result$from_id, "character") - expect_is(result$to_id, "character") - expect_is(result$kappa_tot, "numeric") - expect_is(result$schools, "numeric") - - result <- tester(smaller_matrix[0], detailed_results = TRUE) - expect_true(ncol(result) >= 5) - expect_true(nrow(result) == 0) - expect_is(result$mode, "character") - expect_is(result$from_id, "character") - expect_is(result$to_id, "character") - expect_is(result$kappa_tot, "numeric") - expect_is(result$jobs, "numeric") -}) From e7c07bfb86b100c0ff963468d61f124ae944e794 Mon Sep 17 00:00:00 2001 From: Rafael H M Pereira Date: Fri, 5 Dec 2025 21:38:30 -0300 Subject: [PATCH 4/4] Clean up references by removing abstracts and file paths Removed abstracts and file paths from references. --- inst/REFERENCES.bib | 30 ++++++++---------------------- 1 file changed, 8 insertions(+), 22 deletions(-) diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 2f32593..5ea0171 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -9,9 +9,7 @@ @article{barboza2021balancing pages = {102924}, issn = {09666923}, doi = {10.1016/j.jtrangeo.2020.102924}, - abstract = {One of the most common applications of accessibility is in evaluating inequality in access to jobs. A vital factor to be incorporated by accessibility indicators when analyzing job accessibility is the competition for job positions by job seekers; otherwise, the results may be inaccurate or misleading. Despite efforts by researchers to develop accessibility measures that capture job competition, they fail to ensure that these measures are practical and easily interpretable and communicable, which in turn makes planners and policymakers continue to opt for more straightforward measures. In this paper we aim to fill this gap by providing a simple accessibility measure that accounts for competition effects, while remaining practical, intuitive, and highly communicable. The proposed indicator -- Balancing Time -- is applied to assess the inequality in job accessibility in 160 neighborhoods within the city of Rio de Janeiro, and the results are compared with the most popular indicator used in practice, the cumulative opportunities. The findings suggest that Balancing Time overcomes some of the limitations of cumulative opportunities and that it is a useful tool for planners, particularly in the cities with job opportunities concentrated in central areas. Given its simplicity, Balancing Time is especially relevant in the context of the Global South, where most transport agencies face data limitations and have low skilled technical staff.}, - langid = {english}, - file = {/home/dhersz/Zotero/storage/6NEC3LRI/Barboza et al. - 2021 - Balancing time Using a new accessibility measure .pdf} + langid = {english} } @article{bauer2016measuring, @@ -41,10 +39,8 @@ @article{erreygers2009correcting pages = {504--515}, issn = {0167-6296}, doi = {10.1016/j.jhealeco.2008.02.003}, - abstract = {In recent years attention has been drawn to several shortcomings of the Concentration Index, a frequently used indicator of the socioeconomic inequality of health. Some modifications have been suggested, but these are only partial remedies. This paper proposes a corrected version of the Concentration Index which is superior to the original Concentration Index and its variants, in the sense that it is a rank-dependent indicator which satisfies four key requirements (transfer, level independence, cardinal invariance, and mirror). The paper also shows how the corrected Concentration Index can be decomposed and generalized.}, langid = {english}, - keywords = {Concentration Index,Health inequality,Socioeconomic inequality}, - file = {/home/dhersz/Zotero/storage/56WQ9IGG/S0167629608000076.html} + keywords = {Concentration Index,Health inequality,Socioeconomic inequality} } @article{foster1984class, @@ -75,8 +71,7 @@ @article{kwan1998spacetime issn = {00167363}, doi = {10.1111/j.1538-4632.1998.tb00396.x}, urldate = {2022-06-07}, - langid = {english}, - file = {/home/dhersz/Zotero/storage/KEXX4KKN/Kwan - 2010 - Space-Time and Integral Measures of Individual Acc.pdf;/home/dhersz/Zotero/storage/VX6AZ25Z/Kwan - 2010 - Space-Time and Integral Measures of Individual Acc.pdf} + langid = {english} } @article{luo2003measures, @@ -93,8 +88,7 @@ @article{luo2003measures publisher = {SAGE Publications Sage UK: London, England}, issn = {0265-8135, 1472-3417}, doi = {10.1068/b29120}, - langid = {english}, - file = {/home/dhersz/Zotero/storage/EIXZH5GT/Luo and Wang - 2003 - Measures of Spatial Accessibility to Health Care i.pdf} + langid = {english} } @article{paez2019demand, @@ -110,8 +104,7 @@ @article{paez2019demand issn = {1932-6203}, doi = {10.1371/journal.pone.0218773}, langid = {english}, - keywords = {Body weight,Geography,Health care facilities,Health care policy,Health care providers,Ontario,Physicians,Urban areas}, - file = {/home/dhersz/Zotero/storage/C9UJVHK6/Paez et al. - 2019 - Demand and level of service inflation in Floating .pdf;/home/dhersz/Zotero/storage/FD2GDJ8N/Paez et al. - 2019 - Demand and level of service inflation in Floating .pdf;/home/dhersz/Zotero/storage/ISWY2APP/article.html} + keywords = {Body weight,Geography,Health care facilities,Health care policy,Health care providers,Ontario,Physicians,Urban areas} } @article{pereira2021geographic, @@ -124,10 +117,8 @@ @article{pereira2021geographic pages = {113773}, issn = {0277-9536}, doi = {10.1016/j.socscimed.2021.113773}, - abstract = {The rapid spread of COVID-19 across the world has raised concerns about the responsiveness of cities and healthcare systems during pandemics. Recent studies try to model how the number of COVID-19 infections will likely grow and impact the demand for hospitalization services at national and regional levels. However, less attention has been paid to the geographic access to COVID-19 healthcare services and to hospitals' response capacity at the local level, particularly in urban areas in the Global South. This paper shows how transport accessibility analysis can provide actionable information to help improve healthcare coverage and responsiveness. It analyzes accessibility to COVID-19 healthcare at high spatial resolution in the 20 largest cities of Brazil. Using network-distance metrics, we estimate the vulnerable population living in areas with poor access to healthcare facilities that could either screen or hospitalize COVID-19 patients. We then use a new balanced floating catchment area (BFCA) indicator to estimate spatial, income, and racial inequalities in access to hospitals with intensive care unit (ICU) beds and mechanical ventilators while taking into account congestion effects. Based on this analysis, we identify substantial social and spatial inequalities in access to health services during the pandemic. The availability of ICU equipment varies considerably between cities, and it is substantially lower among black and poor communities. The study maps territorial inequalities in healthcare access and reflects on different policy lessons that can be learned for other countries based on the Brazilian case.}, langid = {english}, - keywords = {Accessibility,Brazil,COVID-19,Equity,Floating catchment area,ICU,Race,Ventilators}, - file = {/home/dhersz/Zotero/storage/5D6YJBBD/Pereira et al. - 2021 - Geographic access to COVID-19 healthcare in Brazil.pdf;/home/dhersz/Zotero/storage/8WZL4NCW/Pereira et al. - 2021 - Geographic access to COVID-19 healthcare in Brazil.pdf;/home/dhersz/Zotero/storage/PE6IRGRZ/S0277953621001052.html} + keywords = {Accessibility,Brazil,COVID-19,Equity,Floating catchment area,ICU,Race,Ventilators} } @article{soukhov2023introducing, @@ -142,9 +133,7 @@ @article{soukhov2023introducing pages = {e0278468}, issn = {1932-6203}, doi = {10.1371/journal.pone.0278468}, - abstract = {Accessibility indicators are widely used in transportation, urban and healthcare planning, among many other applications. These measures are weighted sums of reachable opportunities from a given origin, conditional on the cost of movement, and are estimates of the potential for spatial interaction. Over time, various proposals have been forwarded to improve their interpretability: one of those methodological additions have been the introduction of competition. In this paper we focus on competition, but first demonstrate how a widely used measure of accessibility with congestion fails to properly match the opportunity-seeking population. We then propose an alternative formulation of accessibility with competition, a measure we call spatial availability . This measure relies on proportional allocation balancing factors (friction of distance and population competition) that are equivalent to imposing a single constraint on conventional gravity-based accessibility. In other words, the proportional allocation of opportunities results in a spatially available opportunities value which is assigned to each origin that, when all origin values are summed, equals the total number of opportunities in the region. We also demonstrate how Two-Stage Floating Catchment Area (2SFCA) methods are equivalent to spatial availability and can be reconceptualized as singly-constrained accessibility. To illustrate the application of spatial availability and compare it to other relevant measures, we use data from the 2016 Transportation Tomorrow Survey of the Greater Golden Horseshoe area in southern Ontario, Canada. Spatial availability is an important contribution since it clarifies the interpretation of accessibility with competition and paves the way for future applications in equity analysis (e.g., spatial mismatch, opportunity benchmarking, policy intervention scenario analysis).}, - langid = {english}, - file = {/home/dhersz/Zotero/storage/P9M77AZM/Soukhov et al. - 2023 - Introducing spatial availability, a singly-constra.pdf} + langid = {english} } @article{tomasiello2023time, @@ -157,10 +146,8 @@ @article{tomasiello2023time pages = {103007}, issn = {0143-6228}, doi = {10.1016/j.apgeog.2023.103007}, - abstract = {Cumulative accessibility measures are the most-used accessibility metric and indicate the number of opportunities reached within a given travel time threshold. However, these measures require an ad-hoc choice of a single travel time threshold, which can influence the conclusions of transport project evaluations and equity analyses. In this paper, we introduce the time interval cumulative accessibility measure, a new metric that mitigates the impacts of arbitrary choices of trip duration on cumulative accessibility analyses while preserving computation and communicability advantages. The proposed indicator estimates the average (or median) number of opportunities that can be reached considering multiple minute-by-minute cutoffs within a given travel time interval. We demonstrate the new metric in a case study assessing how a planned subway expansion could impact employment accessibility in Fortaleza, Brazil. Using sensitivity analyses with Monte Carlo simulations, we show that the proposed metric makes the results of accessibility estimates and equity analyses significantly less sensitive to ad-hoc methodological choices while yielding results that are very similar to those obtained with traditional threshold-based measures. Future accessibility-oriented research and planning could benefit from the way in which the proposed time interval cumulative opportunity measure provides more robust accessibility estimates without compromising the communicability of results.}, langid = {english}, - keywords = {Cumulative accessibility,Equity,MTUP,Public transport,Time interval,Transport policy,Travel time}, - file = {/home/dhersz/Zotero/storage/NNWIAXLP/S0143622823001388.html} + keywords = {Cumulative accessibility,Equity,MTUP,Public transport,Time interval,Transport policy,Travel time} } @@ -168,7 +155,6 @@ @misc{soukhovFamilyAccessibilityMeasures2025 title = {A family of accessibility measures derived from spatial interaction principles}, url = {https://osf.io/a9dxb_v1}, doi = {10.31219/osf.io/a9dxb_v1}, - abstract = {Transportation planning has long prioritized the efficiency of movement, or mobility. However, the concept of accessibility represents a more comprehensive evolution, shifting focus from mere movement to the potential to reach (i.e., spatially interact) with desired destinations. Despite growing recognition of accessibility-based planning approaches, the concept remains fragmented, with inconsistent definitions and unclear interpretations. This work's aim is to clarify and unify the concept of accessibility by connecting it into spatial interaction modeling. We demonstrate that widely used mobility and accessibility models, such as gravity-based accessibility and spatial interaction models, share common theoretical roots. From this foundation, this paper offers three contributions: (A) we introduce a family of accessibility measures within the principles of spatial interaction, and (B) formally define four members of the family, namely the 'unconstrained' measure (i.e., Hansen-type accessibility), the 'total constrained' measure (i.e., a constrained version of the Hansen-type accessibility), the 'singly constrained' measure (i.e., related to the popular two step floating catchment approach - 2SFCA), and the 'doubly constrained' measure representing realized interactions or 'access', effectively equal to the doubly constrained spatial interaction model; and (C) we demonstrate the interpretability advantages of the family, as these constrained accessibility measures yield values in units of the number of potential "opportunities for spatial interaction" or "population for spatial interaction" for each zone and zonal flow. The family of accessibility measures proposed here clarifies the concept of 'potential' in accessibility, demonstrates theoretical and formulaic linkages across popular accessibility and spatial interaction models, and reintroduces measurement units into accessibility measures. By doing so, we believe this family of measures can unlock a clearer, more interpretable, and cohesive foundation for accessibility analysis.}, publisher = {{OSF}}, author = {Soukhov, Anastasia and Pereira, Rafael and Higgins, Christopher and Paez, Antonio}, urldate = {2025-09-19},