GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
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) / 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[labs(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  }
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[labs(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[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
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, 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,
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[labs(id)]) {
218  visited[to] = 1;
219  queue[end++] = to;
220  }
221  }
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 = labs(dglEdgeGet_Id(graph, edge));
240  if (!visited[to] && flow[edge_id] != 0) {
241  Vect_list_append(cut, edge_id);
242  total_flow += abs(flow[labs(edge_id)]);
243  }
244  }
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  */
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 };
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)) {
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  }
328  }
329  dglNode_T_Release(&nt);
330  if (dglFlatten(out) < 0)
331  G_fatal_error(_("GngFlatten error"));
332  return edge_cnt;
333 }
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
Definition: path.h:11
dglInt32_t * dglGetNode(dglGraph_s *pGraph, dglInt32_t nNodeId)
int Vect_reset_list(struct ilist *)
Reset ilist structure.
unsigned char dglByte_t
Definition: type.h:36
dglInt32_t dglNodeGet_Id(dglGraph_s *pGraph, dglInt32_t *pnNode)
void dglNode_T_Release(dglNodeTraverser_s *pT)
dglInt32_t dglEdgeGet_Cost(dglGraph_s *pGraph, dglInt32_t *pnEdge)
int n_values
Number of values in the list.
Definition: gis.h:698
void G_free(void *)
Free allocated memory.
Definition: gis/alloc.c:149
dglInt32_t * dglEdgeGet_Head(dglGraph_s *pGraph, dglInt32_t *pnEdge)
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 dglGet_NodeCount(dglGraph_s *pgraph)
dglInt32_t * dglEdgeGet_Tail(dglGraph_s *pGraph, dglInt32_t *pnEdge)
#define NULL
Definition: ccmath.h:32
#define x
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:270
#define G_calloc(m, n)
Definition: defs/gis.h:113
int dglInitialize(dglGraph_s *pGraph, dglByte_t Version, dglInt32_t NodeAttrSize, dglInt32_t EdgeAttrSize, dglInt32_t *pOpaqueSet)
Definition: queue.h:43
dglInt32_t * dglNode_T_First(dglNodeTraverser_s *pT)
dglInt32_t * dglEdgeset_T_First(dglEdgesetTraverser_s *pT)
long dglInt32_t
Definition: type.h:37
void dglEdgeset_T_Release(dglEdgesetTraverser_s *pT)
int dglAddEdge(dglGraph_s *pGraph, dglInt32_t nHead, dglInt32_t nTail, dglInt32_t nCost, dglInt32_t nEdge)
int Vect_list_append(struct ilist *, int)
Append new item to the end of list if not yet present.
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 * dglNodeGet_Attr(dglGraph_s *pGraph, dglInt32_t *pnNode)
int dglGet_EdgeCount(dglGraph_s *pgraph)
#define _(str)
Definition: glocale.h:10
List of integers.
Definition: gis.h:689
int * value
Array of values.
Definition: gis.h:694
dglInt32_t sign(dglInt32_t x)
Definition: flow.c:25
int dglNode_T_Initialize(dglNodeTraverser_s *pT, dglGraph_s *pGraph)
dglInt32_t * dglEdgeset_T_Next(dglEdgesetTraverser_s *pT)
dglInt32_t * dglNodeGet_OutEdgeset(dglGraph_s *pGraph, dglInt32_t *pnNode)
int dglEdgeset_T_Initialize(dglEdgesetTraverser_s *pT, dglGraph_s *pGraph, dglInt32_t *pnEdgeset)
dglInt32_t * dglNode_T_Next(dglNodeTraverser_s *pT)
int dglFlatten(dglGraph_s *pGraph)
dglInt32_t dglEdgeGet_Id(dglGraph_s *pGraph, dglInt32_t *pnEdge)
int dglGet_NodeAttrSize(dglGraph_s *pgraph)