diff --git a/configuration.conf b/configuration.conf index df075557a84..5f3c180d3dd 100644 --- a/configuration.conf +++ b/configuration.conf @@ -31,6 +31,7 @@ vrppdtw | N | Y | Y withPoints | Y | Y | Y lineGraph | Y | Y | Y components | Y | Y | Y +isochrones | Y | Y | N #---------------------- # SQL only directories #---------------------- diff --git a/include/drivers/isochrones/isochrones_driver.h b/include/drivers/isochrones/isochrones_driver.h new file mode 100644 index 00000000000..85f0f213f11 --- /dev/null +++ b/include/drivers/isochrones/isochrones_driver.h @@ -0,0 +1,53 @@ +/*PGR-GNU***************************************************************** +File: isochrones_driver.h + +Copyright (c) 2020 Vjeran Crnjak +vjeran@crnjak.xyz + +------ + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +********************************************************************PGR-GNU*/ + +#ifndef INCLUDE_DRIVERS_ISOCHRONES_MANY_TO_ISOCHRONES_H_ +#define INCLUDE_DRIVERS_ISOCHRONES_MANY_TO_ISOCHRONES_H_ + +/* for size_t */ +#ifdef __cplusplus +#include +#else +#include +#endif + +#include "c_types/general_path_element_t.h" +#include "c_types/pgr_edge_t.h" + +#ifdef __cplusplus +extern "C" { +#endif + +void do_pgr_many_to_isochrones(pgr_edge_t *edges, size_t total_edges, + int64_t *start_vertex, size_t s_len, + double distance, bool directed, bool equicost, + General_path_element_t **return_tuples, + size_t *return_count, char **log_msg, + char **notice_msg, char **err_msg); + +#ifdef __cplusplus +} +#endif + +#endif // INCLUDE_DRIVERS_ISOCHRONES_MANY_TO_ISOCHRONES_H_ diff --git a/sql/isochrones/CMakeLists.txt b/sql/isochrones/CMakeLists.txt new file mode 100644 index 00000000000..3aaf7ca2485 --- /dev/null +++ b/sql/isochrones/CMakeLists.txt @@ -0,0 +1,12 @@ + +SET(LOCAL_FILES + isochrones.sql +) + +foreach (f ${LOCAL_FILES}) +configure_file(${f} ${f}) +list(APPEND PACKAGE_SQL_FILES ${CMAKE_CURRENT_BINARY_DIR}/${f}) +endforeach() + +set(PgRouting_SQL_FILES ${PgRouting_SQL_FILES} ${PACKAGE_SQL_FILES} PARENT_SCOPE) + diff --git a/sql/isochrones/isochrones.sql b/sql/isochrones/isochrones.sql new file mode 100644 index 00000000000..df57fe2afc7 --- /dev/null +++ b/sql/isochrones/isochrones.sql @@ -0,0 +1,61 @@ +/*PGR-GNU***************************************************************** + +Copyright (c) 2020 Vjeran Crnjak +Mail: vjeran@crnjak.xyz + +------ + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +********************************************************************PGR-GNU*/ + +CREATE OR REPLACE FUNCTION pgr_isochrones( + edges_sql text, + start_vids anyarray, + distance FLOAT, + directed BOOLEAN DEFAULT TRUE, + equicost BOOLEAN DEFAULT FALSE, + OUT seq integer, + OUT from_v bigint, + OUT node bigint, + OUT edge bigint, + OUT cost FLOAT, + OUT agg_cost FLOAT) + RETURNS SETOF RECORD AS + '${MODULE_PATHNAME}', 'many_to_isochrones' + LANGUAGE c VOLATILE STRICT; + + +CREATE OR REPLACE FUNCTION pgr_isochrones( + edges_sql text, + start_vid bigint, + distance FLOAT8, + directed BOOLEAN DEFAULT TRUE, + OUT seq integer, + OUT node bigint, + OUT edge bigint, + OUT cost FLOAT, + OUT agg_cost FLOAT) + RETURNS SETOF RECORD AS +$BODY$ + SELECT a.seq, a.node, a.edge, a.cost, a.agg_cost + FROM pgr_isochrones($1, ARRAY[$2]::BIGINT[], $3, $4, false) a; +$BODY$ +LANGUAGE SQL VOLATILE STRICT +COST 100 +ROWS 1000; + + + diff --git a/sql/sigs/pgrouting--2.6.3.sig b/sql/sigs/pgrouting--2.6.3.sig index 069e6eb86ea..d4f5ad33159 100644 --- a/sql/sigs/pgrouting--2.6.3.sig +++ b/sql/sigs/pgrouting--2.6.3.sig @@ -75,6 +75,8 @@ pgr_dijkstravia(text,anyarray,boolean,boolean,boolean) pgr_drivingdistance(text,anyarray,double precision,boolean,boolean) pgr_drivingdistance(text,bigint,double precision,boolean) pgr_drivingdistance(text,bigint,double precision,boolean,boolean) +pgr_isochrones(text,anyarray,double precision,boolean,boolean) +pgr_isochrones(text,bigint,double precision,boolean) pgr_edgedisjointpaths(text,anyarray,anyarray,boolean) pgr_edgedisjointpaths(text,anyarray,bigint,boolean) pgr_edgedisjointpaths(text,bigint,anyarray,boolean) diff --git a/src/isochrones/CMakeLists.txt b/src/isochrones/CMakeLists.txt new file mode 100644 index 00000000000..84c34ed3824 --- /dev/null +++ b/src/isochrones/CMakeLists.txt @@ -0,0 +1,4 @@ +ADD_LIBRARY(isochrones OBJECT + many_to_isochrones.c + isochrones_driver.cpp + ) diff --git a/src/isochrones/isochrones_driver.cpp b/src/isochrones/isochrones_driver.cpp new file mode 100644 index 00000000000..944328c6a7a --- /dev/null +++ b/src/isochrones/isochrones_driver.cpp @@ -0,0 +1,271 @@ +/*PGR-GNU***************************************************************** +File: isochrones_driver.cpp + +Copyright (c) 2020 Vjeran Crnjak +vjeran@crnjak.xyz + +------ + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + ********************************************************************PGR-GNU*/ + +#include +#include +#include +#include +#include +#include + +#include "drivers/isochrones/isochrones_driver.h" + +#include "cpp_common/basic_vertex.h" +#include "cpp_common/pgr_alloc.hpp" +#include "cpp_common/pgr_assert.h" + +std::vector> +construct_adjacency_matrix(size_t n, const pgr_edge_t *edges, + size_t total_edges) { + std::vector> adj(n); + for (size_t i = 0; i < total_edges; ++i) { + if (edges[i].cost >= 0.) { + adj[edges[i].source].push_back(&edges[i]); + } + if (edges[i].reverse_cost >= 0.) { + adj[edges[i].target].push_back(&edges[i]); + } + } + return adj; +} + +struct Pred { + int64_t edge_id; + int64_t pred_id; + double cost; // cost of the edge by going from pred_id to target, cost can be + // partial if this is the leaf (driving_distance is respected). + // target can be either pgr_edge_t.{source,target}, depending on + // which cost is used pgr_edge_t.{cost,reverse_cost}. +}; + +void dijkstra(int64_t start_vertex, double driving_distance, + const std::vector> &adj, + std::vector *predecessors, std::vector *distances) { + size_t n = adj.size(); + predecessors->assign(n, {-1, -1, std::numeric_limits::infinity()}); + distances->assign(n, std::numeric_limits::infinity()); + (*predecessors)[0] = {-1, -1, 0.}; + + typedef std::tuple pq_el; // + std::set q; // priority queue + q.insert({0., start_vertex}); + while (!q.empty()) { + double dist; + int64_t node_id; + std::tie(dist, node_id) = *q.begin(); + if (dist >= driving_distance) { + break; + } + q.erase(q.begin()); + for (auto &&e : adj[node_id]) { + int64_t target = e->target == node_id ? e->source : e->target; + double cost = e->target == node_id ? e->reverse_cost : e->cost; + double agg_cost = dist + cost; + if ((*distances)[target] > agg_cost) { + q.erase({(*distances)[target], target}); + (*distances)[target] = + agg_cost > driving_distance ? driving_distance : agg_cost; + double edge_cost = + agg_cost > driving_distance ? (driving_distance - dist) : cost; + // Storing the partial travel edge cost. + (*predecessors)[target] = {e->id, node_id, edge_cost}; + q.emplace((*distances)[target], target); + } + } + } +} + +std::unordered_map remap_edges(pgr_edge_t *data_edges, + size_t total_edges) { + std::unordered_map mapping; + int64_t id = 0; + for (size_t i = 0; i < total_edges; ++i) { + pgr_edge_t *e = data_edges + i; + int64_t source_id, target_id; + auto it = mapping.find(e->source); + // better with if-initialization-statement + if (it != mapping.end()) { + source_id = it->second; + } else { + source_id = id++; + mapping[e->source] = source_id; + } + it = mapping.find(e->target); + if (it != mapping.end()) { + target_id = it->second; + } else { + target_id = id++; + mapping[e->target] = target_id; + } + e->source = source_id; + e->target = target_id; + } + return mapping; +} + +std::vector +do_many_dijkstras(pgr_edge_t *data_edges, size_t total_edges, + int64_t *start_vertex, size_t s_len, double distance) { + // Extracting vertices and mapping the ids from 0 to N-1. Remapping is done + // so that data structures used can be simpler (arrays instead of maps). + std::unordered_map mapping = + // modifying data_edges source/target fields. + remap_edges(data_edges, total_edges); + size_t nodes_count = mapping.size(); + std::vector results; + std::vector start_vertices(start_vertex, start_vertex + s_len); + auto adj = + construct_adjacency_matrix(mapping.size(), data_edges, total_edges); + // Storing the result of dijkstra call and reusing the memory for each vertex. + std::vector predecessors(nodes_count); + std::vector distances(nodes_count); + for (int64_t start_v : start_vertices) { + auto it = mapping.find(start_v); + // If start_v did not appear in edges then it has no particular mapping but + // pgr_drivingDistance result includes one row for this node. + if (it == mapping.end()) { + General_path_element_t r; + r.seq = 0; + r.start_id = start_v; + r.end_id = start_v; + r.node = start_v; + // -2 tags the unmapped starting vertex and won't use the reverse_mapping + // because mapping does not exist. -2 is changed to -1 later. + r.edge = -2; + r.cost = 0.0; + r.agg_cost = 0.0; + results.push_back(r); + continue; + } + // Calling the dijkstra algorithm and storing the results in predecessors + // and distances. + dijkstra(it->second, + /* driving_distance */ distance, adj, &predecessors, &distances); + // Appending the row results. + for (int64_t node_id = 0; node_id < predecessors.size(); ++node_id) { + const Pred &pred = predecessors[node_id]; + if (std::isinf(pred.cost) || pred.cost > distance) { + // Node not reached or reached too late is skipped. + continue; + } + General_path_element_t r; + // r.seq is filled later. r.node is remapped later. + r.start_id = start_v; + r.end_id = start_v; + r.edge = pred.edge_id; + r.node = node_id; + r.cost = pred.cost; + r.agg_cost = distances[node_id]; + results.push_back(r); + } + } + // The 0 to N-1 ids need to be mapped back to their previous values. + std::vector reverse_mapping(mapping.size()); + for (auto &m : mapping) { + reverse_mapping[m.second] = m.first; + } + int seq_cnt = 0; + int64_t prev_start_id = -1; + for (size_t i = 0; i < results.size(); ++i) { + auto &r = results[i]; + if (r.edge == -2) { + // This is a starting vertex that has no edges in graph. + r.seq = 0; + // Changing back the -2 edge indicator to -1. + r.edge = -1; + // Force seq_cnt reset in next iteration and avoid accidental id clash + // with the mapping. + prev_start_id = -1; + continue; + } + if (r.start_id != prev_start_id) { + prev_start_id = r.start_id; + seq_cnt = 0; + } + r.seq = seq_cnt++; + r.node = reverse_mapping[r.node]; + } + return results; +} + +void do_pgr_many_to_isochrones(pgr_edge_t *data_edges, size_t total_edges, + int64_t *start_vertex, size_t s_len, + double distance, bool directedFlag, + bool equiCostFlag, + General_path_element_t **return_tuples, + size_t *return_count, char **log_msg, + char **notice_msg, char **err_msg) { + std::ostringstream log; + std::ostringstream err; + std::ostringstream notice; + + try { + pgassert(total_edges != 0); + pgassert(!(*log_msg)); + pgassert(!(*notice_msg)); + pgassert(!(*err_msg)); + pgassert(!(*return_tuples)); + pgassert(*return_count == 0); + pgassert((*return_tuples) == NULL); + + auto results = do_many_dijkstras(data_edges, total_edges, start_vertex, + s_len, distance); + + size_t count(results.size()); + if (count == 0) { + log << "\nNo return values were found"; + *notice_msg = pgr_msg(log.str().c_str()); + return; + } + *return_tuples = pgr_alloc(count, (*return_tuples)); + *return_count = count; + + // Moving results to allocated return_tuples array. + for (size_t i = 0; i < count; ++i) { + (*return_tuples)[i] = results[i]; + } + + *log_msg = log.str().empty() ? *log_msg : pgr_msg(log.str().c_str()); + *notice_msg = + notice.str().empty() ? *notice_msg : pgr_msg(notice.str().c_str()); + } catch (AssertFailedException &except) { + (*return_tuples) = pgr_free(*return_tuples); + (*return_count) = 0; + err << except.what(); + *err_msg = pgr_msg(err.str().c_str()); + *log_msg = pgr_msg(log.str().c_str()); + } catch (std::exception &except) { + (*return_tuples) = pgr_free(*return_tuples); + (*return_count) = 0; + err << except.what(); + *err_msg = pgr_msg(err.str().c_str()); + *log_msg = pgr_msg(log.str().c_str()); + } catch (...) { + (*return_tuples) = pgr_free(*return_tuples); + (*return_count) = 0; + err << "Caught unknown exception!"; + *err_msg = pgr_msg(err.str().c_str()); + *log_msg = pgr_msg(log.str().c_str()); + } +} diff --git a/src/isochrones/many_to_isochrones.c b/src/isochrones/many_to_isochrones.c new file mode 100644 index 00000000000..aca7c3c8264 --- /dev/null +++ b/src/isochrones/many_to_isochrones.c @@ -0,0 +1,187 @@ +/*PGR-GNU***************************************************************** +File: many_to_isochrones.c + +Copyright (c) 2020 Vjeran Crnjak +Mail: vjeran@crnjak.xyz + +------ + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License +along with this program; if not, write to the Free Software +Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + ********************************************************************PGR-GNU*/ + +#include "c_common/postgres_connection.h" +#include "utils/array.h" + +#include "c_common/debug_macro.h" +#include "c_common/e_report.h" +#include "c_common/time_msg.h" + +#include "c_common/edges_input.h" +#include "c_common/arrays_input.h" +#include "drivers/isochrones/isochrones_driver.h" + + +PGDLLEXPORT Datum many_to_isochrones(PG_FUNCTION_ARGS); +PG_FUNCTION_INFO_V1(many_to_isochrones); + + +static +void process( + char* sql, + ArrayType *starts, + float8 distance, + bool directed, + bool equicost, + General_path_element_t **result_tuples, + size_t *result_count) { + pgr_SPI_connect(); + + size_t size_start_vidsArr = 0; + int64_t* start_vidsArr = pgr_get_bigIntArray(&size_start_vidsArr, starts); + + pgr_edge_t *edges = NULL; + size_t total_tuples = 0; + pgr_get_edges(sql, &edges, &total_tuples); + + if (total_tuples == 0) { + return; + } + + PGR_DBG("Starting timer"); + clock_t start_t = clock(); + char* log_msg = NULL; + char* notice_msg = NULL; + char* err_msg = NULL; + do_pgr_many_to_isochrones( + edges, total_tuples, + start_vidsArr, size_start_vidsArr, + distance, + directed, + equicost, + result_tuples, result_count, + &log_msg, + ¬ice_msg, + &err_msg); + + time_msg("processing pgr_isochrones()", + start_t, clock()); + + if (err_msg && (*result_tuples)) { + pfree(*result_tuples); + (*result_tuples) = NULL; + (*result_count) = 0; + } + + pgr_global_report(log_msg, notice_msg, err_msg); + + if (log_msg) pfree(log_msg); + if (notice_msg) pfree(notice_msg); + if (err_msg) pfree(err_msg); + if (edges) pfree(edges); + if (start_vidsArr) pfree(start_vidsArr); + + pgr_SPI_finish(); +} + + +PGDLLEXPORT Datum +many_to_isochrones(PG_FUNCTION_ARGS) { + FuncCallContext *funcctx; + TupleDesc tuple_desc; + + /**********************************************************************/ + General_path_element_t* result_tuples = 0; + size_t result_count = 0; + /**********************************************************************/ + + if (SRF_IS_FIRSTCALL()) { + MemoryContext oldcontext; + + funcctx = SRF_FIRSTCALL_INIT(); + oldcontext = MemoryContextSwitchTo(funcctx->multi_call_memory_ctx); + + /**********************************************************************/ + + PGR_DBG("Calling driving_many_to_dist_driver"); + process( + text_to_cstring(PG_GETARG_TEXT_P(0)), + PG_GETARG_ARRAYTYPE_P(1), + PG_GETARG_FLOAT8(2), + PG_GETARG_BOOL(3), + PG_GETARG_BOOL(4), + &result_tuples, &result_count); + + /**********************************************************************/ + + +#if PGSQL_VERSION > 95 + funcctx->max_calls = result_count; +#else + funcctx->max_calls = (uint32_t)result_count; +#endif + funcctx->user_fctx = result_tuples; + if (get_call_result_type(fcinfo, NULL, &tuple_desc) + != TYPEFUNC_COMPOSITE) + ereport(ERROR, + (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), + errmsg("function returning record called in context " + "that cannot accept type record"))); + + funcctx->tuple_desc = tuple_desc; + + MemoryContextSwitchTo(oldcontext); + } + + funcctx = SRF_PERCALL_SETUP(); + + tuple_desc = funcctx->tuple_desc; + result_tuples = (General_path_element_t*) funcctx->user_fctx; + + if (funcctx->call_cntr < funcctx->max_calls) { + HeapTuple tuple; + Datum result; + Datum *values; + bool* nulls; + + /**********************************************************************/ + size_t numb = 6; + values = palloc(numb * sizeof(Datum)); + nulls = palloc(numb * sizeof(bool)); + + size_t i; + for (i = 0; i < numb; ++i) { + nulls[i] = false; + } + values[0] = Int32GetDatum(funcctx->call_cntr + 1); + values[1] = Int64GetDatum(result_tuples[funcctx->call_cntr].start_id); + values[2] = Int64GetDatum(result_tuples[funcctx->call_cntr].node); + values[3] = Int64GetDatum(result_tuples[funcctx->call_cntr].edge); + values[4] = Float8GetDatum(result_tuples[funcctx->call_cntr].cost); + values[5] = Float8GetDatum(result_tuples[funcctx->call_cntr].agg_cost); + + /**********************************************************************/ + tuple = heap_form_tuple(tuple_desc, values, nulls); + result = HeapTupleGetDatum(tuple); + + pfree(values); + pfree(nulls); + + SRF_RETURN_NEXT(funcctx, result); + } else { + SRF_RETURN_DONE(funcctx); + } +} +