1 /*!
2 \file vector/neta/flow.c
3
4 \brief Network Analysis library - flow in graph
5
6 Computes the length of the shortest path between all pairs of nodes
7 in the network.
8
9 (C) 2009-2010 by Daniel Bundala, and the GRASS Development Team
10
11 This program is free software under the GNU General Public License
12 (>=v2). Read the file COPYING that comes with GRASS for details.
13
14 \author Daniel Bundala (Google Summer of Code 2009)
15 */
16
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <grass/gis.h>
20 #include <grass/vector.h>
21 #include <grass/glocale.h>
22 #include <grass/dgl/graph.h>
23 #include <grass/neta.h>
24
sign(dglInt32_t x)25 dglInt32_t sign(dglInt32_t x)
26 {
27 if (x >= 0)
28 return 1;
29 return -1;
30 }
31
32 /*!
33 \brief Get max flow from source to sink.
34
35 Array flow stores flow for each edge. Negative flow corresponds to a
36 flow in opposite direction. The function assumes that the edge costs
37 correspond to edge capacities.
38
39 \param graph input graph
40 \param source_list list of sources
41 \param sink_list list of sinks
42 \param[out] flow max flows
43
44 \return number of flows
45 \return -1 on failure
46 */
NetA_flow(dglGraph_s * graph,struct ilist * source_list,struct ilist * sink_list,int * flow)47 int NetA_flow(dglGraph_s * graph, struct ilist *source_list,
48 struct ilist *sink_list, int *flow)
49 {
50 int nnodes, nlines, i;
51 dglEdgesetTraverser_s et;
52 dglInt32_t *queue;
53 dglInt32_t **prev;
54 char *is_source, *is_sink;
55 int begin, end, total_flow;
56 int have_node_costs;
57 dglInt32_t ncost;
58
59 nnodes = dglGet_NodeCount(graph);
60 nlines = dglGet_EdgeCount(graph) / 2; /*each line corresponds to two edges. One in each direction */
61 queue = (dglInt32_t *) G_calloc(nnodes + 3, sizeof(dglInt32_t));
62 prev = (dglInt32_t **) G_calloc(nnodes + 3, sizeof(dglInt32_t *));
63 is_source = (char *)G_calloc(nnodes + 3, sizeof(char));
64 is_sink = (char *)G_calloc(nnodes + 3, sizeof(char));
65 if (!queue || !prev || !is_source || !is_sink) {
66 G_fatal_error(_("Out of memory"));
67 return -1;
68 }
69
70 for (i = 0; i < source_list->n_values; i++)
71 is_source[source_list->value[i]] = 1;
72 for (i = 0; i < sink_list->n_values; i++)
73 is_sink[sink_list->value[i]] = 1;
74
75 for (i = 0; i <= nlines; i++)
76 flow[i] = 0;
77
78 ncost = 0;
79 have_node_costs = dglGet_NodeAttrSize(graph);
80
81 total_flow = 0;
82 while (1) {
83 dglInt32_t node, edge_id, min_residue;
84 int found = -1;
85
86 begin = end = 0;
87 for (i = 0; i < source_list->n_values; i++)
88 queue[end++] = source_list->value[i];
89
90 for (i = 1; i <= nnodes; i++) {
91 prev[i] = NULL;
92 }
93 while (begin != end && found == -1) {
94 dglInt32_t vertex = queue[begin++];
95 dglInt32_t *edge, *node = dglGetNode(graph, vertex);
96
97 dglEdgeset_T_Initialize(&et, graph,
98 dglNodeGet_OutEdgeset(graph, node));
99 for (edge = dglEdgeset_T_First(&et); edge;
100 edge = dglEdgeset_T_Next(&et)) {
101 dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
102 dglInt32_t id = dglEdgeGet_Id(graph, edge);
103 dglInt32_t to =
104 dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
105 if (!is_source[to] && prev[to] == NULL &&
106 cap > sign(id) * flow[abs(id)]) {
107 prev[to] = edge;
108 if (is_sink[to]) {
109 found = to;
110 break;
111 }
112 /* do not go through closed nodes */
113 if (have_node_costs) {
114 memcpy(&ncost, dglNodeGet_Attr(graph, dglEdgeGet_Tail(graph, edge)),
115 sizeof(ncost));
116 }
117 if (ncost >= 0)
118 queue[end++] = to;
119 }
120 }
121 dglEdgeset_T_Release(&et);
122 }
123 if (found == -1)
124 break; /*no augmenting path */
125 /*find minimum residual capacity along the augmenting path */
126 node = found;
127 edge_id = dglEdgeGet_Id(graph, prev[node]);
128 min_residue =
129 dglEdgeGet_Cost(graph,
130 prev[node]) - sign(edge_id) * flow[abs(edge_id)];
131 while (!is_source[node]) {
132 dglInt32_t residue;
133
134 edge_id = dglEdgeGet_Id(graph, prev[node]);
135 residue =
136 dglEdgeGet_Cost(graph,
137 prev[node]) -
138 sign(edge_id) * flow[abs(edge_id)];
139 if (residue < min_residue)
140 min_residue = residue;
141 node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
142 }
143 total_flow += min_residue;
144 /*update flow along the augmenting path */
145 node = found;
146 while (!is_source[node]) {
147 edge_id = dglEdgeGet_Id(graph, prev[node]);
148 flow[abs(edge_id)] += sign(edge_id) * min_residue;
149 node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
150 }
151 }
152
153 G_free(queue);
154 G_free(prev);
155 G_free(is_source);
156 G_free(is_sink);
157
158 return total_flow;
159 }
160
161 /*!
162 \brief Calculates minimum cut between source(s) and sink(s).
163
164 Flow is the array produced by NetA_flow() method when called with
165 source_list and sink_list as the input. The output of this and
166 NetA_flow() method should be the same.
167
168 \param graph input graph
169 \param source_list list of sources
170 \param sink_list list of sinks
171 \param flow
172 \param[out] cut list of edges (cut)
173
174 \return number of edges
175 \return -1 on failure
176 */
NetA_min_cut(dglGraph_s * graph,struct ilist * source_list,struct ilist * sink_list,int * flow,struct ilist * cut)177 int NetA_min_cut(dglGraph_s * graph, struct ilist *source_list,
178 struct ilist *sink_list, int *flow, struct ilist *cut)
179 {
180 int nnodes, i;
181 dglEdgesetTraverser_s et;
182 dglInt32_t *queue;
183 char *visited;
184 int begin, end, total_flow;
185
186 nnodes = dglGet_NodeCount(graph);
187 queue = (dglInt32_t *) G_calloc(nnodes + 3, sizeof(dglInt32_t));
188 visited = (char *)G_calloc(nnodes + 3, sizeof(char));
189 if (!queue || !visited) {
190 G_fatal_error(_("Out of memory"));
191 return -1;
192 }
193
194 total_flow = begin = end = 0;
195
196 for (i = 1; i <= nnodes; i++)
197 visited[i] = 0;
198
199 for (i = 0; i < source_list->n_values; i++) {
200 queue[end++] = source_list->value[i];
201 visited[source_list->value[i]] = 1;
202 }
203
204 /* find vertices reachable from source(s) using only non-saturated edges */
205 while (begin != end) {
206 dglInt32_t vertex = queue[begin++];
207 dglInt32_t *edge, *node = dglGetNode(graph, vertex);
208
209 dglEdgeset_T_Initialize(&et, graph,
210 dglNodeGet_OutEdgeset(graph, node));
211 for (edge = dglEdgeset_T_First(&et); edge;
212 edge = dglEdgeset_T_Next(&et)) {
213 dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
214 dglInt32_t id = dglEdgeGet_Id(graph, edge);
215 dglInt32_t to =
216 dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
217 if (!visited[to] && cap > sign(id) * flow[abs(id)]) {
218 visited[to] = 1;
219 queue[end++] = to;
220 }
221 }
222 dglEdgeset_T_Release(&et);
223 }
224 /*saturated edges from reachable vertices to non-reachable ones form a minimum cost */
225 Vect_reset_list(cut);
226 for (i = 1; i <= nnodes; i++) {
227 if (!visited[i])
228 continue;
229 dglInt32_t *node, *edgeset, *edge;
230
231 node = dglGetNode(graph, i);
232 edgeset = dglNodeGet_OutEdgeset(graph, node);
233 dglEdgeset_T_Initialize(&et, graph, edgeset);
234 for (edge = dglEdgeset_T_First(&et); edge;
235 edge = dglEdgeset_T_Next(&et)) {
236 dglInt32_t to, edge_id;
237
238 to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
239 edge_id = abs(dglEdgeGet_Id(graph, edge));
240 if (!visited[to] && flow[edge_id] != 0) {
241 Vect_list_append(cut, edge_id);
242 total_flow += abs(flow[abs(edge_id)]);
243 }
244 }
245 dglEdgeset_T_Release(&et);
246 }
247
248 G_free(visited);
249 G_free(queue);
250 return total_flow;
251 }
252
253 /*!
254 \brief Splits each vertex of in graph into two vertices
255
256 The method splits each vertex of in graph into two vertices: in
257 vertex and out vertex. Also, it adds an edge from an in vertex to
258 the corresponding out vertex (capacity=2) and it adds an edge from
259 out vertex to in vertex for each edge present in the in graph
260 (forward capacity=1, backward capacity=0). If the id of a vertex is
261 v then id of in vertex is 2*v-1 and of out vertex 2*v.
262
263 \param in from graph
264 \param out to graph
265 \param node_costs list of node costs
266
267 \return number of undirected edges in the graph
268 \return -1 on failure
269 */
NetA_split_vertices(dglGraph_s * in,dglGraph_s * out,int * node_costs)270 int NetA_split_vertices(dglGraph_s * in, dglGraph_s * out, int *node_costs)
271 {
272 dglInt32_t opaqueset[16] =
273 { 360000, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
274 dglNodeTraverser_s nt;
275 dglInt32_t nnodes, edge_cnt;
276 dglInt32_t *cur_node;
277
278 nnodes = dglGet_NodeCount(in);
279 dglInitialize(out, (dglByte_t) 1, (dglInt32_t) 0, (dglInt32_t) 0,
280 opaqueset);
281 dglNode_T_Initialize(&nt, in);
282 edge_cnt = 0;
283 dglInt32_t max_node_cost = 0;
284
285 for (cur_node = dglNode_T_First(&nt); cur_node;
286 cur_node = dglNode_T_Next(&nt)) {
287 dglInt32_t v = dglNodeGet_Id(in, cur_node);
288
289 edge_cnt++;
290 dglInt32_t cost = 1;
291
292 if (node_costs)
293 cost = node_costs[v];
294 /* skip closed nodes */
295 if (cost < 0)
296 continue;
297 if (cost > max_node_cost)
298 max_node_cost = cost;
299 dglAddEdge(out, 2 * v - 1, 2 * v, cost, edge_cnt);
300 dglAddEdge(out, 2 * v, 2 * v - 1, (dglInt32_t) 0, -edge_cnt);
301 }
302 dglNode_T_Release(&nt);
303 dglNode_T_Initialize(&nt, in);
304 for (cur_node = dglNode_T_First(&nt); cur_node;
305 cur_node = dglNode_T_Next(&nt)) {
306 dglEdgesetTraverser_s et;
307 dglInt32_t *edge;
308 dglInt32_t v = dglNodeGet_Id(in, cur_node);
309 dglInt32_t cost = 1;
310
311 if (node_costs)
312 cost = node_costs[v];
313 /* skip closed nodes */
314 if (cost < 0)
315 continue;
316
317 dglEdgeset_T_Initialize(&et, in, dglNodeGet_OutEdgeset(in, cur_node));
318 for (edge = dglEdgeset_T_First(&et); edge;
319 edge = dglEdgeset_T_Next(&et)) {
320 dglInt32_t to;
321
322 to = dglNodeGet_Id(in, dglEdgeGet_Tail(in, edge));
323 edge_cnt++;
324 dglAddEdge(out, 2 * v, 2 * to - 1, max_node_cost + 1, edge_cnt);
325 dglAddEdge(out, 2 * to - 1, 2 * v, (dglInt32_t) 0, -edge_cnt);
326 }
327 dglEdgeset_T_Release(&et);
328 }
329 dglNode_T_Release(&nt);
330 if (dglFlatten(out) < 0)
331 G_fatal_error(_("GngFlatten error"));
332 return edge_cnt;
333 }
334