1
2 /****************************************************************
3 *
4 * MODULE: v.net.flow
5 *
6 * AUTHOR(S): Daniel Bundala
7 *
8 * PURPOSE: Max flow and min cut between two sets of nodes
9 *
10 * COPYRIGHT: (C) 2002-2005 by the GRASS Development Team
11 *
12 * This program is free software under the
13 * GNU General Public License (>=v2).
14 * Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 ****************************************************************/
18
19 #include <stdio.h>
20 #include <stdlib.h>
21 #include <grass/gis.h>
22 #include <grass/vector.h>
23 #include <grass/glocale.h>
24 #include <grass/dbmi.h>
25 #include <grass/neta.h>
26
27
main(int argc,char * argv[])28 int main(int argc, char *argv[])
29 {
30 struct Map_info In, Out, cut_map;
31 static struct line_pnts *Points;
32 struct line_cats *Cats;
33 struct GModule *module; /* GRASS module for parsing arguments */
34 struct Option *map_in, *map_out, *cut_out;
35 struct Option *afield_opt, *nfield_opt, *abcol, *afcol, *ncol;
36 struct Option *catsource_opt, *wheresource_opt;
37 struct Option *catsink_opt, *wheresink_opt;
38 int with_z;
39 int afield, nfield, mask_type;
40 struct varray *varray_source, *varray_sink;
41 dglGraph_s *graph;
42 int i, nlines, *flow, total_flow;
43 struct ilist *source_list, *sink_list, *cut;
44 int find_cut;
45
46 char buf[2000];
47
48 /* Attribute table */
49 dbString sql;
50 dbDriver *driver;
51 struct field_info *Fi;
52
53 /* initialize GIS environment */
54 G_gisinit(argv[0]); /* reads grass env, stores program name to G_program_name() */
55
56 /* initialize module */
57 module = G_define_module();
58 G_add_keyword(_("vector"));
59 G_add_keyword(_("network"));
60 G_add_keyword(_("flow"));
61 module->description =
62 _("Computes the maximum flow between two sets of nodes in the network.");
63
64 /* Define the different options as defined in gis.h */
65 map_in = G_define_standard_option(G_OPT_V_INPUT);
66
67 afield_opt = G_define_standard_option(G_OPT_V_FIELD);
68 afield_opt->key = "arc_layer";
69 afield_opt->answer = "1";
70 afield_opt->label = _("Arc layer");
71 afield_opt->guisection = _("Cost");
72
73 nfield_opt = G_define_standard_option(G_OPT_V_FIELD);
74 nfield_opt->key = "node_layer";
75 nfield_opt->answer = "2";
76 nfield_opt->label = _("Node layer");
77 nfield_opt->guisection = _("Cost");
78
79 map_out = G_define_standard_option(G_OPT_V_OUTPUT);
80
81 cut_out = G_define_standard_option(G_OPT_V_OUTPUT);
82 cut_out->key = "cut";
83 cut_out->description =
84 _("Name for output vector map containing a minimum cut");
85
86 afcol = G_define_standard_option(G_OPT_DB_COLUMN);
87 afcol->key = "arc_column";
88 afcol->required = NO;
89 afcol->description =
90 _("Arc forward/both direction(s) cost column (number)");
91 afcol->guisection = _("Cost");
92
93 abcol = G_define_standard_option(G_OPT_DB_COLUMN);
94 abcol->key = "arc_backward_column";
95 abcol->required = NO;
96 abcol->description = _("Arc backward direction cost column (number)");
97 abcol->guisection = _("Cost");
98
99 ncol = G_define_standard_option(G_OPT_DB_COLUMN);
100 ncol->key = "node_column";
101 ncol->required = NO;
102 ncol->description = _("Node cost column (number)");
103 ncol->guisection = _("Cost");
104
105 catsource_opt = G_define_standard_option(G_OPT_V_CATS);
106 catsource_opt->key = "source_cats";
107 catsource_opt->label = _("Source category values");
108 catsource_opt->guisection = _("Source");
109
110 wheresource_opt = G_define_standard_option(G_OPT_DB_WHERE);
111 wheresource_opt->key = "source_where";
112 wheresource_opt->label =
113 _("Source WHERE conditions of SQL statement without 'where' keyword");
114 wheresource_opt->guisection = _("Source");
115
116 catsink_opt = G_define_standard_option(G_OPT_V_CATS);
117 catsink_opt->key = "sink_cats";
118 catsink_opt->label = _("Sink category values");
119 catsink_opt->guisection = _("Sink");
120
121 wheresink_opt = G_define_standard_option(G_OPT_DB_WHERE);
122 wheresink_opt->key = "sink_where";
123 wheresink_opt->label =
124 _("Sink WHERE conditions of SQL statement without 'where' keyword");
125 wheresink_opt->guisection = _("Sink");
126
127 /* options and flags parser */
128 if (G_parser(argc, argv))
129 exit(EXIT_FAILURE);
130 find_cut = (cut_out->answer[0]);
131 /* TODO: make an option for this */
132 mask_type = GV_LINE | GV_BOUNDARY;
133
134 Points = Vect_new_line_struct();
135 Cats = Vect_new_cats_struct();
136
137 Vect_check_input_output_name(map_in->answer, map_out->answer,
138 G_FATAL_EXIT);
139
140 Vect_set_open_level(2);
141
142 if (1 > Vect_open_old(&In, map_in->answer, ""))
143 G_fatal_error(_("Unable to open vector map <%s>"), map_in->answer);
144
145 with_z = Vect_is_3d(&In);
146
147 if (0 > Vect_open_new(&Out, map_out->answer, with_z)) {
148 Vect_close(&In);
149 G_fatal_error(_("Unable to create vector map <%s>"), map_out->answer);
150 }
151
152 if (find_cut && 0 > Vect_open_new(&cut_map, cut_out->answer, with_z)) {
153 Vect_close(&In);
154 Vect_close(&Out);
155 G_fatal_error(_("Unable to create vector map <%s>"), cut_out->answer);
156 }
157
158 /* parse filter option and select appropriate lines */
159 afield = Vect_get_field_number(&In, afield_opt->answer);
160 nfield = Vect_get_field_number(&In, nfield_opt->answer);
161
162 /* Create table */
163 Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
164 Vect_map_add_dblink(&Out, 1, NULL, Fi->table, GV_KEY_COLUMN, Fi->database,
165 Fi->driver);
166 db_init_string(&sql);
167 driver = db_start_driver_open_database(Fi->driver, Fi->database);
168 if (driver == NULL)
169 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
170 Fi->database, Fi->driver);
171 db_set_error_handler_driver(driver);
172
173 sprintf(buf, "create table %s (cat integer, flow double precision)",
174 Fi->table);
175
176 db_set_string(&sql, buf);
177 G_debug(2, "%s", db_get_string(&sql));
178
179 if (db_execute_immediate(driver, &sql) != DB_OK) {
180 db_close_database_shutdown_driver(driver);
181 G_fatal_error(_("Unable to create table: '%s'"), db_get_string(&sql));
182 }
183
184 if (db_create_index2(driver, Fi->table, GV_KEY_COLUMN) != DB_OK)
185 G_warning(_("Cannot create index"));
186
187 if (db_grant_on_table
188 (driver, Fi->table, DB_PRIV_SELECT, DB_GROUP | DB_PUBLIC) != DB_OK)
189 G_fatal_error(_("Cannot grant privileges on table <%s>"), Fi->table);
190
191 db_begin_transaction(driver);
192
193 source_list = Vect_new_list();
194 sink_list = Vect_new_list();
195
196 if (NetA_initialise_varray
197 (&In, nfield, GV_POINT,
198 wheresource_opt->answer, catsource_opt->answer, &varray_source) <= 0) {
199 G_fatal_error(_("No source features selected. "
200 "Please check options '%s', '%s'."),
201 catsource_opt->key, wheresource_opt->key);
202 }
203 if (NetA_initialise_varray
204 (&In, nfield, GV_POINT, wheresink_opt->answer,
205 catsink_opt->answer, &varray_sink) <= 0) {
206 G_fatal_error(_("No sink features selected. "
207 "Please check options '%s', '%s'."),
208 catsink_opt->key, wheresink_opt->key);
209 }
210
211 NetA_varray_to_nodes(&In, varray_source, source_list, NULL);
212 NetA_varray_to_nodes(&In, varray_sink, sink_list, NULL);
213
214 if (source_list->n_values == 0)
215 G_fatal_error(_("No sources"));
216
217 if (sink_list->n_values == 0)
218 G_fatal_error(_("No sinks"));
219
220 Vect_copy_head_data(&In, &Out);
221 Vect_hist_copy(&In, &Out);
222 Vect_hist_command(&Out);
223
224 if (0 != Vect_net_build_graph(&In, mask_type, afield, nfield, afcol->answer, abcol->answer,
225 ncol->answer, 0, 0))
226 G_fatal_error(_("Unable to build graph for vector map <%s>"), Vect_get_full_name(&In));
227
228 graph = Vect_net_get_graph(&In);
229 nlines = Vect_get_num_lines(&In);
230 flow = (int *)G_calloc(nlines + 1, sizeof(int));
231 if (!flow)
232 G_fatal_error(_("Out of memory"));
233
234 total_flow = NetA_flow(graph, source_list, sink_list, flow);
235 G_debug(3, "Max flow: %d", total_flow);
236 if (find_cut) {
237 cut = Vect_new_list();
238 total_flow = NetA_min_cut(graph, source_list, sink_list, flow, cut);
239 G_debug(3, "Min cut: %d", total_flow);
240 }
241
242 G_message(_("Writing the output..."));
243 G_percent_reset();
244 for (i = 1; i <= nlines; i++) {
245 G_percent(i, nlines, 1);
246 int type = Vect_read_line(&In, Points, Cats, i);
247
248 Vect_write_line(&Out, type, Points, Cats);
249 if (type == GV_LINE) {
250 int cat;
251
252 Vect_cat_get(Cats, afield, &cat);
253 if (cat == -1)
254 continue; /*TODO: warning? */
255 sprintf(buf, "insert into %s values (%d, %f)", Fi->table, cat,
256 flow[i] / (double)In.dgraph.cost_multip);
257 db_set_string(&sql, buf);
258 G_debug(3, "%s", db_get_string(&sql));
259
260 if (db_execute_immediate(driver, &sql) != DB_OK) {
261 db_close_database_shutdown_driver(driver);
262 G_fatal_error(_("Cannot insert new record: %s"),
263 db_get_string(&sql));
264 };
265 }
266 }
267
268 if (find_cut) {
269 for (i = 0; i < cut->n_values; i++) {
270 int type = Vect_read_line(&In, Points, Cats, cut->value[i]);
271
272 Vect_write_line(&cut_map, type, Points, Cats);
273 }
274 Vect_destroy_list(cut);
275
276 Vect_build(&cut_map);
277 Vect_close(&cut_map);
278 }
279
280 db_commit_transaction(driver);
281 db_close_database_shutdown_driver(driver);
282
283 G_free(flow);
284 Vect_destroy_list(source_list);
285 Vect_destroy_list(sink_list);
286
287 Vect_build(&Out);
288
289 Vect_close(&In);
290 Vect_close(&Out);
291
292 exit(EXIT_SUCCESS);
293 }
294