diff --git a/DESCRIPTION b/DESCRIPTION index 8ed30bc..8638e8f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: cppRouting Type: Package Title: Algorithms for Routing and Solving the Traffic Assignment Problem -Version: 3.2 +Version: 4.0 Date: 2025-11-28 Authors@R: person(given = "Vincent",family = "Larmet",role = c("aut", "cre"),email = "larmet.vincent@gmail.com") Description: Calculation of distances, shortest paths and isochrones on weighted graphs using several variants of Dijkstra algorithm. @@ -15,15 +15,16 @@ Description: Calculation of distances, shortest paths and isochrones on weighted Method of Successive Averages, Frank-Wolfe algorithm (M. Fukushima (1984) ), Conjugate and Bi-Conjugate Frank-Wolfe algorithms (M. Mitradjieva, P. O. Lindberg (2012) ), - Algorithm-B (R. B. Dial (2006) ). + Algorithm-B (R. B. Dial (2006) ), + Customizable Contraction Hierarchies (J. Dibbelt, B. Strasser, D. Wagner (2014) ). License: GPL (>= 2) Encoding: UTF-8 LazyData: true Imports: Rcpp (>= 1.0.7), RcppParallel, RcppProgress, data.table LinkingTo: Rcpp, RcppParallel, RcppProgress SystemRequirements: GNU make -RoxygenNote: 7.2.1 URL: https://github.com/vlarmet/cppRouting Suggests: knitr, rmarkdown, igraph VignetteBuilder: knitr NeedsCompilation: yes +Config/roxygen2/version: 8.0.0 diff --git a/NAMESPACE b/NAMESPACE index c5a2240..56b4fb0 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -11,6 +11,7 @@ export(get_detour) export(cpp_contract) export(assign_traffic) export(get_aon) +export(cpp_customize) importFrom(Rcpp, evalCpp) importFrom(RcppParallel, RcppParallelLibs) importFrom("stats", "setNames") diff --git a/NEWS.md b/NEWS.md index 6012d73..411f682 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,19 @@ +cppRouting v4.0 +=============== +Major changes + +- Customizable Contraction Hierarchies (CCH) support, via `cpp_contract(..., customizable = TRUE)` and `cpp_customize()`. CCH separates preprocessing (topology) from customization (edge weights), so a prepared topology can be reused across many weight updates (e.g. congestion assignment) without a full re-contraction. +- `get_distance_pair`, `get_distance_matrix`, `get_path_pair` and `get_aon` all support CCH graphs, in addition to normal and classically-contracted graphs. +- `makegraph()`'s `aux` argument now accepts a `data.frame` with one or more named columns, in addition to a single vector. `aggregate_aux` in `get_distance_pair`/`get_distance_matrix` accepts `TRUE` (all columns) or a character vector of column names, and is now supported on normal, contracted (CH) and CCH graphs, and on graphs simplified with `cpp_simplify`. Aggregation runs inline in C++ during the shortest-path search itself, no path reconstruction needed. +- `assign_traffic(..., aon_method = "cch", cch = cch)` uses a prepared CCH topology for repeated traffic assignment iterations, avoiding a full re-contraction at each run. + +Minor changes + +- `cpp_simplify()` result is now classed (`cppRouting_simplified`), fixing downstream functions that silently failed on a simplified graph. +- `cpp_simplify()` now correctly aggregates multiple auxiliary columns across merged edges (previously single-column only). +- `cpp_contract(..., customizable = TRUE)` shows progress, like classical contraction. +- Various CCH correctness and performance fixes: `get_path_pair` on CCH graphs, elimination ordering (near-linear instead of quadratic), many-to-many distance matrix on CCH graphs. + cppRouting v3.2 =============== diff --git a/R/RcppExports.R b/R/RcppExports.R index c346789..6a07ba8 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -17,14 +17,42 @@ cpppathmat <- function(gfrom, gto, gw, nb, dict, keep, dep, arr, lim, setdif, al .Call(`_cppRouting_cpppathmat`, gfrom, gto, gw, nb, dict, keep, dep, arr, lim, setdif, algo, reverse) } -cppsimplify <- function(gfrom, gto, gw, nb, keep, rm_loop, iterate, progress) { - .Call(`_cppRouting_cppsimplify`, gfrom, gto, gw, nb, keep, rm_loop, iterate, progress) +cppsimplify <- function(gfrom, gto, gw, nb, keep, rm_loop, iterate, progress, aux) { + .Call(`_cppRouting_cppsimplify`, gfrom, gto, gw, nb, keep, rm_loop, iterate, progress, aux) } cppcontract <- function(gfrom, gto, gw, NbNodes, display_progress) { .Call(`_cppRouting_cppcontract`, gfrom, gto, gw, NbNodes, display_progress) } +cppcchprepare <- function(gfrom, gto, NbNodes, order, display_progress = FALSE) { + .Call(`_cppRouting_cppcchprepare`, gfrom, gto, NbNodes, order, display_progress) +} + +cppcchcustomize <- function(gfrom, gto, gw, NbNodes, rank, tail, head, rank_first_out, rank_adj_head, rank_adj_arc, input_arc, input_forward) { + .Call(`_cppRouting_cppcchcustomize`, gfrom, gto, gw, NbNodes, rank, tail, head, rank_first_out, rank_adj_head, rank_adj_arc, input_arc, input_forward) +} + +cppdistcch <- function(NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, dep, arr) { + .Call(`_cppRouting_cppdistcch`, NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, dep, arr) +} + +cppdistmatcch <- function(NbNodes, first_out, adj_head, adj_arc, forward, backward, dep, arr) { + .Call(`_cppRouting_cppdistmatcch`, NbNodes, first_out, adj_head, adj_arc, forward, backward, dep, arr) +} + +cpppathvaluescch <- function(gfrom, NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, dep, arr, values) { + .Call(`_cppRouting_cpppathvaluescch`, gfrom, NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, dep, arr, values) +} + +cpppathcch <- function(NbNodes, rank, gfrom, gto, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, dict, keep, dep, arr) { + .Call(`_cppRouting_cpppathcch`, NbNodes, rank, gfrom, gto, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, dict, keep, dep, arr) +} + +cppaoncchelim <- function(gfrom, gto, gw, NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, dep, arr, dem) { + .Call(`_cppRouting_cppaoncchelim`, gfrom, gto, gw, NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, dep, arr, dem) +} + cppdistC <- function(gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr, algo) { .Call(`_cppRouting_cppdistC`, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr, algo) } @@ -45,18 +73,34 @@ cpppadd <- function(gfrom, gto, gw, gadd, nb, dep, arr) { .Call(`_cppRouting_cpppadd`, gfrom, gto, gw, gadd, nb, dep, arr) } +cppmataddmulti <- function(gfrom, gto, gw, aux, nb, dep, arr) { + .Call(`_cppRouting_cppmataddmulti`, gfrom, gto, gw, aux, nb, dep, arr) +} + cppdistadd <- function(gfrom, gto, gw, gadd, nb, lat, lon, k, dep, arr, algo) { .Call(`_cppRouting_cppdistadd`, gfrom, gto, gw, gadd, nb, lat, lon, k, dep, arr, algo) } +cppdistaddmulti <- function(gfrom, gto, gw, aux, nb, lat, lon, k, dep, arr, algo) { + .Call(`_cppRouting_cppdistaddmulti`, gfrom, gto, gw, aux, nb, lat, lon, k, dep, arr, algo) +} + cppdistaddC <- function(orfrom, orto, orw, gadd, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr, algo) { .Call(`_cppRouting_cppdistaddC`, orfrom, orto, orw, gadd, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr, algo) } +cppdistaddCmulti <- function(orfrom, orto, orw, aux, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr) { + .Call(`_cppRouting_cppdistaddCmulti`, orfrom, orto, orw, aux, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr) +} + cppaddC <- function(orfrom, orto, orw, gadd, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr, algo) { .Call(`_cppRouting_cppaddC`, orfrom, orto, orw, gadd, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr, algo) } +cppmataddCmulti <- function(orfrom, orto, orw, aux, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr) { + .Call(`_cppRouting_cppmataddCmulti`, orfrom, orto, orw, aux, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr) +} + cppaon <- function(gfrom, gto, gw, nb, lat, lon, k, dep, arr, dem, algo) { .Call(`_cppRouting_cppaon`, gfrom, gto, gw, nb, lat, lon, k, dep, arr, dem, algo) } @@ -65,6 +109,10 @@ cppaonC <- function(orfrom, orto, orw, gfrom, gto, gw, nb, rank, shortf, shortt, .Call(`_cppRouting_cppaonC`, orfrom, orto, orw, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr, dem, algo) } +cpptrafficcch <- function(gfrom, gto, gw, gflow, gaux, gftt, galpha, gbeta, gcap, nb, dep, arr, dem, max_gap, max_it, method, rank, tail, head, first_out, adj_head, adj_arc, rank_first_out, rank_adj_head, rank_adj_arc, input_arc, input_forward, elimination_tree_parent, verbose) { + .Call(`_cppRouting_cpptrafficcch`, gfrom, gto, gw, gflow, gaux, gftt, galpha, gbeta, gcap, nb, dep, arr, dem, max_gap, max_it, method, rank, tail, head, first_out, adj_head, adj_arc, rank_first_out, rank_adj_head, rank_adj_arc, input_arc, input_forward, elimination_tree_parent, verbose) +} + cpptraffic <- function(gfrom, gto, gw, gflow, gaux, gftt, galpha, gbeta, gcap, nb, lat, lon, k, dep, arr, dem, max_gap, max_it, method, aon_method, contract, phast, verbose) { .Call(`_cppRouting_cpptraffic`, gfrom, gto, gw, gflow, gaux, gftt, galpha, gbeta, gcap, nb, lat, lon, k, dep, arr, dem, max_gap, max_it, method, aon_method, contract, phast, verbose) } diff --git a/R/assign_traffic.R b/R/assign_traffic.R index b3fe6fc..706f738 100644 --- a/R/assign_traffic.R +++ b/R/assign_traffic.R @@ -7,7 +7,8 @@ #' @param algorithm character. \code{msa}, \code{fw}, \code{cfw}, \code{bfw} or \code{dial}. Default to \code{bfw}. See details. #' @param max_gap Numeric. Relative gap to achieve. Default to 0.001. #' @param max_it Numeric. Maximum number of iterations. Default to \code{.Machine$integer.max} -#' @param aon_method Character.\code{d}, \code{bi}, \code{nba}, \code{cphast} or \code{cbi}. Default to \code{bi}. See details. +#' @param aon_method Character.\code{d}, \code{bi}, \code{nba}, \code{cphast}, \code{cbi} or \code{cch}. Default to \code{bi}. See details. +#' @param cch Optional CCH topology generated by \code{cpp_contract(Graph, customizable = TRUE)}. Use it with \code{aon_method = "cch"} to prepare the road topology once and reuse it across many congestion runs. #' @param constant numeric. Constant to maintain the heuristic function admissible in NBA* algorithm. Default to 1, when cost is expressed in the same unit than coordinates. See details #' @param dial_params List. Named list of hyperparameters for \code{dial} algorithm. See details. #' @param verbose Logical. If \code{TRUE} (default), progression is displayed. @@ -62,12 +63,20 @@ #' \item \code{d} : Dijkstra algorithm #' \item \code{cbi} : contraction hierarchies + bidirectional search #' \item \code{cphast} : contraction hierarchies + phast algorithm +#' \item \code{cch} : customizable contraction hierarchies with RoutingKit-style elimination-tree grouped queries. The topology is prepared once, then edge costs are customized at each assignment iteration. #' } #' These AON algorithm can be decomposed into two families, depending the sparsity of origin-destination matrix : \itemize{ #' \item recursive pairwise : \code{bi}, \code{nba} and \code{cbi}. Optimal for high sparsity. One-to-one algorithm is called N times, with N being the length of from. -#' \item recursive one-to-many : \code{d} and \code{cphast}. Optimal for dense matrix. One-to-many algorithm is called N times, with N being the number of unique from (or to) nodes +#' \item recursive one-to-many : \code{d}, \code{cphast} and \code{cch}. Optimal for dense or repeated-endpoint matrices. One-to-many algorithm is called N times, with N being the number of unique from (or to) nodes. #' } #' +#' CCH is most useful when the same graph topology is reused with changing edge +#' costs, for example during congestion assignment. In that case, run +#' \code{cch <- cpp_contract(graph, customizable = TRUE)} once, save it if needed, and pass +#' \code{cch = cch} with \code{aon_method = "cch"}. The assignment will +#' then customize the CCH at each iteration instead of rebuilding a full +#' contraction hierarchy. +#' #' For large instance, it may be appropriate to test different \code{aon_method} for few iterations and choose the fastest one for the final estimation. #' #' Hyperparameters for algorithm-b are : \itemize{ @@ -120,7 +129,7 @@ assign_traffic <- function(Graph, from, to, demand, algorithm = "bfw", max_gap = 0.001, max_it = .Machine$integer.max, - aon_method = "bi", constant = 1, dial_params = NULL, verbose = TRUE){ + aon_method = "bi", cch = NULL, constant = 1, dial_params = NULL, verbose = TRUE){ if (length(from)!=length(to) | length(from) != length(demand)) stop("From, to and demand have not the same length") demand <- as.numeric(demand) if (any(is.na(data.frame(from,to, demand)))) stop("NAs are not allowed in origin/destination trips") @@ -129,7 +138,7 @@ assign_traffic <- function(Graph, from, to, demand, algorithm = "bfw", max_gap = to <- to[ind] demand <- demand[ind] - if (length(Graph) != 5) stop("Graph object must be generated by makegraph() function") + if (!inherits(Graph, "cppRouting_graph")) stop("Graph object must be generated by makegraph() function") if (is.null(Graph$attrib$alpha) | is.null(Graph$attrib$beta)| is.null(Graph$attrib$cap)) stop("alpha, beta and capacity must be defined during graph construction") @@ -143,35 +152,30 @@ assign_traffic <- function(Graph, from, to, demand, algorithm = "bfw", max_gap = to_id <- Graph$dict$id[match(to,Graph$dict$ref)] if (!algorithm %in% c("msa", "fw", "cfw", "bfw", "dial")) stop("algorithm should be 'msa', 'fw', 'cfw', 'bfw' or 'dial'") - if (!aon_method %in% c("bi", "nba", "d", "cbi", "cphast")) stop("aon_method should be 'bi', 'nba', 'd', 'cbi' or 'cphast'") + if (!aon_method %in% c("bi", "nba", "d", "cbi", "cphast", "cch")) stop("aon_method should be 'bi', 'nba', 'd', 'cbi', 'cphast' or 'cch'") + if (algorithm == "dial" && aon_method == "cch") stop("CCH aon_method is only implemented for link-based algorithms: 'msa', 'fw', 'cfw' and 'bfw'") + if (!is.null(cch) && aon_method != "cch") stop("cch can only be used with aon_method = 'cch'") + if (aon_method == "cch" && !is.null(cch)){ + validate_cch_for_graph(cch, Graph) + } if (aon_method == "nba" & is.null(Graph$coords)){ aon_method <- "bi" message("nodes coordinates are not provided, using bidirectional Dijkstra") } - phast <- FALSE - - if (aon_method == "d"){ - aon <- ifelse(length(unique(from_id)) <= length(unique(to_id)), 0, 1) - contract <- FALSE - } - if (aon_method == "bi") { - aon <- 2 - contract <- FALSE - } - if (aon_method == "nba") { - aon <- 3 - contract <- FALSE - } - if (aon_method == "cbi") { - aon <- 2 - contract <- TRUE - } - if (aon_method == "cphast"){ - aon <- ifelse(length(unique(from_id)) <= length(unique(to_id)), 0, 1) - contract <- TRUE - phast <- TRUE - } + n_from <- length(unique(from_id)) + n_to <- length(unique(to_id)) + grouped_aon <- ifelse(n_from <= n_to, 0, 1) + aon_config <- switch(aon_method, + d = list(aon = grouped_aon, contract = FALSE, phast = FALSE), + bi = list(aon = 2L, contract = FALSE, phast = FALSE), + nba = list(aon = 3L, contract = FALSE, phast = FALSE), + cbi = list(aon = 2L, contract = TRUE, phast = FALSE), + cphast = list(aon = grouped_aon, contract = TRUE, phast = TRUE), + cch = list(aon = 7L, contract = FALSE, phast = FALSE)) + aon <- aon_config$aon + contract <- aon_config$contract + phast <- aon_config$phast vec <- rep(0, nrow(Graph$data)) @@ -188,9 +192,22 @@ assign_traffic <- function(Graph, from, to, demand, algorithm = "bfw", max_gap = if (algorithm %in% c("msa", "fw", "cfw", "bfw")){ algo <- which(c("msa", "fw", "cfw", "bfw") == algorithm) - 1 - res <- cpptraffic(Graph$data$from, Graph$data$to, Graph$data$dist, vec, vec, Graph$data$dist, - Graph$attrib$alpha, Graph$attrib$beta, Graph$attrib$cap, Graph$nbnode, - lat, lon, constant, from_id, to_id, demand, max_gap, max_it, algo, aon, contract, phast, verbose) + if (aon_method == "cch"){ + if (is.null(cch)) cch <- cpp_contract(Graph, customizable = TRUE) + if (!inherits(cch, "cppRouting_cch")) stop("cch must be generated by cpp_contract() with customizable = TRUE") + res <- cpptrafficcch(Graph$data$from, Graph$data$to, Graph$data$dist, vec, vec, Graph$data$dist, + Graph$attrib$alpha, Graph$attrib$beta, Graph$attrib$cap, Graph$nbnode, + from_id, to_id, demand, max_gap, max_it, algo, + cch$rank, cch$tail, cch$head, + cch$first_out, cch$adj_head, cch$adj_arc, + cch$rank_first_out, cch$rank_adj_head, cch$rank_adj_arc, + cch$input_arc, cch$input_forward, cch$elimination_tree_parent, + verbose) + } else { + res <- cpptraffic(Graph$data$from, Graph$data$to, Graph$data$dist, vec, vec, Graph$data$dist, + Graph$attrib$alpha, Graph$attrib$beta, Graph$attrib$cap, Graph$nbnode, + lat, lon, constant, from_id, to_id, demand, max_gap, max_it, algo, aon, contract, phast, verbose) + } } diff --git a/R/aux_helpers.R b/R/aux_helpers.R new file mode 100644 index 0000000..4c457f8 --- /dev/null +++ b/R/aux_helpers.R @@ -0,0 +1,43 @@ +# Internal helpers shared by get_distance_pair(), get_distance_matrix() (and, +# eventually, other query functions) to support aggregate_aux with either a +# single legacy aux vector or a multi-column aux data.frame declared at +# makegraph() time. + +# Turn (aux attribute, aggregate_aux argument) into a numeric matrix ready for +# the C++ layer, plus enough metadata to shape the result the same way the +# request was expressed. +resolve_aggregate_aux <- function(aux, aggregate_aux){ + if (is.null(aux)) stop("network don't have additional weight") + aux_is_df <- is.data.frame(aux) + + if (isTRUE(aggregate_aux)) { + cols <- if (aux_is_df) colnames(aux) else "aux" + } else { + cols <- as.character(aggregate_aux) + if (length(cols) == 0) stop("aggregate_aux must be TRUE or a character vector of aux column names") + } + + if (aux_is_df) { + missing <- setdiff(cols, colnames(aux)) + if (length(missing) > 0) stop(sprintf("aux column(s) not found: %s", paste(missing, collapse = ", "))) + mat <- as.matrix(aux[, cols, drop = FALSE]) + } else { + unknown <- setdiff(cols, "aux") + if (length(unknown) > 0) stop(sprintf("aux column(s) not found: %s", paste(unknown, collapse = ", "))) + mat <- matrix(as.double(aux), ncol = 1, dimnames = list(NULL, "aux")) + } + storage.mode(mat) <- "double" + + list(matrix = mat, names = cols, single_legacy = isTRUE(aggregate_aux) && !aux_is_df) +} + +# Reshape a (n_pairs x n_cols) numeric matrix of aggregated values back into a +# plain vector for the single-column legacy case (aggregate_aux = TRUE on a +# graph with one aux column, matching the historical return shape), or a +# data.frame with named columns otherwise. +shape_aggregate_result <- function(value_matrix, col_names, single_legacy){ + if (single_legacy) return(as.numeric(value_matrix[, 1])) + out <- as.data.frame(value_matrix) + colnames(out) <- col_names + out +} diff --git a/R/cch.R b/R/cch.R new file mode 100644 index 0000000..cab2954 --- /dev/null +++ b/R/cch.R @@ -0,0 +1,83 @@ +# CCH objects are deliberately kept as plain R lists. They can be saved with +# saveRDS() and reused later, which is important for traffic assignment models +# that run many congestion iterations on the same road topology. + +cch_topology_fields <- function(){ + c("rank", "tail", "head", "first_out", "adj_head", "adj_arc", + "rank_first_out", "rank_adj_head", "rank_adj_arc", + "input_arc", "input_forward", "elimination_tree_parent") +} + +validate_cch_for_graph <- function(cch, Graph){ + if (!inherits(cch, "cppRouting_cch")) stop("cch must be generated by cpp_contract() with customizable = TRUE") + if (cch$nbnode != Graph$nbnode) stop("cch was prepared for a graph with a different number of nodes") + if (nrow(cch$original$data) != nrow(Graph$data)) stop("cch was prepared for a graph with a different number of edges") + if (!all(cch$original$data$from == Graph$data$from) || !all(cch$original$data$to == Graph$data$to)) { + stop("cch must be prepared from the same graph topology") + } + + missing_fields <- setdiff(cch_topology_fields(), names(cch)) + if (length(missing_fields) > 0) { + stop("cch was prepared with an older CCH format. Regenerate it with cpp_contract(customizable = TRUE).") + } + + TRUE +} + +#' Customize a CCH with edge weights +#' +#' @description +#' Update a prepared CCH with the edge weights that should be used for shortest +#' path queries. Customization is much cheaper than preparing the CCH topology, +#' so it is the step to repeat when congestion changes travel times. +#' +#' @param CCH An object generated by \link{cpp_contract} with \code{customizable = TRUE}. +#' @param weights Optional edge weights. Defaults to the original graph cost. +#' @return A customized CCH metric object usable by \link{get_distance_pair}, +#' \link{get_distance_matrix} and \link{get_aon}. +#' @details +#' The \code{weights} vector must follow the same edge order as +#' \code{CCH$original$data}. In a traffic model this is normally the current +#' congested travel time for every edge. +#' @examples +#' RcppParallel::setThreadOptions(numThreads = 1) +#' +#' edges <- data.frame(from_vertex = c(0, 0, 1, 1, 2), +#' to_vertex = c(1, 2, 2, 3, 3), +#' cost = c(1, 4, 1, 2, 1)) +#' graph <- makegraph(edges, directed = TRUE) +#' +#' cch <- cpp_contract(graph, customizable = TRUE) +#' metric <- cpp_customize(cch) +#' get_distance_pair(metric, from = c(0, 1), to = c(3, 3)) +#' @export +cpp_customize <- function(CCH, weights = NULL){ + if (!inherits(CCH, "cppRouting_cch")) stop("Input should be generated by cpp_contract() with customizable = TRUE") + if (length(setdiff(cch_topology_fields(), names(CCH))) > 0) { + stop("CCH was prepared with an older CCH format. Regenerate it with cpp_contract(customizable = TRUE).") + } + + if (is.null(weights)) weights <- CCH$original$data$dist + weights <- as.numeric(weights) + if (length(weights) != nrow(CCH$original$data)) stop("weights must have one value per original edge") + if (any(is.na(weights))) stop("NAs are not allowed in weights") + if (any(weights < 0)) stop("Negative weights are not allowed") + + metric <- cppcchcustomize(CCH$original$data$from, + CCH$original$data$to, + weights, + CCH$nbnode, + CCH$rank, + CCH$tail, + CCH$head, + CCH$rank_first_out, + CCH$rank_adj_head, + CCH$rank_adj_arc, + CCH$input_arc, + CCH$input_forward) + + CCH$weight <- metric + CCH$weights <- weights + structure(CCH, + class = "cppRouting_cch_metric") +} diff --git a/R/contract.R b/R/contract.R index c798cf3..a9936ad 100644 --- a/R/contract.R +++ b/R/contract.R @@ -1,10 +1,32 @@ #' Contraction hierarchies algorithm -#' @description Contract a graph by using contraction hierarchies algorithm +#' @description Contract a graph by using contraction hierarchies algorithm. #' @param Graph An object generated by \link{makegraph} or \link{cpp_simplify} function. #' @param silent Logical. If \code{TRUE}, progress is not displayed. -#' @return A contracted graph. +#' @param customizable Logical. If \code{TRUE}, build a customizable contraction +#' hierarchy (CCH) instead of a classical one: only the graph topology is contracted, so +#' shortcut costs can be recomputed cheaply with \link{cpp_customize} whenever edge weights +#' change. See Details. +#' @param order Optional node order, only used when \code{customizable = TRUE}. Values must +#' be node references from \code{Graph$dict$ref}. If omitted, cppRouting computes a simple +#' local order. Large road networks usually benefit from a stronger nested-dissection order +#' prepared outside this function. +#' @return A contracted graph, or (when \code{customizable = TRUE}) a CCH topology object. +#' Run \link{cpp_customize} before distance queries on a CCH topology object, or pass it to +#' \link{assign_traffic} through the \code{cch} argument to avoid rebuilding the topology at +#' each model run. #' @details Contraction hierarchies is a speed-up technique for finding shortest path in a graph. #' It consist of two steps : preprocessing phase and query. \code{cpp_contract()} preprocess the input graph to later use special query algorithm implemented in \link{get_distance_pair}, \link{get_distance_matrix}, \link{get_aon} and \link{get_path_pair} functions. +#' +#' When \code{customizable = TRUE}, contraction is split into two phases instead: \itemize{ +#' \item \strong{Preparation} (this function): builds shortcuts from the road topology. This +#' step is reusable as long as the set of edges does not change. +#' \item \strong{Customization} (\link{cpp_customize}): updates shortcut costs from the +#' current edge weights. This step is repeated when congestion changes travel times. +#' } +#' For traffic assignment, prepare once and call \code{assign_traffic(..., aon_method = "cch", +#' cch = cch)}. The CCH assignment path uses the elimination-tree query, which is the fast +#' query variant for repeated origins or repeated destinations. +#' #' To see the benefits of using contraction hierarchies, see the package description : \url{https://github.com/vlarmet/cppRouting/blob/master/README.md}. #' @seealso \link{cpp_simplify} #' @examples @@ -18,23 +40,48 @@ #' #' #Contract graph #' contracted_graph<-cpp_contract(graph,silent=TRUE) +#' +#' #Build a customizable contraction hierarchy instead +#' cch <- cpp_contract(graph, customizable = TRUE) +#' metric <- cpp_customize(cch) +#' get_distance_pair(metric, from = 0, to = 3) + + +cpp_contract<-function(Graph,silent=FALSE,customizable=FALSE,order=NULL){ + if (!inherits(Graph, "cppRouting_graph")) stop("Input should be generated by makegraph() or cpp_simplify() function") + + if (customizable){ + if (is.null(order)){ + order_id <- integer(0) + } else { + order <- as.character(order) + if (length(order) != Graph$nbnode) stop("order must contain one entry for each node") + if (sum(order %in% Graph$dict$ref) < length(order)) stop("Some order nodes are not in the graph") + if (any(duplicated(order))) stop("order contains duplicated nodes") + order_id <- Graph$dict$id[match(order, Graph$dict$ref)] + } + cch <- cppcchprepare(Graph$data$from, Graph$data$to, Graph$nbnode, order_id, !silent) -cpp_contract<-function(Graph,silent=FALSE){ - if (length(Graph) != 5) stop("Input should be generated by makegraph() or cpp_simplify() function") + cch$nbnode <- Graph$nbnode + cch$dict <- Graph$dict + cch$original <- Graph[c("data", "attrib")] + return(structure(cch, class = "cppRouting_cch")) + } res <- cppcontract(Graph$data$from,Graph$data$to,Graph$data$dist,Graph$nbnode,!silent) - return(list(data=data.frame(from=res[[1]][[1]], - to=res[[1]][[2]], - dist=res[[1]][[3]]), - rank=res[[2]], - shortcuts=data.frame(shortf=res[[3]][[1]], - shortt=res[[3]][[2]], - shortc=res[[3]][[3]]), - nbnode=Graph$nbnode, - dict=Graph$dict, - original = Graph[c("data", "attrib")])) + structure(list(data=data.frame(from=res[[1]][[1]], + to=res[[1]][[2]], + dist=res[[1]][[3]]), + rank=res[[2]], + shortcuts=data.frame(shortf=res[[3]][[1]], + shortt=res[[3]][[2]], + shortc=res[[3]][[3]]), + nbnode=Graph$nbnode, + dict=Graph$dict, + original = Graph[c("data", "attrib")]), + class = "cppRouting_contracted") } diff --git a/R/get_aon.R b/R/get_aon.R index a4789f3..9b0945b 100644 --- a/R/get_aon.R +++ b/R/get_aon.R @@ -60,7 +60,40 @@ get_aon<-function(Graph, from, to, demand, algorithm = "bi",constant = 1){ from_id <- Graph$dict$id[match(from,Graph$dict$ref)] to_id <- Graph$dict$id[match(to,Graph$dict$ref)] - if (length(Graph) == 5){ + if (inherits(Graph, "cppRouting_cch")) stop("Run cpp_customize() before querying a CCH") + + if (inherits(Graph, "cppRouting_cch_metric")){ + res <- cppaoncchelim(Graph$original$data$from, + Graph$original$data$to, + Graph$weights, + Graph$nbnode, + Graph$rank, + Graph$first_out, + Graph$adj_head, + Graph$adj_arc, + Graph$elimination_tree_parent, + Graph$weight$forward, + Graph$weight$backward, + Graph$weight$forward_first_arc, + Graph$weight$forward_first_dir, + Graph$weight$forward_second_arc, + Graph$weight$forward_second_dir, + Graph$weight$forward_original, + Graph$weight$backward_first_arc, + Graph$weight$backward_first_dir, + Graph$weight$backward_second_arc, + Graph$weight$backward_second_dir, + Graph$weight$backward_original, + from_id, + to_id, + demand) + res <- as.data.frame(res, col.names = c("from", "to", "cost", "flow")) + res$from <- Graph$dict$ref[match(res$from, Graph$dict$id)] + res$to <- Graph$dict$ref[match(res$to, Graph$dict$id)] + return(res) + } + + if (inherits(Graph, "cppRouting_graph")){ if (!algorithm %in% c("bi", "d", "nba")) stop("algorithm argument should be 'd', 'bi' or 'nba'") if (algorithm == "nba" & is.null(Graph$coords)) { algorithm <- "bi" @@ -84,7 +117,7 @@ get_aon<-function(Graph, from, to, demand, algorithm = "bi",constant = 1){ } } - if (length(Graph) == 6){ + if (inherits(Graph, "cppRouting_contracted")){ if (!algorithm %in% c("bi", "phast")) stop("algorithm argument should be 'phast' or 'bi'") if (algorithm == "bi"){ res <- cppaonC(Graph$original$data$from, Graph$original$data$to, Graph$original$data$dist, diff --git a/R/get_detour.R b/R/get_detour.R index 350b6c8..f778d2f 100644 --- a/R/get_detour.R +++ b/R/get_detour.R @@ -61,7 +61,7 @@ #' } get_detour<-function(Graph,from,to,extra=NULL,keep=NULL,long=FALSE){ - if (length(Graph) != 5) stop("Input should be generated by makegraph() or cpp_simplify() function") + if (!inherits(Graph, "cppRouting_graph")) stop("Input should be generated by makegraph() or cpp_simplify() function") if (length(from)!=length(to)) stop("From and to have not the same length") if (is.null(extra)) stop("No additional cost provided") if (length(extra)!=1) stop("extra should be one number") diff --git a/R/get_distance_mat.R b/R/get_distance_mat.R index 69c68a4..f5bd2d3 100644 --- a/R/get_distance_mat.R +++ b/R/get_distance_mat.R @@ -4,10 +4,14 @@ #' @param from A vector of one or more vertices from which distances are calculated (origin). #' @param to A vector of one or more vertices (destination). #' @param algorithm Character. Only for contracted graph, \code{mch} for Many to many CH, \code{phast} for PHAST algorithm -#' @param aggregate_aux Logical. If \code{TRUE}, the additional weight is summed along shortest paths. +#' @param aggregate_aux Logical or character vector. If \code{TRUE}, the additional weight(s) +#' declared in \link{makegraph}'s \code{aux} are summed along shortest paths (all columns, if +#' \code{aux} was a data.frame). A character vector of column names restricts this to specific +#' columns. #' @param allcores Logical (deprecated). If \code{TRUE}, all cores are used. -#' @return Matrix of shortest distances. -#' @note It is not possible to aggregate auxiliary weights on a Graph object coming from \link{cpp_simplify} function. +#' @return Matrix of shortest distances, or of aggregated \code{aux} values when +#' \code{aggregate_aux} is used. Requesting more than one \code{aux} column returns +#' a named list of matrices (one per requested value) instead of a single matrix. #' @details If graph is not contracted, \code{get_distance_matrix()} recursively perform Dijkstra algorithm for each \code{from} nodes. #' If graph is contracted, the user has the choice between : \itemize{ #' \item many to many contraction hierarchies (\code{mch}) : optimal for square matrix. @@ -17,6 +21,12 @@ #' Shortest path is always computed according to the main edge weights, corresponding to the 3rd column of \code{df} argument in \code{makegraph()} function. #' If \code{aggregate_aux} argument is \code{TRUE}, the values returned are the sum of auxiliary weights along shortest paths. #' +#' When several shortest paths tie on the main edge weight, \code{aggregate_aux} is summed along +#' whichever tied path the search happens to settle on. Normal, \link{cpp_contract}-contracted and +#' CCH graphs can break such ties differently, so the aggregated value for a given \code{from}/\code{to} +#' pair may occasionally differ between them, even though the main shortest distance always matches +#' exactly. This is expected: it reflects a genuine choice between equally short paths, not an error. +#' #' All algorithms are \strong{multithreaded.} \code{allcores} argument is deprecated, please use \code{RcppParallel::setThreadOptions()} to set the number of threads. #' #' See details in package website : \url{https://github.com/vlarmet/cppRouting/blob/master/README.md} @@ -59,9 +69,152 @@ get_distance_matrix<-function(Graph, from, to, algorithm="phast", aggregate_aux if (allcores) message("allcores argument is deprecated since v3.0. \nPlease use RcppParallel::setThreadOptions() to set the number of threads") + if (inherits(Graph, "cppRouting_cch")) stop("Run cpp_customize() before querying a CCH") + + if (inherits(Graph, "cppRouting_cch_metric")){ + if (isFALSE(aggregate_aux)){ + res <- cppdistmatcch(Graph$nbnode, + Graph$first_out, + Graph$adj_head, + Graph$adj_arc, + Graph$weight$forward, + Graph$weight$backward, + from_id, + to_id) + rownames(res)<-from + colnames(res)<-to + return(res) + } + + # aggregate_aux on CCH: expand to the full from x to cross product and + # reuse the pair-based path_values_pair_cch engine (grouped by shared + # source/target, so this stays efficient for a dense matrix), then + # reshape each requested column back into a from x to matrix. + resolved <- resolve_aggregate_aux(Graph$original$attrib$aux, aggregate_aux) + from_id_x <- rep(from_id, times = length(to_id)) + to_id_x <- rep(to_id, each = length(from_id)) + result <- cpppathvaluescch(Graph$original$data$from, + Graph$nbnode, + Graph$rank, + Graph$first_out, + Graph$adj_head, + Graph$adj_arc, + Graph$elimination_tree_parent, + Graph$weight$forward, + Graph$weight$backward, + Graph$weight$forward_first_arc, + Graph$weight$forward_first_dir, + Graph$weight$forward_second_arc, + Graph$weight$forward_second_dir, + Graph$weight$forward_original, + Graph$weight$backward_first_arc, + Graph$weight$backward_first_dir, + Graph$weight$backward_second_arc, + Graph$weight$backward_second_dir, + Graph$weight$backward_original, + from_id_x, + to_id_x, + resolved$matrix) + + mats <- lapply(seq_along(resolved$names), function(col) { + m <- matrix(result[, col + 1], nrow = length(from), ncol = length(to)) + rownames(m) <- from; colnames(m) <- to + m + }) + names(mats) <- resolved$names + if (resolved$single_legacy) return(mats[[1]]) + return(mats) + } + + # Resolve aggregate_aux before the algorithm-dispatch blocks below (same + # strategy as get_distance_pair(): single-column requests are substituted + # into a local copy of Graph so the existing fast inline-summing code below + # runs unchanged; multi-column requests are handled once here, inline in + # C++ (see cppmataddmulti/cppmataddCmulti below), no path reconstruction). + aux_multi <- NULL + single_col_resolved <- NULL + if (!isFALSE(aggregate_aux)) { + if (inherits(Graph, "cppRouting_graph")) { + resolved <- resolve_aggregate_aux(Graph$attrib$aux, aggregate_aux) + if (ncol(resolved$matrix) > 1) { + aux_multi <- resolved + } else { + Graph$attrib$aux <- resolved$matrix[, 1] + aggregate_aux <- TRUE + single_col_resolved <- resolved + } + } else if (inherits(Graph, "cppRouting_contracted")) { + resolved <- resolve_aggregate_aux(Graph$original$attrib$aux, aggregate_aux) + if (ncol(resolved$matrix) > 1) { + aux_multi <- resolved + } else { + Graph$original$attrib$aux <- resolved$matrix[, 1] + aggregate_aux <- TRUE + single_col_resolved <- resolved + } + } + } + + if (!is.null(aux_multi)) { + # Normal graph: aux columns summed inline during the same search (no + # path reconstruction, no R round-trip). Mirrors the scalar dispatch + # below: Dijkstra is rooted at the smaller side of from/to, transposing + # once at the end if `to` was the smaller side. + if (inherits(Graph, "cppRouting_graph")) { + if (length(from) < length(to)) { + arr3d <- cppmataddmulti(Graph$data[,1], Graph$data[,2], Graph$data[,3], aux_multi$matrix, Graph$nbnode, from_id, to_id) + transpose_needed <- FALSE + } else { + arr3d <- cppmataddmulti(Graph$data[,2], Graph$data[,1], Graph$data[,3], aux_multi$matrix, Graph$nbnode, to_id, from_id) + transpose_needed <- TRUE + } + mats <- lapply(seq_along(aux_multi$names), function(col) { + m <- matrix(arr3d[,,col], nrow = dim(arr3d)[1], ncol = dim(arr3d)[2]) + if (transpose_needed) m <- t(m) + rownames(m) <- from; colnames(m) <- to + m + }) + names(mats) <- aux_multi$names + if (aux_multi$single_legacy) return(mats[[1]]) + return(mats) + } + + # Contracted (CH) graph: aux columns summed inline via the same + # per-shortcut precompute pass (aggCAdd) used by get_distance_pair(), + # feeding a PHAST-based sweep (the only CH matrix algorithm that + # supports aux, same constraint as the existing single-column path, + # which is why `algorithm` is ignored here too). + if (inherits(Graph, "cppRouting_contracted")) { + invrank <- Graph$nbnode - Graph$rank + if (length(from) < length(to)) { + arr3d <- cppmataddCmulti(invrank[Graph$original$data[,1] + 1], invrank[Graph$original$data[,2] + 1], Graph$original$data[,3], aux_multi$matrix, + invrank[Graph$data[,1] + 1], invrank[Graph$data[,2] + 1], Graph$data[,3], Graph$nbnode, Graph$rank, + invrank[Graph$shortcuts$shortf + 1], invrank[Graph$shortcuts$shortt + 1], invrank[Graph$shortcuts$shortc + 1], + TRUE, + invrank[from_id + 1], invrank[to_id + 1]) + transpose_needed <- FALSE + } else { + arr3d <- cppmataddCmulti(invrank[Graph$original$data[,2] + 1], invrank[Graph$original$data[,1] + 1], Graph$original$data[,3], aux_multi$matrix, + invrank[Graph$data[,2] + 1], invrank[Graph$data[,1] + 1], Graph$data[,3], Graph$nbnode, Graph$rank, + invrank[Graph$shortcuts$shortt + 1], invrank[Graph$shortcuts$shortf + 1], invrank[Graph$shortcuts$shortc + 1], + TRUE, + invrank[to_id + 1], invrank[from_id + 1]) + transpose_needed <- TRUE + } + mats <- lapply(seq_along(aux_multi$names), function(col) { + m <- matrix(arr3d[,,col], nrow = dim(arr3d)[1], ncol = dim(arr3d)[2]) + if (transpose_needed) m <- t(m) + rownames(m) <- from; colnames(m) <- to + m + }) + names(mats) <- aux_multi$names + if (aux_multi$single_legacy) return(mats[[1]]) + return(mats) + } + } # not contracted - if (length(Graph) == 5){ + if (inherits(Graph, "cppRouting_graph")){ if (aggregate_aux & length(Graph$attrib$aux) == 0) stop("network don't have additional weight") if (length(from)< length(to)){ @@ -86,7 +239,7 @@ get_distance_matrix<-function(Graph, from, to, algorithm="phast", aggregate_aux } # contracted - if (length(Graph) == 6){ + if (inherits(Graph, "cppRouting_contracted")){ if (aggregate_aux & length(Graph$original$attrib$aux) == 0) stop("network don't have additional weight") invrank <- Graph$nbnode - Graph$rank # invrank2 <- data.frame(invrank = invrank, id = 0:(length(invrank)-1)) @@ -142,5 +295,12 @@ get_distance_matrix<-function(Graph, from, to, algorithm="phast", aggregate_aux rownames(res)<-from colnames(res)<-to + + if (!is.null(single_col_resolved) && !single_col_resolved$single_legacy) { + mats <- list(res) + names(mats) <- single_col_resolved$names + return(mats) + } + return(res) } diff --git a/R/get_distance_pair.R b/R/get_distance_pair.R index 51417a5..8e36ee6 100644 --- a/R/get_distance_pair.R +++ b/R/get_distance_pair.R @@ -3,14 +3,18 @@ #' @param Graph An object generated by \link{makegraph}, \link{cpp_simplify} or \link{cpp_contract} function. #' @param from A vector of one or more vertices from which distances are calculated (origin). #' @param to A vector of one or more vertices (destination). -#' @param aggregate_aux Logical. If \code{TRUE}, the additional weight is summed along shortest paths. +#' @param aggregate_aux Logical or character vector. If \code{TRUE}, the additional weight(s) +#' declared in \link{makegraph}'s \code{aux} are summed along shortest paths (all columns, if +#' \code{aux} was a data.frame). A character vector of column names restricts this to specific +#' columns. #' @param algorithm character. \code{Dijkstra} for uni-directional Dijkstra, \code{bi} for bi-directional Dijkstra, \code{A*} for A star unidirectional search or \code{NBA} for New bi-directional A star .Default to \code{bi} #' @param constant numeric. Constant to maintain the heuristic function admissible in \code{A*} and \code{NBA} algorithms. #' Default to 1, when cost is expressed in the same unit than coordinates. See details #' @param allcores Logical (deprecated). If \code{TRUE}, all cores are used. -#' @return Vector of shortest distances. +#' @return Vector of shortest distances, or of aggregated \code{aux} values when +#' \code{aggregate_aux} is used. Requesting more than one \code{aux} column returns +#' a data.frame with one column per requested value instead of a vector. #' @note \code{from} and \code{to} must be the same length. -#' It is not possible to aggregate auxiliary weights on a Graph object coming from \link{cpp_simplify} function. #' #' #' @details If graph is not contracted, the user has the choice between : \itemize{ @@ -24,6 +28,12 @@ #' Shortest path is always computed according to the main edge weights, corresponding to the 3rd column of \code{df} argument in \link{makegraph} function. #' If \code{aggregate_aux} argument is \code{TRUE}, the values returned are the sum of auxiliary weights along shortest paths. #' +#' When several shortest paths tie on the main edge weight, \code{aggregate_aux} is summed along +#' whichever tied path the search happens to settle on. Normal, \link{cpp_contract}-contracted and +#' CCH graphs can break such ties differently, so the aggregated value for a given \code{from}/\code{to} +#' pair may occasionally differ between them, even though the main shortest distance always matches +#' exactly. This is expected: it reflects a genuine choice between equally short paths, not an error. +#' #' #' In A* and New Bidirectional A star algorithms, euclidean distance is used as heuristic function. #' @@ -73,8 +83,116 @@ get_distance_pair<-function(Graph, from, to, aggregate_aux = FALSE, algorithm="b from_id <- Graph$dict$id[match(from,Graph$dict$ref)] to_id <- Graph$dict$id[match(to,Graph$dict$ref)] + if (inherits(Graph, "cppRouting_cch")) stop("Run cpp_customize() before querying a CCH") + + if (inherits(Graph, "cppRouting_cch_metric")){ + if (isFALSE(aggregate_aux)){ + res <- cppdistcch(Graph$nbnode, + Graph$rank, + Graph$first_out, + Graph$adj_head, + Graph$adj_arc, + Graph$elimination_tree_parent, + Graph$weight$forward, + Graph$weight$backward, + from_id, + to_id) + return(res) + } + + # aggregate_aux on CCH recovers the shortest path and sums the requested + # aux column(s) over its original edges (path_values_pair_cch), instead of + # the distance-only elimination-tree query used above. + resolved <- resolve_aggregate_aux(Graph$original$attrib$aux, aggregate_aux) + result <- cpppathvaluescch(Graph$original$data$from, + Graph$nbnode, + Graph$rank, + Graph$first_out, + Graph$adj_head, + Graph$adj_arc, + Graph$elimination_tree_parent, + Graph$weight$forward, + Graph$weight$backward, + Graph$weight$forward_first_arc, + Graph$weight$forward_first_dir, + Graph$weight$forward_second_arc, + Graph$weight$forward_second_dir, + Graph$weight$forward_original, + Graph$weight$backward_first_arc, + Graph$weight$backward_first_dir, + Graph$weight$backward_second_arc, + Graph$weight$backward_second_dir, + Graph$weight$backward_original, + from_id, + to_id, + resolved$matrix) + return(shape_aggregate_result(result[, -1, drop = FALSE], resolved$names, resolved$single_legacy)) + } + + # Resolve aggregate_aux before the algorithm-dispatch blocks below, so they + # never need to know about multi-column requests. Single-column requests + # (the historical case) are turned back into a plain TRUE + a one-column + # aux vector substituted into a local copy of Graph, so the existing fast + # inline-summing code path runs completely unchanged. Multi-column requests + # are handled once here, inline in C++ (see cppdistaddmulti/cppdistaddCmulti + # below), no path reconstruction needed. + aux_multi <- NULL + single_col_resolved <- NULL + if (!isFALSE(aggregate_aux)) { + if (inherits(Graph, "cppRouting_graph")) { + resolved <- resolve_aggregate_aux(Graph$attrib$aux, aggregate_aux) + if (ncol(resolved$matrix) > 1) { + aux_multi <- resolved + } else { + Graph$attrib$aux <- resolved$matrix[, 1] + aggregate_aux <- TRUE + single_col_resolved <- resolved + } + } else if (inherits(Graph, "cppRouting_contracted")) { + resolved <- resolve_aggregate_aux(Graph$original$attrib$aux, aggregate_aux) + if (ncol(resolved$matrix) > 1) { + aux_multi <- resolved + } else { + Graph$original$attrib$aux <- resolved$matrix[, 1] + aggregate_aux <- TRUE + single_col_resolved <- resolved + } + } + } + + if (!is.null(aux_multi)) { + # Normal graph: aux columns are summed inline during the same search (no + # path reconstruction, no R round-trip), for all 4 algorithms. + if (inherits(Graph, "cppRouting_graph")) { + has_coords <- !is.null(Graph$coords) + algo_int <- if (algorithm == "bi") 1L + else if (has_coords && algorithm == "A*") 2L + else if (has_coords && algorithm == "NBA") 3L + else 0L + if (algo_int %in% c(2L, 3L) && constant == 1) warning("Are you sure constant is equal to 1 ?") + coords_lat <- if (has_coords) Graph$coords[,2] else c(0,0) + coords_lon <- if (has_coords) Graph$coords[,3] else c(0,0) + value_matrix <- cppdistaddmulti(Graph$data$from, Graph$data$to, Graph$data$dist, + aux_multi$matrix, Graph$nbnode, + coords_lat, coords_lon, constant, + from_id, to_id, algo_int) + return(shape_aggregate_result(value_matrix, aux_multi$names, aux_multi$single_legacy)) + } - if (length(Graph) == 5) { + # Contracted (CH) graph: aux columns are summed inline too, via a + # per-shortcut precompute pass (aggCAdd) so the query itself never needs + # to unpack shortcuts. + if (inherits(Graph, "cppRouting_contracted")) { + value_matrix <- cppdistaddCmulti(Graph$original$data$from, Graph$original$data$to, Graph$original$data$dist, aux_multi$matrix, + Graph$data$from, Graph$data$to, Graph$data$dist, Graph$nbnode, Graph$rank, + Graph$shortcuts$shortf, Graph$shortcuts$shortt, Graph$shortcuts$shortc, FALSE, + from_id, to_id) + return(shape_aggregate_result(value_matrix, aux_multi$names, aux_multi$single_legacy)) + } + + } + + if (inherits(Graph, "cppRouting_graph")) { if (aggregate_aux & length(Graph$attrib$aux) == 0) stop("network don't have additional weight") if (!is.null(Graph$coords)){ if (algorithm %in% c("NBA","A*","bi")){ @@ -161,7 +279,7 @@ get_distance_pair<-function(Graph, from, to, aggregate_aux = FALSE, algorithm="b } - if (length(Graph) == 6){ + if (inherits(Graph, "cppRouting_contracted")){ if (aggregate_aux & length(Graph$original$attrib$aux) == 0) stop("network don't have additional weight") if (aggregate_aux){ @@ -178,6 +296,10 @@ get_distance_pair<-function(Graph, from, to, aggregate_aux = FALSE, algorithm="b } + if (!is.null(single_col_resolved)) { + return(shape_aggregate_result(matrix(res, ncol = 1), single_col_resolved$names, single_col_resolved$single_legacy)) + } + return(res) diff --git a/R/get_isochrone.R b/R/get_isochrone.R index 1c788fb..406160b 100644 --- a/R/get_isochrone.R +++ b/R/get_isochrone.R @@ -35,7 +35,7 @@ get_isochrone<-function(Graph,from,lim,setdif=FALSE,keep=NULL,long=FALSE){ - if (length(Graph) != 5) stop("Input should be generated by makegraph() or cpp_simplify() function") + if (!inherits(Graph, "cppRouting_graph")) stop("Input should be generated by makegraph() or cpp_simplify() function") if (any(is.na(from))) stop("NAs are not allowed in origin nodes") from<-as.character(from) diff --git a/R/get_multi_paths.R b/R/get_multi_paths.R index da6987d..cdd70b7 100644 --- a/R/get_multi_paths.R +++ b/R/get_multi_paths.R @@ -38,7 +38,7 @@ get_multi_paths<-function(Graph, from ,to ,keep=NULL, long=FALSE){ - if (length(Graph) != 5) stop("Input should be generated by makegraph() or cpp_simplify() function") + if (!inherits(Graph, "cppRouting_graph")) stop("Input should be generated by makegraph() or cpp_simplify() function") if (any(is.na(cbind(from,to)))) stop("NAs are not allowed in origin/destination nodes") diff --git a/R/get_path_pair.R b/R/get_path_pair.R index fff2e93..52aad07 100644 --- a/R/get_path_pair.R +++ b/R/get_path_pair.R @@ -19,6 +19,11 @@ #' #' If the input graph has been contracted by \link{cpp_contract} function, the algorithm is a modified bidirectional search. #' +#' When several shortest paths tie on the main edge weight, the exact node sequence returned for a +#' given \code{from}/\code{to} pair can differ between algorithms and graph representations (normal, +#' contracted or CCH), since each may break the tie differently. All of these paths are equally +#' optimal in terms of the main shortest distance. +#' #' In \code{A*} and \code{NBA} algorithms, euclidean distance is used as heuristic function. #' #' All algorithms are \strong{multithreaded.} Please use \code{RcppParallel::setThreadOptions()} to set the number of threads. @@ -75,7 +80,7 @@ get_path_pair<-function(Graph,from,to,algorithm="bi",constant=1,keep=NULL,long=F to_keep<-rep(1,Graph$nbnode) } - if (length(Graph) == 5){ + if (inherits(Graph, "cppRouting_graph")){ if (!is.null(Graph$coords)){ if (algorithm %in% c("NBA","A*","bi")){ if (algorithm=="A*"){ @@ -118,12 +123,21 @@ get_path_pair<-function(Graph,from,to,algorithm="bi",constant=1,keep=NULL,long=F } - if (length(Graph) == 6){ + if (inherits(Graph, "cppRouting_contracted")){ res <- cpppathC(Graph$data$from, Graph$data$to, Graph$data$dist, Graph$nbnode, Graph$rank, Graph$shortcuts$shortf, Graph$shortcuts$shortt, Graph$shortcuts$shortc, FALSE, Graph$dict$ref, to_keep, from_id, to_id, 0) } + if (inherits(Graph, "cppRouting_cch_metric")){ + res <- cpppathcch(Graph$nbnode, Graph$rank, Graph$original$data$from, Graph$original$data$to, Graph$first_out, Graph$adj_head, Graph$adj_arc, + Graph$elimination_tree_parent, + Graph$weight$forward, Graph$weight$backward, + Graph$weight$forward_first_arc, Graph$weight$forward_first_dir, Graph$weight$forward_second_arc, Graph$weight$forward_second_dir, Graph$weight$forward_original, + Graph$weight$backward_first_arc, Graph$weight$backward_first_dir, Graph$weight$backward_second_arc, Graph$weight$backward_second_dir, Graph$weight$backward_original, + Graph$dict$ref, to_keep, from_id, to_id) + } + if (long){ diff --git a/R/graphs.R b/R/graphs.R index ef0d0b8..f480a9c 100644 --- a/R/graphs.R +++ b/R/graphs.R @@ -3,7 +3,8 @@ #' @param df A data.frame or matrix containing 3 columns: from, to, cost. See details. #' @param directed logical. If \code{FALSE}, then all edges are duplicated by inverting 'from' and 'to' nodes. #' @param coords Optional. A data.frame or matrix containing all nodes coordinates. Columns order should be 'node_ID', 'X', 'Y'. -#' @param aux Optional. A vector or a single value describing an additional edge weight. +#' @param aux Optional. A vector, a single value, or a data.frame with one or +#' more columns, describing one or more additional edge weights. #' @param capacity Optional. A vector or a single value describing edge capacity. Used for traffic assignment. #' @param alpha Optional. A vector or a single value describing alpha parameter. Used for traffic assignment. #' @param beta Optional. A vector or a single value describing beta parameter. Used for traffic assignment. @@ -65,7 +66,17 @@ makegraph<-function(df, Nodes=unique(c(df[,1],df[,2])) - if (!is.null(aux)){ + aux_is_df <- is.data.frame(aux) + if (!is.null(aux) && aux_is_df){ + if (ncol(aux) == 0) stop("aux data.frame must have at least one column") + if (nrow(aux) == 1) aux <- aux[rep(1, nrow(df)), , drop = FALSE] + if (nrow(aux) != nrow(df)) stop("nrow(aux) must equal 1 or nrow(df)") + if (is.null(colnames(aux)) || any(colnames(aux) == "")) stop("aux columns must be named") + for (col in colnames(aux)) { + if (!is.numeric(aux[[col]])) stop("aux columns must be numeric") + if (sum(aux[[col]] < 0) > 0) warning(sprintf("aux column '%s' contains negative values, this weight is aggregate-only and cannot be minimized", col)) + } + } else if (!is.null(aux)){ if (length(aux) == 1) aux <- rep(aux, nrow(df)) if (length(aux) != nrow(df)) stop("length(aux) must equal 1 or nrow(df)") if (sum(aux < 0) > 0) warning("aux contains negative values, this weight is aggregate-only and cannot be minimized") @@ -96,7 +107,7 @@ makegraph<-function(df, df2<-df[,c(2,1,3)] colnames(df2)<-colnames(df) df<-rbind(df,df2) - attrib <- lapply(attrib, rep, 2) + attrib <- lapply(attrib, function(x) if (is.data.frame(x)) rbind(x, x) else rep(x, 2)) } @@ -124,10 +135,12 @@ makegraph<-function(df, coords<-coords[match(Nodes,coords[,1]),] - return(list(data=df, - coords=coords, - nbnode=length(Nodes), - dict=dict, - attrib = attrib)) + Graph <- list(data=df, + coords=coords, + nbnode=length(Nodes), + dict=dict, + attrib = attrib) + class(Graph) <- "cppRouting_graph" + return(Graph) } diff --git a/R/simplify.R b/R/simplify.R index c9c35fe..05d89cb 100644 --- a/R/simplify.R +++ b/R/simplify.R @@ -7,7 +7,10 @@ #' @param silent Logical. If \code{TRUE} and iterate set to \code{TRUE}, number of iteration and number of removed nodes are printed to the console. #' @return The simplified cppRouting graph #' @details To understand why process can be iterated, see the package description : \url{https://github.com/vlarmet/cppRouting/blob/master/README.md} -#' @note Additional edge attributes like \code{aux}, \code{alpha}, \code{beta} and \code{capacity} will be removed. +#' @note \code{aux} is preserved through simplification: when a chain of edges is +#' collapsed into one, its \code{aux} value(s) are summed over the collapsed +#' edges, the same way \code{aggregate_aux} sums them over a shortest path. +#' \code{alpha}, \code{beta} and \code{capacity} are still removed. #' The first iteration usually eliminates the majority of non-intersection nodes and is therefore faster. #' @examples #' #Simple directed graph @@ -58,11 +61,33 @@ cpp_simplify<-function(Graph,keep=NULL,rm_loop=TRUE,iterate=FALSE,silent=TRUE){ to_keep[Graph$dict$ref %in% keep]<-1 } + # aux is carried through simplification (summed over collapsed edges), same + # as get_distance_pair(aggregate_aux=TRUE) sums it over a path. Accept both + # today's single vector and a data.frame of several columns. + aux_in <- Graph$attrib$aux + aux_was_df <- is.data.frame(aux_in) + if (is.null(aux_in)) { + aux_mat <- matrix(numeric(0), nrow = nrow(Graph$data), ncol = 0) + } else if (aux_was_df) { + aux_mat <- as.matrix(aux_in) + storage.mode(aux_mat) <- "double" + } else { + aux_mat <- matrix(as.double(aux_in), ncol = 1, dimnames = list(NULL, "aux")) + } + + simp_raw <- cppsimplify(Graph$data$from,Graph$data$to,Graph$data$dist,Graph$nbnode, + keep = to_keep, rm_loop=rm_loop, iterate, !silent, aux = aux_mat) + simp <- data.frame(from = simp_raw[[1]], to = simp_raw[[2]], dist = simp_raw[[3]]) - simp <- cppsimplify(Graph$data$from,Graph$data$to,Graph$data$dist,Graph$nbnode, - keep = to_keep, rm_loop=rm_loop, iterate, !silent) - simp <- data.frame(simp) - colnames(simp) <- c("from","to","dist") + aux_out <- simp_raw[[4]] + if (ncol(aux_out) == 0) { + new_aux <- NULL + } else if (aux_was_df) { + new_aux <- as.data.frame(aux_out) + colnames(new_aux) <- colnames(aux_in) + } else { + new_aux <- as.numeric(aux_out[, 1]) + } if (nrow(simp)==0) stop("All nodes have been removed") @@ -90,10 +115,11 @@ cpp_simplify<-function(Graph,keep=NULL,rm_loop=TRUE,iterate=FALSE,silent=TRUE){ - return (list(data=simp, - coords=coords, - nbnode=length(Nodes), - dict=dict, - attrib = list(aux = NULL, alpha = NULL, beta = NULL, cap = NULL))) + return (structure(list(data=simp, + coords=coords, + nbnode=length(Nodes), + dict=dict, + attrib = list(aux = new_aux, alpha = NULL, beta = NULL, cap = NULL)), + class = c("cppRouting_simplified", "cppRouting_graph"))) } diff --git a/R/to_df.R b/R/to_df.R index 58c8b6a..f087437 100644 --- a/R/to_df.R +++ b/R/to_df.R @@ -19,7 +19,7 @@ to_df<-function(Graph){ - if (length(Graph) != 5) stop("Invalid graph input, must be an object generated by makegraph() or cpp_simplify() function") + if (!inherits(Graph, "cppRouting_graph")) stop("Invalid graph input, must be an object generated by makegraph() or cpp_simplify() function") dat<-Graph$data dat$from<-Graph$dict$ref[match(dat$from,Graph$dict$id)] diff --git a/man/assign_traffic.Rd b/man/assign_traffic.Rd index ed84bf1..6023300 100644 --- a/man/assign_traffic.Rd +++ b/man/assign_traffic.Rd @@ -13,6 +13,7 @@ assign_traffic( max_gap = 0.001, max_it = .Machine$integer.max, aon_method = "bi", + cch = NULL, constant = 1, dial_params = NULL, verbose = TRUE @@ -33,7 +34,9 @@ assign_traffic( \item{max_it}{Numeric. Maximum number of iterations. Default to \code{.Machine$integer.max}} -\item{aon_method}{Character.\code{d}, \code{bi}, \code{nba}, \code{cphast} or \code{cbi}. Default to \code{bi}. See details.} +\item{aon_method}{Character.\code{d}, \code{bi}, \code{nba}, \code{cphast}, \code{cbi} or \code{cch}. Default to \code{bi}. See details.} + +\item{cch}{Optional CCH topology generated by \code{cpp_contract(Graph, customizable = TRUE)}. Use it with \code{aon_method = "cch"} to prepare the road topology once and reuse it across many congestion runs.} \item{constant}{numeric. Constant to maintain the heuristic function admissible in NBA* algorithm. Default to 1, when cost is expressed in the same unit than coordinates. See details} @@ -93,12 +96,20 @@ AON algorithms are : \itemize{ \item \code{d} : Dijkstra algorithm \item \code{cbi} : contraction hierarchies + bidirectional search \item \code{cphast} : contraction hierarchies + phast algorithm +\item \code{cch} : customizable contraction hierarchies with RoutingKit-style elimination-tree grouped queries. The topology is prepared once, then edge costs are customized at each assignment iteration. } These AON algorithm can be decomposed into two families, depending the sparsity of origin-destination matrix : \itemize{ \item recursive pairwise : \code{bi}, \code{nba} and \code{cbi}. Optimal for high sparsity. One-to-one algorithm is called N times, with N being the length of from. -\item recursive one-to-many : \code{d} and \code{cphast}. Optimal for dense matrix. One-to-many algorithm is called N times, with N being the number of unique from (or to) nodes +\item recursive one-to-many : \code{d}, \code{cphast} and \code{cch}. Optimal for dense or repeated-endpoint matrices. One-to-many algorithm is called N times, with N being the number of unique from (or to) nodes. } +CCH is most useful when the same graph topology is reused with changing edge +costs, for example during congestion assignment. In that case, run +\code{cch <- cpp_contract(graph, customizable = TRUE)} once, save it if needed, and pass +\code{cch = cch} with \code{aon_method = "cch"}. The assignment will +then customize the CCH at each iteration instead of rebuilding a full +contraction hierarchy. + For large instance, it may be appropriate to test different \code{aon_method} for few iterations and choose the fastest one for the final estimation. Hyperparameters for algorithm-b are : \itemize{ diff --git a/man/cpp_contract.Rd b/man/cpp_contract.Rd index a0d3ec7..f1d40be 100644 --- a/man/cpp_contract.Rd +++ b/man/cpp_contract.Rd @@ -4,22 +4,46 @@ \alias{cpp_contract} \title{Contraction hierarchies algorithm} \usage{ -cpp_contract(Graph, silent = FALSE) +cpp_contract(Graph, silent = FALSE, customizable = FALSE, order = NULL) } \arguments{ \item{Graph}{An object generated by \link{makegraph} or \link{cpp_simplify} function.} \item{silent}{Logical. If \code{TRUE}, progress is not displayed.} + +\item{customizable}{Logical. If \code{TRUE}, build a customizable contraction +hierarchy (CCH) instead of a classical one: only the graph topology is contracted, so +shortcut costs can be recomputed cheaply with \link{cpp_customize} whenever edge weights +change. See Details.} + +\item{order}{Optional node order, only used when \code{customizable = TRUE}. Values must +be node references from \code{Graph$dict$ref}. If omitted, cppRouting computes a simple +local order. Large road networks usually benefit from a stronger nested-dissection order +prepared outside this function.} } \value{ -A contracted graph. +A contracted graph, or (when \code{customizable = TRUE}) a CCH topology object. +Run \link{cpp_customize} before distance queries on a CCH topology object, or pass it to +\link{assign_traffic} through the \code{cch} argument to avoid rebuilding the topology at +each model run. } \description{ -Contract a graph by using contraction hierarchies algorithm +Contract a graph by using contraction hierarchies algorithm. } \details{ Contraction hierarchies is a speed-up technique for finding shortest path in a graph. It consist of two steps : preprocessing phase and query. \code{cpp_contract()} preprocess the input graph to later use special query algorithm implemented in \link{get_distance_pair}, \link{get_distance_matrix}, \link{get_aon} and \link{get_path_pair} functions. + +When \code{customizable = TRUE}, contraction is split into two phases instead: \itemize{ + \item \strong{Preparation} (this function): builds shortcuts from the road topology. This + step is reusable as long as the set of edges does not change. + \item \strong{Customization} (\link{cpp_customize}): updates shortcut costs from the + current edge weights. This step is repeated when congestion changes travel times. +} +For traffic assignment, prepare once and call \code{assign_traffic(..., aon_method = "cch", +cch = cch)}. The CCH assignment path uses the elimination-tree query, which is the fast +query variant for repeated origins or repeated destinations. + To see the benefits of using contraction hierarchies, see the package description : \url{https://github.com/vlarmet/cppRouting/blob/master/README.md}. } \examples{ @@ -33,6 +57,11 @@ graph<-makegraph(edges,directed=TRUE) #Contract graph contracted_graph<-cpp_contract(graph,silent=TRUE) + +#Build a customizable contraction hierarchy instead +cch <- cpp_contract(graph, customizable = TRUE) +metric <- cpp_customize(cch) +get_distance_pair(metric, from = 0, to = 3) } \seealso{ \link{cpp_simplify} diff --git a/man/cpp_customize.Rd b/man/cpp_customize.Rd new file mode 100644 index 0000000..d242d42 --- /dev/null +++ b/man/cpp_customize.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/cch.R +\name{cpp_customize} +\alias{cpp_customize} +\title{Customize a CCH with edge weights} +\usage{ +cpp_customize(CCH, weights = NULL) +} +\arguments{ +\item{CCH}{An object generated by \link{cpp_contract} with \code{customizable = TRUE}.} + +\item{weights}{Optional edge weights. Defaults to the original graph cost.} +} +\value{ +A customized CCH metric object usable by \link{get_distance_pair}, +\link{get_distance_matrix} and \link{get_aon}. +} +\description{ +Update a prepared CCH with the edge weights that should be used for shortest +path queries. Customization is much cheaper than preparing the CCH topology, +so it is the step to repeat when congestion changes travel times. +} +\details{ +The \code{weights} vector must follow the same edge order as +\code{CCH$original$data}. In a traffic model this is normally the current +congested travel time for every edge. +} +\examples{ +RcppParallel::setThreadOptions(numThreads = 1) + +edges <- data.frame(from_vertex = c(0, 0, 1, 1, 2), + to_vertex = c(1, 2, 2, 3, 3), + cost = c(1, 4, 1, 2, 1)) +graph <- makegraph(edges, directed = TRUE) + +cch <- cpp_contract(graph, customizable = TRUE) +metric <- cpp_customize(cch) +get_distance_pair(metric, from = c(0, 1), to = c(3, 3)) +} diff --git a/man/cpp_simplify.Rd b/man/cpp_simplify.Rd index 9194d20..4a35d33 100644 --- a/man/cpp_simplify.Rd +++ b/man/cpp_simplify.Rd @@ -33,7 +33,10 @@ Reduce the number of edges by removing non-intersection nodes, duplicated edges To understand why process can be iterated, see the package description : \url{https://github.com/vlarmet/cppRouting/blob/master/README.md} } \note{ -Additional edge attributes like \code{aux}, \code{alpha}, \code{beta} and \code{capacity} will be removed. +\code{aux} is preserved through simplification: when a chain of edges is +collapsed into one, its \code{aux} value(s) are summed over the collapsed +edges, the same way \code{aggregate_aux} sums them over a shortest path. +\code{alpha}, \code{beta} and \code{capacity} are still removed. The first iteration usually eliminates the majority of non-intersection nodes and is therefore faster. } \examples{ diff --git a/man/get_distance_matrix.Rd b/man/get_distance_matrix.Rd index 65c20c1..8fdd2f0 100644 --- a/man/get_distance_matrix.Rd +++ b/man/get_distance_matrix.Rd @@ -22,12 +22,17 @@ get_distance_matrix( \item{algorithm}{Character. Only for contracted graph, \code{mch} for Many to many CH, \code{phast} for PHAST algorithm} -\item{aggregate_aux}{Logical. If \code{TRUE}, the additional weight is summed along shortest paths.} +\item{aggregate_aux}{Logical or character vector. If \code{TRUE}, the additional weight(s) +declared in \link{makegraph}'s \code{aux} are summed along shortest paths (all columns, if +\code{aux} was a data.frame). A character vector of column names restricts this to specific +columns.} \item{allcores}{Logical (deprecated). If \code{TRUE}, all cores are used.} } \value{ -Matrix of shortest distances. +Matrix of shortest distances, or of aggregated \code{aux} values when +\code{aggregate_aux} is used. Requesting more than one \code{aux} column returns +a named list of matrices (one per requested value) instead of a single matrix. } \description{ Compute all shortest distance between origin and destination nodes. @@ -42,13 +47,16 @@ If graph is contracted, the user has the choice between : \itemize{ Shortest path is always computed according to the main edge weights, corresponding to the 3rd column of \code{df} argument in \code{makegraph()} function. If \code{aggregate_aux} argument is \code{TRUE}, the values returned are the sum of auxiliary weights along shortest paths. +When several shortest paths tie on the main edge weight, \code{aggregate_aux} is summed along +whichever tied path the search happens to settle on. Normal, \link{cpp_contract}-contracted and +CCH graphs can break such ties differently, so the aggregated value for a given \code{from}/\code{to} +pair may occasionally differ between them, even though the main shortest distance always matches +exactly. This is expected: it reflects a genuine choice between equally short paths, not an error. + All algorithms are \strong{multithreaded.} \code{allcores} argument is deprecated, please use \code{RcppParallel::setThreadOptions()} to set the number of threads. See details in package website : \url{https://github.com/vlarmet/cppRouting/blob/master/README.md} } -\note{ -It is not possible to aggregate auxiliary weights on a Graph object coming from \link{cpp_simplify} function. -} \examples{ #Choose number of cores used by cppRouting RcppParallel::setThreadOptions(numThreads = 1) diff --git a/man/get_distance_pair.Rd b/man/get_distance_pair.Rd index b3dba76..6b67845 100644 --- a/man/get_distance_pair.Rd +++ b/man/get_distance_pair.Rd @@ -21,7 +21,10 @@ get_distance_pair( \item{to}{A vector of one or more vertices (destination).} -\item{aggregate_aux}{Logical. If \code{TRUE}, the additional weight is summed along shortest paths.} +\item{aggregate_aux}{Logical or character vector. If \code{TRUE}, the additional weight(s) +declared in \link{makegraph}'s \code{aux} are summed along shortest paths (all columns, if +\code{aux} was a data.frame). A character vector of column names restricts this to specific +columns.} \item{algorithm}{character. \code{Dijkstra} for uni-directional Dijkstra, \code{bi} for bi-directional Dijkstra, \code{A*} for A star unidirectional search or \code{NBA} for New bi-directional A star .Default to \code{bi}} @@ -31,7 +34,9 @@ Default to 1, when cost is expressed in the same unit than coordinates. See deta \item{allcores}{Logical (deprecated). If \code{TRUE}, all cores are used.} } \value{ -Vector of shortest distances. +Vector of shortest distances, or of aggregated \code{aux} values when +\code{aggregate_aux} is used. Requesting more than one \code{aux} column returns +a data.frame with one column per requested value instead of a vector. } \description{ Compute shortest distance between origin and destination nodes. @@ -48,6 +53,12 @@ If the input graph has been contracted by \link{cpp_contract} function, the algo Shortest path is always computed according to the main edge weights, corresponding to the 3rd column of \code{df} argument in \link{makegraph} function. If \code{aggregate_aux} argument is \code{TRUE}, the values returned are the sum of auxiliary weights along shortest paths. +When several shortest paths tie on the main edge weight, \code{aggregate_aux} is summed along +whichever tied path the search happens to settle on. Normal, \link{cpp_contract}-contracted and +CCH graphs can break such ties differently, so the aggregated value for a given \code{from}/\code{to} +pair may occasionally differ between them, even though the main shortest distance always matches +exactly. This is expected: it reflects a genuine choice between equally short paths, not an error. + In A* and New Bidirectional A star algorithms, euclidean distance is used as heuristic function. @@ -58,7 +69,6 @@ To understand the importance of constant parameter, see the package description } \note{ \code{from} and \code{to} must be the same length. -It is not possible to aggregate auxiliary weights on a Graph object coming from \link{cpp_simplify} function. } \examples{ #Choose number of cores used by cppRouting diff --git a/man/get_path_pair.Rd b/man/get_path_pair.Rd index 7e70bf9..df39128 100644 --- a/man/get_path_pair.Rd +++ b/man/get_path_pair.Rd @@ -46,6 +46,11 @@ If graph is not contracted, the user has the choice between : \itemize{ If the input graph has been contracted by \link{cpp_contract} function, the algorithm is a modified bidirectional search. +When several shortest paths tie on the main edge weight, the exact node sequence returned for a +given \code{from}/\code{to} pair can differ between algorithms and graph representations (normal, +contracted or CCH), since each may break the tie differently. All of these paths are equally +optimal in terms of the main shortest distance. + In \code{A*} and \code{NBA} algorithms, euclidean distance is used as heuristic function. All algorithms are \strong{multithreaded.} Please use \code{RcppParallel::setThreadOptions()} to set the number of threads. diff --git a/man/makegraph.Rd b/man/makegraph.Rd index 2bc094a..3aa9a37 100644 --- a/man/makegraph.Rd +++ b/man/makegraph.Rd @@ -21,7 +21,8 @@ makegraph( \item{coords}{Optional. A data.frame or matrix containing all nodes coordinates. Columns order should be 'node_ID', 'X', 'Y'.} -\item{aux}{Optional. A vector or a single value describing an additional edge weight.} +\item{aux}{Optional. A vector, a single value, or a data.frame with one or +more columns, describing one or more additional edge weights.} \item{capacity}{Optional. A vector or a single value describing edge capacity. Used for traffic assignment.} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 074df7a..d3fc6bf 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -92,8 +92,8 @@ BEGIN_RCPP END_RCPP } // cppsimplify -Rcpp::List cppsimplify(std::vector& gfrom, std::vector& gto, std::vector& gw, int nb, std::vector keep, bool rm_loop, bool iterate, bool progress); -RcppExport SEXP _cppRouting_cppsimplify(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP nbSEXP, SEXP keepSEXP, SEXP rm_loopSEXP, SEXP iterateSEXP, SEXP progressSEXP) { +Rcpp::List cppsimplify(std::vector& gfrom, std::vector& gto, std::vector& gw, int nb, std::vector keep, bool rm_loop, bool iterate, bool progress, Rcpp::NumericMatrix aux); +RcppExport SEXP _cppRouting_cppsimplify(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP nbSEXP, SEXP keepSEXP, SEXP rm_loopSEXP, SEXP iterateSEXP, SEXP progressSEXP, SEXP auxSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; @@ -105,7 +105,8 @@ BEGIN_RCPP Rcpp::traits::input_parameter< bool >::type rm_loop(rm_loopSEXP); Rcpp::traits::input_parameter< bool >::type iterate(iterateSEXP); Rcpp::traits::input_parameter< bool >::type progress(progressSEXP); - rcpp_result_gen = Rcpp::wrap(cppsimplify(gfrom, gto, gw, nb, keep, rm_loop, iterate, progress)); + Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type aux(auxSEXP); + rcpp_result_gen = Rcpp::wrap(cppsimplify(gfrom, gto, gw, nb, keep, rm_loop, iterate, progress, aux)); return rcpp_result_gen; END_RCPP } @@ -124,6 +125,181 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// cppcchprepare +Rcpp::List cppcchprepare(std::vector& gfrom, std::vector& gto, int NbNodes, std::vector order, bool display_progress); +RcppExport SEXP _cppRouting_cppcchprepare(SEXP gfromSEXP, SEXP gtoSEXP, SEXP NbNodesSEXP, SEXP orderSEXP, SEXP display_progressSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gto(gtoSEXP); + Rcpp::traits::input_parameter< int >::type NbNodes(NbNodesSEXP); + Rcpp::traits::input_parameter< std::vector >::type order(orderSEXP); + Rcpp::traits::input_parameter< bool >::type display_progress(display_progressSEXP); + rcpp_result_gen = Rcpp::wrap(cppcchprepare(gfrom, gto, NbNodes, order, display_progress)); + return rcpp_result_gen; +END_RCPP +} +// cppcchcustomize +Rcpp::List cppcchcustomize(std::vector& gfrom, std::vector& gto, std::vector& gw, int NbNodes, std::vector& rank, std::vector& tail, std::vector& head, std::vector& rank_first_out, std::vector& rank_adj_head, std::vector& rank_adj_arc, std::vector& input_arc, std::vector& input_forward); +RcppExport SEXP _cppRouting_cppcchcustomize(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP NbNodesSEXP, SEXP rankSEXP, SEXP tailSEXP, SEXP headSEXP, SEXP rank_first_outSEXP, SEXP rank_adj_headSEXP, SEXP rank_adj_arcSEXP, SEXP input_arcSEXP, SEXP input_forwardSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gto(gtoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gw(gwSEXP); + Rcpp::traits::input_parameter< int >::type NbNodes(NbNodesSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank(rankSEXP); + Rcpp::traits::input_parameter< std::vector& >::type tail(tailSEXP); + Rcpp::traits::input_parameter< std::vector& >::type head(headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank_first_out(rank_first_outSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank_adj_head(rank_adj_headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank_adj_arc(rank_adj_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type input_arc(input_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type input_forward(input_forwardSEXP); + rcpp_result_gen = Rcpp::wrap(cppcchcustomize(gfrom, gto, gw, NbNodes, rank, tail, head, rank_first_out, rank_adj_head, rank_adj_arc, input_arc, input_forward)); + return rcpp_result_gen; +END_RCPP +} +// cppdistcch +Rcpp::NumericVector cppdistcch(int NbNodes, std::vector& rank, std::vector& first_out, std::vector& adj_head, std::vector& adj_arc, std::vector& elimination_tree_parent, std::vector& forward, std::vector& backward, std::vector dep, std::vector arr); +RcppExport SEXP _cppRouting_cppdistcch(SEXP NbNodesSEXP, SEXP rankSEXP, SEXP first_outSEXP, SEXP adj_headSEXP, SEXP adj_arcSEXP, SEXP elimination_tree_parentSEXP, SEXP forwardSEXP, SEXP backwardSEXP, SEXP depSEXP, SEXP arrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type NbNodes(NbNodesSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank(rankSEXP); + Rcpp::traits::input_parameter< std::vector& >::type first_out(first_outSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_head(adj_headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_arc(adj_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type elimination_tree_parent(elimination_tree_parentSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward(forwardSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward(backwardSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + rcpp_result_gen = Rcpp::wrap(cppdistcch(NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, dep, arr)); + return rcpp_result_gen; +END_RCPP +} +// cppdistmatcch +Rcpp::NumericMatrix cppdistmatcch(int NbNodes, std::vector& first_out, std::vector& adj_head, std::vector& adj_arc, std::vector& forward, std::vector& backward, std::vector dep, std::vector arr); +RcppExport SEXP _cppRouting_cppdistmatcch(SEXP NbNodesSEXP, SEXP first_outSEXP, SEXP adj_headSEXP, SEXP adj_arcSEXP, SEXP forwardSEXP, SEXP backwardSEXP, SEXP depSEXP, SEXP arrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type NbNodes(NbNodesSEXP); + Rcpp::traits::input_parameter< std::vector& >::type first_out(first_outSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_head(adj_headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_arc(adj_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward(forwardSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward(backwardSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + rcpp_result_gen = Rcpp::wrap(cppdistmatcch(NbNodes, first_out, adj_head, adj_arc, forward, backward, dep, arr)); + return rcpp_result_gen; +END_RCPP +} +// cpppathvaluescch +Rcpp::NumericMatrix cpppathvaluescch(std::vector& gfrom, int NbNodes, std::vector& rank, std::vector& first_out, std::vector& adj_head, std::vector& adj_arc, std::vector& elimination_tree_parent, std::vector& forward, std::vector& backward, std::vector& forward_first_arc, std::vector& forward_first_dir, std::vector& forward_second_arc, std::vector& forward_second_dir, std::vector& forward_original, std::vector& backward_first_arc, std::vector& backward_first_dir, std::vector& backward_second_arc, std::vector& backward_second_dir, std::vector& backward_original, std::vector dep, std::vector arr, Rcpp::NumericMatrix values); +RcppExport SEXP _cppRouting_cpppathvaluescch(SEXP gfromSEXP, SEXP NbNodesSEXP, SEXP rankSEXP, SEXP first_outSEXP, SEXP adj_headSEXP, SEXP adj_arcSEXP, SEXP elimination_tree_parentSEXP, SEXP forwardSEXP, SEXP backwardSEXP, SEXP forward_first_arcSEXP, SEXP forward_first_dirSEXP, SEXP forward_second_arcSEXP, SEXP forward_second_dirSEXP, SEXP forward_originalSEXP, SEXP backward_first_arcSEXP, SEXP backward_first_dirSEXP, SEXP backward_second_arcSEXP, SEXP backward_second_dirSEXP, SEXP backward_originalSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP valuesSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< int >::type NbNodes(NbNodesSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank(rankSEXP); + Rcpp::traits::input_parameter< std::vector& >::type first_out(first_outSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_head(adj_headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_arc(adj_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type elimination_tree_parent(elimination_tree_parentSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward(forwardSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward(backwardSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_first_arc(forward_first_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_first_dir(forward_first_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_second_arc(forward_second_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_second_dir(forward_second_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_original(forward_originalSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_first_arc(backward_first_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_first_dir(backward_first_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_second_arc(backward_second_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_second_dir(backward_second_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_original(backward_originalSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type values(valuesSEXP); + rcpp_result_gen = Rcpp::wrap(cpppathvaluescch(gfrom, NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, dep, arr, values)); + return rcpp_result_gen; +END_RCPP +} +// cpppathcch +Rcpp::List cpppathcch(int NbNodes, std::vector& rank, std::vector& gfrom, std::vector& gto, std::vector& first_out, std::vector& adj_head, std::vector& adj_arc, std::vector& elimination_tree_parent, std::vector& forward, std::vector& backward, std::vector& forward_first_arc, std::vector& forward_first_dir, std::vector& forward_second_arc, std::vector& forward_second_dir, std::vector& forward_original, std::vector& backward_first_arc, std::vector& backward_first_dir, std::vector& backward_second_arc, std::vector& backward_second_dir, std::vector& backward_original, std::vector dict, std::vector keep, std::vector dep, std::vector arr); +RcppExport SEXP _cppRouting_cpppathcch(SEXP NbNodesSEXP, SEXP rankSEXP, SEXP gfromSEXP, SEXP gtoSEXP, SEXP first_outSEXP, SEXP adj_headSEXP, SEXP adj_arcSEXP, SEXP elimination_tree_parentSEXP, SEXP forwardSEXP, SEXP backwardSEXP, SEXP forward_first_arcSEXP, SEXP forward_first_dirSEXP, SEXP forward_second_arcSEXP, SEXP forward_second_dirSEXP, SEXP forward_originalSEXP, SEXP backward_first_arcSEXP, SEXP backward_first_dirSEXP, SEXP backward_second_arcSEXP, SEXP backward_second_dirSEXP, SEXP backward_originalSEXP, SEXP dictSEXP, SEXP keepSEXP, SEXP depSEXP, SEXP arrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< int >::type NbNodes(NbNodesSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank(rankSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gto(gtoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type first_out(first_outSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_head(adj_headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_arc(adj_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type elimination_tree_parent(elimination_tree_parentSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward(forwardSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward(backwardSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_first_arc(forward_first_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_first_dir(forward_first_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_second_arc(forward_second_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_second_dir(forward_second_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_original(forward_originalSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_first_arc(backward_first_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_first_dir(backward_first_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_second_arc(backward_second_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_second_dir(backward_second_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_original(backward_originalSEXP); + Rcpp::traits::input_parameter< std::vector >::type dict(dictSEXP); + Rcpp::traits::input_parameter< std::vector >::type keep(keepSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + rcpp_result_gen = Rcpp::wrap(cpppathcch(NbNodes, rank, gfrom, gto, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, dict, keep, dep, arr)); + return rcpp_result_gen; +END_RCPP +} +// cppaoncchelim +Rcpp::List cppaoncchelim(std::vector& gfrom, std::vector& gto, std::vector& gw, int NbNodes, std::vector& rank, std::vector& first_out, std::vector& adj_head, std::vector& adj_arc, std::vector& elimination_tree_parent, std::vector& forward, std::vector& backward, std::vector& forward_first_arc, std::vector& forward_first_dir, std::vector& forward_second_arc, std::vector& forward_second_dir, std::vector& forward_original, std::vector& backward_first_arc, std::vector& backward_first_dir, std::vector& backward_second_arc, std::vector& backward_second_dir, std::vector& backward_original, std::vector dep, std::vector arr, std::vector dem); +RcppExport SEXP _cppRouting_cppaoncchelim(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP NbNodesSEXP, SEXP rankSEXP, SEXP first_outSEXP, SEXP adj_headSEXP, SEXP adj_arcSEXP, SEXP elimination_tree_parentSEXP, SEXP forwardSEXP, SEXP backwardSEXP, SEXP forward_first_arcSEXP, SEXP forward_first_dirSEXP, SEXP forward_second_arcSEXP, SEXP forward_second_dirSEXP, SEXP forward_originalSEXP, SEXP backward_first_arcSEXP, SEXP backward_first_dirSEXP, SEXP backward_second_arcSEXP, SEXP backward_second_dirSEXP, SEXP backward_originalSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP demSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gto(gtoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gw(gwSEXP); + Rcpp::traits::input_parameter< int >::type NbNodes(NbNodesSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank(rankSEXP); + Rcpp::traits::input_parameter< std::vector& >::type first_out(first_outSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_head(adj_headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_arc(adj_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type elimination_tree_parent(elimination_tree_parentSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward(forwardSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward(backwardSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_first_arc(forward_first_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_first_dir(forward_first_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_second_arc(forward_second_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_second_dir(forward_second_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type forward_original(forward_originalSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_first_arc(backward_first_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_first_dir(backward_first_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_second_arc(backward_second_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_second_dir(backward_second_dirSEXP); + Rcpp::traits::input_parameter< std::vector& >::type backward_original(backward_originalSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + Rcpp::traits::input_parameter< std::vector >::type dem(demSEXP); + rcpp_result_gen = Rcpp::wrap(cppaoncchelim(gfrom, gto, gw, NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, dep, arr, dem)); + return rcpp_result_gen; +END_RCPP +} // cppdistC Rcpp::NumericVector cppdistC(std::vector& gfrom, std::vector& gto, std::vector& gw, int nb, std::vector& rank, std::vector& shortf, std::vector& shortt, std::vector& shortc, bool phast, std::vector dep, std::vector arr, int algo); RcppExport SEXP _cppRouting_cppdistC(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP nbSEXP, SEXP rankSEXP, SEXP shortfSEXP, SEXP shorttSEXP, SEXP shortcSEXP, SEXP phastSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP algoSEXP) { @@ -232,6 +408,23 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// cppmataddmulti +Rcpp::NumericVector cppmataddmulti(std::vector& gfrom, std::vector& gto, std::vector& gw, Rcpp::NumericMatrix aux, int nb, std::vector dep, std::vector arr); +RcppExport SEXP _cppRouting_cppmataddmulti(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP auxSEXP, SEXP nbSEXP, SEXP depSEXP, SEXP arrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gto(gtoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gw(gwSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type aux(auxSEXP); + Rcpp::traits::input_parameter< int >::type nb(nbSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + rcpp_result_gen = Rcpp::wrap(cppmataddmulti(gfrom, gto, gw, aux, nb, dep, arr)); + return rcpp_result_gen; +END_RCPP +} // cppdistadd Rcpp::NumericVector cppdistadd(std::vector& gfrom, std::vector& gto, std::vector& gw, std::vector& gadd, int nb, std::vector lat, std::vector lon, double k, std::vector dep, std::vector arr, int algo); RcppExport SEXP _cppRouting_cppdistadd(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP gaddSEXP, SEXP nbSEXP, SEXP latSEXP, SEXP lonSEXP, SEXP kSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP algoSEXP) { @@ -253,6 +446,27 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// cppdistaddmulti +Rcpp::NumericMatrix cppdistaddmulti(std::vector& gfrom, std::vector& gto, std::vector& gw, Rcpp::NumericMatrix aux, int nb, std::vector lat, std::vector lon, double k, std::vector dep, std::vector arr, int algo); +RcppExport SEXP _cppRouting_cppdistaddmulti(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP auxSEXP, SEXP nbSEXP, SEXP latSEXP, SEXP lonSEXP, SEXP kSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP algoSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gto(gtoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gw(gwSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type aux(auxSEXP); + Rcpp::traits::input_parameter< int >::type nb(nbSEXP); + Rcpp::traits::input_parameter< std::vector >::type lat(latSEXP); + Rcpp::traits::input_parameter< std::vector >::type lon(lonSEXP); + Rcpp::traits::input_parameter< double >::type k(kSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + Rcpp::traits::input_parameter< int >::type algo(algoSEXP); + rcpp_result_gen = Rcpp::wrap(cppdistaddmulti(gfrom, gto, gw, aux, nb, lat, lon, k, dep, arr, algo)); + return rcpp_result_gen; +END_RCPP +} // cppdistaddC Rcpp::NumericVector cppdistaddC(std::vector& orfrom, std::vector& orto, std::vector& orw, std::vector& gadd, std::vector& gfrom, std::vector& gto, std::vector& gw, int nb, std::vector& rank, std::vector& shortf, std::vector& shortt, std::vector& shortc, bool phast, std::vector dep, std::vector arr, int algo); RcppExport SEXP _cppRouting_cppdistaddC(SEXP orfromSEXP, SEXP ortoSEXP, SEXP orwSEXP, SEXP gaddSEXP, SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP nbSEXP, SEXP rankSEXP, SEXP shortfSEXP, SEXP shorttSEXP, SEXP shortcSEXP, SEXP phastSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP algoSEXP) { @@ -279,6 +493,31 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// cppdistaddCmulti +Rcpp::NumericMatrix cppdistaddCmulti(std::vector& orfrom, std::vector& orto, std::vector& orw, Rcpp::NumericMatrix aux, std::vector& gfrom, std::vector& gto, std::vector& gw, int nb, std::vector& rank, std::vector& shortf, std::vector& shortt, std::vector& shortc, bool phast, std::vector dep, std::vector arr); +RcppExport SEXP _cppRouting_cppdistaddCmulti(SEXP orfromSEXP, SEXP ortoSEXP, SEXP orwSEXP, SEXP auxSEXP, SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP nbSEXP, SEXP rankSEXP, SEXP shortfSEXP, SEXP shorttSEXP, SEXP shortcSEXP, SEXP phastSEXP, SEXP depSEXP, SEXP arrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::vector& >::type orfrom(orfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type orto(ortoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type orw(orwSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type aux(auxSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gto(gtoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gw(gwSEXP); + Rcpp::traits::input_parameter< int >::type nb(nbSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank(rankSEXP); + Rcpp::traits::input_parameter< std::vector& >::type shortf(shortfSEXP); + Rcpp::traits::input_parameter< std::vector& >::type shortt(shorttSEXP); + Rcpp::traits::input_parameter< std::vector& >::type shortc(shortcSEXP); + Rcpp::traits::input_parameter< bool >::type phast(phastSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + rcpp_result_gen = Rcpp::wrap(cppdistaddCmulti(orfrom, orto, orw, aux, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr)); + return rcpp_result_gen; +END_RCPP +} // cppaddC Rcpp::NumericMatrix cppaddC(std::vector& orfrom, std::vector& orto, std::vector& orw, std::vector& gadd, std::vector& gfrom, std::vector& gto, std::vector& gw, int nb, std::vector& rank, std::vector& shortf, std::vector& shortt, std::vector& shortc, bool phast, std::vector dep, std::vector arr, int algo); RcppExport SEXP _cppRouting_cppaddC(SEXP orfromSEXP, SEXP ortoSEXP, SEXP orwSEXP, SEXP gaddSEXP, SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP nbSEXP, SEXP rankSEXP, SEXP shortfSEXP, SEXP shorttSEXP, SEXP shortcSEXP, SEXP phastSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP algoSEXP) { @@ -305,6 +544,31 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// cppmataddCmulti +Rcpp::NumericVector cppmataddCmulti(std::vector& orfrom, std::vector& orto, std::vector& orw, Rcpp::NumericMatrix aux, std::vector& gfrom, std::vector& gto, std::vector& gw, int nb, std::vector& rank, std::vector& shortf, std::vector& shortt, std::vector& shortc, bool phast, std::vector dep, std::vector arr); +RcppExport SEXP _cppRouting_cppmataddCmulti(SEXP orfromSEXP, SEXP ortoSEXP, SEXP orwSEXP, SEXP auxSEXP, SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP nbSEXP, SEXP rankSEXP, SEXP shortfSEXP, SEXP shorttSEXP, SEXP shortcSEXP, SEXP phastSEXP, SEXP depSEXP, SEXP arrSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::vector& >::type orfrom(orfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type orto(ortoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type orw(orwSEXP); + Rcpp::traits::input_parameter< Rcpp::NumericMatrix >::type aux(auxSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gto(gtoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gw(gwSEXP); + Rcpp::traits::input_parameter< int >::type nb(nbSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank(rankSEXP); + Rcpp::traits::input_parameter< std::vector& >::type shortf(shortfSEXP); + Rcpp::traits::input_parameter< std::vector& >::type shortt(shorttSEXP); + Rcpp::traits::input_parameter< std::vector& >::type shortc(shortcSEXP); + Rcpp::traits::input_parameter< bool >::type phast(phastSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + rcpp_result_gen = Rcpp::wrap(cppmataddCmulti(orfrom, orto, orw, aux, gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast, dep, arr)); + return rcpp_result_gen; +END_RCPP +} // cppaon Rcpp::List cppaon(std::vector& gfrom, std::vector& gto, std::vector& gw, int nb, std::vector lat, std::vector lon, double k, std::vector dep, std::vector arr, std::vector dem, int algo); RcppExport SEXP _cppRouting_cppaon(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP nbSEXP, SEXP latSEXP, SEXP lonSEXP, SEXP kSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP demSEXP, SEXP algoSEXP) { @@ -352,6 +616,45 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } +// cpptrafficcch +Rcpp::List cpptrafficcch(std::vector& gfrom, std::vector& gto, std::vector& gw, std::vector& gflow, std::vector& gaux, std::vector& gftt, std::vector& galpha, std::vector& gbeta, std::vector& gcap, int nb, std::vector dep, std::vector arr, std::vector dem, double max_gap, int max_it, int method, std::vector& rank, std::vector& tail, std::vector& head, std::vector& first_out, std::vector& adj_head, std::vector& adj_arc, std::vector& rank_first_out, std::vector& rank_adj_head, std::vector& rank_adj_arc, std::vector& input_arc, std::vector& input_forward, std::vector& elimination_tree_parent, bool verbose); +RcppExport SEXP _cppRouting_cpptrafficcch(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP gflowSEXP, SEXP gauxSEXP, SEXP gfttSEXP, SEXP galphaSEXP, SEXP gbetaSEXP, SEXP gcapSEXP, SEXP nbSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP demSEXP, SEXP max_gapSEXP, SEXP max_itSEXP, SEXP methodSEXP, SEXP rankSEXP, SEXP tailSEXP, SEXP headSEXP, SEXP first_outSEXP, SEXP adj_headSEXP, SEXP adj_arcSEXP, SEXP rank_first_outSEXP, SEXP rank_adj_headSEXP, SEXP rank_adj_arcSEXP, SEXP input_arcSEXP, SEXP input_forwardSEXP, SEXP elimination_tree_parentSEXP, SEXP verboseSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< std::vector& >::type gfrom(gfromSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gto(gtoSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gw(gwSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gflow(gflowSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gaux(gauxSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gftt(gfttSEXP); + Rcpp::traits::input_parameter< std::vector& >::type galpha(galphaSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gbeta(gbetaSEXP); + Rcpp::traits::input_parameter< std::vector& >::type gcap(gcapSEXP); + Rcpp::traits::input_parameter< int >::type nb(nbSEXP); + Rcpp::traits::input_parameter< std::vector >::type dep(depSEXP); + Rcpp::traits::input_parameter< std::vector >::type arr(arrSEXP); + Rcpp::traits::input_parameter< std::vector >::type dem(demSEXP); + Rcpp::traits::input_parameter< double >::type max_gap(max_gapSEXP); + Rcpp::traits::input_parameter< int >::type max_it(max_itSEXP); + Rcpp::traits::input_parameter< int >::type method(methodSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank(rankSEXP); + Rcpp::traits::input_parameter< std::vector& >::type tail(tailSEXP); + Rcpp::traits::input_parameter< std::vector& >::type head(headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type first_out(first_outSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_head(adj_headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type adj_arc(adj_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank_first_out(rank_first_outSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank_adj_head(rank_adj_headSEXP); + Rcpp::traits::input_parameter< std::vector& >::type rank_adj_arc(rank_adj_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type input_arc(input_arcSEXP); + Rcpp::traits::input_parameter< std::vector& >::type input_forward(input_forwardSEXP); + Rcpp::traits::input_parameter< std::vector& >::type elimination_tree_parent(elimination_tree_parentSEXP); + Rcpp::traits::input_parameter< bool >::type verbose(verboseSEXP); + rcpp_result_gen = Rcpp::wrap(cpptrafficcch(gfrom, gto, gw, gflow, gaux, gftt, galpha, gbeta, gcap, nb, dep, arr, dem, max_gap, max_it, method, rank, tail, head, first_out, adj_head, adj_arc, rank_first_out, rank_adj_head, rank_adj_arc, input_arc, input_forward, elimination_tree_parent, verbose)); + return rcpp_result_gen; +END_RCPP +} // cpptraffic Rcpp::List cpptraffic(std::vector& gfrom, std::vector& gto, std::vector& gw, std::vector& gflow, std::vector& gaux, std::vector& gftt, std::vector& galpha, std::vector& gbeta, std::vector& gcap, int nb, std::vector lat, std::vector lon, double k, std::vector dep, std::vector arr, std::vector dem, double max_gap, int max_it, int method, int aon_method, bool contract, bool phast, bool verbose); RcppExport SEXP _cppRouting_cpptraffic(SEXP gfromSEXP, SEXP gtoSEXP, SEXP gwSEXP, SEXP gflowSEXP, SEXP gauxSEXP, SEXP gfttSEXP, SEXP galphaSEXP, SEXP gbetaSEXP, SEXP gcapSEXP, SEXP nbSEXP, SEXP latSEXP, SEXP lonSEXP, SEXP kSEXP, SEXP depSEXP, SEXP arrSEXP, SEXP demSEXP, SEXP max_gapSEXP, SEXP max_itSEXP, SEXP methodSEXP, SEXP aon_methodSEXP, SEXP contractSEXP, SEXP phastSEXP, SEXP verboseSEXP) { @@ -428,18 +731,30 @@ static const R_CallMethodDef CallEntries[] = { {"_cppRouting_cpppath", (DL_FUNC) &_cppRouting_cpppath, 13}, {"_cppRouting_cppdistmat", (DL_FUNC) &_cppRouting_cppdistmat, 6}, {"_cppRouting_cpppathmat", (DL_FUNC) &_cppRouting_cpppathmat, 12}, - {"_cppRouting_cppsimplify", (DL_FUNC) &_cppRouting_cppsimplify, 8}, + {"_cppRouting_cppsimplify", (DL_FUNC) &_cppRouting_cppsimplify, 9}, {"_cppRouting_cppcontract", (DL_FUNC) &_cppRouting_cppcontract, 5}, + {"_cppRouting_cppcchprepare", (DL_FUNC) &_cppRouting_cppcchprepare, 5}, + {"_cppRouting_cppcchcustomize", (DL_FUNC) &_cppRouting_cppcchcustomize, 12}, + {"_cppRouting_cppdistcch", (DL_FUNC) &_cppRouting_cppdistcch, 10}, + {"_cppRouting_cppdistmatcch", (DL_FUNC) &_cppRouting_cppdistmatcch, 8}, + {"_cppRouting_cpppathvaluescch", (DL_FUNC) &_cppRouting_cpppathvaluescch, 22}, + {"_cppRouting_cpppathcch", (DL_FUNC) &_cppRouting_cpppathcch, 24}, + {"_cppRouting_cppaoncchelim", (DL_FUNC) &_cppRouting_cppaoncchelim, 24}, {"_cppRouting_cppdistC", (DL_FUNC) &_cppRouting_cppdistC, 12}, {"_cppRouting_cpppathC", (DL_FUNC) &_cppRouting_cpppathC, 14}, {"_cppRouting_cppdistmatC", (DL_FUNC) &_cppRouting_cppdistmatC, 12}, {"_cppRouting_cpppathmatC", (DL_FUNC) &_cppRouting_cpppathmatC, 13}, {"_cppRouting_cpppadd", (DL_FUNC) &_cppRouting_cpppadd, 7}, + {"_cppRouting_cppmataddmulti", (DL_FUNC) &_cppRouting_cppmataddmulti, 7}, {"_cppRouting_cppdistadd", (DL_FUNC) &_cppRouting_cppdistadd, 11}, + {"_cppRouting_cppdistaddmulti", (DL_FUNC) &_cppRouting_cppdistaddmulti, 11}, {"_cppRouting_cppdistaddC", (DL_FUNC) &_cppRouting_cppdistaddC, 16}, + {"_cppRouting_cppdistaddCmulti", (DL_FUNC) &_cppRouting_cppdistaddCmulti, 15}, {"_cppRouting_cppaddC", (DL_FUNC) &_cppRouting_cppaddC, 16}, + {"_cppRouting_cppmataddCmulti", (DL_FUNC) &_cppRouting_cppmataddCmulti, 15}, {"_cppRouting_cppaon", (DL_FUNC) &_cppRouting_cppaon, 11}, {"_cppRouting_cppaonC", (DL_FUNC) &_cppRouting_cppaonC, 16}, + {"_cppRouting_cpptrafficcch", (DL_FUNC) &_cppRouting_cpptrafficcch, 29}, {"_cppRouting_cpptraffic", (DL_FUNC) &_cppRouting_cpptraffic, 23}, {"_cppRouting_cppalgB", (DL_FUNC) &_cppRouting_cppalgB, 27}, {NULL, NULL, 0} diff --git a/src/aggc.cpp b/src/aggc.cpp index 480dc1e..67b86ad 100644 --- a/src/aggc.cpp +++ b/src/aggc.cpp @@ -45,7 +45,7 @@ void aggC::operator()(std::size_t begin, std::size_t end){ } } - if (edge_index == -1) Rcpp::Rcout << dep << "->"<%d", dep, second); weight += m_or->add[edge_index]; dep = second; @@ -84,7 +84,7 @@ void aggC::operator()(std::size_t begin, std::size_t end){ } } - if (edge_index == -1) Rcpp::Rcout << dep << "->"<%d", dep, second); weight += m_or->add[edge_index]; dep = second; @@ -96,3 +96,107 @@ void aggC::operator()(std::size_t begin, std::size_t end){ } } + +aggCAdd::aggCAdd(CGraph* gr, Graph* original): + m_gr(gr), m_or(original){ + + n_add = m_or->n_add; + m_result.assign((size_t)m_gr->nodeG.size() * n_add, 0.0); + m_result2.assign((size_t)m_gr->nodeGr.size() * n_add, 0.0); +} + +void aggCAdd::operator()(std::size_t begin, std::size_t end){ + + int total = (int)m_or->nodeG.size(); // stride of m_or->add + int strideup = (int)m_gr->nodeG.size(); + int stridedown = (int)m_gr->nodeGr.size(); + + for (size_t k = begin; k != end; k++){ + // upgraph + for (int i = m_gr->indG[k]; i < m_gr->indG[k+1]; i++){ + + DVec weight(n_add, 0.0); + IVec result2(2); + result2[0] = k; + result2[1] = m_gr->nodeG[i]; + m_gr->unpack(result2); + + int pos = 0; + int dep = result2[pos]; + int second = result2[pos + 1]; + int arr = result2.back(); + + + while (dep != arr){ + + second = result2[pos + 1]; + int edge_index = -1; + double w = std::numeric_limits::max(); + for (int j = m_or->indG[dep]; j < m_or->indG[dep + 1]; j++){ + + if (m_or->nodeG[j] == second && m_or->wG[j] < w){ + w = m_or->wG[j]; + edge_index = j; + + } + } + + if (edge_index == -1) Rcpp::stop("aggCAdd: no matching original edge found for %d->%d", dep, second); + + for (int col = 0; col < n_add; col++){ + weight[col] += m_or->add[(size_t)col * total + edge_index]; + } + dep = second; + pos += 1; + } + + for (int col = 0; col < n_add; col++){ + m_result[(size_t)col * strideup + i] = weight[col]; + } + } + + // downgraph + for (int i = m_gr->indGr[k]; i < m_gr->indGr[k+1]; i++){ + + DVec weight(n_add, 0.0); + IVec result2(2); + result2[0] = m_gr->nodeGr[i];; + result2[1] = k; + m_gr->unpack(result2); + + int pos = 0; + int dep = result2[pos]; + int second = result2[pos + 1]; + int arr = result2.back(); + + + while (dep != arr){ + + second = result2[pos + 1]; + int edge_index = -1; + double w = std::numeric_limits::max(); + for (int j = m_or->indG[dep]; j < m_or->indG[dep + 1]; j++){ + + if (m_or->nodeG[j] == second && m_or->wG[j] < w){ + w = m_or->wG[j]; + edge_index = j; + + } + } + + if (edge_index == -1) Rcpp::stop("aggCAdd: no matching original edge found for %d->%d", dep, second); + + for (int col = 0; col < n_add; col++){ + weight[col] += m_or->add[(size_t)col * total + edge_index]; + } + dep = second; + pos += 1; + } + + for (int col = 0; col < n_add; col++){ + m_result2[(size_t)col * stridedown + i] = weight[col]; + } + } + + } +} diff --git a/src/aggc.h b/src/aggc.h index b68426d..2327af4 100644 --- a/src/aggc.h +++ b/src/aggc.h @@ -21,4 +21,21 @@ struct aggC : public Worker{ void operator()(std::size_t begin, std::size_t end); }; +// Same as aggC, generalized to N aux columns at once: for every shortcut +// edge in the contracted graph, sums each of the original graph's n_add +// columns over the unpacked original edges (instead of a single running +// sum), so distancePairCAdd never needs to unpack shortcuts during queries. +struct aggCAdd : public Worker{ + CGraph* m_gr; + Graph* m_or; + int n_add; + DVec m_result; // flattened (m_gr->nodeG.size() x n_add), column-major, upgraph + DVec m_result2; // flattened (m_gr->nodeGr.size() x n_add), column-major, downgraph + + aggCAdd(CGraph* gr, Graph* original); + + // iterate over m_gr->nbnode + void operator()(std::size_t begin, std::size_t end); +}; + #endif diff --git a/src/cch.h b/src/cch.h new file mode 100644 index 0000000..32007c8 --- /dev/null +++ b/src/cch.h @@ -0,0 +1,190 @@ +#ifndef CCH_H +#define CCH_H + +#include "graph.h" + +struct CCHPrepared { + int nbnode; + + // rank[node] is the contraction order position. Lower rank means the node is + // contracted earlier. CCH arcs always point from lower rank to higher rank. + IVec rank; + + // Topology of the prepared CCH, stored in original node ids. These arcs are + // independent from edge weights and can be saved/reused while costs change. + IVec tail; + IVec head; + IVec first_out; + IVec adj_head; + IVec adj_arc; + + // The same upward graph indexed by rank. Customization needs this form to + // scan triangles in contraction order. + IVec rank_first_out; + IVec rank_adj_head; + IVec rank_adj_arc; + + // Maps each input edge to the CCH arc that represents the same undirected + // connection, plus the direction of the edge on that arc. + IVec input_arc; + IVec input_forward; + + // Parent pointers for the RoutingKit-style elimination-tree query. Walking + // these parents visits the ancestors that can appear on an upward CCH search. + IVec elimination_tree_parent; +}; + +CCHPrepared build_cch(IVec &gfrom, IVec >o, int nb, IVec order, bool display_progress = false); + +void remap_cch_input_arcs(IVec &gfrom, IVec >o, CCHPrepared &prepared); + +void remap_cch_input_arcs(IVec &gfrom, + IVec >o, + CCHPrepared &prepared, + IVec &input_arc, + IVec &input_forward); + +void customize_cch(IVec &gfrom, + IVec >o, + DVec &gw, + int nb, + IVec &rank, + IVec &tail, + IVec &head, + IVec &rank_first_out, + IVec &rank_adj_head, + IVec &rank_adj_arc, + IVec &input_arc, + IVec &input_forward, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original); + +DVec distance_pair_cch(int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec dep, + IVec arr); + +Rcpp::NumericMatrix distance_matrix_cch(int nb, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + DVec &forward, + DVec &backward, + IVec dep, + IVec arr); + +vector path_pair_cch(int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &gfrom, + IVec >o, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec dep, + IVec arr, + IVec keep); + +Rcpp::NumericMatrix path_values_pair_cch(IVec &gfrom, + int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec dep, + IVec arr, + DVec &values, + int n_values); + +DVec aon_flow_cch_elimination_grouped(IVec &gfrom, + IVec >o, + DVec &gw, + int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec &dep, + IVec &arr, + DVec &demand); + +Rcpp::List aon_cch_elimination_grouped(IVec &gfrom, + IVec >o, + DVec &gw, + int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec dep, + IVec arr, + DVec demand); + +#endif diff --git a/src/cch_assign.cpp b/src/cch_assign.cpp new file mode 100644 index 0000000..daeeb0e --- /dev/null +++ b/src/cch_assign.cpp @@ -0,0 +1,319 @@ +#include "cch_internal.h" + +using namespace cch_internal; + +struct CCHEliminationGroupedAonWorker : public Worker { + // RoutingKit-style grouped elimination-tree AON query. + // + // Instead of running one independent shortest-path search per OD pair, this + // worker groups OD pairs by the side with fewer unique endpoints. A group then + // shares one search from that repeated endpoint and only performs the smaller + // pinned-endpoint work needed for the individual OD pairs. + // + // The recovered paths are unpacked to original graph edges immediately and + // accumulated in flow. That is why this worker is the traffic-assignment CCH + // path, while the distance query above is not. + IVec &gfrom; + int nb; + IVec &rank; + IVec &first_out; + IVec &adj_head; + IVec &adj_arc; + IVec &elimination_tree_parent; + DVec &forward; + DVec &backward; + IVec &forward_first_arc; + IVec &forward_first_dir; + IVec &forward_second_arc; + IVec &forward_second_dir; + IVec &forward_original; + IVec &backward_first_arc; + IVec &backward_first_dir; + IVec &backward_second_arc; + IVec &backward_second_dir; + IVec &backward_original; + IVec &dep; + IVec &arr; + DVec &demand; + vector &groups; + bool group_sources; + DVec flow; + + CCHEliminationGroupedAonWorker(IVec &gfrom, + int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec &dep, + IVec &arr, + DVec &demand, + vector &groups, + bool group_sources) + : gfrom(gfrom), nb(nb), rank(rank), first_out(first_out), adj_head(adj_head), adj_arc(adj_arc), + elimination_tree_parent(elimination_tree_parent), forward(forward), backward(backward), + forward_first_arc(forward_first_arc), forward_first_dir(forward_first_dir), + forward_second_arc(forward_second_arc), forward_second_dir(forward_second_dir), + forward_original(forward_original), + backward_first_arc(backward_first_arc), backward_first_dir(backward_first_dir), + backward_second_arc(backward_second_arc), backward_second_dir(backward_second_dir), + backward_original(backward_original), dep(dep), arr(arr), demand(demand), + groups(groups), group_sources(group_sources), flow(gfrom.size(), 0.0) {} + + CCHEliminationGroupedAonWorker(CCHEliminationGroupedAonWorker &other, Split) + : gfrom(other.gfrom), nb(other.nb), rank(other.rank), + first_out(other.first_out), adj_head(other.adj_head), adj_arc(other.adj_arc), + elimination_tree_parent(other.elimination_tree_parent), + forward(other.forward), backward(other.backward), + forward_first_arc(other.forward_first_arc), forward_first_dir(other.forward_first_dir), + forward_second_arc(other.forward_second_arc), forward_second_dir(other.forward_second_dir), + forward_original(other.forward_original), + backward_first_arc(other.backward_first_arc), backward_first_dir(other.backward_first_dir), + backward_second_arc(other.backward_second_arc), backward_second_dir(other.backward_second_dir), + backward_original(other.backward_original), dep(other.dep), arr(other.arr), demand(other.demand), + groups(other.groups), group_sources(other.group_sources), flow(other.gfrom.size(), 0.0) {} + + void join(CCHEliminationGroupedAonWorker &other){ + for (int i = 0; i < flow.size(); i++) flow[i] += other.flow[i]; + } + + void add_recovered_path(int start, + int root, + bool reverse_path, + double volume, + IVec &pred_node, + IVec &pred_arc, + IVec &pred_dir, + IVec &path_arcs, + IVec &path_dirs, + IVec &path_edges){ + int node = start; + while (node != root && pred_node[node] != -1){ + path_arcs.push_back(pred_arc[node]); + path_dirs.push_back(pred_dir[node]); + node = pred_node[node]; + } + if (node != root) { + path_arcs.clear(); + path_dirs.clear(); + return; + } + if (reverse_path){ + reverse(path_arcs.begin(), path_arcs.end()); + reverse(path_dirs.begin(), path_dirs.end()); + } + + for (int i = 0; i < path_arcs.size(); i++){ + unpack_metric_arc(path_arcs[i], path_dirs[i], + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + path_edges); + } + for (int i = 0; i < path_edges.size(); i++){ + if (path_edges[i] >= 0) flow[path_edges[i]] += volume; + } + } + + void operator()(std::size_t begin, std::size_t end){ + DVec dist(nb, INF); + IVec pred_node(nb, -1); + IVec pred_arc(nb, -1); + IVec pred_dir(nb, -1); + IVec is_touched(nb, 0); + IVec touched; + IVec is_pinned_ancestor(nb, 0); + IVec pinned_stack; + IVec pinned_touched; + IVec path_arcs; + IVec path_dirs; + IVec path_edges; + + for (std::size_t group_id = begin; group_id != end; group_id++){ + if (groups[group_id].size() == 0) continue; + + int root = static_cast(group_id); + touch_distance_node(root, dist, pred_node, pred_arc, pred_dir, is_touched, touched); + dist[root] = 0.0; + + // Shared phase: one upward search from the repeated endpoint. For source + // groups this is a forward search from the origin; for destination groups + // it is the symmetric backward search from the destination. + if (group_sources){ + for_ancestors(elimination_tree_parent, root, [&](int node){ + relax_outgoing_elimination(node, first_out, adj_head, adj_arc, forward, 1, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + return true; + }); + } else { + for_ancestors(elimination_tree_parent, root, [&](int node){ + relax_outgoing_elimination(node, first_out, adj_head, adj_arc, backward, 0, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + return true; + }); + } + + // Pinned phase: only visit ancestors of the OD endpoints in this group. + // Processing high ranks first is essential. It makes the dynamic program + // propagate distances from common ancestors back down toward each pinned + // endpoint before we recover paths. + for (int pos = 0; pos < groups[group_id].size(); pos++){ + int k = groups[group_id][pos]; + int pinned = group_sources ? arr[k] : dep[k]; + for_ancestors(elimination_tree_parent, pinned, [&](int node){ + if (is_pinned_ancestor[node] == 0){ + is_pinned_ancestor[node] = 1; + pinned_touched.push_back(node); + pinned_stack.push_back(node); + } + return true; + }); + } + + sort(pinned_stack.begin(), pinned_stack.end(), [&](int a, int b){ + return rank[a] > rank[b]; + }); + for (int stack_pos = 0; stack_pos < pinned_stack.size(); stack_pos++){ + int node = pinned_stack[stack_pos]; + if (group_sources){ + relax_incoming_elimination(node, first_out, adj_head, adj_arc, backward, 0, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + } else { + relax_incoming_elimination(node, first_out, adj_head, adj_arc, forward, 1, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + } + } + + // Recover and unpack one shortest path for every OD pair in the group. + // Shortcut unpacking writes flow back to original input edges, which is + // what traffic assignment needs for the next cost update. + for (int pos = 0; pos < groups[group_id].size(); pos++){ + int k = groups[group_id][pos]; + int start = group_sources ? arr[k] : dep[k]; + if (start == root || dist[start] == INF) continue; + + path_arcs.clear(); + path_dirs.clear(); + path_edges.clear(); + add_recovered_path(start, root, group_sources, demand[k], + pred_node, pred_arc, pred_dir, path_arcs, path_dirs, path_edges); + } + + clear_elimination_state(dist, pred_node, pred_arc, pred_dir, is_touched, touched); + for (int i = 0; i < pinned_touched.size(); i++) is_pinned_ancestor[pinned_touched[i]] = 0; + pinned_touched.clear(); + pinned_stack.clear(); + } + } +}; + +DVec aon_flow_cch_elimination_grouped(IVec &gfrom, + IVec >o, + DVec &gw, + int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec &dep, + IVec &arr, + DVec &demand){ + // Public native AON entry point for the current CCH implementation. The R API + // exposes this as aon_method = "cch". Older pairwise/grouped CCH query modes + // were removed so assignment always uses this elimination-tree path. + IVec source_count(nb, 0); + IVec target_count(nb, 0); + int unique_sources = 0; + int unique_targets = 0; + for (int i = 0; i < dep.size(); i++){ + if (source_count[dep[i]] == 0) unique_sources++; + if (target_count[arr[i]] == 0) unique_targets++; + source_count[dep[i]]++; + target_count[arr[i]]++; + } + + bool group_sources = unique_sources <= unique_targets; + vector groups(nb); + for (int i = 0; i < dep.size(); i++){ + int group_id = group_sources ? dep[i] : arr[i]; + groups[group_id].push_back(i); + } + + CCHEliminationGroupedAonWorker worker(gfrom, nb, rank, first_out, adj_head, adj_arc, elimination_tree_parent, + forward, backward, + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + dep, arr, demand, groups, group_sources); + if (dep.size() <= 1000 || nb <= 1000){ + worker(0, nb); + } else { + parallelReduce(0, nb, worker); + } + + return worker.flow; +} + +Rcpp::List aon_cch_elimination_grouped(IVec &gfrom, + IVec >o, + DVec &gw, + int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec dep, + IVec arr, + DVec demand){ + // Thin Rcpp-friendly wrapper around the native flow vector. It returns the + // same four-column shape as the other cppRouting AON backends. + DVec flow = aon_flow_cch_elimination_grouped(gfrom, gto, gw, nb, rank, first_out, adj_head, adj_arc, + elimination_tree_parent, forward, backward, + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + dep, arr, demand); + + Rcpp::List result(4); + result[0] = gfrom; + result[1] = gto; + result[2] = gw; + result[3] = flow; + return result; +} diff --git a/src/cch_customize.cpp b/src/cch_customize.cpp new file mode 100644 index 0000000..7e82e65 --- /dev/null +++ b/src/cch_customize.cpp @@ -0,0 +1,121 @@ +#include "cch_internal.h" + +using namespace cch_internal; + +void customize_cch(IVec &gfrom, + IVec >o, + DVec &gw, + int nb, + IVec &rank, + IVec &tail, + IVec &head, + IVec &rank_first_out, + IVec &rank_adj_head, + IVec &rank_adj_arc, + IVec &input_arc, + IVec &input_forward, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original){ + // Customization phase. + // + // The caller supplies one current cost per original edge. This function + // updates the prepared shortcut graph so every CCH arc stores the best known + // cost in each travel direction: + // - forward: lower-ranked node -> higher-ranked node; + // - backward: higher-ranked node -> lower-ranked node. + // + // For original arcs, *_original points to the input edge. For shortcuts, + // *_first_arc and *_second_arc describe how to unpack the shortcut later. + // CCH customization has two jobs: + // 1. Load current edge weights onto the prepared CCH arcs. + // 2. Relax all upward triangles so each shortcut stores the cheapest unpackable + // path for the current congestion state. + forward.assign(tail.size(), INF); + backward.assign(tail.size(), INF); + forward_first_arc.assign(tail.size(), -1); + forward_first_dir.assign(tail.size(), -1); + forward_second_arc.assign(tail.size(), -1); + forward_second_dir.assign(tail.size(), -1); + forward_original.assign(tail.size(), -1); + backward_first_arc.assign(tail.size(), -1); + backward_first_dir.assign(tail.size(), -1); + backward_second_arc.assign(tail.size(), -1); + backward_second_dir.assign(tail.size(), -1); + backward_original.assign(tail.size(), -1); + + for (int i = 0; i < gfrom.size(); i++){ + int arc = input_arc[i]; + if (arc < 0) continue; + + if (input_forward[i] == 1){ + if (gw[i] < forward[arc]){ + forward[arc] = gw[i]; + forward_original[arc] = i; + forward_first_arc[arc] = -1; + forward_second_arc[arc] = -1; + } + } else { + if (gw[i] < backward[arc]){ + backward[arc] = gw[i]; + backward_original[arc] = i; + backward_first_arc[arc] = -1; + backward_second_arc[arc] = -1; + } + } + } + + for (int lower_rank = 0; lower_rank < nb; lower_rank++){ + int begin = rank_first_out[lower_rank]; + int end = rank_first_out[lower_rank + 1]; + + for (int i = begin; i < end; i++){ + int mid_rank = rank_adj_head[i]; + int left = rank_adj_arc[i]; + int mid_cursor = rank_first_out[mid_rank]; + + for (int j = i + 1; j < end; j++){ + int upper_rank = rank_adj_head[j]; + int right = rank_adj_arc[j]; + int upper = find_rank_arc(mid_rank, upper_rank, rank_first_out, rank_adj_head, rank_adj_arc, mid_cursor); + if (upper < 0) continue; + + // left and right form lower -> mid -> upper in rank space. Depending on + // the original travel direction, this can improve either metric + // direction stored on the upper shortcut. + if (backward[left] != INF && forward[right] != INF){ + double candidate = backward[left] + forward[right]; + if (candidate < forward[upper]){ + forward[upper] = candidate; + forward_original[upper] = -1; + forward_first_arc[upper] = left; + forward_first_dir[upper] = 0; + forward_second_arc[upper] = right; + forward_second_dir[upper] = 1; + } + } + + if (backward[right] != INF && forward[left] != INF){ + double candidate = backward[right] + forward[left]; + if (candidate < backward[upper]){ + backward[upper] = candidate; + backward_original[upper] = -1; + backward_first_arc[upper] = right; + backward_first_dir[upper] = 0; + backward_second_arc[upper] = left; + backward_second_dir[upper] = 1; + } + } + } + } + } +} diff --git a/src/cch_internal.h b/src/cch_internal.h new file mode 100644 index 0000000..ed0cad9 --- /dev/null +++ b/src/cch_internal.h @@ -0,0 +1,342 @@ +#ifndef CCH_INTERNAL_H +#define CCH_INTERNAL_H + +#include "cch.h" +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; +using namespace RcppParallel; + +namespace cch_internal { + +const double INF = numeric_limits::max(); + +inline long long edge_key(int a, int b){ + if (a > b) std::swap(a, b); + return (static_cast(a) << 32) | static_cast(b); +} + +inline IVec identity_order(int nb){ + IVec order(nb); + for (int i = 0; i < nb; i++) order[i] = i; + return order; +} + +// [[Rcpp::depends(RcppProgress)]] +inline IVec min_degree_order(IVec &gfrom, IVec >o, int nb, bool display_progress = false){ + vector > graph(nb); + for (int i = 0; i < gfrom.size(); i++){ + if (gfrom[i] == gto[i]) continue; + graph[gfrom[i]].insert(gto[i]); + graph[gto[i]].insert(gfrom[i]); + } + + IVec order; + order.reserve(nb); + IVec removed(nb, 0); + + // Bucket-by-degree queue: buckets[d] holds, in ascending node-id order, the + // unremoved nodes whose current degree is d. This finds the minimum-degree + // node (ties broken by smallest id, same as a linear scan) in amortized + // near-O(1) instead of rescanning every node at each of the nb elimination + // steps, which made preparation O(nb^2) and dominated cpp_cch_prepare's + // runtime on large road networks. + IVec degree(nb); + int max_degree = 0; + for (int node = 0; node < nb; node++){ + degree[node] = (int)graph[node].size(); + max_degree = max(max_degree, degree[node]); + } + + vector > buckets(max_degree + 1); + for (int node = 0; node < nb; node++) buckets[degree[node]].insert(node); + + int min_ptr = 0; + + Progress p(nb, display_progress); + for (int step = 0; step < nb; step++){ + Rcpp::checkUserInterrupt(); + p.increment(); + + while (min_ptr < (int)buckets.size() && buckets[min_ptr].empty()) min_ptr++; + if (min_ptr >= (int)buckets.size()) break; + + int best = *buckets[min_ptr].begin(); + buckets[min_ptr].erase(buckets[min_ptr].begin()); + + IVec neighbors; + for (set::iterator it = graph[best].begin(); it != graph[best].end(); ++it){ + if (removed[*it] == 0) neighbors.push_back(*it); + } + + for (int i = 0; i < neighbors.size(); i++){ + for (int j = i + 1; j < neighbors.size(); j++){ + graph[neighbors[i]].insert(neighbors[j]); + graph[neighbors[j]].insert(neighbors[i]); + } + } + + for (int i = 0; i < neighbors.size(); i++) graph[neighbors[i]].erase(best); + graph[best].clear(); + removed[best] = 1; + order.push_back(best); + + // Fill-in only changes the degree of best's former neighbors: they lose + // the arc to best and gain any new clique edges among themselves. + for (int i = 0; i < neighbors.size(); i++){ + int node = neighbors[i]; + int new_degree = (int)graph[node].size(); + if (new_degree == degree[node]) continue; + + buckets[degree[node]].erase(node); + if (new_degree >= (int)buckets.size()) buckets.resize(new_degree + 1); + buckets[new_degree].insert(node); + degree[node] = new_degree; + if (new_degree < min_ptr) min_ptr = new_degree; + } + } + + return order; +} + +inline IVec rank_from_order(IVec order, int nb){ + if (order.size() == 0) order = identity_order(nb); + if (order.size() != nb) throw invalid_argument("order length must be equal to the number of nodes"); + + IVec rank(nb, -1); + for (int i = 0; i < nb; i++){ + int node = order[i]; + if (node < 0 || node >= nb) throw invalid_argument("order contains an invalid node id"); + if (rank[node] != -1) throw invalid_argument("order contains duplicated node ids"); + rank[node] = i; + } + return rank; +} + +inline void dijkstra_up(int start, + int nb, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + DVec &weight, + DVec &distances, + IVec &touched){ + // Unidirectional Dijkstra restricted to the upward CCH graph. Also records + // which nodes were actually reached, so callers can combine many such + // searches (e.g. a many-to-many distance matrix) without rescanning the + // full node set for every pair. + distances[start] = 0.0; + touched.push_back(start); + PQ queue; + queue.push(make_pair(start, 0.0)); + + while (!queue.empty()){ + int node = queue.top().first; + double dist = queue.top().second; + queue.pop(); + + if (dist > distances[node]) continue; + + for (int i = first_out[node]; i < first_out[node + 1]; i++){ + int arc = adj_arc[i]; + int next = adj_head[i]; + double edge_weight = weight[arc]; + if (edge_weight == INF) continue; + + double tentative = dist + edge_weight; + if (tentative < distances[next]){ + if (distances[next] == INF) touched.push_back(next); + distances[next] = tentative; + queue.push(make_pair(next, tentative)); + } + } + } +} + +inline void touch_node(int node, IVec &is_touched, IVec &touched){ + if (is_touched[node] == 0){ + is_touched[node] = 1; + touched.push_back(node); + } +} + +inline void sort_unique(IVec &values){ + sort(values.begin(), values.end()); + values.erase(unique(values.begin(), values.end()), values.end()); +} + +inline IVec merge_unique_tail(IVec &base, IVec &added){ + IVec result; + result.reserve(base.size() + added.size()); + merge(base.begin(), base.end(), added.begin(), added.end(), back_inserter(result)); + result.erase(unique(result.begin(), result.end()), result.end()); + return result; +} + +inline int find_rank_arc(int tail_rank, + int head_rank, + IVec &rank_first_out, + IVec &rank_adj_head, + IVec &rank_adj_arc, + int &cursor){ + int end = rank_first_out[tail_rank + 1]; + while (cursor < end && rank_adj_head[cursor] < head_rank) cursor++; + if (cursor < end && rank_adj_head[cursor] == head_rank) return rank_adj_arc[cursor]; + return -1; +} + +inline IVec elimination_tree_from_up_graph(int nb, IVec &first_out, IVec &adj_head){ + IVec parent(nb, -1); + for (int node = 0; node < nb; node++){ + if (first_out[node] != first_out[node + 1]) parent[node] = adj_head[first_out[node]]; + } + return parent; +} + +template +inline void for_ancestors(IVec &parent, int node, const F &f){ + while (node != -1){ + if (!f(node)) return; + node = parent[node]; + } +} + +inline void touch_distance_node(int node, DVec &dist, IVec &pred_node, IVec &pred_arc, IVec &pred_dir, + IVec &is_touched, IVec &touched){ + if (is_touched[node] == 0){ + is_touched[node] = 1; + touched.push_back(node); + } +} + +inline void relax_outgoing_elimination(int node, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + DVec &weight, + int arc_dir, + DVec &dist, + IVec &pred_node, + IVec &pred_arc, + IVec &pred_dir, + IVec &is_touched, + IVec &touched){ + if (dist[node] == INF) return; + for (int i = first_out[node]; i < first_out[node + 1]; i++){ + int arc = adj_arc[i]; + double edge_weight = weight[arc]; + if (edge_weight == INF) continue; + + int next = adj_head[i]; + double tentative = dist[node] + edge_weight; + if (tentative < dist[next]){ + touch_distance_node(next, dist, pred_node, pred_arc, pred_dir, is_touched, touched); + dist[next] = tentative; + pred_node[next] = node; + pred_arc[next] = arc; + pred_dir[next] = arc_dir; + } + } +} + +inline void relax_incoming_elimination(int node, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + DVec &weight, + int arc_dir, + DVec &dist, + IVec &pred_node, + IVec &pred_arc, + IVec &pred_dir, + IVec &is_touched, + IVec &touched){ + for (int i = first_out[node]; i < first_out[node + 1]; i++){ + int arc = adj_arc[i]; + double edge_weight = weight[arc]; + if (edge_weight == INF) continue; + + int next = adj_head[i]; + if (dist[next] == INF) continue; + double tentative = dist[next] + edge_weight; + if (tentative < dist[node]){ + touch_distance_node(node, dist, pred_node, pred_arc, pred_dir, is_touched, touched); + dist[node] = tentative; + pred_node[node] = next; + pred_arc[node] = arc; + pred_dir[node] = arc_dir; + } + } +} + +inline void clear_elimination_state(DVec &dist, + IVec &pred_node, + IVec &pred_arc, + IVec &pred_dir, + IVec &is_touched, + IVec &touched){ + for (int i = 0; i < touched.size(); i++){ + int node = touched[i]; + dist[node] = INF; + pred_node[node] = -1; + pred_arc[node] = -1; + pred_dir[node] = -1; + is_touched[node] = 0; + } + touched.clear(); +} + +inline void unpack_metric_arc(int arc, + int direction, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec &path_edges){ + if (arc < 0) return; + + if (direction == 1){ + if (forward_original[arc] >= 0){ + path_edges.push_back(forward_original[arc]); + return; + } + unpack_metric_arc(forward_first_arc[arc], forward_first_dir[arc], + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + path_edges); + unpack_metric_arc(forward_second_arc[arc], forward_second_dir[arc], + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + path_edges); + } else { + if (backward_original[arc] >= 0){ + path_edges.push_back(backward_original[arc]); + return; + } + unpack_metric_arc(backward_first_arc[arc], backward_first_dir[arc], + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + path_edges); + unpack_metric_arc(backward_second_arc[arc], backward_second_dir[arc], + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + path_edges); + } +} + +} // namespace cch_internal + +#endif diff --git a/src/cch_prepare.cpp b/src/cch_prepare.cpp new file mode 100644 index 0000000..ded4301 --- /dev/null +++ b/src/cch_prepare.cpp @@ -0,0 +1,176 @@ +#include "cch_internal.h" + +using namespace cch_internal; + +// [[Rcpp::depends(RcppProgress)]] +CCHPrepared build_cch(IVec &gfrom, IVec >o, int nb, IVec order, bool display_progress){ + // Preparation phase. + // + // Input: + // - original directed road edges, using integer node ids; + // - an optional node order. + // + // Output: + // - a directed upward shortcut graph; + // - the elimination tree used by fast grouped queries; + // - a map from each original road edge to the CCH arc where its cost enters. + // + // No travel-time/congestion weight is stored here. That is intentional: the + // result can be saved and reused as long as the road topology is unchanged. + if (order.size() == 0){ + if (display_progress) Rcpp::Rcout << "Computing elimination order..." << endl; + order = min_degree_order(gfrom, gto, nb, display_progress); + } + + CCHPrepared result; + result.nbnode = nb; + result.rank = rank_from_order(order, nb); + + // Build the upward graph induced by the order. Each input edge is first + // treated as undirected; direction is restored during customization. + vector upward(nb); + for (int i = 0; i < gfrom.size(); i++){ + int from_rank = result.rank[gfrom[i]]; + int to_rank = result.rank[gto[i]]; + if (from_rank == to_rank) continue; + if (from_rank < to_rank){ + upward[from_rank].push_back(to_rank); + } else { + upward[to_rank].push_back(from_rank); + } + } + + for (int node_rank = 0; node_rank < nb; node_rank++) sort_unique(upward[node_rank]); + + vector > rank_arcs; + rank_arcs.reserve(gfrom.size()); + + if (display_progress) Rcpp::Rcout << "Eliminating nodes..." << endl; + Progress p(nb, display_progress); + for (int node_rank = 0; node_rank < nb; node_rank++){ + Rcpp::checkUserInterrupt(); + p.increment(); + + IVec &neighbors = upward[node_rank]; + if (neighbors.size() == 0) continue; + + for (int i = 0; i < neighbors.size(); i++) rank_arcs.push_back(make_pair(node_rank, neighbors[i])); + + // Eliminate this rank by making all higher-ranked neighbors adjacent. The + // fill edges are what let later customization evaluate shortcut triangles. + int lowest_upward_neighbor = neighbors[0]; + if (neighbors.size() > 1){ + IVec fill_edges(neighbors.begin() + 1, neighbors.end()); + upward[lowest_upward_neighbor] = merge_unique_tail(upward[lowest_upward_neighbor], fill_edges); + } + } + + sort(rank_arcs.begin(), rank_arcs.end()); + rank_arcs.erase(unique(rank_arcs.begin(), rank_arcs.end()), rank_arcs.end()); + + result.tail.reserve(rank_arcs.size()); + result.head.reserve(rank_arcs.size()); + + vector cch_upward(nb); + for (int arc = 0; arc < rank_arcs.size(); arc++){ + int lower_rank = rank_arcs[arc].first; + int higher_rank = rank_arcs[arc].second; + int lower_node = order[lower_rank]; + int higher_node = order[higher_rank]; + + result.tail.push_back(lower_node); + result.head.push_back(higher_node); + cch_upward[lower_rank].push_back(higher_rank); + } + + result.rank_first_out.assign(nb + 1, 0); + for (int rank_id = 0; rank_id < nb; rank_id++) result.rank_first_out[rank_id + 1] = cch_upward[rank_id].size(); + for (int rank_id = 1; rank_id <= nb; rank_id++) result.rank_first_out[rank_id] += result.rank_first_out[rank_id - 1]; + + result.rank_adj_head.resize(result.tail.size()); + result.rank_adj_arc.resize(result.tail.size()); + IVec rank_cursor = result.rank_first_out; + for (int arc = 0; arc < rank_arcs.size(); arc++){ + int lower_rank = rank_arcs[arc].first; + int higher_rank = rank_arcs[arc].second; + int pos = rank_cursor[lower_rank]++; + result.rank_adj_head[pos] = higher_rank; + result.rank_adj_arc[pos] = arc; + } + + result.first_out.assign(nb + 1, 0); + for (int i = 0; i < result.tail.size(); i++) result.first_out[result.tail[i] + 1]++; + for (int i = 1; i <= nb; i++) result.first_out[i] += result.first_out[i - 1]; + + result.adj_head.resize(result.tail.size()); + result.adj_arc.resize(result.tail.size()); + IVec cursor = result.first_out; + for (int arc = 0; arc < result.tail.size(); arc++){ + int pos = cursor[result.tail[arc]]++; + result.adj_head[pos] = result.head[arc]; + result.adj_arc[pos] = arc; + } + result.elimination_tree_parent = elimination_tree_from_up_graph(nb, result.first_out, result.adj_head); + + // Remember where each original edge enters the CCH. Congestion assignment + // changes edge costs many times, so customization must not search for this + // mapping at every iteration. + unordered_map edge_to_arc; + edge_to_arc.reserve(result.tail.size() * 2); + for (int arc = 0; arc < result.tail.size(); arc++){ + edge_to_arc[edge_key(result.tail[arc], result.head[arc])] = arc; + } + + result.input_arc.assign(gfrom.size(), -1); + result.input_forward.assign(gfrom.size(), 0); + for (int i = 0; i < gfrom.size(); i++){ + int from = gfrom[i]; + int to = gto[i]; + if (from == to) continue; + + int lower = result.rank[from] < result.rank[to] ? from : to; + int higher = result.rank[from] < result.rank[to] ? to : from; + unordered_map::iterator found = edge_to_arc.find(edge_key(lower, higher)); + if (found == edge_to_arc.end()) continue; + + result.input_arc[i] = found->second; + result.input_forward[i] = (from == lower && to == higher) ? 1 : 0; + } + + return result; +} + +void remap_cch_input_arcs(IVec &gfrom, IVec >o, CCHPrepared &prepared){ + remap_cch_input_arcs(gfrom, gto, prepared, prepared.input_arc, prepared.input_forward); +} + +void remap_cch_input_arcs(IVec &gfrom, + IVec >o, + CCHPrepared &prepared, + IVec &input_arc, + IVec &input_forward){ + // Prepared CCH objects can be reused by assign_traffic(). The assignment + // graph has the same topology but its current adjacency edge order may differ, + // so only this cheap input-edge map is rebuilt. + unordered_map edge_to_arc; + edge_to_arc.reserve(prepared.tail.size() * 2); + for (int arc = 0; arc < prepared.tail.size(); arc++){ + edge_to_arc[edge_key(prepared.tail[arc], prepared.head[arc])] = arc; + } + + input_arc.assign(gfrom.size(), -1); + input_forward.assign(gfrom.size(), 0); + for (int i = 0; i < gfrom.size(); i++){ + int from = gfrom[i]; + int to = gto[i]; + if (from == to) continue; + + int lower = prepared.rank[from] < prepared.rank[to] ? from : to; + int higher = prepared.rank[from] < prepared.rank[to] ? to : from; + unordered_map::iterator found = edge_to_arc.find(edge_key(lower, higher)); + if (found == edge_to_arc.end()) continue; + + input_arc[i] = found->second; + input_forward[i] = (from == lower && to == higher) ? 1 : 0; + } +} diff --git a/src/cch_query.cpp b/src/cch_query.cpp new file mode 100644 index 0000000..c2b5545 --- /dev/null +++ b/src/cch_query.cpp @@ -0,0 +1,762 @@ +#include "cch_internal.h" + +using namespace cch_internal; + +struct CCHDistancePairWorker : public Worker { + // One-to-one distance queries on a customized CCH metric, using the + // elimination-tree query (RoutingKit-style): pairs are grouped by shared + // source or shared target so the ancestor-chain relaxation is shared + // across a whole group, instead of running a full bidirectional Dijkstra + // per pair. This is useful for distance APIs and correctness checks, but + // it does not recover paths and is not the traffic-assignment path. + int nb; + IVec &rank; + IVec &first_out; + IVec &adj_head; + IVec &adj_arc; + IVec &elimination_tree_parent; + DVec &forward; + DVec &backward; + IVec &dep; + IVec &arr; + vector &groups; + bool group_sources; + DVec &result; + + CCHDistancePairWorker(int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &dep, + IVec &arr, + vector &groups, + bool group_sources, + DVec &result) + : nb(nb), rank(rank), first_out(first_out), adj_head(adj_head), adj_arc(adj_arc), + elimination_tree_parent(elimination_tree_parent), forward(forward), backward(backward), + dep(dep), arr(arr), groups(groups), group_sources(group_sources), result(result) {} + + CCHDistancePairWorker(CCHDistancePairWorker &other, Split) + : nb(other.nb), rank(other.rank), first_out(other.first_out), adj_head(other.adj_head), adj_arc(other.adj_arc), + elimination_tree_parent(other.elimination_tree_parent), forward(other.forward), backward(other.backward), + dep(other.dep), arr(other.arr), groups(other.groups), group_sources(other.group_sources), result(other.result) {} + + void operator()(std::size_t begin, std::size_t end){ + DVec dist(nb, INF); + IVec pred_node(nb, -1); + IVec pred_arc(nb, -1); + IVec pred_dir(nb, -1); + IVec is_touched(nb, 0); + IVec touched; + IVec is_pinned_ancestor(nb, 0); + IVec pinned_stack; + IVec pinned_touched; + + for (std::size_t group_id = begin; group_id != end; group_id++){ + if (groups[group_id].size() == 0) continue; + + int root = static_cast(group_id); + touch_distance_node(root, dist, pred_node, pred_arc, pred_dir, is_touched, touched); + dist[root] = 0.0; + + if (group_sources){ + for_ancestors(elimination_tree_parent, root, [&](int node){ + relax_outgoing_elimination(node, first_out, adj_head, adj_arc, forward, 1, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + return true; + }); + } else { + for_ancestors(elimination_tree_parent, root, [&](int node){ + relax_outgoing_elimination(node, first_out, adj_head, adj_arc, backward, 0, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + return true; + }); + } + + for (int pos = 0; pos < groups[group_id].size(); pos++){ + int k = groups[group_id][pos]; + int pinned = group_sources ? arr[k] : dep[k]; + for_ancestors(elimination_tree_parent, pinned, [&](int node){ + if (is_pinned_ancestor[node] == 0){ + is_pinned_ancestor[node] = 1; + pinned_touched.push_back(node); + pinned_stack.push_back(node); + } + return true; + }); + } + + sort(pinned_stack.begin(), pinned_stack.end(), [&](int a, int b){ + return rank[a] > rank[b]; + }); + for (int stack_pos = 0; stack_pos < pinned_stack.size(); stack_pos++){ + int node = pinned_stack[stack_pos]; + if (group_sources){ + relax_incoming_elimination(node, first_out, adj_head, adj_arc, backward, 0, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + } else { + relax_incoming_elimination(node, first_out, adj_head, adj_arc, forward, 1, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + } + } + + for (int pos = 0; pos < groups[group_id].size(); pos++){ + int k = groups[group_id][pos]; + int start = group_sources ? arr[k] : dep[k]; + result[k] = dist[start]; + } + + clear_elimination_state(dist, pred_node, pred_arc, pred_dir, is_touched, touched); + for (int i = 0; i < pinned_touched.size(); i++) is_pinned_ancestor[pinned_touched[i]] = 0; + pinned_touched.clear(); + pinned_stack.clear(); + } + } +}; + +DVec distance_pair_cch(int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec dep, + IVec arr){ + // Distance-only pair query, using the elimination-tree query grouped by + // shared source or shared target. It computes shortest path costs but + // intentionally does not recover original road edges, so it cannot be + // used for AON flow loading. + IVec source_count(nb, 0); + IVec target_count(nb, 0); + int unique_sources = 0; + int unique_targets = 0; + for (int i = 0; i < dep.size(); i++){ + if (source_count[dep[i]] == 0) unique_sources++; + if (target_count[arr[i]] == 0) unique_targets++; + source_count[dep[i]]++; + target_count[arr[i]]++; + } + + bool group_sources = unique_sources <= unique_targets; + vector groups(nb); + for (int i = 0; i < dep.size(); i++){ + int group_id = group_sources ? dep[i] : arr[i]; + groups[group_id].push_back(i); + } + + DVec result(dep.size(), INF); + CCHDistancePairWorker worker(nb, rank, first_out, adj_head, adj_arc, elimination_tree_parent, + forward, backward, dep, arr, groups, group_sources, result); + if (dep.size() <= 1000 || nb <= 1000){ + worker(0, nb); + } else { + parallelFor(0, nb, worker); + } + + return result; +} + +Rcpp::NumericMatrix distance_matrix_cch(int nb, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + DVec &forward, + DVec &backward, + IVec dep, + IVec arr){ + // Distance matrix helper (RPHAST-style many-to-many). It precomputes one + // backward upward search per destination and, for every node it touches, + // drops (destination column, distance) into that node's bucket. Each + // origin then runs a single forward upward search and only combines with + // the buckets of the nodes it actually touches, instead of rescanning + // every node in the graph for every (origin, destination) pair. + Rcpp::NumericMatrix result(dep.size(), arr.size()); + fill(result.begin(), result.end(), INF); + + vector bucket_col(nb); + vector bucket_dist(nb); + DVec dist_backward(nb, INF); + IVec touched; + + for (int j = 0; j < arr.size(); j++){ + dijkstra_up(arr[j], nb, first_out, adj_head, adj_arc, backward, dist_backward, touched); + for (int t = 0; t < touched.size(); t++){ + int node = touched[t]; + bucket_col[node].push_back(j); + bucket_dist[node].push_back(dist_backward[node]); + dist_backward[node] = INF; + } + touched.clear(); + } + + DVec dist_forward(nb, INF); + for (int i = 0; i < dep.size(); i++){ + dijkstra_up(dep[i], nb, first_out, adj_head, adj_arc, forward, dist_forward, touched); + for (int t = 0; t < touched.size(); t++){ + int node = touched[t]; + double from_dist = dist_forward[node]; + IVec &cols = bucket_col[node]; + DVec &dists = bucket_dist[node]; + for (int k = 0; k < cols.size(); k++){ + double candidate = from_dist + dists[k]; + if (candidate < result(i, cols[k])) result(i, cols[k]) = candidate; + } + dist_forward[node] = INF; + } + touched.clear(); + } + + return result; +} + +struct CCHPathValuesPairWorker : public Worker { + // Pair query that routes on the customized CCH cost and accumulates one or + // more original-edge value columns along the selected shortest path. + IVec &gfrom; + int nb; + int n_values; + IVec &rank; + IVec &first_out; + IVec &adj_head; + IVec &adj_arc; + IVec &elimination_tree_parent; + DVec &forward; + DVec &backward; + IVec &forward_first_arc; + IVec &forward_first_dir; + IVec &forward_second_arc; + IVec &forward_second_dir; + IVec &forward_original; + IVec &backward_first_arc; + IVec &backward_first_dir; + IVec &backward_second_arc; + IVec &backward_second_dir; + IVec &backward_original; + IVec &dep; + IVec &arr; + DVec &values; + vector &groups; + bool group_sources; + DVec &result; + + CCHPathValuesPairWorker(IVec &gfrom, + int nb, + int n_values, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec &dep, + IVec &arr, + DVec &values, + vector &groups, + bool group_sources, + DVec &result) + : gfrom(gfrom), nb(nb), n_values(n_values), rank(rank), first_out(first_out), adj_head(adj_head), adj_arc(adj_arc), + elimination_tree_parent(elimination_tree_parent), forward(forward), backward(backward), + forward_first_arc(forward_first_arc), forward_first_dir(forward_first_dir), + forward_second_arc(forward_second_arc), forward_second_dir(forward_second_dir), + forward_original(forward_original), + backward_first_arc(backward_first_arc), backward_first_dir(backward_first_dir), + backward_second_arc(backward_second_arc), backward_second_dir(backward_second_dir), + backward_original(backward_original), dep(dep), arr(arr), values(values), + groups(groups), group_sources(group_sources), result(result) {} + + CCHPathValuesPairWorker(CCHPathValuesPairWorker &other, Split) + : gfrom(other.gfrom), nb(other.nb), n_values(other.n_values), rank(other.rank), + first_out(other.first_out), adj_head(other.adj_head), adj_arc(other.adj_arc), + elimination_tree_parent(other.elimination_tree_parent), + forward(other.forward), backward(other.backward), + forward_first_arc(other.forward_first_arc), forward_first_dir(other.forward_first_dir), + forward_second_arc(other.forward_second_arc), forward_second_dir(other.forward_second_dir), + forward_original(other.forward_original), + backward_first_arc(other.backward_first_arc), backward_first_dir(other.backward_first_dir), + backward_second_arc(other.backward_second_arc), backward_second_dir(other.backward_second_dir), + backward_original(other.backward_original), dep(other.dep), arr(other.arr), values(other.values), + groups(other.groups), group_sources(other.group_sources), result(other.result) {} + + void add_recovered_values(int pair_index, + int start, + int root, + bool reverse_path, + double cost, + IVec &pred_node, + IVec &pred_arc, + IVec &pred_dir, + IVec &path_arcs, + IVec &path_dirs, + IVec &path_edges){ + int n_pairs = dep.size(); + result[pair_index] = cost; + for (int col = 0; col < n_values; col++) result[(col + 1) * n_pairs + pair_index] = 0.0; + + if (start == root) return; + + int node = start; + while (node != root && pred_node[node] != -1){ + path_arcs.push_back(pred_arc[node]); + path_dirs.push_back(pred_dir[node]); + node = pred_node[node]; + } + if (node != root){ + result[pair_index] = INF; + for (int col = 0; col < n_values; col++) result[(col + 1) * n_pairs + pair_index] = INF; + return; + } + if (reverse_path){ + reverse(path_arcs.begin(), path_arcs.end()); + reverse(path_dirs.begin(), path_dirs.end()); + } + + for (int i = 0; i < path_arcs.size(); i++){ + unpack_metric_arc(path_arcs[i], path_dirs[i], + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + path_edges); + } + for (int i = 0; i < path_edges.size(); i++){ + int edge = path_edges[i]; + if (edge < 0) continue; + for (int col = 0; col < n_values; col++){ + result[(col + 1) * n_pairs + pair_index] += values[col * gfrom.size() + edge]; + } + } + } + + void operator()(std::size_t begin, std::size_t end){ + DVec dist(nb, INF); + IVec pred_node(nb, -1); + IVec pred_arc(nb, -1); + IVec pred_dir(nb, -1); + IVec is_touched(nb, 0); + IVec touched; + IVec is_pinned_ancestor(nb, 0); + IVec pinned_stack; + IVec pinned_touched; + IVec path_arcs; + IVec path_dirs; + IVec path_edges; + + for (std::size_t group_id = begin; group_id != end; group_id++){ + if (groups[group_id].size() == 0) continue; + + int root = static_cast(group_id); + touch_distance_node(root, dist, pred_node, pred_arc, pred_dir, is_touched, touched); + dist[root] = 0.0; + + if (group_sources){ + for_ancestors(elimination_tree_parent, root, [&](int node){ + relax_outgoing_elimination(node, first_out, adj_head, adj_arc, forward, 1, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + return true; + }); + } else { + for_ancestors(elimination_tree_parent, root, [&](int node){ + relax_outgoing_elimination(node, first_out, adj_head, adj_arc, backward, 0, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + return true; + }); + } + + for (int pos = 0; pos < groups[group_id].size(); pos++){ + int k = groups[group_id][pos]; + int pinned = group_sources ? arr[k] : dep[k]; + for_ancestors(elimination_tree_parent, pinned, [&](int node){ + if (is_pinned_ancestor[node] == 0){ + is_pinned_ancestor[node] = 1; + pinned_touched.push_back(node); + pinned_stack.push_back(node); + } + return true; + }); + } + + sort(pinned_stack.begin(), pinned_stack.end(), [&](int a, int b){ + return rank[a] > rank[b]; + }); + for (int stack_pos = 0; stack_pos < pinned_stack.size(); stack_pos++){ + int node = pinned_stack[stack_pos]; + if (group_sources){ + relax_incoming_elimination(node, first_out, adj_head, adj_arc, backward, 0, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + } else { + relax_incoming_elimination(node, first_out, adj_head, adj_arc, forward, 1, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + } + } + + for (int pos = 0; pos < groups[group_id].size(); pos++){ + int k = groups[group_id][pos]; + int start = group_sources ? arr[k] : dep[k]; + if (dist[start] == INF) continue; + + path_arcs.clear(); + path_dirs.clear(); + path_edges.clear(); + add_recovered_values(k, start, root, group_sources, dist[start], + pred_node, pred_arc, pred_dir, path_arcs, path_dirs, path_edges); + } + + clear_elimination_state(dist, pred_node, pred_arc, pred_dir, is_touched, touched); + for (int i = 0; i < pinned_touched.size(); i++) is_pinned_ancestor[pinned_touched[i]] = 0; + pinned_touched.clear(); + pinned_stack.clear(); + } + } +}; + +Rcpp::NumericMatrix path_values_pair_cch(IVec &gfrom, + int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec dep, + IVec arr, + DVec &values, + int n_values){ + IVec source_count(nb, 0); + IVec target_count(nb, 0); + int unique_sources = 0; + int unique_targets = 0; + for (int i = 0; i < dep.size(); i++){ + if (source_count[dep[i]] == 0) unique_sources++; + if (target_count[arr[i]] == 0) unique_targets++; + source_count[dep[i]]++; + target_count[arr[i]]++; + } + + bool group_sources = unique_sources <= unique_targets; + vector groups(nb); + for (int i = 0; i < dep.size(); i++){ + int group_id = group_sources ? dep[i] : arr[i]; + groups[group_id].push_back(i); + } + + DVec result((n_values + 1) * dep.size(), INF); + CCHPathValuesPairWorker worker(gfrom, nb, n_values, rank, first_out, adj_head, adj_arc, elimination_tree_parent, + forward, backward, + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + dep, arr, values, groups, group_sources, result); + if (dep.size() <= 1000 || nb <= 1000){ + worker(0, nb); + } else { + parallelFor(0, nb, worker); + } + + Rcpp::NumericMatrix out(dep.size(), n_values + 1); + for (int col = 0; col < n_values + 1; col++){ + for (int row = 0; row < dep.size(); row++){ + out(row, col) = result[col * dep.size() + row]; + } + } + return out; +} + +struct CCHPathPairWorker : public Worker { + // Pair query that recovers the original node path on a customized CCH + // metric, using the elimination-tree query grouped by shared source or + // shared target (same algorithm as CCHPathValuesPairWorker), then unpacks + // each CCH arc on the resulting path to original edges and converts them + // to nodes via gfrom/gto. + int nb; + IVec &rank; + IVec &first_out; + IVec &adj_head; + IVec &adj_arc; + IVec &elimination_tree_parent; + DVec &forward; + DVec &backward; + IVec &dep; + IVec &arr; + IVec &keep; + vector &groups; + bool group_sources; + vector &result; + // For unpacking CCH arcs to original edges and then to nodes + IVec &gfrom; + IVec >o; + IVec &forward_first_arc; + IVec &forward_first_dir; + IVec &forward_second_arc; + IVec &forward_second_dir; + IVec &forward_original; + IVec &backward_first_arc; + IVec &backward_first_dir; + IVec &backward_second_arc; + IVec &backward_second_dir; + IVec &backward_original; + + CCHPathPairWorker(int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &dep, + IVec &arr, + IVec &keep, + vector &groups, + bool group_sources, + vector &result, + IVec &gfrom, + IVec >o, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original) + : nb(nb), rank(rank), first_out(first_out), adj_head(adj_head), adj_arc(adj_arc), + elimination_tree_parent(elimination_tree_parent), + forward(forward), backward(backward), dep(dep), arr(arr), keep(keep), + groups(groups), group_sources(group_sources), result(result), + gfrom(gfrom), gto(gto), + forward_first_arc(forward_first_arc), forward_first_dir(forward_first_dir), + forward_second_arc(forward_second_arc), forward_second_dir(forward_second_dir), + forward_original(forward_original), + backward_first_arc(backward_first_arc), backward_first_dir(backward_first_dir), + backward_second_arc(backward_second_arc), backward_second_dir(backward_second_dir), + backward_original(backward_original) {} + + CCHPathPairWorker(CCHPathPairWorker &other, Split) + : nb(other.nb), rank(other.rank), first_out(other.first_out), adj_head(other.adj_head), adj_arc(other.adj_arc), + elimination_tree_parent(other.elimination_tree_parent), + forward(other.forward), backward(other.backward), dep(other.dep), arr(other.arr), keep(other.keep), + groups(other.groups), group_sources(other.group_sources), result(other.result), + gfrom(other.gfrom), gto(other.gto), + forward_first_arc(other.forward_first_arc), forward_first_dir(other.forward_first_dir), + forward_second_arc(other.forward_second_arc), forward_second_dir(other.forward_second_dir), + forward_original(other.forward_original), + backward_first_arc(other.backward_first_arc), backward_first_dir(other.backward_first_dir), + backward_second_arc(other.backward_second_arc), backward_second_dir(other.backward_second_dir), + backward_original(other.backward_original) {} + + void unpack_arcs_to_nodes(const IVec &path_arcs, + const IVec &path_dirs, + IVec &out_nodes){ + // Unpack CCH arcs to original edges, then convert to nodes. + IVec path_edges; + for (int i = 0; i < path_arcs.size(); i++){ + unpack_metric_arc(path_arcs[i], path_dirs[i], + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + path_edges); + } + + // Convert original edges to nodes. For each edge (u,v), add v if it's not + // the same as the last added node. + for (int i = 0; i < path_edges.size(); i++){ + int edge = path_edges[i]; + if (edge < 0) continue; + int u = gfrom[edge]; + int v = gto[edge]; + // We need to determine direction. Since unpack_metric_arc may have + // reversed the edge, we check which endpoint matches the previous node. + if (out_nodes.empty()){ + out_nodes.push_back(u); + out_nodes.push_back(v); + } else { + int last = out_nodes.back(); + if (last == u){ + out_nodes.push_back(v); + } else if (last == v){ + // Edge was traversed in reverse, add u + out_nodes.push_back(u); + } + } + } + } + + void operator()(std::size_t begin, std::size_t end){ + DVec dist(nb, INF); + IVec pred_node(nb, -1); + IVec pred_arc(nb, -1); + IVec pred_dir(nb, -1); + IVec is_touched(nb, 0); + IVec touched; + IVec is_pinned_ancestor(nb, 0); + IVec pinned_stack; + IVec pinned_touched; + IVec path_arcs; + IVec path_dirs; + + for (std::size_t group_id = begin; group_id != end; group_id++){ + if (groups[group_id].size() == 0) continue; + + int root = static_cast(group_id); + touch_distance_node(root, dist, pred_node, pred_arc, pred_dir, is_touched, touched); + dist[root] = 0.0; + + if (group_sources){ + for_ancestors(elimination_tree_parent, root, [&](int node){ + relax_outgoing_elimination(node, first_out, adj_head, adj_arc, forward, 1, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + return true; + }); + } else { + for_ancestors(elimination_tree_parent, root, [&](int node){ + relax_outgoing_elimination(node, first_out, adj_head, adj_arc, backward, 0, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + return true; + }); + } + + for (int pos = 0; pos < groups[group_id].size(); pos++){ + int k = groups[group_id][pos]; + int pinned = group_sources ? arr[k] : dep[k]; + for_ancestors(elimination_tree_parent, pinned, [&](int node){ + if (is_pinned_ancestor[node] == 0){ + is_pinned_ancestor[node] = 1; + pinned_touched.push_back(node); + pinned_stack.push_back(node); + } + return true; + }); + } + + sort(pinned_stack.begin(), pinned_stack.end(), [&](int a, int b){ + return rank[a] > rank[b]; + }); + for (int stack_pos = 0; stack_pos < pinned_stack.size(); stack_pos++){ + int node = pinned_stack[stack_pos]; + if (group_sources){ + relax_incoming_elimination(node, first_out, adj_head, adj_arc, backward, 0, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + } else { + relax_incoming_elimination(node, first_out, adj_head, adj_arc, forward, 1, + dist, pred_node, pred_arc, pred_dir, is_touched, touched); + } + } + + for (int pos = 0; pos < groups[group_id].size(); pos++){ + int k = groups[group_id][pos]; + int start = group_sources ? arr[k] : dep[k]; + + IVec path_result; + if (dist[start] != INF){ + if (start != root){ + path_arcs.clear(); + path_dirs.clear(); + + int node = start; + while (node != root && pred_node[node] != -1){ + path_arcs.push_back(pred_arc[node]); + path_dirs.push_back(pred_dir[node]); + node = pred_node[node]; + } + if (node == root){ + // group_sources: start=target, root=source, so the chain above + // was built target->source and needs reversing to source->target. + if (group_sources){ + reverse(path_arcs.begin(), path_arcs.end()); + reverse(path_dirs.begin(), path_dirs.end()); + } + unpack_arcs_to_nodes(path_arcs, path_dirs, path_result); + } + } + if (path_result.empty()) path_result.push_back(start); + } + result[k] = path_result; + } + + clear_elimination_state(dist, pred_node, pred_arc, pred_dir, is_touched, touched); + for (int i = 0; i < pinned_touched.size(); i++) is_pinned_ancestor[pinned_touched[i]] = 0; + pinned_touched.clear(); + pinned_stack.clear(); + } + } +}; + +vector path_pair_cch(int nb, + IVec &rank, + IVec &first_out, + IVec &adj_head, + IVec &adj_arc, + IVec &elimination_tree_parent, + DVec &forward, + DVec &backward, + IVec &gfrom, + IVec >o, + IVec &forward_first_arc, + IVec &forward_first_dir, + IVec &forward_second_arc, + IVec &forward_second_dir, + IVec &forward_original, + IVec &backward_first_arc, + IVec &backward_first_dir, + IVec &backward_second_arc, + IVec &backward_second_dir, + IVec &backward_original, + IVec dep, + IVec arr, + IVec keep){ + IVec source_count(nb, 0); + IVec target_count(nb, 0); + int unique_sources = 0; + int unique_targets = 0; + for (int i = 0; i < dep.size(); i++){ + if (source_count[dep[i]] == 0) unique_sources++; + if (target_count[arr[i]] == 0) unique_targets++; + source_count[dep[i]]++; + target_count[arr[i]]++; + } + + bool group_sources = unique_sources <= unique_targets; + vector groups(nb); + for (int i = 0; i < dep.size(); i++){ + int group_id = group_sources ? dep[i] : arr[i]; + groups[group_id].push_back(i); + } + + vector result(dep.size()); + CCHPathPairWorker worker(nb, rank, first_out, adj_head, adj_arc, elimination_tree_parent, + forward, backward, dep, arr, keep, groups, group_sources, result, + gfrom, gto, forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original); + if (dep.size() <= 1000 || nb <= 1000){ + worker(0, nb); + } else { + parallelFor(0, nb, worker); + } + return result; +} diff --git a/src/cgraph.h b/src/cgraph.h index 7ccceec..12ce835 100644 --- a/src/cgraph.h +++ b/src/cgraph.h @@ -36,6 +36,7 @@ class CGraph{ // additional weight DVec add; DVec addr; + int n_add = 1; // number of columns flattened into add/addr (multi-column aggregate_aux path only) // constructors CGraph(Graph* graph); //from a c++ Graph object diff --git a/src/distance_mat.cpp b/src/distance_mat.cpp index cebc99f..c8e1986 100644 --- a/src/distance_mat.cpp +++ b/src/distance_mat.cpp @@ -79,6 +79,87 @@ void distanceMat::dijkstra_mat(std::size_t begin, std::size_t end){ } } +// constructor +distanceMatAdd::distanceMatAdd(Graph *gr, IVec dep, IVec arr): + m_gr(gr), m_dep(dep), m_arr(arr) +{ + n_add = m_gr->n_add; + m_result.assign((size_t)m_dep.size() * m_arr.size() * n_add, 0.0); +} + +void distanceMatAdd::operator()(std::size_t begin, std::size_t end){ + dijkstra_mat(begin, end); +} + +// same as distanceMat::dijkstra_mat, generalized to n_add running sums, with +// touched-array resets instead of full fills +void distanceMatAdd::dijkstra_mat(std::size_t begin, std::size_t end){ + + int n_dep = (int)m_dep.size(); + int n_arr = (int)m_arr.size(); + int total = (int)m_gr->nodeG.size(); + + DVec distances(m_gr->nbnode, numeric_limits::max()); + DVec distadd((size_t)m_gr->nbnode * n_add, numeric_limits::max()); + IVec touched; + + for (std::size_t k = begin; k != end; k++){ + + int start = m_dep[k]; + distances[start] = 0.0; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + start] = 0.0; + } + touched.push_back(start); + + PQ Q; + Q.push(make_pair(start, 0.0)); + + while (!Q.empty()) { + int v = Q.top().first; + double w = Q.top().second; + Q.pop(); + + if (w <= distances[v]) { + + for (int i = m_gr->indG[v]; i< m_gr->indG[v+1]; i++) { + + int v2 = m_gr->nodeG[i]; + double w2 = m_gr->wG[i]; + + if (distances[v] + w2 < distances[v2]) { + distances[v2] = distances[v] + w2; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + v2] = distadd[(size_t)col * m_gr->nbnode + v] + m_gr->add[(size_t)col * total + i]; + } + + Q.push(make_pair(v2, distances[v2])); + touched.push_back(v2); + } + } + } + } + + + for (int i = 0; i != n_arr; i++){ + int node = m_arr[i]; + for (int col = 0; col < n_add; col++){ + m_result[(size_t)col * n_dep * n_arr + (size_t)i * n_dep + k] = distadd[(size_t)col * m_gr->nbnode + node]; + } + } + + //Reinitialize (touched-only) + for (int t = 0; t < (int)touched.size(); t++){ + int node = touched[t]; + distances[node] = numeric_limits::max(); + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + node] = numeric_limits::max(); + } + } + touched.clear(); + } +} + ///////////////////////////////////////////////////// contracted graph @@ -261,3 +342,98 @@ void phastC::operator()(std::size_t begin, std::size_t end){ } } + +// constructor +phastCAdd::phastCAdd(CGraph *gr, IVec dep, IVec arr): + m_gr(gr), m_dep(dep), m_arr(arr) +{ + n_add = m_gr->n_add; + m_result.assign((size_t)m_dep.size() * m_arr.size() * n_add, 0.0); +} + +// same as phastC, generalized to n_add running sums; add/addr already hold +// per-shortcut sums precomputed by aggCAdd, so no unpacking happens here +void phastCAdd::operator()(std::size_t begin, std::size_t end){ + + int n_dep = (int)m_dep.size(); + int n_arr = (int)m_arr.size(); + int total = (int)m_gr->nodeG.size(); + int totalr = (int)m_gr->nodeGr.size(); + + DVec Distances(m_gr->nbnode, std::numeric_limits::max()); + DVec distadd((size_t)m_gr->nbnode * n_add, std::numeric_limits::max()); + + for (std::size_t k=begin; k!=end;k++){ + + int StartNode=m_dep[k]; + + Distances[StartNode] = 0.0; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + StartNode] = 0.0; + } + PQ Q; + Q.push(std::make_pair(StartNode, 0.0)); + + + while (true) { + + if (Q.empty()){ + break; + } + + if (!Q.empty()){ + + int v=Q.top().first; + double w=Q.top().second; + Q.pop(); + + if (w <= Distances[v]) { + if (Stall_par(v, Distances, m_gr->nodeGr, m_gr->wGr, m_gr->indGr)) continue; + for (int i = m_gr->indG[v]; i< m_gr->indG[v+1]; i++){ + int v2 = m_gr->nodeG[i]; + double w2 = m_gr->wG[i]; + + + if (Distances[v] + w2 < Distances[v2]) { + Distances[v2] = Distances[v] + w2; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + v2] = distadd[(size_t)col * m_gr->nbnode + v] + m_gr->add[(size_t)col * total + i]; + } + Q.push(std::make_pair(v2, Distances[v2])); + + } + } + } + + } + } + + //Backward + + for (int i=0; i < (m_gr->indGr.size()-1); i++){ + + for (int j=m_gr->indGr[i]; j < m_gr->indGr[i+1]; j++){ + + + if (Distances[m_gr->nodeGr[j]]+m_gr->wGr[j] < Distances[i]) { + Distances[i] = Distances[m_gr->nodeGr[j]]+m_gr->wGr[j]; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + i] = distadd[(size_t)col * m_gr->nbnode + m_gr->nodeGr[j]] + m_gr->addr[(size_t)col * totalr + j]; + } + } + + } + } + + for (int i=0; i != n_arr; i++){ + int node = m_arr[i]; + for (int col = 0; col < n_add; col++){ + m_result[(size_t)col * n_dep * n_arr + (size_t)i * n_dep + k] = distadd[(size_t)col * m_gr->nbnode + node]; + } + } + + std::fill(Distances.begin(),Distances.end(),std::numeric_limits::max()); + std::fill(distadd.begin(),distadd.end(),std::numeric_limits::max()); + + } +} diff --git a/src/distance_mat.h b/src/distance_mat.h index 762e0b8..29be509 100644 --- a/src/distance_mat.h +++ b/src/distance_mat.h @@ -32,6 +32,26 @@ struct distanceMat : public Worker{ }; +// Same as distanceMat, generalized to n_add aux columns summed inline. +// m_result is flattened n_dep x n_arr x n_add, column-major (R array +// convention: dep varies fastest, then arr, then column) so the wrapper can +// hand it back to R as a 3-D array with dim = c(n_dep, n_arr, n_add). +struct distanceMatAdd : public Worker{ + + const Graph *m_gr; + IVec m_dep; + IVec m_arr; + int n_add; + DVec m_result; + + distanceMatAdd(Graph *gr, IVec dep, IVec arr); + + void dijkstra_mat(std::size_t begin, std::size_t end); + + void operator()(std::size_t begin, std::size_t end); + +}; + struct distanceMatC : public Worker{ @@ -72,4 +92,22 @@ struct phastC : public Worker{ }; +// Same as phastC, generalized to n_add aux columns summed inline. m_gr->add +// / m_gr->addr are expected to already hold per-shortcut N-column sums +// (precomputed by aggCAdd), same contract as distancePairCAdd. m_result is +// flattened n_dep x n_arr x n_add, column-major (R array convention). +struct phastCAdd : public Worker{ + + CGraph *m_gr; + IVec m_dep; + IVec m_arr; + int n_add; + DVec m_result; + + phastCAdd(CGraph *gr, IVec dep, IVec arr); + + void operator()(std::size_t begin, std::size_t end); + +}; + #endif diff --git a/src/distance_pair.cpp b/src/distance_pair.cpp index 7769523..20f3bce 100644 --- a/src/distance_pair.cpp +++ b/src/distance_pair.cpp @@ -480,6 +480,532 @@ void distancePair::nba(std::size_t begin, std::size_t end){ } +//////////////////////////////////////////////////////////// Normal graph, N aux columns + +// constructor +distancePairAdd::distancePairAdd(Graph *gr, IVec dep, IVec arr, int algo) : + m_gr(gr), m_dep(dep), m_arr(arr), algorithm(algo) +{ + n_add = m_gr->n_add; + m_result.assign((size_t)m_dep.size() * n_add, 0.0); +} + +void distancePairAdd::operator()(std::size_t begin, std::size_t end){ + if (algorithm == 0) dijkstra_early_stop(begin, end); + if (algorithm == 1) bidir(begin, end); + if (algorithm == 2) astar(begin, end); + if (algorithm == 3) nba(begin, end); +} + +// dijkstra early stop, N aux columns summed inline (same structure as +// distancePair::dijkstra_early_stop, generalized to n_add running sums, with +// touched-array resets instead of full fills) +void distancePairAdd::dijkstra_early_stop(std::size_t begin, std::size_t end){ + + int n_pairs = (int)m_dep.size(); + int total = (int)m_gr->nodeG.size(); + + std::vector Distances(m_gr->nbnode, std::numeric_limits::max()); + DVec distadd((size_t)m_gr->nbnode * n_add, std::numeric_limits::max()); + IVec touched; + + for (std::size_t k = begin; k != end; k++){ + + int StartNode=m_dep[k]; + int EndNode=m_arr[k]; + + Distances[StartNode] = 0.0; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + StartNode] = 0.0; + } + touched.push_back(StartNode); + + PQ Q; + Q.push(std::make_pair(StartNode, 0.0)); + + while (!Q.empty()) { + int v = Q.top().first; + double w = Q.top().second; + Q.pop(); + + if (w <= Distances[v]) { + + for (int i = m_gr->indG[v]; i< m_gr->indG[v+1]; i++){ + int v2 = m_gr->nodeG[i]; + double w2 = m_gr->wG[i]; + + if (Distances[v] + w2 < Distances[v2]) { + Distances[v2] = Distances[v] + w2; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + v2] = distadd[(size_t)col * m_gr->nbnode + v] + m_gr->add[(size_t)col * total + i]; + } + + Q.push(std::make_pair(v2, Distances[v2])); + touched.push_back(v2); + } + } + } + if (v == EndNode){ + break; + } + } + + for (int col = 0; col < n_add; col++){ + m_result[(size_t)col * n_pairs + k] = distadd[(size_t)col * m_gr->nbnode + EndNode]; + } + + for (int i = 0; i < (int)touched.size(); i++){ + int node = touched[i]; + Distances[node] = std::numeric_limits::max(); + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + node] = std::numeric_limits::max(); + } + } + touched.clear(); + + } + +} + +// bidirectional dijkstra, N aux columns summed inline (same structure as +// distancePair::bidir, generalized from one running sum to n_add of them) +void distancePairAdd::bidir(std::size_t begin, std::size_t end){ + + int n_pairs = (int)m_dep.size(); + int total = (int)m_gr->nodeG.size(); // stride of m_gr->add + int totalr = (int)m_gr->nodeGr.size(); // stride of m_gr->addr + + std::vector Distances(m_gr->nbnode, std::numeric_limits::max()); + std::vector Distances2(m_gr->nbnode, std::numeric_limits::max()); + DVec distadd((size_t)m_gr->nbnode * n_add, std::numeric_limits::max()); + DVec distadd2((size_t)m_gr->nbnode * n_add, std::numeric_limits::max()); + DVec mu2(n_add, std::numeric_limits::max()); + + std::vector Visited(m_gr->nbnode,0); + std::vector Visiting(m_gr->nbnode,0); + std::vector Visited2(m_gr->nbnode,0); + std::vector Visiting2(m_gr->nbnode,0); + IVec touched; + + for (std::size_t k=begin; k!=end;k++){ + + int StartNode = m_dep[k]; + int EndNode=m_arr[k]; + Distances[StartNode] = 0.0; + Distances2[EndNode] = 0.0; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + StartNode] = 0.0; + distadd2[(size_t)col * m_gr->nbnode + EndNode] = 0.0; + } + touched.push_back(StartNode); + touched.push_back(EndNode); + + PQ Q; + PQ Qr; + Q.push(std::make_pair(StartNode, 0.0)); + Qr.push(std::make_pair(EndNode, 0.0)); + Visiting[StartNode]=1; + Visiting2[EndNode]=1; + + double mu=std::numeric_limits::max(); + std::fill(mu2.begin(), mu2.end(), std::numeric_limits::max()); + + while (!Q.empty() && !Qr.empty()) { + if (Q.top().second+Qr.top().second >= mu){ + break; + } + + + if (!Q.empty()){ + int v=Q.top().first; + double w=Q.top().second; + Q.pop(); + + + Visited[v]=1; + + + if (w <= Distances[v]) { + for (int i = m_gr->indG[v]; i< m_gr->indG[v+1]; i++){ + int v2 = m_gr->nodeG[i]; + double w2 = m_gr->wG[i]; + + + if (Distances[v] + w2 < Distances[v2]) { + Distances[v2] = Distances[v] + w2; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + v2] = distadd[(size_t)col * m_gr->nbnode + v] + m_gr->add[(size_t)col * total + i]; + } + + Q.push(std::make_pair(v2, Distances[v2])); + Visiting[v2]=1; + touched.push_back(v2); + + } + } + } + if ((Visited2[v]==1 || Visiting2[v]==1) && (Distances[v]+Distances2[v]) < mu){ + + mu=Distances[v]+Distances2[v]; + for (int col = 0; col < n_add; col++){ + mu2[col] = distadd[(size_t)col * m_gr->nbnode + v] + distadd2[(size_t)col * m_gr->nbnode + v]; + } + + } + } + + if (!Qr.empty()){ + int vv=Qr.top().first; + double ww=Qr.top().second; + Qr.pop(); + + + Visited2[vv]=1; + + if (ww <= Distances2[vv]) { + for (int i=m_gr->indGr[vv]; i< m_gr->indGr[vv+1]; i++){ + int vv2 = m_gr->nodeGr[i]; + double ww2 = m_gr->wGr[i]; + + + if (Distances2[vv] + ww2 < Distances2[vv2]) { + Distances2[vv2] = Distances2[vv] + ww2; + for (int col = 0; col < n_add; col++){ + distadd2[(size_t)col * m_gr->nbnode + vv2] = distadd2[(size_t)col * m_gr->nbnode + vv] + m_gr->addr[(size_t)col * totalr + i]; + } + + Qr.push(std::make_pair(vv2, Distances2[vv2])); + Visiting2[vv2]=1; + touched.push_back(vv2); + } + } + } + if ((Visited[vv]==1 || Visiting[vv]==1) && (Distances[vv]+Distances2[vv]) < mu){ + + mu=Distances[vv]+Distances2[vv]; + for (int col = 0; col < n_add; col++){ + mu2[col] = distadd[(size_t)col * m_gr->nbnode + vv] + distadd2[(size_t)col * m_gr->nbnode + vv]; + } + + } + } + + } + + for (int col = 0; col < n_add; col++){ + m_result[(size_t)col * n_pairs + k] = mu2[col]; + } + + for (int i = 0; i < touched.size(); i++){ + int node = touched[i]; + Distances[node] = std::numeric_limits::max(); + Distances2[node] = std::numeric_limits::max(); + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + node] = std::numeric_limits::max(); + distadd2[(size_t)col * m_gr->nbnode + node] = std::numeric_limits::max(); + } + Visited[node] = 0; + Visited2[node] = 0; + Visiting[node] = 0; + Visiting2[node] = 0; + } + touched.clear(); + + } +} + +// A*, N aux columns summed inline (same structure as distancePair::astar, +// generalized to n_add running sums, with touched-array resets instead of +// full fills) +void distancePairAdd::astar(std::size_t begin, std::size_t end){ + + int n_pairs = (int)m_dep.size(); + int total = (int)m_gr->nodeG.size(); + + std::vector Distances(m_gr->nbnode, std::numeric_limits::max()); + DVec distadd((size_t)m_gr->nbnode * n_add, std::numeric_limits::max()); + std::vector Distances2(m_gr->nbnode, std::numeric_limits::max()); + std::vector closedList(m_gr->nbnode,0); + std::vector openList(m_gr->nbnode,0); + IVec touched; + + for (std::size_t k=begin; k!=end;k++){ + + int StartNode=m_dep[k]; + int endNode=m_arr[k]; + double lata=m_gr->lat[endNode]; + double lona=m_gr->lon[endNode]; + + Distances[StartNode] = 0.0; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + StartNode] = 0.0; + } + + Distances2[StartNode] = sqrt(pow(m_gr->lat[StartNode]-lata,2)+pow(m_gr->lon[StartNode]-lona,2))/m_gr->k; + PQ Q; + Q.push(std::make_pair(StartNode,Distances2[StartNode])); + openList[StartNode]=1; + touched.push_back(StartNode); + + while (!Q.empty()) { + int v = Q.top().first; + Q.pop(); + if (closedList[v]==1){ + continue; + } + openList[v]=0; + closedList[v]=1; + + for (int i=m_gr->indG[v]; i< m_gr->indG[v+1]; i++){ + int v2 = m_gr->nodeG[i]; + double w2 = m_gr->wG[i]; + if (closedList[v2]==1) { + continue; + } + + double temp; + temp = Distances[v] + w2; + if (openList[v2]==0){ + + Q.push(std::make_pair(v2,Distances2[v2])); + openList[v2]=1; + touched.push_back(v2); + } + + + else if (temp>=Distances[v2]){ + continue; + } + + + Distances[v2]=temp; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + v2] = distadd[(size_t)col * m_gr->nbnode + v] + m_gr->add[(size_t)col * total + i]; + } + + Distances2[v2]=Distances[v2]+sqrt(pow(m_gr->lat[v2]-lata,2)+pow(m_gr->lon[v2]-lona,2))/m_gr->k; + Q.push(std::make_pair(v2,Distances2[v2])); + openList[v2]=1; + } + + + if (v==endNode){ + break; + } + + } + + for (int col = 0; col < n_add; col++){ + m_result[(size_t)col * n_pairs + k] = distadd[(size_t)col * m_gr->nbnode + endNode]; + } + + for (int i = 0; i < (int)touched.size(); i++){ + int node = touched[i]; + Distances[node] = std::numeric_limits::max(); + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + node] = std::numeric_limits::max(); + } + Distances2[node] = std::numeric_limits::max(); + closedList[node] = 0; + openList[node] = 0; + } + touched.clear(); + + } +} + +// bidirectional A*, N aux columns summed inline (same structure as +// distancePair::nba, generalized to n_add running sums, with touched-array +// resets instead of full fills) +void distancePairAdd::nba(std::size_t begin, std::size_t end){ + + int n_pairs = (int)m_dep.size(); + int total = (int)m_gr->nodeG.size(); + int totalr = (int)m_gr->nodeGr.size(); + + std::vector Distances(m_gr->nbnode, std::numeric_limits::max()); + std::vector Distances2(m_gr->nbnode, std::numeric_limits::max()); + DVec distadd((size_t)m_gr->nbnode * n_add, std::numeric_limits::max()); + DVec distadd2((size_t)m_gr->nbnode * n_add, std::numeric_limits::max()); + DVec mu2(n_add, std::numeric_limits::max()); + std::vector Visited(m_gr->nbnode,0); + std::vector Visited1check(m_gr->nbnode,0); + std::vector Visited2check(m_gr->nbnode,0); + IVec touched; + + for (std::size_t k=begin; k!=end;k++){ + + int StartNode=m_dep[k]; + int EndNode=m_arr[k]; + + if (StartNode==EndNode){ + for (int col = 0; col < n_add; col++){ + m_result[(size_t)col * n_pairs + k] = 0; + } + continue; + } + + double lata=m_gr->lat[EndNode]; + double lona=m_gr->lon[EndNode]; + double lata2=m_gr->lat[StartNode]; + double lona2=m_gr->lon[StartNode]; + + Distances[StartNode] = 0.0; + Visited1check[StartNode]=1; + Distances2[EndNode] = 0.0; + Visited2check[EndNode]=1; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + StartNode] = 0.0; + distadd2[(size_t)col * m_gr->nbnode + EndNode] = 0.0; + } + touched.push_back(StartNode); + touched.push_back(EndNode); + + PQ Q; + PQ Qr; + Q.push(std::make_pair(StartNode, sqrt(pow(m_gr->lat[StartNode]-lata,2)+pow(m_gr->lon[StartNode]-lona,2))/m_gr->k)); + Qr.push(std::make_pair(EndNode, sqrt(pow(m_gr->lat[EndNode]-lata2,2)+pow(m_gr->lon[EndNode]-lona2,2))/m_gr->k)); + + + + double total1=sqrt(pow(lata2-lata,2)+pow(lona2-lona,2))/m_gr->k; + double total2=sqrt(pow(lata2-lata,2)+pow(lona2-lona,2))/m_gr->k; + double mu=std::numeric_limits::max(); + std::fill(mu2.begin(), mu2.end(), std::numeric_limits::max()); + + while (!Q.empty() && !Qr.empty()) { + //Forward + if (Q.size() < Qr.size()){ + int v=Q.top().first; + Q.pop(); + if (Visited[v]==0){ + Visited[v]=1; + + if ((Distances[v] + sqrt(pow(m_gr->lat[v]-lata,2)+pow(m_gr->lon[v]-lona,2))/m_gr->k) >= mu || (Distances[v] + total2 - sqrt(pow(m_gr->lat[v]-lata2,2)+pow(m_gr->lon[v]-lona2,2))/m_gr->k) >= mu){} + + else { + for (int i=m_gr->indG[v]; i< m_gr->indG[v+1]; i++){ + int v2 = m_gr->nodeG[i]; + double w2 = m_gr->wG[i]; + if (Visited[v2]==1){ + continue; + } + double tentative=Distances[v]+w2; + + if (Visited1check[v2]==0 || Distances[v2] > tentative){ + Distances[v2]=tentative; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + v2] = distadd[(size_t)col * m_gr->nbnode + v] + m_gr->add[(size_t)col * total + i]; + } + if (Visited1check[v2]==0) touched.push_back(v2); + Visited1check[v2]=1; + + Q.push(std::make_pair(v2, tentative + sqrt(pow(m_gr->lat[v2]-lata,2)+pow(m_gr->lon[v2]-lona,2))/m_gr->k)); + + if (Visited2check[v2]==1){ + double temp=tentative + Distances2[v2]; + if (mu > temp){ + mu=temp; + for (int col = 0; col < n_add; col++){ + mu2[col] = distadd[(size_t)col * m_gr->nbnode + v] + m_gr->add[(size_t)col * total + i] + distadd2[(size_t)col * m_gr->nbnode + v2]; + } + } + } + + } + + } + } + + if (!Q.empty()){ + total1=Q.top().second; + } + + + } + } + + //Backward + else { + + int vv=Qr.top().first; + + Qr.pop(); + + if (Visited[vv]==0){ + + Visited[vv]=1; + + if ((Distances2[vv] + sqrt(pow(m_gr->lat[vv]-lata2,2)+pow(m_gr->lon[vv]-lona2,2))/m_gr->k) >= mu || (Distances2[vv] + total1 - sqrt(pow(m_gr->lat[vv]-lata,2)+pow(m_gr->lon[vv]-lona,2))/m_gr->k) >= mu){} + + else { + for (int i=m_gr->indGr[vv]; i< m_gr->indGr[vv+1]; i++){ + int vv2 = m_gr->nodeGr[i]; + double ww2 = m_gr->wGr[i]; + if (Visited[vv2]==1){ + continue; + } + double tentative=Distances2[vv]+ww2; + + if (Visited2check[vv2]==0 || Distances2[vv2] > tentative){ + Distances2[vv2]=tentative; + for (int col = 0; col < n_add; col++){ + distadd2[(size_t)col * m_gr->nbnode + vv2] = distadd2[(size_t)col * m_gr->nbnode + vv] + m_gr->addr[(size_t)col * totalr + i]; + } + if (Visited2check[vv2]==0) touched.push_back(vv2); + Visited2check[vv2]=1; + + Qr.push(std::make_pair(vv2, tentative + sqrt(pow(m_gr->lat[vv2]-lata2,2)+pow(m_gr->lon[vv2]-lona2,2))/m_gr->k)); + + if (Visited1check[vv2]==1){ + double temp=tentative + Distances[vv2]; + if (mu > temp){ + mu=temp; + for (int col = 0; col < n_add; col++){ + mu2[col] = distadd2[(size_t)col * m_gr->nbnode + vv] + m_gr->addr[(size_t)col * totalr + i] + distadd[(size_t)col * m_gr->nbnode + vv2]; + } + } + } + + } + + } + } + + if (!Qr.empty()){ + total2=Qr.top().second; + } + + } + } + + + } + + + for (int col = 0; col < n_add; col++){ + m_result[(size_t)col * n_pairs + k] = mu2[col]; + } + + for (int i = 0; i < (int)touched.size(); i++){ + int node = touched[i]; + Distances[node] = std::numeric_limits::max(); + Distances2[node] = std::numeric_limits::max(); + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + node] = std::numeric_limits::max(); + distadd2[(size_t)col * m_gr->nbnode + node] = std::numeric_limits::max(); + } + Visited[node] = 0; + Visited1check[node] = 0; + Visited2check[node] = 0; + } + touched.clear(); + + } + +} + + //////////////////////////////////////////////////////////////// Contracted graph // constructor @@ -656,3 +1182,178 @@ void distancePairC::bidirmod(std::size_t begin, std::size_t end){ } +//////////////////////////////////////////////////////////// Contracted graph, N aux columns + +// constructor +distancePairCAdd::distancePairCAdd(CGraph *gr, IVec dep, IVec arr) : + m_gr(gr), m_dep(dep), m_arr(arr) +{ + n_add = m_gr->n_add; + m_result.assign((size_t)m_dep.size() * n_add, 0.0); +} + +void distancePairCAdd::operator()(std::size_t begin, std::size_t end){ + bidirmod(begin, end); +} + +// modified bidirectional dijkstra, N aux columns summed inline (same +// structure as distancePairC::bidirmod, generalized from one running sum to +// n_add of them; add/addr already hold per-shortcut sums precomputed by +// aggCAdd, so no unpacking happens here) +void distancePairCAdd::bidirmod(std::size_t begin, std::size_t end){ + + int n_pairs = (int)m_dep.size(); + int total = (int)m_gr->nodeG.size(); // stride of m_gr->add + int totalr = (int)m_gr->nodeGr.size(); // stride of m_gr->addr + + std::vector Distances(m_gr->nbnode, std::numeric_limits::max()); + std::vector Distances2(m_gr->nbnode, std::numeric_limits::max()); + DVec distadd((size_t)m_gr->nbnode * n_add, std::numeric_limits::max()); + DVec distadd2((size_t)m_gr->nbnode * n_add, std::numeric_limits::max()); + DVec mu2(n_add, std::numeric_limits::max()); + + std::vector Visited1(m_gr->nbnode, 0); + std::vector Visited2(m_gr->nbnode, 0); + std::vector touched; + + double mu=std::numeric_limits::max(); + + for (std::size_t k=begin; k!=end;k++){ + + int StartNode = m_dep[k]; + int EndNode = m_arr[k]; + + Distances[StartNode] = 0.0; + Distances2[EndNode] = 0.0; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + StartNode] = 0.0; + distadd2[(size_t)col * m_gr->nbnode + EndNode] = 0.0; + } + touched.push_back(StartNode); + touched.push_back(EndNode); + + PQ Q; + PQ Qr; + Q.push(std::make_pair(StartNode, 0.0)); + Qr.push(std::make_pair(EndNode, 0.0)); + + + while (true) { + + if (Q.top().second > mu && Qr.top().second > mu){ + break; + } + if (Q.empty() && Qr.empty()){ + break; + } + + if (!Q.empty()){ + + int v=Q.top().first; + double w=Q.top().second; + Q.pop(); + + touched.push_back(v); + + Visited1[v] = 1; + + if ((Visited2[v]==1) && (Distances[v]+Distances2[v]) < mu){ + + mu=Distances[v]+Distances2[v]; + for (int col = 0; col < n_add; col++){ + mu2[col] = distadd[(size_t)col * m_gr->nbnode + v] + distadd2[(size_t)col * m_gr->nbnode + v]; + } + + } + + if (w <= Distances[v]) { + if (Stall_par(v, Distances, m_gr->nodeGr, m_gr->wGr, m_gr->indGr)) continue; + for (int i=m_gr->indG[v]; i< m_gr->indG[v+1]; i++){ + int v2 = m_gr->nodeG[i]; + double w2 = m_gr->wG[i]; + + if (Distances[v] + w2 < Distances[v2]) { + Distances[v2] = Distances[v] + w2; + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + v2] = distadd[(size_t)col * m_gr->nbnode + v] + m_gr->add[(size_t)col * total + i]; + } + Q.push(std::make_pair(v2, Distances[v2])); + Visited1[v2]=1; + + touched.push_back(v2); + + } + } + } + + } + + if (!Qr.empty()){ + int vv = Qr.top().first; + double ww = Qr.top().second; + Qr.pop(); + + touched.push_back(vv); + + Visited2[vv]=1; + + if ((Visited1[vv]== 1) && (Distances[vv]+Distances2[vv]) < mu){ + + mu=Distances[vv]+Distances2[vv]; + for (int col = 0; col < n_add; col++){ + mu2[col] = distadd[(size_t)col * m_gr->nbnode + vv] + distadd2[(size_t)col * m_gr->nbnode + vv]; + } + + } + + if (ww <= Distances2[vv]) { + if (Stall_par(vv, Distances2, m_gr->nodeG, m_gr->wG, m_gr->indG)) continue; + + for (int i=m_gr->indGr[vv]; i< m_gr->indGr[vv+1]; i++){ + int vv2 = m_gr->nodeGr[i]; + double ww2 = m_gr->wGr[i]; + + if (Distances2[vv] + ww2 < Distances2[vv2]) { + Distances2[vv2] = Distances2[vv] + ww2; + for (int col = 0; col < n_add; col++){ + distadd2[(size_t)col * m_gr->nbnode + vv2] = distadd2[(size_t)col * m_gr->nbnode + vv] + m_gr->addr[(size_t)col * totalr + i]; + } + Qr.push(std::make_pair(vv2, Distances2[vv2])); + Visited2[vv2]=1; + + touched.push_back(vv2); + + } + } + } + + } + + } + + + for (int col = 0; col < n_add; col++){ + m_result[(size_t)col * n_pairs + k] = mu2[col]; + } + + for (int i=0; i<(int)touched.size();i++) { + int node = touched[i]; + Distances[node] = std::numeric_limits::max(); + Distances2[node] = std::numeric_limits::max(); + for (int col = 0; col < n_add; col++){ + distadd[(size_t)col * m_gr->nbnode + node] = std::numeric_limits::max(); + distadd2[(size_t)col * m_gr->nbnode + node] = std::numeric_limits::max(); + } + Visited1[node] = 0; + Visited2[node] = 0; + } + + touched.clear(); + + mu=std::numeric_limits::max(); + std::fill(mu2.begin(), mu2.end(), std::numeric_limits::max()); + } + +} + + diff --git a/src/distance_pair.h b/src/distance_pair.h index 211ee06..db62155 100644 --- a/src/distance_pair.h +++ b/src/distance_pair.h @@ -34,6 +34,29 @@ struct distancePair : public Worker{ }; +// Distance-only pair query that sums N aux columns inline during the same +// search (no path reconstruction), generalizing distancePair's single-column +// add/addr mechanism to n_add columns at once. +struct distancePairAdd : public Worker{ + + Graph *m_gr; + IVec m_dep; + IVec m_arr; + int n_add; + int algorithm; // 0 : dijkstra with early stop, 1 : bidirectional dijkstra, 2 : A*, 3 : NBA* + DVec m_result; // flattened n_pairs x n_add, column-major + + distancePairAdd(Graph *gr, IVec dep, IVec arr, int algo); + + void dijkstra_early_stop(std::size_t begin, std::size_t end); + void bidir(std::size_t begin, std::size_t end); + void astar(std::size_t begin, std::size_t end); + void nba(std::size_t begin, std::size_t end); + + void operator()(std::size_t begin, std::size_t end); + +}; + struct distancePairC : public Worker{ // graph pointer to avoid copying input @@ -57,5 +80,25 @@ struct distancePairC : public Worker{ }; +// Same as distancePairC (contracted graph, modified bidirectional dijkstra), +// generalized to n_add aux columns summed inline. Aux values on m_gr->add / +// m_gr->addr are per-shortcut sums already precomputed by aggCAdd, so no +// unpacking happens during the query itself. +struct distancePairCAdd : public Worker{ + + CGraph *m_gr; + IVec m_dep; + IVec m_arr; + int n_add; + DVec m_result; // flattened n_pairs x n_add, column-major + + distancePairCAdd(CGraph *gr, IVec dep, IVec arr); + + void bidirmod(std::size_t begin, std::size_t end); + + void operator()(std::size_t begin, std::size_t end); + +}; + #endif diff --git a/src/graph.cpp b/src/graph.cpp index 9b02bab..3b15966 100644 --- a/src/graph.cpp +++ b/src/graph.cpp @@ -26,6 +26,7 @@ Graph::Graph(IVec &gfrom, IVec >o, DVec &gw, int nb, bool rm_dup){ nbedge = gfrom.size(); nbnode = nb; data = G(nbnode); + edge_origin = vector >(nbnode); add.resize(0); addr.resize(0); @@ -36,11 +37,15 @@ Graph::Graph(IVec &gfrom, IVec >o, DVec &gw, int nb, bool rm_dup){ exist = true; if (data[gfrom[i]][j].second > gw[i]){ data[gfrom[i]][j].second = gw[i]; + edge_origin[gfrom[i]][j] = IVec(1, i); } break; } } - if (!exist) data[gfrom[i]].push_back(make_pair(gto[i], gw[i])); + if (!exist) { + data[gfrom[i]].push_back(make_pair(gto[i], gw[i])); + edge_origin[gfrom[i]].push_back(IVec(1, i)); + } } } @@ -109,6 +114,74 @@ Graph::Graph(IVec &gfrom, IVec >o, DVec &gw, DVec &gadd, int nb){ } +Graph::Graph(IVec &gfrom, IVec >o, DVec &gw, DVec &gauxflat, int n_aux, int nb){ + nbedge = gfrom.size(); + nbnode = nb; + n_add = n_aux; + data = G(nbnode); + vector data2(nbnode); + vector data3(nbnode); // reversed + for (int i = 0; i != nbedge; ++i) { + data[gfrom[i]].push_back(std::make_pair(gto[i], gw[i])); + data2[gfrom[i]].push_back(i); + data3[gto[i]].push_back(i); + } + + // adj list + int count = 0; + for (int i = 0; i < nbnode; i++){ + count += data[i].size(); + } + int total = count; // == nbedge here (this constructor never merges edges) + + nodeG.resize(count); + wG.resize(count); + add.assign((size_t)total * n_aux, 0.0); + indG.resize(nbnode + 1); + + count=0; + + for (int i=0; i < data2.size();i++){ + indG[i] = count; + + for (int j = 0; j < data2[i].size(); j++){ + int orig_row = data2[i][j]; + nodeG[count] = gto[orig_row]; + wG[count] = gw[orig_row]; + for (int col = 0; col < n_aux; col++){ + add[(size_t)col * total + count] = gauxflat[(size_t)col * nbedge + orig_row]; + } + count += 1; + } + } + indG[nbnode]=count; + + // reversed + + nodeGr.resize(nodeG.size()); + wGr.resize(nodeG.size()); + addr.assign((size_t)total * n_aux, 0.0); + indGr.resize(nbnode + 1); + + count=0; + + for (int i=0; i < data3.size();i++){ + indGr[i] = count; + + for (int j = 0; j < data3[i].size(); j++){ + int orig_row = data3[i][j]; + nodeGr[count] = gfrom[orig_row]; + wGr[count] = gw[orig_row]; + for (int col = 0; col < n_aux; col++){ + addr[(size_t)col * total + count] = gauxflat[(size_t)col * nbedge + orig_row]; + } + count += 1; + } + } + indGr[nbnode]=count; + +} + Graph::Graph(IVec &gfrom, IVec >o, DVec &gw, DVec &gflow, DVec &gaux, DVec &gftt, DVec &galpha, DVec &gbeta, DVec &gcap, int nb){ @@ -220,6 +293,39 @@ Rcpp::List Graph::getEdges(){ return edges; } +Rcpp::List Graph::getEdgesWithOrigin(){ + nbedge = 0; + for (int i=0; i < data.size(); i++){ + nbedge += data[i].size(); + } + + IVec Newfrom(nbedge); + IVec Newto(nbedge); + DVec Neww(nbedge); + Rcpp::List origin(nbedge); + + int index=0; + for (int i=0; i < data.size();i++){ + for (int j=0; j < data[i].size();j++){ + + Newfrom[index]=i; + Newto[index]= data[i][j].first; + Neww[index]= data[i][j].second; + origin[index] = edge_origin[i][j]; + index+=1; + + } + } + + Rcpp::List edges(4); + edges[0] = Newfrom; + edges[1] = Newto; + edges[2] = Neww; + edges[3] = origin; + + return edges; +} + void Graph::to_adj_list(bool reversed){ int count = 0; diff --git a/src/graph.h b/src/graph.h index 6ab4e2f..3f40492 100644 --- a/src/graph.h +++ b/src/graph.h @@ -5,6 +5,8 @@ #include #include #include + +struct CCHPrepared; using namespace std; // Priority queue comparator @@ -46,6 +48,13 @@ class Graph{ SVec dict; G dataR; + // For each stored edge in `data`, the row indices into the original input + // edge table it represents. A fresh edge maps to itself; simplify() folds + // chains of edges into one and concatenates their origin lists so + // auxiliary edge weights can be summed once, after simplification, instead + // of being carried arithmetically through every intermediate merge. + vector > edge_origin; + // Parallel graph : 3 vectors (adjacency list) IVec nodeG; IVec indG; @@ -62,6 +71,7 @@ class Graph{ // additional weight DVec add; DVec addr; + int n_add = 1; // number of columns flattened into add/addr (multi-column aggregate_aux path only) /// reversed IVec nodeGr; @@ -82,6 +92,13 @@ class Graph{ // Constructor with additional weight Graph(IVec &gfrom, IVec >o, DVec &gw, DVec &gadd, int nb); + // Constructor with N additional weight columns (flattened, nbedge x n_aux, + // column-major). add/addr become flattened nbedge x n_aux matrices instead + // of a single column; used only by the multi-column aggregate_aux path so + // the single-column constructor above (and everything that depends on add + // being one value per edge) stays untouched. + Graph(IVec &gfrom, IVec >o, DVec &gw, DVec &gauxflat, int n_aux, int nb); + // constructor with full edges description Graph(IVec &gfrom, IVec >o, DVec &gw, DVec &gflow, DVec &gaux, DVec &gftt, @@ -96,6 +113,11 @@ class Graph{ // return edge description from current graph Rcpp::List getEdges(); + // same as getEdges(), plus the origin edge-index list for every surviving + // edge (see edge_origin above). Used by cpp_simplify() to sum auxiliary + // weights over the original edges once simplification is done. + Rcpp::List getEdgesWithOrigin(); + // routing algorithms, grouped by return type DVec routing_dvec(IVec dep, IVec arr, int algo); @@ -133,7 +155,8 @@ class Graph{ void assign_traffic(int method, double max_gap, int max_iter, IVec dep, IVec arr, DVec dem, int mode, - bool contract, bool phast, bool verbose); + bool contract, bool phast, bool verbose, + CCHPrepared* prepared_cch = NULL); void algorithmB(int batch_size, int n_batch, string file_path, double max_gap, int max_iter, IVec dep, IVec arr, DVec dem, int mode, int inner_iter, double NUM_TOL, diff --git a/src/simplify.cpp b/src/simplify.cpp index b5d461c..3fbd7f3 100644 --- a/src/simplify.cpp +++ b/src/simplify.cpp @@ -16,6 +16,12 @@ void quickDelete2( int idx ,vector &vec) vec.pop_back(); } +void quickDeleteVec( int idx ,vector &vec) +{ + vec[idx] = vec.back(); + vec.pop_back(); +} + void Graph::simp(vector &Gr, IVec &Junction, vector &Edges, @@ -234,6 +240,8 @@ void Graph::simp(vector &Gr, DVec Acc(Edges.size()); DVec Acc2(Edges.size()); + vector AccOrigin(Edges.size()); + vector AccOrigin2(Edges.size()); for (int i=0; i < Edges.size(); ++i){ IVec path=Edges[i]; double w=0.0; @@ -246,6 +254,7 @@ void Graph::simp(vector &Gr, for (int m=0; m < data[node1].size(); ++m){ if (data[node1][m].first==node2){ w+=data[node1][m].second; + AccOrigin[i].insert(AccOrigin[i].end(), edge_origin[node1][m].begin(), edge_origin[node1][m].end()); } } @@ -266,6 +275,7 @@ void Graph::simp(vector &Gr, for (int m=0; m < data[node1].size(); ++m){ if (data[node1][m].first==node2){ w2+=data[node1][m].second; + AccOrigin2[i].insert(AccOrigin2[i].end(), edge_origin[node1][m].begin(), edge_origin[node1][m].end()); } } @@ -284,9 +294,11 @@ void Graph::simp(vector &Gr, vector > Gr_remov=data[Edges[i][j]]; data[Edges[i][j]].erase(data[Edges[i][j]].begin(),data[Edges[i][j]].end()); + edge_origin[Edges[i][j]].erase(edge_origin[Edges[i][j]].begin(),edge_origin[Edges[i][j]].end()); for (int k=0; k < Gr[Edges[i][j]].size();k++){ for (int m=0; m < data[Gr[Edges[i][j]][k]].size();m++){ if (data[Gr[Edges[i][j]][k]][m].first==Edges[i][j]) { + quickDeleteVec(m,edge_origin[Gr[Edges[i][j]][k]]); quickDelete(m,data[Gr[Edges[i][j]][k]]); } } @@ -314,9 +326,11 @@ void Graph::simp(vector &Gr, vector > Gr_remov=data[i]; data[i].erase(data[i].begin(),data[i].end()); + edge_origin[i].erase(edge_origin[i].begin(),edge_origin[i].end()); for (int k=0; k < Gr[i].size();k++){ for (int m=0; m < data[Gr[i][k]].size();m++){ if (data[Gr[i][k]][m].first==i) { + quickDeleteVec(m,edge_origin[Gr[i][k]]); quickDelete(m,data[Gr[i][k]]); } } @@ -355,6 +369,7 @@ void Graph::simp(vector &Gr, dup=true; if (data[node1][k].second >= Acc[i]){ data[node1][k].second=Acc[i]; + edge_origin[node1][k]=AccOrigin[i]; break; } @@ -364,6 +379,7 @@ void Graph::simp(vector &Gr, if (dup==false) { data[node1].push_back(make_pair(node2,Acc[i])); + edge_origin[node1].push_back(AccOrigin[i]); Gr[node2].push_back(node1); } @@ -375,6 +391,7 @@ void Graph::simp(vector &Gr, dup=true; if (data[node2][k].second >= Acc2[i]){ data[node2][k].second=Acc2[i]; + edge_origin[node2][k]=AccOrigin2[i]; break; } @@ -384,6 +401,7 @@ void Graph::simp(vector &Gr, if (dup==false) { data[node2].push_back(make_pair(node1,Acc2[i])); + edge_origin[node2].push_back(AccOrigin2[i]); Gr[node1].push_back(node2); } } diff --git a/src/traffic.cpp b/src/traffic.cpp index 7aa34fb..d4ca6b2 100644 --- a/src/traffic.cpp +++ b/src/traffic.cpp @@ -1,5 +1,6 @@ #include "graph.h" #include "cgraph.h" +#include "cch.h" #include "aon.h" #include "bush.h" #include @@ -10,6 +11,14 @@ using namespace std; using namespace RcppParallel; +namespace { + +const int AON_DIJKSTRA_GROUPED = 0; +const int AON_BIDIRECTIONAL = 2; +const int AON_NBA = 3; +const int AON_CCH_ELIMINATION = 7; + +} // namespace double Graph::bissection(double tol){ double inf = 0.0; @@ -226,26 +235,126 @@ DVec ch_aon(Graph* ptr, IVec dep, IVec arr, DVec &dem, bool phast, int algo){ } +void graph_adjacency_edges(Graph* ptr, IVec &gfrom, IVec >o){ + gfrom.resize(ptr->nbedge); + gto.resize(ptr->nbedge); + + int count = 0; + for (int i = 0; i < ptr->indG.size() - 1; i++){ + for (int j = ptr->indG[i]; j < ptr->indG[i + 1]; j++){ + gfrom[count] = i; + gto[count] = ptr->nodeG[j]; + count++; + } + } +} + +struct CCHAssignmentScratch { + DVec forward; + DVec backward; + IVec forward_first_arc; + IVec forward_first_dir; + IVec forward_second_arc; + IVec forward_second_dir; + IVec forward_original; + IVec backward_first_arc; + IVec backward_first_dir; + IVec backward_second_arc; + IVec backward_second_dir; + IVec backward_original; +}; + +DVec cch_aon_current(Graph* ptr, + CCHPrepared &prepared, + IVec &gfrom, + IVec >o, + IVec &dep, + IVec &arr, + DVec &dem, + IVec &input_arc, + IVec &input_forward, + CCHAssignmentScratch &scratch){ + + customize_cch(gfrom, gto, ptr->wG, ptr->nbnode, prepared.rank, prepared.tail, prepared.head, + prepared.rank_first_out, prepared.rank_adj_head, prepared.rank_adj_arc, + input_arc, input_forward, scratch.forward, scratch.backward, + scratch.forward_first_arc, scratch.forward_first_dir, scratch.forward_second_arc, scratch.forward_second_dir, scratch.forward_original, + scratch.backward_first_arc, scratch.backward_first_dir, scratch.backward_second_arc, scratch.backward_second_dir, scratch.backward_original); + + return aon_flow_cch_elimination_grouped(gfrom, gto, ptr->wG, ptr->nbnode, + prepared.rank, prepared.first_out, prepared.adj_head, prepared.adj_arc, + prepared.elimination_tree_parent, scratch.forward, scratch.backward, + scratch.forward_first_arc, scratch.forward_first_dir, scratch.forward_second_arc, scratch.forward_second_dir, scratch.forward_original, + scratch.backward_first_arc, scratch.backward_first_dir, scratch.backward_second_arc, scratch.backward_second_dir, scratch.backward_original, + dep, arr, dem); +} + +DVec assignment_aon(Graph* ptr, + aonGraph &gr, + IVec &dep, + IVec &arr, + DVec &dem, + int mode, + bool contract, + bool phast, + bool use_cch, + CCHPrepared &prepared, + IVec &input_arc, + IVec &input_forward, + IVec &cch_from, + IVec &cch_to, + CCHAssignmentScratch &cch_scratch){ + if (use_cch) { + return cch_aon_current(ptr, prepared, cch_from, cch_to, dep, arr, dem, input_arc, input_forward, cch_scratch); + } + + if (contract) return ch_aon(ptr, dep, arr, dem, phast, mode); + + parallelReduce(0, gr.n_iter, gr); + DVec result = gr.m_result; + fill(gr.m_result.begin(), gr.m_result.end(), 0.0); + return result; +} + void Graph::assign_traffic(int method, double max_gap, int max_iter, IVec dep, IVec arr, DVec dem, int mode, bool contract, bool phast, - bool verbose){ // 0 : msa, 1 : fw, 2 : cfw, 3 : bfw + bool verbose, + CCHPrepared* prepared_cch){ // 0 : msa, 1 : fw, 2 : cfw, 3 : bfw setReverse(); to_adj_list(true); - aonGraph gr(this, dep, arr, dem, mode); + bool use_cch = mode == AON_CCH_ELIMINATION; + int graph_mode = use_cch ? AON_BIDIRECTIONAL : mode; + aonGraph gr(this, dep, arr, dem, graph_mode); + IVec cch_from; + IVec cch_to; + CCHPrepared cch_prepared; + CCHPrepared* active_cch = &cch_prepared; + IVec* cch_input_arc = &cch_prepared.input_arc; + IVec* cch_input_forward = &cch_prepared.input_forward; + IVec remapped_input_arc; + IVec remapped_input_forward; + CCHAssignmentScratch cch_scratch; + + if (use_cch){ + graph_adjacency_edges(this, cch_from, cch_to); + if (prepared_cch == NULL){ + IVec empty_order; + cch_prepared = build_cch(cch_from, cch_to, nbnode, empty_order); + } else { + active_cch = prepared_cch; + remap_cch_input_arcs(cch_from, cch_to, *active_cch, remapped_input_arc, remapped_input_forward); + cch_input_arc = &remapped_input_arc; + cch_input_forward = &remapped_input_forward; + } + } // MSA if (method == 0) { gap = numeric_limits::max(); - if (contract) { - aux = ch_aon(this, dep, arr, dem, phast, mode); - } else{ - parallelReduce(0, gr.n_iter, gr); - aux = gr.m_result; - fill(gr.m_result.begin(), gr.m_result.end(), 0.0); - } + aux = assignment_aon(this, gr, dep, arr, dem, graph_mode, contract, phast, use_cch, *active_cch, *cch_input_arc, *cch_input_forward, cch_from, cch_to, cch_scratch); it = 1; while (gap > max_gap){ @@ -256,13 +365,7 @@ void Graph::assign_traffic(int method, double max_gap, int max_iter, update_flow(theta); update_cost(); - if (contract) { - aux = ch_aon(this, dep, arr, dem, phast, mode); - } else{ - parallelReduce(0, gr.n_iter, gr); - aux = gr.m_result; - fill(gr.m_result.begin(), gr.m_result.end(), 0.0); - } + aux = assignment_aon(this, gr, dep, arr, dem, graph_mode, contract, phast, use_cch, *active_cch, *cch_input_arc, *cch_input_forward, cch_from, cch_to, cch_scratch); update_tstt(); update_sptt(); @@ -280,13 +383,7 @@ void Graph::assign_traffic(int method, double max_gap, int max_iter, { gap = numeric_limits::max(); - if (contract) { - aux = ch_aon(this, dep, arr, dem, phast, mode); - } else{ - parallelReduce(0, gr.n_iter, gr); - aux = gr.m_result; - fill(gr.m_result.begin(), gr.m_result.end(), 0.0); - } + aux = assignment_aon(this, gr, dep, arr, dem, graph_mode, contract, phast, use_cch, *active_cch, *cch_input_arc, *cch_input_forward, cch_from, cch_to, cch_scratch); it = 1; double theta = 0.0; @@ -300,13 +397,7 @@ void Graph::assign_traffic(int method, double max_gap, int max_iter, update_flow(theta); update_cost(); - if (contract) { - aux = ch_aon(this, dep, arr, dem, phast, mode); - } else{ - parallelReduce(0, gr.n_iter, gr); - aux = gr.m_result; - fill(gr.m_result.begin(), gr.m_result.end(), 0.0); - } + aux = assignment_aon(this, gr, dep, arr, dem, graph_mode, contract, phast, use_cch, *active_cch, *cch_input_arc, *cch_input_forward, cch_from, cch_to, cch_scratch); update_tstt(); update_sptt(); @@ -323,13 +414,7 @@ void Graph::assign_traffic(int method, double max_gap, int max_iter, { gap = numeric_limits::max(); - if (contract) { - aux = ch_aon(this, dep, arr, dem, phast, mode); - } else{ - parallelReduce(0, gr.n_iter, gr); - aux = gr.m_result; - fill(gr.m_result.begin(), gr.m_result.end(), 0.0); - } + aux = assignment_aon(this, gr, dep, arr, dem, graph_mode, contract, phast, use_cch, *active_cch, *cch_input_arc, *cch_input_forward, cch_from, cch_to, cch_scratch); it = 1; double theta = 0.0; @@ -349,13 +434,7 @@ void Graph::assign_traffic(int method, double max_gap, int max_iter, update_cost(); aux_prev = aux; - if (contract) { - aux = ch_aon(this, dep, arr, dem, phast, mode); - } else{ - parallelReduce(0, gr.n_iter, gr); - aux = gr.m_result; - fill(gr.m_result.begin(), gr.m_result.end(), 0.0); - } + aux = assignment_aon(this, gr, dep, arr, dem, graph_mode, contract, phast, use_cch, *active_cch, *cch_input_arc, *cch_input_forward, cch_from, cch_to, cch_scratch); update_tstt(); update_sptt(); @@ -373,13 +452,7 @@ void Graph::assign_traffic(int method, double max_gap, int max_iter, { gap = numeric_limits::max(); - if (contract) { - aux = ch_aon(this, dep, arr, dem, phast, mode); - } else{ - parallelReduce(0, gr.n_iter, gr); - aux = gr.m_result; - fill(gr.m_result.begin(), gr.m_result.end(), 0.0); - } + aux = assignment_aon(this, gr, dep, arr, dem, graph_mode, contract, phast, use_cch, *active_cch, *cch_input_arc, *cch_input_forward, cch_from, cch_to, cch_scratch); it = 1; @@ -403,13 +476,7 @@ void Graph::assign_traffic(int method, double max_gap, int max_iter, aux_prev2 = aux_prev; aux_prev = aux; - if (contract) { - aux = ch_aon(this, dep, arr, dem, phast, mode); - } else{ - parallelReduce(0, gr.n_iter, gr); - aux = gr.m_result; - fill(gr.m_result.begin(), gr.m_result.end(), 0.0); - } + aux = assignment_aon(this, gr, dep, arr, dem, graph_mode, contract, phast, use_cch, *active_cch, *cch_input_arc, *cch_input_forward, cch_from, cch_to, cch_scratch); update_tstt(); update_sptt(); diff --git a/src/wrapper.cpp b/src/wrapper.cpp index 731b565..0dbd1ff 100644 --- a/src/wrapper.cpp +++ b/src/wrapper.cpp @@ -1,7 +1,10 @@ #include "graph.h" #include "cgraph.h" +#include "cch.h" #include "aon.h" #include "aggc.h" +#include "distance_pair.h" +#include "distance_mat.h" #include #include #include @@ -162,14 +165,37 @@ Rcpp::List cpppathmat(std::vector &gfrom, std::vector >o, std::vecto // [[Rcpp::export]] Rcpp::List cppsimplify(std::vector &gfrom, std::vector >o, std::vector &gw, int nb, - std::vector keep, bool rm_loop, bool iterate, bool progress){ + std::vector keep, bool rm_loop, bool iterate, bool progress, + Rcpp::NumericMatrix aux){ Graph network(gfrom, gto, gw, nb, true); network.simplify(rm_loop, keep, iterate, progress); - Rcpp::List res = network.getEdges(); + Rcpp::List res = network.getEdgesWithOrigin(); + + Rcpp::List result(4); + result[0] = res[0]; + result[1] = res[1]; + result[2] = res[2]; + + int n_aux = aux.ncol(); + if (n_aux > 0){ + Rcpp::List origin = res[3]; + Rcpp::NumericMatrix summed_aux(origin.size(), n_aux); + for (int i = 0; i < origin.size(); i++){ + std::vector rows = origin[i]; + for (int col = 0; col < n_aux; col++){ + double total = 0.0; + for (int r = 0; r < rows.size(); r++) total += aux(rows[r], col); + summed_aux(i, col) = total; + } + } + result[3] = summed_aux; + } else { + result[3] = Rcpp::NumericMatrix(0, 0); + } - return res; + return result; } @@ -194,6 +220,192 @@ Rcpp::List cppcontract(std::vector &gfrom,std::vector >o,std::vector +// [[Rcpp::export]] +Rcpp::List cppcchprepare(std::vector &gfrom, std::vector >o, int NbNodes, std::vector order, bool display_progress = false){ + CCHPrepared prepared = build_cch(gfrom, gto, NbNodes, order, display_progress); + + Rcpp::List result(12); + result[0] = prepared.rank; + result[1] = prepared.tail; + result[2] = prepared.head; + result[3] = prepared.first_out; + result[4] = prepared.adj_head; + result[5] = prepared.adj_arc; + result[6] = prepared.rank_first_out; + result[7] = prepared.rank_adj_head; + result[8] = prepared.rank_adj_arc; + result[9] = prepared.input_arc; + result[10] = prepared.input_forward; + result[11] = prepared.elimination_tree_parent; + result.attr("names") = Rcpp::CharacterVector::create("rank", "tail", "head", "first_out", "adj_head", "adj_arc", + "rank_first_out", "rank_adj_head", "rank_adj_arc", + "input_arc", "input_forward", "elimination_tree_parent"); + return result; +} + +// [[Rcpp::export]] +Rcpp::List cppcchcustomize(std::vector &gfrom, std::vector >o, std::vector &gw, int NbNodes, + std::vector &rank, std::vector &tail, std::vector &head, + std::vector &rank_first_out, std::vector &rank_adj_head, std::vector &rank_adj_arc, + std::vector &input_arc, std::vector &input_forward){ + DVec forward; + DVec backward; + IVec forward_first_arc; + IVec forward_first_dir; + IVec forward_second_arc; + IVec forward_second_dir; + IVec forward_original; + IVec backward_first_arc; + IVec backward_first_dir; + IVec backward_second_arc; + IVec backward_second_dir; + IVec backward_original; + customize_cch(gfrom, gto, gw, NbNodes, rank, tail, head, rank_first_out, rank_adj_head, rank_adj_arc, + input_arc, input_forward, forward, backward, + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original); + + Rcpp::List result(12); + result[0] = forward; + result[1] = backward; + result[2] = forward_first_arc; + result[3] = forward_first_dir; + result[4] = forward_second_arc; + result[5] = forward_second_dir; + result[6] = forward_original; + result[7] = backward_first_arc; + result[8] = backward_first_dir; + result[9] = backward_second_arc; + result[10] = backward_second_dir; + result[11] = backward_original; + result.attr("names") = Rcpp::CharacterVector::create("forward", "backward", + "forward_first_arc", "forward_first_dir", + "forward_second_arc", "forward_second_dir", + "forward_original", + "backward_first_arc", "backward_first_dir", + "backward_second_arc", "backward_second_dir", + "backward_original"); + return result; +} + +// [[Rcpp::export]] +Rcpp::NumericVector cppdistcch(int NbNodes, + std::vector &rank, + std::vector &first_out, std::vector &adj_head, std::vector &adj_arc, + std::vector &elimination_tree_parent, + std::vector &forward, std::vector &backward, + std::vector dep, std::vector arr){ + Rcpp::NumericVector result = Rcpp::wrap(distance_pair_cch(NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, dep, arr)); + check_nas_vec(result); + return result; +} + +// [[Rcpp::export]] +Rcpp::NumericMatrix cppdistmatcch(int NbNodes, + std::vector &first_out, std::vector &adj_head, std::vector &adj_arc, + std::vector &forward, std::vector &backward, + std::vector dep, std::vector arr){ + Rcpp::NumericMatrix result = distance_matrix_cch(NbNodes, first_out, adj_head, adj_arc, forward, backward, dep, arr); + check_nas_mat(result); + return result; +} + +// [[Rcpp::export]] +Rcpp::NumericMatrix cpppathvaluescch(std::vector &gfrom, + int NbNodes, + std::vector &rank, + std::vector &first_out, + std::vector &adj_head, + std::vector &adj_arc, + std::vector &elimination_tree_parent, + std::vector &forward, + std::vector &backward, + std::vector &forward_first_arc, + std::vector &forward_first_dir, + std::vector &forward_second_arc, + std::vector &forward_second_dir, + std::vector &forward_original, + std::vector &backward_first_arc, + std::vector &backward_first_dir, + std::vector &backward_second_arc, + std::vector &backward_second_dir, + std::vector &backward_original, + std::vector dep, + std::vector arr, + Rcpp::NumericMatrix values){ + std::vector flat_values(values.begin(), values.end()); + Rcpp::NumericMatrix result = path_values_pair_cch(gfrom, NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, + forward, backward, + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + dep, arr, flat_values, values.ncol()); + check_nas_mat(result); + return result; +} + +// [[Rcpp::export]] +Rcpp::List cpppathcch(int NbNodes, + std::vector &rank, + std::vector &gfrom, + std::vector >o, + std::vector &first_out, + std::vector &adj_head, + std::vector &adj_arc, + std::vector &elimination_tree_parent, + std::vector &forward, + std::vector &backward, + std::vector &forward_first_arc, + std::vector &forward_first_dir, + std::vector &forward_second_arc, + std::vector &forward_second_dir, + std::vector &forward_original, + std::vector &backward_first_arc, + std::vector &backward_first_dir, + std::vector &backward_second_arc, + std::vector &backward_second_dir, + std::vector &backward_original, + std::vector dict, + std::vector keep, + std::vector dep, + std::vector arr){ + vector result = path_pair_cch(NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, forward, backward, + gfrom, gto, + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + dep, arr, keep); + vector dict_str(dict.begin(), dict.end()); + vector result_str(result.size()); + for (size_t i = 0; i < result.size(); i++){ + result_str[i].reserve(result[i].size()); + for (size_t j = 0; j < result[i].size(); j++){ + if (result[i][j] >= 0 && result[i][j] < (int)dict_str.size()){ + result_str[i].push_back(dict_str[result[i][j]]); + } + } + } + return Rcpp::wrap(result_str); +} + +// [[Rcpp::export]] +Rcpp::List cppaoncchelim(std::vector &gfrom, std::vector >o, std::vector &gw, int NbNodes, + std::vector &rank, + std::vector &first_out, std::vector &adj_head, std::vector &adj_arc, + std::vector &elimination_tree_parent, + std::vector &forward, std::vector &backward, + std::vector &forward_first_arc, std::vector &forward_first_dir, + std::vector &forward_second_arc, std::vector &forward_second_dir, + std::vector &forward_original, + std::vector &backward_first_arc, std::vector &backward_first_dir, + std::vector &backward_second_arc, std::vector &backward_second_dir, + std::vector &backward_original, + std::vector dep, std::vector arr, std::vector dem){ + return aon_cch_elimination_grouped(gfrom, gto, gw, NbNodes, rank, first_out, adj_head, adj_arc, elimination_tree_parent, + forward, backward, + forward_first_arc, forward_first_dir, forward_second_arc, forward_second_dir, forward_original, + backward_first_arc, backward_first_dir, backward_second_arc, backward_second_dir, backward_original, + dep, arr, dem); +} + // [[Rcpp::export]] Rcpp::NumericVector cppdistC(std::vector &gfrom, std::vector >o, std::vector &gw, int nb, @@ -287,6 +499,23 @@ Rcpp::NumericMatrix cpppadd(std::vector &gfrom, std::vector >o, std: return result; } +// [[Rcpp::export]] +Rcpp::NumericVector cppmataddmulti(std::vector &gfrom, std::vector >o, std::vector &gw, + Rcpp::NumericMatrix aux, int nb, + std::vector dep, std::vector arr){ + int n_aux = aux.ncol(); + std::vector aux_flat(aux.begin(), aux.end()); + Graph network(gfrom, gto, gw, aux_flat, n_aux, nb); + + distanceMatAdd dijfunc(&network, dep, arr); + parallelFor(0, dep.size(), dijfunc, 1, 12); + + Rcpp::NumericVector result = Rcpp::wrap(dijfunc.m_result); + result.attr("dim") = Rcpp::IntegerVector::create((int)dep.size(), (int)arr.size(), n_aux); + check_nas_vec(result); + return result; +} + // [[Rcpp::export]] @@ -310,6 +539,34 @@ Rcpp::NumericVector cppdistadd(std::vector &gfrom, std::vector >o, s +// [[Rcpp::export]] + +Rcpp::NumericMatrix cppdistaddmulti(std::vector &gfrom, std::vector >o, std::vector &gw, + Rcpp::NumericMatrix aux, int nb, + std::vector lat, std::vector lon, double k, + std::vector dep, std::vector arr, int algo){ + int n_aux = aux.ncol(); + std::vector aux_flat(aux.begin(), aux.end()); + Graph network(gfrom, gto, gw, aux_flat, n_aux, nb); + + if (algo == 2 || algo == 3){ + network.setLatLon(lat, lon); + network.k = k; + } + + distancePairAdd dijfunc(&network, dep, arr, algo); + parallelFor(0, dep.size(), dijfunc, 1, 12); + + Rcpp::NumericMatrix result(dep.size(), n_aux); + for (int col = 0; col < n_aux; col++){ + for (int row = 0; row < (int)dep.size(); row++){ + result(row, col) = dijfunc.m_result[(size_t)col * dep.size() + row]; + } + } + check_nas_mat(result); + return result; +} + // [[Rcpp::export]] Rcpp::NumericVector cppdistaddC(std::vector &orfrom, std::vector &orto, std::vector &orw, std::vector &gadd, @@ -336,6 +593,42 @@ Rcpp::NumericVector cppdistaddC(std::vector &orfrom, std::vector &orto return result; } +// [[Rcpp::export]] + +Rcpp::NumericMatrix cppdistaddCmulti(std::vector &orfrom, std::vector &orto, std::vector &orw, Rcpp::NumericMatrix aux, + std::vector &gfrom, std::vector >o, std::vector &gw, int nb, + std::vector &rank, std::vector &shortf, std::vector &shortt, std::vector &shortc, + bool phast, + std::vector dep, std::vector arr){ + int n_aux = aux.ncol(); + std::vector aux_flat(aux.begin(), aux.end()); + Graph network(orfrom, orto, orw, aux_flat, n_aux, nb); + + CGraph networkc(gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast); + networkc.n_add = n_aux; + networkc.construct_shortcuts(); + networkc.to_adj_list(false, phast); + networkc.to_adj_list(true, phast); + + aggCAdd dijfunc(&networkc, &network); + parallelFor(0, networkc.nbnode, dijfunc); + + networkc.add = dijfunc.m_result; + networkc.addr = dijfunc.m_result2; + + distancePairCAdd distfunc(&networkc, dep, arr); + parallelFor(0, dep.size(), distfunc, 1, 12); + + Rcpp::NumericMatrix result(dep.size(), n_aux); + for (int col = 0; col < n_aux; col++){ + for (int row = 0; row < (int)dep.size(); row++){ + result(row, col) = distfunc.m_result[(size_t)col * dep.size() + row]; + } + } + check_nas_mat(result); + return result; +} + // [[Rcpp::export]] @@ -363,6 +656,38 @@ Rcpp::NumericMatrix cppaddC(std::vector &orfrom, std::vector &orto, st return result; } +// [[Rcpp::export]] + +Rcpp::NumericVector cppmataddCmulti(std::vector &orfrom, std::vector &orto, std::vector &orw, Rcpp::NumericMatrix aux, + std::vector &gfrom, std::vector >o, std::vector &gw, int nb, + std::vector &rank, std::vector &shortf, std::vector &shortt, std::vector &shortc, + bool phast, + std::vector dep, std::vector arr){ + int n_aux = aux.ncol(); + std::vector aux_flat(aux.begin(), aux.end()); + Graph network(orfrom, orto, orw, aux_flat, n_aux, nb); + + CGraph networkc(gfrom, gto, gw, nb, rank, shortf, shortt, shortc, phast); + networkc.n_add = n_aux; + networkc.construct_shortcuts(); + networkc.to_adj_list(false, phast); + networkc.to_adj_list(true, phast); + + aggCAdd dijfunc(&networkc, &network); + parallelFor(0, networkc.nbnode, dijfunc); + + networkc.add = dijfunc.m_result; + networkc.addr = dijfunc.m_result2; + + phastCAdd distfunc(&networkc, dep, arr); + parallelFor(0, dep.size(), distfunc, 1, 12); + + Rcpp::NumericVector result = Rcpp::wrap(distfunc.m_result); + result.attr("dim") = Rcpp::IntegerVector::create((int)dep.size(), (int)arr.size(), n_aux); + check_nas_vec(result); + return result; +} + // [[Rcpp::export]] @@ -449,6 +774,62 @@ Rcpp::List cppaonC(std::vector &orfrom, std::vector &orto, std::vector } +// [[Rcpp::export]] + +Rcpp::List cpptrafficcch(std::vector &gfrom, std::vector >o, std::vector &gw, + std::vector &gflow, std::vector &gaux, std::vector &gftt, + std::vector &galpha, std::vector &gbeta, std::vector &gcap, + int nb, + std::vector dep, std::vector arr, std::vector dem, + double max_gap, int max_it, int method, + std::vector &rank, + std::vector &tail, + std::vector &head, + std::vector &first_out, + std::vector &adj_head, + std::vector &adj_arc, + std::vector &rank_first_out, + std::vector &rank_adj_head, + std::vector &rank_adj_arc, + std::vector &input_arc, + std::vector &input_forward, + std::vector &elimination_tree_parent, + bool verbose){ + + Graph network(gfrom,gto, gw, gflow,gaux, gftt, galpha, gbeta, gcap, nb); + + CCHPrepared prepared; + prepared.nbnode = nb; + prepared.rank = rank; + prepared.tail = tail; + prepared.head = head; + prepared.first_out = first_out; + prepared.adj_head = adj_head; + prepared.adj_arc = adj_arc; + prepared.rank_first_out = rank_first_out; + prepared.rank_adj_head = rank_adj_head; + prepared.rank_adj_arc = rank_adj_arc; + prepared.input_arc = input_arc; + prepared.input_forward = input_forward; + prepared.elimination_tree_parent = elimination_tree_parent; + + network.assign_traffic(method, max_gap, max_it, dep, arr, dem, 7, false, false, verbose, &prepared); + + Rcpp::List result(10); + result[0] = network.indG2; + result[1] = network.nodeG; + result[2] = network.ftt; + result[3] = network.wG; + result[4] = network.flow; + result[5] = network.cap; + result[6] = network.alpha; + result[7] = network.beta; + result[8] = network.gap; + result[9] = network.it; + return (result); + +} + // [[Rcpp::export]] Rcpp::List cpptraffic(std::vector &gfrom, std::vector >o, std::vector &gw,