1 #include "engine/routing_algorithms/many_to_many.hpp"
2 #include "engine/routing_algorithms/routing_base_mld.hpp"
3 
4 #include <boost/assert.hpp>
5 #include <boost/range/iterator_range_core.hpp>
6 
7 #include <limits>
8 #include <memory>
9 #include <unordered_map>
10 #include <vector>
11 
12 namespace osrm
13 {
14 namespace engine
15 {
16 namespace routing_algorithms
17 {
18 
19 namespace mld
20 {
21 
22 using PackedEdge = std::tuple</*from*/ NodeID, /*to*/ NodeID, /*from_clique_arc*/ bool>;
23 using PackedPath = std::vector<PackedEdge>;
24 
25 template <typename MultiLevelPartition>
getNodeQueryLevel(const MultiLevelPartition & partition,const NodeID node,const PhantomNode & phantom_node,const LevelID maximal_level)26 inline LevelID getNodeQueryLevel(const MultiLevelPartition &partition,
27                                  const NodeID node,
28                                  const PhantomNode &phantom_node,
29                                  const LevelID maximal_level)
30 {
31     const auto node_level = getNodeQueryLevel(partition, node, phantom_node);
32 
33     if (node_level >= maximal_level)
34         return INVALID_LEVEL_ID;
35 
36     return node_level;
37 }
38 
39 template <bool DIRECTION>
relaxBorderEdges(const DataFacade<mld::Algorithm> & facade,const NodeID node,const EdgeWeight weight,const EdgeDuration duration,const EdgeDistance distance,typename SearchEngineData<mld::Algorithm>::ManyToManyQueryHeap & query_heap,LevelID level)40 void relaxBorderEdges(const DataFacade<mld::Algorithm> &facade,
41                       const NodeID node,
42                       const EdgeWeight weight,
43                       const EdgeDuration duration,
44                       const EdgeDistance distance,
45                       typename SearchEngineData<mld::Algorithm>::ManyToManyQueryHeap &query_heap,
46                       LevelID level)
47 {
48     for (const auto edge : facade.GetBorderEdgeRange(level, node))
49     {
50         const auto &data = facade.GetEdgeData(edge);
51         if ((DIRECTION == FORWARD_DIRECTION) ? facade.IsForwardEdge(edge)
52                                              : facade.IsBackwardEdge(edge))
53         {
54             const NodeID to = facade.GetTarget(edge);
55             if (facade.ExcludeNode(to))
56             {
57                 continue;
58             }
59 
60             const auto turn_id = data.turn_id;
61             const auto node_id = DIRECTION == FORWARD_DIRECTION ? node : facade.GetTarget(edge);
62             const auto node_weight = facade.GetNodeWeight(node_id);
63             const auto node_duration = facade.GetNodeDuration(node_id);
64             const auto node_distance = facade.GetNodeDistance(node_id);
65             const auto turn_weight = node_weight + facade.GetWeightPenaltyForEdgeID(turn_id);
66             const auto turn_duration = node_duration + facade.GetDurationPenaltyForEdgeID(turn_id);
67 
68             BOOST_ASSERT_MSG(node_weight + turn_weight > 0, "edge weight is invalid");
69             const auto to_weight = weight + turn_weight;
70             const auto to_duration = duration + turn_duration;
71             const auto to_distance = distance + node_distance;
72 
73             // New Node discovered -> Add to Heap + Node Info Storage
74             const auto toHeapNode = query_heap.GetHeapNodeIfWasInserted(to);
75             if (!toHeapNode)
76             {
77                 query_heap.Insert(to, to_weight, {node, false, to_duration, to_distance});
78             }
79             // Found a shorter Path -> Update weight and set new parent
80             else if (std::tie(to_weight, to_duration, to_distance, node) <
81                      std::tie(toHeapNode->weight,
82                               toHeapNode->data.duration,
83                               toHeapNode->data.distance,
84                               toHeapNode->data.parent))
85             {
86                 toHeapNode->data = {node, false, to_duration, to_distance};
87                 toHeapNode->weight = to_weight;
88                 query_heap.DecreaseKey(*toHeapNode);
89             }
90         }
91     }
92 }
93 
94 template <bool DIRECTION, typename... Args>
relaxOutgoingEdges(const DataFacade<mld::Algorithm> & facade,const typename SearchEngineData<mld::Algorithm>::ManyToManyQueryHeap::HeapNode & heapNode,typename SearchEngineData<mld::Algorithm>::ManyToManyQueryHeap & query_heap,Args...args)95 void relaxOutgoingEdges(
96     const DataFacade<mld::Algorithm> &facade,
97     const typename SearchEngineData<mld::Algorithm>::ManyToManyQueryHeap::HeapNode &heapNode,
98     typename SearchEngineData<mld::Algorithm>::ManyToManyQueryHeap &query_heap,
99     Args... args)
100 {
101     BOOST_ASSERT(!facade.ExcludeNode(heapNode.node));
102 
103     const auto &partition = facade.GetMultiLevelPartition();
104 
105     const auto level = getNodeQueryLevel(partition, heapNode.node, args...);
106 
107     // Break outgoing edges relaxation if node at the restricted level
108     if (level == INVALID_LEVEL_ID)
109         return;
110 
111     const auto &cells = facade.GetCellStorage();
112     const auto &metric = facade.GetCellMetric();
113 
114     if (level >= 1 && !heapNode.data.from_clique_arc)
115     {
116         const auto &cell = cells.GetCell(metric, level, partition.GetCell(level, heapNode.node));
117         if (DIRECTION == FORWARD_DIRECTION)
118         { // Shortcuts in forward direction
119             auto destination = cell.GetDestinationNodes().begin();
120             auto shortcut_durations = cell.GetOutDuration(heapNode.node);
121             auto shortcut_distances = cell.GetOutDistance(heapNode.node);
122             for (auto shortcut_weight : cell.GetOutWeight(heapNode.node))
123             {
124                 BOOST_ASSERT(destination != cell.GetDestinationNodes().end());
125                 BOOST_ASSERT(!shortcut_durations.empty());
126                 BOOST_ASSERT(!shortcut_distances.empty());
127                 const NodeID to = *destination;
128 
129                 if (shortcut_weight != INVALID_EDGE_WEIGHT && heapNode.node != to)
130                 {
131                     const auto to_weight = heapNode.weight + shortcut_weight;
132                     const auto to_duration = heapNode.data.duration + shortcut_durations.front();
133                     const auto to_distance = heapNode.data.distance + shortcut_distances.front();
134                     const auto toHeapNode = query_heap.GetHeapNodeIfWasInserted(to);
135                     if (!toHeapNode)
136                     {
137                         query_heap.Insert(
138                             to, to_weight, {heapNode.node, true, to_duration, to_distance});
139                     }
140                     else if (std::tie(to_weight, to_duration, to_distance, heapNode.node) <
141                              std::tie(toHeapNode->weight,
142                                       toHeapNode->data.duration,
143                                       toHeapNode->data.distance,
144                                       toHeapNode->data.parent))
145                     {
146                         toHeapNode->data = {heapNode.node, true, to_duration, to_distance};
147                         toHeapNode->weight = to_weight;
148                         query_heap.DecreaseKey(*toHeapNode);
149                     }
150                 }
151                 ++destination;
152                 shortcut_durations.advance_begin(1);
153                 shortcut_distances.advance_begin(1);
154             }
155             BOOST_ASSERT(shortcut_durations.empty());
156             BOOST_ASSERT(shortcut_distances.empty());
157         }
158         else
159         { // Shortcuts in backward direction
160             auto source = cell.GetSourceNodes().begin();
161             auto shortcut_durations = cell.GetInDuration(heapNode.node);
162             auto shortcut_distances = cell.GetInDistance(heapNode.node);
163             for (auto shortcut_weight : cell.GetInWeight(heapNode.node))
164             {
165                 BOOST_ASSERT(source != cell.GetSourceNodes().end());
166                 BOOST_ASSERT(!shortcut_durations.empty());
167                 BOOST_ASSERT(!shortcut_distances.empty());
168                 const NodeID to = *source;
169 
170                 if (shortcut_weight != INVALID_EDGE_WEIGHT && heapNode.node != to)
171                 {
172                     const auto to_weight = heapNode.weight + shortcut_weight;
173                     const auto to_duration = heapNode.data.duration + shortcut_durations.front();
174                     const auto to_distance = heapNode.data.distance + shortcut_distances.front();
175                     const auto toHeapNode = query_heap.GetHeapNodeIfWasInserted(to);
176                     if (!toHeapNode)
177                     {
178                         query_heap.Insert(
179                             to, to_weight, {heapNode.node, true, to_duration, to_distance});
180                     }
181                     else if (std::tie(to_weight, to_duration, to_distance, heapNode.node) <
182                              std::tie(toHeapNode->weight,
183                                       toHeapNode->data.duration,
184                                       toHeapNode->data.distance,
185                                       toHeapNode->data.parent))
186                     {
187                         toHeapNode->data = {heapNode.node, true, to_duration, to_distance};
188                         toHeapNode->weight = to_weight;
189                         query_heap.DecreaseKey(*toHeapNode);
190                     }
191                 }
192                 ++source;
193                 shortcut_durations.advance_begin(1);
194                 shortcut_distances.advance_begin(1);
195             }
196             BOOST_ASSERT(shortcut_durations.empty());
197             BOOST_ASSERT(shortcut_distances.empty());
198         }
199     }
200 
201     relaxBorderEdges<DIRECTION>(facade,
202                                 heapNode.node,
203                                 heapNode.weight,
204                                 heapNode.data.duration,
205                                 heapNode.data.distance,
206                                 query_heap,
207                                 level);
208 }
209 
210 //
211 // Unidirectional multi-layer Dijkstra search for 1-to-N and N-to-1 matrices
212 //
213 template <bool DIRECTION>
214 std::pair<std::vector<EdgeDuration>, std::vector<EdgeDistance>>
oneToManySearch(SearchEngineData<Algorithm> & engine_working_data,const DataFacade<Algorithm> & facade,const std::vector<PhantomNode> & phantom_nodes,std::size_t phantom_index,const std::vector<std::size_t> & phantom_indices,const bool calculate_distance)215 oneToManySearch(SearchEngineData<Algorithm> &engine_working_data,
216                 const DataFacade<Algorithm> &facade,
217                 const std::vector<PhantomNode> &phantom_nodes,
218                 std::size_t phantom_index,
219                 const std::vector<std::size_t> &phantom_indices,
220                 const bool calculate_distance)
221 {
222     std::vector<EdgeWeight> weights_table(phantom_indices.size(), INVALID_EDGE_WEIGHT);
223     std::vector<EdgeDuration> durations_table(phantom_indices.size(), MAXIMAL_EDGE_DURATION);
224     std::vector<EdgeDistance> distances_table(calculate_distance ? phantom_indices.size() : 0,
225                                               MAXIMAL_EDGE_DISTANCE);
226     std::vector<NodeID> middle_nodes_table(phantom_indices.size(), SPECIAL_NODEID);
227 
228     // Collect destination (source) nodes into a map
229     std::unordered_multimap<NodeID, std::tuple<std::size_t, EdgeWeight, EdgeDuration, EdgeDistance>>
230         target_nodes_index;
231     target_nodes_index.reserve(phantom_indices.size());
232     for (std::size_t index = 0; index < phantom_indices.size(); ++index)
233     {
234         const auto &phantom_index = phantom_indices[index];
235         const auto &phantom_node = phantom_nodes[phantom_index];
236 
237         if (DIRECTION == FORWARD_DIRECTION)
238         {
239             if (phantom_node.IsValidForwardTarget())
240                 target_nodes_index.insert(
241                     {phantom_node.forward_segment_id.id,
242                      std::make_tuple(index,
243                                      phantom_node.GetForwardWeightPlusOffset(),
244                                      phantom_node.GetForwardDuration(),
245                                      phantom_node.GetForwardDistance())});
246             if (phantom_node.IsValidReverseTarget())
247                 target_nodes_index.insert(
248                     {phantom_node.reverse_segment_id.id,
249                      std::make_tuple(index,
250                                      phantom_node.GetReverseWeightPlusOffset(),
251                                      phantom_node.GetReverseDuration(),
252                                      phantom_node.GetReverseDistance())});
253         }
254         else if (DIRECTION == REVERSE_DIRECTION)
255         {
256             if (phantom_node.IsValidForwardSource())
257                 target_nodes_index.insert(
258                     {phantom_node.forward_segment_id.id,
259                      std::make_tuple(index,
260                                      -phantom_node.GetForwardWeightPlusOffset(),
261                                      -phantom_node.GetForwardDuration(),
262                                      -phantom_node.GetForwardDistance())});
263             if (phantom_node.IsValidReverseSource())
264                 target_nodes_index.insert(
265                     {phantom_node.reverse_segment_id.id,
266                      std::make_tuple(index,
267                                      -phantom_node.GetReverseWeightPlusOffset(),
268                                      -phantom_node.GetReverseDuration(),
269                                      -phantom_node.GetReverseDistance())});
270         }
271     }
272 
273     // Initialize query heap
274     engine_working_data.InitializeOrClearManyToManyThreadLocalStorage(
275         facade.GetNumberOfNodes(), facade.GetMaxBorderNodeID() + 1);
276     auto &query_heap = *(engine_working_data.many_to_many_heap);
277 
278     // Check if node is in the destinations list and update weights/durations
279     auto update_values =
280         [&](NodeID node, EdgeWeight weight, EdgeDuration duration, EdgeDistance distance) {
281             auto candidates = target_nodes_index.equal_range(node);
282             for (auto it = candidates.first; it != candidates.second;)
283             {
284                 std::size_t index;
285                 EdgeWeight target_weight;
286                 EdgeDuration target_duration;
287                 EdgeDistance target_distance;
288                 std::tie(index, target_weight, target_duration, target_distance) = it->second;
289 
290                 const auto path_weight = weight + target_weight;
291                 if (path_weight >= 0)
292                 {
293                     const auto path_duration = duration + target_duration;
294                     const auto path_distance = distance + target_distance;
295 
296                     EdgeDistance nulldistance = 0;
297                     auto &current_distance =
298                         distances_table.empty() ? nulldistance : distances_table[index];
299 
300                     if (std::tie(path_weight, path_duration, path_distance) <
301                         std::tie(weights_table[index], durations_table[index], current_distance))
302                     {
303                         weights_table[index] = path_weight;
304                         durations_table[index] = path_duration;
305                         current_distance = path_distance;
306                         middle_nodes_table[index] = node;
307                     }
308 
309                     // Remove node from destinations list
310                     it = target_nodes_index.erase(it);
311                 }
312                 else
313                 {
314                     ++it;
315                 }
316             }
317         };
318 
319     auto insert_node = [&](NodeID node,
320                            EdgeWeight initial_weight,
321                            EdgeDuration initial_duration,
322                            EdgeDistance initial_distance) {
323         if (target_nodes_index.count(node))
324         {
325             // Source and target on the same edge node. If target is not reachable directly via
326             // the node (e.g destination is before source on oneway segment) we want to allow
327             // node to be visited later in the search along a reachable path.
328             // Therefore, we manually run first step of search without marking node as visited.
329             update_values(node, initial_weight, initial_duration, initial_distance);
330             relaxBorderEdges<DIRECTION>(
331                 facade, node, initial_weight, initial_duration, initial_distance, query_heap, 0);
332         }
333         else
334         {
335             query_heap.Insert(node, initial_weight, {node, initial_duration, initial_distance});
336         }
337     };
338 
339     { // Place source (destination) adjacent nodes into the heap
340         const auto &phantom_node = phantom_nodes[phantom_index];
341 
342         if (DIRECTION == FORWARD_DIRECTION)
343         {
344             if (phantom_node.IsValidForwardSource())
345             {
346                 insert_node(phantom_node.forward_segment_id.id,
347                             -phantom_node.GetForwardWeightPlusOffset(),
348                             -phantom_node.GetForwardDuration(),
349                             -phantom_node.GetForwardDistance());
350             }
351 
352             if (phantom_node.IsValidReverseSource())
353             {
354                 insert_node(phantom_node.reverse_segment_id.id,
355                             -phantom_node.GetReverseWeightPlusOffset(),
356                             -phantom_node.GetReverseDuration(),
357                             -phantom_node.GetReverseDistance());
358             }
359         }
360         else if (DIRECTION == REVERSE_DIRECTION)
361         {
362             if (phantom_node.IsValidForwardTarget())
363             {
364                 insert_node(phantom_node.forward_segment_id.id,
365                             phantom_node.GetForwardWeightPlusOffset(),
366                             phantom_node.GetForwardDuration(),
367                             phantom_node.GetForwardDistance());
368             }
369 
370             if (phantom_node.IsValidReverseTarget())
371             {
372                 insert_node(phantom_node.reverse_segment_id.id,
373                             phantom_node.GetReverseWeightPlusOffset(),
374                             phantom_node.GetReverseDuration(),
375                             phantom_node.GetReverseDistance());
376             }
377         }
378     }
379 
380     while (!query_heap.Empty() && !target_nodes_index.empty())
381     {
382         // Extract node from the heap. Take a copy (no ref) because otherwise can be modified later
383         // if toHeapNode is the same
384         const auto heapNode = query_heap.DeleteMinGetHeapNode();
385 
386         // Update values
387         update_values(
388             heapNode.node, heapNode.weight, heapNode.data.duration, heapNode.data.distance);
389 
390         // Relax outgoing edges
391         relaxOutgoingEdges<DIRECTION>(
392             facade, heapNode, query_heap, phantom_nodes, phantom_index, phantom_indices);
393     }
394 
395     return std::make_pair(std::move(durations_table), std::move(distances_table));
396 }
397 
398 //
399 // Bidirectional multi-layer Dijkstra search for M-to-N matrices
400 //
401 template <bool DIRECTION>
forwardRoutingStep(const DataFacade<Algorithm> & facade,const unsigned row_idx,const unsigned number_of_sources,const unsigned number_of_targets,typename SearchEngineData<Algorithm>::ManyToManyQueryHeap & query_heap,const std::vector<NodeBucket> & search_space_with_buckets,std::vector<EdgeWeight> & weights_table,std::vector<EdgeDuration> & durations_table,std::vector<EdgeDistance> & distances_table,std::vector<NodeID> & middle_nodes_table,const PhantomNode & phantom_node)402 void forwardRoutingStep(const DataFacade<Algorithm> &facade,
403                         const unsigned row_idx,
404                         const unsigned number_of_sources,
405                         const unsigned number_of_targets,
406                         typename SearchEngineData<Algorithm>::ManyToManyQueryHeap &query_heap,
407                         const std::vector<NodeBucket> &search_space_with_buckets,
408                         std::vector<EdgeWeight> &weights_table,
409                         std::vector<EdgeDuration> &durations_table,
410                         std::vector<EdgeDistance> &distances_table,
411                         std::vector<NodeID> &middle_nodes_table,
412                         const PhantomNode &phantom_node)
413 {
414     // Take a copy of the extracted node because otherwise could be modified later if toHeapNode is
415     // the same
416     const auto heapNode = query_heap.DeleteMinGetHeapNode();
417 
418     // Check if each encountered node has an entry
419     const auto &bucket_list = std::equal_range(search_space_with_buckets.begin(),
420                                                search_space_with_buckets.end(),
421                                                heapNode.node,
422                                                NodeBucket::Compare());
423     for (const auto &current_bucket : boost::make_iterator_range(bucket_list))
424     {
425         // Get target id from bucket entry
426         const auto column_idx = current_bucket.column_index;
427         const auto target_weight = current_bucket.weight;
428         const auto target_duration = current_bucket.duration;
429         const auto target_distance = current_bucket.distance;
430 
431         // Get the value location in the results tables:
432         //  * row-major direct (row_idx, column_idx) index for forward direction
433         //  * row-major transposed (column_idx, row_idx) for reversed direction
434         const auto location = DIRECTION == FORWARD_DIRECTION
435                                   ? row_idx * number_of_targets + column_idx
436                                   : row_idx + column_idx * number_of_sources;
437         auto &current_weight = weights_table[location];
438         auto &current_duration = durations_table[location];
439 
440         EdgeDistance nulldistance = 0;
441         auto &current_distance = distances_table.empty() ? nulldistance : distances_table[location];
442 
443         // Check if new weight is better
444         auto new_weight = heapNode.weight + target_weight;
445         auto new_duration = heapNode.data.duration + target_duration;
446         auto new_distance = heapNode.data.distance + target_distance;
447 
448         if (new_weight >= 0 && std::tie(new_weight, new_duration, new_distance) <
449                                    std::tie(current_weight, current_duration, current_distance))
450         {
451             current_weight = new_weight;
452             current_duration = new_duration;
453             current_distance = new_distance;
454             middle_nodes_table[location] = heapNode.node;
455         }
456     }
457 
458     relaxOutgoingEdges<DIRECTION>(facade, heapNode, query_heap, phantom_node);
459 }
460 
461 template <bool DIRECTION>
backwardRoutingStep(const DataFacade<Algorithm> & facade,const unsigned column_idx,typename SearchEngineData<Algorithm>::ManyToManyQueryHeap & query_heap,std::vector<NodeBucket> & search_space_with_buckets,const PhantomNode & phantom_node)462 void backwardRoutingStep(const DataFacade<Algorithm> &facade,
463                          const unsigned column_idx,
464                          typename SearchEngineData<Algorithm>::ManyToManyQueryHeap &query_heap,
465                          std::vector<NodeBucket> &search_space_with_buckets,
466                          const PhantomNode &phantom_node)
467 {
468     // Take a copy of the extracted node because otherwise could be modified later if toHeapNode is
469     // the same
470     const auto heapNode = query_heap.DeleteMinGetHeapNode();
471 
472     // Store settled nodes in search space bucket
473     search_space_with_buckets.emplace_back(heapNode.node,
474                                            heapNode.data.parent,
475                                            heapNode.data.from_clique_arc,
476                                            column_idx,
477                                            heapNode.weight,
478                                            heapNode.data.duration,
479                                            heapNode.data.distance);
480 
481     const auto &partition = facade.GetMultiLevelPartition();
482     const auto maximal_level = partition.GetNumberOfLevels() - 1;
483 
484     relaxOutgoingEdges<!DIRECTION>(facade, heapNode, query_heap, phantom_node, maximal_level);
485 }
486 
487 template <bool DIRECTION>
retrievePackedPathFromSearchSpace(NodeID middle_node_id,const unsigned column_idx,const std::vector<NodeBucket> & search_space_with_buckets,PackedPath & path)488 void retrievePackedPathFromSearchSpace(NodeID middle_node_id,
489                                        const unsigned column_idx,
490                                        const std::vector<NodeBucket> &search_space_with_buckets,
491                                        PackedPath &path)
492 {
493     auto bucket_list = std::equal_range(search_space_with_buckets.begin(),
494                                         search_space_with_buckets.end(),
495                                         middle_node_id,
496                                         NodeBucket::ColumnCompare(column_idx));
497 
498     BOOST_ASSERT_MSG(std::distance(bucket_list.first, bucket_list.second) == 1,
499                      "The pointers are not pointing to the same element.");
500 
501     NodeID current_node_id = middle_node_id;
502 
503     while (bucket_list.first->parent_node != current_node_id &&
504            bucket_list.first != search_space_with_buckets.end())
505     {
506         const auto parent_node_id = bucket_list.first->parent_node;
507 
508         const auto from = DIRECTION == FORWARD_DIRECTION ? current_node_id : parent_node_id;
509         const auto to = DIRECTION == FORWARD_DIRECTION ? parent_node_id : current_node_id;
510         path.emplace_back(std::make_tuple(from, to, bucket_list.first->from_clique_arc));
511 
512         current_node_id = parent_node_id;
513         bucket_list = std::equal_range(search_space_with_buckets.begin(),
514                                        search_space_with_buckets.end(),
515                                        current_node_id,
516                                        NodeBucket::ColumnCompare(column_idx));
517 
518         BOOST_ASSERT_MSG(std::distance(bucket_list.first, bucket_list.second) == 1,
519                          "The pointers are not pointing to the same element.");
520     }
521 }
522 
523 template <bool DIRECTION>
524 std::pair<std::vector<EdgeDuration>, std::vector<EdgeDistance>>
manyToManySearch(SearchEngineData<Algorithm> & engine_working_data,const DataFacade<Algorithm> & facade,const std::vector<PhantomNode> & phantom_nodes,const std::vector<std::size_t> & source_indices,const std::vector<std::size_t> & target_indices,const bool calculate_distance)525 manyToManySearch(SearchEngineData<Algorithm> &engine_working_data,
526                  const DataFacade<Algorithm> &facade,
527                  const std::vector<PhantomNode> &phantom_nodes,
528                  const std::vector<std::size_t> &source_indices,
529                  const std::vector<std::size_t> &target_indices,
530                  const bool calculate_distance)
531 {
532     const auto number_of_sources = source_indices.size();
533     const auto number_of_targets = target_indices.size();
534     const auto number_of_entries = number_of_sources * number_of_targets;
535 
536     std::vector<EdgeWeight> weights_table(number_of_entries, INVALID_EDGE_WEIGHT);
537     std::vector<EdgeDuration> durations_table(number_of_entries, MAXIMAL_EDGE_DURATION);
538     std::vector<EdgeDistance> distances_table(calculate_distance ? number_of_entries : 0,
539                                               INVALID_EDGE_DISTANCE);
540     std::vector<NodeID> middle_nodes_table(number_of_entries, SPECIAL_NODEID);
541 
542     std::vector<NodeBucket> search_space_with_buckets;
543 
544     // Populate buckets with paths from all accessible nodes to destinations via backward searches
545     for (std::uint32_t column_idx = 0; column_idx < target_indices.size(); ++column_idx)
546     {
547         const auto index = target_indices[column_idx];
548         const auto &target_phantom = phantom_nodes[index];
549 
550         engine_working_data.InitializeOrClearManyToManyThreadLocalStorage(
551             facade.GetNumberOfNodes(), facade.GetMaxBorderNodeID() + 1);
552         auto &query_heap = *(engine_working_data.many_to_many_heap);
553 
554         if (DIRECTION == FORWARD_DIRECTION)
555             insertTargetInHeap(query_heap, target_phantom);
556         else
557             insertSourceInHeap(query_heap, target_phantom);
558 
559         // explore search space
560         while (!query_heap.Empty())
561         {
562             backwardRoutingStep<DIRECTION>(
563                 facade, column_idx, query_heap, search_space_with_buckets, target_phantom);
564         }
565     }
566 
567     // Order lookup buckets
568     std::sort(search_space_with_buckets.begin(), search_space_with_buckets.end());
569 
570     // Find shortest paths from sources to all accessible nodes
571     for (std::uint32_t row_idx = 0; row_idx < source_indices.size(); ++row_idx)
572     {
573         const auto source_index = source_indices[row_idx];
574         const auto &source_phantom = phantom_nodes[source_index];
575 
576         // Clear heap and insert source nodes
577         engine_working_data.InitializeOrClearManyToManyThreadLocalStorage(
578             facade.GetNumberOfNodes(), facade.GetMaxBorderNodeID() + 1);
579 
580         auto &query_heap = *(engine_working_data.many_to_many_heap);
581 
582         if (DIRECTION == FORWARD_DIRECTION)
583             insertSourceInHeap(query_heap, source_phantom);
584         else
585             insertTargetInHeap(query_heap, source_phantom);
586 
587         // Explore search space
588         while (!query_heap.Empty())
589         {
590             forwardRoutingStep<DIRECTION>(facade,
591                                           row_idx,
592                                           number_of_sources,
593                                           number_of_targets,
594                                           query_heap,
595                                           search_space_with_buckets,
596                                           weights_table,
597                                           durations_table,
598                                           distances_table,
599                                           middle_nodes_table,
600                                           source_phantom);
601         }
602     }
603 
604     return std::make_pair(std::move(durations_table), std::move(distances_table));
605 }
606 
607 } // namespace mld
608 
609 // Dispatcher function for one-to-many and many-to-one tasks that can be handled by MLD differently:
610 //
611 // * one-to-many (many-to-one) tasks use a unidirectional forward (backward) Dijkstra search
612 //   with the candidate node level `min(GetQueryLevel(phantom_node, node, phantom_nodes)`
613 //   for all destination (source) phantom nodes
614 //
615 // * many-to-many search tasks use a bidirectional Dijkstra search
616 //   with the candidate node level `min(GetHighestDifferentLevel(phantom_node, node))`
617 //   Due to pruned backward search space it is always better to compute the durations matrix
618 //   when number of sources is less than targets. If number of targets is less than sources
619 //   then search is performed on a reversed graph with phantom nodes with flipped roles and
620 //   returning a transposed matrix.
621 template <>
622 std::pair<std::vector<EdgeDuration>, std::vector<EdgeDistance>>
manyToManySearch(SearchEngineData<mld::Algorithm> & engine_working_data,const DataFacade<mld::Algorithm> & facade,const std::vector<PhantomNode> & phantom_nodes,const std::vector<std::size_t> & source_indices,const std::vector<std::size_t> & target_indices,const bool calculate_distance)623 manyToManySearch(SearchEngineData<mld::Algorithm> &engine_working_data,
624                  const DataFacade<mld::Algorithm> &facade,
625                  const std::vector<PhantomNode> &phantom_nodes,
626                  const std::vector<std::size_t> &source_indices,
627                  const std::vector<std::size_t> &target_indices,
628                  const bool calculate_distance)
629 {
630     if (source_indices.size() == 1)
631     { // TODO: check if target_indices.size() == 1 and do a bi-directional search
632         return mld::oneToManySearch<FORWARD_DIRECTION>(engine_working_data,
633                                                        facade,
634                                                        phantom_nodes,
635                                                        source_indices.front(),
636                                                        target_indices,
637                                                        calculate_distance);
638     }
639 
640     if (target_indices.size() == 1)
641     {
642         return mld::oneToManySearch<REVERSE_DIRECTION>(engine_working_data,
643                                                        facade,
644                                                        phantom_nodes,
645                                                        target_indices.front(),
646                                                        source_indices,
647                                                        calculate_distance);
648     }
649 
650     if (target_indices.size() < source_indices.size())
651     {
652         return mld::manyToManySearch<REVERSE_DIRECTION>(engine_working_data,
653                                                         facade,
654                                                         phantom_nodes,
655                                                         target_indices,
656                                                         source_indices,
657                                                         calculate_distance);
658     }
659 
660     return mld::manyToManySearch<FORWARD_DIRECTION>(engine_working_data,
661                                                     facade,
662                                                     phantom_nodes,
663                                                     source_indices,
664                                                     target_indices,
665                                                     calculate_distance);
666 }
667 
668 } // namespace routing_algorithms
669 } // namespace engine
670 } // namespace osrm
671