GRASS 8 Programmer's Manual 8.6.0dev(2026)-f6f2c534ea
Loading...
Searching...
No Matches
flow.c
Go to the documentation of this file.
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
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 */
48 struct ilist *sink_list, int *flow)
49{
50 int nnodes, nlines, i;
53 dglInt32_t **prev;
54 char *is_source, *is_sink;
55 int begin, end, total_flow;
58
60 nlines = dglGet_EdgeCount(graph) /
61 2; /*each line corresponds to two edges. One in each direction */
62 queue = (dglInt32_t *)G_calloc(nnodes + 3, sizeof(dglInt32_t));
63 prev = (dglInt32_t **)G_calloc(nnodes + 3, sizeof(dglInt32_t *));
64 is_source = (char *)G_calloc(nnodes + 3, sizeof(char));
65 is_sink = (char *)G_calloc(nnodes + 3, sizeof(char));
66 if (!queue || !prev || !is_source || !is_sink) {
67 G_fatal_error(_("Out of memory"));
68 return -1;
69 }
70
71 for (i = 0; i < source_list->n_values; i++)
72 is_source[source_list->value[i]] = 1;
73 for (i = 0; i < sink_list->n_values; i++)
74 is_sink[sink_list->value[i]] = 1;
75
76 for (i = 0; i <= nlines; i++)
77 flow[i] = 0;
78
79 ncost = 0;
81
82 total_flow = 0;
83 while (1) {
85 int found = -1;
86
87 begin = end = 0;
88 for (i = 0; i < source_list->n_values; i++)
89 queue[end++] = source_list->value[i];
90
91 for (i = 1; i <= nnodes; i++) {
92 prev[i] = NULL;
93 }
94 while (begin != end && found == -1) {
96 dglInt32_t *edge, *node = dglGetNode(graph, vertex);
97
100 for (edge = dglEdgeset_T_First(&et); edge;
101 edge = dglEdgeset_T_Next(&et)) {
103 dglInt32_t id = dglEdgeGet_Id(graph, edge);
104 dglInt32_t to =
106 if (!is_source[to] && prev[to] == NULL &&
107 cap > sign(id) * flow[labs(id)]) {
108 prev[to] = edge;
109 if (is_sink[to]) {
110 found = to;
111 break;
112 }
113 /* do not go through closed nodes */
114 if (have_node_costs) {
115 memcpy(&ncost,
117 dglEdgeGet_Tail(graph, edge)),
118 sizeof(ncost));
119 }
120 if (ncost >= 0)
121 queue[end++] = to;
122 }
123 }
125 }
126 if (found == -1)
127 break; /*no augmenting path */
128 /*find minimum residual capacity along the augmenting path */
129 node = found;
130 edge_id = dglEdgeGet_Id(graph, prev[node]);
131 min_residue = dglEdgeGet_Cost(graph, prev[node]) -
133 while (!is_source[node]) {
135
136 edge_id = dglEdgeGet_Id(graph, prev[node]);
137 residue = dglEdgeGet_Cost(graph, prev[node]) -
139 if (residue < min_residue)
141 node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
142 }
144 /*update flow along the augmenting path */
145 node = found;
146 while (!is_source[node]) {
147 edge_id = dglEdgeGet_Id(graph, prev[node]);
149 node = dglNodeGet_Id(graph, dglEdgeGet_Head(graph, prev[node]));
150 }
151 }
152
153 G_free(queue);
154 G_free(prev);
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 (unused)
171 \param flow
172 \param[out] cut list of edges (cut)
173
174 \return number of edges
175 \return -1 on failure
176 */
178 struct ilist *sink_list UNUSED, int *flow, struct ilist *cut)
179{
180 int nnodes, i;
183 char *visited;
184 int begin, end, total_flow;
185
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) {
207 dglInt32_t *edge, *node = dglGetNode(graph, vertex);
208
210 for (edge = dglEdgeset_T_First(&et); edge;
211 edge = dglEdgeset_T_Next(&et)) {
213 dglInt32_t id = dglEdgeGet_Id(graph, edge);
215 if (!visited[to] && cap > sign(id) * flow[labs(id)]) {
216 visited[to] = 1;
217 queue[end++] = to;
218 }
219 }
221 }
222 /*saturated edges from reachable vertices to non-reachable ones form a
223 * minimum cost */
225 for (i = 1; i <= nnodes; i++) {
226 if (!visited[i])
227 continue;
228 dglInt32_t *node, *edgeset, *edge;
229
230 node = dglGetNode(graph, i);
233 for (edge = dglEdgeset_T_First(&et); edge;
234 edge = dglEdgeset_T_Next(&et)) {
236
239 if (!visited[to] && flow[edge_id] != 0) {
242 }
243 }
245 }
246
248 G_free(queue);
249 return total_flow;
250}
251
252/*!
253 \brief Splits each vertex of in graph into two vertices
254
255 The method splits each vertex of in graph into two vertices: in
256 vertex and out vertex. Also, it adds an edge from an in vertex to
257 the corresponding out vertex (capacity=2) and it adds an edge from
258 out vertex to in vertex for each edge present in the in graph
259 (forward capacity=1, backward capacity=0). If the id of a vertex is
260 v then id of in vertex is 2*v-1 and of out vertex 2*v.
261
262 \param in from graph
263 \param out to graph
264 \param node_costs list of node costs
265
266 \return number of undirected edges in the graph
267 \return -1 on failure
268 */
269int NetA_split_vertices(dglGraph_s *in, dglGraph_s *out, int *node_costs)
270{
271 dglInt32_t opaqueset[16] = {360000, 0, 0, 0, 0, 0, 0, 0,
272 0, 0, 0, 0, 0, 0, 0, 0};
276
279 edge_cnt = 0;
281
285
286 edge_cnt++;
287 dglInt32_t cost = 1;
288
289 if (node_costs)
290 cost = node_costs[v];
291 /* skip closed nodes */
292 if (cost < 0)
293 continue;
294 if (cost > max_node_cost)
295 max_node_cost = cost;
296 dglAddEdge(out, 2 * v - 1, 2 * v, cost, edge_cnt);
297 dglAddEdge(out, 2 * v, 2 * v - 1, (dglInt32_t)0, -edge_cnt);
298 }
304 dglInt32_t *edge;
306 dglInt32_t cost = 1;
307
308 if (node_costs)
309 cost = node_costs[v];
310 /* skip closed nodes */
311 if (cost < 0)
312 continue;
313
315 for (edge = dglEdgeset_T_First(&et); edge;
316 edge = dglEdgeset_T_Next(&et)) {
317 dglInt32_t to;
318
319 to = dglNodeGet_Id(in, dglEdgeGet_Tail(in, edge));
320 edge_cnt++;
321 dglAddEdge(out, 2 * v, 2 * to - 1, max_node_cost + 1, edge_cnt);
322 dglAddEdge(out, 2 * to - 1, 2 * v, (dglInt32_t)0, -edge_cnt);
323 }
325 }
327 if (dglFlatten(out) < 0)
328 G_fatal_error(_("GngFlatten error"));
329 return edge_cnt;
330}
#define NULL
Definition ccmath.h:32
Definition queue.h:43
void G_free(void *)
Free allocated memory.
Definition gis/alloc.c:147
#define G_calloc(m, n)
Definition defs/gis.h:140
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
int Vect_list_append(struct ilist *, int)
Append new item to the end of list if not yet present.
int Vect_reset_list(struct ilist *)
Reset ilist structure.
int NetA_min_cut(dglGraph_s *graph, struct ilist *source_list, struct ilist *sink_list, int *flow, struct ilist *cut)
Calculates minimum cut between source(s) and sink(s).
Definition flow.c:177
dglInt32_t sign(dglInt32_t x)
Definition flow.c:25
int NetA_flow(dglGraph_s *graph, struct ilist *source_list, struct ilist *sink_list, int *flow)
Get max flow from source to sink.
Definition flow.c:47
int NetA_split_vertices(dglGraph_s *in, dglGraph_s *out, int *node_costs)
Splits each vertex of in graph into two vertices.
Definition flow.c:269
#define UNUSED
A macro for an attribute, if attached to a variable, indicating that the variable is not used.
Definition gis.h:46
#define _(str)
Definition glocale.h:10
List of integers.
Definition gis.h:715
Definition path.h:10
unsigned char dglByte_t
Definition type.h:35
long dglInt32_t
Definition type.h:36
dglInt32_t * dglNode_T_Next(dglNodeTraverser_s *pT)
dglInt32_t * dglNode_T_First(dglNodeTraverser_s *pT)
dglInt32_t * dglGetNode(dglGraph_s *pGraph, dglInt32_t nNodeId)
int dglGet_NodeAttrSize(dglGraph_s *pgraph)
dglInt32_t * dglNodeGet_OutEdgeset(dglGraph_s *pGraph, dglInt32_t *pnNode)
int dglAddEdge(dglGraph_s *pGraph, dglInt32_t nHead, dglInt32_t nTail, dglInt32_t nCost, dglInt32_t nEdge)
int dglEdgeset_T_Initialize(dglEdgesetTraverser_s *pT, dglGraph_s *pGraph, dglInt32_t *pnEdgeset)
int dglNode_T_Initialize(dglNodeTraverser_s *pT, dglGraph_s *pGraph)
dglInt32_t dglEdgeGet_Id(dglGraph_s *pGraph, dglInt32_t *pnEdge)
dglInt32_t * dglEdgeset_T_First(dglEdgesetTraverser_s *pT)
int dglInitialize(dglGraph_s *pGraph, dglByte_t Version, dglInt32_t NodeAttrSize, dglInt32_t EdgeAttrSize, dglInt32_t *pOpaqueSet)
int dglGet_EdgeCount(dglGraph_s *pgraph)
dglInt32_t * dglEdgeset_T_Next(dglEdgesetTraverser_s *pT)
dglInt32_t * dglNodeGet_Attr(dglGraph_s *pGraph, dglInt32_t *pnNode)
int dglGet_NodeCount(dglGraph_s *pgraph)
int dglFlatten(dglGraph_s *pGraph)
dglInt32_t dglEdgeGet_Cost(dglGraph_s *pGraph, dglInt32_t *pnEdge)
void dglNode_T_Release(dglNodeTraverser_s *pT)
void dglEdgeset_T_Release(dglEdgesetTraverser_s *pT)
dglInt32_t * dglEdgeGet_Head(dglGraph_s *pGraph, dglInt32_t *pnEdge)
dglInt32_t * dglEdgeGet_Tail(dglGraph_s *pGraph, dglInt32_t *pnEdge)
dglInt32_t dglNodeGet_Id(dglGraph_s *pGraph, dglInt32_t *pnNode)
#define x