1 /***************************************************************
2 *
3 * MODULE: v.normal
4 *
5 * AUTHOR(S): James Darrell McCauley darrell@mccauley-usa.com
6 * OGR support by Martin Landa <landa.martin gmail.com> (2009)
7 *
8 * PURPOSE: GRASS program for distributional testing (based on s.normal)
9 *
10 * COPYRIGHT: (C) 2001-2014 by the GRASS Development Team
11 *
12 * This program is free software under the GNU General
13 * Public License (>=v2). Read the file COPYING that
14 * comes with GRASS for details.
15 * Modification History:
16 * <23 Jan 2001> - added field parameter, fixed reading of sites (MN)
17 * <27 Aug 1994> - began coding. Adapted cdh.f from statlib (jdm)
18 * <30 Sep 1994> - finished alpha version of cdh-c (jdm)
19 * <10 Oct 1994> - announced version 0.1B on pasture.ecn.purdue.edu (jdm)
20 * <02 Jan 1995> - cleaned man page, added html page (jdm)
21 * <25 Feb 1995> - cleaned 'gcc -Wall' warnings (jdm)
22 * <21 Jun 1995> - adapted to use new sites API (jdm)
23 * <13 Sep 2000> - released under GPL
24 *
25 **************************************************************/
26
27 #include <stdlib.h>
28 #include <math.h>
29
30 #include <grass/gis.h>
31 #include <grass/dbmi.h>
32 #include <grass/vector.h>
33 #include <grass/cdhc.h>
34 #include <grass/glocale.h>
35
36 int scan_cats(char *, long *, long *);
37
main(int argc,char ** argv)38 int main(int argc, char **argv)
39 {
40 int i, nsites, warn_once = 0;
41 int all;
42 long x, y;
43 struct Cell_head window;
44 struct GModule *module;
45 struct
46 {
47 struct Option *input, *tests, *dfield, *layer;
48 } parm;
49 struct
50 {
51 struct Flag *q, *l, *region;
52 } flag;
53 double *w, *z;
54
55 struct Map_info Map;
56 int line, nlines, npoints;
57 int field;
58 struct line_pnts *Points;
59 struct line_cats *Cats;
60 struct bound_box box;
61
62 /* Attributes */
63 int nrecords;
64 int ctype;
65 struct field_info *Fi;
66 dbDriver *Driver;
67 dbCatValArray cvarr;
68
69 G_gisinit(argv[0]);
70
71 module = G_define_module();
72 G_add_keyword(_("vector"));
73 G_add_keyword(_("statistics"));
74 G_add_keyword(_("points"));
75 G_add_keyword(_("point pattern"));
76 module->description = _("Tests for normality for vector points.");
77
78 parm.input = G_define_standard_option(G_OPT_V_MAP);
79
80 parm.layer = G_define_standard_option(G_OPT_V_FIELD);
81
82 parm.tests = G_define_option();
83 parm.tests->key = "tests";
84 parm.tests->key_desc = "range";
85 parm.tests->type = TYPE_STRING;
86 parm.tests->multiple = YES;
87 parm.tests->required = YES;
88 parm.tests->label = _("Lists of tests (1-15)");
89 parm.tests->description = _("E.g. 1,3-8,13");
90
91 parm.dfield = G_define_standard_option(G_OPT_DB_COLUMN);
92 parm.dfield->required = YES;
93
94 flag.region = G_define_flag();
95 flag.region->key = 'r';
96 flag.region->description = _("Use only points in current region");
97
98 flag.l = G_define_flag();
99 flag.l->key = 'l';
100 flag.l->description = _("Lognormality instead of normality");
101
102 if (G_parser(argc, argv))
103 exit(EXIT_FAILURE);
104
105 all = flag.region->answer ? 0 : 1;
106
107 /* Open input */
108 Vect_set_open_level(2);
109 if (Vect_open_old2(&Map, parm.input->answer, "", parm.layer->answer) < 0)
110 G_fatal_error(_("Unable to open vector map <%s>"), parm.input->answer);
111
112 field = Vect_get_field_number(&Map, parm.layer->answer);
113
114 /* Read attributes */
115 Fi = Vect_get_field(&Map, field);
116 if (Fi == NULL) {
117 G_fatal_error("Database connection not defined for layer %d", field);
118 }
119
120 Driver = db_start_driver_open_database(Fi->driver, Fi->database);
121 if (Driver == NULL)
122 G_fatal_error(_("Unable to open database <%s> by driver <%s>"),
123 Fi->database, Fi->driver);
124
125 nrecords = db_select_CatValArray(Driver, Fi->table, Fi->key, parm.dfield->answer,
126 NULL, &cvarr);
127 G_debug(1, "nrecords = %d", nrecords);
128
129 ctype = cvarr.ctype;
130 if (ctype != DB_C_TYPE_INT && ctype != DB_C_TYPE_DOUBLE)
131 G_fatal_error(_("Only numeric column type supported"));
132
133 if (nrecords < 0)
134 G_fatal_error(_("Unable to select data from table"));
135 G_verbose_message(_("%d records selected from table"), nrecords);
136
137 db_close_database_shutdown_driver(Driver);
138
139 /* Read points */
140 npoints = Vect_get_num_primitives(&Map, GV_POINT);
141 z = (double *)G_malloc(npoints * sizeof(double));
142
143 G_get_window(&window);
144 Vect_region_box(&window, &box);
145
146 Points = Vect_new_line_struct();
147 Cats = Vect_new_cats_struct();
148
149 nlines = Vect_get_num_lines(&Map);
150 nsites = 0;
151 for (line = 1; line <= nlines; line++) {
152 int type, cat, ret, cval;
153 double dval;
154
155 G_debug(3, "line = %d", line);
156
157 type = Vect_read_line(&Map, Points, Cats, line);
158 if (!(type & GV_POINT))
159 continue;
160
161 if (!all) {
162 if (!Vect_point_in_box(Points->x[0], Points->y[0], 0.0, &box))
163 continue;
164 }
165
166 Vect_cat_get(Cats, 1, &cat);
167
168 G_debug(3, "cat = %d", cat);
169
170 /* find actual value */
171 if (ctype == DB_C_TYPE_INT) {
172 ret = db_CatValArray_get_value_int(&cvarr, cat, &cval);
173 if (ret != DB_OK) {
174 G_warning(_("No record for cat %d"), cat);
175 continue;
176 }
177 dval = cval;
178 }
179 else if (ctype == DB_C_TYPE_DOUBLE) {
180 ret = db_CatValArray_get_value_double(&cvarr, cat, &dval);
181 if (ret != DB_OK) {
182 G_warning(_("No record for cat %d"), cat);
183 continue;
184 }
185 }
186
187 G_debug(3, "dval = %e", dval);
188 z[nsites] = dval;
189 nsites++;
190 }
191
192 G_verbose_message(_("Number of points: %d"), nsites);
193
194 if (nsites <= 0)
195 G_fatal_error(_("No points found"));
196
197 if (nsites < 4)
198 G_warning(_("Too small sample"));
199
200 if (flag.l->answer) {
201 warn_once = 0;
202 for (i = 0; i < nsites; ++i) {
203 if (z[i] > 1.0e-10)
204 z[i] = log10(z[i]);
205 else if (!warn_once) {
206 G_warning(_("Negative or very small point values set to -10.0"));
207 z[i] = -10.0;
208 warn_once = 1;
209 }
210 }
211 }
212
213 for (i = 0; parm.tests->answers[i]; i++)
214 if (!scan_cats(parm.tests->answers[i], &x, &y)) {
215 G_usage();
216 exit(EXIT_FAILURE);
217 }
218 for (i = 0; parm.tests->answers[i]; i++) {
219 scan_cats(parm.tests->answers[i], &x, &y);
220 while (x <= y)
221 switch (x++) {
222 case 1: /* moments */
223 fprintf(stdout, _("Moments \\sqrt{b_1} and b_2: "));
224 w = Cdhc_omnibus_moments(z, nsites);
225 fprintf(stdout, "%g %g\n", w[0], w[1]);
226 break;
227 case 2: /* geary */
228 fprintf(stdout, _("Geary's a-statistic & an approx. normal: "));
229 w = Cdhc_geary_test(z, nsites);
230 fprintf(stdout, "%g %g\n", w[0], w[1]);
231 break;
232 case 3: /* extreme deviates */
233 fprintf(stdout, _("Extreme normal deviates: "));
234 w = Cdhc_extreme(z, nsites);
235 fprintf(stdout, "%g %g\n", w[0], w[1]);
236 break;
237 case 4: /* D'Agostino */
238 fprintf(stdout, _("D'Agostino's D & an approx. normal: "));
239 w = Cdhc_dagostino_d(z, nsites);
240 fprintf(stdout, "%g %g\n", w[0], w[1]);
241 break;
242 case 5: /* Kuiper */
243 fprintf(stdout,
244 _("Kuiper's V (regular & modified for normality): "));
245 w = Cdhc_kuipers_v(z, nsites);
246 fprintf(stdout, "%g %g\n", w[1], w[0]);
247 break;
248 case 6: /* Watson */
249 fprintf(stdout,
250 _("Watson's U^2 (regular & modified for normality): "));
251 w = Cdhc_watson_u2(z, nsites);
252 fprintf(stdout, "%g %g\n", w[1], w[0]);
253 break;
254 case 7: /* Durbin */
255 fprintf(stdout,
256 _("Durbin's Exact Test (modified Kolmogorov): "));
257 w = Cdhc_durbins_exact(z, nsites);
258 fprintf(stdout, "%g\n", w[0]);
259 break;
260 case 8: /* Anderson-Darling */
261 fprintf(stdout,
262 _("Anderson-Darling's A^2 (regular & modified for normality): "));
263 w = Cdhc_anderson_darling(z, nsites);
264 fprintf(stdout, "%g %g\n", w[1], w[0]);
265 break;
266 case 9: /* Cramer-Von Mises */
267 fprintf(stdout,
268 _("Cramer-Von Mises W^2(regular & modified for normality): "));
269 w = Cdhc_cramer_von_mises(z, nsites);
270 fprintf(stdout, "%g %g\n", w[1], w[0]);
271 break;
272 case 10: /* Kolmogorov-Smirnov */
273 fprintf(stdout,
274 _("Kolmogorov-Smirnov's D (regular & modified for normality): "));
275 w = Cdhc_kolmogorov_smirnov(z, nsites);
276 fprintf(stdout, "%g %g\n", w[1], w[0]);
277 break;
278 case 11: /* chi-square */
279 fprintf(stdout,
280 _("Chi-Square stat (equal probability classes) and d.f.: "));
281 w = Cdhc_chi_square(z, nsites);
282 fprintf(stdout, "%g %d\n", w[0], (int)w[1]);
283 break;
284 case 12: /* Shapiro-Wilk */
285 if (nsites > 50) {
286 G_warning(_("Shapiro-Wilk's W cannot be used for n > 50"));
287 if (nsites < 99)
288 G_message(_("Use Weisberg-Binghams's W''"));
289 }
290 else {
291 fprintf(stdout, _("Shapiro-Wilk W: "));
292 w = Cdhc_shapiro_wilk(z, nsites);
293 fprintf(stdout, "%g\n", w[0]);
294 }
295 break;
296 case 13: /* Weisberg-Bingham */
297 if (nsites > 99 || nsites < 50)
298 G_warning(_("Weisberg-Bingham's W'' cannot be used for n < 50 or n > 99"));
299 else {
300 fprintf(stdout, _("Weisberg-Bingham's W'': "));
301 w = Cdhc_weisberg_bingham(z, nsites);
302 fprintf(stdout, "%g\n", w[0]);
303 }
304 break;
305 case 14: /* Royston */
306 if (nsites > 2000)
307 G_warning(_("Royston only extended Shapiro-Wilk's W up to n = 2000"));
308 else {
309 fprintf(stdout, _("Shapiro-Wilk W'': "));
310 w = Cdhc_royston(z, nsites);
311 fprintf(stdout, "%g\n", w[0]);
312 }
313 break;
314 case 15: /* Kotz */
315 fprintf(stdout, _("Kotz' T'_f (Lognormality vs. Normality): "));
316 w = Cdhc_kotz_families(z, nsites);
317 fprintf(stdout, "%g\n", w[0]);
318 break;
319 default:
320 break;
321 }
322 }
323 exit(EXIT_SUCCESS);
324 }
325