[GRASS-SVN] r38293 - in grass-addons/vector/net.analyze: . netalib
v.net.allpairs v.net.centrality v.net.connectivity
v.net.distance v.net.flow
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Jul 7 09:21:27 EDT 2009
Author: dano
Date: 2009-07-07 09:21:27 -0400 (Tue, 07 Jul 2009)
New Revision: 38293
Added:
grass-addons/vector/net.analyze/netalib/centrality.c
grass-addons/vector/net.analyze/netalib/path.c
grass-addons/vector/net.analyze/v.net.centrality/
grass-addons/vector/net.analyze/v.net.centrality/Makefile
grass-addons/vector/net.analyze/v.net.centrality/main.c
grass-addons/vector/net.analyze/v.net.connectivity/
grass-addons/vector/net.analyze/v.net.connectivity/Makefile
grass-addons/vector/net.analyze/v.net.connectivity/main.c
grass-addons/vector/net.analyze/v.net.distance/
grass-addons/vector/net.analyze/v.net.distance/Makefile
grass-addons/vector/net.analyze/v.net.distance/main.c
Modified:
grass-addons/vector/net.analyze/Makefile
grass-addons/vector/net.analyze/netalib/flow.c
grass-addons/vector/net.analyze/netalib/neta.h
grass-addons/vector/net.analyze/netalib/utils.c
grass-addons/vector/net.analyze/v.net.allpairs/main.c
grass-addons/vector/net.analyze/v.net.flow/main.c
Log:
New modules v.net.distance/connectivity/centrality and several bug fixes/updates of other modules as well
Modified: grass-addons/vector/net.analyze/Makefile
===================================================================
--- grass-addons/vector/net.analyze/Makefile 2009-07-07 11:38:23 UTC (rev 38292)
+++ grass-addons/vector/net.analyze/Makefile 2009-07-07 13:21:27 UTC (rev 38293)
@@ -5,7 +5,10 @@
v.net.bridge \
v.net.spanningtree \
v.net.allpairs \
- v.net.flow
+ v.net.flow \
+ v.net.connectivity \
+ v.net.centrality \
+ v.net.distance
SUBDIRS = netalib $(SUBDIRS1)
Added: grass-addons/vector/net.analyze/netalib/centrality.c
===================================================================
--- grass-addons/vector/net.analyze/netalib/centrality.c (rev 0)
+++ grass-addons/vector/net.analyze/netalib/centrality.c 2009-07-07 13:21:27 UTC (rev 38293)
@@ -0,0 +1,212 @@
+
+/****************************************************************
+ *
+ * MODULE: netalib
+ *
+ * AUTHOR(S): Daniel Bundala
+ *
+ * PURPOSE: Centrality measures
+ *
+ *
+ * COPYRIGHT: (C) 2002-2005 by 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.
+ *
+ ****************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+#include <grass/neta.h>
+
+/* Computes degree centrality measure. Array degree has to be properly initialised to nnodes+1 elements */
+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;
+}
+
+/* Computes eigenvector centrality using edge costs as weights. Returns 0 on success. */
+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;
+}
+
+/*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.
+ * Returns 0 on success */
+int neta_betweenness_closeness(dglGraph_s * graph, double *betweenness,
+ double *closeness)
+{
+ int i, j, nnodes, stack_size;
+ 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;
+ }
+
+
+ dglNode_T_Initialize(&nt, graph);
+ for (node = dglNode_T_First(&nt); node; node = dglNode_T_Next(&nt)) {
+ dglInt32_t s = dglNodeGet_Id(graph, node);
+ dglHeapData_u heap_data;
+ dglHeapNode_s heap_node;
+
+ stack_size = 0;
+ for (i = 1; i <= nnodes; i++)
+ Vect_reset_list(prev[i]);
+ for (i = 1; i <= nnodes; i++) {
+ cnt[i] = 0;
+ dst[i] = -1;
+ }
+ dst[s] = 0;
+ cnt[s] = 1;
+ dglHeapInit(&heap);
+ heap_data.ul = s;
+ dglHeapInsertMin(&heap, 0, ' ', heap_data);
+ while (1) {
+ dglInt32_t v, dist;
+ if (!dglHeapExtractMin(&heap, &heap_node))
+ break;
+ v = heap_node.value.ul;
+ dist = heap_node.key;
+ if (dst[v] < dist)
+ continue;
+ stack[stack_size++] = v;
+
+ dglInt32_t *edge;
+ dglEdgeset_T_Initialize(&et, graph,
+ dglNodeGet_OutEdgeset(graph,
+ dglGetNode(graph,
+ v)));
+ for (edge = dglEdgeset_T_First(&et); edge;
+ edge = dglEdgeset_T_Next(&et)) {
+ dglInt32_t *to = dglEdgeGet_Tail(graph, edge);
+ dglInt32_t to_id = dglNodeGet_Id(graph, to);
+ dglInt32_t d = dglEdgeGet_Cost(graph, edge);
+ if (dst[to_id] == -1 || dst[to_id] > dist + d) {
+ dst[to_id] = dist + d;
+ Vect_reset_list(prev[to_id]);
+ heap_data.ul = to_id;
+ dglHeapInsertMin(&heap, dist + d, ' ', heap_data);
+ }
+ if (dst[to_id] == dist + d) {
+ cnt[to_id] += cnt[v];
+ Vect_list_append(prev[to_id], v);
+ }
+ }
+
+ dglEdgeset_T_Release(&et);
+ }
+ dglHeapFree(&heap, NULL);
+ for (i = 1; i <= nnodes; i++)
+ delta[i] = 0;
+ for (i = stack_size - 1; i >= 0; i--) {
+ dglInt32_t w = stack[i];
+ if (closeness)
+ closeness[s] += dst[w];
+
+ for (j = 0; j < prev[w]->n_values; j++) {
+ dglInt32_t v = prev[w]->value[j];
+ delta[v] += (cnt[v] / (double)cnt[w]) * (1.0 + delta[w]);
+ }
+ if (w != s && betweenness)
+ betweenness[w] += delta[w];
+
+ }
+ if (closeness)
+ closeness[s] /= (double)stack_size;
+ }
+ dglNode_T_Release(&nt);
+
+ for (i = 1; i <= nnodes; i++)
+ Vect_destroy_list(prev[i]);
+ G_free(delta);
+ G_free(cnt);
+ G_free(stack);
+ G_free(prev);
+ G_free(dst);
+
+ return 0;
+};
Property changes on: grass-addons/vector/net.analyze/netalib/centrality.c
___________________________________________________________________
Added: svn:executable
+ *
Modified: grass-addons/vector/net.analyze/netalib/flow.c
===================================================================
--- grass-addons/vector/net.analyze/netalib/flow.c 2009-07-07 11:38:23 UTC (rev 38292)
+++ grass-addons/vector/net.analyze/netalib/flow.c 2009-07-07 13:21:27 UTC (rev 38293)
@@ -23,6 +23,7 @@
#include <grass/Vect.h>
#include <grass/glocale.h>
#include <grass/dgl/graph.h>
+#include <grass/neta.h>
dglInt32_t sign(dglInt32_t x)
{
@@ -36,35 +37,46 @@
* each edge. Negative flow corresponds to a flow in opposite direction
* The function assumes that the edge costs correspond to edge capacities.
*/
-int neta_flow(dglGraph_s * graph, int source, int sink, int *flow)
+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 *));
- if (!queue || !prev) {
+ 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 = 0;
- begin = 0;
- end = 1;
- queue[0] = source;
+ 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) {
+ while (begin != end && found == -1) {
dglInt32_t vertex = queue[begin++];
dglInt32_t *edge, *node = dglGetNode(graph, vertex);
dglEdgeset_T_Initialize(&et, graph,
@@ -75,11 +87,11 @@
dglInt32_t id = dglEdgeGet_Id(graph, edge);
dglInt32_t to =
dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
- if (to != source && prev[to] == NULL &&
+ if (!is_source[to] && prev[to] == NULL &&
cap > sign(id) * flow[abs(id)]) {
prev[to] = edge;
- if (to == sink) {
- found = 1;
+ if (is_sink[to]) {
+ found = to;
break;
}
queue[end++] = to;
@@ -87,15 +99,15 @@
}
dglEdgeset_T_Release(&et);
}
- if (!found)
+ if (found == -1)
break; /*no augmenting path */
/*find minimum residual capacity along the augmenting path */
- node = sink;
+ node = found;
edge_id = dglEdgeGet_Id(graph, prev[node]);
min_residue =
dglEdgeGet_Cost(graph,
prev[node]) - sign(edge_id) * flow[abs(edge_id)];
- while (node != source) {
+ while (!is_source[node]) {
dglInt32_t residue;
edge_id = dglEdgeGet_Id(graph, prev[node]);
residue =
@@ -108,8 +120,8 @@
}
total_flow += min_residue;
/*update flow along the augmenting path */
- node = sink;
- while (node != source) {
+ 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]));
@@ -118,6 +130,141 @@
G_free(queue);
G_free(prev);
+ G_free(is_source);
+ G_free(is_sink);
return total_flow;
}
+
+/*
+ * This methods calculates minimum cut between source(s) and sink(s) and returns the list of edges (cut).
+ * 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.
+ */
+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;
+}
+
+/* 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.
+ *
+ * Method returns the number of undirected edges in the graph.
+ */
+
+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;
+}
Modified: grass-addons/vector/net.analyze/netalib/neta.h
===================================================================
--- grass-addons/vector/net.analyze/netalib/neta.h 2009-07-07 11:38:23 UTC (rev 38292)
+++ grass-addons/vector/net.analyze/netalib/neta.h 2009-07-07 13:21:27 UTC (rev 38293)
@@ -23,6 +23,7 @@
#include <stdlib.h>
#include <grass/gis.h>
#include <grass/Vect.h>
+#include <grass/dbmi.h>
#include <grass/glocale.h>
#include <grass/dgl/graph.h>
@@ -44,5 +45,21 @@
void neta_add_point_on_node(struct Map_info *In, struct Map_info *Out, int node, struct line_cats *Cats);
/*neta_flow.c*/
-int neta_flow(dglGraph_s *graph, int source, int sink, int *flow);
+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, VARRAY *varray, struct ilist *nodes, int *nodes_to_features);
+/*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);
#endif
Added: grass-addons/vector/net.analyze/netalib/path.c
===================================================================
--- grass-addons/vector/net.analyze/netalib/path.c (rev 0)
+++ grass-addons/vector/net.analyze/netalib/path.c 2009-07-07 13:21:27 UTC (rev 38293)
@@ -0,0 +1,87 @@
+
+/****************************************************************
+ *
+ * MODULE: netalib
+ *
+ * AUTHOR(S): Daniel Bundala
+ *
+ * PURPOSE: Shortest paths from a set of nodes
+ *
+ *
+ * COPYRIGHT: (C) 2002-2005 by 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.
+ *
+ ****************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+#include <grass/dgl/graph.h>
+#include <grass/neta.h>
+
+/* 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 */
+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;
+}
Property changes on: grass-addons/vector/net.analyze/netalib/path.c
___________________________________________________________________
Added: svn:executable
+ *
Modified: grass-addons/vector/net.analyze/netalib/utils.c
===================================================================
--- grass-addons/vector/net.analyze/netalib/utils.c 2009-07-07 11:38:23 UTC (rev 38292)
+++ grass-addons/vector/net.analyze/netalib/utils.c 2009-07-07 13:21:27 UTC (rev 38293)
@@ -25,6 +25,7 @@
#include <grass/neta.h>
+/* Writes GV_POINT to Out at the position of the node in In */
void neta_add_point_on_node(struct Map_info *In, struct Map_info *Out, int node,
struct line_cats *Cats)
{
@@ -37,3 +38,114 @@
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);
+ * }
+ */
+
+/* Finds the node corresponding to each point in the point_list */
+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;
+ }
+}
+
+/* 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)
+ *
+ * Returns 1 on success, 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;
+}
+
+/*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 */
+void neta_varray_to_nodes(struct Map_info *map, 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);
+ 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);
+ 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);
+ nodes_to_features[node1] = nodes_to_features[node2] = i;
+ }
+ }
+}
Modified: grass-addons/vector/net.analyze/v.net.allpairs/main.c
===================================================================
--- grass-addons/vector/net.analyze/v.net.allpairs/main.c 2009-07-07 11:38:23 UTC (rev 38292)
+++ grass-addons/vector/net.analyze/v.net.allpairs/main.c 2009-07-07 13:21:27 UTC (rev 38293)
@@ -234,7 +234,8 @@
for (j = 1; j <= nnodes; j++)
if (cats[j] != -1) {
sprintf(buf, "insert into %s values (%d, %d, %f)",
- Fi->table, i, j, dist[i][j] / 1000.0);
+ Fi->table, i, j,
+ dist[i][j] / (double)In.cost_multip);
db_set_string(&sql, buf);
G_debug(3, db_get_string(&sql));
Added: grass-addons/vector/net.analyze/v.net.centrality/Makefile
===================================================================
--- grass-addons/vector/net.analyze/v.net.centrality/Makefile (rev 0)
+++ grass-addons/vector/net.analyze/v.net.centrality/Makefile 2009-07-07 13:21:27 UTC (rev 38293)
@@ -0,0 +1,15 @@
+MODULE_TOPDIR = ../../..
+
+PGM=v.net.centrality
+
+include ../Netalib.make
+
+LIBES = $(NETALIB) $(VECTLIB) $(VECTLIB_REAL) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(NETADEP) $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Property changes on: grass-addons/vector/net.analyze/v.net.centrality/Makefile
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/vector/net.analyze/v.net.centrality/main.c
===================================================================
--- grass-addons/vector/net.analyze/v.net.centrality/main.c (rev 0)
+++ grass-addons/vector/net.analyze/v.net.centrality/main.c 2009-07-07 13:21:27 UTC (rev 38293)
@@ -0,0 +1,297 @@
+
+/****************************************************************
+ *
+ * MODULE: v.net.centrality
+ *
+ * AUTHOR(S): Daniel Bundala
+ *
+ * PURPOSE: This module computes various cetrality measures
+ *
+ * COPYRIGHT: (C) 2002-2005 by 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.
+ *
+ ****************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+#include <grass/dbmi.h>
+#include <grass/neta.h>
+
+void update_record(dbDriver * driver, dbString * sql, char *table, char *key,
+ char *column, double value, int cat)
+{
+ char buf[2000];
+ sprintf(buf, "update %s set %s=%f where %s=%d", table, column, value, key,
+ cat);
+ db_set_string(sql, buf);
+ G_debug(3, db_get_string(sql));
+ if (db_execute_immediate(driver, sql) != DB_OK) {
+ db_close_database_shutdown_driver(driver);
+ G_fatal_error(_("Cannot insert new record: %s"), db_get_string(sql));
+ };
+}
+
+int main(int argc, char *argv[])
+{
+ struct Map_info In;
+ struct line_cats *Cats;
+ char *mapset;
+ struct GModule *module; /* GRASS module for parsing arguments */
+ struct Option *map_in;
+ struct Option *cat_opt, *field_opt, *where_opt, *abcol, *afcol;
+ struct Option *deg_opt, *close_opt, *betw_opt, *eigen_opt;
+ struct Option *iter_opt, *error_opt;
+ struct Flag *geo_f;
+ int chcat, with_z;
+ int layer, mask_type;
+ VARRAY *varray;
+ dglGraph_s *graph;
+ int i, geo, nnodes, nlines;
+ double *deg, *close, *betw, *eigen;
+
+ /* Attribute table */
+ dbString sql;
+ dbDriver *driver;
+ struct field_info *Fi;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ module->keywords =
+ _
+ ("network, centrality measures, degree, betweeness, closeness, eigenvector");
+ module->description =
+ _
+ ("Computes degree, centrality, betweeness, closeness and eigenvector cetrality measures");
+
+ /* Define the different options as defined in gis.h */
+ map_in = G_define_standard_option(G_OPT_V_INPUT);
+
+ field_opt = G_define_standard_option(G_OPT_V_FIELD);
+ cat_opt = G_define_standard_option(G_OPT_V_CATS);
+ where_opt = G_define_standard_option(G_OPT_WHERE);
+
+ afcol = G_define_option();
+ afcol->key = "afcolumn";
+ afcol->type = TYPE_STRING;
+ afcol->required = NO;
+ afcol->description = _("Arc forward/both direction(s) cost column");
+
+ abcol = G_define_option();
+ abcol->key = "abcolumn";
+ abcol->type = TYPE_STRING;
+ abcol->required = NO;
+ abcol->description = _("Arc backward direction cost column");
+
+ deg_opt = G_define_option();
+ deg_opt->key = "degree";
+ deg_opt->type = TYPE_STRING;
+ deg_opt->required = NO;
+ deg_opt->description = _("Degree centrality column");
+
+ close_opt = G_define_option();
+ close_opt->key = "closeness";
+ close_opt->type = TYPE_STRING;
+ close_opt->required = NO;
+ close_opt->description = _("Closeness centrality column");
+
+ betw_opt = G_define_option();
+ betw_opt->key = "betweenness";
+ betw_opt->type = TYPE_STRING;
+ betw_opt->required = NO;
+ betw_opt->description = _("Betweenness centrality column");
+
+ eigen_opt = G_define_option();
+ eigen_opt->key = "eigenvector";
+ eigen_opt->type = TYPE_STRING;
+ eigen_opt->required = NO;
+ eigen_opt->description = _("Eigenvector centrality column");
+
+ iter_opt = G_define_option();
+ iter_opt->key = "iterations";
+ iter_opt->answer = "1000";
+ iter_opt->type = TYPE_INTEGER;
+ iter_opt->required = NO;
+ iter_opt->description =
+ _("Maximum number of iterations to compute eigenvector centrality");
+
+ error_opt = G_define_option();
+ error_opt->key = "error";
+ error_opt->answer = "0.1";
+ error_opt->type = TYPE_DOUBLE;
+ error_opt->required = NO;
+ error_opt->description =
+ _("Cuumulative error tolerance for eigenvector centrality");
+
+
+ geo_f = G_define_flag();
+ geo_f->key = 'g';
+ geo_f->description =
+ _("Use geodesic calculation for longitude-latitude locations");
+
+ /* options and flags parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+ /* TODO: make an option for this */
+ mask_type = GV_LINE | GV_BOUNDARY;
+
+ Cats = Vect_new_cats_struct();
+
+ if ((mapset = G_find_vector2(map_in->answer, "")) == NULL)
+ G_fatal_error(_("Vector map <%s> not found"), map_in->answer);
+
+ Vect_set_open_level(2);
+
+ if (1 > Vect_open_old(&In, map_in->answer, mapset))
+ G_fatal_error(_("Unable to open vector map <%s>"),
+ G_fully_qualified_name(map_in->answer, mapset));
+
+ with_z = Vect_is_3d(&In);
+
+ if (geo_f->answer) {
+ geo = 1;
+ if (G_projection() != PROJECTION_LL)
+ G_warning(_("The current projection is not longitude-latitude"));
+ }
+ else
+ geo = 0;
+
+ /* parse filter option and select appropriate lines */
+ layer = atoi(field_opt->answer);
+ if (where_opt->answer) {
+ if (layer < 1)
+ G_fatal_error(_("'%s' must be > 0 for '%s'"), "layer", "where");
+ if (cat_opt->answer)
+ G_warning(_
+ ("'where' and 'cats' parameters were supplied, cat will be ignored"));
+ chcat = 1;
+ varray = Vect_new_varray(Vect_get_num_lines(&In));
+ if (Vect_set_varray_from_db
+ (&In, layer, where_opt->answer, mask_type, 1, varray) == -1) {
+ G_warning(_("Unable to load data from database"));
+ }
+ }
+ else if (cat_opt->answer) {
+ if (layer < 1)
+ G_fatal_error(_("'%s' must be > 0 for '%s'"), "layer", "cat");
+ varray = Vect_new_varray(Vect_get_num_lines(&In));
+ chcat = 1;
+ if (Vect_set_varray_from_cat_string
+ (&In, layer, cat_opt->answer, mask_type, 1, varray) == -1) {
+ G_warning(_("Problem loading category values"));
+ }
+ }
+ else {
+ chcat = 0;
+ varray = NULL;
+ }
+
+ /* Open database */
+ 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);
+ db_begin_transaction(driver);
+ db_init_string(&sql);
+
+ Vect_net_build_graph(&In, mask_type, atoi(field_opt->answer), 0,
+ afcol->answer, abcol->answer, NULL, geo, 0);
+ graph = &(In.graph);
+ nnodes = dglGet_NodeCount(graph);
+
+ deg = close = betw = eigen = NULL;
+
+ if (deg_opt->answer) {
+ deg = (double *)G_calloc(nnodes + 1, sizeof(double));
+ if (!deg)
+ G_fatal_error(_("Out of memory"));
+ }
+
+ if (close_opt->answer) {
+ close = (double *)G_calloc(nnodes + 1, sizeof(double));
+ if (!close)
+ G_fatal_error(_("Out of memory"));
+ }
+
+ if (betw_opt->answer) {
+ betw = (double *)G_calloc(nnodes + 1, sizeof(double));
+ if (!betw)
+ G_fatal_error(_("Out of memory"));
+ }
+
+ if (eigen_opt->answer) {
+ eigen = (double *)G_calloc(nnodes + 1, sizeof(double));
+ if (!eigen)
+ G_fatal_error(_("Out of memory"));
+ }
+
+
+ if (deg_opt->answer) {
+ G_message(_("Computing degree centrality measure"));
+ neta_degree_centrality(graph, deg);
+ }
+ if (betw_opt->answer || close_opt->answer) {
+ G_message(_
+ ("Computing betweenness and/or closeness centrality measure"));
+ neta_betweenness_closeness(graph, betw, close);
+ }
+ if (eigen_opt->answer) {
+ G_message(_("Computing eigenvector centrality measure"));
+ neta_eigenvector_centrality(graph, atoi(iter_opt->answer),
+ atof(error_opt->answer), eigen);
+ }
+
+ nlines = Vect_get_num_lines(&In);
+ G_message(_("Writing data into the table..."));
+ G_percent_reset();
+ for (i = 1; i <= nlines; i++) {
+ G_percent(i, nlines, 1);
+ // if(!varray->c[i])continue;
+ int type = Vect_read_line(&In, NULL, Cats, i);
+ if (type == GV_POINT) {
+ int cat, node;
+ if (!Vect_cat_get(Cats, layer, &cat))
+ continue;
+ Vect_get_line_nodes(&In, i, &node, NULL);
+ if (deg_opt->answer)
+ update_record(driver, &sql, Fi->table, Fi->key, deg_opt->answer,
+ deg[node], cat);
+ if (close_opt->answer)
+ update_record(driver, &sql, Fi->table, Fi->key,
+ close_opt->answer, close[node], cat);
+ if (betw_opt->answer)
+ update_record(driver, &sql, Fi->table, Fi->key,
+ betw_opt->answer, betw[node], cat);
+ if (eigen_opt->answer)
+ update_record(driver, &sql, Fi->table, Fi->key,
+ eigen_opt->answer, eigen[node], cat);
+ }
+
+
+ }
+ db_commit_transaction(driver);
+ db_close_database_shutdown_driver(driver);
+
+ if (deg)
+ G_free(deg);
+ if (close)
+ G_free(close);
+ if (betw)
+ G_free(betw);
+ if (eigen)
+ G_free(eigen);
+
+ Vect_close(&In);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/vector/net.analyze/v.net.centrality/main.c
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/vector/net.analyze/v.net.connectivity/Makefile
===================================================================
--- grass-addons/vector/net.analyze/v.net.connectivity/Makefile (rev 0)
+++ grass-addons/vector/net.analyze/v.net.connectivity/Makefile 2009-07-07 13:21:27 UTC (rev 38293)
@@ -0,0 +1,15 @@
+MODULE_TOPDIR = ../../..
+
+PGM=v.net.connectivity
+
+include ../Netalib.make
+
+LIBES = $(NETALIB) $(VECTLIB) $(VECTLIB_REAL) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(NETADEP) $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Property changes on: grass-addons/vector/net.analyze/v.net.connectivity/Makefile
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/vector/net.analyze/v.net.connectivity/main.c
===================================================================
--- grass-addons/vector/net.analyze/v.net.connectivity/main.c (rev 0)
+++ grass-addons/vector/net.analyze/v.net.connectivity/main.c 2009-07-07 13:21:27 UTC (rev 38293)
@@ -0,0 +1,221 @@
+
+/****************************************************************
+ *
+ * MODULE: v.net.connectivity
+ *
+ * AUTHOR(S): Daniel Bundala
+ *
+ * PURPOSE: Vertex connectivity between two sets of nodes
+ *
+ * COPYRIGHT: (C) 2002-2005 by 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.
+ *
+ ****************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+#include <grass/dbmi.h>
+#include <grass/neta.h>
+
+
+int main(int argc, char *argv[])
+{
+ struct Map_info In, Out;
+ static struct line_pnts *Points;
+ struct line_cats *Cats;
+ char *mapset;
+ struct GModule *module; /* GRASS module for parsing arguments */
+ struct Option *map_in, *map_out;
+ struct Option *cat_opt, *field_opt, *where_opt, *ncol, *set1_opt, *set2_opt;
+ int chcat, with_z;
+ int layer, mask_type;
+ VARRAY *varray;
+ dglGraph_s *graph;
+ int i, nnodes, nlines, *flow, total_flow, nedges;
+ struct ilist *set1_list, *set2_list, *cut;
+ struct cat_list *set1_cats, *set2_cats;
+ int *node_costs;
+
+ dglGraph_s vg;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ module->keywords = _("network, connectivity, vertex connectivity");
+ module->description =
+ _("Computes vertex connectivity between two sets of nodes.");
+
+ /* Define the different options as defined in gis.h */
+ map_in = G_define_standard_option(G_OPT_V_INPUT);
+ map_out = G_define_standard_option(G_OPT_V_OUTPUT);
+
+ field_opt = G_define_standard_option(G_OPT_V_FIELD);
+ cat_opt = G_define_standard_option(G_OPT_V_CATS);
+ where_opt = G_define_standard_option(G_OPT_WHERE);
+
+ ncol = G_define_option();
+ ncol->key = "ncolumn";
+ ncol->type = TYPE_STRING;
+ ncol->required = NO;
+ ncol->description = _("Node capacity column");
+
+ set1_opt = G_define_standard_option(G_OPT_V_CATS);
+ set1_opt->key = "set1";
+ set1_opt->required = YES;
+ set1_opt->description = _("Categories of the first set");
+
+ set2_opt = G_define_standard_option(G_OPT_V_CATS);
+ set2_opt->key = "set2";
+ set2_opt->required = YES;
+ set2_opt->description = _("Categories of the second set");
+
+ /* options and flags parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+ /* TODO: make an option for this */
+ mask_type = GV_LINE | GV_BOUNDARY;
+
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ Vect_check_input_output_name(map_in->answer, map_out->answer,
+ GV_FATAL_EXIT);
+
+ if ((mapset = G_find_vector2(map_in->answer, "")) == NULL)
+ G_fatal_error(_("Vector map <%s> not found"), map_in->answer);
+
+ Vect_set_open_level(2);
+
+ if (1 > Vect_open_old(&In, map_in->answer, mapset))
+ G_fatal_error(_("Unable to open vector map <%s>"),
+ G_fully_qualified_name(map_in->answer, mapset));
+
+ with_z = Vect_is_3d(&In);
+
+ if (0 > Vect_open_new(&Out, map_out->answer, with_z)) {
+ Vect_close(&In);
+ G_fatal_error(_("Unable to create vector map <%s>"), map_out->answer);
+ }
+
+ /* parse filter option and select appropriate lines */
+ layer = atoi(field_opt->answer);
+ if (where_opt->answer) {
+ if (layer < 1)
+ G_fatal_error(_("'%s' must be > 0 for '%s'"), "layer", "where");
+ if (cat_opt->answer)
+ G_warning(_
+ ("'where' and 'cats' parameters were supplied, cat will be ignored"));
+ chcat = 1;
+ varray = Vect_new_varray(Vect_get_num_lines(&In));
+ if (Vect_set_varray_from_db
+ (&In, layer, where_opt->answer, mask_type, 1, varray) == -1) {
+ G_warning(_("Unable to load data from database"));
+ }
+ }
+ else if (cat_opt->answer) {
+ if (layer < 1)
+ G_fatal_error(_("'%s' must be > 0 for '%s'"), "layer", "cat");
+ varray = Vect_new_varray(Vect_get_num_lines(&In));
+ chcat = 1;
+ if (Vect_set_varray_from_cat_string
+ (&In, layer, cat_opt->answer, mask_type, 1, varray) == -1) {
+ G_warning(_("Problem loading category values"));
+ }
+ }
+ else {
+ chcat = 0;
+ varray = NULL;
+ }
+
+ set1_cats = Vect_new_cat_list();
+ set2_cats = Vect_new_cat_list();
+ Vect_str_to_cat_list(set1_opt->answer, set1_cats);
+ Vect_str_to_cat_list(set2_opt->answer, set2_cats);
+ set1_list = Vect_new_list();
+ set2_list = Vect_new_list();
+
+ nlines = Vect_get_num_lines(&In);
+ nnodes = Vect_get_num_nodes(&In);
+ for (i = 1; i <= nlines; i++) {
+ int type, cat;
+ type = Vect_read_line(&In, NULL, Cats, i);
+ if (type != GV_POINT)
+ continue;
+ Vect_cat_get(Cats, layer, &cat);
+ if (Vect_cat_in_cat_list(cat, set1_cats))
+ Vect_list_append(set1_list, i);
+ if (Vect_cat_in_cat_list(cat, set2_cats))
+ Vect_list_append(set2_list, i);
+ }
+
+ if (set1_list->n_values == 0)
+ G_fatal_error(_("No points with categories [%s]"), set1_opt->answer);
+
+ if (set2_list->n_values == 0)
+ G_fatal_error(_("No points with categories [%s]"), set2_opt->answer);
+
+ neta_points_to_nodes(&In, set1_list);
+ neta_points_to_nodes(&In, set2_list);
+
+ Vect_copy_head_data(&In, &Out);
+ Vect_hist_copy(&In, &Out);
+
+ Vect_hist_command(&Out);
+
+ Vect_net_build_graph(&In, mask_type, 0, atoi(field_opt->answer),
+ NULL, NULL, NULL, 0, 0);
+ graph = &(In.graph);
+
+ /*build new graph */
+ if (ncol->answer) {
+ node_costs = (int *)G_calloc(nnodes + 1, sizeof(int));
+ if (!node_costs)
+ G_fatal_error(_("Out of memory"));
+ neta_get_node_costs(&In, layer, ncol->answer, node_costs);
+ nedges = neta_split_vertices(graph, &vg, node_costs);
+ G_free(node_costs);
+ }
+ else
+ nedges = neta_split_vertices(graph, &vg, NULL);
+ graph = &vg;
+
+ for (i = 0; i < set1_list->n_values; i++)
+ set1_list->value[i] = set1_list->value[i] * 2; /*out vertex */
+ for (i = 0; i < set2_list->n_values; i++)
+ set2_list->value[i] = set2_list->value[i] * 2 - 1; /*in vertex */
+
+ flow = (int *)G_calloc(nedges + 1, sizeof(int));
+ if (!flow)
+ G_fatal_error(_("Out of memory"));
+
+ total_flow = neta_flow(graph, set1_list, set2_list, flow);
+ G_debug(3, "Connectivity: %d", total_flow);
+ cut = Vect_new_list();
+ total_flow = neta_min_cut(graph, set1_list, set2_list, flow, cut);
+
+ /*TODO: copy old points */
+ for (i = 0; i < cut->n_values; i++)
+ neta_add_point_on_node(&In, &Out, cut->value[i], Cats);
+
+ Vect_destroy_list(cut);
+
+ G_free(flow);
+ Vect_destroy_list(set1_list);
+ Vect_destroy_list(set2_list);
+
+ Vect_build(&Out);
+
+ Vect_close(&In);
+ Vect_close(&Out);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/vector/net.analyze/v.net.connectivity/main.c
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/vector/net.analyze/v.net.distance/Makefile
===================================================================
--- grass-addons/vector/net.analyze/v.net.distance/Makefile (rev 0)
+++ grass-addons/vector/net.analyze/v.net.distance/Makefile 2009-07-07 13:21:27 UTC (rev 38293)
@@ -0,0 +1,16 @@
+
+MODULE_TOPDIR = ../../..
+
+PGM=v.net.distance
+
+include ../Netalib.make
+
+LIBES = $(NETALIB) $(VECTLIB) $(VECTLIB_REAL) $(DBMILIB) $(GISLIB)
+DEPENDENCIES = $(NETADEP) $(VECTDEP) $(DBMIDEP) $(GISDEP)
+EXTRA_INC = $(VECT_INC)
+EXTRA_CFLAGS = $(VECT_CFLAGS)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+default: cmd
+
Property changes on: grass-addons/vector/net.analyze/v.net.distance/Makefile
___________________________________________________________________
Added: svn:executable
+ *
Added: grass-addons/vector/net.analyze/v.net.distance/main.c
===================================================================
--- grass-addons/vector/net.analyze/v.net.distance/main.c (rev 0)
+++ grass-addons/vector/net.analyze/v.net.distance/main.c 2009-07-07 13:21:27 UTC (rev 38293)
@@ -0,0 +1,308 @@
+
+/****************************************************************
+ *
+ * MODULE: v.net.distance
+ *
+ * AUTHOR(S): Daniel Bundala
+ *
+ * PURPOSE: Computes shortest distance via the network between
+ * the given sets of features.
+ *
+ * COPYRIGHT: (C) 2002-2005 by 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.
+ *
+ ****************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <grass/gis.h>
+#include <grass/Vect.h>
+#include <grass/glocale.h>
+#include <grass/dbmi.h>
+#include <grass/neta.h>
+
+int initiliase_varray(struct Map_info *In, int layer, int mask_type,
+ char *where, char *cat, 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 {
+ G_fatal_error(_("Neither 'where' nor 'cat' were specified"));
+ return 0;
+ }
+
+
+}
+
+int main(int argc, char *argv[])
+{
+ struct Map_info In, Out;
+ static struct line_pnts *Points;
+ struct line_cats *Cats;
+ char *mapset;
+ struct GModule *module; /* GRASS module for parsing arguments */
+ struct Option *map_in, *map_out, *abcol, *afcol;
+ struct Option *catf_opt, *fieldf_opt, *wheref_opt;
+ struct Option *catt_opt, *fieldt_opt, *wheret_opt;
+ struct Flag *geo_f, *newpoints_f;
+ int with_z, geo;
+ int mask_type;
+ VARRAY *varrayf, *varrayt;
+ int flayer, tlayer;
+ dglGraph_s *graph;
+ struct ilist *nodest;
+ int i, nnodes, nlines;
+ int *dst, *nodes_to_features;
+ dglInt32_t **prev;
+ struct line_cats **on_path;
+ char buf[2000];
+
+ /* Attribute table */
+ dbString sql;
+ dbDriver *driver;
+ struct field_info *Fi;
+
+ /* initialize GIS environment */
+ G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
+
+ /* initialize module */
+ module = G_define_module();
+ module->keywords = _("network, shortest path");
+ module->description =
+ _
+ ("Find the shortest paths from a feature 'to' to every feature 'from' and various information about this realtion are uploaded to the attribute table.");
+
+ /* Define the different options as defined in gis.h */
+ map_in = G_define_standard_option(G_OPT_V_INPUT);
+ map_out = G_define_standard_option(G_OPT_V_OUTPUT);
+
+ fieldf_opt = G_define_standard_option(G_OPT_V_FIELD);
+ fieldf_opt->key = "from_layer";
+ fieldf_opt->description = _("From layer");
+ catf_opt = G_define_standard_option(G_OPT_V_CATS);
+ catf_opt->key = "from_cats";
+ catf_opt->description = _("From category values");
+ wheref_opt = G_define_standard_option(G_OPT_WHERE);
+ wheref_opt->key = "from_where";
+ wheref_opt->description =
+ _("From WHERE conditions of SQL statement without 'where' keyword");
+
+ fieldt_opt = G_define_standard_option(G_OPT_V_FIELD);
+ fieldt_opt->key = "to_layer";
+ fieldt_opt->description = _("To layer");
+ catt_opt = G_define_standard_option(G_OPT_V_CATS);
+ catt_opt->key = "to_cats";
+ catt_opt->description = _("To category values");
+ wheret_opt = G_define_standard_option(G_OPT_WHERE);
+ wheret_opt->key = "to_where";
+ wheret_opt->description =
+ _("To WHERE conditions of SQL statement without 'where' keyword");
+
+ afcol = G_define_option();
+ afcol->key = "afcolumn";
+ afcol->type = TYPE_STRING;
+ afcol->required = NO;
+ afcol->description = _("Arc forward/both direction(s) cost column");
+
+ abcol = G_define_option();
+ abcol->key = "abcolumn";
+ abcol->type = TYPE_STRING;
+ abcol->required = NO;
+ abcol->description = _("Arc backward direction cost column");
+
+ geo_f = G_define_flag();
+ geo_f->key = 'g';
+ geo_f->description =
+ _("Use geodesic calculation for longitude-latitude locations");
+
+
+ /* options and flags parser */
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+ /* TODO: make an option for this */
+ mask_type = GV_POINT | GV_LINE | GV_BOUNDARY;
+
+ Points = Vect_new_line_struct();
+ Cats = Vect_new_cats_struct();
+
+ Vect_check_input_output_name(map_in->answer, map_out->answer,
+ GV_FATAL_EXIT);
+
+ if ((mapset = G_find_vector2(map_in->answer, "")) == NULL)
+ G_fatal_error(_("Vector map <%s> not found"), map_in->answer);
+
+ Vect_set_open_level(2);
+
+ if (1 > Vect_open_old(&In, map_in->answer, mapset))
+ G_fatal_error(_("Unable to open vector map <%s>"),
+ G_fully_qualified_name(map_in->answer, mapset));
+
+ with_z = Vect_is_3d(&In);
+
+ if (0 > Vect_open_new(&Out, map_out->answer, with_z)) {
+ Vect_close(&In);
+ G_fatal_error(_("Unable to create vector map <%s>"), map_out->answer);
+ }
+
+
+ if (geo_f->answer) {
+ geo = 1;
+ if (G_projection() != PROJECTION_LL)
+ G_warning(_("The current projection is not longitude-latitude"));
+ }
+ else
+ geo = 0;
+
+
+ nnodes = Vect_get_num_nodes(&In);
+ nlines = Vect_get_num_lines(&In);
+
+ dst = (int *)G_calloc(nnodes + 1, sizeof(int));
+ prev = (dglInt32_t **) G_calloc(nnodes + 1, sizeof(dglInt32_t *));
+ nodes_to_features = (int *)G_calloc(nnodes + 1, sizeof(int));
+ on_path =
+ (struct line_cats **)G_calloc(nlines + 1, sizeof(struct line_cats *));
+ if (!dst || !prev || !nodes_to_features || !on_path)
+ G_fatal_error(_("Out of memory"));
+
+ for (i = 1; i <= nlines; i++)
+ on_path[i] = Vect_new_cats_struct();
+
+ /*initialise varrays and nodes list appropriatelly */
+ flayer = atoi(fieldf_opt->answer);
+ tlayer = atoi(fieldt_opt->answer);
+ initiliase_varray(&In, flayer, GV_POINT, wheref_opt->answer,
+ catf_opt->answer, &varrayf);
+ initiliase_varray(&In, tlayer, mask_type, wheret_opt->answer,
+ catt_opt->answer, &varrayt);
+
+ nodest = Vect_new_list();
+ neta_varray_to_nodes(&In, varrayt, nodest, nodes_to_features);
+
+ Vect_net_build_graph(&In, mask_type, 1, 0, afcol->answer, abcol->answer,
+ NULL, geo, 0);
+ graph = &(In.graph);
+ neta_distance_from_points(graph, nodest, dst, prev);
+
+ /* Create table */
+ Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
+ Vect_map_add_dblink(&Out, 1, NULL, Fi->table, "cat", Fi->database,
+ Fi->driver);
+
+ 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);
+
+ sprintf(buf,
+ "create table %s ( cat integer, tcat integer, dist double precision)",
+ Fi->table);
+
+ db_set_string(&sql, buf);
+ G_debug(2, db_get_string(&sql));
+
+ if (db_execute_immediate(driver, &sql) != DB_OK) {
+ db_close_database_shutdown_driver(driver);
+ G_fatal_error(_("Unable to create table: '%s'"), db_get_string(&sql));
+ }
+
+ if (db_create_index2(driver, Fi->table, "cat") != DB_OK)
+ G_warning(_("Cannot create index"));
+
+ if (db_grant_on_table
+ (driver, Fi->table, DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK)
+ G_fatal_error(_("Cannot grant privileges on table <%s>"), Fi->table);
+
+ db_begin_transaction(driver);
+
+ Vect_copy_head_data(&In, &Out);
+ Vect_hist_copy(&In, &Out);
+ Vect_hist_command(&Out);
+
+
+ for (i = 1; i <= nlines; i++)
+ if (varrayf->c[i]) {
+ int type = Vect_read_line(&In, Points, Cats, i);
+ int node, tcat, cat;
+ double cost;
+ dglInt32_t *vertex, vertex_id;
+ if (!Vect_cat_get(Cats, flayer, &cat))
+ continue;
+ Vect_get_line_nodes(&In, i, &node, NULL);
+ Vect_write_line(&Out, type, Points, Cats);
+ cost = dst[node] / (double)In.cost_multip;
+ vertex = dglGetNode(graph, node);
+ vertex_id = node;
+ while (prev[vertex_id] != NULL) {
+ Vect_cat_set(on_path
+ [abs(dglEdgeGet_Id(graph, prev[vertex_id]))], 1,
+ cat);
+ vertex = dglEdgeGet_Head(graph, prev[vertex_id]);
+ vertex_id = dglNodeGet_Id(graph, vertex);
+ }
+ Vect_read_line(&In, NULL, Cats, nodes_to_features[vertex_id]);
+ if (!Vect_cat_get(Cats, tlayer, &tcat))
+ continue;
+ sprintf(buf, "insert into %s values (%d, %d, %f)", Fi->table, cat,
+ tcat, cost);
+
+ db_set_string(&sql, buf);
+ G_debug(3, db_get_string(&sql));
+ if (db_execute_immediate(driver, &sql) != DB_OK) {
+ db_close_database_shutdown_driver(driver);
+ G_fatal_error(_("Cannot insert new record: %s"),
+ db_get_string(&sql));
+ };
+ }
+
+ for (i = 1; i <= nlines; i++)
+ if (on_path[i]->n_cats > 0) {
+ int type = Vect_read_line(&In, Points, NULL, i);
+ Vect_write_line(&Out, type, Points, on_path[i]);
+ }
+
+ db_commit_transaction(driver);
+ db_close_database_shutdown_driver(driver);
+
+ Vect_build(&Out);
+
+ Vect_close(&In);
+ Vect_close(&Out);
+
+ for (i = 1; i <= nlines; i++)
+ Vect_destroy_cats_struct(on_path[i]);
+ G_free(on_path);
+ G_free(nodes_to_features);
+ G_free(dst);
+ G_free(prev);
+
+ exit(EXIT_SUCCESS);
+}
Property changes on: grass-addons/vector/net.analyze/v.net.distance/main.c
___________________________________________________________________
Added: svn:executable
+ *
Modified: grass-addons/vector/net.analyze/v.net.flow/main.c
===================================================================
--- grass-addons/vector/net.analyze/v.net.flow/main.c 2009-07-07 11:38:23 UTC (rev 38292)
+++ grass-addons/vector/net.analyze/v.net.flow/main.c 2009-07-07 13:21:27 UTC (rev 38293)
@@ -5,7 +5,7 @@
*
* AUTHOR(S): Daniel Bundala
*
- * PURPOSE: Max flow between two nodes
+ * PURPOSE: Max flow and min cut between two sets of nodes
*
* COPYRIGHT: (C) 2002-2005 by the GRASS Development Team
*
@@ -27,18 +27,23 @@
int main(int argc, char *argv[])
{
- struct Map_info In, Out;
+ struct Map_info In, Out, cut_map;
static struct line_pnts *Points;
struct line_cats *Cats;
char *mapset;
struct GModule *module; /* GRASS module for parsing arguments */
- struct Option *map_in, *map_out;
- struct Option *cat_opt, *field_opt, *where_opt, *abcol, *afcol;
+ struct Option *map_in, *map_out, *cut_out;
+ struct Option *cat_opt, *field_opt, *where_opt, *abcol, *afcol, *source_opt,
+ *sink_opt;
int chcat, with_z;
int layer, mask_type;
VARRAY *varray;
dglGraph_s *graph;
int i, nlines, *flow, total_flow;
+ struct ilist *source_list, *sink_list, *cut;
+ struct cat_list *source_cats, *sink_cats;
+ int find_cut;
+
char buf[2000];
/* Attribute table */
@@ -52,12 +57,19 @@
/* initialize module */
module = G_define_module();
module->keywords = _("network, max flow");
- module->description = _("Computes the maximum flow between two nodes.");
+ module->description =
+ _("Computes the maximum flow between two sets of nodes.");
/* Define the different options as defined in gis.h */
map_in = G_define_standard_option(G_OPT_V_INPUT);
map_out = G_define_standard_option(G_OPT_V_OUTPUT);
+ cut_out = G_define_option();
+ cut_out->type = TYPE_STRING;
+ cut_out->required = YES;
+ cut_out->key = "cut";
+ cut_out->description = _("Name of output map containing a minimum cut");
+
field_opt = G_define_standard_option(G_OPT_V_FIELD);
cat_opt = G_define_standard_option(G_OPT_V_CATS);
where_opt = G_define_standard_option(G_OPT_WHERE);
@@ -74,9 +86,20 @@
abcol->required = NO;
abcol->description = _("Arc backward direction capacity column");
+ source_opt = G_define_standard_option(G_OPT_V_CATS);
+ source_opt->key = "source";
+ source_opt->required = YES;
+ source_opt->description = _("Categories of source(s)");
+
+ sink_opt = G_define_standard_option(G_OPT_V_CATS);
+ sink_opt->key = "sink";
+ sink_opt->required = YES;
+ sink_opt->description = _("Categories of sink(s)");
+
/* options and flags parser */
if (G_parser(argc, argv))
exit(EXIT_FAILURE);
+ find_cut = (cut_out->answer[0]);
/* TODO: make an option for this */
mask_type = GV_LINE | GV_BOUNDARY;
@@ -102,6 +125,12 @@
G_fatal_error(_("Unable to create vector map <%s>"), map_out->answer);
}
+ if (find_cut && 0 > Vect_open_new(&cut_map, cut_out->answer, with_z)) {
+ Vect_close(&In);
+ Vect_close(&Out);
+ G_fatal_error(_("Unable to create vector map <%s>"), cut_out->answer);
+ }
+
/* parse filter option and select appropriate lines */
layer = atoi(field_opt->answer);
if (where_opt->answer) {
@@ -161,6 +190,35 @@
db_begin_transaction(driver);
+ source_cats = Vect_new_cat_list();
+ sink_cats = Vect_new_cat_list();
+ Vect_str_to_cat_list(source_opt->answer, source_cats);
+ Vect_str_to_cat_list(sink_opt->answer, sink_cats);
+ source_list = Vect_new_list();
+ sink_list = Vect_new_list();
+
+ nlines = Vect_get_num_lines(&In);
+ for (i = 1; i <= nlines; i++) {
+ int type, cat;
+ type = Vect_read_line(&In, NULL, Cats, i);
+ if (type != GV_POINT)
+ continue;
+ Vect_cat_get(Cats, layer, &cat);
+ if (Vect_cat_in_cat_list(cat, source_cats))
+ Vect_list_append(source_list, i);
+ if (Vect_cat_in_cat_list(cat, sink_cats))
+ Vect_list_append(sink_list, i);
+ }
+
+ if (source_list->n_values == 0)
+ G_fatal_error(_("No points with categories [%s]"), source_opt->answer);
+
+ if (sink_list->n_values == 0)
+ G_fatal_error(_("No points with categories [%s]"), sink_opt->answer);
+
+ neta_points_to_nodes(&In, source_list);
+ neta_points_to_nodes(&In, sink_list);
+
Vect_copy_head_data(&In, &Out);
Vect_hist_copy(&In, &Out);
@@ -169,14 +227,18 @@
Vect_net_build_graph(&In, mask_type, atoi(field_opt->answer), 0,
afcol->answer, abcol->answer, NULL, 0, 0);
graph = &(In.graph);
- nlines = Vect_get_num_lines(&In);
flow = (int *)G_calloc(nlines + 1, sizeof(int));
if (!flow)
G_fatal_error(_("Out of memory"));
- total_flow = neta_flow(graph, 215, 219, flow);
+ total_flow = neta_flow(graph, source_list, sink_list, flow);
G_debug(3, "Max flow: %d", total_flow);
+ if (find_cut) {
+ cut = Vect_new_list();
+ total_flow = neta_min_cut(graph, source_list, sink_list, flow, cut);
+ G_debug(3, "Min cut: %d", total_flow);
+ }
G_message(_("Writing the output..."));
G_percent_reset();
@@ -190,7 +252,7 @@
if (cat == -1)
continue; /*TODO: warning? */
sprintf(buf, "insert into %s values (%d, %d)", Fi->table, cat,
- flow[i] / 1000);
+ flow[i] / In.cost_multip);
db_set_string(&sql, buf);
G_debug(3, db_get_string(&sql));
@@ -201,12 +263,24 @@
};
}
}
- neta_add_point_on_node(&In, &Out, 215, Cats);
- neta_add_point_on_node(&In, &Out, 219, Cats);
+
+ if (find_cut) {
+ for (i = 0; i < cut->n_values; i++) {
+ int type = Vect_read_line(&In, Points, Cats, cut->value[i]);
+ Vect_write_line(&cut_map, type, Points, Cats);
+ }
+ Vect_destroy_list(cut);
+
+ Vect_build(&cut_map);
+ Vect_close(&cut_map);
+ }
+
db_commit_transaction(driver);
db_close_database_shutdown_driver(driver);
G_free(flow);
+ Vect_destroy_list(source_list);
+ Vect_destroy_list(sink_list);
Vect_build(&Out);
More information about the grass-commit
mailing list