netalib: fix network components

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

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,42 +235,47 @@ 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;
 }