GRASS GIS 8 Programmer's Manual  8.5.0dev(2024)-3f124c5c78
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  */
47 int NetA_flow(dglGraph_s *graph, struct ilist *source_list,
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;
56  int have_node_costs;
57  dglInt32_t ncost;
58 
59  nnodes = dglGet_NodeCount(graph);
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;
80  have_node_costs = dglGet_NodeAttrSize(graph);
81 
82  total_flow = 0;
83  while (1) {
84  dglInt32_t node, edge_id, min_residue;
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) {
95  dglInt32_t vertex = queue[begin++];
96  dglInt32_t *edge, *node = dglGetNode(graph, vertex);
97 
98  dglEdgeset_T_Initialize(&et, graph,
99  dglNodeGet_OutEdgeset(graph, node));
100  for (edge = dglEdgeset_T_First(&et); edge;
101  edge = dglEdgeset_T_Next(&et)) {
102  dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
103  dglInt32_t id = dglEdgeGet_Id(graph, edge);
104  dglInt32_t to =
105  dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
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,
116  dglNodeGet_Attr(graph,
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]) -
132  sign(edge_id) * flow[labs(edge_id)];
133  while (!is_source[node]) {
134  dglInt32_t residue;
135 
136  edge_id = dglEdgeGet_Id(graph, prev[node]);
137  residue = dglEdgeGet_Cost(graph, prev[node]) -
138  sign(edge_id) * flow[labs(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[labs(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 (unused)
171  \param flow
172  \param[out] cut list of edges (cut)
173 
174  \return number of edges
175  \return -1 on failure
176  */
177 int NetA_min_cut(dglGraph_s *graph, struct ilist *source_list,
178  struct ilist *sink_list UNUSED, int *flow, struct ilist *cut)
179 {
180  int nnodes, i;
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, dglNodeGet_OutEdgeset(graph, node));
210  for (edge = dglEdgeset_T_First(&et); edge;
211  edge = dglEdgeset_T_Next(&et)) {
212  dglInt32_t cap = dglEdgeGet_Cost(graph, edge);
213  dglInt32_t id = dglEdgeGet_Id(graph, edge);
214  dglInt32_t to = dglNodeGet_Id(graph, dglEdgeGet_Tail(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 */
224  Vect_reset_list(cut);
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);
231  edgeset = dglNodeGet_OutEdgeset(graph, node);
232  dglEdgeset_T_Initialize(&et, graph, edgeset);
233  for (edge = dglEdgeset_T_First(&et); edge;
234  edge = dglEdgeset_T_Next(&et)) {
235  dglInt32_t to, edge_id;
236 
237  to = dglNodeGet_Id(graph, dglEdgeGet_Tail(graph, edge));
238  edge_id = labs(dglEdgeGet_Id(graph, edge));
239  if (!visited[to] && flow[edge_id] != 0) {
240  Vect_list_append(cut, edge_id);
241  total_flow += abs(flow[labs(edge_id)]);
242  }
243  }
245  }
246 
247  G_free(visited);
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  */
269 int 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};
274  dglInt32_t edge_cnt;
275  dglInt32_t *cur_node;
276 
277  dglInitialize(out, (dglByte_t)1, (dglInt32_t)0, (dglInt32_t)0, opaqueset);
278  dglNode_T_Initialize(&nt, in);
279  edge_cnt = 0;
280  dglInt32_t max_node_cost = 0;
281 
282  for (cur_node = dglNode_T_First(&nt); cur_node;
283  cur_node = dglNode_T_Next(&nt)) {
284  dglInt32_t v = dglNodeGet_Id(in, cur_node);
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  }
299  dglNode_T_Release(&nt);
300  dglNode_T_Initialize(&nt, in);
301  for (cur_node = dglNode_T_First(&nt); cur_node;
302  cur_node = dglNode_T_Next(&nt)) {
304  dglInt32_t *edge;
305  dglInt32_t v = dglNodeGet_Id(in, cur_node);
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 
314  dglEdgeset_T_Initialize(&et, in, dglNodeGet_OutEdgeset(in, cur_node));
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  }
326  dglNode_T_Release(&nt);
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:150
#define G_calloc(m, n)
Definition: defs/gis.h:95
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 UNUSED, 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:47
#define _(str)
Definition: glocale.h:10
List of integers.
Definition: gis.h:708
int n_values
Number of values in the list.
Definition: gis.h:716
int * value
Array of values.
Definition: gis.h:712
Definition: path.h:10
unsigned char dglByte_t
Definition: type.h:35
long dglInt32_t
Definition: type.h:36
int dglGet_NodeAttrSize(dglGraph_s *pgraph)
dglInt32_t * dglNode_T_First(dglNodeTraverser_s *pT)
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_Next(dglEdgesetTraverser_s *pT)
dglInt32_t * dglEdgeGet_Tail(dglGraph_s *pGraph, dglInt32_t *pnEdge)
dglInt32_t * dglEdgeGet_Head(dglGraph_s *pGraph, dglInt32_t *pnEdge)
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 * dglNodeGet_Attr(dglGraph_s *pGraph, dglInt32_t *pnNode)
dglInt32_t * dglNode_T_Next(dglNodeTraverser_s *pT)
int dglGet_NodeCount(dglGraph_s *pgraph)
void dglEdgeset_T_Release(dglEdgesetTraverser_s *pT UNUSED)
int dglFlatten(dglGraph_s *pGraph)
dglInt32_t dglEdgeGet_Cost(dglGraph_s *pGraph, dglInt32_t *pnEdge)
dglInt32_t * dglEdgeset_T_First(dglEdgesetTraverser_s *pT)
void dglNode_T_Release(dglNodeTraverser_s *pT)
dglInt32_t dglNodeGet_Id(dglGraph_s *pGraph, dglInt32_t *pnNode)
dglInt32_t * dglGetNode(dglGraph_s *pGraph, dglInt32_t nNodeId)
#define x