← Back to branch summary

~ubuntu-branches/ubuntu/vivid/grass/vivid-proposed

« back to all changes in this revision

Viewing changes to vector/v.kcv/main.c

  • Committer: Package Import Robot
  • Author(s): Bas Couwenberg
  • Date: 2015-02-20 23:12:08 UTC
  • mfrom: (8.2.6 experimental)
  • Revision ID: package-import@ubuntu.com-20150220231208-1u6qvqm84v430b10
Tags: 7.0.0-1~exp1
* New upstream release.
* Update python-ctypes-ternary.patch to use if/else instead of and/or.
* Drop check4dev patch, rely on upstream check.
* Add build dependency on libpq-dev to grass-dev for libpq-fe.h.
* Drop patches applied upstream, refresh remaining patches.
* Update symlinks for images switched from jpg to png.

Show diffs side-by-side

added added

removed removed

Lines of Context:
vector/v.kcv/main.c
1
 
/*
2
 
 * from s.kcv
3
 
 * Copyright (C) 1993-1994. James Darrell McCauley.
4
 
 *
5
 
 * Author: James Darrell McCauley darrell@mccauley-usa.com
6
 
 *                                http://mccauley-usa.com/
7
 
 *
8
 
 * This program is free software; you can redistribute it and/or
9
 
 * modify it under the terms of the GNU General Public License
10
 
 * as published by the Free Software Foundation; either version 2
11
 
 * of the License, or (at your option) any later version.
12
 
 *
13
 
 * This program is distributed in the hope that it will be useful,
14
 
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15
 
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16
 
 * GNU General Public License for more details.
17
 
 *
18
 
 * You should have received a copy of the GNU General Public License
19
 
 * along with this program; if not, write to the Free Software
20
 
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
21
 
 *
22
 
 * Modification History:
23
 
 * 4.2B <27 Jan 1994>  fixed RAND_MAX for Solaris 2.3
24
 
 * <13 Sep 2000> released under GPL
25
 
 * 10/2004 Upgrade to 5.7 Radim Blazek
26
 
 */
 
1
 
 
2
/****************************************************************
 
3
 *
 
4
 * MODULE:     v.kcv
 
5
 *
 
6
 * AUTHOR(S):  James Darrell McCauley darrell@mccauley-usa.com (1993-1994)
 
7
 *             OGR support by Martin Landa <landa.martin gmail.com> (2009)
 
8
 *             Speed-up by Jan Vandrol and Jan Ruzicka (2013)
 
9
 *
 
10
 * PURPOSE:    Randomly partition points into test/train sets.
 
11
 *
 
12
 * COPYRIGHT:  (C) 2001-2014 by James Darrell McCauley, and the GRASS Development Team
 
13
 *
 
14
 *             This program is free software under the GNU General
 
15
 *             Public License (>=v2).  Read the file COPYING that
 
16
 *             comes with GRASS for details.
 
17
 *
 
18
 ****************************************************************/
 
19
 
 
20
 
27
21
#include <stdlib.h>
28
 
#include <unistd.h>
29
22
#include <math.h>
30
23
#include <grass/gis.h>
31
24
#include <grass/dbmi.h>
32
 
#include <grass/Vect.h>
 
25
#include <grass/vector.h>
33
26
#include <grass/glocale.h>
34
 
#include "kcv.h"
35
 
 
36
 
#ifndef RAND_MAX
37
 
#define RAND_MAX (pow(2.0,31.0)-1)
38
 
#endif
39
 
#if defined(__CYGWIN__) || defined(__APPLE__) || defined(__MINGW32__)
40
 
double drand48()
41
 
{
42
 
    return (rand() / 32767.0);
43
 
}
44
 
 
45
 
#define srand48(sv) (srand((unsigned)(sv)))
46
 
#else
47
 
double drand48();
48
 
void srand48();
49
 
#endif
50
27
 
51
28
struct Cell_head window;
52
29
 
53
 
int main(int argc, char *argv[])
 
30
int main(int argc, char *argv[]) 
54
31
{
55
 
    int line, nlines, nlinks;
56
 
    double east, north, (*rng) (), max, myrand();
57
 
    int i, j, n, nsites, verbose, np, *p, dcmp();
58
 
    int *pnt_part;
59
 
    struct Map_info In, Out;
60
 
    static struct line_pnts *Points;
 
32
    int i, line, nlines, nlinks;
 
33
    int idx, *line_idx, lines_left;
 
34
    double (*rng)(void) = G_drand48;
 
35
    int nsites, p, np, min_count, spill;
 
36
    struct partition {
 
37
        int id, count, max;
 
38
    } *part;
 
39
    int last_pidx;
 
40
    struct Map_info Map;
 
41
    struct line_pnts *Points;
61
42
    struct line_cats *Cats;
62
 
    char *mapset;
63
43
    struct GModule *module;
64
 
    struct Option *in_opt, *out_opt, *col_opt, *npart_opt;
65
 
    struct Flag *drand48_flag, *q_flag;
66
 
    BOUND_BOX box;
67
 
    double maxdist;
68
 
 
 
44
    struct Option *map_opt, *col_opt, *npart_opt, *field_opt;
 
45
 
69
46
    /* Attributes */
70
47
    struct field_info *Fi;
71
 
    dbDriver *Driver;
 
48
    int layer;
 
49
    dbDriver * Driver;
72
50
    char buf[2000];
73
51
    dbString sql;
74
52
 
 
53
 
75
54
    module = G_define_module();
76
 
    module->keywords = _("vector, statistics");
 
55
    G_add_keyword(_("vector"));
 
56
    G_add_keyword(_("statistics"));
 
57
    G_add_keyword(_("points"));
 
58
    G_add_keyword(_("point pattern"));
77
59
    module->description =
78
60
        _("Randomly partition points into test/train sets.");
79
61
 
80
 
    in_opt = G_define_standard_option(G_OPT_V_INPUT);
81
 
    out_opt = G_define_standard_option(G_OPT_V_OUTPUT);
 
62
    map_opt = G_define_standard_option(G_OPT_V_MAP);
 
63
 
 
64
    field_opt = G_define_standard_option(G_OPT_V_FIELD);
82
65
 
83
66
    npart_opt = G_define_option();
84
 
    npart_opt->key = "k";
 
67
    npart_opt->key = "npartitions";
85
68
    npart_opt->type = TYPE_INTEGER;
86
69
    npart_opt->required = YES;
87
 
    npart_opt->description = _("Number of partitions");
88
 
    npart_opt->options = "1-32767";
 
70
    npart_opt->label = _("Number of partitions");
 
71
    npart_opt->description = _("Must be > 1");
89
72
 
90
 
    col_opt = G_define_option();
91
 
    col_opt->key = "column";
92
 
    col_opt->type = TYPE_STRING;
93
 
    col_opt->required = YES;
94
 
    col_opt->multiple = NO;
 
73
    col_opt = G_define_standard_option(G_OPT_DB_COLUMN);
95
74
    col_opt->answer = "part";
96
75
    col_opt->description =
97
76
        _("Name for new column to which partition number is written");
98
77
 
99
 
    drand48_flag = G_define_flag();
100
 
    drand48_flag->key = 'd';
101
 
    drand48_flag->description = _("Use drand48()");
102
 
 
103
 
    /* please remove in GRASS 7 */
104
 
    q_flag = G_define_flag();
105
 
    q_flag->key = 'q';
106
 
    q_flag->description = _("Quiet");
107
 
 
108
78
    G_gisinit(argv[0]);
109
79
 
110
80
    if (G_parser(argc, argv))
111
81
        exit(EXIT_FAILURE);
112
82
 
113
 
    verbose = (!q_flag->answer);
114
83
    np = atoi(npart_opt->answer);
 
84
    if (np < 2)
 
85
        G_fatal_error(_("'%s' must be > 1"), npart_opt->key);
115
86
 
116
 
    if (drand48_flag->answer) {
117
 
        rng = drand48;
118
 
        max = 1.0;
119
 
        srand48((long)getpid());
120
 
    }
121
 
    else {
122
 
        rng = myrand;
123
 
        max = RAND_MAX;
124
 
        srand(getpid());
125
 
    }
 
87
    /* FIXME - allow seed to be specified for repeatability */
 
88
    G_srand48_auto();
126
89
 
127
90
    Points = Vect_new_line_struct();
128
91
    Cats = Vect_new_cats_struct();
129
92
 
130
93
    /* open input vector */
131
 
    if ((mapset = G_find_vector2(in_opt->answer, "")) == NULL) {
132
 
        G_fatal_error(_("Vector map <%s> not found"), in_opt->answer);
133
 
    }
134
 
 
135
94
    Vect_set_open_level(2);
136
 
    if (Vect_open_old(&In, in_opt->answer, mapset) < 2) {
 
95
    if (Vect_open_old2(&Map, map_opt->answer, G_mapset(), field_opt->answer) < 2) {
137
96
        G_fatal_error(_("Unable to open vector map <%s> at topological level %d"),
138
 
                      in_opt->answer, 2);
 
97
                      map_opt->answer, 2);
139
98
    }
140
99
 
141
 
    Vect_get_map_box(&In, &box);
 
100
    layer = Vect_get_field_number(&Map, field_opt->answer);
 
101
    if (layer <= 0)
 
102
        G_fatal_error(_("Layer number must be positive"));
142
103
 
143
 
    nsites = Vect_get_num_primitives(&In, GV_POINT);
 
104
    nsites = Vect_get_num_primitives(&Map, GV_POINT);
144
105
 
145
106
    if (nsites < np) {
146
 
        G_fatal_error("Sites found: %i\nMore partitions than sites.", nsites);
 
107
        G_fatal_error(_("More partitions than points (%d)"), nsites);
147
108
    }
148
109
 
149
 
    Vect_set_fatal_error(GV_FATAL_PRINT);
150
 
    Vect_open_new(&Out, out_opt->answer, Vect_is_3d(&In));
151
 
 
152
 
    Vect_copy_head_data(&In, &Out);
153
 
    Vect_hist_copy(&In, &Out);
154
 
    Vect_hist_command(&Out);
155
 
 
156
 
    /* Copy vector lines */
157
 
    Vect_copy_map_lines(&In, &Out);
158
 
 
159
 
    /* Copy tables */
160
 
    if (Vect_copy_tables(&In, &Out, 0))
161
 
        G_warning(_("Failed to copy attribute table to output map"));
162
 
 
163
110
    /* Add column */
164
111
    db_init_string(&sql);
165
112
 
166
113
    /* Check if there is a database for output */
167
 
    nlinks = Vect_get_num_dblinks(&Out);
 
114
    nlinks = Vect_get_num_dblinks(&Map);
168
115
    if (nlinks < 1)
169
 
        Fi = Vect_default_field_info(&Out, 1, NULL, GV_1TABLE);
 
116
        Fi = Vect_default_field_info(&Map, layer, NULL, GV_1TABLE);
170
117
    else
171
 
        Fi = Vect_get_field(&Out, 1);
 
118
        Fi = Vect_get_field(&Map, layer);
172
119
    if (Fi == NULL) {
173
 
        Vect_delete(Out.name);
174
120
        G_fatal_error(_("Unable to get layer info for vector map <%s>"),
175
 
                      out_opt->answer);
 
121
                      map_opt->answer);
176
122
    }
177
123
 
178
124
    Driver = db_start_driver_open_database(Fi->driver, Fi->database);
179
125
    if (Driver == NULL) {
180
 
        Vect_delete(Out.name);
181
126
        G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
182
127
                      Fi->database, Fi->driver);
183
128
    }
 
129
    db_set_error_handler_driver(Driver);
184
130
 
 
131
    buf[0] = '\0';
185
132
    if (nlinks < 1)
186
133
        sprintf(buf, "create table %s (%s integer, %s integer)", Fi->table,
187
134
                Fi->key, col_opt->answer);
188
 
    else
189
 
        sprintf(buf, "alter table %s add column %s integer", Fi->table,
190
 
                col_opt->answer);
191
 
 
192
 
    db_set_string(&sql, buf);
193
 
 
194
 
    G_debug(3, "SQL: %s", db_get_string(&sql));
195
 
 
196
 
    if (db_execute_immediate(Driver, &sql) != DB_OK) {
197
 
        Vect_delete(Out.name);
198
 
        G_fatal_error(_("Cannot alter table: %s"), db_get_string(&sql));
199
 
    }
200
 
    if (nlinks < 1)
201
 
        Vect_map_add_dblink(&Out, Fi->number, Fi->name, Fi->table, Fi->key,
 
135
    else {
 
136
        dbColumn *column = NULL;
 
137
 
 
138
        /* check if column exists */
 
139
        db_get_column(Driver, Fi->table, col_opt->answer, &column);
 
140
        if (column) {
 
141
            int sqltype = db_get_column_sqltype(column);
 
142
            int ctype = db_sqltype_to_Ctype(sqltype);
 
143
 
 
144
            if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE) {
 
145
                G_fatal_error(_("Column <%s> already exists but is not numeric"),
 
146
                              col_opt->answer);
 
147
            }
 
148
        }
 
149
        else
 
150
            sprintf(buf, "alter table %s add column %s integer", Fi->table,
 
151
                    col_opt->answer);
 
152
    }
 
153
 
 
154
    if (*buf) {
 
155
        db_set_string(&sql, buf);
 
156
 
 
157
        G_debug(3, "SQL: %s", db_get_string(&sql));
 
158
 
 
159
        if (db_execute_immediate(Driver, &sql) != DB_OK) {
 
160
            G_fatal_error(_("Unable to alter table: %s"), db_get_string(&sql));
 
161
        }
 
162
    }
 
163
 
 
164
    if (nlinks < 1) {
 
165
        Vect_set_open_level(1);
 
166
        Vect_close(&Map);
 
167
        if (Vect_open_update_head(&Map, map_opt->answer, G_mapset()) < 1)
 
168
            G_fatal_error(_("Unable to modify vector map stored in other mapset"));
 
169
        Vect_map_add_dblink(&Map, layer, Fi->name, Fi->table, Fi->key,
202
170
                            Fi->database, Fi->driver);
203
 
 
204
 
    /*
205
 
     * make histogram of number sites in each test partition since the
206
 
     * number of sites will not always be a multiple of the number of
207
 
     * partitions. make_histo() returns the maximum bin size.
208
 
     */
209
 
    n = make_histo(&p, np, nsites);
210
 
 
211
 
    nlines = Vect_get_num_lines(&In);
212
 
    pnt_part = (int *)G_calloc(nlines + 1, sizeof(int));
213
 
 
214
 
    maxdist = 1.1 * sqrt(pow(box.E - box.W, 2.0) + pow(box.N - box.S, 2.0));
215
 
 
216
 
    /* Assign fold to each point */
217
 
    for (i = 0; i < np; ++i) {  /* for each partition */
218
 
        for (j = 0; j < p[i]; ++j) {    /* create p[i] random points */
219
 
            int nearest = 0;
220
 
            double dist;
221
 
 
222
 
            east = rng() / max * (box.E - box.W) + box.W;
223
 
            north = rng() / max * (box.N - box.S) + box.S;
224
 
 
225
 
            G_debug(3, "east = %f north = %f", east, north);
226
 
 
227
 
            /* find nearest */
228
 
            for (line = 1; line <= nlines; line++) {
229
 
                int type;
230
 
                double cur_dist;
231
 
 
232
 
                if (pnt_part[line] > 0)
233
 
                    continue;
234
 
 
235
 
                type = Vect_read_line(&In, Points, Cats, line);
236
 
 
237
 
                if (!(type & GV_POINT))
238
 
                    continue;
239
 
 
240
 
                cur_dist = hypot(Points->x[0] - east, Points->y[0] - north);
241
 
 
242
 
                if (nearest < 1 || cur_dist < dist) {
243
 
                    nearest = line;
244
 
                    dist = cur_dist;
245
 
                }
246
 
            }
247
 
 
248
 
            G_debug(3, "nearest = %d", nearest);
249
 
 
250
 
            /* Update */
251
 
            if (nearest > 0) {  /* shopuld be always */
252
 
                int cat;
253
 
 
254
 
                Vect_read_line(&In, Points, Cats, nearest);
255
 
                Vect_cat_get(Cats, 1, &cat);
256
 
 
257
 
                if (nlinks < 1)
258
 
                    sprintf(buf, "insert into %s (%s, %s) values (%d, %d)",
259
 
                            Fi->table, Fi->key, col_opt->answer, cat, i + 1);
260
 
                else
261
 
                    sprintf(buf, "update %s set %s = %d where %s = %d",
262
 
                            Fi->table, col_opt->answer, i + 1, Fi->key, cat);
263
 
 
264
 
                db_set_string(&sql, buf);
265
 
 
266
 
                G_debug(3, "SQL: %s", db_get_string(&sql));
267
 
 
268
 
                if (db_execute_immediate(Driver, &sql) != DB_OK) {
269
 
                    Vect_delete(Out.name);
270
 
                    G_fatal_error(_("Unable to insert row: %s"),
271
 
                                  db_get_string(&sql));
272
 
                }
273
 
                pnt_part[nearest] = i + 1;
274
 
 
275
 
            }
276
 
        }
277
 
    }
278
 
 
279
 
    Vect_close(&In);
280
 
 
 
171
        Vect_close(&Map);
 
172
 
 
173
        if (db_create_index2(Driver, Fi->table, Fi->key) !=
 
174
            DB_OK)
 
175
            G_warning(_("Cannot create index"));
 
176
 
 
177
        if (db_grant_on_table(Driver, Fi->table, DB_PRIV_SELECT,
 
178
                              DB_GROUP | DB_PUBLIC) != DB_OK)
 
179
            G_warning(_("Cannot grant privileges on table %s"),
 
180
                      Fi->table);
 
181
 
 
182
        G_important_message(_("Select privileges were granted on the table"));
 
183
 
 
184
        Vect_set_open_level(2);
 
185
        if (Vect_open_old2(&Map, map_opt->answer, G_mapset(), field_opt->answer) < 2) {
 
186
            G_fatal_error(_("Unable to open vector map <%s> at topological level %d"),
 
187
                          map_opt->answer, 2);
 
188
        }
 
189
    }
 
190
 
 
191
    /* nlines - count of records */ 
 
192
    nlines = Vect_get_num_lines(&Map);
 
193
 
 
194
    /* last_pidx - current maximal index of partition array */ 
 
195
    last_pidx = np - 1;
 
196
 
 
197
    /* minimum number of features in one partition */
 
198
    min_count = nsites / np;
 
199
    G_debug(1, "min count: %d", min_count);
 
200
 
 
201
    /* number of partitions that need min_count + 1 features */
 
202
    spill = nsites - np * min_count;  /* nsites % np, but modulus is evil */
 
203
    G_debug(1, "spill: %d", spill);
 
204
 
 
205
    /* array of available partitions */ 
 
206
    part = (struct partition *)G_calloc(np, sizeof(struct partition));
 
207
    if (!part)
 
208
        G_fatal_error(_("Out of memory"));
 
209
 
 
210
    /* line index */ 
 
211
    line_idx = (int *)G_calloc(nlines, sizeof(int));
 
212
 
 
213
    /* initialization of arrays */ 
 
214
    for (p = 0; p < np; p++) {
 
215
        part[p].id = p + 1;     /* partition ids */
 
216
        part[p].count = 0;
 
217
        part[p].max = min_count;
 
218
    }
 
219
    for (p = 0; p < spill; p++) {
 
220
        part[p].max++;
 
221
    }
 
222
 
 
223
    for (i = 0; i < nlines; i++)
 
224
        line_idx[i] = i + 1;
 
225
 
 
226
    /* proper randomization requires a 2-step randomization
 
227
     * randomize in space
 
228
     * randomize partition assignment
 
229
     * looping sequentially through all points creates a bias */
 
230
 
 
231
    /* process all points */
 
232
    db_begin_transaction(Driver);
 
233
    lines_left = nlines;
 
234
    while (lines_left) {
 
235
        int type, cat;
 
236
 
 
237
        G_percent(nlines - lines_left, nlines, 4);
 
238
 
 
239
        /* random point */
 
240
        do {
 
241
            idx = rng() * lines_left;
 
242
            /* in case rng() returns 1.0: */
 
243
        } while (idx >= lines_left);
 
244
 
 
245
        line = line_idx[idx];
 
246
        lines_left--;
 
247
        line_idx[idx] = line_idx[lines_left];
 
248
 
 
249
        if (!Vect_line_alive(&Map, line))
 
250
            continue;
 
251
 
 
252
        type = Vect_read_line(&Map, Points, Cats, line);
 
253
 
 
254
        if (!(type & GV_POINT))
 
255
            continue;
 
256
 
 
257
        /* random partition for current point */
 
258
        do {
 
259
            p = rng() * (last_pidx + 1);
 
260
            /* in case rng() returns 1.0: */
 
261
        } while (p > last_pidx);
 
262
 
 
263
        G_debug(3, "partition id = %d", part[p].id);
 
264
 
 
265
        Vect_cat_get(Cats, 1, &cat);
 
266
        if (cat < 0) {
 
267
            G_warning(_("No category for line %d in layer %d"), line, layer);
 
268
            continue;
 
269
        }
 
270
        if (nlinks < 1)
 
271
            sprintf(buf, "insert into %s (%s, %s) values (%d, %d)",
 
272
                         Fi->table, Fi->key, col_opt->answer, cat,
 
273
                         part[p].id);
 
274
        else
 
275
            sprintf(buf, "update %s set %s = %d where %s = %d", 
 
276
                         Fi->table, col_opt->answer, part[p].id,
 
277
                         Fi->key, cat);
 
278
 
 
279
        db_set_string(&sql, buf);
 
280
 
 
281
        G_debug(3, "SQL: %s", db_get_string(&sql));
 
282
 
 
283
        if (db_execute_immediate(Driver, &sql) != DB_OK) {
 
284
            G_fatal_error(_("Unable to insert row: %s"),
 
285
            db_get_string(&sql));
 
286
        }
 
287
 
 
288
        /* increase count of features in partition */
 
289
        part[p].count++;
 
290
        /* if the partition is full */
 
291
        if (part[p].count >= part[p].max) {
 
292
            /* move last partition to current position */
 
293
            if (p != last_pidx)
 
294
                part[p] = part[last_pidx];
 
295
 
 
296
            /* disable last partition */
 
297
            last_pidx--;
 
298
 
 
299
            if (last_pidx < 0 && lines_left) {
 
300
                G_fatal_error(_("internal error"));
 
301
            }
 
302
        }
 
303
    }
 
304
    G_percent(1, 1, 1);
 
305
 
 
306
    db_commit_transaction(Driver);
281
307
    db_close_database_shutdown_driver(Driver);
282
308
 
283
 
    Vect_build(&Out);
284
 
    Vect_close(&Out);
 
309
    Vect_set_db_updated(&Map);
 
310
    Vect_close(&Map);
285
311
 
286
312
    exit(EXIT_SUCCESS);
287
313
}
 
314
 
 
315