netalib: backport bugfixes from trunk

git-svn-id: https://svn.osgeo.org/grass/grass/branches/releasebranch_7_0@68024 15284696-431f-4ddb-bdfa-cd5b030d7da7

include/defs/neta.h

@@ -40,6 +40,8 @@ int NetA_betweenness_closeness(dglGraph_s * graph, double *betweenness,
 /*path.c */
 int NetA_distance_from_points(dglGraph_s * graph, struct ilist *from, int *dst,
 			      dglInt32_t ** prev);
+int NetA_distance_to_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);
 

lib/vector/neta/components.c

@@ -11,6 +11,26 @@
    (>=v2). Read the file COPYING that comes with GRASS for details.
 
    \author Daniel Bundala (Google Summer of Code 2009)
+   \author Markus Metz
+ */
+
+/* example:
+ * 
+ * X -->-- X ---- X --<-- X ---- X
+ * N1      N2     N3      N4     N5
+ * 
+ * -->--, --<-- one-way
+ * ---- both ways
+ * 
+ * weakly connected:
+ * all 5 nodes, even though there is no direct path from N1 to N4, 5
+ * but N1 connects to N2, 3, and N4, 5 also connect to N2, 3
+ * 
+ * strongly connected:
+ * no path from N2 to N1, no path from N3 to N4
+ * component 1: N1
+ * component 2: N2, 3
+ * Component3: N4, 5
  */
 
 #include <stdio.h>
@@ -21,43 +41,53 @@
 #include <grass/dgl/graph.h>
 
 /*!
-   \brief Computes weekly connected components
+   \brief Computes weakly connected components
 
    \param graph input graph
-   \param[out] component list of components
+   \param[out] component array of component ids
 
    \return number of components
    \return -1 on failure
  */
 int NetA_weakly_connected_components(dglGraph_s * graph, int *component)
 {
-    int nnodes;
+    int nnodes, i;
     dglInt32_t *stack;
-    int *visited;
     int stack_size, components;
     dglInt32_t *cur_node;
     dglNodeTraverser_s nt;
+    int have_node_costs;
+    dglInt32_t ncost;
+
+    if (graph->Version < 2) {
+	G_warning("Directed graph must be version 2 or 3 for NetA_weakly_connected_components()");
+	return -1;
+    }
 
     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) {
+    if (!stack) {
 	G_fatal_error(_("Out of memory"));
 	return -1;
     }
 
+    for (i = 1; i <= nnodes; i++)
+	component[i] = 0;
+
+    ncost = 0;
+    have_node_costs = dglGet_NodeAttrSize(graph);
+
     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);
+	dglInt32_t cur_node_id = dglNodeGet_Id(graph, cur_node);
 
-	if (!visited[node_id]) {
-	    visited[node_id] = 1;
-	    stack[0] = node_id;
+	if (!component[cur_node_id]) {
+	    stack[0] = cur_node_id;
 	    stack_size = 1;
-	    component[node_id] = ++components;
+	    component[cur_node_id] = ++components;
 	    while (stack_size) {
 		dglInt32_t *node, *edgeset, *edge;
 		dglEdgesetTraverser_s et;
@@ -70,10 +100,35 @@ int NetA_weakly_connected_components(dglGraph_s * graph, int *component)
 		    dglInt32_t to;
 
 		    to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
-		    if (!visited[to]) {
-			visited[to] = 1;
+		    if (!component[to]) {
 			component[to] = components;
-			stack[stack_size++] = to;
+			/* do not go through closed nodes */
+			if (have_node_costs) {
+			    memcpy(&ncost, dglNodeGet_Attr(graph, dglEdgeGet_Tail(graph, edge)),
+				   sizeof(ncost));
+			}
+			if (ncost >= 0)
+			    stack[stack_size++] = to;
+		    }
+		}
+		dglEdgeset_T_Release(&et);
+
+		edgeset = dglNodeGet_InEdgeset(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_Head(graph, edge));
+		    if (!component[to]) {
+			component[to] = components;
+			/* do not go through closed nodes */
+			if (have_node_costs) {
+			    memcpy(&ncost, dglNodeGet_Attr(graph, dglEdgeGet_Tail(graph, edge)),
+				   sizeof(ncost));
+			}
+			if (ncost >= 0)
+			    stack[stack_size++] = to;
 		    }
 		}
 		dglEdgeset_T_Release(&et);
@@ -81,73 +136,95 @@ int NetA_weakly_connected_components(dglGraph_s * graph, int *component)
 	}
     }
     dglNode_T_Release(&nt);
-    G_free(visited);
+
+    G_free(stack);
     return components;
 }
 
 /*!
-   \brief Computes strongly connected components
+   \brief Computes strongly connected components with Kosaraju's 
+   two-pass algorithm
 
    \param graph input graph
-   \param[out] component list of components
+   \param[out] component array of component ids
 
    \return number of components
    \return -1 on failure
  */
 int NetA_strongly_connected_components(dglGraph_s * graph, int *component)
 {
-    int nnodes;
+    int nnodes, i;
     dglInt32_t *stack, *order;
-    int *visited, *processed;
+    int *processed;
     int stack_size, order_size, components;
-    dglInt32_t *node;
+    dglInt32_t *cur_node;
     dglNodeTraverser_s nt;
+    int have_node_costs;
+    dglInt32_t ncost;
+
+    if (graph->Version < 2) {
+	G_warning("Directed graph must be version 2 or 3 for NetA_strongly_connected_components()");
+	return -1;
+    }
 
     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) {
+    if (!stack || !order || !processed) {
 	G_fatal_error(_("Out of memory"));
 	return -1;
     }
 
+    for (i = 1; i <= nnodes; i++) {
+	component[i] = 0;
+    }
+
+    ncost = 0;
+    have_node_costs = dglGet_NodeAttrSize(graph);
+
     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);
+    for (cur_node = dglNode_T_First(&nt); cur_node;
+	 cur_node = dglNode_T_Next(&nt)) {
+	dglInt32_t cur_node_id = dglNodeGet_Id(graph, cur_node);
 
-	component[node_id] = 0;
-	if (!visited[node_id]) {
-	    visited[node_id] = 1;
-	    stack[0] = node_id;
+	if (!component[cur_node_id]) {
+	    component[cur_node_id] = --components;
+	    stack[0] = cur_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];
+		dglInt32_t node_id = stack[stack_size - 1];
 
-		if (processed[cur_node_id]) {
+		if (processed[node_id]) {
 		    stack_size--;
-		    order[order_size++] = cur_node_id;
+		    order[order_size++] = node_id;
 		    continue;
 		}
-		processed[cur_node_id] = 1;
-		node = dglGetNode(graph, cur_node_id);
+		processed[node_id] = 1;
+
+		node = dglGetNode(graph, 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;
+		    if (!component[to]) {
+			component[to] = components;
+			/* do not go through closed nodes */
+			if (have_node_costs) {
+			    memcpy(&ncost, dglNodeGet_Attr(graph, dglEdgeGet_Tail(graph, edge)),
+				   sizeof(ncost));
+			}
+			if (ncost < 0)
+			    processed[to] = 1;
+
 			stack[stack_size++] = to;
 		    }
 		}
@@ -158,41 +235,48 @@ int NetA_strongly_connected_components(dglGraph_s * graph, int *component)
 
     dglNode_T_Release(&nt);
 
+    components = 0;
+    dglNode_T_Initialize(&nt, graph);
+
     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;
+	dglInt32_t cur_node_id = order[--order_size];
+	int cur_comp = component[cur_node_id];
+
+	if (cur_comp < 1) {
+	    component[cur_node_id] = ++components;
+	    stack[0] = cur_node_id;
+	    stack_size = 1;
+	    while (stack_size) {
+		dglInt32_t *node, *edgeset, *edge;
+		dglEdgesetTraverser_s et;
+		dglInt32_t node_id = stack[--stack_size];
+
+		node = dglGetNode(graph, node_id);
+		edgeset = dglNodeGet_InEdgeset(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_Head(graph, edge));
+		    if (component[to] == cur_comp) {
+			component[to] = components;
+			/* do not go through closed nodes */
+			if (have_node_costs) {
+			    memcpy(&ncost, dglNodeGet_Attr(graph, dglEdgeGet_Head(graph, edge)),
+				   sizeof(ncost));
+			}
+			if (ncost >= 0)
+			    stack[stack_size++] = to;
+		    }
 		}
+		dglEdgeset_T_Release(&et);
 	    }
-	    dglEdgeset_T_Release(&et);
 	}
     }
+    dglNode_T_Release(&nt);
 
     G_free(stack);
-    G_free(visited);
     G_free(order);
     G_free(processed);
     return components;

lib/vector/neta/path.c

@@ -11,6 +11,7 @@
    (>=v2). Read the file COPYING that comes with GRASS for details.
 
    \author Daniel Bundala (Google Summer of Code 2009)
+   \author Markus Metz
  */
 
 #include <stdio.h>
@@ -22,16 +23,16 @@
 #include <grass/neta.h>
 
 /*!
-   \brief Computes shortests paths to every node from nodes in "from".
+   \brief Computes shortest paths to every node from nodes in "from".
 
-   Array "dst" contains the length of the path or -1 if the node is not
+   Array "dst" contains the cost 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
+   \param[out] dst array of costs to reach nodes
+   \param[out] prev array of edges from predecessor along the shortest path
 
    \return 0 on success
    \return -1 on failure
@@ -41,6 +42,8 @@ int NetA_distance_from_points(dglGraph_s *graph, struct ilist *from,
 {
     int i, nnodes;
     dglHeap_s heap;
+    int have_node_costs;
+    dglInt32_t ncost;
 
     nnodes = dglGet_NodeCount(graph);
     dglEdgesetTraverser_s et;
@@ -51,13 +54,16 @@ int NetA_distance_from_points(dglGraph_s *graph, struct ilist *from,
 	prev[i] = NULL;
     }
 
+    ncost = 0;
+    have_node_costs = dglGet_NodeAttrSize(graph);
+
     dglHeapInit(&heap);
 
     for (i = 0; i < from->n_values; i++) {
 	int v = from->value[i];
 
 	if (dst[v] == 0)
-	    continue;		/*ingore duplicates */
+	    continue;		/* ignore duplicates */
 	dst[v] = 0;		/* make sure all from nodes are processed first */
 	dglHeapData_u heap_data;
 
@@ -68,6 +74,8 @@ int NetA_distance_from_points(dglGraph_s *graph, struct ilist *from,
 	dglInt32_t v, dist;
 	dglHeapNode_s heap_node;
 	dglHeapData_u heap_data;
+	dglInt32_t *edge;
+	dglInt32_t *node;
 
 	if (!dglHeapExtractMin(&heap, &heap_node))
 	    break;
@@ -76,11 +84,20 @@ int NetA_distance_from_points(dglGraph_s *graph, struct ilist *from,
 	if (dst[v] < dist)
 	    continue;
 
-	dglInt32_t *edge;
+	node = dglGetNode(graph, v);
+
+	if (have_node_costs && prev[v]) {
+	    memcpy(&ncost, dglNodeGet_Attr(graph, node),
+		   sizeof(ncost));
+	    if (ncost > 0)
+		dist += ncost;
+	    /* do not go through closed nodes */
+	    if (ncost < 0)
+		continue;
+	}
 
 	dglEdgeset_T_Initialize(&et, graph,
-				dglNodeGet_OutEdgeset(graph,
-						      dglGetNode(graph, v)));
+				dglNodeGet_OutEdgeset(graph, node));
 
 	for (edge = dglEdgeset_T_First(&et); edge;
 	     edge = dglEdgeset_T_Next(&et)) {
@@ -105,16 +122,123 @@ int NetA_distance_from_points(dglGraph_s *graph, struct ilist *from,
 }
 
 /*!
-   \brief Find a path (minimum number of edges) from 'from' to 'to' using only edges in 'edges'.
+   \brief Computes shortest paths from every node to nodes in "to".
+
+   Array "dst" contains the cost of the path or -1 if the node is not
+   reachable. Nxt contains edges from successor along the shortest
+   path. This method does reverse search starting with "to" nodes and 
+   going backward.
+
+   \param graph input graph
+   \param to list of 'to' positions
+   \param[out] dst array of costs to reach nodes
+   \param[out] nxt array of edges from successor along the shortest path
+
+   \return 0 on success
+   \return -1 on failure
+ */
+int NetA_distance_to_points(dglGraph_s *graph, struct ilist *to,
+			      int *dst, dglInt32_t **nxt)
+{
+    int i, nnodes;
+    dglHeap_s heap;
+    dglEdgesetTraverser_s et;
+    int have_node_costs;
+    dglInt32_t ncost;
+
+    nnodes = dglGet_NodeCount(graph);
+
+    /* initialize costs and edge list */
+    for (i = 1; i <= nnodes; i++) {
+	dst[i] = -1;
+	nxt[i] = NULL;
+    }
+
+    if (graph->Version < 2) {
+	G_warning("Directed graph must be version 2 or 3 for NetA_distance_to_points()");
+	return -1;
+    }
+
+    ncost = 0;
+    have_node_costs = dglGet_NodeAttrSize(graph);
+
+    dglHeapInit(&heap);
+
+    for (i = 0; i < to->n_values; i++) {
+	int v = to->value[i];
+
+	if (dst[v] == 0)
+	    continue;		/* ignore duplicates */
+	dst[v] = 0;		/* make sure all to nodes are processed first */
+	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;
+	dglInt32_t *edge;
+	dglInt32_t *node;
+
+	if (!dglHeapExtractMin(&heap, &heap_node))
+	    break;
+	v = heap_node.value.ul;
+	dist = heap_node.key;
+	if (dst[v] < dist)
+	    continue;
+
+	node = dglGetNode(graph, v);
+
+	if (have_node_costs && nxt[v]) {
+	    memcpy(&ncost, dglNodeGet_Attr(graph, node),
+		   sizeof(ncost));
+	    if (ncost > 0)
+		dist += ncost;
+	    /* do not go through closed nodes */
+	    if (ncost < 0)
+		continue;
+	}
+
+	dglEdgeset_T_Initialize(&et, graph,
+				dglNodeGet_InEdgeset(graph, node));
+
+	for (edge = dglEdgeset_T_First(&et); edge;
+	     edge = dglEdgeset_T_Next(&et)) {
+	    dglInt32_t *from = dglEdgeGet_Head(graph, edge);
+	    dglInt32_t from_id = dglNodeGet_Id(graph, from);
+	    dglInt32_t d = dglEdgeGet_Cost(graph, edge);
+
+	    if (dst[from_id] < 0 || dst[from_id] > dist + d) {
+		dst[from_id] = dist + d;
+		nxt[from_id] = edge;
+		heap_data.ul = from_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 flagged as valid in 'edges'. Edge costs are not 
+   considered. Closed nodes are not traversed.
 
    Precisely, edge with id I is used if edges[abs(i)] == 1. List
-   stores the indices of lines on the path. Method return number of
-   edges or -1 if no path exist.
+   stores the indices of lines on the path. The method returns the 
+   number of edges or -1 if no path exists.
 
    \param graph input graph
    \param from 'from' position
    \param to 'to' position
-   \param edges list of available edges
+   \param edges array of edges indicating wether an edge should be used
    \param[out] list list of edges
 
    \return number of edges
@@ -127,6 +251,8 @@ int NetA_find_path(dglGraph_s * graph, int from, int to, int *edges,
     dglEdgesetTraverser_s et;
     char *vis;
     int begin, end, cur, nnodes;
+    int have_node_costs;
+    dglInt32_t ncost;
 
     nnodes = dglGet_NodeCount(graph);
     prev = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
@@ -138,6 +264,9 @@ int NetA_find_path(dglGraph_s * graph, int from, int to, int *edges,
     }
     Vect_reset_list(list);
 
+    ncost = 0;
+    have_node_costs = dglGet_NodeAttrSize(graph);
+
     begin = 0;
     end = 1;
     vis[from] = 'y';
@@ -145,22 +274,32 @@ int NetA_find_path(dglGraph_s * graph, int from, int to, int *edges,
     prev[from] = NULL;
     while (begin != end) {
 	dglInt32_t vertex = queue[begin++];
+	dglInt32_t *edge, *node;
 
 	if (vertex == to)
 	    break;
-	dglInt32_t *edge, *node = dglGetNode(graph, vertex);
+
+	/* do not go through closed nodes */
+	if (have_node_costs && prev[vertex]) {
+	    memcpy(&ncost, dglNodeGet_Attr(graph, dglEdgeGet_Tail(graph, edge)),
+		   sizeof(ncost));
+	    if (ncost < 0)
+		continue;
+	}
+
+	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 =
+	    dglInt32_t edge_id = abs(dglEdgeGet_Id(graph, edge));
+	    dglInt32_t node_id =
 		dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
-	    if (edges[id] && !vis[to]) {
-		vis[to] = 'y';
-		prev[to] = edge;
-		queue[end++] = to;
+	    if (edges[edge_id] && !vis[node_id]) {
+		vis[node_id] = 'y';
+		prev[node_id] = edge;
+		queue[end++] = node_id;
 	    }
 	}
 	dglEdgeset_T_Release(&et);

lib/vector/neta/spanningtree.c

@@ -136,6 +136,7 @@ int NetA_spanning_tree(dglGraph_s * graph, struct ilist *tree_list)
 	    Vect_list_append(tree_list, dglEdgeGet_Id(graph, perm[i].edge));
 	}
     }
+    G_percent(index, index, 1);
     G_free(perm);
     uf_release(&uf);
     return edges;

lib/vector/neta/utils.c

@@ -1,5 +1,5 @@
 /*!
-   \file vector/neta/timetables.c
+   \file vector/neta/utils.c
 
    \brief Network Analysis library - utils
 
@@ -91,8 +91,8 @@ void NetA_points_to_nodes(struct Map_info *In, struct ilist *point_list)
    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)
+   node_costs are multiplied by the graph's cost multiplier and  
+   truncated to integers (as is done in Vect_net_build_graph)
 
    \param In pointer to Map_info structure
    \param layer layer number
@@ -141,8 +141,12 @@ int NetA_get_node_costs(struct Map_info *In, int layer, char *column,
 	    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;
+	    if (db_CatValArray_get_value_double(&vals, cat, &value) == DB_OK) {
+		if (value < 0)
+		    node_costs[node] = -1;
+		else
+		    node_costs[node] = value * In->dgraph.cost_multip;
+	    }
 	}
     }
 
@@ -159,12 +163,13 @@ int NetA_get_node_costs(struct Map_info *In, int layer, char *column,
    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.
+   unitialised. Nodes_to_features will be wrong if several lines connect 
+   to the same node.
 
    \param map pointer to Map_info structure
    \param varray pointer to varray structure
    \param[out] nodes list of node ids
-   \param node_to_features ?
+   \param[out] nodes_to_features maps nodes to features
  */
 void NetA_varray_to_nodes(struct Map_info *map, struct varray *varray,
 			  struct ilist *nodes, int *nodes_to_features)