Merge pull request #32 from vlarmet/feat/cch

feat(CCH) : support for customizable contraction hierarchies

Author: vlarmet

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) <doi:10.1016/0191-2615(84)90029-8>),
     Conjugate and Bi-Conjugate Frank-Wolfe algorithms (M. Mitradjieva, P. O. Lindberg (2012) <doi:10.1287/trsc.1120.0409>),
-    Algorithm-B (R. B. Dial (2006) <doi:10.1016/j.trb.2006.02.008>).
+    Algorithm-B (R. B. Dial (2006) <doi:10.1016/j.trb.2006.02.008>),
+    Customizable Contraction Hierarchies (J. Dibbelt, B. Strasser, D. Wagner (2014) <doi:10.1145/2886843>).
 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

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")

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
 ===============
 

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)
 }

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)
+    }
 
   }