[GRASS-SVN] r41777 - in grass/trunk: include include/Make
lib/vector lib/vector/neta
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Apr 10 11:12:32 EDT 2010
Author: martinl
Date: 2010-04-10 11:12:32 -0400 (Sat, 10 Apr 2010)
New Revision: 41777
Added:
grass/trunk/include/neta.h
grass/trunk/lib/vector/neta/
grass/trunk/lib/vector/neta/Makefile
grass/trunk/lib/vector/neta/allpairs.c
grass/trunk/lib/vector/neta/articulation_point.c
grass/trunk/lib/vector/neta/bridge.c
grass/trunk/lib/vector/neta/centrality.c
grass/trunk/lib/vector/neta/components.c
grass/trunk/lib/vector/neta/flow.c
grass/trunk/lib/vector/neta/path.c
grass/trunk/lib/vector/neta/spanningtree.c
grass/trunk/lib/vector/neta/timetables.c
grass/trunk/lib/vector/neta/utils.c
Modified:
grass/trunk/include/Make/Grass.make
grass/trunk/include/raster.h
Log:
added neta lib from grass-addons
Modified: grass/trunk/include/Make/Grass.make
===================================================================
--- grass/trunk/include/Make/Grass.make 2010-04-10 14:51:15 UTC (rev 41776)
+++ grass/trunk/include/Make/Grass.make 2010-04-10 15:12:32 UTC (rev 41777)
@@ -171,6 +171,7 @@
TRANS:trans \
VECTOR:vector \
VEDIT:vedit \
+ NETA:neta \
XDISPLAY:xdisplay \
XGD:xgd \
XGI:xgi \
@@ -220,6 +221,7 @@
TRANSDEPS = $(MATHLIB)
VECTORDEPS = $(DBMILIB) $(GRAPHLIB) $(DIG2LIB) $(LINKMLIB) $(RTREELIB) $(GISLIB) $(GEOSLIBS) $(GDALLIBS) $(MATHLIB)
VEDITDEPS = $(VECTORLIB) $(DBMILIB) $(GISLIB) $(MATHLIB)
+NETADEPS = $(VECTORDEP) $(DBMIDEP) $(GISDEP)
ifneq ($(USE_X11),)
CAIRODRIVERDEPS += $(XLIBPATH) $(XLIB) $(XEXTRALIBS)
Added: grass/trunk/include/neta.h
===================================================================
--- grass/trunk/include/neta.h (rev 0)
+++ grass/trunk/include/neta.h 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,114 @@
+
+/****************************************************************
+ *
+ * MODULE: netalib
+ *
+ * AUTHOR(S): Daniel Bundala (Google Summer of Code 2009)
+ *
+ * PURPOSE: NETwork Analysis Library
+ *
+ * COPYRIGHT: (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+ *
+ * This program is free software under the GNU General
+ * Public License (>=v2). Read the file COPYING that
+ * comes with GRASS for details.
+ *
+ ****************************************************************/
+
+#ifndef _NETA_H_
+#define _NETA_H_
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+
+/*bridge.c */
+int NetA_compute_bridges(dglGraph_s * graph, struct ilist *bridge_list);
+int NetA_articulation_points(dglGraph_s * graph,
+ struct ilist *articulation_list);
+
+/*components.c */
+int NetA_weakly_connected_components(dglGraph_s * graph, int *component);
+int NetA_strongly_connected_components(dglGraph_s * graph, int *component);
+
+/*spanningtree.c */
+int NetA_spanning_tree(dglGraph_s * graph, struct ilist *tree_list);
+
+/*allpairs.c */
+int NetA_allpairs(dglGraph_s * graph, dglInt32_t ** dist);
+
+/*neta_flow.c */
+int NetA_flow(dglGraph_s * graph, struct ilist *source_list,
+ struct ilist *sink_list, int *flow);
+int NetA_min_cut(dglGraph_s * graph, struct ilist *source_list,
+ struct ilist *sink_list, int *flow, struct ilist *cut);
+int NetA_split_vertices(dglGraph_s * in, dglGraph_s * out, int *node_costs);
+
+/*utils.c */
+void NetA_add_point_on_node(struct Map_info *In, struct Map_info *Out, int node,
+ struct line_cats *Cats);
+void NetA_points_to_nodes(struct Map_info *In, struct ilist *point_list);
+int NetA_get_node_costs(struct Map_info *In, int layer, char *column,
+ int *node_costs);
+void NetA_varray_to_nodes(struct Map_info *map, struct varray * varray,
+ struct ilist *nodes, int *nodes_to_features);
+int NetA_initialise_varray(struct Map_info *In, int layer, int mask_type,
+ char *where, char *cat, struct varray ** varray);
+/*centrality.c */
+void NetA_degree_centrality(dglGraph_s * graph, double *degree);
+int NetA_eigenvector_centrality(dglGraph_s * graph, int iterations,
+ double error, double *eigenvector);
+int NetA_betweenness_closeness(dglGraph_s * graph, double *betweenness,
+ double *closeness);
+
+/*path.c */
+int NetA_distance_from_points(dglGraph_s * graph, struct ilist *from, int *dst,
+ dglInt32_t ** prev);
+int NetA_find_path(dglGraph_s * graph, int from, int to, int *edges,
+ struct ilist *list);
+
+/*timetables.c */
+
+/*Structure containing all information about a timetable.
+ * Everything in indexed from 0.
+ */
+typedef struct
+{
+ int routes; /*Number of different routes. Two routes are different even if they differ only in time. */
+ int *route_length; /*Length of each route, i.e., number of stops */
+ int **route_stops; /*list of stops on each route in order (time increases) */
+ int **route_times; /*stop arrival times on overy route. Stops are given in the same order as above */
+ int stops; /*number of stops */
+ int *stop_length; /*Number of routes stopping at each stop */
+ int **stop_routes; /*List of routes for each stop. Routes are in increasing order */
+ int **stop_times; /*arrival times of routes for each stop. Routes are given in the same order as above */
+ int *walk_length; /*number of stops with "walking connection" for each stop */
+ int **walk_stops; /*list of stops within walking distance for each stop */
+ int **walk_times; /*walking times between stops as given above */
+} neta_timetable;
+
+typedef struct
+{
+ int **dst;
+ int **prev_stop;
+ int **prev_route;
+ int **prev_conn;
+ int rows, routes;
+} neta_timetable_result;
+int NetA_init_timetable_from_db(struct Map_info *In, int route_layer,
+ int walk_layer, char *route_id, char *times,
+ char *to_stop, char *walk_length,
+ neta_timetable * timetable, int **route_ids,
+ int **stop_ids);
+int NetA_timetable_shortest_path(neta_timetable * timetable, int from_stop,
+ int to_stop, int start_time, int min_change,
+ int max_changes, int walking_change,
+ neta_timetable_result * result);
+int NetA_timetable_get_route_time(neta_timetable * timetable, int stop,
+ int route);
+void NetA_timetable_result_release(neta_timetable_result * result);
+#endif
Property changes on: grass/trunk/include/neta.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Property changes on: grass/trunk/include/raster.h
___________________________________________________________________
Added: svn:mime-type
+ text/x-chdr
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Property changes on: grass/trunk/lib/vector/neta
___________________________________________________________________
Added: svn:ignore
+ OBJ.*
Added: grass/trunk/lib/vector/neta/Makefile
===================================================================
--- grass/trunk/lib/vector/neta/Makefile (rev 0)
+++ grass/trunk/lib/vector/neta/Makefile 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,13 @@
+MODULE_TOPDIR = ../../..
+
+LIB = NETA
+
+LIBES = $(VECTLIB) $(DBMILIB) $(GISLIB) $(GRAPHLIB)
+DEPENDENCIES= $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Lib.make
+include $(MODULE_TOPDIR)/include/Make/Doxygen.make
+
+default: lib
Property changes on: grass/trunk/lib/vector/neta/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-sh
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/allpairs.c
===================================================================
--- grass/trunk/lib/vector/neta/allpairs.c (rev 0)
+++ grass/trunk/lib/vector/neta/allpairs.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,93 @@
+/*!
+ \file vector/neta/allpairs.c
+
+ \brief Network Analysis library - shortest path between all pairs
+
+ Computes the length of the shortest path between all pairs of nodes
+ in the network.
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+
+/*!
+ \brief Stores the directed distance between every pair
+
+ \todo Use only O(W*W) memory where W is the number of nodes present in the graph
+
+ Upon the completion, dist stores the directed distance between every pair (i,j)
+ or -1 if the nodes are unreachable. It must be an array of dimension [nodes+1]*[nodes+1]
+
+ \param graph input graph
+ \param[out] dist list of directed distance
+
+ \return -1 on error
+ \return 0 on success
+ */
+int NetA_allpairs(dglGraph_s * graph, dglInt32_t ** dist)
+{
+ int nnodes, i, j, k, indices;
+ dglEdgesetTraverser_s et;
+ dglNodeTraverser_s nt;
+ dglInt32_t *node;
+ nnodes = dglGet_NodeCount(graph);
+ dglInt32_t *node_indices;
+ node_indices = (dglInt32_t *) G_calloc(nnodes, sizeof(dglInt32_t));
+ if (!node_indices) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ G_message(_("Computing all pairs shortest paths..."));
+ G_percent_reset();
+ for (i = 0; i <= nnodes; i++)
+ for (j = 0; j <= nnodes; j++)
+ dist[i][j] = -1;
+ dglNode_T_Initialize(&nt, graph);
+ indices = 0;
+ for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) {
+ dglInt32_t node_id = dglNodeGet_Id(graph, node);
+ node_indices[indices++] = node_id;
+ dglInt32_t *edge;
+ dglEdgeset_T_Initialize(&et, graph, dglNodeGet_OutEdgeset(graph, node));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et))
+ if (dglEdgeGet_Id(graph, edge) < 0) /*ignore backward edges */
+ dist[node_id][dglNodeGet_Id
+ (graph, dglEdgeGet_Tail(graph, edge))] =
+ dglEdgeGet_Cost(graph, edge);
+ dglEdgeset_T_Release(&et);
+ }
+ dglNode_T_Release(&nt);
+ for (k = 0; k < indices; k++) {
+ dglInt32_t k_index = node_indices[k];
+ G_percent(k + 1, indices, 1);
+ for (i = 0; i < indices; i++) {
+ dglInt32_t i_index = node_indices[i];
+ if (dist[i_index][k_index] == -1)
+ continue; /*no reason to proceed along infinite path */
+ for (j = 0; j < indices; j++) {
+ dglInt32_t j_index = node_indices[j];
+ if (dist[k_index][j_index] != -1 &&
+ (dist[i_index][k_index] + dist[k_index][j_index] <
+ dist[i_index][j_index] || dist[i_index][j_index] == -1)) {
+ dist[i_index][j_index] =
+ dist[i_index][k_index] + dist[k_index][j_index];
+ }
+ }
+ }
+ }
+
+ G_free(node_indices);
+ return 0;
+}
Property changes on: grass/trunk/lib/vector/neta/allpairs.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/articulation_point.c
===================================================================
--- grass/trunk/lib/vector/neta/articulation_point.c (rev 0)
+++ grass/trunk/lib/vector/neta/articulation_point.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,143 @@
+/*!
+ \file vector/neta/articulation_point.c
+
+ \brief Network Analysis library - connected components
+
+ Computes strongly and weakly connected components.
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+
+/*!
+ \brief Get number of articulation points in the graph
+
+ \param graph input graph
+ \param[out] articulation_list list of articulation points
+
+ \return number of points
+ \return -1 on error
+ */
+int NetA_articulation_points(dglGraph_s * graph,
+ struct ilist *articulation_list)
+{
+ int nnodes;
+ int points = 0;
+
+ dglEdgesetTraverser_s *current; /*edge to be processed when the node is visited */
+ int *tin, *min_tin; /*time in, and smallest tin over all successors. 0 if not yet visited */
+ dglInt32_t **parent; /*parents of the nodes */
+ dglInt32_t **stack; /*stack of nodes */
+ dglInt32_t **current_edge; /*current edge for each node */
+ int *mark; /*marked articulation points */
+ dglNodeTraverser_s nt;
+ dglInt32_t *current_node;
+ int stack_size;
+ int i, time;
+ nnodes = dglGet_NodeCount(graph);
+ current =
+ (dglEdgesetTraverser_s *) G_calloc(nnodes + 1,
+ sizeof(dglEdgesetTraverser_s));
+ tin = (int *)G_calloc(nnodes + 1, sizeof(int));
+ min_tin = (int *)G_calloc(nnodes + 1, sizeof(int));
+ parent = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
+ stack = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
+ current_edge = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
+ mark = (int *)G_calloc(nnodes + 1, sizeof(int));
+ if (!tin || !min_tin || !parent || !stack || !current || !mark) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+
+ for (i = 1; i <= nnodes; i++) {
+ dglEdgeset_T_Initialize(¤t[i], graph,
+ dglNodeGet_OutEdgeset(graph,
+ dglGetNode(graph, i)));
+ current_edge[i] = dglEdgeset_T_First(¤t[i]);
+ tin[i] = mark[i] = 0;
+ }
+
+ dglNode_T_Initialize(&nt, graph);
+
+ time = 0;
+ for (current_node = dglNode_T_First(&nt); current_node;
+ current_node = dglNode_T_Next(&nt)) {
+ dglInt32_t current_id = dglNodeGet_Id(graph, current_node);
+ if (tin[current_id] == 0) {
+ int children = 0; /*number of subtrees rooted at the root/current_node */
+ stack[0] = current_node;
+ stack_size = 1;
+ parent[current_id] = NULL;
+ while (stack_size) {
+ dglInt32_t *node = stack[stack_size - 1];
+ dglInt32_t node_id = dglNodeGet_Id(graph, node);
+ if (tin[node_id] == 0) /*vertex visited for the first time */
+ min_tin[node_id] = tin[node_id] = ++time;
+ else { /*return from the recursion */
+ dglInt32_t to = dglNodeGet_Id(graph,
+ dglEdgeGet_Tail(graph,
+ current_edge
+ [node_id]));
+ if (min_tin[to] >= tin[node_id]) /*no path from the subtree above the current node */
+ mark[node_id] = 1; /*so the current node must be an articulation point */
+
+ if (min_tin[to] < min_tin[node_id])
+ min_tin[node_id] = min_tin[to];
+ current_edge[node_id] = dglEdgeset_T_Next(¤t[node_id]); /*proceed to the next edge */
+ }
+ for (; current_edge[node_id]; current_edge[node_id] = dglEdgeset_T_Next(¤t[node_id])) { //try next edges
+ dglInt32_t *to =
+ dglEdgeGet_Tail(graph, current_edge[node_id]);
+ if (to == parent[node_id])
+ continue; /*skip parrent */
+ int to_id = dglNodeGet_Id(graph, to);
+ if (tin[to_id]) { /*back edge, cannot be a bridge/articualtion point */
+ if (tin[to_id] < min_tin[node_id])
+ min_tin[node_id] = tin[to_id];
+ }
+ else { /*forward edge */
+ if (node_id == current_id)
+ children++; /*if root, increase number of children */
+ parent[to_id] = node;
+ stack[stack_size++] = to;
+ break;
+ }
+ }
+ if (!current_edge[node_id])
+ stack_size--; /*current node completely processed */
+ }
+ if (children > 1)
+ mark[current_id] = 1; /*if the root has more than 1 subtrees rooted at it, then it is an
+ * articulation point */
+ }
+ }
+
+ for (i = 1; i <= nnodes; i++)
+ if (mark[i]) {
+ points++;
+ Vect_list_append(articulation_list, i);
+ }
+
+ dglNode_T_Release(&nt);
+ for (i = 1; i <= nnodes; i++)
+ dglEdgeset_T_Release(¤t[i]);
+
+ G_free(current);
+ G_free(tin);
+ G_free(min_tin);
+ G_free(parent);
+ G_free(stack);
+ G_free(current_edge);
+ return points;
+}
Property changes on: grass/trunk/lib/vector/neta/articulation_point.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/bridge.c
===================================================================
--- grass/trunk/lib/vector/neta/bridge.c (rev 0)
+++ grass/trunk/lib/vector/neta/bridge.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,133 @@
+/*!
+ \file vector/neta/bridges.c
+
+ \brief Network Analysis library - bridges
+
+ Computes number of bridges in the graph.
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+
+/*!
+ \brief Get number of bridges in the graph.
+
+ Bridge is an array containing the indices of the bridges.
+
+ \param graph input graph
+ \param[out] bridges_list list of bridges
+
+ \return number of bridges
+ \return -1 on error
+*/
+int NetA_compute_bridges(dglGraph_s * graph, struct ilist *bridge_list)
+{
+ int nnodes;
+ int bridges = 0;
+
+ dglEdgesetTraverser_s *current; /*edge to be processed when the node is visited */
+ int *tin, *min_tin; /*time in, and smallest tin over all successors. 0 if not yet visited */
+ dglInt32_t *parent; /*edge from parent to the node */
+ dglInt32_t **stack; /*stack of nodes */
+ dglInt32_t **current_edge; /*current edge for each node */
+ dglNodeTraverser_s nt;
+ dglInt32_t *current_node;
+ int stack_size;
+ int i, time;
+ nnodes = dglGet_NodeCount(graph);
+ current =
+ (dglEdgesetTraverser_s *) G_calloc(nnodes + 1,
+ sizeof(dglEdgesetTraverser_s));
+ tin = (int *)G_calloc(nnodes + 1, sizeof(int));
+ min_tin = (int *)G_calloc(nnodes + 1, sizeof(int));
+ parent = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
+ stack = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
+ current_edge = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
+ if (!tin || !min_tin || !parent || !stack || !current) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+
+ for (i = 1; i <= nnodes; i++) {
+ dglEdgeset_T_Initialize(¤t[i], graph,
+ dglNodeGet_OutEdgeset(graph,
+ dglGetNode(graph, i)));
+ current_edge[i] = dglEdgeset_T_First(¤t[i]);
+ tin[i] = 0;
+ }
+
+ dglNode_T_Initialize(&nt, graph);
+
+ time = 0;
+ for (current_node = dglNode_T_First(&nt); current_node;
+ current_node = dglNode_T_Next(&nt)) {
+ dglInt32_t current_id = dglNodeGet_Id(graph, current_node);
+ if (tin[current_id] == 0) {
+ stack[0] = current_node;
+ stack_size = 1;
+ parent[current_id] = 0;
+ while (stack_size) {
+ dglInt32_t *node = stack[stack_size - 1];
+ dglInt32_t node_id = dglNodeGet_Id(graph, node);
+ if (tin[node_id] == 0) /*vertex visited for the first time */
+ min_tin[node_id] = tin[node_id] = ++time;
+ else { /*return from the recursion */
+ dglInt32_t to = dglNodeGet_Id(graph,
+ dglEdgeGet_Tail(graph,
+ current_edge
+ [node_id]));
+ if (min_tin[to] > tin[node_id]) { /*no path from the subtree above the current node */
+ Vect_list_append(bridge_list, dglEdgeGet_Id(graph, current_edge[node_id])); /*so it must be a bridge */
+ bridges++;
+ }
+ if (min_tin[to] < min_tin[node_id])
+ min_tin[node_id] = min_tin[to];
+ current_edge[node_id] = dglEdgeset_T_Next(¤t[node_id]); /*proceed to the next edge */
+ }
+ for (; current_edge[node_id]; current_edge[node_id] = dglEdgeset_T_Next(¤t[node_id])) { //try next edges
+ dglInt32_t *to =
+ dglEdgeGet_Tail(graph, current_edge[node_id]);
+ dglInt32_t edge_id =
+ dglEdgeGet_Id(graph, current_edge[node_id]);
+ if (abs(edge_id) == parent[node_id])
+ continue; /*skip edge we used to travel to this node */
+ int to_id = dglNodeGet_Id(graph, to);
+ if (tin[to_id]) { /*back edge, cannot be a bridge/articualtion point */
+ if (tin[to_id] < min_tin[node_id])
+ min_tin[node_id] = tin[to_id];
+ }
+ else { /*forward edge */
+ parent[to_id] = abs(edge_id);
+ stack[stack_size++] = to;
+ break;
+ }
+ }
+ if (!current_edge[node_id])
+ stack_size--; /*current node completely processed */
+ }
+ }
+ }
+
+ dglNode_T_Release(&nt);
+ for (i = 1; i <= nnodes; i++)
+ dglEdgeset_T_Release(¤t[i]);
+
+ G_free(current);
+ G_free(tin);
+ G_free(min_tin);
+ G_free(parent);
+ G_free(stack);
+ G_free(current_edge);
+ return bridges;
+}
Property changes on: grass/trunk/lib/vector/neta/bridge.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/centrality.c
===================================================================
--- grass/trunk/lib/vector/neta/centrality.c (rev 0)
+++ grass/trunk/lib/vector/neta/centrality.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,237 @@
+/*!
+ \file vector/neta/centality.c
+
+ \brief Network Analysis library - centrality
+
+ Centrality measures
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+#include <grass/neta.h>
+
+/*!
+ \brief Computes degree centrality measure.
+
+ Array degree has to be properly initialised to nnodes+1 elements
+
+ \param graph input graph
+ \param[out] array of degrees
+*/
+void NetA_degree_centrality(dglGraph_s * graph, double *degree)
+{
+ int i;
+ double nnodes = dglGet_NodeCount(graph);
+ for (i = 1; i <= nnodes; i++)
+ degree[i] =
+ dglNodeGet_OutDegree(graph,
+ dglGetNode(graph, (dglInt32_t) i)) / nnodes;
+}
+
+/*!
+ \brief Computes eigenvector centrality using edge costs as weights.
+
+ \param graph input graph
+ \param iterations number of iterations
+ \param error ?
+ \param[out] eigenvector eigen vector value
+
+ \return 0 on success
+ \return -1 on failure
+*/
+int NetA_eigenvector_centrality(dglGraph_s * graph, int iterations,
+ double error, double *eigenvector)
+{
+ int i, iter, nnodes;
+ double *tmp;
+ nnodes = dglGet_NodeCount(graph);
+ tmp = (double *)G_calloc(nnodes + 1, sizeof(double));
+ if (!tmp) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+
+ error *= error;
+ for (i = 1; i <= nnodes; i++)
+ eigenvector[i] = 1;
+ for (iter = 0; iter < iterations; iter++) {
+ for (i = 1; i <= nnodes; i++)
+ tmp[i] = 0;
+ dglInt32_t *node;
+ dglNodeTraverser_s nt;
+ dglEdgesetTraverser_s et;
+ dglNode_T_Initialize(&nt, graph);
+ for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) {
+ dglInt32_t node_id = dglNodeGet_Id(graph, node);
+ double cur_value = eigenvector[node_id];
+ dglInt32_t *edge;
+ dglEdgeset_T_Initialize(&et, graph,
+ dglNodeGet_OutEdgeset(graph, node));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et))
+ tmp[dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge))] +=
+ cur_value * dglEdgeGet_Cost(graph, edge);
+
+ dglEdgeset_T_Release(&et);
+ }
+ dglNode_T_Release(&nt);
+ double cum_error = 0, max_value = tmp[1];
+ for (i = 2; i <= nnodes; i++)
+ if (tmp[i] > max_value)
+ max_value = tmp[i];
+ for (i = 1; i <= nnodes; i++) {
+ tmp[i] /= max_value;
+ cum_error += (tmp[i] - eigenvector[i]) * (tmp[i] - eigenvector[i]);
+ eigenvector[i] = tmp[i];
+ }
+ if (cum_error < error)
+ break;
+
+ }
+
+ G_free(tmp);
+ return 0;
+}
+
+/*!
+ \brief Computes betweenness and closeness centrality measure using Brandes algorithm.
+
+ Edge costs must be nonnegative. If some edge costs are negative then
+ the behaviour of this method is undefined.
+
+ \param graph input graph
+ \param[out] betweenness betweeness values
+ \param[out] closeness cloneness values
+
+ \return 0 on success
+ \return -1 on failure
+*/
+int NetA_betweenness_closeness(dglGraph_s * graph, double *betweenness,
+ double *closeness)
+{
+ int i, j, nnodes, stack_size, count;
+ dglInt32_t *dst, *node, *stack, *cnt, *delta;
+ dglNodeTraverser_s nt;
+ dglEdgesetTraverser_s et;
+ dglHeap_s heap;
+ struct ilist **prev;
+
+ nnodes = dglGet_NodeCount(graph);
+
+ dst = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
+ prev = (struct ilist **)G_calloc(nnodes + 1, sizeof(struct ilist *));
+ stack = (dglInt32_t *) G_calloc(nnodes, sizeof(dglInt32_t));
+ cnt = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
+ delta = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
+
+ if (!dst || !prev || !stack || !cnt || !delta) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+
+
+ for (i = 1; i <= nnodes; i++) {
+ prev[i] = Vect_new_list();
+ if (closeness)
+ closeness[i] = 0;
+ if (betweenness)
+ betweenness[i] = 0;
+ }
+
+ count = 0;
+ G_percent_reset();
+ dglNode_T_Initialize(&nt, graph);
+ for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) {
+ G_percent(count++, nnodes, 1);
+ dglInt32_t s = dglNodeGet_Id(graph, node);
+ dglHeapData_u heap_data;
+ dglHeapNode_s heap_node;
+
+ stack_size = 0;
+ for (i = 1; i <= nnodes; i++)
+ Vect_reset_list(prev[i]);
+ for (i = 1; i <= nnodes; i++) {
+ cnt[i] = 0;
+ dst[i] = -1;
+ }
+ dst[s] = 0;
+ cnt[s] = 1;
+ dglHeapInit(&heap);
+ heap_data.ul = s;
+ dglHeapInsertMin(&heap, 0, ' ', heap_data);
+ while (1) {
+ dglInt32_t v, dist;
+ if (!dglHeapExtractMin(&heap, &heap_node))
+ break;
+ v = heap_node.value.ul;
+ dist = heap_node.key;
+ if (dst[v] < dist)
+ continue;
+ stack[stack_size++] = v;
+
+ dglInt32_t *edge;
+ dglEdgeset_T_Initialize(&et, graph,
+ dglNodeGet_OutEdgeset(graph,
+ dglGetNode(graph,
+ v)));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t *to = dglEdgeGet_Tail(graph, edge);
+ dglInt32_t to_id = dglNodeGet_Id(graph, to);
+ dglInt32_t d = dglEdgeGet_Cost(graph, edge);
+ if (dst[to_id] == -1 || dst[to_id] > dist + d) {
+ dst[to_id] = dist + d;
+ Vect_reset_list(prev[to_id]);
+ heap_data.ul = to_id;
+ dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
+ }
+ if (dst[to_id] == dist + d) {
+ cnt[to_id] += cnt[v];
+ Vect_list_append(prev[to_id], v);
+ }
+ }
+
+ dglEdgeset_T_Release(&et);
+ }
+ dglHeapFree(&heap, NULL);
+ for (i = 1; i <= nnodes; i++)
+ delta[i] = 0;
+ for (i = stack_size - 1; i >= 0; i--) {
+ dglInt32_t w = stack[i];
+ if (closeness)
+ closeness[s] += dst[w];
+
+ for (j = 0; j < prev[w]->n_values; j++) {
+ dglInt32_t v = prev[w]->value[j];
+ delta[v] += (cnt[v] / (double)cnt[w]) * (1.0 + delta[w]);
+ }
+ if (w != s && betweenness)
+ betweenness[w] += delta[w];
+
+ }
+ if (closeness)
+ closeness[s] /= (double)stack_size;
+ }
+ dglNode_T_Release(&nt);
+
+ for (i = 1; i <= nnodes; i++)
+ Vect_destroy_list(prev[i]);
+ G_free(delta);
+ G_free(cnt);
+ G_free(stack);
+ G_free(prev);
+ G_free(dst);
+
+ return 0;
+};
Property changes on: grass/trunk/lib/vector/neta/centrality.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/components.c
===================================================================
--- grass/trunk/lib/vector/neta/components.c (rev 0)
+++ grass/trunk/lib/vector/neta/components.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,190 @@
+/*!
+ \file vector/neta/components.c
+
+ \brief Network Analysis library - graph componets
+
+ Computes strongly and weakly connected components.
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+
+/*!
+ \brief Computes weekly connected components
+
+ \param graph input graph
+ \param[out] component list of components
+
+ \return number of components
+ \return -1 on failure
+*/
+int NetA_weakly_connected_components(dglGraph_s * graph, int *component)
+{
+ int nnodes;
+ dglInt32_t *stack;
+ int *visited;
+ int stack_size, components;
+ dglInt32_t *cur_node;
+ dglNodeTraverser_s nt;
+
+ components = 0;
+ nnodes = dglGet_NodeCount(graph);
+ stack = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
+ visited = (int *)G_calloc(nnodes + 1, sizeof(int));
+ if (!stack || !visited) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+
+ dglNode_T_Initialize(&nt, graph);
+
+ for (cur_node = dglNode_T_First(&nt); cur_node;
+ cur_node = dglNode_T_Next(&nt)) {
+ dglInt32_t node_id = dglNodeGet_Id(graph, cur_node);
+ if (!visited[node_id]) {
+ visited[node_id] = 1;
+ stack[0] = node_id;
+ stack_size = 1;
+ component[node_id] = ++components;
+ while (stack_size) {
+ dglInt32_t *node, *edgeset, *edge;
+ dglEdgesetTraverser_s et;
+ node = dglGetNode(graph, stack[--stack_size]);
+ edgeset = dglNodeGet_OutEdgeset(graph, node);
+ dglEdgeset_T_Initialize(&et, graph, edgeset);
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t to;
+ to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
+ if (!visited[to]) {
+ visited[to] = 1;
+ component[to] = components;
+ stack[stack_size++] = to;
+ }
+ }
+ dglEdgeset_T_Release(&et);
+ }
+ }
+ }
+ dglNode_T_Release(&nt);
+ G_free(visited);
+ return components;
+}
+
+/*!
+ \brief Computes strongly connected components
+
+ \param graph input graph
+ \param[out] component list of components
+
+ \return number of components
+ \return -1 on failure
+*/
+int NetA_strongly_connected_components(dglGraph_s * graph, int *component)
+{
+ int nnodes;
+ dglInt32_t *stack, *order;
+ int *visited, *processed;
+ int stack_size, order_size, components;
+ dglInt32_t *node;
+ dglNodeTraverser_s nt;
+
+ components = 0;
+ nnodes = dglGet_NodeCount(graph);
+ stack = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
+ order = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
+ visited = (int *)G_calloc(nnodes + 1, sizeof(int));
+ processed = (int *)G_calloc(nnodes + 1, sizeof(int));
+ if (!stack || !visited || !order || !processed) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+
+ order_size = 0;
+ dglNode_T_Initialize(&nt, graph);
+
+ for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) {
+ dglInt32_t node_id = dglNodeGet_Id(graph, node);
+ component[node_id] = 0;
+ if (!visited[node_id]) {
+ visited[node_id] = 1;
+ stack[0] = node_id;
+ stack_size = 1;
+ while (stack_size) {
+ dglInt32_t *node, *edgeset, *edge;
+ dglEdgesetTraverser_s et;
+ dglInt32_t cur_node_id = stack[stack_size - 1];
+ if (processed[cur_node_id]) {
+ stack_size--;
+ order[order_size++] = cur_node_id;
+ continue;
+ }
+ processed[cur_node_id] = 1;
+ node = dglGetNode(graph, cur_node_id);
+ edgeset = dglNodeGet_OutEdgeset(graph, node);
+ dglEdgeset_T_Initialize(&et, graph, edgeset);
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t to;
+ if (dglEdgeGet_Id(graph, edge) < 0)
+ continue; /*ignore backward edges */
+ to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
+ if (!visited[to]) {
+ visited[to] = 1;
+ stack[stack_size++] = to;
+ }
+ }
+ dglEdgeset_T_Release(&et);
+ }
+ }
+ }
+
+ dglNode_T_Release(&nt);
+
+ while (order_size) {
+ dglInt32_t node_id = order[--order_size];
+ if (component[node_id])
+ continue;
+ components++;
+ component[node_id] = components;
+ stack[0] = node_id;
+ stack_size = 1;
+ while (stack_size) {
+ dglInt32_t *node, *edgeset, *edge;
+ dglEdgesetTraverser_s et;
+ dglInt32_t cur_node_id = stack[--stack_size];
+ node = dglGetNode(graph, cur_node_id);
+ edgeset = dglNodeGet_OutEdgeset(graph, node);
+ dglEdgeset_T_Initialize(&et, graph, edgeset);
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t to;
+ if (dglEdgeGet_Id(graph, edge) > 0)
+ continue; /*ignore forward edges */
+ to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
+ if (!component[to]) {
+ component[to] = components;
+ stack[stack_size++] = to;
+ }
+ }
+ dglEdgeset_T_Release(&et);
+ }
+ }
+
+ G_free(stack);
+ G_free(visited);
+ G_free(order);
+ G_free(processed);
+ return components;
+}
Property changes on: grass/trunk/lib/vector/neta/components.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/flow.c
===================================================================
--- grass/trunk/lib/vector/neta/flow.c (rev 0)
+++ grass/trunk/lib/vector/neta/flow.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,297 @@
+/*!
+ \file vector/neta/flow.c
+
+ \brief Network Analysis library - flow in graph
+
+ Computes the length of the shortest path between all pairs of nodes
+ in the network.
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+#include <grass/neta.h>
+
+dglInt32_t sign(dglInt32_t x)
+{
+ if (x >= 0)
+ return 1;
+ return -1;
+}
+
+/*!
+ \brief Get max flow from source to sink.
+
+ Array flow stores flow for each edge. Negative flow corresponds to a
+ flow in opposite direction The function assumes that the edge costs
+ correspond to edge capacities.
+
+ \param graph input graph
+ \param source_list list of sources
+ \param sink_list list of sinks
+ \param[out] flow max flows
+
+ \return number of flows
+ \return -1 on failure
+*/
+int NetA_flow(dglGraph_s * graph, struct ilist *source_list,
+ struct ilist *sink_list, int *flow)
+{
+ int nnodes, nlines, i;
+ dglEdgesetTraverser_s et;
+ dglInt32_t *queue;
+ dglInt32_t **prev;
+ char *is_source, *is_sink;
+ int begin, end, total_flow;
+
+ nnodes = dglGet_NodeCount(graph);
+ nlines = dglGet_EdgeCount(graph) / 2; /*each line corresponds to two edges. One in each direction */
+ queue = (dglInt32_t *) G_calloc(nnodes + 3, sizeof(dglInt32_t));
+ prev = (dglInt32_t **) G_calloc(nnodes + 3, sizeof(dglInt32_t *));
+ is_source = (char *)G_calloc(nnodes + 3, sizeof(char));
+ is_sink = (char *)G_calloc(nnodes + 3, sizeof(char));
+ if (!queue || !prev || !is_source || !is_sink) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+
+ for (i = 0; i < source_list->n_values; i++)
+ is_source[source_list->value[i]] = 1;
+ for (i = 0; i < sink_list->n_values; i++)
+ is_sink[sink_list->value[i]] = 1;
+
+ for (i = 0; i <= nlines; i++)
+ flow[i] = 0;
+
+ total_flow = 0;
+ while (1) {
+ dglInt32_t node, edge_id, min_residue;
+ int found = -1;
+ begin = end = 0;
+ for (i = 0; i < source_list->n_values; i++)
+ queue[end++] = source_list->value[i];
+
+ for (i = 1; i <= nnodes; i++) {
+ prev[i] = NULL;
+ }
+ while (begin != end && found == -1) {
+ dglInt32_t vertex = queue[begin++];
+ dglInt32_t *edge, *node = dglGetNode(graph, vertex);
+ dglEdgeset_T_Initialize(&et, graph,
+ dglNodeGet_OutEdgeset(graph, node));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
+ dglInt32_t id = dglEdgeGet_Id(graph, edge);
+ dglInt32_t to =
+ dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
+ if (!is_source[to] && prev[to] == NULL &&
+ cap > sign(id) * flow[abs(id)]) {
+ prev[to] = edge;
+ if (is_sink[to]) {
+ found = to;
+ break;
+ }
+ queue[end++] = to;
+ }
+ }
+ dglEdgeset_T_Release(&et);
+ }
+ if (found == -1)
+ break; /*no augmenting path */
+ /*find minimum residual capacity along the augmenting path */
+ node = found;
+ edge_id = dglEdgeGet_Id(graph, prev[node]);
+ min_residue =
+ dglEdgeGet_Cost(graph,
+ prev[node]) - sign(edge_id) * flow[abs(edge_id)];
+ while (!is_source[node]) {
+ dglInt32_t residue;
+ edge_id = dglEdgeGet_Id(graph, prev[node]);
+ residue =
+ dglEdgeGet_Cost(graph,
+ prev[node]) -
+ sign(edge_id) * flow[abs(edge_id)];
+ if (residue < min_residue)
+ min_residue = residue;
+ node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
+ }
+ total_flow += min_residue;
+ /*update flow along the augmenting path */
+ node = found;
+ while (!is_source[node]) {
+ edge_id = dglEdgeGet_Id(graph, prev[node]);
+ flow[abs(edge_id)] += sign(edge_id) * min_residue;
+ node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
+ }
+ }
+
+ G_free(queue);
+ G_free(prev);
+ G_free(is_source);
+ G_free(is_sink);
+
+ return total_flow;
+}
+
+/*!
+ \brief Calculates minimum cut between source(s) and sink(s).
+
+ Flow is the array produced by NetA_flow() method when called with
+ source_list and sink_list as the input. The output of this and
+ NetA_flow() method should be the same.
+
+ \param graph input graph
+ \param source_list list of sources
+ \param sink_list list of sinks
+ \param[out] cut list of edges (cut)
+
+ \return number of edges
+ \return -1 on failure
+*/
+int NetA_min_cut(dglGraph_s * graph, struct ilist *source_list,
+ struct ilist *sink_list, int *flow, struct ilist *cut)
+{
+ int nnodes, i;
+ dglEdgesetTraverser_s et;
+ dglInt32_t *queue;
+ char *visited;
+ int begin, end, total_flow;
+
+ nnodes = dglGet_NodeCount(graph);
+ queue = (dglInt32_t *) G_calloc(nnodes + 3, sizeof(dglInt32_t));
+ visited = (char *)G_calloc(nnodes + 3, sizeof(char));
+ if (!queue || !visited) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+
+ total_flow = begin = end = 0;
+
+ for (i = 1; i <= nnodes; i++)
+ visited[i] = 0;
+
+ for (i = 0; i < source_list->n_values; i++) {
+ queue[end++] = source_list->value[i];
+ visited[source_list->value[i]] = 1;
+ }
+
+ /* find vertices reachable from source(s) using only non-saturated edges */
+ while (begin != end) {
+ dglInt32_t vertex = queue[begin++];
+ dglInt32_t *edge, *node = dglGetNode(graph, vertex);
+ dglEdgeset_T_Initialize(&et, graph, dglNodeGet_OutEdgeset(graph, node));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
+ dglInt32_t id = dglEdgeGet_Id(graph, edge);
+ dglInt32_t to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
+ if (!visited[to] && cap > sign(id) * flow[abs(id)]) {
+ visited[to] = 1;
+ queue[end++] = to;
+ }
+ }
+ dglEdgeset_T_Release(&et);
+ }
+ /*saturated edges from reachable vertices to non-reachable ones form a minimum cost */
+ Vect_reset_list(cut);
+ for (i = 1; i <= nnodes; i++) {
+ if (!visited[i])
+ continue;
+ dglInt32_t *node, *edgeset, *edge;
+ node = dglGetNode(graph, i);
+ edgeset = dglNodeGet_OutEdgeset(graph, node);
+ dglEdgeset_T_Initialize(&et, graph, edgeset);
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t to, edge_id;
+ to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
+ edge_id = abs(dglEdgeGet_Id(graph, edge));
+ if (!visited[to] && flow[edge_id] != 0) {
+ Vect_list_append(cut, edge_id);
+ total_flow += abs(flow[abs(edge_id)]);
+ }
+ }
+ dglEdgeset_T_Release(&et);
+ }
+
+ G_free(visited);
+ G_free(queue);
+ return total_flow;
+}
+
+/*!
+ \brief Splits each vertex of in graph into two vertices
+
+ The method splits each vertex of in graph into two vertices: in
+ vertex and out vertex. Also, it adds an edge from an in vertex to
+ the corresponding out vertex (capacity=2) and it adds an edge from
+ out vertex to in vertex for each edge present in the in graph
+ (forward capacity=1, backward capacity=0). If the id of a vertex is
+ v then id of in vertex is 2*v-1 and of out vertex 2*v.
+
+ \param in from graph
+ \param out to graph
+ \param node_cost list of node costs
+
+ \return number of undirected edges in the graph
+ \return -1 on failure
+*/
+int NetA_split_vertices(dglGraph_s * in, dglGraph_s * out, int *node_costs)
+{
+ dglInt32_t opaqueset[16] =
+ { 360000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
+ dglNodeTraverser_s nt;
+ dglInt32_t nnodes, edge_cnt;
+ dglInt32_t *cur_node;
+ nnodes = dglGet_NodeCount(in);
+ dglInitialize(out, (dglByte_t) 1, (dglInt32_t) 0, (dglInt32_t) 0,
+ opaqueset);
+ dglNode_T_Initialize(&nt, in);
+ edge_cnt = 0;
+ dglInt32_t max_node_cost = 0;
+ for (cur_node = dglNode_T_First(&nt); cur_node;
+ cur_node = dglNode_T_Next(&nt)) {
+ dglInt32_t v = dglNodeGet_Id(in, cur_node);
+ edge_cnt++;
+ dglInt32_t cost = 1;
+ if (node_costs)
+ cost = node_costs[v];
+ if (cost > max_node_cost)
+ max_node_cost = cost;
+ dglAddEdge(out, 2 * v - 1, 2 * v, cost, edge_cnt);
+ dglAddEdge(out, 2 * v, 2 * v - 1, (dglInt32_t) 0, -edge_cnt);
+ }
+ dglNode_T_Release(&nt);
+ dglNode_T_Initialize(&nt, in);
+ for (cur_node = dglNode_T_First(&nt); cur_node;
+ cur_node = dglNode_T_Next(&nt)) {
+ dglEdgesetTraverser_s et;
+ dglInt32_t *edge;
+ dglInt32_t v = dglNodeGet_Id(in, cur_node);
+ dglEdgeset_T_Initialize(&et, in, dglNodeGet_OutEdgeset(in, cur_node));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t to;
+ to = dglNodeGet_Id(in, dglEdgeGet_Tail(in, edge));
+ edge_cnt++;
+ dglAddEdge(out, 2 * v, 2 * to - 1, max_node_cost + 1, edge_cnt);
+ dglAddEdge(out, 2 * to - 1, 2 * v, (dglInt32_t) 0, -edge_cnt);
+ }
+ dglEdgeset_T_Release(&et);
+ }
+ dglNode_T_Release(&nt);
+ if (dglFlatten(out) < 0)
+ G_fatal_error(_("GngFlatten error"));
+ return edge_cnt;
+}
Property changes on: grass/trunk/lib/vector/neta/flow.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/path.c
===================================================================
--- grass/trunk/lib/vector/neta/path.c (rev 0)
+++ grass/trunk/lib/vector/neta/path.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,171 @@
+/*!
+ \file vector/neta/path.c
+
+ \brief Network Analysis library - shortest path
+
+ Shortest paths from a set of nodes.
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+#include <grass/neta.h>
+
+/*!
+ \brief Computes shortests paths to every node from nodes in "from".
+
+ Array "dst" contains the length of the path or -1 if the node is not
+ reachable. Prev contains edges from predecessor along the shortest
+ path.
+
+ \param graph input graph
+ \param from list of 'from' positions
+ \param dst list of 'to' positions
+ \param[out] prev list of edges from predecessor along the shortest path
+
+ \return 0 on success
+ \return -1 on failure
+*/
+int NetA_distance_from_points(dglGraph_s * graph, struct ilist *from, int *dst,
+ dglInt32_t ** prev)
+{
+ int i, nnodes;
+ dglHeap_s heap;
+ nnodes = dglGet_NodeCount(graph);
+ dglEdgesetTraverser_s et;
+ for (i = 1; i <= nnodes; i++) {
+ dst[i] = -1;
+ prev[i] = NULL;
+ }
+
+ dglHeapInit(&heap);
+
+ for (i = 0; i < from->n_values; i++) {
+ int v = from->value[i];
+ if (dst[v] == 0)
+ continue; /*ingore duplicates */
+ dst[v] = 0;
+ dglHeapData_u heap_data;
+ heap_data.ul = v;
+ dglHeapInsertMin(&heap, 0, ' ', heap_data);
+ }
+ while (1) {
+ dglInt32_t v, dist;
+ dglHeapNode_s heap_node;
+ dglHeapData_u heap_data;
+ if (!dglHeapExtractMin(&heap, &heap_node))
+ break;
+ v = heap_node.value.ul;
+ dist = heap_node.key;
+ if (dst[v] < dist)
+ continue;
+
+ dglInt32_t *edge;
+ dglEdgeset_T_Initialize(&et, graph,
+ dglNodeGet_OutEdgeset(graph,
+ dglGetNode(graph, v)));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t *to = dglEdgeGet_Tail(graph, edge);
+ dglInt32_t to_id = dglNodeGet_Id(graph, to);
+ dglInt32_t d = dglEdgeGet_Cost(graph, edge);
+ if (dst[to_id] == -1 || dst[to_id] > dist + d) {
+ dst[to_id] = dist + d;
+ prev[to_id] = edge;
+ heap_data.ul = to_id;
+ dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
+ }
+ }
+
+ dglEdgeset_T_Release(&et);
+ }
+
+ dglHeapFree(&heap, NULL);
+
+ return 0;
+}
+
+/*!
+ \brief Find a path (minimum number of edges) from 'from' to 'to' using only edges in 'edges'.
+
+ Precisely, edge with id I is used iff edges[abs(i)] == 1. List
+ stores the indices of lines on the path. Method return number of
+ edges or -1 if no path exist.
+
+ \param graph input graph
+ \param from 'from' position
+ \param to 'to' position
+ \param edges list of available edges
+ \param[out] list list of edges
+
+ \return number of edges
+ \return -1 on failure
+*/
+int NetA_find_path(dglGraph_s * graph, int from, int to, int *edges,
+ struct ilist *list)
+{
+ dglInt32_t **prev, *queue;
+ dglEdgesetTraverser_s et;
+ char *vis;
+ int begin, end, cur, nnodes;
+
+ nnodes = dglGet_NodeCount(graph);
+ prev = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
+ queue = (dglInt32_t *) G_calloc(nnodes + 1, sizeof(dglInt32_t));
+ vis = (char *)G_calloc(nnodes + 1, sizeof(char));
+ if (!prev || !queue || !vis) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ Vect_reset_list(list);
+
+ begin = 0;
+ end = 1;
+ vis[from] = 'y';
+ queue[0] = from;
+ prev[from] = NULL;
+ while (begin != end) {
+ dglInt32_t vertex = queue[begin++];
+ if (vertex == to)
+ break;
+ dglInt32_t *edge, *node = dglGetNode(graph, vertex);
+ dglEdgeset_T_Initialize(&et, graph, dglNodeGet_OutEdgeset(graph, node));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t id = abs(dglEdgeGet_Id(graph, edge));
+ dglInt32_t to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
+ if (edges[id] && !vis[to]) {
+ vis[to] = 'y';
+ prev[to] = edge;
+ queue[end++] = to;
+ }
+ }
+ dglEdgeset_T_Release(&et);
+ }
+ G_free(queue);
+ if (!vis[to]) {
+ G_free(prev);
+ G_free(vis);
+ return -1;
+ }
+
+ cur = to;
+ while (prev[cur] != NULL) {
+ Vect_list_append(list, abs(dglEdgeGet_Id(graph, prev[cur])));
+ cur = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[cur]));
+ }
+
+ G_free(prev);
+ G_free(vis);
+ return list->n_values;
+}
Property changes on: grass/trunk/lib/vector/neta/path.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/spanningtree.c
===================================================================
--- grass/trunk/lib/vector/neta/spanningtree.c (rev 0)
+++ grass/trunk/lib/vector/neta/spanningtree.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,138 @@
+/*!
+ \file vector/neta/spanningtree.c
+
+ \brief Network Analysis library - spanning tree
+
+ Computes minimum spanning tree in the network.
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+
+struct union_find
+{
+ int *parent;
+};
+
+static int uf_initialize(struct union_find *uf, int size)
+{
+ int i;
+ uf->parent = (int *)G_calloc(size, sizeof(int));
+ if (!uf->parent)
+ return 0;
+ for (i = 0; i < size; i++)
+ uf->parent[i] = i;
+ return 1;
+}
+
+static void uf_release(struct union_find *uf)
+{
+ G_free(uf->parent);
+}
+
+static int uf_find(struct union_find *uf, int v)
+{
+ int cur = v, tmp;
+ while (uf->parent[cur] != cur)
+ cur = uf->parent[cur];
+ while (uf->parent[v] != v) {
+ tmp = uf->parent[v];
+ uf->parent[v] = cur;
+ v = tmp;
+ }
+ return cur;
+}
+
+/*TODO: union by rank */
+static void uf_union(struct union_find *uf, int u, int v)
+{
+ int parent_u = uf_find(uf, u);
+ int parent_v = uf_find(uf, v);
+ if (parent_u != parent_v)
+ uf->parent[parent_u] = parent_v;
+}
+
+typedef struct
+{
+ dglInt32_t cost;
+ dglInt32_t *edge;
+} edge_cost_pair;
+
+static int cmp_edge(const void *pa, const void *pb)
+{
+ return ((edge_cost_pair *) pa)->cost - ((edge_cost_pair *) pb)->cost;
+}
+
+/*!
+ \brief Get number of edges in the spanning forest
+
+ \param graph input graph
+ \param[out] list of edges
+
+ \return number of edges
+ \return -1 on failure
+*/
+int NetA_spanning_tree(dglGraph_s * graph, struct ilist *tree_list)
+{
+ int nnodes, edges, nedges, i, index;
+ edge_cost_pair *perm; /*permutaion of edges in ascending order */
+ struct union_find uf;
+ dglEdgesetTraverser_s et;
+ nnodes = dglGet_NodeCount(graph);
+ nedges = dglGet_EdgeCount(graph);
+ perm = (edge_cost_pair *) G_calloc(nedges, sizeof(edge_cost_pair));
+ if (!perm || !uf_initialize(&uf, nnodes + 1)) {
+ G_fatal_error(_("Out of memory"));
+ return -1;
+ }
+ /*for some obscure reasons, dglGetEdge always returns NULL. Therefore this complicated enumeration of the edges... */
+ index = 0;
+ G_message(_("Computing minimum spanning tree..."));
+ G_percent_reset();
+ for (i = 1; i <= nnodes; i++) {
+ G_percent(i, nnodes + nedges, 1);
+ dglInt32_t *edge;
+ dglEdgeset_T_Initialize(&et, graph,
+ dglNodeGet_OutEdgeset(graph,
+ dglGetNode(graph,
+ (dglInt32_t)
+ i)));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et))
+ if (dglEdgeGet_Id(graph, edge) > 0) {
+ perm[index].edge = edge;
+ perm[index].cost = dglEdgeGet_Cost(graph, edge);
+ index++;
+ }
+
+ dglEdgeset_T_Release(&et);
+ }
+ edges = 0;
+ qsort((void *)perm, index, sizeof(edge_cost_pair), cmp_edge);
+ for (i = 0; i < index; i++) {
+ G_percent(i + nnodes, nnodes + nedges, 1);
+ dglInt32_t head =
+ dglNodeGet_Id(graph, dglEdgeGet_Head(graph, perm[i].edge));
+ dglInt32_t tail =
+ dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, perm[i].edge));
+ if (uf_find(&uf, head) != uf_find(&uf, tail)) {
+ uf_union(&uf, head, tail);
+ edges++;
+ Vect_list_append(tree_list, dglEdgeGet_Id(graph, perm[i].edge));
+ }
+ }
+ G_free(perm);
+ uf_release(&uf);
+ return edges;
+}
Property changes on: grass/trunk/lib/vector/neta/spanningtree.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/timetables.c
===================================================================
--- grass/trunk/lib/vector/neta/timetables.c (rev 0)
+++ grass/trunk/lib/vector/neta/timetables.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,544 @@
+/*!
+ \file vector/neta/timetables.c
+
+ \brief Network Analysis library - timetables
+
+ Shortest path using timetables.
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/dbmi.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+#include <grass/neta.h>
+
+/*!
+ \brief Get number of distinct elements
+
+ \param driver DB driver
+ \param sql SQl string
+ \param[out] list of lengths
+ \param[out] list of ids
+
+ \return number of distinct elements
+ \return -1 on failure
+*/
+int NetA_init_distinct(dbDriver * driver, dbString * sql, int **lengths,
+ int **ids)
+{
+ int count, last, cur, result, index, more;
+ dbCursor cursor;
+ dbTable *table;
+ dbColumn *column;
+ dbValue *value;
+
+ if (db_open_select_cursor(driver, sql, &cursor, DB_SEQUENTIAL) != DB_OK) {
+ G_warning(_("Unable to open select cursor: %s"), db_get_string(sql));
+ return -1;
+ }
+ /*TODO: check column types */
+
+ count = last = 0;
+ /*count number of distinct routes */
+ table = db_get_cursor_table(&cursor);
+ column = db_get_table_column(table, 0);
+ while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
+ value = db_get_column_value(column);
+ cur = db_get_value_int(value);
+ if (count == 0 || cur != last) {
+ last = cur;
+ count++;
+ }
+ }
+ result = count;
+ db_close_cursor(&cursor);
+
+ *lengths = (int *)G_calloc(count, sizeof(int));
+ *ids = (int *)G_calloc(count, sizeof(int));
+ if (!*lengths || !*ids) {
+ G_warning(_("Out of memory"));
+ return -1;
+ }
+ db_open_select_cursor(driver, sql, &cursor, DB_SEQUENTIAL);
+ count = index = 0;
+ /*calculate the lengths of the routes */
+ table = db_get_cursor_table(&cursor);
+ column = db_get_table_column(table, 0);
+ while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
+ value = db_get_column_value(column);
+ cur = db_get_value_int(value);
+ if (count != 0 && cur != last)
+ index++;
+ if (count == 0 || cur != last)
+ (*ids)[index] = cur;
+ (*lengths)[index]++;
+ last = cur;
+ count++;
+ }
+ return result;
+}
+
+static int cmp_int(const void *a, const void *b)
+{
+ return *(int *)a - *(int *)b;
+}
+
+/*!
+ \brief Initialises timetable from a database
+
+ \param In pointer to Map_info structure
+ \param route_layer layer number of routes
+ \param walk_layer layer number of walkers
+ \param route_id id of route
+ \param times list of timestamps
+ \param to_stop ?
+ \param walk_length walk length as string
+ \param timetable pointer to neta_timetable
+ \param route_ids list of route ids
+ \param stop_ids lits of stop ids
+
+ \return 0 on success
+ \return non-zero value on failure
+*/
+int NetA_init_timetable_from_db(struct Map_info *In, int route_layer,
+ int walk_layer, char *route_id, char *times,
+ char *to_stop, char *walk_length,
+ neta_timetable * timetable, int **route_ids,
+ int **stop_ids)
+{
+ int more, i, stop, route, time, *stop_pnt, stop1, stop2;
+ dbString sql;
+ dbCursor cursor;
+ dbTable *table;
+ dbColumn *column1, *column2, *column3;
+ dbValue *value;
+ char buf[2000];
+
+ dbDriver *driver;
+ struct field_info *Fi;
+ Fi = Vect_get_field(In, route_layer);
+ driver = db_start_driver_open_database(Fi->driver, Fi->database);
+ if (driver == NULL)
+ G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+ Fi->database, Fi->driver);
+
+
+ db_init_string(&sql);
+ sprintf(buf, "select %s from %s order by %s", route_id, Fi->table,
+ route_id);
+ db_set_string(&sql, buf);
+ timetable->routes =
+ NetA_init_distinct(driver, &sql, &(timetable->route_length), route_ids);
+ if (timetable->routes < 0)
+ return 1;
+
+ sprintf(buf, "select %s from %s order by %s", Fi->key, Fi->table, Fi->key);
+ db_set_string(&sql, buf);
+ timetable->stops =
+ NetA_init_distinct(driver, &sql, &(timetable->stop_length), stop_ids);
+ if (timetable->stops < 0)
+ return 1;
+
+ timetable->route_stops = (int **)G_calloc(timetable->routes, sizeof(int *));
+ timetable->route_times = (int **)G_calloc(timetable->routes, sizeof(int *));
+ timetable->stop_routes = (int **)G_calloc(timetable->stops, sizeof(int *));
+ timetable->stop_times = (int **)G_calloc(timetable->stops, sizeof(int *));
+ timetable->walk_length = (int *)G_calloc(timetable->stops, sizeof(int));
+ timetable->walk_stops = (int **)G_calloc(timetable->stops, sizeof(int *));
+ timetable->walk_times = (int **)G_calloc(timetable->stops, sizeof(int *));
+ if (!timetable->route_stops || !timetable->route_times ||
+ !timetable->stop_routes || !timetable->stop_times ||
+ !timetable->walk_length) {
+ G_warning(_("Out of memory"));
+ return 2;
+ }
+
+ for (i = 0; i < timetable->routes; i++) {
+ timetable->route_stops[i] =
+ (int *)G_calloc(timetable->route_length[i], sizeof(int));
+ timetable->route_times[i] =
+ (int *)G_calloc(timetable->route_length[i], sizeof(int));
+ if (!timetable->route_stops[i] || !timetable->route_times[i]) {
+ G_warning(_("Out of memory"));
+ return 2;
+ }
+
+ timetable->route_length[i] = 0;
+ }
+
+ for (i = 0; i < timetable->stops; i++) {
+ timetable->stop_routes[i] =
+ (int *)G_calloc(timetable->stop_length[i], sizeof(int));
+ timetable->stop_times[i] =
+ (int *)G_calloc(timetable->stop_length[i], sizeof(int));
+ if (!timetable->stop_routes[i] || !timetable->stop_times[i]) {
+ G_warning(_("Out of memory"));
+ return 2;
+ }
+ timetable->walk_length[i] = 0;
+ timetable->stop_length[i] = 0;
+ }
+
+ sprintf(buf, "select %s, %s, %s from %s order by %s", Fi->key, route_id,
+ times, Fi->table, times);
+ db_set_string(&sql, buf);
+
+ if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) != DB_OK) {
+ G_warning(_("Unable to open select cursor: %s"), db_get_string(&sql));
+ return 1;
+ }
+
+
+ table = db_get_cursor_table(&cursor);
+ column1 = db_get_table_column(table, 0);
+ column2 = db_get_table_column(table, 1);
+ column3 = db_get_table_column(table, 2);
+ while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
+ value = db_get_column_value(column1);
+ stop = db_get_value_int(value);
+ value = db_get_column_value(column2);
+ route = db_get_value_int(value);
+ value = db_get_column_value(column3);
+ time = db_get_value_int(value);
+ stop =
+ (int *)bsearch(&stop, *stop_ids, timetable->stops, sizeof(int),
+ cmp_int) - (*stop_ids);
+ route =
+ (int *)bsearch(&route, *route_ids, timetable->routes, sizeof(int),
+ cmp_int) - (*route_ids);
+
+ timetable->stop_routes[stop][timetable->stop_length[stop]] = route;
+ timetable->stop_times[stop][timetable->stop_length[stop]++] = time;
+
+ timetable->route_stops[route][timetable->route_length[route]] = stop;
+ timetable->route_times[route][timetable->route_length[route]++] = time;
+ }
+ db_close_cursor(&cursor);
+
+ if (walk_layer != -1) {
+
+ Fi = Vect_get_field(In, walk_layer);
+ sprintf(buf, "select %s, %s, %s from %s", Fi->key, to_stop, walk_length,
+ Fi->table);
+ db_set_string(&sql, buf);
+
+ if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) !=
+ DB_OK) {
+ G_warning(_("Unable to open select cursor: %s"),
+ db_get_string(&sql));
+ return 1;
+ }
+
+ table = db_get_cursor_table(&cursor);
+ column1 = db_get_table_column(table, 0);
+ column2 = db_get_table_column(table, 1);
+ while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
+ value = db_get_column_value(column2);
+ stop = db_get_value_int(value);
+ stop_pnt =
+ (int *)bsearch(&stop, *stop_ids, timetable->stops, sizeof(int),
+ cmp_int);
+ if (stop_pnt) {
+ value = db_get_column_value(column1);
+ stop = db_get_value_int(value);
+ stop_pnt =
+ (int *)bsearch(&stop, *stop_ids, timetable->stops,
+ sizeof(int), cmp_int);
+ if (stop_pnt) {
+ stop = stop_pnt - (*stop_ids);
+ timetable->walk_length[stop]++;
+ }
+ }
+ }
+ db_close_cursor(&cursor);
+
+ for (i = 0; i < timetable->stops; i++) {
+ timetable->walk_stops[i] =
+ (int *)G_calloc(timetable->walk_length[i], sizeof(int));
+ timetable->walk_times[i] =
+ (int *)G_calloc(timetable->walk_length[i], sizeof(int));
+ if (!timetable->walk_stops[i] || !timetable->walk_times[i]) {
+ G_warning(_("Out of memory"));
+ return 2;
+ }
+ timetable->walk_length[i] = 0;
+ }
+
+ if (db_open_select_cursor(driver, &sql, &cursor, DB_SEQUENTIAL) !=
+ DB_OK) {
+ G_warning(_("Unable to open select cursor: %s"),
+ db_get_string(&sql));
+ return 1;
+ }
+
+ table = db_get_cursor_table(&cursor);
+ column1 = db_get_table_column(table, 0);
+ column2 = db_get_table_column(table, 1);
+ column3 = db_get_table_column(table, 2);
+ while (db_fetch(&cursor, DB_NEXT, &more) == DB_OK && more) {
+ value = db_get_column_value(column2);
+ stop = db_get_value_int(value);
+ stop_pnt =
+ (int *)bsearch(&stop, *stop_ids, timetable->stops, sizeof(int),
+ cmp_int);
+ if (stop_pnt) {
+ stop2 = stop_pnt - (*stop_ids);
+ value = db_get_column_value(column1);
+ stop = db_get_value_int(value);
+ stop_pnt =
+ (int *)bsearch(&stop, *stop_ids, timetable->stops,
+ sizeof(int), cmp_int);
+ if (stop_pnt) {
+ stop1 = stop_pnt - (*stop_ids);
+ value = db_get_column_value(column3);
+ time = db_get_value_int(value);
+ timetable->walk_stops[stop1][timetable->
+ walk_length[stop1]] = stop2;
+ timetable->walk_times[stop1][timetable->
+ walk_length[stop1]++] = time;
+ }
+ }
+ }
+ db_close_cursor(&cursor);
+ }
+ db_close_database_shutdown_driver(driver);
+
+ return 0;
+}
+
+typedef struct
+{
+ int v;
+ int conns;
+} neta_heap_data;
+
+static neta_heap_data *new_heap_data(int conns, int v)
+{
+ neta_heap_data *d = (neta_heap_data *) G_calloc(1, sizeof(neta_heap_data));
+ d->v = v;
+ d->conns = conns;
+ return d;
+}
+
+/*!
+ \brief Update Dijkstra structures
+
+ \param olc_conns old connection
+ \param new_conns new connection
+ \param to old 'to' node
+ \param new_dst new 'to' node
+ \param v ?
+ \param route id of route
+ \param rows ?
+ \param update ?
+ \param[out] result pointer to neta_timetable_result structure
+ \param heap ?
+*/
+void NetA_update_dijkstra(int old_conns, int new_conns, int to, int new_dst,
+ int v, int route, int rows, int update,
+ neta_timetable_result * result, dglHeap_s * heap)
+{
+ if (result->dst[new_conns][to] == -1 ||
+ result->dst[new_conns][to] > new_dst) {
+ result->dst[new_conns][to] = new_dst;
+ result->prev_stop[new_conns][to] = v;
+ result->prev_route[new_conns][to] = route;
+ result->prev_conn[new_conns][to] = old_conns;
+ if (update) {
+ dglHeapData_u heap_data;
+ heap_data.pv = (void *)new_heap_data(new_conns, to);
+ dglHeapInsertMin(heap, new_dst, ' ', heap_data);
+ }
+ }
+}
+
+/*!
+ \brief Computes the earliest arrival time.
+
+ Computes the earliest arrival time to to_stop from from_stop
+ starting at start_time, or -1 if no path exists
+
+ \param timetable pointer to neta_timetable structure
+ \param from_stop 'from' node
+ \param to_stop 'to' stop
+ \param start_time start timestamp
+ \param min_change ?
+ \param max_changes ?
+ \param walking_change ?
+ \param[out] result pointer to neta_timetable_result
+
+ \return ?
+ \return -1 on error
+*/
+int NetA_timetable_shortest_path(neta_timetable * timetable, int from_stop,
+ int to_stop, int start_time, int min_change,
+ int max_changes, int walking_change,
+ neta_timetable_result * result)
+{
+ int i, j;
+ dglHeap_s heap;
+
+ int opt_conns, rows = 1;
+ if (max_changes != -1)
+ rows = max_changes + 2;
+
+ result->rows = rows;
+ result->dst = (int **)G_calloc(rows, sizeof(int *));
+ result->prev_stop = (int **)G_calloc(rows, sizeof(int *));
+ result->prev_route = (int **)G_calloc(rows, sizeof(int *));
+ result->prev_conn = (int **)G_calloc(rows, sizeof(int *));
+
+ if (!result->dst || !result->prev_stop || !result->prev_route ||
+ !result->prev_conn) {
+ G_warning(_("Out of memory"));
+ return -1;
+ }
+
+ for (i = 0; i < rows; i++) {
+ result->dst[i] = (int *)G_calloc(timetable->stops, sizeof(int));
+ result->prev_stop[i] = (int *)G_calloc(timetable->stops, sizeof(int));
+ result->prev_route[i] = (int *)G_calloc(timetable->stops, sizeof(int));
+ result->prev_conn[i] = (int *)G_calloc(timetable->stops, sizeof(int));
+ if (!result->dst[i] || !result->prev_stop[i] || !result->prev_route[i]
+ || !result->prev_conn[i]) {
+ G_warning(_("Out of memory"));
+ return -1;
+ }
+ }
+
+ if (from_stop == to_stop) {
+ result->dst[0][to_stop] = start_time;
+ result->prev_route[0][to_stop] = result->prev_stop[0][to_stop] = -1;
+ result->routes = 0;
+ return start_time;
+ }
+
+ result->routes = -1;
+ if (walking_change > 1)
+ walking_change = 1;
+ if (walking_change < 0 || max_changes == -1)
+ walking_change = 0;
+ dglHeapInit(&heap);
+
+ for (i = 0; i < rows; i++)
+ for (j = 0; j < timetable->stops; j++)
+ result->dst[i][j] = result->prev_stop[i][j] =
+ result->prev_route[i][j] = -1;
+
+ result->dst[0][from_stop] = start_time - min_change;
+ result->prev_stop[0][from_stop] = result->prev_route[0][from_stop] = -1;
+ dglHeapData_u heap_data;
+ heap_data.pv = (void *)new_heap_data(0, from_stop);
+ dglHeapInsertMin(&heap, start_time - min_change, ' ', heap_data);
+
+ while (1) {
+ dglInt32_t v, dist, conns;
+ dglHeapNode_s heap_node;
+ int new_conns, walk_conns, update;
+ if (!dglHeapExtractMin(&heap, &heap_node))
+ break;
+ v = ((neta_heap_data *) (heap_node.value.pv))->v;
+ conns = ((neta_heap_data *) (heap_node.value.pv))->conns;
+ dist = heap_node.key;
+
+ if (dist > result->dst[conns][v])
+ continue;
+ if (v == to_stop)
+ break;
+ new_conns = (max_changes == -1) ? 0 : (conns + 1);
+ walk_conns = conns + walking_change;
+
+ /*walking */
+ if (walk_conns < rows) {
+ /* update = (max_changes == -1 || walk_conns <= max_changes); */
+ update = 1;
+ for (i = 0; i < timetable->walk_length[v]; i++) {
+ int to = timetable->walk_stops[v][i];
+ int new_dst = dist + timetable->walk_times[v][i];
+ NetA_update_dijkstra(conns, walk_conns, to, new_dst, v, -2,
+ rows, update, result, &heap);
+ }
+ }
+
+ if (new_conns >= rows)
+ continue;
+ /*process all routes arriving after dist+min_change */
+ for (i = 0; i < timetable->stop_length[v]; i++)
+ if (timetable->stop_times[v][i] >= dist + min_change) {
+ int route = timetable->stop_routes[v][i];
+ /*find the index of v on the route */
+ for (j = 0; j < timetable->route_length[route]; j++)
+ if (timetable->route_stops[route][j] == v)
+ break;
+ j++;
+ for (; j < timetable->route_length[route]; j++) {
+ int to = timetable->route_stops[route][j];
+ NetA_update_dijkstra(conns, new_conns, to,
+ timetable->route_times[route][j], v,
+ route, rows, 1, result, &heap);
+ }
+ }
+ }
+ dglHeapFree(&heap, NULL);
+ opt_conns = -1;
+ for (i = 0; i < rows; i++)
+ if (result->dst[i][to_stop] != -1 &&
+ (opt_conns == -1 ||
+ result->dst[opt_conns][to_stop] > result->dst[i][to_stop]))
+ opt_conns = i;
+ result->routes = opt_conns;
+
+ if (opt_conns != -1)
+ return result->dst[opt_conns][to_stop];
+ return -1;
+}
+
+/*!
+ \brief Get time
+
+ Get time when route "route" arrives at stop "stop" or -1.
+
+ \param timetable pointer to neta_timetable structure
+ \param stop 'stop' node id
+ \param route route id
+
+ \return time
+ \return -1 if not found
+*/
+int NetA_timetable_get_route_time(neta_timetable * timetable, int stop,
+ int route)
+{
+ int i;
+ for (i = 0; i < timetable->route_length[route]; i++)
+ if (timetable->route_stops[route][i] == stop)
+ return timetable->route_times[route][i];
+ return -1;
+}
+
+/*!
+ \brief Free neta_timetable_result structure
+
+ \param result pointer to neta_timetable_result structure
+*/
+void NetA_timetable_result_release(neta_timetable_result * result)
+{
+ int i;
+ for (i = 0; i < result->rows; i++) {
+ G_free(result->dst[i]);
+ G_free(result->prev_stop[i]);
+ G_free(result->prev_route[i]);
+ }
+ G_free(result->dst);
+ G_free(result->prev_stop);
+ G_free(result->prev_route);
+}
Property changes on: grass/trunk/lib/vector/neta/timetables.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass/trunk/lib/vector/neta/utils.c
===================================================================
--- grass/trunk/lib/vector/neta/utils.c (rev 0)
+++ grass/trunk/lib/vector/neta/utils.c 2010-04-10 15:12:32 UTC (rev 41777)
@@ -0,0 +1,238 @@
+/*!
+ \file vector/neta/timetables.c
+
+ \brief Network Analysis library - utils
+
+ Utils subroutines.
+
+ (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
+
+ This program is free software under the GNU General Public License
+ (>=v2). Read the file COPYING that comes with GRASS for details.
+
+ \author Daniel Bundala (Google Summer of Code 2009)
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/vector.h>
+#include <grass/glocale.h>
+#include <grass/dbmi.h>
+#include <grass/neta.h>
+
+
+/*!
+ \brief Writes point
+
+ Writes GV_POINT to Out at the position of the node in <em>In</em>.
+
+ \param In pointer to Map_info structure (input vector map)
+ \param[in,out] Out pointer to Map_info structure (output vector map)
+ \param node node id
+ \param Cats pointer to line_cats structures
+*/
+void NetA_add_point_on_node(struct Map_info *In, struct Map_info *Out, int node,
+ struct line_cats *Cats)
+{
+ static struct line_pnts *Points;
+ double x, y, z;
+ Points = Vect_new_line_struct();
+ Vect_get_node_coor(In, node, &x, &y, &z);
+ Vect_reset_line(Points);
+ Vect_append_point(Points, x, y, z);
+ Vect_write_line(Out, GV_POINT, Points, Cats);
+ Vect_destroy_line_struct(Points);
+}
+
+/* Returns the list of all points with the given category and field */
+/*void NetA_get_points_by_category(struct Map_info *In, int field, int cat, struct ilist *point_list)
+ * {
+ * int i, nlines;
+ * struct line_cats *Cats;
+ * Cats = Vect_new_cats_struct();
+ * Vect_get_num_lines(In);
+ * for(i=1;i<=nlines;i++){
+ * int type = Vect_read_line(In, NULL, Cats, i);
+ * if(type!=GV_POINT)continue;
+ * }
+ *
+ * Vect_destroy_cats_struct(Cats);
+ * }
+ */
+
+/*!
+ \brief Finds node
+
+ Find the node corresponding to each point in the point_list
+
+ \param In pointer to Map_info structure
+ \param point_list list of points (their ids)
+*/
+void NetA_points_to_nodes(struct Map_info *In, struct ilist *point_list)
+{
+ int i, node;
+ for (i = 0; i < point_list->n_values; i++) {
+ Vect_get_line_nodes(In, point_list->value[i], &node, NULL);
+ point_list->value[i] = node;
+ }
+}
+
+/*!
+ \brief Get node cost
+
+ For each node in the map, finds the category of the point on it (if
+ there is any) and stores the value associated with this category in
+ the array node_costs. If there is no point with a category,
+ node_costs=0.
+
+ node_costs are multiplied by 1000000 and truncated to integers (as
+ is done in Vect_net_build_graph)
+
+ \param In pointer to Map_info structure
+ \param layer layer number
+ \param column name of column
+ \param[out] node_costs list of node costs
+
+ \returns 1 on success
+ \return 0 on failure
+ */
+int NetA_get_node_costs(struct Map_info *In, int layer, char *column,
+ int *node_costs)
+{
+ int i, nlines, nnodes;
+ dbCatValArray vals;
+ struct line_cats *Cats;
+ struct line_pnts *Points;
+
+ dbDriver *driver;
+ struct field_info *Fi;
+ Fi = Vect_get_field(In, layer);
+ driver = db_start_driver_open_database(Fi->driver, Fi->database);
+ if (driver == NULL)
+ G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
+ Fi->database, Fi->driver);
+
+ nlines = Vect_get_num_lines(In);
+ nnodes = Vect_get_num_nodes(In);
+ Cats = Vect_new_cats_struct();
+ Points = Vect_new_line_struct();
+ for (i = 1; i <= nnodes; i++)
+ node_costs[i] = 0;
+
+ db_CatValArray_init(&vals);
+
+ if (db_select_CatValArray(driver, Fi->table, Fi->key, column, NULL, &vals)
+ == -1)
+ return 0;
+ for (i = 1; i <= nlines; i++) {
+ int type = Vect_read_line(In, Points, Cats, i);
+ if (type == GV_POINT) {
+ int node, cat;
+ double value;
+ if (!Vect_cat_get(Cats, layer, &cat))
+ continue;
+ Vect_get_line_nodes(In, i, &node, NULL);
+ if (db_CatValArray_get_value_double(&vals, cat, &value) == DB_OK)
+ node_costs[node] = value * 1000000.0;
+ }
+ }
+
+ Vect_destroy_cats_struct(Cats);
+ db_CatValArray_free(&vals);
+ db_close_database_shutdown_driver(driver);
+ return 1;
+}
+
+/*!
+ \brief Get list of nodes from varray
+
+ Returns the list of all nodes on features selected by varray.
+ nodes_to_features conains the index of a feature adjecent to each
+ node or -1 if no such feature specified by varray
+ exists. Nodes_to_features might be NULL, in which case it is left
+ unitialised.
+
+ \param map pointer to Map_info structure
+ \param varray pointer to varray structure
+ \param[out] nodes list of node ids
+ \param node_to_features ?
+*/
+void NetA_varray_to_nodes(struct Map_info *map, struct varray * varray,
+ struct ilist *nodes, int *nodes_to_features)
+{
+ int nlines, nnodes, i;
+ nlines = Vect_get_num_lines(map);
+ nnodes = Vect_get_num_nodes(map);
+ if (nodes_to_features)
+ for (i = 1; i <= nnodes; i++)
+ nodes_to_features[i] = -1;
+
+ for (i = 1; i <= nlines; i++)
+ if (varray->c[i]) {
+ int type = Vect_read_line(map, NULL, NULL, i);
+ if (type == GV_POINT) {
+ int node;
+ Vect_get_line_nodes(map, i, &node, NULL);
+ Vect_list_append(nodes, node);
+ if (nodes_to_features)
+ nodes_to_features[node] = i;
+ }
+ else {
+ int node1, node2;
+ Vect_get_line_nodes(map, i, &node1, &node2);
+ Vect_list_append(nodes, node1);
+ Vect_list_append(nodes, node2);
+ if (nodes_to_features)
+ nodes_to_features[node1] = nodes_to_features[node2] = i;
+ }
+ }
+}
+
+/*!
+ \brief Initialize varray
+
+ \param In pointer to Map_info structure
+ \param layer layer number
+ \param mask_type ?
+ \param where where statement
+ \param cat ?
+ \param[out] pointer to varray structure
+
+ \return ?
+*/
+int NetA_initialise_varray(struct Map_info *In, int layer, int mask_type,
+ char *where, char *cat, struct varray ** varray)
+{
+ /* parse filter option and select appropriate lines */
+ if (where) {
+ if (layer < 1)
+ G_fatal_error(_("'%s' must be > 0 for '%s'"), "layer", "where");
+ if (cat)
+ G_warning(_
+ ("'where' and 'cats' parameters were supplied, cat will be ignored"));
+ *varray = Vect_new_varray(Vect_get_num_lines(In));
+ if (Vect_set_varray_from_db
+ (In, layer, where, mask_type, 1, *varray) == -1) {
+ G_warning(_("Unable to load data from database"));
+ return 0;
+ }
+ return 1;
+ }
+ else if (cat) {
+ if (layer < 1)
+ G_fatal_error(_("'%s' must be > 0 for '%s'"), "layer", "cat");
+ *varray = Vect_new_varray(Vect_get_num_lines(In));
+ if (Vect_set_varray_from_cat_string
+ (In, layer, cat, mask_type, 1, *varray) == -1) {
+ G_warning(_("Problem loading category values"));
+ return 0;
+ }
+ return 1;
+ }
+ else {
+ return 2;
+ }
+
+
+}
Property changes on: grass/trunk/lib/vector/neta/utils.c
___________________________________________________________________
Added: svn:mime-type
+ text/x-csrc
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list