added neta lib from grass-addons

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@41777 15284696-431f-4ddb-bdfa-cd5b030d7da7

include/Make/Grass.make

@@ -171,6 +171,7 @@ libs = \
 	TRANS:trans \
 	VECTOR:vector \
 	VEDIT:vedit \
+	NETA:neta \
 	XDISPLAY:xdisplay \
 	XGD:xgd \
 	XGI:xgi \
@@ -220,6 +221,7 @@ SYMBDEPS         = $(GISLIB) $(MATHLIB)
 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)

include/neta.h

@@ -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

lib/vector/neta/Makefile

@@ -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

lib/vector/neta/allpairs.c

@@ -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;
+}

lib/vector/neta/articulation_point.c

@@ -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(&current[i], graph,
+				dglNodeGet_OutEdgeset(graph,
+						      dglGetNode(graph, i)));
+	current_edge[i] = dglEdgeset_T_First(&current[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(&current[node_id]);	/*proceed to the next edge */
+		}
+		for (; current_edge[node_id]; current_edge[node_id] = dglEdgeset_T_Next(&current[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(&current[i]);
+
+    G_free(current);
+    G_free(tin);
+    G_free(min_tin);
+    G_free(parent);
+    G_free(stack);
+    G_free(current_edge);
+    return points;
+}

lib/vector/neta/bridge.c

@@ -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(&current[i], graph,
+				dglNodeGet_OutEdgeset(graph,
+						      dglGetNode(graph, i)));
+	current_edge[i] = dglEdgeset_T_First(&current[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(&current[node_id]);	/*proceed to the next edge */
+		}
+		for (; current_edge[node_id]; current_edge[node_id] = dglEdgeset_T_Next(&current[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(&current[i]);
+
+    G_free(current);
+    G_free(tin);
+    G_free(min_tin);
+    G_free(parent);
+    G_free(stack);
+    G_free(current_edge);
+    return bridges;
+}

lib/vector/neta/centrality.c

@@ -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;
+};

lib/vector/neta/components.c

@@ -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;
+}

lib/vector/neta/flow.c

@@ -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;
+}

lib/vector/neta/path.c

@@ -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;
+}

lib/vector/neta/spanningtree.c

@@ -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;
+}

lib/vector/neta/timetables.c

@@ -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);
+}

lib/vector/neta/utils.c

@@ -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;
+    }
+
+
+}