Space-Efficient, Fast and Exact Routing in Time-Dependent Road Networks ^†

Abstract

We study the problem of quickly computing point-to-point shortest paths in massive road networks with traffic predictions. Incorporating traffic predictions into routing allows, for example, to avoid commuter traffic congestions. Existing techniques follow a two-phase approach: In a preprocessing step, an index is built. The index depends on the road network and the traffic patterns but not on the path start and end. The latter are the input of the query phase, in which shortest paths are computed. All existing techniques have large index size, slow query running times or may compute suboptimal paths. In this work, we introduce CATCHUp (Customizable Approximated Time-dependent Contraction Hierarchies through Unpacking), the first algorithm that simultaneously achieves all three objectives. The core idea of CATCHUp is to store paths instead of travel times at shortcuts. Shortcut travel times are derived lazily from the stored paths. We perform an experimental study on a set of real world instances and compare our approach with state-of-the-art techniques. Our approach achieves the fastest preprocessing, competitive query running times and up to 38 times smaller indexes than competing approaches.

1. Introduction

Routing in road networks is a well-studied topic with a plethora of real world applications. Services such as Google, Baidu, Yandex, Bing, Apple or HERE Maps are ubiquitous and used by millions of users on a daily basis. The core problem is to compute the fastest route between a source and a target. The idealized problem can be formalized as the classic point-to-point shortest path problem. Streets are modeled as arcs. Street intersections are modeled as nodes. Travel times are modeled as scalar arc weights. Unfortunately, this idealized view does not model certain important real world effects. An example are recurring commuter congestions. In this article, we consider an extended problem in which travel times are time-dependent. The travel time of an arc is a function of the moment where a car enters the arc. Figure 1 depicts an example.

Computing shortest-paths using Dijkstra’s [1] algorithm is possible both in the classical and in the time-dependent setting. However, for many applications, its running time is too large. To achieve fast running times, a two-phase approach is used. In the first phase, the preprocessing phase, an index is constructed. The index only depends on the road networks and the arc travel times. In the second phase, the query phase, shortest paths are computed utilizing this index.

An important ingredient for many two-phase techniques [2,3,4,5,6] are shortcuts. Shortcuts are additional arcs introduced during preprocessing, which bypass parts of the input graph such as in Figure 2a. The weight of a shortcut is set to the length of the shortest path between its endpoints. When computing shortest paths, only few shortcuts are explored instead of many arcs in the input graph. The path represented by a shortcut can be obtained lazily, for example by running local Dijkstra searches [7] or iterating over possible middle nodes when shortcuts always represent two other (shortcut) arcs [5,6].

This approach has been extended to the time-dependent setting [8,9]. Shortcuts are no longer associated with scalar weights. Instead, travel time functions are used that map the entry time into a shortcut to the travel time through it. Typically, these functions are represented as piecewise linear functions. They are stored as a sequence of breakpoints. Unfortunately, these functions can become very complex. Computing and storing them is expensive. The number of breakpoints in a shortcut’s function practically corresponds to the accumulated number of breakpoints of the functions of bypassed arcs. Contrary to the classic setting, shortcuts aggregate the complexity of paths they represent, rather than skipping it. This leads to slow preprocessing and prohibitive memory consumption.

In this paper, we explore an alternative approach to shortcut travel time functions. Rather than explicitly storing them and obtaining paths lazily, we store paths and obtain travel times lazily. We expect that the shortest path between two nodes changes less frequently than the travel time. Intuitively, going via a highway may be slower due to congestion but is usually still the fastest option. Consider the functions f and g in Figure 2b. These functions are travel time functions of two paths between the same endpoints and have many breakpoints. If we want to store the travel time function of a shortcut between these endpoints, we need to store the function $h=min(f,g)$. Storing h explicitly requires roughly a number of breakpoints proportional to the number of breakpoints in f and g. However, if we only store which path is the fastest, we only need to store the points in time when the faster path switches. We expect significantly fewer switches than breakpoints. In this paper, we employ this alternative approach to adapt an existing speed-up technique to the time-dependent setting, describe engineering techniques employed in our implementation and present experimental results demonstrating that our approach significantly reduces memory consumption while achieving competitive query times.

Related Work Route planning in road networks has been extensively studied in the past decade. An overview over the field can be found in [2]. Here, we focus on speed-up techniques for time-dependent road networks. Several time-independent speed-up techniques have been generalized to the time-dependent setting. ALT [10], an approach using landmarks to obtain good A* [11] potentials, has been generalized to TD-ALT [12] and successively extended with node contraction to TD-CALT [4]. Even when combined with approximation, TD-CALT queries may take longer than 10 ms on continental sized graphs. SHARC [3], a combination of ARC-Flags [13] with shortcuts which allows unidirectional queries, was also extended to the time-dependent scenario [14]. It can additionally be combined with ALT yielding L-SHARC [14]. SHARC can find short paths in less than a millisecond but does not always find a shortest path. MLD/CRP [7,15] has been extended to TD-CRP [9] which can be used in a time-dependent setting. TD-CRP requires approximation to achieve reasonable memory consumption. It may find suboptimal paths. Another approach is FLAT [16] and its extension CFLAT [17]. CFLAT features sublinear query running time after subquadratic preprocessing and guarantees on the approximation error. Similar to our approach, CFLAT uses timestamped combinatoric structures to represent the changes in shortest paths over time and computes travel times lazily. Unfortunately, preprocessing takes long in practice and generates a prohibitively large index size.

There are several approaches based on Contraction Hierarchies [6]. Three were introduced by Batz et al. [8]: Time-dependent CH (TCH), inexact TCH and Approximated TCH (ATCH). TCH achieve great query performance but at the cost of a huge index size on state-of-the-art continental sized instances. The index size can be reduced at the cost of exactness (inexact TCH) or query performance (ATCH). An open-source reimplementation of the method in [8] named KaTCH (https://github.com/GVeitBatz/KaTCH (accessed on 15 March 2021)) exists. A simple heuristic named Time-Dependent Sampling (TD-S) was introduced by Strasser [18]. It samples a fixed set of scalar values from the time-dependent functions. It has manageable index sizes and fast query times but does not always find shortest paths.

Time-dependent shortest path algorithms are sometimes used as a subroutine in time-dependent vehicle routing problems (VRP) [19]. VRP research is a large field and discussing it beyond the scope of this article. For an overview over time-dependent VRP and other time-dependent routing problem variants, we refer to the work by Gendreau et al. [20].

Contribution In this work, we explore a variant of time-dependent Contraction Hierarchies, where shortcuts store paths instead of travel times. We introduce CATCHUp (Customizable Approximated Time-dependent Contraction Hierarchies through Unpacking), a time-dependent generalization of Customizable Contraction Hierarchies [5] and a thoroughly engineered implementation. Preprocessing takes only minutes even on modern production-grade continental sized instances with tens of millions of nodes. We also present algorithms which allow us to employ approximation to accelerate preprocessing without sacrificing exactness for the queries. Our implementation achieves fast and exact queries with performance competitive to TCH queries while requiring up to 38 times less memory.

This paper is an extended version of a conference paper [21]. In addition to the previously reported results, we describe our algorithms in greater depth and provide additional important engineering details. We also introduce new algorithms for profile queries. The experimental evaluation has been significantly extended. We perform experiments with more graphs and provide a deeper analysis on the performance of our preprocessing algorithms.

2. Materials and Methods

In this section, we describe our algorithms, data structures and implementation. After introducing preliminaries, we describe existing algorithms we build upon, most importantly Customizable Contraction Hierarchies. Section 2.2 introduces our shortcut unpacking data structure which we use to efficiently reconstruct the paths represented by a shortcut. We continue by presenting our our algorithms. The preprocessing, which computes auxiliary data from a road network with traffic predictions, is discussed in Section 2.3. In Section 2.4, we present query algorithms which utilize the auxiliary data to efficiently compute shortest travel times and paths between two given location.

2.1. Preliminaries

We model road networks as directed graphs $G=(V,A)$. A node $v\in V$ represents an intersection and an arc $a=uv\in A$ with $u,v\in V$ represents a road segment. Every arc a has a travel time function$f_a:\mathbb{R}\to {\mathbb{R}}^{>0}$ mapping departure time to travel time. These functions are also referred to as travel time profiles. We assume that travel time functions fulfill the First-In-First-Out (FIFO) property, that is, for any $\sigma ,\tau \in \mathbb{R}$ with $\sigma \le \tau$, $\sigma +f\left(\sigma \right)\le \tau +f\left(\tau \right)$ has to hold. Informally, this means that it is not possible to arrive earlier by starting later. If there are arcs that do not fulfill the FIFO property, the shortest path problem becomes $\mathcal{NP}$-hard [22] if waiting is not allowed. In our implementation, travel time functions are periodic piecewise linear functions represented by a sequence of breakpoints. We denote the number of breakpoints in a function f by its complexity$\left|f\right|$. A path is a sequence of nodes $[v_1,\dots ,v_k]$ such that $v_i{v}_{i+1}\in A$. We denote the concatenation of two paths by $[v_1,\dots ,v_k]·[v_k,\dots ,v_l]=[v_1,\dots ,v_k,\dots v_l]$. The travel time to traverse a path $[v_1,\dots ,v_k]$ can be evaluated by successively evaluating each link’s travel time: EVAL$([v_1,\dots ,v_k],\tau )=f_{v_1{v}_2}\left(\tau \right)+\mathrm{EVAL}([v_2,\dots ,v_k],\tau +f_{v_1{v}_2}\left(\tau \right))$.

Given two travel time functions f and g for arcs $uv$ and $vw$, we are often interested in the travel time function of traversing first $uv$ and then $vw$, that is $f\left(\tau \right)+g\left(f\right(\tau )+\tau )$. Computing this function is called linking. In a slight abuse of notation, we write $g\circ f$ for this linked function. When combining two travel time functions f and g for different paths $[u,\dots ,v]$ with the same start and end, we often want to know the travel time of the best path between u and v, that is $min(f,g)$. Computing this function is called merging. Both linking and merging can be implemented with coordinated linear sweeps over the breakpoints of both functions.

Given a departure time $\tau$ and nodes s and t, an earliest-arrival query asks for earliest point in time one can arrive at t when starting from s at $\tau$. Such a query can be handled by Dijkstra’s algorithm [1] with minor modifications [23]. The algorithm keeps track of the earliest known arrival time ${ea}_v$ at each node v. These labels are initialized with $\tau$ for s and ∞ for all other nodes. A priority queue is initialized with $(s,\tau )$. In each step, the node u with minimal earliest arrival ${ea}_u$ is popped from the queue and outgoing arcs are relaxed. To relax an arc $uv$, the algorithm checks if ${ea}_u+f_{uv}\left({ea}_u\right)$ improves ${ea}_v$ and updates label and queue position of v accordingly.

When nodes are popped from the queue, their earliest arrival is final. This property is denoted as label-setting. Once t is extracted from the queue, the earliest arrival at t is known. To retrieve the shortest path, one can use parent pointers which store the previous node on the shortest path from s for each node. We refer to this algorithm as TD-Dijkstra.

A profile query asks for the shortest travel time function between two nodes s and t for a time interval T. This query type can also be solved by a variant of Dijkstra’s algorithm. For each node v, instead of the earliest arrival time a tentative travel time function $f_v$ from s to v is maintained. Initially, $f_v\left(\tau \right)$ is set to ∞ for all $\tau \in T$ and to zero for s. In the priority queue, nodes are ordered by the current lower bound of their travel time function. The queue is initialized with $(s,0)$. When a node is popped from the queue, outgoing arcs are relaxed. Here, relaxing an arc $uv$ means linking $f_u$ with $f_{uv}$ and merging the result with the travel time function of v: $f_v=min(f_v,f_{uv}\circ f_u)$. If the travel time to v can be improved, its label and queue position will be updated accordingly.

This algorithm is not label-setting. Nodes may be popped several times from the queue. The algorithm can terminate as soon as the minimum key in the queue is greater than ${max}_{\tau \in T}{f}_t\left(\tau \right)$. We refer to this algorithm as TD-Profile-Dijkstra.

The A* algorithm [11] is an extension to Dijkstra’s algorithm. It reduces the number of explored nodes by guiding the search towards t. Each node u has a potential ${\rho }_t\left(u\right)$ which is an estimate of the distance to t. The priority queue is then ordered by ${ea}_u+{\rho }_t\left(u\right)$.

2.1.1. Contraction Hierarchies

Contraction Hierarchies (CH) [6] is a speed-up technique exploiting the inherent hierarchy in road networks. It was initially developed for networks with scalar edge weights. Nodes are heuristically ranked by their importance. Nodes with higher rank should cover more shortest paths. During preprocessing, all nodes are contracted in order of ascending importance. Contracting a node v means removing it from the network but preserving all shortest distances among remaining higher ranked nodes. For this, shortcut arcs are inserted between the neighbors of v if a shortest path goes through v. A shortcut is only necessary if it represents the only remaining shortest path between its endpoints. This can be checked with a local Dijkstra search (called witness search) between the endpoints. The result of the preprocessing is called an augmented graph.

Queries can be answered by performing a bidirectional Dijkstra search on the augmented graph. The forward search starts at s and relaxes only forward arcs to higher ranked nodes. The backward search starts at t and traverses arcs in reverse direction and also only searches to higher ranked nodes. The construction of the augmented graph guarantees that the searches will meet and find a path that has the same length as shortest paths in the original graph. The higher ranked nodes reachable from a node are referred to as the node’s CH search space.

2.1.2. Customizable Contraction Hierarchies

Customizable Contraction Hierarchies (CCH) [5] is a CH extension. It splits CH preprocessing into two phases where only the second uses weights. In the first phase, a separator decomposition and an associated nested dissection order [24,25] are computed. This order determines the node ranks. Nodes in the top-level separators have the highest ranks, followed by the nodes of each cell, recursively ordered by the same method. Since all shortest paths between different cells have to use separator nodes, these nodes cover many shortest paths. Thus, a nested dissection order is a good CH order.

Then, nodes are contracted iteratively ordered by ascending rank. Because weights are not considered in this phase, no witness search can be performed. All potential shortcuts between the higher ranked neighbors of the current node will be inserted. The upward neighbors become a clique. This phase is performed on the bidirected input graph with arcs $A^{\prime }=A\cup \{vu:uv\in A\}$ where each arc exists in both directions.

In the second phase called customization, arc weights are computed. Arcs in the augmented graph corresponding to an input graph arc are initialized with the corresponding weight. All other arc weights are set to ∞. Then, all arcs are processed in ascending order of their lower ranked endpoint. To process an arc $uv$, all lower triangles$[u,w,v]$, where w has lower rank than u and v are enumerated, checking if the path $[u,w,v]$ is shorter than $uv$. If so, the weight of $uv$ is set to the length of the $[u,w,v]$. We denote this as lower triangle relaxation.

The result of this basic customization fulfills all necessary properties for the CH query algorithm to find correct shortest distances. However, it contains many unnecessary arcs. These can optionally be removed using the perfect customization algorithm. Here, all arcs are processed in descending order of their higher ranked endpoint. For each arc $uv$, upper and intermediate triangles are enumerated, i.e., the third node w has greater rank than either u or v. The weight of $uv$ will then be decreased if possible to the weight of the path $[u,w,v]$. Once all arcs have been processed, all arcs where the weight changed can be removed. The augmented graph is now as small as possible.

The CH query algorithm can be reused without modifications. Another query algorithm is described in [5] which does not need priority queues. It is based on the elimination tree. A node’s parent in the elimination tree is its lowest-ranked upward neighbor in the augmented graph. The ancestors of a node in the elimination tree are exactly the set of nodes that are reachable in a CH search from this node [24]. Thus, instead of exploring the search space through Dijkstra searches, the elimination tree query searches the same nodes by traversing the path to the root in the elimination tree.

2.2. Shortcut Unpacking Data

The key element of our approach is the information we store with each arc of the augmented graph. We store time-dependent unpacking information, which allows us to efficiently reconstruct the original path represented by a shortcut for a point in time. In (C)CH, shortcuts $uv$ are inserted when a node w is contracted and the arcs $uw$ and $wv$ exist. Thus, a shortcut $uv$ always skips over a triangle $[u,w,v]$. However, there may be several triangles and which one is the fastest may change over time. There may also be an arc $uv$ in the input graph which might sometimes be the fastest path. This is the information our unpacking data structure has to capture.

For each arc $uv$, we store a set of time-dependent expansions$X_{uv}$ for unpacking. Figure 3 presents an example. For an expansion $x\in X_{uv}$, we denote the time during which x represents the shortest path as the validity interval ${\Pi }_x$ of x. When formally referring to the path represented by an expansion, we use the expand function $E:X\to V\cup A$. E either maps to an original arc or to the middle node $w_x$ of the lower triangle $[u,w_x,v]$. Knowing the middle node for each expansion is also sufficient to obtain longer paths. These can be computed by unpacking shortcuts recursively.

In our implementation, the expansion information is represented as an array of triples $(\pi ,u{w}_x,w_{x}v)$. $\pi$ is the beginning of the validity interval and $uw$ and $wv$ are arc ids. This information can be stored in 16 bytes for each entry—8 bytes for the timestamp and 4 bytes for each arc id. An expansion can also represent an original arc or no arc at all during a certain time interval. Both these cases are represented as special arc id values. Two invalid arc ids indicate no arc at all. One invalid id indicates that the other arc id represents an original arc. Beside the expansions, we also maintain a scalar lower bound ${\underset{¯}{b}}_{uv}$ and an upper bound ${\overline{b}}_{uv}$ on the travel time for each arc.

During preprocessing, we have to compute the bounds and expansion sets for each arc in the augmented graph. This is done using the same schema as in CCH. We iterate over all arcs and relax their lower triangles. Algorithm 1 depicts the routine for each triangle. The routine requires travel time functions for all involved arcs. We maintain these functions during preprocessing but discard them later. To relax the lower triangle $[u,w,v]$, the functions $f_{uw}$ and $f_{wv}$ are linked and the result is merged with the current function of $f_{uv}$. Where $[u,v,w]$ is faster, new expansions are inserted into $X_{uv}$. Where the current $uv$ travel time is faster, the current expansions are kept. The bounds are updated with the new minimum and maximum of the merged function.

Once the unpacking information for an arc is complete, we can use it to compute the arc’s travel time, unpack it to the path in the original graph or compute the travel time function for the arc. All these operations follow the same basic schema: Determine the relevant expansions and apply the operation recursively until arcs from the input graph are reached. The simplest case is the travel time evaluation. Algorithm 2 depicts this operation. First, the relevant expansion is determined using binary search. If it points to an original arc, this arc’s travel time can be evaluated and returned. If the expansion points to a lower triangle, we first recursively evaluate the first arc of the triangle. Then, the second arc can be evaluated at the entry time plus the travel time of the first arc.

Algorithm 1:LowerTriangleRelax

Input: Preliminary unpacking data for arc $uv$: $X_{uv},{\underset{¯}{b}}_{uv},{\overline{b}}_{uv}$ and travel time function $f_{uv}$. Travel time functions of lower triangle arcs $f_{uw}$ and $f_{wv}$.

Output: Improved data $X_{uv},{\underset{¯}{b}}_{uv},{\overline{b}}_{uv}$ and function $f_{uv}$.

1$g\leftarrow f_{wv}\circ f_{uw}$

2${\Pi }_{uv}\leftarrow \{\tau \mid f_{uw}\left(\tau \right)\le g\left(\tau \right)\}$

3${\Pi }_{\left[uwv\right]}\leftarrow \{\tau \mid g\left(\tau \right)<f_{uw}\left(\tau \right)\}$

4Xuv ← {NewExpansion(x, Πuv ∩ Πx)∣ x ∈   Xuv} ∪ {NewExpansion([u, w, v], Π[uwv])}

5$f_{uv}\leftarrow min(f_{uv},g)$

6${\underset{¯}{b}}_{uv}\leftarrow min({\underset{¯}{b}}_{uv},{min}_{\tau }\left(f_{uv}\left(\tau \right)\right))$

7${\overline{b}}_{uv}\leftarrow min({\overline{b}}_{uv},{max}_{\tau }\left(f_{uv}\left(\tau \right)\right))$

8return$X_{uv},{\underset{¯}{b}}_{uv},{\overline{b}}_{uv},f_{uv}$

Algorithm 2:Eval

Algorithm 3 depicts the procedure for determining the path represented by an expansion set for a given time. The recursive schema is the same as for Eval but the result is a path instead of a travel time. Nevertheless, when unpacking lower triangles, we still need to evaluate the first arcs travel time to determine the time for unpacking the second arc.

Algorithm 3:Unpack

Constructing the travel time function is also similar and shown in Algorithm 4. We recursively unpack expansions until we reach arcs of the original graph where exact travel time functions are available. We may need to unpack several expansions for different times and combine them. For each expansion we check if its validity overlaps with the time range for which we want to construct the travel time function. If so, we recursively retrieve the function for the first arc during this overlap. Then, we calculate the function for the second arc during the overlap. For the second arc, the time interval must be shifted by the travel time of the first arc at the start and end of the time interval. Both functions are then linked and appended to the final function.

Implementing this algorithm naively may cause performance issues since many memory allocations are performed for intermediate results. We avoid this by keeping all intermediate results in two buffers which are reused for all invocations of this algorithm. The buffers are stored as dynamically sized arrays (C++ vectors) and can grow dynamically but will never shrink. Once they have grown to an appropriate size, no more memory allocations will be necessary. Each buffer can contain many travel time functions stored consecutively. The link operation will read the last two functions from one buffer and append the result to the other buffer. Then, the two input functions will be truncated from the first buffer. After swapping, the buffers can be used again for the next link operation. Swapping is necessary, because it is not possible to read from and write to the same buffer during the same operation. The same schema can be employed for joining partial functions (see Figure 4 for a visualization).

Algorithm 4:ReconstructTravelTimeFunction

2.3. Preprocessing

In this section, we present our preprocessing algorithms. The first phase of CCH preprocessing is performed only on the topology of the graph. Since no travel time functions are involved, we can adapt the algorithms of Dibbelt et al. [5] without modification. We use InertialFlowCutter [26] to obtain the nested dissection order. To generate the augmented graph, we implement an improved contraction algorithm first presented in [27]. When contracting a node, we insert all upward neighbors of the current node only into the neighborhood of its lowest ranked upward neighbor. This algorithm can be implemented to run in linear time in the size of the output graph.

The goal of the second phase of preprocessing—the customization—for classical CCH is to compute the shortcut weights. For our approach, we have to compute the travel time bounds and time-dependent expansions for all arcs in the augmented graph. Recall that a shortcut $uv$ always bypasses one or many lower triangles $[u,w_i,v]$ for different nodes $w_i$, where $w_i$ has lower rank than u and v. For the bounds, we want to find the minimum and maximum travel time of the fastest travel time function between u and v over any $w_i$. For the expansions, we need to determine for each point in time which lower triangle is the fastest. Assuming we know the final travel time functions of all $u{w}_i$ and $w_{i}v$, we can compute this using the LowerTriangleRelax algorithm (see Algorithm 1). This leads to the following algorithmic schema: Maintain a set of necessary travel time functions in memory starting with the functions from the input graph. Iterate over all arcs in the augmented graph in a bottom-up fashion. For each arc enumerate lower triangles. Link and merge their functions to compute the function, bounds and expansions of the current arc. Keep the current arc’s travel time function in memory until it is no longer needed.

We implement this schema as follows: We process all arcs $uv$ ordered ascending by their lower ranked endpoint. Since the middle node w of a lower triangle $[u,w,v]$ has always lower rank than u and v, the arcs $uw$ and $wv$ will have been processed already. To process an arc $uv$ we enumerate lower triangles $[u,w,v]$. Perform LowerTriangleRelax for each triangle in both directions. Once all arcs $uv$ have been processed where u is the lower ranked endpoint, we drop the travel time functions of all arcs $wu$ where u is the higher ranked endpoint. Algorithm 5 depicts this in pseudo code.

Algorithm 5:Customization

When enumerating triangles, we order them ascending by ${\underset{¯}{b}}_{uw}+{\underset{¯}{b}}_{wv}$. This way, we process triangles first, which are likely faster. This gives us preliminary bounds on the travel time of $uv$. Before linking the functions of another triangle $f_{uw}$ and $f_{wv}$, we check if ${\overline{b}}_{uv}\le {\underset{¯}{b}}_{uw}+{\underset{¯}{b}}_{wv}$. If so, the linked path would be dominated by the shortcut, and we can skip linking and merging completely. If not, we link $f_{uw}$ and $f_{wv}$ and obtain $f_{[u,w,v]}$. We still can skip merging if one function is strictly smaller than the other, that is either ${\overline{b}}_{uv}\le min\left(f_{[u,w,v]}\right)$ or $max\left(f_{[u,w,v]}\right)\le {\underset{¯}{b}}_{uv}$. Even if the bounds overlap, one function might still dominate the other. To check for this case, we simultaneously sweep over the breakpoints of both functions, determining the value of the respectively other function by linear interpolation. Only when this check fails, we perform the merge operation.

Before the time-dependent customization, we first use the basic and perfect customization algorithms from [5] to compute preliminary scalar upper and lower bounds for all arcs. With these bounds, we can skip additional linking and merging operations. Employing perfect customization, we can remove some arcs completely, when a dominating path through higher ranked nodes exists.

Parallelization We employ both loop based and task based parallelism. The original CCH publication [5] suggests processing arcs with their lower ranked endpoint on the same level in parallel. The level of a node is the maximum level of its downward neighbors increased by one, or zero if the node does not have downward neighbors. We use this approach to process arcs in the top-level separators. However, this approach requires a synchronization step on each level. This is detrimental to load balancing. Thus, we use a different strategy when possible.

In [28], a task based parallelism approach utilizing the separator decomposition of the graph is proposed. Each task is responsible for a subgraph $G^{\prime }$. Removing the top-level separator in $G^{\prime }$ decomposes the subgraph into two or more disconnected components. For each component, a new task is spawned to process the arcs the component. After all child tasks are completed, the arcs in the separator are processed utilizing the loop based parallelization schema. If the size of subgraph $G^{\prime }$ is below a certain threshold, the task processes the arcs in $G^{\prime }$ sequentially without spawning subtasks. We use $n/(\alpha ·c)$ as the threshold with c being the number of cores and the tuning parameter $\alpha$ set to 32, as suggested by Buchhold et al. [28].

Approximation As we process increasingly higher ranked arcs, the associated travel time functions become increasingly complex. This leads to two problems. First, linking and merging becomes very time-consuming as running times scale with the complexity of the functions. Second, storing these functions—even though it is only temporary—requires a lot of memory. We employ approximation to mitigate these issues. However, for exact queries, we need exact unpacking information. We achieve this by lazily reconstructing parts of exact travel time functions during merging.

When approximating, we do not store one approximated function but two—a lower bound function and an upper bound function with maximum difference $ϵ$ where $ϵ$ is a configurable parameter. These approximations replace the exact function stored for later merge operations and will also be dropped when no longer needed. To obtain the bound functions, we first compute an approximation using the algorithm of Douglas and Peucker [29]. Previous works [8,9] have reported using the algorithm of Imai and Iri [30] for approximation. Given a maximum error bound $ϵ$, this algorithm can compute in linear time the piecewise linear function with the minimum number of breakpoints within the given bound. The Douglas–Peucker algorithm has a quadratic worst case running time and no such guarantees on the number of breakpoints in the approximation. However, the theoretic guarantees of the Imai–Iri algorithm come at the cost of considerable implementation complexity and high constant runtime factors. Preliminary experiments showed that, compared to Imai–Iri, our Douglas–Peucker implementation actually produces insignificantly more breakpoints and also runs faster due to better constants. In addition, the implementation needs 30 instead of 400 lines of code, so we use the Douglas–Peucker variant. Then, we add or subtract $ϵ$ to the value of each breakpoint to obtain an upper or lower bound, respectively. This bounds are valid, but they may not be as tight as possible. Therefore, we iterate over all approximated points and move each point back towards the original function. Both adjacent segments in the approximated functions have a minimum absolute error to the original function. We move the breakpoint by the smaller of the two errors. This yields sufficiently good bounds.

When linking approximated functions, we link both lower and both upper bound functions. Linking two lower bounds yields a valid lower bound of the linked exact functions because of the FIFO property. The same argument holds for upper bounds.

Merging approximated functions is slightly more involved. Our goal is to determine the exact expansions for each arc. We use the approximated bounds to narrow down the time ranges where intersections are possible. To identify these parts, we merge the first function’s lower bound with the second function’s upper bound and vice versa. Where the bounds overlap, an intersection might occur. We then obtain the exact functions in the overlapping parts using Algorithm 4 and perform the exact merging. To obtain approximated upper and lower bounds of the merged function, we merge both lower bounds and both upper bounds (see Figure 5 for a visualization).

We approximate whenever a function has more than $\beta$ breakpoints after merging. This includes already approximated functions. Both $\beta$ and the maximum difference $ϵ$ are tuning parameters which influence the performance (but not the correctness). We evaluate different choices in Section 3.3.1.

2.4. Queries

2.4.1. Earliest Arrival Queries

Recall that, for an earliest arrival query, we are given a source node s, a target t and a departure time $\tau$ from s. The goal is to obtain the earliest arrival at t. Compared to a standard (C)CH query, our query algorithm has to deal with two challenges. First, we cannot simply perform a backwards search, as we do not know the arrival time at the target node. Second, to evaluate the travel time of a shortcut, we need to obtain the path in the original graph which is an expensive operation. To address the first challenge, the query is split in two phases. In the first phase we obtain a subgraph on which we can run a forward-only Dijkstra-like search in the second phase. We now present the basic query algorithm and later introduce optimizations to address the second challenge.

In the first phase, the union of the subgraphs reachable through arcs to higher ranked neighbors from s and t is obtained. We construct these subgraphs by traversing the elimination tree starting from both s and t to the root and marking all encountered arcs as part of the search space. The backward search from t maintains parent pointers to represent the subgraph: For each encountered arc $uv$ (where v has the higher rank), we store the arc id and the tail u at v. This allows efficiently traversing these downward arcs in the forward-only Dijkstra in the second phase. In the second phase, we run Dijkstra’s algorithm on the combined search spaces. Shortcut travel times are evaluated with Algorithm 2.

By the construction of CH, the search space contains the shortest path, thus Dijkstra’s algorithm will find it and this algorithm will determine the correct earliest arrival at t. However, the search space is bigger than strictly necessary. This slows down the query. In the next paragraph, we discuss how to construct smaller subgraphs using an elimination tree interval query.

Elimination Tree Interval Query The elimination tree interval query is a bidirectional search starting from both the source node s and the target node t. It constructs a smaller subgraph for the second phase. We denote this subgraph as a shortest path corridor. Node labels contain an upper ${\overline{t}}_v$ and a lower bound ${\underset{¯}{t}}_v$ on the travel time to/from the search origin and a parent pointer to the previous node and the respective arc id. The bounds ${\overline{t}}_s$, ${\overline{t}}_t$, ${\underset{¯}{t}}_s$, ${\underset{¯}{t}}_t$ are all initialized to zero in their respective direction, all other bounds to infinity. We also track tentative travel time bounds for the total travel time from s to t. For both directions, the path from the start node to the root of the elimination tree is traversed. For each node u, all arcs $uv$ to higher ranked neighbors are relaxed, that is checking if ${\overline{t}}_u+{\overline{b}}_{uv}<{\overline{t}}_v$ or ${\underset{¯}{t}}_u+{\underset{¯}{b}}_{uv}<{\underset{¯}{t}}_v$ and improving the bounds of v if possible. When the new travel time bounds from an arc relaxation overlap with the current bounds, more than one label has to be stored. As an optimization [28], nodes can be skipped if the lower bound on the travel time to it is already greater than the tentative upper bound on the total travel time between s and t. After both directions are finished, we have several nodes in the intersection of the search spaces. Where the sum of the forward and backward distance bounds of such a node overlaps with the total travel time bounds, the parent pointers are traversed to the search origin and all arcs on the paths are marked as part of the shortest path corridor.

Lazy Shortcut Unpacking In the second query phase, we perform Dijkstra’s algorithm on the corridor obtained in the first phase. In the basic query, shortcuts are unpacked completely to evaluate their travel time. However, this may cause unnecessary and duplicate unpacking work. We now describe an optimized version of the algorithm which performs unpacking lazily. The algorithm starts with the same initialization as a regular TD-Dijkstra. All earliest arrivals are set to infinity, except for the start node which is set to the departure time. The start node is inserted into the queue. Then, nodes are popped from the queue until it is empty or the target node is reached. For each node, all outgoing arcs within the shortest path corridor are relaxed. When an arc is from the input graph, its travel time function can be evaluated directly. Shortcut arcs, however, need to be unpacked. The lazy unpacking algorithm defers as much work as possible: Only the first arc of the triangle of each shortcut will be recursively unpacked until an input arc is reached, the second arc will be added to the corridor. Figure 6 shows an example. This way, we unpack only the necessary parts and avoid relaxing arcs multiple times when shortcuts share the same paths.

Corridor A* The query can be accelerated further, by using the lower bounds obtained during the elimination tree interval query as potentials for an A*-search. For nodes in the CH search space of t, the lower bounds from the backward search can be used. For nodes in the CH search space of s, we start at the meeting nodes from the corridor search and propagate the bounds backwards down along the parent pointers. This yields potentials for all nodes in the initial corridor. However, we also need potentials for nodes added to the corridor through unpacking. These potentials are computed during the shortcut unpacking. When unpacking a shortcut $uv$ into the arcs $uw$ and $wv$, then $\rho \left(w\right)$ will be set to $min(\rho \left(w\right),\rho \left(v\right)+{\underset{¯}{b}}_{wv})$.

Justifying that A* with these potentials will always find the correct earliest arrival is surprisingly non-trivial. In fact, these potentials are not feasible in the sense that $\forall \tau \in T,uv\in A:f_{uv}\left(\tau \right)-\rho \left(u\right)+\rho \left(v\right)\ge 0$ [10]. Figure 7 shows an example where the term becomes negative and the same node has to be popped several times from the queue. Assume that all arcs have a constant travel time for the departure time of this query and lower bounds are equal to the travel time. The exception is $v_{3}t$ which has constant travel time 100 during this query but the lower bound is zero. We use zero weights to simplify the example. They are not strictly necessary for such an example. The shortest path from s to t is $[s,w_2,w_1,v_2,t]$ and has length 22. After s is settled, the queue will contain $v_3$ with key $1+0$ (distance plus lower bound to t), $v_1$ with key $1+2$ and $w_2$ with key $1+21$. Then, $v_3$ will be settled which will insert t with key $101+0$ into the queue. Then, $v_1$ will be settled and $w_1$ will be inserted into the queue with key $3+0$. Then, $w_1$ will be settled even though the current distance of 3 is greater than the actual shortest distance of 2. This will insert $v_2$ with key $13+10$ into the queue. Now, $w_2$ will be popped and the distance to $w_1$ will be improved and it will be reinserted into the queue with key $2+20$. $w_1$ will be popped immediately afterwards which improves the distance and key of $v_2$ which is the next node to be popped from the queue. After it has been processed, the final distance to t (22) is known, and t is the final node to be settled.

Nevertheless, we claim that once t is popped from the queue the algorithm always has found the correct earliest arrival. The reason is that for all nodes v on the shortest path ${ea}_v+\rho \left(v\right)\le {ea}_t$ holds. Since this is the queue key, all nodes on the shortest path will have been popped (possibly several times) before t. Let $P=[s,\dots ,t]$ be the desired shortest path from s to t when departing from s at $\tau$ and $P[v..]\subset P$ the subpath from a node $v\in P$ to t. We denote the global lower bound on the travel time between a node v and t by ${\underset{¯}{f}}_{vt}$ and the lower bound travel time on $P[v..]$ by ${\underset{¯}{f}}_P[v..]$. For nodes in the initial corridor obtained by the interval query, ${ea}_v+\rho \left(v\right)\le {ea}_t$ always holds because the potential $\rho \left(v\right)$ is set to the global lower bound ${\underset{¯}{f}}_{vt}$. However, nodes u added later to the corridor through unpacking may have a greater potential than ${\underset{¯}{f}}_{ut}$, as depicted in the example. However, their final potential cannot be greater than the lower bound of the travel time on the shortest path ${\underset{¯}{f}}_P[u..]$. This is enough to satisfy ${ea}_u+\rho \left(u\right)\le {ea}_t$. In addition, the potential value will be set to this final value before t is settled. Assume for contraction that $p_i$ is the first node on P for which this statement does not hold. $p_i$ cannot be a node from the initial corridor. However, $\rho \left(p_i\right)$ will be set at most to ${\underset{¯}{f}}_P[p_i..]$ once $p_{i-1}$ is settled which by assumption happens before t is settled. This is a contraction. Thus, ${ea}_v+\rho \left(v\right)\le {ea}_t$ holds for all nodes $v\in P$ and the query algorithm always finds the correct earliest arrival when terminating once t is popped from the queue.

2.4.2. Profile Queries

A profile query asks for the function of the fastest travel time between two nodes over a given time period T. Without loss of generality, we assume that the T equals the entire time period covered by the input network. As discussed in Section 2.1, such a query can be answered with TD-Profile-Dijkstra. However, TD-Profile-Dijkstra exhibits both prohibitive running time and memory consumption. Consider a path $[v_0,\dots ,v_k]$ where the travel time function of each arc has b breakpoints. In general, linking two functions f and g creates a new function with $\mathsf{\Theta }\left(\right|f|+|g\left|\right)$ breakpoints. Thus, the travel time function from $v_0$ to node $v_i$ contains $\mathsf{\Theta }(i·b)$ breakpoints. Thus, when applying TD-Profile-Dijkstra to compute the function between $v_0$ and $v_k$, the total memory consumption and the running time grows quadratically with the length of the path. We conclude that Dijkstra-based approaches to profile queries are not a promising direction. The experiments with Dijkstra-based TCH profile queries in [8] support this conclusion. We also performed preliminary experiments with a proof-of-concept implementation where we adapted our earliest arrival query algorithm to the profile query setting. The Dijkstra-based approach was more than an order of magnitude slower than the approach described in the following. Instead of a Dijkstra-like search, our algorithm uses contraction and methods from the preprocessing to construct the desired profile.

Our algorithm has four phases. The first phase uses the elimination tree interval query to obtain a shortest path corridor and is the same as for earliest arrival queries. In the second phase, we obtain travel time functions for all arcs in the shortest path corridor. During the third phase, additional shortcuts will be inserted and their unpacking data will be computed, reusing the preprocessing algorithms. The result is a new $st$ shortcut with exact expansions $X_{st}$. From this unpacking information, the exact travel time function (among other things) can be obtained in the fourth phase. We now describe Phases 2–4 in detail.

Reconstruction In this phase, we obtain travel time functions for all arcs in the shortest path corridor. Similar to the customization, the obtained travel time functions may be either exact, or approximated upper and lower bound functions. This keeps the memory consumption low and linking and merging operations fast. We obtain these functions by first recursively reconstructing the travel time functions of all arcs referenced by expansions. The reconstructed functions will be saved for both arcs in the corridor and unpacked arcs, in case another reconstruction operation might reuse an arc. If all arcs referenced by the expansions have an exact function available, we can compute an exact travel time function for the current arc. If not, we end up with an approximation. After reconstruction, we always check if a function has more than $\beta$ breakpoints. If that is the case, we approximate it to reduce the complexity as described in Section 2.3.

Contraction In the third phase, we insert additional shortcuts and compute their unpacking data by simulating the contraction of the nodes in the shortest path corridor. We reuse the existing nested dissection order. The ranks of s and t are increased such that they are higher in the hierarchy than all other nodes. This construction leads to a shortcut between s and t, one between s and each node in the corridor on the path from s to the elimination tree root and one from each node in the corridor on the path from t to the root to t. Some of these shortcuts may already exist (see Figure 8 for an illustration).

These new shortcuts are now processed as in Algorithm 5 to compute their unpacking data. We initialize the shortcut bounds with the bounds obtained from the elimination tree query. This allows to prune unnecessary operations. We process shortcuts ordered by their lower ranked endpoint. For each shortcut we enumerate and relax lower triangles using Algorithm 1. We can enumerate these triangles efficiently using the parent pointer from the interval query. Each shortcut has an endpoint node in the corridor. The parent pointers of this node correspond to the triangles that need to be relaxed. The shortcuts from s and the shortcuts to t are independent of each other and can be processed in parallel. The lower triangles of the $st$ shortcut can be enumerated by iterating over the meeting nodes from the interval query. We also employ the triangle sorting optimization.

Extraction In this final phase, we can use the unpacking information of the $st$ shortcut to efficiently compute the final result. The shortcut already contains a possibly approximated travel time function from the contraction phase. This may suffice for some applications. If the shortcut contains only an approximation, but we need an exact travel time profile, we can use Algorithm 4 to compute it. For some practical applications, the different shortest paths over the day may be more useful than the travel time profile. Algorithm 6 depicts a routine to compute path switches and the associated shortest paths. The algorithm follows the same schema as all unpacking algorithms. The operation is recursively applied to all expansions limited to the validity time of the expansion. Only the Combine operation is more involved. It performs a coordinated linear sweep over the path sets from $u{w}_x$ and $w_{x}v$ and appends the paths where the validity intervals overlap. For the paths from $w_x$ to v, we only know the validity times with respect to departure at $w_x$. To obtain the corresponding departure time at u, we reverse evaluate the current $u{w}_x$ path, i.e., we successively evaluate the inverted arrival time function of all arcs on the path in reverse order.

Algorithm 6:UnpackPaths

3. Results

In this section, we present our experimental results. We first discuss the experimental setup and the input road networks. Then, we discuss the performance of each of our presented algorithms in turn. Finally, we compare our approach to related work.

3.1. Experimental Setup

Our benchmark machine runs openSUSE Leap 15.2 (kernel 5.3.18) and has 192 GB of DDR4-2666 RAM and two Intel Xeon Gold 6144 CPUs, each of which has eight cores clocked at 3.5 GHz and 8 × 64 KB of L1, 8 × 1 MB of L2 and 24.75 MB of shared L3 cache. Hyperthreading was disabled and parallel experiments use 16 threads. We implemented our algorithms in Rust (the code for this paper and all experiments is available at https://github.com/kit-algo/catchup (accessed on 15 March 2021)) and compiled them with rustc 1.49.0-nightly (cf9cf7c92 2020-11-10) in the release profile with the target-cpu = native option (we disabled AVX512 instructions, as they caused misoptimizations). To compile competing implementations written in C++, we use GCC 9.3.1 using optimization level 3.

We investigated the performance of our preprocessing and query algorithms and compare dit to competing algorithms. Our experiments were focused on but not limited to space consumption and running times. We performed preprocessing five times for each input network and report arithmetic means of the running times. Unless reported otherwise, preprocessing utilized all 16 cores. For queries, we generated 100,000 source, target, departure time triples chosen uniformly at random for each graph. These were executed in bulk. Competing algorithms were evaluated with the same query set. For profile queries, we only used 1000 queries (and discarded the departure time). We report arithmetic means of query running times and machine independent measures such as number of nodes popped from the queue and number of evaluated travel time functions.

3.2. Input Road Networks

In this section, we report results for five road networks: an instance of Berlin and the surrounding area, an old instance of Germany from 2006 which is used in many related works, two production-grade instances for Germany and Europe and with traffic predictions from 2017 and a recent production-grade instance of Europe. All instances include traffic predictions as piecewise linear functions. We use traffic predictions for a car on a typical midweek day or a Tuesday.

Table 1. Characteristics of test instances used. The third column contains the percentage of arcs with a non-constant travel time function. The fourth column the average number of breakpoints among those. The fifth and sixth columns report the relative total delay for all/only non-constant arcs. The final column contains the graph size in a compact representation in memory.

	Nodes [·103]	Arcs [·103]	TD Arcs [%]	Avg. |f| per TD Arc	Rel. Delay [%]	Rel. Delay TD [%]	Size [GB]
Ber	443.2	988.5	27.4	75.0	3.1	17.6	0.2
Ger06	4688.2	10,795.8	7.2	19.5	1.7	33.1	0.3
Ger17	7247.6	15,752.1	29.2	31.6	3.5	20.8	1.3
Eur17	25,758.0	55,503.8	27.2	29.5	2.7	19.0	4.2
Eur20	28,510.0	60,898.8	76.3	22.5	21.0	34.9	8.7

lists characteristics of each graph.

For time-independent routing, the performance of algorithms primarily depends on the size of the network. In our case, however, the amount and complexity of time-dependent information also has a significant impact on the performance. To measure this, we report the fraction of arcs which have a non-constant travel time and the average number of breakpoints among all non-constant travel time functions. We also report the relative total delay $\frac{{\sum }_{a}max{f}_a-min{f}_a}{{\sum }_{a}min{f}_a}$ as a measure for the degree of time-dependency of the predictions. (Variants of this measure have been used in previous works. Delling reported the average relative delay of time-dependent earliest arrival queries over the result of a time-independent query with lower bound travel times in [14]. Batz reported the average relative delay $\frac{max{f}_a-min{f}_a}{min{f}_a}$ in [8]. We use the total delay because averages of ratios have hard to interpret semantics [31]. For example, a short arc with a large relative delay could have a much bigger influence on the average relative delay than it has on the actual shortest path structure.). The smaller is the relative delay, the greater is the effectiveness of pruning with upper and lower bounds.

The Ber instance was provided by TomTom (https://www.tomtom.com (accessed on 15 March 2021)) in 2013. It contains the northern eastern part of Germany, but in the literature it is referred to by the largest city, Berlin. This instance features the most complex input functions and 27.4% of the arcs have a non-constant travel time function. However, the road network is quite small which makes the instance comparatively easy. All other instances were provided by PTV (https://ptvgroup.com (accessed on 15 March 2021)).

The Ger06 instance is the easiest with respect to all time-dependency measures except the total relative delay limited to non-constant travel time functions. That means that the degree of time-dependency in the non-constant functions is significant. However, the overall influence on the shortest path structure is limited, because only 7.2% of the arcs have a non-constant travel time. Ger06 has also the smallest average complexity of non-constant travel time functions.

The newer graphs are not only bigger in terms of number of nodes and arcs but have also significantly more non-constant travel time functions. Ger17 has four times as many time-dependent arcs as Ger06 and about 1.5 times as many breakpoints per time-dependent function. Even though the relative delay among non-constant functions is not as high as for Ger06, the relative delay among all arcs is twice as high. Eur17 exhibits similar characteristics but additionally has 3.5 times as many nodes and arcs.

Our newest instance is Eur20 with 28.5 M nodes and 61 M arcs. Three quarters of the arcs have a non-constant travel time. With around 35%, the delay among non-constant functions is the greatest among all instances. The total delay among all arcs is more than an order of magnitude higher than on Ger06. This makes it the hardest instance.

We also performed experiments with predictions for different weekdays and the Western Europe graph provided for the 9th DIMACs implementation challenge [32] with synthetic travel time functions [12]. The results for these networks and prediction sets did not provide much additional insight. We report them in Appendix A for the sake of comparability and completeness.

3.3. Preprocessing

Table 2. Preprocessing statistics. Running times are for parallel execution with 16 cores.

	CCH Arcs	Expansions per Arc			Index	Running Time [s]
	[·103]	Avg.	Max.	=1 [%]	[GB]	Phase 1	Phase 2
Ber	1977	1.039	31	98.6	0.09	1.5	6.2
Ger06	22,519	1.075	44	98.4	1.06	30.1	21.6
Ger17	31,488	1.090	107	98.5	1.50	35.0	107.4
Eur17	114,857	1.099	115	98.4	5.47	189.6	557.0
Eur20	128,921	1.191	109	96.9	6.32	209.6	1039.5

reports the results for our preprocessing. On Ger06, the first preprocessing step takes longer than the second. However, for the newer instances with more time-dependent arcs and more breakpoints per function this changes and the second step becomes more expensive. Despite that, the size of the final index corresponds only to the number of arcs in the augmented graph and does not grow as much for the newer instances. The high complexity of the input function on Ber also does not depict any negative influence on the index size or the number of expansions per shortcut. The augmented graphs have about twice as many arcs as the original graphs. On average, only 1.1 expansions per arc need to be stored for all graphs (1.2 for Eur20). About 98% of all arcs have only a single expansion. The maximum number of expansions per arc is only 115. This is two orders of magnitude less than the number of breakpoints in the travel time function of that arc. On Eur20, our hardest instance, the total preprocessing time is about 20 min, roughly twice as much as for Eur17. However, the index size grows by less than 1 GB and is in fact smaller than the input graph. This clearly demonstrates the advantage in space efficiency of expansion information over explicitly storing travel time functions.

3.3.1. Customization

We now focus on the second preprocessing phase. The first preprocessing phase does not use any time-dependent information and only reuses existing algorithms which have been evaluated in greater detail in other works [26,27]. For the scope of this subsection, we only use the Eur20 graph.

To evaluate their impact, we selectively disable the triangle sorting and precustomization optimizations. Even though both optimizations speed-up the customization by improving bounds, both have a significant impact on their own. Disabling the precustomization increases the overall customization running time by about 5 min to 1311 s. The effect of triangle sorting is even stronger. Disabling it roughly doubles customization running time to 2156 s.

Parallelization We evaluate the effectiveness of our parallelization schema and run the customization with a varying number of threads. Figure 9 depicts the results. As a baseline, we run the experiment with all parallelization code disabled. The baseline running time is indicated by the dashed line. Enabling parallelization but running with only one thread causes only little overhead. Running with more threads introduces more overhead due to synchronization. With 16 threads, parallel efficiency is still around 0.9. We conclude that our parallelization schema scales well and that customization times could be reduced further by utilizing additional cores.

Approximation We perform customization experiments with different approximation parameters. Over the course of the customization, we track the progress over time, the memory consumption and the average travel time function complexity. Figure 10 displays these measurements. We use the number of processed triangles to measure the progress because it corresponds roughly linearly to the time passed (though different parameters lead to different slopes). After about 60% of the triangles, the slope changes slightly. At this point only high-level separator nodes/arcs remain. These have complex travel time functions so linking and merging becomes more expensive. In addition, we switch from task based to loop based parallelization which is less effective. Measuring progress by processed nodes (the for loop in Line 1 in Algorithm 5) or processed arcs (for loop on Line 2) is also insightful. However, for these, the correspondence is not linear. The last couple of thousand nodes and the last million arcs take almost half the total time.

We observe that the choice of approximation parameters has a huge influence on running time and memory consumption. The best running time of around 1000 s is achieved with $\beta =1000$ and $ϵ=1.0$ s. Thus, we use it as our default configuration. The worst running time is over six times higher. In the best configuration, we use only around 25 GB for travel time functions while a bad parameter choice or no approximation leads to crashes with out-of-memory errors. Generally, a larger $ϵ$ leads to looser approximation, lower travel time function complexity and thus less memory consumption. Conversely, a larger $\beta$ causes the approximation to be executed less often and the memory consumption and function complexity increases. The average function complexity is usually well below $\beta$ except for very small values of $ϵ$ or $\beta$. In that case the complexity cannot be reduced sufficiently to keep the complexity below $\beta$. If $\beta$ is large, approximation is performed seldom and the influence of $ϵ$ becomes smaller. Similarly, when $ϵ$ is small, the influence of $\beta$ is limited because the complexity cannot be reduced enough no matter how often approximation is performed. Clearly, the right choice of approximation parameters is essential to the performance of the preprocessing. When travel time functions are too complex, too much memory is used, and linking and merging are very expensive. However, when functions are approximated too loosely, a lot of time is spent in the reconstruction of exact functions for the times when bounds overlap. Thus, extreme parameter choice for $ϵ$ and $\beta$ are detrimental to the running time, even though they may reduce memory consumption.

Through all configurations, the memory usage peaks after around 60% of the triangles have been processed. The reason is the way we maintain travel time functions in memory during the customization. An arc’s travel time function is stored once its lower ranked endpoint has been processed until its higher ranked endpoint has been processed. In the beginning, we process nodes with low rank and store many travel time functions. Only once we reach the higher ranked nodes do we start dropping a significant amount of the stored functions. This causes the observed peak.

3.4. Queries

In this section, we investigate the performance of our query algorithms.

Table 3. Query performance with different optimizations. We report the number of nodes popped from the queue, the number of evaluated travel time functions (TTFs) and the running time. All values are arithmetic means over 100,000 queries executed in bulk with source, target and departure times drawn uniformly at random.

		Queue Pops	Evaluated TTFs	Running Time [ms]
Ber	Basic	167.4	100,820.5	8.8
	+ Corridor	38.1	5224.1	0.6
	+ Lazy	1603.6	1747.4	0.6
	+ A*	635.2	691.5	0.3
Ger06	Basic	492.3	818,721.3	46.4
	+ Corridor	79.7	31,740.8	2.3
	+ Lazy	3323.2	3838.0	1.7
	+ A*	831.0	995.1	0.6
Ger17	Basic	510.3	2,100,731.8	169.7
	+ Corridor	143.4	164,372.5	13.7
	+ Lazy	18,450.0	19,910.5	9.1
	+ A*	3,099.2	3495.5	1.7
Eur17	Basic	861.6	9,951,623.1	808.6
	+ Corridor	229.3	806,727.8	62.3
	+ Lazy	39,714.8	43,581.1	20.8
	+ A*	6876.5	7911.0	4.1
Eur20	Basic	871.0	10,527,072.7	813.2
	+ Corridor	335.6	1,222,655.6	92.9
	+ Lazy	62,677.7	70,145.4	33.7
	+ A*	7231.9	8844.7	4.7

depicts the influence of the query optimizations. The basic approach does not achieve competitive running times. Queries take almost a second on average on the newer Europe graphs. Naively evaluating shortcut travel times with Algorithm 2 is too slow.

Limiting the search space to a shortest path corridor using the elimination tree interval query significantly reduces running times. The speedup is between a factor of 20 on Ger06 and 8 and Eur20. The effectiveness of this optimization corresponds inversely to the relative delay (see Table 1). Greater relative delays lead to bigger corridors and thus smaller speedups.

The lazy evaluation optimization has limited impact on the running time (speedups between 1.3 and 3). However, it drastically shifts the balance between queue operations and travel time function evaluations. On the newer graphs, the number of queue pops increases by more than two orders of magnitude, while the number of travel time function evaluations decreases by up to a factor of 20. The additional queue operations introduce some overhead. However, this is mitigated by the avoided unnecessary and duplicate evaluations. Reusing the lower bounds from the corridor search for goal directed search yields an additional speedup of factor two to seven.

3.4.1. Local Queries

We generate another set of queries to investigate the performance of our algorithms depending on the distance of source and target. We draw 10,000 start nodes and departure times uniformly at random and perform time-dependent Dijkstra without a specific target. For every $2^i$th settled node, we store it as the target of a query of Dijkstra rank i. This methodology was introduced by Sanders and Schultes [33]. Figure 11 shows query running times by rank for the query algorithm with the various query optimizations enabled successively.

Obviously, query running times scale with the distance. The fully optimized algorithm takes only fractions of milliseconds for short range queries, except for some outliers which take up to a millisecond. For long range queries, we usually achieve query times within a couple of milliseconds and the maximum query time was 39 ms. The basic query algorithm is around two orders of magnitude slower across all ranks. The impact of lazy optimization appears to depend on the rank of the query. For lower ranks, it introduces some overhead but reduces outliers compared to only the corridor optimization. This is due to the overhead of the queue operations. For long range queries, this is completely amortized by the reduction in arc relaxations.

3.4.2. Profile Queries

We perform experiments for profile queries and report the results in

Table 4. Running times of profile queries and characteristics of the obtained profiles. We report the total running time and the running time of each phase (Corridor, Reconstruction, Contraction, Extraction) of the query. The total running time is slightly larger than the sum of all phases as it includes some additional initialization and cleanup work. In addition, we report the number of breakpoints in the obtained travel time profile (Column $\left|f\right|$). Column $\left|X\right|$ contains the number of times the shortest path changes during the day. Since the same path may be the fastest for several times, we also report the number of distinct paths in the last column. All values are arithmetic means over 1000 queries executed in bulk with source and target nodes drawn uniformly at random.

	Running Time [ms]
	I	II	III	IV		Total	|f|	|X|	Distinct
				TTF	Paths				Paths
Ber	0.1	37.8	13.9	2.8	0.4	55.6	30,974.9	2.7	2.3
Ger06	0.4	56.5	23.0	0.7	0.8	83.6	9359.6	6.9	3.3
Ger17	0.8	452.2	189.0	6.1	2.4	660.1	66,146.0	9.7	3.7
Eur17	1.7	1135.8	732.9	13.0	7.8	1913.2	122,192.1	16.5	6.8
Eur20	2.8	3166.6	1507.7	12.3	10.3	4747.5	107,690.7	24.4	11.6

. The total running time depends on the amount of time-dependent information in the instance. From Ger06 to Ger17, the total running time increases by a factor of eight even though the network grows only little. The same can be observed between Eur17 and Eur20.

The total running time is dominated by the reconstruction and contraction phases. Reconstruction of the travel time functions of the existing shortcuts takes roughly twice as long as computing the functions for the new shortcuts in the contraction phase. The corridor phase with the elimination tree interval query takes a negligible amount of time. Surprisingly, the time required to compute a final exact travel time profile is not much greater than to compute a path profile (on Ger06, the path profile is even slower than the travel time profile).

The average complexity of the final travel time function varies by an order of magnitude across the different instances. This confirms that Ger06 is a relatively simple instance. Surprisingly, the average complexity on Eur17 is higher than on Eur20. We suspect that this is caused by the higher average complexity in the input graph and that many important arcs that cover many shortest paths already have a non-constant travel time function in Eur17. The number of path switches is magnitudes smaller. It ranges from 3 path switches on Ber to 24 on Eur20. The number of distinct paths during the day are roughly half of that. The arithmetic mean of these numbers is slightly skewed upwards by a few very high values. However, the median is still fairly close: For example, the median number of distinct paths is two on Ger06 and ten on Eur20.

3.5. Comparison with Related Work

Table 5. Comparison with related work. We list unscaled numbers as reported in the respective publications for algorithms we could not run ourselves. Values not reported are indicated as n/r. OOM means that the program crashed while trying to allocate more memory than available. A similar overview with scaled numbers can be found in [34].

		Preprocessing		Index	Query
		Time	Cores	Size	Time	Rel. Error
		[s]		[GB]	[ms]	Avg. [%]	Max. [%]
Ger06	TD-Dijkstra	-	-	-	719.26	-	-
	TDCALT [4]	540	1	0.23	5.36	-	-
	TDCALT-K1.15 [4]	540	1	0.23	1.87	0.050	13.840
	eco L-SHARC [14]	4680	1	1.03	6.31	-	-
	heu SHARC [14]	12,360	1	0.64	0.69	n/r	0.610
	KaTCH	169	16	4.66	0.64	-	-
	TCH [8]	378	8	4.66	0.75	-	-
	ATCH (1.0) [8]	378	8	1.12	1.24	-	-
	ATCH (∞) [8]	378	8	0.55	1.66	-	-
	inex. TCH (0.1) [8]	378	8	1.34	0.70	0.020	0.100
	inex. TCH (1.0) [8]	378	8	1.00	0.69	0.270	1.010
	TD-CRP (0.1) [9]	289	16	0.78	1.92	0.050	0.250
	TD-CRP (1.0) [9]	281	16	0.36	1.66	0.680	2.850
	FLAT [17]	158,760	6	54.63	1.27	0.015	n/r
	CFLAT [17]	104,220	6	34.63	0.58	0.008	0.918
	TD-S+9	542	1	3.61	2.07	0.001	1.523
	CATCHUp	52	16	1.06	0.72	-	-
Ger17	TD-Dijkstra	-	-	-	814.60	-	-
	KaTCH	859	16	42.81	1.26	-	-
	TD-S+9	601	1	5.28	2.61	0.001	0.963
	CATCHUp	142	16	1.50	2.02	-	-
Eur17	TD-Dijkstra	-	-	-	2929.72	-	-
	KaTCH	3066	16	146.97	OOM	-	-
	TD-S+9	3149	1	18.84	4.70	0.002	1.159
	CATCHUp	747	16	5.47	4.92	-	-
Eur20	TD-Dijkstra	-	-	-	3784.11	-	-
	KaTCH	7149	16	239.78	OOM	-	-
	TD-S+9	3352	1	20.65	4.23	0.006	1.733
	CATCHUp	1249	16	6.32	5.60	-	-

provides an overview over different techniques, their preprocessing and query times, space overhead of the index data structures and average query errors where approximation is used. Where possible, we obtained the code of competing algorithms (KaTCH: https://github.com/GVeitBatz/KaTCH (accessed on 15 March 2021); and TD-S: https://github.com/ben-strasser/td_p (accessed on 15 March 2021 ) and evaluated them with same methodology, instances and queries as our algorithms. For other competitors, we report available numbers from the respective publications.

In our comparison, KaTCH, heu SHARC, CFLAT and CATCHUp all achieve query times around 0.6 ms on Ger06. The original research implementation TCH reports slightly slower times than KaTCH. This may be because experiments were run on an older machine, but also because according to the KaTCH documentation, the newer query is somewhat more efficient. TCH pays for this speed with 4.7 GB index data. Reducing the KaTCH memory consumption while keeping exactness (ATCH) brings query times up to 1.24 ms. ATCH also feature a configuration where they only keep upper and lower bounds for each travel time function (ATCH ∞). This configuration uses even less memory than CATCHUp because the optimized order results in fewer shortcuts. However, query running times degrade to 1.66 ms. Giving up on exactness allows keeping the query times at 0.7 ms (inex. TCH) but introduces noticeable errors.

While achieving competitive query times for acceptable memory consumption, heu SHARC suffers from huge preprocessing times of several hours. The original publication does not report average query errors, only a maximum error of 0.61%. TDCALT has the smallest memory consumption but does not achieve competitive query times, even when approximating. FLAT and CFLAT both suffer from extreme preprocessing times and memory consumption despite having no exact queries. CATCHUp offers competitive query times for exact results while keeping memory consumption reasonable. TD-CRP offers even lower memory consumption. However, this is only possible through the use of approximation. TD-CRP queries depict a noticeable error and perform somewhat worse than KaTCH or CATCHUp queries. TD-S+9 depicts the smallest average error of all non-exact approaches (Kontogiannis et al. [17] reported another CFLAT configuration with even smaller errors but significantly slower queries).

Path retrieval in the time-dependent scenario is not as easy as in the static setting. Table 5 reports running times for the earliest arrival query and the path retrieval combined. We only have separate numbers for KaTCH and CATCHUp. For CFLAT, Kontogiannis et al. [17] reported that path retrieval takes up to a third of the total query time. Our experiments show a similar amount for KaTCH. For CATCHUp, path retrieval takes up less than 10% of the query time. TD-CRP and FLAT do not support path retrieval.

On Ger17, KaTCH query times increase by a factor of about two. However, memory usage grows by almost an order of magnitude. For TD-S, both the growth in space consumption and query times corresponds roughly to the growth of the graph size, but not to the increased number of breakpoints. The index of CATCHUp grows by a similar factor. Query times get about 2.7 times slower.

On Eur17, the memory consumption of KaTCH becomes prohibitive. While KaTCH is still able to finish preprocessing and output 150 GB of data, queries crash since the 192 GB RAM of our machine are not enough. Using ATCH or inexact TCH, the memory consumption could likely be reduced sufficiently to perform queries. However, this would either introduce errors or slow down queries significantly. On the other hand, with only 5.5 GB of index data, CATCHUp is still able to perform exact queries in less than 5 ms on average. This is fast enough to enable interactive applications. Total preprocessing for CATCHUp takes less than a quarter of the time KaTCH needs. TD-S+9 is also able to handle this instance with similar query times but only with a small average error.

On Eur20, the behavior is similar, only more extreme. KaTCH preprocessing time increases by more than a factor of two and index data grows to 240 GB. The TD-S+9 numbers remain relatively stable. Query times get slightly faster but errors become larger. CATCHUp preprocessing times also become slower but by less than a factor of two. Query times increase to 5.6 ms. The index takes only 6.3 GB, which is smaller than the input network.

4. Discussion

We introduce CATCHUp, a speed-up technique for routing in time-dependent road networks. It features a small index size and fast, exact queries. To the best of our knowledge, our approach is the first to simultaneously achieve all three objectives. We perform an extensive experimental study to evaluate the performance of CATCHUp and compare it to competing approaches. Our approach achieves the fastest preprocessing, competitive query running times and up to 38 times smaller indexes than other approaches. This clearly demonstrates the advantage of storing expansion information instead of travel time functions.

Revisiting ATCH, TCH and TD-CRP with the insights gained in this work could be fruitful. Combining ATCH with our A* query extension could reduce ATCH query running times. CATCHUp makes use of travel time independent node orders. Combining CATCHUp with TCH-like node orders could result in even smaller index sizes and query running times. We further expect that some of our optimizations to the preprocessing can also be applied in a TD-CRP context. Another possible direction for future research would be to support partial updates to further accelerate the customization. This could enable the integration of live traffic information.