OmniSciDB  c1a53651b2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GraphTableFunctions.hpp
Go to the documentation of this file.
1 /*
2  * Copyright 2022 HEAVY.AI, Inc., Inc.
3  *
4  * Licensed under the Apache License, Version 2.0 (the "License");
5  * you may not use this file except in compliance with the License.
6  * You may obtain a copy of the License at
7  *
8  * http://www.apache.org/licenses/LICENSE-2.0
9  *
10  * Unless required by applicable law or agreed to in writing, software
11  * distributed under the License is distributed on an "AS IS" BASIS,
12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13  * See the License for the specific language governing permissions and
14  * limitations under the License.
15  */
16 
17 #pragma once
18 
19 #ifndef __CUDACC__
20 #ifdef HAVE_TBB
21 
22 #include <tbb/parallel_for.h>
23 
24 #include <boost/config.hpp>
25 #include <boost/graph/adjacency_list.hpp>
26 #include <boost/graph/dijkstra_shortest_paths.hpp>
27 #include <boost/graph/graph_traits.hpp>
28 #include <boost/property_map/property_map.hpp>
29 
30 #include "ThirdParty/robin_hood/robin_hood.h"
31 
33 
34 template <class T>
35 struct get_mapped;
36 
37 template <>
38 struct get_mapped<int32_t> {
39  using type = int32_t;
40 };
41 
42 template <>
43 struct get_mapped<int64_t> {
44  using type = int64_t;
45 };
46 
47 template <>
48 struct get_mapped<TextEncodingDict> {
49  using type = int32_t;
50 };
51 
52 template <typename T>
53 struct AttrIdxMap {
54  using T2 = typename get_mapped<T>::type;
55  bool must_be_unique;
56  std::vector<T2> idx_to_attr_map;
57  using AttrMap = robin_hood::unordered_flat_map<T2, int32_t>;
58  // std::map<T, int32_t> attr_to_idx_map;
59  AttrMap attr_to_idx_map;
60  int32_t unique_idx = 0;
61 
62  inline int32_t size() const { return static_cast<int64_t>(idx_to_attr_map.size()); }
63 
64  inline int32_t get_idx_for_key(const T2& key) const {
65  // Note this is not safe if key doesn't exist
66  return attr_to_idx_map.at(key);
67  }
68 
69  inline int64_t get_idx_for_key_safe(const T2& key) const {
70  // Note this is not safe if key doesn't exist
71  const auto key_itr = attr_to_idx_map.find(key);
72  if (key_itr == attr_to_idx_map.end()) {
73  return -1; // Invalid sentinel
74  }
75  return key_itr->second;
76  }
77 
78  inline T get_key_for_idx(const int32_t idx) const {
79  // Note this is not safe if key doesn't exist
80  return idx_to_attr_map[idx];
81  }
82 
83  AttrIdxMap(const bool must_be_unique) : must_be_unique(must_be_unique) {}
84 
85  void add_column(const Column<T>& col) {
86  auto timer = DEBUG_TIMER(__func__);
87  const int64_t num_rows = col.size();
88  if (idx_to_attr_map.empty()) {
89  // Only reserve space for the first batch, with the logic being
90  // that most of the unique values will be seen in the first
91  // column, excepting extreme edge cases
92  idx_to_attr_map.reserve(num_rows);
93  attr_to_idx_map.reserve(num_rows);
94  }
95  for (int64_t row_idx = 0; row_idx < num_rows; ++row_idx) {
96  if (!col.isNull(row_idx)) {
97  const bool is_new_elem = attr_to_idx_map.emplace(col[row_idx], unique_idx).second;
98  if (is_new_elem) {
99  idx_to_attr_map.emplace_back(col[row_idx]);
100  unique_idx++;
101  } else {
102  if (must_be_unique) {
103  throw std::runtime_error(
104  "Expected unique elements but found duplicates for column.");
105  }
106  }
107  }
108  }
109  CHECK_EQ(attr_to_idx_map.size(), idx_to_attr_map.size());
110  CHECK_EQ(idx_to_attr_map.size(), static_cast<size_t>(unique_idx));
111  }
112 };
113 
114 template <typename N>
115 std::vector<std::pair<N, N>> construct_key_normalized_edge_list(
116  const Column<N>& node1,
117  const Column<N>& node2,
118  const AttrIdxMap<N>& attr_idx_map) {
119  auto timer = DEBUG_TIMER(__func__);
120  const int32_t num_edges = node1.size();
121  std::vector<std::pair<N, N>> key_normalized_edge_list(num_edges);
123  tbb::blocked_range<int32_t>(0, num_edges),
124  [&](const tbb::blocked_range<int32_t>& r) {
125  const int32_t r_end = r.end();
126  for (int32_t e = r.begin(); e < r_end; ++e) {
127  key_normalized_edge_list[e].first = attr_idx_map.get_idx_for_key(node1[e]);
128  key_normalized_edge_list[e].second = attr_idx_map.get_idx_for_key(node2[e]);
129  }
130  });
131  return key_normalized_edge_list;
132 }
133 
134 template <typename S>
135 struct TerminalNodes {
136  S start_node;
137  S end_node;
138  bool end_node_is_valid;
139 
140  // In the first constructor we set end_node to start_node to squash a
141  // maybe_uninitialized warning (harmless as we have the end_node_is_valid sentinel)
142  TerminalNodes(const S start_node)
143  : start_node(start_node), end_node(start_node), end_node_is_valid(false) {}
144  TerminalNodes(const S start_node, const S end_node)
145  : start_node(start_node), end_node(end_node), end_node_is_valid(true) {}
146 };
147 
148 template <typename N, typename D>
149 struct GraphTraversalResults {
150  using graph_t = boost::adjacency_list<boost::listS,
151  boost::vecS,
152  boost::directedS,
153  boost::no_property,
154  boost::property<boost::edge_weight_t, int32_t>>;
155  using vertex_descriptor = boost::graph_traits<graph_t>::vertex_descriptor;
156  AttrIdxMap<N> attr_idx_map;
157  int32_t start_node_idx, end_node_idx;
158  graph_t edge_graph;
159  std::vector<vertex_descriptor> parents;
160  std::vector<D> graph_distances;
161  int64_t num_vertices{0};
162  boost::graph_traits<graph_t>::vertex_iterator vi, vend;
163 
164  GraphTraversalResults() : attr_idx_map(false) {}
165 
166  void buildEdgeGraph(const std::vector<std::pair<N, N>>& edge_list,
167  const Column<D>& distance) {
168  auto timer = DEBUG_TIMER(__func__);
169  edge_graph = graph_t(edge_list.data(),
170  edge_list.data() + edge_list.size(),
171  distance.ptr_,
172  attr_idx_map.size());
173  num_vertices = boost::num_vertices(edge_graph);
174  parents.resize(num_vertices);
175  graph_distances.resize(num_vertices);
176  boost::tie(vi, vend) = boost::vertices(edge_graph);
177  }
178 };
179 
180 template <typename N, typename D, typename S>
181 GraphTraversalResults<N, D> graph_shortest_path_impl(
182  const Column<N>& node1,
183  const Column<N>& node2,
184  const Column<D>& distance,
185  const TerminalNodes<S>& terminal_nodes) {
186  auto func_timer = DEBUG_TIMER(__func__);
187  using graph_t = boost::adjacency_list<boost::listS,
188  boost::vecS,
189  boost::directedS,
190  boost::no_property,
191  boost::property<boost::edge_weight_t, int32_t>>;
192  using vertex_descriptor = boost::graph_traits<graph_t>::vertex_descriptor;
193  const int64_t num_edges = node1.size();
194  auto new_node2 = Column<N>(node2);
195  std::vector<N> translated_node2_data;
196  N start_node_translated;
197  N end_node_translated;
198  if constexpr (std::is_same_v<N, TextEncodingDict>) {
199  const auto node1_sdp = node1.string_dict_proxy_;
200  const auto node2_sdp = node2.string_dict_proxy_;
201  if (node1_sdp->getDictKey() != node2_sdp->getDictKey()) {
202  auto translation_timer = DEBUG_TIMER("Dictionary Translation");
203  const auto translation_map =
204  node2_sdp->buildUnionTranslationMapToOtherProxy(node1_sdp, {});
205  translated_node2_data.resize(num_edges);
206  for (int64_t edge_idx = 0; edge_idx < num_edges; ++edge_idx) {
207  const N original_val = node2[edge_idx];
208  translated_node2_data[edge_idx] =
209  original_val == inline_null_value<N>()
210  ? original_val
211  : static_cast<N>(translation_map[original_val]);
212  }
213  new_node2 = Column<N>(
214  translated_node2_data.data(), node2.num_rows_, node1.string_dict_proxy_);
215  }
216  if constexpr (!std::is_same_v<S, TextEncodingNone>) {
217  throw std::runtime_error(
218  "Starting node should be of text type to match graph nodes.");
219  } else {
220  start_node_translated =
221  node1.string_dict_proxy_->getIdOfString(terminal_nodes.start_node);
222  if (terminal_nodes.end_node_is_valid) {
223  end_node_translated =
224  node1.string_dict_proxy_->getIdOfString(terminal_nodes.end_node);
225  }
226  }
227  if (start_node_translated == StringDictionary::INVALID_STR_ID) {
228  throw std::runtime_error("Starting node not found.");
229  }
230  } else {
231  if constexpr (std::is_same_v<S, TextEncodingNone>) {
232  throw std::runtime_error(
233  "Starting node should be of integer type to match graph nodes.");
234  } else {
235  start_node_translated = terminal_nodes.start_node;
236  end_node_translated = terminal_nodes.end_node;
237  }
238  }
239 
240  GraphTraversalResults<N, D> graph_traversal_results;
241  graph_traversal_results.attr_idx_map.add_column(node1);
242  graph_traversal_results.attr_idx_map.add_column(new_node2);
243 
244  const auto edge_list = construct_key_normalized_edge_list(
245  node1, new_node2, graph_traversal_results.attr_idx_map);
246  graph_traversal_results.buildEdgeGraph(edge_list, distance);
247  graph_traversal_results.start_node_idx =
248  graph_traversal_results.attr_idx_map.get_idx_for_key_safe(start_node_translated);
249  graph_traversal_results.end_node_idx =
250  terminal_nodes.end_node_is_valid
251  ? graph_traversal_results.attr_idx_map.get_idx_for_key_safe(end_node_translated)
252  : -1;
253 
254  if (graph_traversal_results.start_node_idx < 0) {
255  throw std::runtime_error("Starting node not found.");
256  }
257  if (terminal_nodes.end_node_is_valid && graph_traversal_results.end_node_idx < 0) {
258  throw std::runtime_error("End node not found.");
259  }
260  vertex_descriptor start_vertex = boost::vertex(graph_traversal_results.start_node_idx,
261  graph_traversal_results.edge_graph);
262  {
263  auto shortest_paths_timer = DEBUG_TIMER("Djikstra Shortest Paths");
264  boost::dijkstra_shortest_paths(
265  graph_traversal_results.edge_graph,
266  start_vertex,
267  boost::predecessor_map(
268  boost::make_iterator_property_map(
269  graph_traversal_results.parents.begin(),
270  boost::get(boost::vertex_index, graph_traversal_results.edge_graph)))
271  .distance_map(boost::make_iterator_property_map(
272  graph_traversal_results.graph_distances.begin(),
273  boost::get(boost::vertex_index, graph_traversal_results.edge_graph))));
274  }
275  return graph_traversal_results;
276 }
277 
278 // clang-format off
279 /*
280  UDTF: tf_graph_shortest_path__cpu_template(TableFunctionManager,
281  Cursor<Column<N> node1, Column<N> node2, Column<D> distance> edge_list,
282  S origin_node, S destination_node) ->
283  Column<int32_t> path_step, Column<N> node | input_id=args<0>, Column<D> cume_distance,
284  N=[int32_t, int64_t, TextEncodingDict], D=[int32_t, int64_t, float, double],
285  S=[int64_t, TextEncodingNone]
286 */
287 // clang-format on
288 
289 template <typename N, typename D, typename S>
290 int64_t tf_graph_shortest_path__cpu_template(TableFunctionManager& mgr,
291  const Column<N>& node1,
292  const Column<N>& node2,
293  const Column<D>& distance,
294  const S& start_node,
295  const S& end_node,
296  Column<int32_t>& output_path_step,
297  Column<N>& output_node,
298  Column<D>& output_distance) {
299  auto func_timer = DEBUG_TIMER(__func__);
300  GraphTraversalResults<N, D> graph_traversal_results;
301  TerminalNodes<S> terminal_nodes(start_node, end_node);
302  try {
303  graph_traversal_results =
304  graph_shortest_path_impl(node1, node2, distance, terminal_nodes);
305  } catch (std::exception& e) {
306  return mgr.ERROR_MESSAGE(e.what());
307  }
308  auto output_timer = DEBUG_TIMER("Output shortest path results");
309  int32_t current_node_idx = graph_traversal_results.end_node_idx;
310  // Doing this first to determine how large the array is
311  int32_t num_hops = 0;
312  const D total_distance_to_origin_node =
313  graph_traversal_results.graph_distances[current_node_idx];
314  // Per the boost graph docs infinity distances are by default represented
315  // by std::numeric_limits<D>::max()
316  // https://www.boost.org/doc/libs/1_79_0/libs/graph/doc/dijkstra_shortest_paths.html
317  if (total_distance_to_origin_node != std::numeric_limits<D>::max()) {
318  while (current_node_idx != graph_traversal_results.start_node_idx) {
319  current_node_idx = graph_traversal_results.parents[current_node_idx];
320  num_hops++;
321  }
322  mgr.set_output_row_size(num_hops + 1);
323  current_node_idx = graph_traversal_results.end_node_idx;
324  int32_t path_step_idx = num_hops;
325  const auto last_node_key =
326  graph_traversal_results.attr_idx_map.get_key_for_idx(current_node_idx);
327  output_path_step[path_step_idx] = path_step_idx + 1;
328  output_node[path_step_idx] = last_node_key;
329  output_distance[path_step_idx] = total_distance_to_origin_node;
330  while (current_node_idx != graph_traversal_results.start_node_idx) {
331  current_node_idx = graph_traversal_results.parents[current_node_idx];
332  path_step_idx--;
333  output_path_step[path_step_idx] = path_step_idx + 1;
334  output_node[path_step_idx] =
335  graph_traversal_results.attr_idx_map.get_key_for_idx(current_node_idx);
336  output_distance[path_step_idx] =
337  graph_traversal_results.graph_distances[current_node_idx];
338  }
339  return num_hops + 1;
340  } else {
341  mgr.set_output_row_size(0);
342  return 0;
343  }
344 }
345 
346 // clang-format off
347 /*
348  UDTF: tf_graph_shortest_paths_distances__cpu_template(TableFunctionManager,
349  Cursor<Column<N> node1, Column<N> node2, Column<D> distance> edge_list,
350  S origin_node) ->
351  Column<N> origin_node | input_id=args<0>, Column<N> destination_node | input_id=args<0>,
352  Column<D> distance, Column<int32_t> num_edges_traversed, N=[int32_t, int64_t, TextEncodingDict],
353  D=[int32_t, int64_t, float, double], S=[int64_t, TextEncodingNone]
354 */
355 // clang-format on
356 
357 template <typename N, typename D, typename S>
358 int64_t tf_graph_shortest_paths_distances__cpu_template(TableFunctionManager& mgr,
359  const Column<N>& node1,
360  const Column<N>& node2,
361  const Column<D>& distance,
362  const S& start_node,
363  Column<N>& out_node1,
364  Column<N>& out_node2,
365  Column<D>& out_distance,
366  Column<int32_t>& out_num_hops) {
367  auto func_timer = DEBUG_TIMER(__func__);
368  GraphTraversalResults<N, D> graph_traversal_results;
369  TerminalNodes<S> terminal_nodes(start_node);
370  try {
371  graph_traversal_results =
372  graph_shortest_path_impl(node1, node2, distance, terminal_nodes);
373  } catch (std::exception& e) {
374  return mgr.ERROR_MESSAGE(e.what());
375  }
376  auto output_timer = DEBUG_TIMER("Output shortest paths results");
377  mgr.set_output_row_size(graph_traversal_results.num_vertices);
378  const N node1_val = graph_traversal_results.attr_idx_map.get_key_for_idx(
379  graph_traversal_results.start_node_idx);
381  tbb::blocked_range<int32_t>(0, graph_traversal_results.num_vertices),
382  [&](const tbb::blocked_range<int32_t>& r) {
383  const int32_t r_end = r.end();
384  for (int32_t vertex_idx = r.begin(); vertex_idx < r_end; ++vertex_idx) {
385  out_node1[vertex_idx] = node1_val;
386  out_node2[vertex_idx] =
387  graph_traversal_results.attr_idx_map.get_key_for_idx(vertex_idx);
388  out_distance[vertex_idx] = graph_traversal_results.graph_distances[vertex_idx];
389  if (out_distance[vertex_idx] == std::numeric_limits<D>::max()) {
390  out_distance.setNull(vertex_idx);
391  out_num_hops.setNull(vertex_idx);
392  } else {
393  int32_t num_hops = 0;
394  int32_t current_node_idx = vertex_idx;
395  while (current_node_idx != graph_traversal_results.start_node_idx) {
396  current_node_idx = graph_traversal_results.parents[current_node_idx];
397  num_hops++;
398  }
399  out_num_hops[vertex_idx] = num_hops;
400  }
401  }
402  });
403 
404  return graph_traversal_results.num_vertices;
405 }
406 
407 template <typename Z>
408 struct RasterGraphEdges {
409  std::vector<int32_t> bins_1;
410  std::vector<int32_t> bins_2;
411  std::vector<Z> distances;
412 };
413 
414 template <typename Z>
415 struct RasterGraphEdgeGenerator {
416  const Z slope_weight_exponent;
417  const Z slope_pct_max;
418  RasterGraphEdges<Z> raster_graph_edges;
419  std::atomic<int32_t> num_edges{0};
420 
421  RasterGraphEdgeGenerator(const int64_t max_size,
422  const Z slope_weight_exponent,
423  const Z slope_pct_max)
424  : slope_weight_exponent(slope_weight_exponent), slope_pct_max(slope_pct_max) {
425  raster_graph_edges.bins_1.resize(max_size);
426  raster_graph_edges.bins_2.resize(max_size);
427  raster_graph_edges.distances.resize(max_size);
428  }
429 
430  int64_t size() const { return num_edges; }
431 
432  inline void conditionally_write_edge_slope(const std::pair<int64_t, Z>& bin_1_idx_and_z,
433  const std::pair<int64_t, Z>& bin_2_idx_and_z,
434  const Z x_distance) {
435  // Only checks bin2's z val for nullness, so its up to the caller
436  // to enforce that bin1 is not null
437  if (bin_2_idx_and_z.second == inline_null_value<Z>()) {
438  return;
439  }
440  const Z slope = (bin_2_idx_and_z.second - bin_1_idx_and_z.second) / x_distance;
441  if (slope * 100.0 > slope_pct_max) {
442  return;
443  }
444  const int32_t edge_idx = num_edges++;
445  raster_graph_edges.bins_1[edge_idx] = static_cast<int32_t>(bin_1_idx_and_z.first);
446  raster_graph_edges.bins_2[edge_idx] = static_cast<int32_t>(bin_2_idx_and_z.first);
447  const Z abs_slope = 1.0 + abs(slope);
448  raster_graph_edges.distances[edge_idx] =
449  x_distance * pow(abs_slope, slope_weight_exponent);
450  }
451 
452  void trim_to_size() {
453  raster_graph_edges.bins_1.resize(num_edges);
454  raster_graph_edges.bins_2.resize(num_edges);
455  raster_graph_edges.distances.resize(num_edges);
456  }
457 };
458 
459 template <typename T, typename Z>
460 RasterGraphEdges<Z> generate_raster_graph_edges(const GeoRaster<T, Z>& geo_raster,
461  const Z slope_weight_exponent,
462  const Z slope_pct_max) {
463  RasterGraphEdgeGenerator<Z> raster_graph_edge_generator(
464  geo_raster.num_bins_ * 8, slope_weight_exponent, slope_pct_max);
465  const Z x_distance = static_cast<Z>(geo_raster.bin_dim_meters_);
466  const Z x_distance_diag = sqrt(2.0 * x_distance * x_distance);
468  tbb::blocked_range<int64_t>(0, geo_raster.num_y_bins_),
469  [&](const tbb::blocked_range<int64_t>& r) {
470  const int64_t end_y_bin = r.end();
471  for (int64_t y_bin = r.begin(); y_bin != end_y_bin; ++y_bin) {
472  for (int64_t x_bin = 0; x_bin < geo_raster.num_x_bins_; ++x_bin) {
473  const auto bin_1_idx_and_z =
474  geo_raster.get_bin_idx_and_z_val_for_xy_bin(x_bin, y_bin);
475  if (bin_1_idx_and_z.second == inline_null_value<Z>()) {
476  continue;
477  }
478  raster_graph_edge_generator.conditionally_write_edge_slope(
479  bin_1_idx_and_z,
480  geo_raster.get_bin_idx_and_z_val_for_xy_bin(x_bin, y_bin - 1),
481  x_distance);
482  raster_graph_edge_generator.conditionally_write_edge_slope(
483  bin_1_idx_and_z,
484  geo_raster.get_bin_idx_and_z_val_for_xy_bin(x_bin, y_bin + 1),
485  x_distance);
486  raster_graph_edge_generator.conditionally_write_edge_slope(
487  bin_1_idx_and_z,
488  geo_raster.get_bin_idx_and_z_val_for_xy_bin(x_bin - 1, y_bin),
489  x_distance);
490  raster_graph_edge_generator.conditionally_write_edge_slope(
491  bin_1_idx_and_z,
492  geo_raster.get_bin_idx_and_z_val_for_xy_bin(x_bin + 1, y_bin),
493  x_distance);
494  raster_graph_edge_generator.conditionally_write_edge_slope(
495  bin_1_idx_and_z,
496  geo_raster.get_bin_idx_and_z_val_for_xy_bin(x_bin - 1, y_bin - 1),
497  x_distance_diag);
498  raster_graph_edge_generator.conditionally_write_edge_slope(
499  bin_1_idx_and_z,
500  geo_raster.get_bin_idx_and_z_val_for_xy_bin(x_bin - 1, y_bin + 1),
501  x_distance_diag);
502  raster_graph_edge_generator.conditionally_write_edge_slope(
503  bin_1_idx_and_z,
504  geo_raster.get_bin_idx_and_z_val_for_xy_bin(x_bin + 1, y_bin + 1),
505  x_distance_diag);
506  raster_graph_edge_generator.conditionally_write_edge_slope(
507  bin_1_idx_and_z,
508  geo_raster.get_bin_idx_and_z_val_for_xy_bin(x_bin + 1, y_bin - 1),
509  x_distance_diag);
510  }
511  }
512  });
513  raster_graph_edge_generator.trim_to_size();
514  return raster_graph_edge_generator.raster_graph_edges;
515 }
516 
517 // clang-format off
518 /*
519  UDTF: tf_raster_graph_shortest_slope_weighted_path__cpu_template(TableFunctionManager,
520  Cursor<Column<T> x, Column<T> y, Column<Z> z> raster, TextEncodingNone agg_type,
521  T bin_dim | require="bin_dim > 0", bool geographic_coords,
522  int64_t neighborhood_fill_radius | require="neighborhood_fill_radius >= 0",
523  bool fill_only_nulls, T origin_x, T origin_y, T destination_x, T destination_y,
524  Z slope_weight_exponent, Z slope_pct_max) -> Column<int32_t> path_step, Column<T> x,
525  Column<T> y, T=[float, double], Z=[float, double]
526 */
527 // clang-format on
528 template <typename T, typename Z>
529 TEMPLATE_NOINLINE int32_t tf_raster_graph_shortest_slope_weighted_path__cpu_template(
531  const Column<T>& input_x,
532  const Column<T>& input_y,
533  const Column<Z>& input_z,
534  const TextEncodingNone& agg_type_str,
535  const T bin_dim_meters,
536  const bool geographic_coords,
537  const int64_t neighborhood_fill_radius,
538  const bool fill_only_nulls,
539  const T origin_x,
540  const T origin_y,
541  const T dest_x,
542  const T dest_y,
543  const Z slope_weight_exponent,
544  const Z slope_pct_max,
545  Column<int32_t>& output_path_step,
546  Column<T>& output_x,
547  Column<T>& output_y) {
548  const auto raster_agg_type = get_raster_agg_type(agg_type_str, false);
549  if (raster_agg_type == RasterAggType::INVALID) {
550  const std::string error_msg =
551  "Invalid Raster Aggregate Type: " + agg_type_str.getString();
552  return mgr.ERROR_MESSAGE(error_msg);
553  }
554 
555  GeoRaster<T, Z> geo_raster(input_x,
556  input_y,
557  input_z,
558  raster_agg_type,
559  bin_dim_meters,
560  geographic_coords,
561  true);
562 
563  // Check that origin and dest bin fall inside raster early
564  // to avoid needless work if we're out of bounds
565  const auto origin_x_bin = geo_raster.get_x_bin(origin_x);
566  const auto origin_y_bin = geo_raster.get_y_bin(origin_y);
567  const auto dest_x_bin = geo_raster.get_x_bin(dest_x);
568  const auto dest_y_bin = geo_raster.get_y_bin(dest_y);
569  if (geo_raster.is_bin_out_of_bounds(origin_x_bin, origin_y_bin)) {
570  return mgr.ERROR_MESSAGE("Origin coordinates are out of bounds.");
571  }
572  if (geo_raster.is_bin_out_of_bounds(dest_x_bin, dest_y_bin)) {
573  return mgr.ERROR_MESSAGE("Destination coordinates are out of bounds.");
574  }
575 
576  const auto origin_bin =
577  static_cast<int32_t>(geo_raster.get_bin_idx_for_xy_coords(origin_x, origin_y));
578  const auto dest_bin =
579  static_cast<int32_t>(geo_raster.get_bin_idx_for_xy_coords(dest_x, dest_y));
580  // Our CHECK macros don't parse the templated class properly, so fetch
581  // out of bounds value before calling CHECK_NE macro
582  const auto bin_out_of_bounds_idx = GeoRaster<T, Z>::BIN_OUT_OF_BOUNDS;
583  CHECK_NE(origin_bin, bin_out_of_bounds_idx);
584  CHECK_NE(dest_bin, bin_out_of_bounds_idx);
585 
586  if (neighborhood_fill_radius > 0) {
587  geo_raster.fill_bins_from_neighbors(
588  neighborhood_fill_radius, fill_only_nulls, RasterAggType::GAUSS_AVG);
589  }
590  auto raster_graph_edges =
591  generate_raster_graph_edges(geo_raster, slope_weight_exponent, slope_pct_max);
592  const Column<int32_t> node1(raster_graph_edges.bins_1);
593  const Column<int32_t> node2(raster_graph_edges.bins_2);
594  const Column<Z> distance(raster_graph_edges.distances);
595  GraphTraversalResults<int32_t, Z> graph_traversal_results;
596  TerminalNodes<int32_t> terminal_nodes(origin_bin, dest_bin);
597  try {
598  graph_traversal_results =
599  graph_shortest_path_impl(node1, node2, distance, terminal_nodes);
600  } catch (std::exception& e) {
601  return mgr.ERROR_MESSAGE(e.what());
602  }
603  auto output_timer = DEBUG_TIMER("Output shortest path results");
604  // const auto node1_val = graph_traversal_results.attr_idx_map.get_key_for_idx(
605  // graph_traversal_results.start_node_idx);
606  const auto node2_val = graph_traversal_results.attr_idx_map.get_key_for_idx(
607  graph_traversal_results.end_node_idx);
608  const int64_t vertex_idx = graph_traversal_results.end_node_idx;
609  int32_t current_node_idx = vertex_idx;
610  // Doing this first to determine how large the results are
611  int32_t num_hops = 0;
612  const auto distance_to_origin_node =
613  graph_traversal_results.graph_distances[vertex_idx];
614  if (distance_to_origin_node != std::numeric_limits<Z>::max()) {
615  while (current_node_idx != graph_traversal_results.start_node_idx) {
616  current_node_idx = graph_traversal_results.parents[current_node_idx];
617  num_hops++;
618  }
619  mgr.set_output_row_size(num_hops + 1);
620  current_node_idx = vertex_idx;
621  int32_t path_step_idx = num_hops;
622  const auto end_bin_idx =
623  graph_traversal_results.attr_idx_map.get_key_for_idx(current_node_idx);
624  CHECK_EQ(end_bin_idx, node2_val);
625  output_path_step[path_step_idx] = path_step_idx + 1;
626  const auto [end_path_x, end_path_y] =
627  geo_raster.get_xy_coords_for_bin_idx(end_bin_idx);
628  output_x[path_step_idx] = end_path_x;
629  output_y[path_step_idx] = end_path_y;
630  while (current_node_idx != graph_traversal_results.start_node_idx) {
631  current_node_idx = graph_traversal_results.parents[current_node_idx];
632  path_step_idx--;
633  const auto bin_idx =
634  graph_traversal_results.attr_idx_map.get_key_for_idx(current_node_idx);
635  output_path_step[path_step_idx] = path_step_idx + 1;
636  const auto [path_x, path_y] = geo_raster.get_xy_coords_for_bin_idx(bin_idx);
637  output_x[path_step_idx] = path_x;
638  output_y[path_step_idx] = path_y;
639  }
640  return num_hops + 1;
641  } else {
642  mgr.set_output_row_size(0);
643  return 0;
644  }
645 }
646 
647 #endif // #ifdef HAVE_TBB
648 #endif // #ifdef __CUDACC__
void set_output_row_size(int64_t num_rows)
#define CHECK_EQ(x, y)
Definition: Logger.h:301
int64_t get_bin_idx_for_xy_coords(const T x, const T y) const
std::pair< T, T > get_xy_coords_for_bin_idx(const int64_t bin_idx) const
std::string getString() const
Definition: heavydbTypes.h:311
DEVICE int64_t size() const
Definition: heavydbTypes.h:751
T * ptr_
Definition: heavydbTypes.h:725
void fill_bins_from_neighbors(const int64_t neighborhood_fill_radius, const bool fill_only_nulls, const RasterAggType raster_agg_type=RasterAggType::GAUSS_AVG)
T get_y_bin(const T input) const
T get_x_bin(const T input) const
bool is_bin_out_of_bounds(const int64_t x_bin, const int64_t y_bin) const
static constexpr int32_t INVALID_STR_ID
#define CHECK_NE(x, y)
Definition: Logger.h:302
int64_t num_rows_
Definition: heavydbTypes.h:726
bool g_enable_smem_group_by true
std::pair< int64_t, Z > get_bin_idx_and_z_val_for_xy_bin(const int64_t x_bin, const int64_t y_bin) const
DEVICE bool isNull(int64_t index) const
Definition: heavydbTypes.h:754
DEVICE void setNull(int64_t index)
Definition: heavydbTypes.h:755
void parallel_for(const blocked_range< Int > &range, const Body &body, const Partitioner &p=Partitioner())
constexpr unsigned N
Definition: Utm.h:110
bool g_enable_watchdog false
Definition: Execute.cpp:79
#define DEBUG_TIMER(name)
Definition: Logger.h:411
RasterAggType get_raster_agg_type(const std::string &agg_type_str, const bool is_fill_agg)
#define TEMPLATE_NOINLINE
Definition: heavydbTypes.h:54