← Back to branch summary

~ubuntu-branches/ubuntu/wily/grass/wily

« back to all changes in this revision

Viewing changes to vector/v.net.salesman/main.c

Tags: 7.0.0~rc1+ds1-1~exp1
http://bugs.debian.org/776349
* New upstream release candidate.
* Repack upstream tarball, remove precompiled Python objects.
* Add upstream metadata.
* Update gbp.conf and Vcs-Git URL to use the experimental branch.
* Update watch file for GRASS 7.0.
* Drop build dependencies for Tcl/Tk, add build dependencies:
  python-numpy, libnetcdf-dev, netcdf-bin, libblas-dev, liblapack-dev
* Update Vcs-Browser URL to use cgit instead of gitweb.
* Update paths to use grass70.
* Add configure options: --with-netcdf, --with-blas, --with-lapack,
  remove --with-tcltk-includes.
* Update patches for GRASS 7.
* Update copyright file, changes:
  - Update copyright years
  - Group files by license
  - Remove unused license sections
* Add patches for various typos.
* Fix desktop file with patch instead of d/rules.
* Use minimal dh rules.
* Bump Standards-Version to 3.9.6, no changes.
* Use dpkg-maintscript-helper to replace directories with symlinks.
  (closes: #776349)
* Update my email to use @debian.org address.

Show diffs side-by-side

added added

removed removed

Lines of Context:
vector/v.net.salesman/main.c
3
3
 * 
4
4
 *  MODULE:       v.net.salesman
5
5
 *  
6
 
 *  AUTHOR(S):    Radim Blazek
 
6
 *  AUTHOR(S):    Radim Blazek, Markus Metz
7
7
 *                
8
8
 *  PURPOSE:      Create a cycle connecting given nodes.
9
 
 *                
10
 
 *  COPYRIGHT:    (C) 2001 by the GRASS Development Team
 
9
 *
 
10
 *  COPYRIGHT:    (C) 2001-2011,2014 by the GRASS Development Team
11
11
 * 
12
12
 *                This program is free software under the 
13
13
 *                GNU General Public License (>=v2). 
14
14
 *                Read the file COPYING that comes with GRASS
15
15
 *                for details.
16
 
 * 
17
 
 **************************************************************/
 
16
 *
 
17
 ****************************************************************/
18
18
#include <stdlib.h>
19
19
#include <string.h>
20
20
#include <time.h>
21
21
#include <grass/gis.h>
22
 
#include <grass/Vect.h>
 
22
#include <grass/vector.h>
23
23
#include <grass/dbmi.h>
24
24
#include <grass/glocale.h>
25
25
 
35
35
int nnodes;                     /* number of nodes */
36
36
int *cities;                    /* array of cities */
37
37
int *cused;                     /* city is in cycle */
38
 
COST **costs;                   /* pointer to array of pointers to arrays of sorted costs */
 
38
COST **costs;                   /* pointer to array of pointers to arrays of sorted forward costs */
 
39
COST **bcosts;                  /* pointer to array of pointers to arrays of sorted backward costs */
39
40
int *cycle;                     /* path */
40
41
int ncyc = 0;                   /* number of cities in cycle */
 
42
int debug_level;
41
43
 
42
44
int cmp(const void *, const void *);
43
45
 
54
56
        cycle[0] = city;
55
57
    }
56
58
    else {
 
59
        /* for a large number of cities this will become slow */
57
60
        for (j = ncyc - 1; j > after; j--)
58
61
            cycle[j + 1] = cycle[j];
59
62
 
62
65
    cused[city] = 1;
63
66
    ncyc++;
64
67
 
65
 
    G_debug(2, "Cycle:\n");
66
 
    for (i = 0; i < ncyc; i++) {
67
 
        G_debug(2, "%d: %d: %d\n", i, cycle[i], cities[cycle[i]]);
68
 
    }
69
 
 
 
68
    if (debug_level >= 2) {
 
69
        G_debug(2, "Cycle:");
 
70
        for (i = 0; i < ncyc; i++) {
 
71
            G_debug(2, "%d: %d: %d", i, cycle[i], cities[cycle[i]]);
 
72
        }
 
73
    }
 
74
 
 
75
}
 
76
 
 
77
/* like Vect_list_append(), but allows duplicates */
 
78
int tsp_list_append(struct ilist *list, int val)
 
79
{
 
80
    size_t size;
 
81
 
 
82
    if (list == NULL)
 
83
        return 1;
 
84
 
 
85
    if (list->n_values == list->alloc_values) {
 
86
        size = (list->n_values + 1000) * sizeof(int);
 
87
        list->value = (int *)G_realloc((void *)list->value, size);
 
88
        list->alloc_values = list->n_values + 1000;
 
89
    }
 
90
 
 
91
    list->value[list->n_values] = val;
 
92
    list->n_values++;
 
93
 
 
94
    return 0;
70
95
}
71
96
 
72
97
 
75
100
    int i, j, k, ret, city, city1;
76
101
    int nlines, type, ltype, afield, tfield, geo, cat;
77
102
    int node, node1, node2, line;
78
 
    struct Option *map, *output, *afield_opt, *tfield_opt, *afcol, *type_opt,
79
 
        *term_opt;
 
103
    double **cost_cache;                        /* pointer to array of pointers to arrays of cached costs */
 
104
    struct Option *map, *output, *afield_opt, *tfield_opt, *afcol, *abcol,
 
105
        *seq, *type_opt, *term_opt;
80
106
    struct Flag *geo_f;
81
107
    struct GModule *module;
82
 
    char *mapset;
83
108
    struct Map_info Map, Out;
84
109
    struct ilist *TList;        /* list of terminal nodes */
85
110
    struct ilist *List;
89
114
    struct cat_list *Clist;
90
115
    struct line_cats *Cats;
91
116
    struct line_pnts *Points;
 
117
    const char *dstr;
 
118
    const char *seqname;
 
119
    int seq2stdout;
 
120
    FILE *fp;
92
121
 
93
122
    /* Initialize the GIS calls */
94
123
    G_gisinit(argv[0]);
95
124
 
96
125
    module = G_define_module();
97
 
    module->keywords = _("vector, networking");
 
126
    G_add_keyword(_("vector"));
 
127
    G_add_keyword(_("network"));
 
128
    G_add_keyword(_("salesman"));
98
129
    module->label =
99
130
        _("Creates a cycle connecting given nodes (Traveling salesman problem).");
100
131
    module->description =
104
135
    map = G_define_standard_option(G_OPT_V_INPUT);
105
136
    output = G_define_standard_option(G_OPT_V_OUTPUT);
106
137
 
 
138
    afield_opt = G_define_standard_option(G_OPT_V_FIELD);
 
139
    afield_opt->key = "arc_layer";
 
140
    afield_opt->label = _("Arc layer");
 
141
 
107
142
    type_opt = G_define_standard_option(G_OPT_V_TYPE);
 
143
    type_opt->key = "arc_type";
108
144
    type_opt->options = "line,boundary";
109
145
    type_opt->answer = "line,boundary";
110
146
    type_opt->description = _("Arc type");
111
147
 
112
 
    afield_opt = G_define_standard_option(G_OPT_V_FIELD);
113
 
    afield_opt->key = "alayer";
114
 
    afield_opt->description = _("Arc layer");
115
 
 
116
148
    tfield_opt = G_define_standard_option(G_OPT_V_FIELD);
117
 
    tfield_opt->key = "nlayer";
 
149
    tfield_opt->key = "node_layer";
118
150
    tfield_opt->answer = "2";
119
 
    tfield_opt->description = _("Node layer (used for cities)");
 
151
    tfield_opt->label = _("Node layer (used for cities)");
120
152
 
121
153
    afcol = G_define_option();
122
 
    afcol->key = "acolumn";
 
154
    afcol->key = "arc_column";
123
155
    afcol->type = TYPE_STRING;
124
156
    afcol->required = NO;
125
 
    afcol->description = _("Arcs' cost column (for both directions)");
 
157
    afcol->description = _("Arc forward/both direction(s) cost column (number)");
 
158
 
 
159
    abcol = G_define_option();
 
160
    abcol->key = "arc_backward_column";
 
161
    abcol->type = TYPE_STRING;
 
162
    abcol->required = NO;
 
163
    abcol->description = _("EXPERIMENTAL: Arc backward direction cost column (number)");
 
164
 
 
165
    seq = G_define_standard_option(G_OPT_F_OUTPUT);
 
166
    seq->key = "sequence";
 
167
    seq->type = TYPE_STRING;
 
168
    seq->required = NO;
 
169
    seq->description = _("Name for output file holding node sequence (\"-\" for stdout)");
126
170
 
127
171
    term_opt = G_define_standard_option(G_OPT_V_CATS);
128
 
    term_opt->key = "ccats";
 
172
    term_opt->key = "center_cats";
129
173
    term_opt->required = YES;
130
174
    term_opt->description = _("Categories of points ('cities') on nodes "
131
175
                              "(layer is specified by nlayer)");
153
197
    tfield = atoi(tfield_opt->answer);
154
198
    Vect_str_to_cat_list(term_opt->answer, Clist);
155
199
 
156
 
    G_debug(1, "Imput categories:\n");
157
 
    for (i = 0; i < Clist->n_ranges; i++) {
158
 
        G_debug(1, "%d - %d\n", Clist->min[i], Clist->max[i]);
 
200
    dstr = G_getenv_nofatal("DEBUG");
 
201
 
 
202
    if (dstr != NULL)
 
203
        debug_level = atoi(dstr);
 
204
    else
 
205
        debug_level = 0;
 
206
 
 
207
    if (debug_level >= 1) {
 
208
        G_debug(1, "Input categories:");
 
209
        for (i = 0; i < Clist->n_ranges; i++) {
 
210
            G_debug(1, "%d - %d", Clist->min[i], Clist->max[i]);
 
211
        }
159
212
    }
160
213
 
161
214
    if (geo_f->answer)
163
216
    else
164
217
        geo = 0;
165
218
 
166
 
    Vect_check_input_output_name(map->answer, output->answer, GV_FATAL_EXIT);
167
 
 
168
 
    mapset = G_find_vector2(map->answer, NULL);
169
 
 
170
 
    if (mapset == NULL)
171
 
        G_fatal_error(_("Vector map <%s> not found"), map->answer);
 
219
    Vect_check_input_output_name(map->answer, output->answer, G_FATAL_EXIT);
172
220
 
173
221
    Vect_set_open_level(2);
174
 
    Vect_open_old(&Map, map->answer, mapset);
 
222
 
 
223
    if (Vect_open_old(&Map, map->answer, "") < 0)
 
224
        G_fatal_error(_("Unable to open vector map <%s>"), map->answer);
 
225
 
175
226
    nnodes = Vect_get_num_nodes(&Map);
 
227
    nlines = Vect_get_num_lines(&Map);
176
228
 
177
229
    /* Create list of terminals based on list of categories */
178
 
    for (i = 1; i <= nnodes; i++) {
179
 
        nlines = Vect_get_node_n_lines(&Map, i);
180
 
        for (j = 0; j < nlines; j++) {
181
 
            line = abs(Vect_get_node_line(&Map, i, j));
182
 
            ltype = Vect_read_line(&Map, NULL, Cats, line);
183
 
            if (!(ltype & GV_POINT))
184
 
                continue;
185
 
            if (!(Vect_cat_get(Cats, tfield, &cat)))
186
 
                continue;
187
 
            if (Vect_cat_in_cat_list(cat, Clist)) {
188
 
                Vect_list_append(TList, i);
 
230
    for (i = 1; i <= nlines; i++) {
 
231
 
 
232
        ltype = Vect_get_line_type(&Map, i);
 
233
        if (!(ltype & GV_POINT))
 
234
            continue;
 
235
 
 
236
        Vect_read_line(&Map, Points, Cats, i);
 
237
        if (!(Vect_cat_get(Cats, tfield, &cat)))
 
238
            continue;
 
239
        if (Vect_cat_in_cat_list(cat, Clist)) {
 
240
            node = Vect_find_node(&Map, Points->x[0], Points->y[0], Points->z[0], 0, 0);
 
241
            if (!node) {
 
242
                G_warning(_("Point is not connected to the network"));
189
243
            }
 
244
            else
 
245
                tsp_list_append(TList, node);
190
246
        }
191
247
    }
 
248
 
192
249
    ncities = TList->n_values;
193
 
    G_message(_("Number of cities: [%d]"), ncities);
 
250
    G_message(_("Number of cities: %d"), ncities);
194
251
    if (ncities < 2)
195
252
        G_fatal_error(_("Not enough cities (< 2)"));
196
253
 
198
255
    cities = (int *)G_malloc(ncities * sizeof(int));
199
256
    cused = (int *)G_malloc(ncities * sizeof(int));
200
257
    for (i = 0; i < ncities; i++) {
201
 
        G_debug(1, "%d\n", TList->value[i]);
 
258
        G_debug(1, "%d", TList->value[i]);
202
259
        cities[i] = TList->value[i];
203
260
        cused[i] = 0;           /* not in cycle */
204
261
    }
205
262
 
206
263
    costs = (COST **) G_malloc(ncities * sizeof(COST *));
207
 
    costs = (COST **) G_malloc((ncities - 1) * sizeof(COST *));
208
264
    for (i = 0; i < ncities; i++) {
209
265
        costs[i] = (COST *) G_malloc(ncities * sizeof(COST));
210
266
    }
 
267
    cost_cache = (double **) G_malloc(ncities * sizeof(double *));
 
268
    for (i = 0; i < ncities; i++) {
 
269
        cost_cache[i] = (double *) G_malloc(ncities * sizeof(double));
 
270
    }
 
271
    if (abcol->answer) {
 
272
        bcosts = (COST **) G_malloc(ncities * sizeof(COST *));
 
273
        for (i = 0; i < ncities; i++) {
 
274
            bcosts[i] = (COST *) G_malloc(ncities * sizeof(COST));
 
275
        }
 
276
    }
 
277
    else
 
278
        bcosts = NULL;
211
279
 
212
280
    cycle = (int *)G_malloc((ncities + 1) * sizeof(int));       /* + 1 is for output cycle */
213
281
 
214
282
    /* Build graph */
215
 
    Vect_net_build_graph(&Map, type, afield, 0, afcol->answer, NULL, NULL,
 
283
    Vect_net_build_graph(&Map, type, afield, 0, afcol->answer, abcol->answer, NULL,
216
284
                         geo, 0);
217
285
 
218
286
    /* Create sorted lists of costs */
 
287
    /* for a large number of cities this will become very slow, can not be fixed */
 
288
    G_message(_("Creating cost cache..."));
219
289
    for (i = 0; i < ncities; i++) {
 
290
        G_percent(i, ncities, 2);
220
291
        k = 0;
221
292
        for (j = 0; j < ncities; j++) {
 
293
            cost_cache[i][j] = 0.0;
222
294
            if (i == j)
223
295
                continue;
 
296
 
224
297
            ret =
225
298
                Vect_net_shortest_path(&Map, cities[i], cities[j], NULL,
226
299
                                       &cost);
227
 
            if (ret == -1)
228
 
                G_fatal_error(_("Destination node [%d] is unreachable "
229
 
                                "from node [%d]"), cities[i], cities[j]);
230
 
 
 
300
 
 
301
            if (ret == -1) {
 
302
                double coor_x, coor_y, coor_z;
 
303
                int cat1, cat2;
 
304
                
 
305
                Vect_get_node_coor(&Map, cities[i], &coor_x, &coor_y, &coor_z);
 
306
                line = Vect_find_line(&Map, coor_x, coor_y, coor_z, GV_POINT, 0, 0, 0);
 
307
                
 
308
                if (!line)
 
309
                    G_fatal_error(_("No point at node %d"), cities[i]);
 
310
 
 
311
                Vect_read_line(&Map, Points, Cats, line);
 
312
                if (!(Vect_cat_get(Cats, tfield, &cat1)))
 
313
                    G_fatal_error(_("No category for point at node %d"), cities[i]);
 
314
 
 
315
                Vect_get_node_coor(&Map, cities[j], &coor_x, &coor_y, &coor_z);
 
316
                line = Vect_find_line(&Map, coor_x, coor_y, coor_z, GV_POINT, 0, 0, 0);
 
317
                
 
318
                if (!line)
 
319
                    G_fatal_error(_("No point at node %d"), cities[j]);
 
320
 
 
321
                Vect_read_line(&Map, Points, Cats, line);
 
322
                if (!(Vect_cat_get(Cats, tfield, &cat2)))
 
323
                    G_fatal_error(_("No category for point at node %d"), cities[j]);
 
324
 
 
325
                G_fatal_error(_("Destination node [cat %d] is unreachable "
 
326
                                "from node [cat %d]"), cat1, cat2);
 
327
            }
 
328
 
 
329
            /* add to directional cost cache: from, to, cost */
231
330
            costs[i][k].city = j;
232
331
            costs[i][k].cost = cost;
 
332
            cost_cache[i][j] = cost;
 
333
 
233
334
            k++;
234
335
        }
235
336
        qsort((void *)costs[i], k, sizeof(COST), cmp);
236
337
    }
237
 
    /* debug: print sorted */
238
 
    for (i = 0; i < ncities; i++) {
239
 
        for (j = 0; j < ncities - 1; j++) {
240
 
            city = costs[i][j].city;
241
 
            G_debug(2, "%d -> %d = %f\n", cities[i], cities[city],
242
 
                    costs[i][j].cost);
243
 
        }
244
 
    }
245
 
 
 
338
    G_percent(1, 1, 2);
 
339
    
 
340
    if (bcosts) {
 
341
        for (i = 0; i < ncities; i++) {
 
342
            /* this should be fast, no need for G_percent() */
 
343
            k = 0;
 
344
            for (j = 0; j < ncities; j++) {
 
345
                if (i == j)
 
346
                    continue;
 
347
                    
 
348
                bcosts[i][k].city = j;
 
349
                bcosts[i][k].cost = cost_cache[j][i];
 
350
 
 
351
                k++;
 
352
            }
 
353
            qsort((void *)bcosts[i], k, sizeof(COST), cmp);
 
354
        }
 
355
    }
 
356
    
 
357
    if (debug_level >= 2) {
 
358
        /* debug: print sorted */
 
359
        for (i = 0; i < ncities; i++) {
 
360
            for (j = 0; j < ncities - 1; j++) {
 
361
                city = costs[i][j].city;
 
362
                G_debug(2, "%d -> %d = %f", cities[i], cities[city],
 
363
                        costs[i][j].cost);
 
364
            }
 
365
        }
 
366
    }
 
367
 
 
368
    G_message(_("Searching for the shortest cycle..."));
246
369
    /* find 2 cities with largest distance */
247
 
    cost = -1;
 
370
    cost = city = -1;
248
371
    for (i = 0; i < ncities; i++) {
249
372
        tmpcost = costs[i][ncities - 2].cost;
250
373
        if (tmpcost > cost) {
252
375
            city = i;
253
376
        }
254
377
    }
255
 
    G_debug(2, "biggest costs %d - %d\n", city,
 
378
    G_debug(2, "biggest costs %d - %d", city,
256
379
            costs[city][ncities - 2].city);
257
380
 
258
 
    /* add this 2 cities to array */
 
381
    /* add these 2 cities to array */
259
382
    add_city(city, -1);
260
383
    add_city(costs[city][ncities - 2].city, 0);
261
384
 
262
385
    /* In each step, find not used city, with biggest cost to any used city, and insert 
263
386
     *  into cycle between 2 nearest nodes */
 
387
    /* for a large number of cities this will become very slow, can be fixed */
264
388
    for (i = 0; i < ncities - 2; i++) {
 
389
        G_percent(i, ncities - 3, 1);
265
390
        cost = -1;
266
 
        G_debug(2, "---- %d ----\n", i);
 
391
        G_debug(2, "---- city %d ----", i);
267
392
        for (j = 0; j < ncities; j++) {
268
393
            if (cused[j] == 1)
269
394
                continue;
270
395
            tmpcost = 0;
271
396
            for (k = 0; k < ncities - 1; k++) {
272
 
                G_debug(2, "? %d (%d) - %d (%d) \n", j, cnode(j),
 
397
                G_debug(2, "forward? %d (%d) - %d (%d)", j, cnode(j),
273
398
                        costs[j][k].city, cnode(costs[j][k].city));
274
399
                if (!cused[costs[j][k].city])
275
400
                    continue;   /* only used */
 
401
                /* directional costs j -> k */
276
402
                tmpcost += costs[j][k].cost;
277
403
                break;          /* first nearest */
278
404
            }
279
 
            G_debug(2, "    cost = %f x %f\n", tmpcost, cost);
 
405
            /* forward/backward: tmpcost = min(fcost) + min(bcost) */
 
406
            if (bcosts) {
 
407
                for (k = 0; k < ncities - 1; k++) {
 
408
                    G_debug(2, "backward? %d (%d) - %d (%d)", j, cnode(j),
 
409
                            bcosts[j][k].city, cnode(bcosts[j][k].city));
 
410
                    if (!cused[bcosts[j][k].city])
 
411
                        continue;       /* only used */
 
412
                    /* directional costs k -> j */
 
413
                    tmpcost += bcosts[j][k].cost;
 
414
                    break;              /* first nearest */
 
415
                }
 
416
            }
 
417
 
 
418
            G_debug(2, "    cost = %f x %f", tmpcost, cost);
280
419
            if (tmpcost > cost) {
281
420
                cost = tmpcost;
282
421
                city = j;
283
422
            }
284
423
        }
285
 
        G_debug(2, "add %d\n", city);
 
424
        G_debug(2, "add city %d", city);
286
425
 
287
 
        /* add to cycle on lovest costs */
288
 
        cycle[ncyc] = cycle[0]; /* tmp for cycle */
 
426
        /* add to cycle on lowest costs */
 
427
        cycle[ncyc] = cycle[0]; /* temporarily close the cycle */
289
428
        cost = PORT_DOUBLE_MAX;
 
429
        city1 = 0;
290
430
        for (j = 0; j < ncyc; j++) {
291
 
            node1 = cities[cycle[j]];
292
 
            node2 = cities[cycle[j + 1]];
293
 
            ret = Vect_net_shortest_path(&Map, node1, node2, NULL, &tcost);
 
431
            /* cost from j to j + 1 (directional) */
 
432
            /* get cost from directional cost cache */
 
433
            tcost = cost_cache[cycle[j]][cycle[j + 1]];
294
434
            tmpcost = -tcost;
295
 
            node1 = cities[cycle[j]];
296
 
            node2 = cities[city];
297
 
            ret = Vect_net_shortest_path(&Map, node1, node2, NULL, &tcost);
298
 
            tmpcost += tcost;
299
 
            node1 = cities[cycle[j + 1]];
300
 
            node2 = cities[city];
301
 
            ret = Vect_net_shortest_path(&Map, node1, node2, NULL, &tcost);
302
 
            tmpcost += tcost;
303
 
 
304
 
            G_debug(2, "? %d - %d cost = %f x %f\n", node1, node2, tmpcost,
 
435
 
 
436
            /* check insertion of city between j and j + 1 */
 
437
 
 
438
            /* cost from j to city (directional) */
 
439
            /* get cost from directional cost cache */
 
440
            tcost = cost_cache[cycle[j]][city];
 
441
            tmpcost += tcost;
 
442
            /* cost from city to j + 1 (directional) */
 
443
            /* get cost from directional cost cache */
 
444
            tcost = cost_cache[city][cycle[j + 1]];
 
445
            tmpcost += tcost;
 
446
            
 
447
            /* tmpcost must always be > 0 */
 
448
 
 
449
            G_debug(2, "? %d - %d cost = %f x %f", node1, node2, tmpcost,
305
450
                    cost);
 
451
            /* always true for j = 0 */
306
452
            if (tmpcost < cost) {
307
453
                city1 = j;
308
454
                cost = tmpcost;
309
455
            }
310
456
        }
311
 
 
312
457
        add_city(city, city1);
313
 
 
314
458
    }
 
459
    
 
460
    /* TODO: optimize tour (some Lin-Kernighan method) */
315
461
 
316
 
    /* Print */
317
 
    G_debug(2, "Cycle:\n");
318
 
    for (i = 0; i < ncities; i++) {
319
 
        G_debug(2, "%d: %d: %d\n", i, cycle[i], cities[cycle[i]]);
 
462
    if (debug_level >= 2) {
 
463
        /* debug print */
 
464
        G_debug(2, "Cycle:");
 
465
        for (i = 0; i < ncities; i++) {
 
466
            G_debug(2, "%d: %d: %d", i, cycle[i], cities[cycle[i]]);
 
467
        }
320
468
    }
321
469
 
322
470
    /* Create list of arcs */
323
 
    cycle[ncities] = cycle[0];
 
471
    cycle[ncities] = cycle[0];  /* close the cycle */
 
472
    cost = 0.0;
324
473
    for (i = 0; i < ncities; i++) {
325
474
        node1 = cities[cycle[i]];
326
475
        node2 = cities[cycle[i + 1]];
327
 
        G_debug(2, " %d -> %d\n", node1, node2);
 
476
        G_debug(2, " %d -> %d", node1, node2);
328
477
        ret = Vect_net_shortest_path(&Map, node1, node2, List, NULL);
 
478
        cost += cost_cache[cycle[i]][cycle[i + 1]];
329
479
        for (j = 0; j < List->n_values; j++) {
330
480
            line = abs(List->value[j]);
331
 
            Vect_list_append(StArcs, line);
 
481
            /* Vect_list_append() appends only if value not yet present !!! 
 
482
             * this breaks the correct sequence */
 
483
            tsp_list_append(StArcs, line);
332
484
            Vect_get_line_nodes(&Map, line, &node1, &node2);
333
 
            Vect_list_append(StNodes, node1);
334
 
            Vect_list_append(StNodes, node2);
 
485
            tsp_list_append(StNodes, node1);
 
486
            tsp_list_append(StNodes, node2);
335
487
        }
336
488
    }
337
489
 
338
 
 
339
 
 
340
490
    /* Write arcs to new map */
341
 
    Vect_open_new(&Out, output->answer, Vect_is_3d(&Map));
 
491
    if (Vect_open_new(&Out, output->answer, Vect_is_3d(&Map)) < 0)
 
492
        G_fatal_error(_("Unable to create vector map <%s>"), output->answer);
 
493
 
342
494
    Vect_hist_command(&Out);
343
495
 
344
 
    fprintf(stdout, "\nCycle:\n");
345
 
    fprintf(stdout, "Arcs' categories (layer %d, %d arcs):\n", afield,
 
496
    G_verbose_message(_("Cycle with total cost %.3f"), cost);
 
497
    G_debug(2, "Arcs' categories (layer %d, %d arcs):", afield,
346
498
            StArcs->n_values);
 
499
 
347
500
    for (i = 0; i < StArcs->n_values; i++) {
348
501
        line = StArcs->value[i];
349
502
        ltype = Vect_read_line(&Map, Points, Cats, line);
350
503
        Vect_write_line(&Out, ltype, Points, Cats);
351
504
        Vect_cat_get(Cats, afield, &cat);
352
 
        if (i > 0)
353
 
            printf(",");
354
 
        fprintf(stdout, "%d", cat);
355
 
    }
356
 
    fprintf(stdout, "\n\n");
357
 
 
358
 
    fprintf(stdout, "Nodes' categories (layer %d, %d nodes):\n", tfield,
359
 
            StNodes->n_values);
 
505
        G_debug(2, "%d. arc: cat %d", i + 1, cat);
 
506
    }
 
507
    
 
508
    seq2stdout = 0;
 
509
    seqname = NULL;
 
510
    if (seq->answer) {
 
511
        if (strcmp(seq->answer, "-")) {
 
512
            seqname = seq->answer;
 
513
        }
 
514
        else {
 
515
            seqname = G_tempfile();
 
516
            seq2stdout = 1;
 
517
        }
 
518
 
 
519
        fp = fopen(seqname, "w");
 
520
        if (!fp)
 
521
            G_fatal_error(_("Unable to open file '%s' for writing"),
 
522
                          seqname);
 
523
 
 
524
        fprintf(fp, "sequence;category;cost_to_next\n");
 
525
    }
 
526
    else
 
527
        fp = NULL;
 
528
 
360
529
    k = 0;
361
 
    for (i = 0; i < StNodes->n_values; i++) {
362
 
        node = StNodes->value[i];
363
 
        nlines = Vect_get_node_n_lines(&Map, node);
364
 
        for (j = 0; j < nlines; j++) {
365
 
            line = abs(Vect_get_node_line(&Map, node, j));
366
 
            ltype = Vect_read_line(&Map, Points, Cats, line);
367
 
            if (!(ltype & GV_POINT))
368
 
                continue;
369
 
            if (!(Vect_cat_get(Cats, tfield, &cat)))
370
 
                continue;
371
 
            Vect_write_line(&Out, ltype, Points, Cats);
372
 
            if (k > 0)
373
 
                fprintf(stdout, ",");
374
 
            fprintf(stdout, "%d", cat);
375
 
            k++;
 
530
    /* this writes out only user-selected nodes, not all visited nodes */
 
531
    G_debug(2, "Nodes' categories (layer %d, %d nodes):", tfield,
 
532
            ncities);
 
533
    for (i = 0; i < ncities; i++) {
 
534
        double coor_x, coor_y, coor_z;
 
535
        
 
536
        node = cities[cycle[i]];
 
537
        Vect_get_node_coor(&Map, node, &coor_x, &coor_y, &coor_z);
 
538
        line = Vect_find_line(&Map, coor_x, coor_y, coor_z, GV_POINT, 0, 0, 0);
 
539
        
 
540
        if (!line)
 
541
            continue;
 
542
 
 
543
        ltype = Vect_read_line(&Map, Points, Cats, line);
 
544
        if (!(ltype & GV_POINT))
 
545
            continue;
 
546
        if (!(Vect_cat_get(Cats, tfield, &cat)))
 
547
            continue;
 
548
        Vect_write_line(&Out, ltype, Points, Cats);
 
549
        k++;
 
550
        if (fp) {
 
551
            fprintf(fp, "%d;%d;%.3f\n", k, cat, cost_cache[cycle[i]][cycle[i + 1]]);
376
552
        }
 
553
 
 
554
        G_debug(2, "%d. node: cat %d", k, cat);
377
555
    }
378
 
    fprintf(stdout, "\n\n");
379
556
 
380
557
    Vect_build(&Out);
381
558
 
385
562
    Vect_close(&Map);
386
563
    Vect_close(&Out);
387
564
 
 
565
    if (fp) {
 
566
        fclose(fp);
 
567
        if (seq2stdout) {
 
568
            char buf[2000];
 
569
 
 
570
            /* spacer to previous output to stderr */
 
571
            G_message(" ");
 
572
 
 
573
            fp = fopen(seqname, "r");
 
574
            while (G_getl2(buf, 2000, fp) != 0)
 
575
                fprintf(stdout, "%s\n", buf);
 
576
 
 
577
            fclose(fp);
 
578
            remove(seqname);
 
579
        }
 
580
    }
 
581
 
388
582
    exit(EXIT_SUCCESS);
389
583
}
390
584