GRASS Programmer's Manual  6.4.2(2012)
flow.c
Go to the documentation of this file.
00001 
00017 #include <stdio.h>
00018 #include <stdlib.h>
00019 #include <grass/gis.h>
00020 #include <grass/Vect.h>
00021 #include <grass/glocale.h>
00022 #include <grass/dgl/graph.h>
00023 #include <grass/neta.h>
00024 
00025 dglInt32_t sign(dglInt32_t x)
00026 {
00027     if (x >= 0)
00028         return 1;
00029     return -1;
00030 }
00031 
00047 int NetA_flow(dglGraph_s * graph, struct ilist *source_list,
00048               struct ilist *sink_list, int *flow)
00049 {
00050     int nnodes, nlines, i;
00051     dglEdgesetTraverser_s et;
00052     dglInt32_t *queue;
00053     dglInt32_t **prev;
00054     char *is_source, *is_sink;
00055     int begin, end, total_flow;
00056 
00057     nnodes = dglGet_NodeCount(graph);
00058     nlines = dglGet_EdgeCount(graph) / 2;       /*each line corresponds to two edges. One in each direction */
00059     queue = (dglInt32_t *) G_calloc(nnodes + 3, sizeof(dglInt32_t));
00060     prev = (dglInt32_t **) G_calloc(nnodes + 3, sizeof(dglInt32_t *));
00061     is_source = (char *)G_calloc(nnodes + 3, sizeof(char));
00062     is_sink = (char *)G_calloc(nnodes + 3, sizeof(char));
00063     if (!queue || !prev || !is_source || !is_sink) {
00064         G_fatal_error(_("Out of memory"));
00065         return -1;
00066     }
00067 
00068     for (i = 0; i < source_list->n_values; i++)
00069         is_source[source_list->value[i]] = 1;
00070     for (i = 0; i < sink_list->n_values; i++)
00071         is_sink[sink_list->value[i]] = 1;
00072 
00073     for (i = 0; i <= nlines; i++)
00074         flow[i] = 0;
00075 
00076     total_flow = 0;
00077     while (1) {
00078         dglInt32_t node, edge_id, min_residue;
00079         int found = -1;
00080 
00081         begin = end = 0;
00082         for (i = 0; i < source_list->n_values; i++)
00083             queue[end++] = source_list->value[i];
00084 
00085         for (i = 1; i <= nnodes; i++) {
00086             prev[i] = NULL;
00087         }
00088         while (begin != end && found == -1) {
00089             dglInt32_t vertex = queue[begin++];
00090             dglInt32_t *edge, *node = dglGetNode(graph, vertex);
00091 
00092             dglEdgeset_T_Initialize(&et, graph,
00093                                     dglNodeGet_OutEdgeset(graph, node));
00094             for (edge = dglEdgeset_T_First(&et); edge;
00095                  edge = dglEdgeset_T_Next(&et)) {
00096                 dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
00097                 dglInt32_t id = dglEdgeGet_Id(graph, edge);
00098                 dglInt32_t to =
00099                     dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
00100                 if (!is_source[to] && prev[to] == NULL &&
00101                     cap > sign(id) * flow[abs(id)]) {
00102                     prev[to] = edge;
00103                     if (is_sink[to]) {
00104                         found = to;
00105                         break;
00106                     }
00107                     queue[end++] = to;
00108                 }
00109             }
00110             dglEdgeset_T_Release(&et);
00111         }
00112         if (found == -1)
00113             break;              /*no augmenting path */
00114         /*find minimum residual capacity along the augmenting path */
00115         node = found;
00116         edge_id = dglEdgeGet_Id(graph, prev[node]);
00117         min_residue =
00118             dglEdgeGet_Cost(graph,
00119                             prev[node]) - sign(edge_id) * flow[abs(edge_id)];
00120         while (!is_source[node]) {
00121             dglInt32_t residue;
00122 
00123             edge_id = dglEdgeGet_Id(graph, prev[node]);
00124             residue =
00125                 dglEdgeGet_Cost(graph,
00126                                 prev[node]) -
00127                 sign(edge_id) * flow[abs(edge_id)];
00128             if (residue < min_residue)
00129                 min_residue = residue;
00130             node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
00131         }
00132         total_flow += min_residue;
00133         /*update flow along the augmenting path */
00134         node = found;
00135         while (!is_source[node]) {
00136             edge_id = dglEdgeGet_Id(graph, prev[node]);
00137             flow[abs(edge_id)] += sign(edge_id) * min_residue;
00138             node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
00139         }
00140     }
00141 
00142     G_free(queue);
00143     G_free(prev);
00144     G_free(is_source);
00145     G_free(is_sink);
00146 
00147     return total_flow;
00148 }
00149 
00165 int NetA_min_cut(dglGraph_s * graph, struct ilist *source_list,
00166                  struct ilist *sink_list, int *flow, struct ilist *cut)
00167 {
00168     int nnodes, i;
00169     dglEdgesetTraverser_s et;
00170     dglInt32_t *queue;
00171     char *visited;
00172     int begin, end, total_flow;
00173 
00174     nnodes = dglGet_NodeCount(graph);
00175     queue = (dglInt32_t *) G_calloc(nnodes + 3, sizeof(dglInt32_t));
00176     visited = (char *)G_calloc(nnodes + 3, sizeof(char));
00177     if (!queue || !visited) {
00178         G_fatal_error(_("Out of memory"));
00179         return -1;
00180     }
00181 
00182     total_flow = begin = end = 0;
00183 
00184     for (i = 1; i <= nnodes; i++)
00185         visited[i] = 0;
00186 
00187     for (i = 0; i < source_list->n_values; i++) {
00188         queue[end++] = source_list->value[i];
00189         visited[source_list->value[i]] = 1;
00190     }
00191 
00192     /* find vertices reachable from source(s) using only non-saturated edges */
00193     while (begin != end) {
00194         dglInt32_t vertex = queue[begin++];
00195         dglInt32_t *edge, *node = dglGetNode(graph, vertex);
00196 
00197         dglEdgeset_T_Initialize(&et, graph,
00198                                 dglNodeGet_OutEdgeset(graph, node));
00199         for (edge = dglEdgeset_T_First(&et); edge;
00200              edge = dglEdgeset_T_Next(&et)) {
00201             dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
00202             dglInt32_t id = dglEdgeGet_Id(graph, edge);
00203             dglInt32_t to =
00204                 dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
00205             if (!visited[to] && cap > sign(id) * flow[abs(id)]) {
00206                 visited[to] = 1;
00207                 queue[end++] = to;
00208             }
00209         }
00210         dglEdgeset_T_Release(&et);
00211     }
00212     /*saturated edges from reachable vertices to non-reachable ones form a minimum cost */
00213     Vect_reset_list(cut);
00214     for (i = 1; i <= nnodes; i++) {
00215         if (!visited[i])
00216             continue;
00217         dglInt32_t *node, *edgeset, *edge;
00218 
00219         node = dglGetNode(graph, i);
00220         edgeset = dglNodeGet_OutEdgeset(graph, node);
00221         dglEdgeset_T_Initialize(&et, graph, edgeset);
00222         for (edge = dglEdgeset_T_First(&et); edge;
00223              edge = dglEdgeset_T_Next(&et)) {
00224             dglInt32_t to, edge_id;
00225 
00226             to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
00227             edge_id = abs(dglEdgeGet_Id(graph, edge));
00228             if (!visited[to] && flow[edge_id] != 0) {
00229                 Vect_list_append(cut, edge_id);
00230                 total_flow += abs(flow[abs(edge_id)]);
00231             }
00232         }
00233         dglEdgeset_T_Release(&et);
00234     }
00235 
00236     G_free(visited);
00237     G_free(queue);
00238     return total_flow;
00239 }
00240 
00258 int NetA_split_vertices(dglGraph_s * in, dglGraph_s * out, int *node_costs)
00259 {
00260     dglInt32_t opaqueset[16] =
00261         { 360000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
00262     dglNodeTraverser_s nt;
00263     dglInt32_t nnodes, edge_cnt;
00264     dglInt32_t *cur_node;
00265 
00266     nnodes = dglGet_NodeCount(in);
00267     dglInitialize(out, (dglByte_t) 1, (dglInt32_t) 0, (dglInt32_t) 0,
00268                   opaqueset);
00269     dglNode_T_Initialize(&nt, in);
00270     edge_cnt = 0;
00271     dglInt32_t max_node_cost = 0;
00272 
00273     for (cur_node = dglNode_T_First(&nt); cur_node;
00274          cur_node = dglNode_T_Next(&nt)) {
00275         dglInt32_t v = dglNodeGet_Id(in, cur_node);
00276 
00277         edge_cnt++;
00278         dglInt32_t cost = 1;
00279 
00280         if (node_costs)
00281             cost = node_costs[v];
00282         if (cost > max_node_cost)
00283             max_node_cost = cost;
00284         dglAddEdge(out, 2 * v - 1, 2 * v, cost, edge_cnt);
00285         dglAddEdge(out, 2 * v, 2 * v - 1, (dglInt32_t) 0, -edge_cnt);
00286     }
00287     dglNode_T_Release(&nt);
00288     dglNode_T_Initialize(&nt, in);
00289     for (cur_node = dglNode_T_First(&nt); cur_node;
00290          cur_node = dglNode_T_Next(&nt)) {
00291         dglEdgesetTraverser_s et;
00292         dglInt32_t *edge;
00293         dglInt32_t v = dglNodeGet_Id(in, cur_node);
00294 
00295         dglEdgeset_T_Initialize(&et, in, dglNodeGet_OutEdgeset(in, cur_node));
00296         for (edge = dglEdgeset_T_First(&et); edge;
00297              edge = dglEdgeset_T_Next(&et)) {
00298             dglInt32_t to;
00299 
00300             to = dglNodeGet_Id(in, dglEdgeGet_Tail(in, edge));
00301             edge_cnt++;
00302             dglAddEdge(out, 2 * v, 2 * to - 1, max_node_cost + 1, edge_cnt);
00303             dglAddEdge(out, 2 * to - 1, 2 * v, (dglInt32_t) 0, -edge_cnt);
00304         }
00305         dglEdgeset_T_Release(&et);
00306     }
00307     dglNode_T_Release(&nt);
00308     if (dglFlatten(out) < 0)
00309         G_fatal_error(_("GngFlatten error"));
00310     return edge_cnt;
00311 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines