From 45444a1e6510f03909bad1336bebdb05f9a1a5d4 Mon Sep 17 00:00:00 2001 From: soukhova Date: Sat, 20 Dec 2025 21:43:00 -0500 Subject: [PATCH] feat: flip active logic, fix marginal sums for total constrained, and improved documentation --- R/constrained_accessibility.R | 295 +++++++--- R/doubly_constrained.R | 8 +- R/singly_constrained.R | 20 +- R/total_constrained.R | 55 +- man/constrained_accessibility.Rd | 264 +++++++-- man/doubly_constrained.Rd | 4 +- man/singly_constrained.Rd | 6 +- man/total_constrained.Rd | 10 +- tests/test_rafa/testing_constrained.Rmd | 543 ++++++++++++++++--- tests/testthat/test-constrained-essentials.R | 20 +- tests/testthat/test-constrained-wrapper.R | 4 +- 11 files changed, 986 insertions(+), 243 deletions(-) diff --git a/R/constrained_accessibility.R b/R/constrained_accessibility.R index 0e36556..fe544d7 100644 --- a/R/constrained_accessibility.R +++ b/R/constrained_accessibility.R @@ -1,71 +1,213 @@ #' Constrained accessibility #' -#' Calculates accessibility using constrained gravity models, as proposed in -#' \insertCite{soukhov2025family;textual}{accessibility}: -#' - `"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. -#' Please see the Details section for more information. +#' Calculates accessibility using constraints, as proposed in +#' \insertCite{soukhov2025family;textual}{accessibility}. Accessibility is +#' conceptualised as potential' spatial interaction. The results are in units of +#' opportunities in the region (i.e., active accessibility `active = TRUE`) or +#' in units of population in the region (if `active = FALSE` reflecting passive +#' accessibility). Results are presented as a sum of a zone or presented as an +#' origin-destination flow (if `detailed_results = TRUE`). #' -#' @template description_generic_cost +#' The following three constraint cases (along with their `active` and passive +#' variants) are included: +#' +#' - `"total"`: Allocates the system-wide total proportionally based on travel +#' impedance. +#' - If `active = TRUE`, results are in units of **opportunities** (supply) +#' accessible by the origin. Only `supply` can be passed, `demand` must be NULL. +#' - If `active = FALSE`, results are in units of **population** (demand) +#' accessible by the destination. Only `demand` can be passed, `supply` must +#' be NULL. +#' - If `detailed_results = TRUE`, accessibility (in units corresponding to the +#' logical for `active`) is presented as a flow along with intermediates. +#' +#' - `"singly"`: Applies a single constraint, allocating one side of the marginal +#' proportionally based on the other marginal and travel impedance: +#' - If `active = TRUE`, returns origin-side results (how much supply is +#' accessible from each origin). Outputs are units of supply. `demand` and +#' `supply` must be passed. +#' - If `active = FALSE`, returns destination-side results (how much demand is +#' accessible from each destination). Outputs are units of demand. `demand` and +#' `supply` must be passed. +#' - If `detailed_results = TRUE`, accessibility (in units corresponding to the +#' logical for `active`) is presented as a flow along with intermediates. +#' +#' - `"doubly"`: Allocates flows so supply at each destination matches demand at +#' each origin. +#' - OD flows are calibrated to both marginals. OD flows are calibrated to both +#' marginals using iterative proportional fitting. The sum of `demand` and +#' `supply` must match; otherwise, the function will not converge. +#' - `active` and `detailed_results` must be NULL. Since supply must match +#' demand, their units are the same and there is no distinction between 'active' +#' and 'passive' notions. +#' +#' Please see the Details section for more information. #' -#' @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"`. See #' Details section for more information. #' @param active A logical. When `TRUE`, the function calculates active #' accessibility (the quantity of opportunities that can be reached from -#' a given origin). when `FALSE`, it calculates passive accessibility (by +#' a given origin). When `FALSE`, it calculates passive accessibility (by #' how many people each destination can be reached), which is equivalent #' to the notion of market potential. This parameter only works for #' `constraint` types `"total"` and `"singly"`. Ignored for #' `constraint = "doubly"`. #' @param error_threshold Numeric. Convergence criterion used only for -#' doubly-constrained case. +#' calibration in the doubly-constrained case (`constraint = "doubly"`). #' @param improvement_threshold Numeric. Convergence criterion for improvement -#' used only for doubly-constrained case. -#' @param max_iterations Integer. Maximum iterations used only for -#' doubly-constrained calibration. +#' used only for calibration in the doubly-constrained case +#' (`constraint = "doubly"`). +#' @param max_iterations Integer. Maximum iterations used only for calibration +#' in the doubly-constrained case (`constraint = "doubly"`). #' @template group_by #' @template fill_missing_ids_combinations #' @param detailed_results Logical. Whether to return detailed OD-level results. #' +#' @template description_generic_cost +#' @template travel_matrix +#' @template land_use_data +#' @template travel_cost +#' @template decay_function +#' @template demand +#' @template supply +#' #' @section Details: #' This function covers the family of constrained accessibility measures #' proposed in \insertCite{soukhov2025family;textual}{accessibility}. #' -#' ## Total constrained accessibility +#' ## Total Constrained Accessibility +#' +#' Allocates the total system-wide quantity proportionally based on travel +#' impedance between origins and destinations. This measure uses the logic of a +#' total ~(or 'unconstrained' by Wilson's terms)~ constraint. +#' +#' Use this measure when the total quantity of **supply** OR **demand** in the +#' system is known and representing accessibility as a proportion of this total +#' is meaningful. +#' +#' **Requirements**: +#' - Either `demand` or `supply` must be provided (cannot provide both). +#' +#' **Interpretation**: +#' - `active = TRUE` (*active accessibility*): Results represent the total +#' number of **opportunities** (supply) accessible from each origin based on +#' region-relative travel impedance. The units are in 'supply' (e.g., jobs, +#' school seats). +#' - If `detailed_results = FALSE`, outputs are aggregated and returned by +#' origin. +#' - If `detailed_results = TRUE`, OD-level flows are returned. Summing flows +#' by origin equals the aggregated result. +#' +#' - `active = FALSE` (*passive accessibility*, the notion of market potential): +#' Results represent the total number of **population** (demand) that can reach +#' each destination based on region-relative travel impedance. The units are in +#' 'demand' (e.g., population). +#' - If `detailed_results = FALSE`, outputs are aggregated by destination. +#' - If `detailed_results = TRUE`, OD-level flows are returned. Summing flows +#' by destination equals the aggregated result. +#' +#' **Use cases**: +#' - Active accessibility (aggregated): +#' "How many jobs can be reached from origin zone A given its region-relative +#' travel impedance?" +#' +#' - Active accessibility (flow-level): +#' "How many jobs can be reached by flow A->1 given A->1's region-relative +#' travel impedance?" +#' +#' - Passive accessibility (aggregated): +#' "How many people can reach destination zone 1 given its region-relative +#' travel impedance?" +#' +#' - Passive accessibility (flow-level): +#' "How many people are reached by flow 1->A given 1->A's region-relative +#' travel impedance?" +#' +#' +#' ## Singly Constrained Accessibility +#' +#' Allocates opportunities at each destination (or population at each origin) +#' proportionally based on travel impedance and the opposite marginal. This +#' measure uses the logic of single constraint from +#' \insertCite{wilson1971family;textual}{accessibility}. +#' +#' Use this measure when modeling **competition**, where both demand and supply +#' are conceptualised to influence accessibility but only one side is fixed. +#' The measure distributes flows so that totals match the constrained side +#' while weighting by travel impedance and the unconstrained side. +#' +#' **Requirements**: +#' - Both `demand` and `supply` must be provided (the logical for `active` +#' determines if either demand or supply is constrained). +#' +#' **Interpretation**: +#' - `active = TRUE` (*active accessibility*): constrains supply. Results +#' represent the total number of **opportunities** (supply) accessible from each +#' origin based on region-relative travel impedance and population at the origin. +#' The units are in 'supply' (e.g., jobs, school seats). +#' - If `detailed_results = FALSE`, outputs are aggregated and returned by origin. +#' - If `detailed_results = TRUE`, OD-level flows are returned. Summing flows +#' by origin equals the aggregated result. +#' +#' - `active = FALSE` (*passive accessibility*, the notion of market potential): +#' constrains demand. Results represent the total number of +#' **population** (demand) that can reach each destination based on +#' region-relative travel impedance and opportunities at the destination. The +#' units are in 'demand' (e.g., population). +#' - If `detailed_results = FALSE`, outputs are aggregated by destination. +#' - If `detailed_results = TRUE`, OD-level flows are returned. Summing flows +#' by destination equals the aggregated result. +#' +#' **Use cases**: +#' - Active accessibility (aggregated): +#' "How many jobs can be reached from origin zone A given its region-relative +#' travel impedance and demand?" +#' +#' - Active accessibility (flow-level): +#' "How many jobs can be reached by flow A->1 given A->1's region-relative +#' travel impedance and demand?" +#' +#' - Passive accessibility (aggregated): +#' "How many people can reach destination zone 1 given its region-relative +#' travel impedance and supply?" +#' +#' - Passive accessibility (flow-level): +#' "How many people are reached by flow 1->A given 1->A's region-relative +#' travel impedance and supply?" +#' +#' **NOTE:** the active form of this measure yields equivalent results to the +#' `spatial_availability()` function, through different logic. #' -#' Sum of accessibility equals total opportunities (supply) in the region. -#' It allocates total opportunities in the region proportionally based on travel -#' impedance. Uses the logic of a total ~(or unconstrained by Wilson's terms)~ -#' constraint. Returns values as either `demand` or `supply`. When -#' `active = FALSE` (market potential variant) is also available. +#' ## Doubly Constrained Accessibility #' -#' ## Singly constrained accessibility +#' Allocates flows between origins and destinations using Wilson's *doubly-constrained* +#' gravity model \insertCite{wilson1971family;textual}{accessibility}. #' -#' Allocates opportunities at each destination proportionally based on travel -#' impedance and population at the origin. Uses the logic of single constraint -#' from \insertCite{wilson1971family;textual}{accessibility}. Returns values as -#' either 'demand' or 'supply'. Supply-constrained (destination totals fixed) -#' when `market_potential = FALSE`. In either case, totals match either the -#' demand at each origin or supply at each destination, depending on variant. -#' This is equivalent to the `spatial_availability()` function. +#' The model uses iterative proportional fitting to update balancing factors +#' (`A_i` for origins and `B_j` for destinations) until convergence. This guarantees +#' that flows satisfy both marginals while being weighted by travel impedance. #' +#' **Requirements**: +#' - Both `demand` and `supply` must be provided. +#' - Unlike `total` and `singly`, `doubly` requires the sum of demand and supply +#' to match; otherwise, the model will not converge. #' -#' ## Doubly constrained accessibility +#' **Interpretation**: +#' - When `detailed_results = TRUE`, results include OD-level flows (`flow`) +#' along with balancing factors (`A_i`, `B_j`) and travel impedance weights. +#' The resulting flows represent the distribution of demand and supply across +#' all origin-destination pairs. NOTE: OD flows are in flow units (jointly +#' determined by demand and supply). #' -#' Calculates accessibility using doubly-constrained gravity model of -#' \insertCite{wilson1971family;textual}{accessibility}. This measure allocates -#' flows between origins and destinations such that origin totals equal demand -#' and destination totals equal supply. Iterative proportional fitting updates -#' (A_i) and (B_j) until convergence. This ensures that row sums equal (O_i) -#' (demand) and column sums equal (D_j) (supply). Note, only OD-level outputs -#' are available (as aggregate outputs just match inputs). +#' - When `detailed_results = FALSE`, results are not returned. As the +#' aggregated outputs simply return the supply at destinations or demand at +#' origin that was fed into `demand` and `supply` parameters. +#' `detailed_results = TRUE` should only be used. +#' +#' **Use cases**: +#' - flow-level: +#' "What is the count of A->1 flows given A->1's region-relative travel impedance, +#' demand and supply?" #' #' @references #' \insertAllCited{} @@ -78,7 +220,8 @@ #' 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) +#' # Total-constrained (active accessibility, aggregated): returns units of +#' # accessible supply by origin (requires supply) #' constrained_accessibility( #' constraint = "total", #' travel_matrix = travel_matrix, @@ -87,10 +230,26 @@ #' decay_function = decay_exponential(0.1), #' demand = NULL, #' supply = "jobs", -#' active = FALSE +#' active = TRUE, +#' detailed_results = FALSE #' ) #' -#' # Singly-constrained (demand-side) +#' # Total-constrained (passive accessibility, aggregated): returns units of +#' # accessible demand by destination (requires demand) +#' constrained_accessibility( +#' constraint = "total", +#' travel_matrix = travel_matrix, +#' land_use_data = land_use_data, +#' travel_cost = "travel_time", +#' decay_function = decay_exponential(0.1), +#' demand = "population", +#' supply = NULL, +#' active = FALSE, +#' detailed_results = FALSE +#' ) +#' +#' # Singly-constrained (active accessibility, aggregated): returns units of +#' # accessible supply by origin (requires supply and demand) #' constrained_accessibility( #' constraint = "singly", #' travel_matrix = travel_matrix, @@ -99,10 +258,14 @@ #' decay_function = decay_exponential(0.1), #' demand = "population", #' supply = "jobs", -#' active = TRUE +#' active = TRUE, +#' detailed_results = FALSE #' ) #' -#' # Doubly-constrained: use a small toy dataset with matching totals +#' # Doubly-constrained: returns units of flow (requires both demand and supply +#' # (totals that match) and `detailed_results = TRUE`) +#' +#' # Using 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")) #' ) @@ -120,7 +283,8 @@ #' travel_cost = "travel_time", #' decay_function = decay_exponential(0.1), #' demand = "population", -#' supply = "jobs" +#' supply = "jobs", +#' detailed_results = TRUE #' ) #' #' @export @@ -131,7 +295,7 @@ constrained_accessibility <- function(constraint, decay_function, demand = NULL, supply = NULL, - active = TRUE, + active = NULL, error_threshold = 0.001, improvement_threshold = 1e-6, max_iterations = 1000, @@ -150,9 +314,6 @@ constrained_accessibility <- function(constraint, if (constraint == "doubly") { - if (!is.null(active)) { - stop("For 'doubly', the 'active' parameter is not used.") - } return(doubly_constrained( travel_matrix = travel_matrix, land_use_data = land_use_data, @@ -168,38 +329,46 @@ constrained_accessibility <- function(constraint, detailed_results = detailed_results )) } - + # For 'total' and 'singly' we require active to be TRUE or FALSE if (is.null(active)) { stop(sprintf("For '%s', active must be TRUE or FALSE.", constraint)) } - if (constraint == "total") { - return(total_constrained( + # Validate the relevant side (opportunity column) only + # 'active = TRUE' -> supply-side is returned; 'active = FALSE' -> demand-side is returned + opportunity_col <- if (active) supply else demand + assert_land_use_data(land_use_data, travel_matrix, opportunity = opportunity_col) + + + 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, + active = active, group_by = group_by, fill_missing_ids = fill_missing_ids, - detailed_results = detailed_results, - active = active, - demand = if (!active) demand else NULL, - supply = if (active) supply else NULL + detailed_results = detailed_results )) } - if (constraint == "singly") { - return(singly_constrained( + if (constraint == "total") { + return(total_constrained( travel_matrix = travel_matrix, land_use_data = land_use_data, travel_cost = travel_cost, decay_function = decay_function, - demand = demand, - supply = supply, - active = active, group_by = group_by, fill_missing_ids = fill_missing_ids, - detailed_results = detailed_results + detailed_results = detailed_results, + active = active, + demand = if (!active) demand else NULL, + supply = if (active) supply else NULL )) } + + } diff --git a/R/doubly_constrained.R b/R/doubly_constrained.R index 659acaa..59c9cd1 100644 --- a/R/doubly_constrained.R +++ b/R/doubly_constrained.R @@ -2,8 +2,8 @@ #' #' 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. This is an -#' internal helper function used by [constrained_accessibility()] when +#' origin totals equal demand and destination totals equal supply. Note: this is +#' an internal helper function used by [constrained_accessibility()] when #' `constraint = "doubly"`. #' #' @name doubly_constrained @@ -31,9 +31,11 @@ doubly_constrained <- function(travel_matrix, error_threshold = 0.001, improvement_threshold = 1e-6, max_iterations = 1000, + active, group_by = character(0), fill_missing_ids = TRUE, - detailed_results = FALSE) { + detailed_results = FALSE + ) { # Validate inputs checkmate::assert_string(demand) checkmate::assert_string(supply) diff --git a/R/singly_constrained.R b/R/singly_constrained.R index 4b1719c..61ba322 100644 --- a/R/singly_constrained.R +++ b/R/singly_constrained.R @@ -1,9 +1,9 @@ #' 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'. This is an internal -#' helper function used by [constrained_accessibility()] when `constraint = "singly"`. +#' impedance and population at the origin. Uses the logic of Wilson's single +#' constraint. Returns values in the unit of 'demand'. +#' Internal helper function used by [constrained_accessibility()] when `constraint = "singly"`. #' #' @name singly_constrained #' @keywords internal @@ -26,7 +26,7 @@ singly_constrained <- function(travel_matrix, decay_function, demand, supply, - active = TRUE, + active, group_by = character(0), fill_missing_ids = TRUE, detailed_results = FALSE) { @@ -54,10 +54,10 @@ singly_constrained <- function(travel_matrix, data <- apply_gravity_measure(data, decay_function, travel_cost) if (active) { - # Supply-constrained (returns the number of supply or JOBS) + # Supply-constrained (returns units of supply).. V_ij = (O_i f(c_ij) / sum_i O_i f(c_ij)) * D_j 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[, kappa_singly := data.table::fifelse(denom_j > 0, weighted_demand / denom_j, 0)] data[, singly_access := kappa_singly * get(supply)] if (detailed_results) { @@ -66,7 +66,7 @@ singly_constrained <- function(travel_matrix, to_id, kappa_singly, supply = singly_access, - B_j = 1 / denom_j + B_j = data.table::fifelse(denom_j > 0, 1 / denom_j, 0) )] } else { access <- data[, .(supply = sum(singly_access)), by = c("from_id", group_by)] @@ -74,10 +74,10 @@ singly_constrained <- function(travel_matrix, } } else { - # Demand-constrained (returns the number of demand or POPULATION) + # Demand-constrained (returns units of demand).... M_ij = (D_j f(c_ij) / sum_j D_j f(c_ij)) * O_i 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[, hatkappa_singly := data.table::fifelse(denom_i > 0, weighted_supply / denom_i, 0)] data[, singly_access := hatkappa_singly * get(demand)] if (detailed_results) { @@ -86,7 +86,7 @@ singly_constrained <- function(travel_matrix, to_id, hatkappa_singly, demand = singly_access, - A_i = 1 / denom_i + A_i = data.table::fifelse(denom_i > 0, 1 / denom_i, 0) )] } else { access <- data[, .(demand = sum(singly_access)), by = c("to_id", group_by)] diff --git a/R/total_constrained.R b/R/total_constrained.R index 6d4b134..69869b3 100644 --- a/R/total_constrained.R +++ b/R/total_constrained.R @@ -1,9 +1,12 @@ #' Total constrained accessibility #' #' 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'. This is an internal -#' helper function used by [constrained_accessibility()] when `constraint = "total"`. +#' impedance. Uses the logic of a total (or unconstrained by Wilson's terms) +#' constraint. Returns values in units of 'supply' (i.e., opportunities) if +#' `active = TRUE` and returns values in units of demand' (i.e., population) +#' if `active = FALSE`. +#' +#' Internal helper used by [constrained_accessibility()] when `constraint = "total"`. #' #' @name total_constrained #' @keywords internal @@ -29,25 +32,17 @@ total_constrained <- function(travel_matrix, group_by = character(0), fill_missing_ids = TRUE, detailed_results = FALSE, - active = TRUE, + active, 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(active, len = 1, any.missing = FALSE) - if (isFALSE(active)) { - if (is.null(demand) || !is.null(supply)) { - stop("For active = FALSE, 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 (active) { + #active accessibility if (is.null(supply) || !is.null(demand)) { stop("For active = TRUE, supply must be specified and demand must be NULL.") } @@ -58,6 +53,18 @@ total_constrained <- function(travel_matrix, total_col <- "K_total" result_col <- "supply" opportunity_col <- supply + } else { + #passive accessibility + if (is.null(demand) || !is.null(supply)) { + stop("For active = FALSE, 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 } assert_decay_function(decay_function) @@ -87,7 +94,10 @@ total_constrained <- function(travel_matrix, 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]]) + + ids_for_total <- unique(data[[merge_id]]) # e.g., A,B,C when active=TRUE + total_opportunities_region <- land_use_data[id %in% ids_for_total, sum(get(opportunity_col), na.rm = TRUE)] + data[, constrained_opportunity := get(kappa_col) * total_opportunities_region] if (detailed_results) { @@ -101,6 +111,7 @@ total_constrained <- function(travel_matrix, hatK_total = total_opportunities_region / total_weighted_system )] } else { + # supply-side details (active accessibility) access <- data[, .( from_id, to_id, @@ -111,12 +122,14 @@ total_constrained <- function(travel_matrix, )] } } else { - if (!active) { - 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 { + if (active) { + # supply aggregated by origin (from_id) 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)) + } else { + # demand aggregated by destination (to_id) + 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)) } } diff --git a/man/constrained_accessibility.Rd b/man/constrained_accessibility.Rd index 6258311..da66969 100644 --- a/man/constrained_accessibility.Rd +++ b/man/constrained_accessibility.Rd @@ -12,7 +12,7 @@ constrained_accessibility( decay_function, demand = NULL, supply = NULL, - active = TRUE, + active = NULL, error_threshold = 0.001, improvement_threshold = 1e-06, max_iterations = 1000, @@ -55,20 +55,21 @@ that will be considered.} \item{active}{A logical. When \code{TRUE}, the function calculates active accessibility (the quantity of opportunities that can be reached from -a given origin). when \code{FALSE}, it calculates passive accessibility (by +a given origin). When \code{FALSE}, it calculates passive accessibility (by how many people each destination can be reached), which is equivalent to the notion of market potential. This parameter only works for \code{constraint} types \code{"total"} and \code{"singly"}. Ignored for \code{constraint = "doubly"}.} \item{error_threshold}{Numeric. Convergence criterion used only for -doubly-constrained case.} +calibration in the doubly-constrained case (\code{constraint = "doubly"}).} \item{improvement_threshold}{Numeric. Convergence criterion for improvement -used only for doubly-constrained case.} +used only for calibration in the doubly-constrained case +(\code{constraint = "doubly"}).} -\item{max_iterations}{Integer. Maximum iterations used only for -doubly-constrained calibration.} +\item{max_iterations}{Integer. Maximum iterations used only for calibration +in the doubly-constrained case (\code{constraint = "doubly"}).} \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 @@ -90,51 +91,208 @@ performance penalty.} \item{detailed_results}{Logical. Whether to return detailed OD-level results.} } \description{ -Calculates accessibility using constrained gravity models, as proposed in -\insertCite{soukhov2025family;textual}{accessibility}: -\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. -Please see the Details section for more information. -} +Calculates accessibility using constraints, as proposed in +\insertCite{soukhov2025family;textual}{accessibility}. Accessibility is +conceptualised as potential' spatial interaction. The results are in units of +opportunities in the region (i.e., active accessibility \code{active = TRUE}) or +in units of population in the region (if \code{active = FALSE} reflecting passive +accessibility). Results are presented as a sum of a zone or presented as an +origin-destination flow (if \code{detailed_results = TRUE}). Please see Details +section for more information. This function is generic over any kind of numeric travel cost, such as distance, time and money. } +\details{ +The following three constraint cases (along with their \code{active} and passive +variants) are included: +\itemize{ +\item \code{"total"}: Allocates the system-wide total proportionally based on travel +impedance. +\item If \code{active = TRUE}, results are in units of \strong{opportunities} (supply) +accessible by the origin. Only \code{supply} can be passed, \code{demand} must be NULL. +\item If \code{active = FALSE}, results are in units of \strong{population} (demand) +accessible by the destination. Only \code{demand} can be passed, \code{supply} must +be NULL. +\item If \code{detailed_results = TRUE}, accessibility (in units corresponding to the +logical for \code{active}) is presented as a flow along with intermediates. +\item \code{"singly"}: Applies a single constraint, allocating one side of the marginal +proportionally based on the other marginal and travel impedance: +\item If \code{active = TRUE}, returns origin-side results (how much supply is +accessible from each origin). Outputs are units of supply. \code{demand} and +\code{supply} must be passed. +\item If \code{active = FALSE}, returns destination-side results (how much demand is +accessible from each destination). Outputs are units of demand. \code{demand} and +\code{supply} must be passed. +\item If \code{detailed_results = TRUE}, accessibility (in units corresponding to the +logical for \code{active}) is presented as a flow along with intermediates. +\item \code{"doubly"}: Allocates flows so supply at each destination matches demand at +each origin. +\item OD flows are calibrated to both marginals. OD flows are calibrated to both +marginals using iterative proportional fitting. The sum of \code{demand} and +\code{supply} must match; otherwise, the function will not converge. +\item \code{active} and \code{detailed_results} must be NULL. Since supply must match +demand, their units are the same and there is no distinction between 'active' +and 'passive' notions. +} + +Please see the Details section for more information. +} \section{Details}{ This function covers the family of constrained accessibility measures proposed in \insertCite{soukhov2025family;textual}{accessibility}. -\subsection{Total constrained accessibility}{ +\subsection{Total Constrained Accessibility}{ + +Allocates the total system-wide quantity proportionally based on travel +impedance between origins and destinations. This measure uses the logic of a +total ~(or 'unconstrained' by Wilson's terms)~ constraint. + +Use this measure when the total quantity of \strong{supply} OR \strong{demand} in the +system is known and representing accessibility as a proportion of this total +is meaningful. + +\strong{Requirements}: +\itemize{ +\item Either \code{demand} or \code{supply} must be provided (cannot provide both). +} + +\strong{Interpretation}: +\itemize{ +\item \code{active = TRUE} (\emph{active accessibility}): Results represent the total +number of \strong{opportunities} (supply) accessible from each origin based on +region-relative travel impedance. The units are in 'supply' (e.g., jobs, +school seats). +\itemize{ +\item If \code{detailed_results = FALSE}, outputs are aggregated and returned by +origin. +\item If \code{detailed_results = TRUE}, OD-level flows are returned. Summing flows +by origin equals the aggregated result. +} +\item \code{active = FALSE} (\emph{passive accessibility}, the notion of market potential): +Results represent the total number of \strong{population} (demand) that can reach +each destination based on region-relative travel impedance. The units are in +'demand' (e.g., population). +\itemize{ +\item If \code{detailed_results = FALSE}, outputs are aggregated by destination. +\item If \code{detailed_results = TRUE}, OD-level flows are returned. Summing flows +by destination equals the aggregated result. +} +} + +\strong{Use cases}: +\itemize{ +\item Active accessibility (aggregated): +"How many jobs can be reached from origin zone A given its region-relative +travel impedance?" +\item Active accessibility (flow-level): +"How many jobs can be reached by flow A->1 given A->1's region-relative +travel impedance?" +\item Passive accessibility (aggregated): +"How many people can reach destination zone 1 given its region-relative +travel impedance?" +\item Passive accessibility (flow-level): +"How many people are reached by flow 1->A given 1->A's region-relative +travel impedance?" +} +} + +\subsection{Singly Constrained Accessibility}{ + +Allocates opportunities at each destination (or population at each origin) +proportionally based on travel impedance and the opposite marginal. This +measure uses the logic of single constraint from +\insertCite{wilson1971family;textual}{accessibility}. + +Use this measure when modeling \strong{competition}, where both demand and supply +are conceptualised to influence accessibility but only one side is fixed. +The measure distributes flows so that totals match the constrained side +while weighting by travel impedance and the unconstrained side. + +\strong{Requirements}: +\itemize{ +\item Both \code{demand} and \code{supply} must be provided (the logical for \code{active} +determines if either demand or supply is constrained). +} + +\strong{Interpretation}: +\itemize{ +\item \code{active = TRUE} (\emph{active accessibility}): constrains supply. Results +represent the total number of \strong{opportunities} (supply) accessible from each +origin based on region-relative travel impedance and population at the origin. +The units are in 'supply' (e.g., jobs, school seats). +\itemize{ +\item If \code{detailed_results = FALSE}, outputs are aggregated and returned by origin. +\item If \code{detailed_results = TRUE}, OD-level flows are returned. Summing flows +by origin equals the aggregated result. +} +\item \code{active = FALSE} (\emph{passive accessibility}, the notion of market potential): +constrains demand. Results represent the total number of +\strong{population} (demand) that can reach each destination based on +region-relative travel impedance and opportunities at the destination. The +units are in 'demand' (e.g., population). +\itemize{ +\item If \code{detailed_results = FALSE}, outputs are aggregated by destination. +\item If \code{detailed_results = TRUE}, OD-level flows are returned. Summing flows +by destination equals the aggregated result. +} +} + +\strong{Use cases}: +\itemize{ +\item Active accessibility (aggregated): +"How many jobs can be reached from origin zone A given its region-relative +travel impedance and demand?" +\item Active accessibility (flow-level): +"How many jobs can be reached by flow A->1 given A->1's region-relative +travel impedance and demand?" +\item Passive accessibility (aggregated): +"How many people can reach destination zone 1 given its region-relative +travel impedance and supply?" +\item Passive accessibility (flow-level): +"How many people are reached by flow 1->A given 1->A's region-relative +travel impedance and supply?" +} -Sum of accessibility equals total opportunities (supply) in the region. -It allocates total opportunities in the region proportionally based on travel -impedance. Uses the logic of a total ~(or unconstrained by Wilson's terms)~ -constraint. Returns values as either \code{demand} or \code{supply}. When -\code{active = FALSE} (market potential variant) is also available. +\strong{NOTE:} the active form of this measure yields equivalent results to the +\code{spatial_availability()} function, through different logic. } -\subsection{Singly constrained accessibility}{ +\subsection{Doubly Constrained Accessibility}{ + +Allocates flows between origins and destinations using Wilson's \emph{doubly-constrained} +gravity model \insertCite{wilson1971family;textual}{accessibility}. -Allocates opportunities at each destination proportionally based on travel -impedance and population at the origin. Uses the logic of single constraint -from \insertCite{wilson1971family;textual}{accessibility}. Returns values as -either 'demand' or 'supply'. Supply-constrained (destination totals fixed) -when \code{market_potential = FALSE}. In either case, totals match either the -demand at each origin or supply at each destination, depending on variant. -This is equivalent to the \code{spatial_availability()} function. +The model uses iterative proportional fitting to update balancing factors +(\code{A_i} for origins and \code{B_j} for destinations) until convergence. This guarantees +that flows satisfy both marginals while being weighted by travel impedance. + +\strong{Requirements}: +\itemize{ +\item Both \code{demand} and \code{supply} must be provided. +\item Unlike \code{total} and \code{singly}, \code{doubly} requires the sum of demand and supply +to match; otherwise, the model will not converge. } -\subsection{Doubly constrained accessibility}{ +\strong{Interpretation}: +\itemize{ +\item When \code{detailed_results = TRUE}, results include OD-level flows (\code{flow}) +along with balancing factors (\code{A_i}, \code{B_j}) and travel impedance weights. +The resulting flows represent the distribution of demand and supply across +all origin-destination pairs. NOTE: OD flows are in flow units (jointly +determined by demand and supply). +\item When \code{detailed_results = FALSE}, results are not returned. As the +aggregated outputs simply return the supply at destinations or demand at +origin that was fed into \code{demand} and \code{supply} parameters. +\code{detailed_results = TRUE} should only be used. +} -Calculates accessibility using doubly-constrained gravity model of -\insertCite{wilson1971family;textual}{accessibility}. This measure allocates -flows between origins and destinations such that origin totals equal demand -and destination totals equal supply. Iterative proportional fitting updates -(A_i) and (B_j) until convergence. This ensures that row sums equal (O_i) -(demand) and column sums equal (D_j) (supply). Note, only OD-level outputs -are available (as aggregate outputs just match inputs). +\strong{Use cases}: +\itemize{ +\item flow-level: +"What is the count of A->1 flows given A->1's region-relative travel impedance, +demand and supply?" +} } } @@ -144,7 +302,8 @@ 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) +# Total-constrained (active accessibility, aggregated): returns units of +# accessible supply by origin (requires supply) constrained_accessibility( constraint = "total", travel_matrix = travel_matrix, @@ -153,10 +312,26 @@ constrained_accessibility( decay_function = decay_exponential(0.1), demand = NULL, supply = "jobs", - active = FALSE + active = TRUE, + detailed_results = FALSE +) + +# Total-constrained (passive accessibility, aggregated): returns units of +# accessible demand by destination (requires demand) +constrained_accessibility( + constraint = "total", + travel_matrix = travel_matrix, + land_use_data = land_use_data, + travel_cost = "travel_time", + decay_function = decay_exponential(0.1), + demand = "population", + supply = NULL, + active = FALSE, + detailed_results = FALSE ) -# Singly-constrained (demand-side) +# Singly-constrained (active accessibility, aggregated): returns units of +# accessible supply by origin (requires supply and demand) constrained_accessibility( constraint = "singly", travel_matrix = travel_matrix, @@ -165,10 +340,14 @@ constrained_accessibility( decay_function = decay_exponential(0.1), demand = "population", supply = "jobs", - active = TRUE + active = TRUE, + detailed_results = FALSE ) -# Doubly-constrained: use a small toy dataset with matching totals +# Doubly-constrained: returns units of flow (requires both demand and supply +# (totals that match) and `detailed_results = TRUE`) + +# Using 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")) ) @@ -186,7 +365,8 @@ constrained_accessibility( travel_cost = "travel_time", decay_function = decay_exponential(0.1), demand = "population", - supply = "jobs" + supply = "jobs", + detailed_results = TRUE ) } diff --git a/man/doubly_constrained.Rd b/man/doubly_constrained.Rd index d0c4268..1a62da9 100644 --- a/man/doubly_constrained.Rd +++ b/man/doubly_constrained.Rd @@ -10,8 +10,8 @@ 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. This is an -internal helper function used by \code{\link[=constrained_accessibility]{constrained_accessibility()}} when +origin totals equal demand and destination totals equal supply. Note: this is +an internal helper function used by \code{\link[=constrained_accessibility]{constrained_accessibility()}} when \code{constraint = "doubly"}. } \examples{ diff --git a/man/singly_constrained.Rd b/man/singly_constrained.Rd index 8dcd36e..6924310 100644 --- a/man/singly_constrained.Rd +++ b/man/singly_constrained.Rd @@ -8,9 +8,9 @@ A \code{data.table}/\code{data.frame} with results (structure mirrors the wrappe } \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'. This is an internal -helper function used by \code{\link[=constrained_accessibility]{constrained_accessibility()}} when \code{constraint = "singly"}. +impedance and population at the origin. Uses the logic of Wilson's single +constraint. Returns values in the unit of 'demand'. +Internal helper function used by \code{\link[=constrained_accessibility]{constrained_accessibility()}} when \code{constraint = "singly"}. } \examples{ NULL diff --git a/man/total_constrained.Rd b/man/total_constrained.Rd index 7535103..d731b0c 100644 --- a/man/total_constrained.Rd +++ b/man/total_constrained.Rd @@ -8,9 +8,13 @@ A \code{data.table}/\code{data.frame} with results (structure mirrors the wrappe } \description{ 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'. This is an internal -helper function used by \code{\link[=constrained_accessibility]{constrained_accessibility()}} when \code{constraint = "total"}. +impedance. Uses the logic of a total (or unconstrained by Wilson's terms) +constraint. Returns values in units of 'supply' (i.e., opportunities) if +\code{active = TRUE} and returns values in units of demand' (i.e., population) +if \code{active = FALSE}. +} +\details{ +Internal helper used by \code{\link[=constrained_accessibility]{constrained_accessibility()}} when \code{constraint = "total"}. } \examples{ NULL diff --git a/tests/test_rafa/testing_constrained.Rmd b/tests/test_rafa/testing_constrained.Rmd index e1b70d0..f2a23da 100644 --- a/tests/test_rafa/testing_constrained.Rmd +++ b/tests/test_rafa/testing_constrained.Rmd @@ -13,15 +13,15 @@ devtools::load_all(".") ## TOT CONST ```{r} # Toy travel matrix and land use data -travel_matrix <- expand.grid(from_id = c("1", "2", "3"), +travel_matrix <- expand.grid(from_id = c("A", "B", "C"), to_id = c("1", "2", "3")) |> - dplyr::mutate(travel_time = c(10, 30, 15, 30, 10, 25, 15, 25, 10)) |> + dplyr::mutate(travel_time = c(15, 30, 100, 30, 15, 100, 100, 100, 15)) |> 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 + id = c("A", "B", "C","1", "2", "3"), + population = c(50000, 150000, 10000), # supply + jobs = c(100000, 100000, 10000) # demand ) # Decay function @@ -30,7 +30,7 @@ decay <- decay_exponential(decay_value = 0.1) ```{r} -agg_false <- constrained_accessibility( +active_det_false <- constrained_accessibility( constraint = "total", travel_matrix, land_use_data, travel_cost = "travel_time", @@ -38,11 +38,14 @@ agg_false <- constrained_accessibility( demand = NULL, supply = "jobs", detailed_results = FALSE, - return_demand_side = FALSE) + active = TRUE) -agg_false +active_det_false +``` -agg_true <- constrained_accessibility( + +```{r} +passive_det_false <- constrained_accessibility( constraint = "total", travel_matrix = travel_matrix, land_use_data = land_use_data, @@ -51,12 +54,15 @@ agg_true <- constrained_accessibility( demand = "population", supply = NULL, detailed_results = FALSE, - return_demand_side = TRUE + active = FALSE ) -agg_true +passive_det_false +``` + -det_false <- constrained_accessibility( +```{r} +active_det_true <- constrained_accessibility( constraint = "total", travel_matrix, land_use_data, @@ -65,85 +71,106 @@ det_false <- constrained_accessibility( demand = NULL, supply = "jobs", detailed_results = TRUE, - return_demand_side = FALSE + active = TRUE ) -det_false -det_true <- constrained_accessibility(constraint = "total", travel_matrix, land_use_data, +active_det_true +``` + + +```{r} +passive_det_true <- constrained_accessibility(constraint = "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 + active = FALSE) +passive_det_true ``` ```{r} -det_true$demand |> sum() -det_false$supply |> sum() -agg_true$demand |> sum() -agg_false$supply |> sum() +passive_det_true$demand |> sum() #POPULATION 210,000 +active_det_true$supply |> sum() #jobs 210,000 +passive_det_false$demand |> sum() #POPULATION 210,000 +active_det_false$supply |> sum() #jobs 210,000 -det_false$kappa_total |> sum() -det_true$hatkappa_total |> sum() +active_det_true$kappa_total |> sum() +passive_det_true$hatkappa_total |> sum() ``` ## SINGLE CONST ```{r} -agg_false <- constrained_accessibility(constraint = "singly", travel_matrix, land_use_data, +active_det_false <- constrained_accessibility(constraint = "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 + active = TRUE) +active_det_false +``` +Should be: +A 66833.466 +B 133203.363 +C 9963.171 -agg_true <- constrained_accessibility(constraint = "singly", travel_matrix, land_use_data, +```{r} +passive_det_false <- constrained_accessibility(constraint = "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 + active = FALSE) +passive_det_false +``` +Should be: +1 68261.682 +2 131775.520 +3 9962.798 -det_false <- constrained_accessibility(constraint = "singly", travel_matrix, land_use_data, +```{r} +active_det_true <- constrained_accessibility(constraint = "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 - + active = TRUE) +active_det_true +``` +e.g., for summing As... 59900.642536+6922.691048+10.132187 = 66833.47 -det_true <- constrained_accessibility(constraint = "singly", travel_matrix, land_use_data, +```{r} +passive_det_true <- constrained_accessibility(constraint = "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 + active = FALSE) +passive_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) +passive_det_true$demand |> sum() #POPULATION 210,000 +active_det_true$supply |> sum() #jobs 210,000 +passive_det_false$demand |> sum() #POPULATION 210,000 +active_det_false$supply |> sum() #jobs 210,000 + +active_det_true$kappa_singly |> sum() +passive_det_true$hatkappa_singly |> sum() ``` ```{r} -det_false$kappa_singly |> sum() #should be the number of destinations -det_false[, sum(kappa_singly), by = to_id] #should all be 1 +active_det_true$kappa_singly |> sum() #should be the number of destinations +active_det_true[, 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 +passive_det_true$hatkappa_singly |> sum() #should be the number of origins +passive_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 @@ -160,70 +187,387 @@ land_use_data <- data.table( ) ``` -11/12 --> total_constrained works! I tested it for all decay scenarios (3, 2, 0.1). +## Total constrained--ACTIVE +Checking for decay_power scenario 3: ```{r} -constrained_accessibility(constraint="total", +active_det_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 = TRUE, - return_demand_side = FALSE) + detailed_results = FALSE, + active = TRUE) -constrained_accessibility(constraint="total", +active_det_true <- 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, + active = TRUE) +active_det_false +active_det_true + +active_det_false$supply |> sum() +active_det_true$supply |> sum() +active_det_true$kappa_total |> sum() +``` +Should be: +1 172.0653 +2 131.6267 +3 186.3080 + +Checking for decay_power scenario 2: +```{r} +active_det_false <- constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = NULL, + supply = "jobs", detailed_results = FALSE, - return_demand_side = FALSE) + active = TRUE) + +active_det_true <- constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = NULL, + supply = "jobs", + detailed_results = TRUE, + active = TRUE) +active_det_false +active_det_true + +active_det_false$supply |> sum() +active_det_true$supply |> sum() +active_det_true$kappa_total |> sum() +``` +Should be: +1 172.6721 +2 132.2474 +3 185.0805 + +Checking for decay_power scenario 0.1: +```{r} +active_det_false <- constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = NULL, + supply = "jobs", + detailed_results = FALSE, + active = TRUE) + +active_det_true <- constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = NULL, + supply = "jobs", + detailed_results = TRUE, + active = TRUE) +active_det_false +active_det_true +active_det_false$supply |> sum() +active_det_true$supply |> sum() +active_det_true$kappa_total |> sum() +``` +Should be: +1 164.0802 +2 160.6921 +3 165.2277 -constrained_accessibility(constraint="total", +## Total constrained--PASSIVE +Checking for decay_power scenario 3: +```{r} +passive_det_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) + detailed_results = FALSE, + active = FALSE) -constrained_accessibility(constraint="total", +passive_det_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 = TRUE, + active = FALSE) +passive_det_false +passive_det_true + +passive_det_false$demand |> sum() +passive_det_true$demand |> sum() +passive_det_true$hatkappa_total |> sum() +``` +Should be: +1 5.017774 +2 8.595749 +3 6.386477 + +Checking for decay_power scenario 2: +```{r} +passive_det_false <- constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = "population", + supply = NULL, detailed_results = FALSE, - return_demand_side = TRUE) + active = FALSE) + +passive_det_true <- constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = "population", + supply = NULL, + detailed_results = TRUE, + active = FALSE) +passive_det_false +passive_det_true + +passive_det_false$demand |> sum() +passive_det_true$demand |> sum() +passive_det_true$hatkappa_total |> sum() ``` +Should be: +1 5.446623 +2 7.986306 +3 6.567071 -11/12 --> Singly-constrained works! I tested it for all decay scenarios (3, 2, 0.1)! +Checking for decay_power scenario 0.1: ```{r} -constrained_accessibility(constraint = "singly", +passive_det_false <- constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = "population", + supply = NULL, + detailed_results = FALSE, + active = FALSE) + +passive_det_true <- constrained_accessibility(constraint="total", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = "population", + supply = NULL, + detailed_results = TRUE, + active = FALSE) +passive_det_false +passive_det_true + +passive_det_false$demand |> sum() +passive_det_true$demand |> sum() +passive_det_true$hatkappa_total |> sum() +``` +Should be: +1 6.598332 +2 6.717223 +3 6.684444 + +## Singly constrained--ACTIVE +Checking for decay_power scenario 3: +```{r} +active_det_false <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 3), + demand = "population", + supply = "jobs", + detailed_results = FALSE, + active = TRUE) + +active_det_true <- constrained_accessibility(constraint="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) + active = TRUE) +active_det_false +active_det_true -constrained_accessibility(constraint = "singly", +active_det_false$supply |> sum() +active_det_true$supply |> sum() +active_det_true$kappa_singly |> sum() #should be the number of destinations or origins + +round(active_det_false$supply / land_use_data$population, 6) +``` +Should be: +1 133.4687 +2 166.7813 +3 189.7499 + +And the per capita should be: +33.36718 16.67813 31.62499 + +Checking for decay_power scenario 2: +```{r} +active_det_false <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = "population", + supply = "jobs", + detailed_results = FALSE, + active = TRUE) + +active_det_true <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = "population", + supply = "jobs", + detailed_results = TRUE, + active = TRUE) +active_det_false +active_det_true + +active_det_false$supply |> sum() +active_det_true$supply |> sum() +active_det_true$kappa_singly|> sum() #should be the number of destinations or origins + +round(active_det_false$supply / land_use_data$population, 6) +``` +Should be: +1 122.2546 +2 185.0957 +3 182.6497 + +And the per capita should be: +30.56365 18.50957 30.44161 + +Checking for decay_power scenario 0.1: +```{r} +active_det_false <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = "population", + supply = "jobs", + detailed_results = FALSE, + active = TRUE) + +active_det_true <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = "population", + supply = "jobs", + detailed_results = TRUE, + active = TRUE) +active_det_false +active_det_true + +active_det_false$supply |> sum() +active_det_true$supply |> sum() +active_det_true$kappa_singly |> sum() #should be the number of destinations or origins + +round(active_det_false$supply / land_use_data$population, 6) +``` +Should be: +1 98.84766 +2 241.87721 +3 149.27513 + +And the per capita should be: +24.71192 24.18772 24.87919 + +## Singly constrained--PASSIVE +Checking for decay_power scenario 3: +```{r} +passive_det_false <- constrained_accessibility(constraint="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) -# + active = FALSE) + +passive_det_true <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 3), + demand = "population", + supply = "jobs", + detailed_results = TRUE, + active = FALSE) +passive_det_false +passive_det_true + +passive_det_false$demand |> sum() +passive_det_true$demand |> sum() +passive_det_true$hatkappa_singly |> sum() ``` -# Test 3. Testing results with doubly constrained toy example in 2025 paper +Checking for decay_power scenario 2: +```{r} +passive_det_false <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = "population", + supply = "jobs", + detailed_results = FALSE, + active = FALSE) + +passive_det_true <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 2), + demand = "population", + supply = "jobs", + detailed_results = TRUE, + active = FALSE) +passive_det_false +passive_det_true + +passive_det_false$demand |> sum() +passive_det_true$demand |> sum() +passive_det_true$hatkappa_singly |> sum() +``` + +Checking for decay_power scenario 0.1: +```{r} +passive_det_false <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = "population", + supply = "jobs", + detailed_results = FALSE, + active = FALSE) + +passive_det_true <- constrained_accessibility(constraint="singly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = "population", + supply = "jobs", + detailed_results = TRUE, + active = FALSE) +passive_det_false +passive_det_true + +passive_det_false$demand |> sum() +passive_det_true$demand |> sum() +passive_det_true$hatkappa_singly |> sum() +``` + +# 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"), @@ -237,42 +581,73 @@ land_use_data <- data.table( 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. + + +Checking for decay_power scenario 3: ```{r} -det <- - constrained_accessibility(constraint = "doubly", - #doubly_constrained( - travel_matrix, land_use_data, +active_det_true <- constrained_accessibility(constraint="doubly", + travel_matrix, land_use_data, travel_cost = "travel_time", - decay_function = decay_power(decay_value = 0.1), + decay_function = decay_power(decay_value = 3), demand = "population", supply = "jobs", detailed_results = TRUE, - return_demand_side = NULL) + active = TRUE) -det +active_det_true -agg <- - constrained_accessibility(constraint = "doubly", - #doubly_constrained( - travel_matrix, land_use_data, +active_det_true$flow |> sum() +active_det_true$kappa_doubly |> sum() #should be the number of destinations or origins +``` +Should be: +1 1 3.23585899 +1 2 0.01032226 +1 3 0.75565683 +2 1 2.13260167 +2 2 4.95932483 +2 3 2.90443915 +3 1 1.63153933 +3 2 0.03035291 +3 3 4.33990401 + +So... this matches table 10, but I accidentally flipped the columns. like 2->1 in the table should be 1->2. + + +Checking for decay_power scenario 2: +```{r} +active_det_true <- constrained_accessibility(constraint="doubly", + travel_matrix, land_use_data, travel_cost = "travel_time", - decay_function = decay_power(decay_value = 0.1), + decay_function = decay_power(decay_value = 2), demand = "population", supply = "jobs", - detailed_results = FALSE, - return_demand_side = NULL) -``` + detailed_results = TRUE, + active = TRUE) +active_det_true +active_det_true$flow |> sum() +active_det_true$kappa_doubly |> sum() #should be the number of destinations or origins +``` + +Checking for decay_power scenario 0.1: ```{r} -agg$flow |> sum() #should be double, -det$flow |> sum() +active_det_true <- constrained_accessibility(constraint="doubly", + travel_matrix, land_use_data, + travel_cost = "travel_time", + decay_function = decay_power(decay_value = 0.1), + demand = "population", + supply = "jobs", + detailed_results = TRUE, + active = TRUE) + +active_det_true -det$kappa_doubly |> sum() #should be the number of destinations or origins +active_det_true$flow |> sum() +active_det_true$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"), @@ -295,7 +670,7 @@ constrained_accessibility(constraint = "singly", demand = "population", supply = "jobs", detailed_results = FALSE, - return_demand_side = FALSE) + active = TRUE) spatial_availability(travel_matrix, land_use_data, @@ -316,7 +691,7 @@ result_singly <- constrained_accessibility(constraint = "singly", demand = "population", supply = "jobs", detailed_results = FALSE, - return_demand_side = FALSE) + active = TRUE) round(result_singly$supply / land_use_data$population, 6) diff --git a/tests/testthat/test-constrained-essentials.R b/tests/testthat/test-constrained-essentials.R index 405c382..ed0ab9d 100644 --- a/tests/testthat/test-constrained-essentials.R +++ b/tests/testthat/test-constrained-essentials.R @@ -33,24 +33,24 @@ test_that("total-constrained: preserves totals on supply/demand sides", { 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 + detailed_results = TRUE, active = TRUE ) 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 + detailed_results = FALSE, active = TRUE ) # 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 + detailed_results = TRUE, active = FALSE ) 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 + detailed_results = FALSE, active = FALSE ) expect_equal(sum(det_supply$supply), sum(lu$jobs)) @@ -68,24 +68,24 @@ test_that("singly-constrained: preserves totals and kappa properties", { 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 + detailed_results = TRUE, active = TRUE ) 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 + detailed_results = FALSE, active = TRUE ) # 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 + detailed_results = TRUE, active = FALSE ) 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 + detailed_results = FALSE, active = FALSE ) # Totals @@ -113,7 +113,7 @@ test_that("doubly-constrained: errors when totals mismatch", { 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 + detailed_results = TRUE, active = NULL ), "sum of origins must equal the sum of destinations" ) @@ -129,7 +129,7 @@ test_that("doubly-constrained: OD flows match marginals", { det <- constrained_accessibility("doubly", tm, lu, travel_cost = "travel_time", decay_function = decay, demand = "population", supply = "jobs", - detailed_results = TRUE, return_demand_side = NULL + detailed_results = TRUE, active = NULL ) # Robustly pick the numeric flow column diff --git a/tests/testthat/test-constrained-wrapper.R b/tests/testthat/test-constrained-wrapper.R index 22959e0..d78d961 100644 --- a/tests/testthat/test-constrained-wrapper.R +++ b/tests/testthat/test-constrained-wrapper.R @@ -20,7 +20,7 @@ tester <- function( decay_function = decay, demand = "population", supply = "jobs", - return_demand_side = TRUE, + active = FALSE, error_threshold = 0.001, improvement_threshold = 1e-6, max_iterations = 1000, @@ -62,7 +62,7 @@ test_that("constrained_accessibility: wrapper argument rules are enforced", { # doubly requires NULL (error if TRUE) expect_error( - tester(constraint = "doubly", return_demand_side = TRUE) + tester(constraint = "doubly", active = FALSE) ) })