17e03130923bab5888580f0493d39f8decd712e3
[pspp] / src / language / stats / frequencies.c
1 /*
2   PSPP - a program for statistical analysis.
3   Copyright (C) 1997-9, 2000, 2007, 2009, 2010, 2011, 2014 Free Software Foundation, Inc.
4    
5   This program is free software: you can redistribute it and/or modify
6   it under the terms of the GNU General Public License as published by
7   the Free Software Foundation, either version 3 of the License, or
8   (at your option) any later version.
9
10   This program is distributed in the hope that it will be useful,
11   but WITHOUT ANY WARRANTY; without even the implied warranty of
12   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13   GNU General Public License for more details.
14   
15   You should have received a copy of the GNU General Public License
16   along with this program.  If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 #include <config.h>
20 #include <stdlib.h>
21 #include <gsl/gsl_histogram.h>
22
23
24 #include "data/case.h"
25 #include "data/casegrouper.h"
26 #include "data/casereader.h"
27 #include "data/dataset.h"
28 #include "data/dictionary.h"
29 #include "data/format.h"
30 #include "data/missing-values.h"
31 #include "data/settings.h"
32 #include "data/value-labels.h"
33 #include "data/variable.h"
34
35 #include "language/dictionary/split-file.h"
36
37 #include "language/command.h"
38 #include "language/lexer/lexer.h"
39 #include "language/lexer/variable-parser.h"
40 #include "language/stats/freq.h"
41
42 #include "libpspp/array.h"
43 #include "libpspp/bit-vector.h"
44 #include "libpspp/compiler.h"
45 #include "libpspp/hmap.h"
46 #include "libpspp/message.h"
47 #include "libpspp/misc.h"
48 #include "libpspp/pool.h"
49
50 #include "math/histogram.h"
51 #include "math/moments.h"
52 #include "math/chart-geometry.h"
53
54
55 #include "output/chart-item.h"
56 #include "output/charts/piechart.h"
57 #include "output/charts/plot-hist.h"
58 #include "output/tab.h"
59
60 #include "gl/minmax.h"
61 #include "gl/xalloc.h"
62
63 #include "gettext.h"
64 #define _(msgid) gettext (msgid)
65 #define N_(msgid) msgid
66
67 /* Percentiles to calculate. */
68
69 struct percentile
70 {
71   double p;        /* the %ile to be calculated */
72   double value;    /* the %ile's value */
73   bool show;       /* True to show this percentile in the statistics box. */
74 };
75
76 static int
77 ptile_3way (const void *_p1, const void *_p2)
78 {
79   const struct percentile *p1 = _p1;
80   const struct percentile *p2 = _p2;
81
82   if (p1->p < p2->p)
83     return -1;
84
85   if (p1->p == p2->p)
86     {
87       if (p1->show > p2->show)
88         return -1;
89
90       return (p1->show < p2->show);
91     }
92
93   return (p1->p > p2->p);
94 }
95
96
97 enum
98   {
99     FRQ_NONORMAL,
100     FRQ_NORMAL
101   };
102
103 enum
104   {
105     FRQ_FREQ,
106     FRQ_PERCENT
107   };
108
109 enum sortprops 
110   {
111     FRQ_AFREQ,
112     FRQ_DFREQ,
113     FRQ_AVALUE,
114     FRQ_DVALUE
115   };
116
117 /* Array indices for STATISTICS subcommand. */
118 enum
119   {
120     FRQ_ST_MEAN,
121     FRQ_ST_SEMEAN,
122     FRQ_ST_MEDIAN,
123     FRQ_ST_MODE,
124     FRQ_ST_STDDEV,
125     FRQ_ST_VARIANCE,
126     FRQ_ST_KURTOSIS,
127     FRQ_ST_SEKURTOSIS,
128     FRQ_ST_SKEWNESS,
129     FRQ_ST_SESKEWNESS,
130     FRQ_ST_RANGE,
131     FRQ_ST_MINIMUM,
132     FRQ_ST_MAXIMUM,
133     FRQ_ST_SUM,
134     FRQ_ST_count
135   };
136
137 /* Description of statistics. */
138 static const char *st_name[FRQ_ST_count] =
139 {
140    N_("Mean"),
141    N_("S.E. Mean"),
142    N_("Median"),
143    N_("Mode"),
144    N_("Std Dev"),
145    N_("Variance"),
146    N_("Kurtosis"),
147    N_("S.E. Kurt"),
148    N_("Skewness"),
149    N_("S.E. Skew"),
150    N_("Range"),
151    N_("Minimum"),
152    N_("Maximum"),
153    N_("Sum")
154 };
155
156 struct freq_tab
157   {
158     struct hmap data;           /* Hash table for accumulating counts. */
159     struct freq *valid;         /* Valid freqs. */
160     int n_valid;                /* Number of total freqs. */
161     const struct dictionary *dict; /* Source of entries in the table. */
162
163     struct freq *missing;       /* Missing freqs. */
164     int n_missing;              /* Number of missing freqs. */
165
166     /* Statistics. */
167     double total_cases;         /* Sum of weights of all cases. */
168     double valid_cases;         /* Sum of weights of valid cases. */
169   };
170
171 struct frq_chart
172   {
173     double x_min;               /* X axis minimum value. */
174     double x_max;               /* X axis maximum value. */
175     int y_scale;                /* Y axis scale: FRQ_FREQ or FRQ_PERCENT. */
176
177     /* Histograms only. */
178     double y_max;               /* Y axis maximum value. */
179     bool draw_normal;           /* Whether to draw normal curve. */
180
181     /* Pie charts only. */
182     bool include_missing;       /* Whether to include missing values. */
183   };
184
185 /* Per-variable frequency data. */
186 struct var_freqs
187   {
188     const struct variable *var;
189
190     /* Freqency table. */
191     struct freq_tab tab;        /* Frequencies table to use. */
192
193     /* Percentiles. */
194     int n_groups;               /* Number of groups. */
195     double *groups;             /* Groups. */
196
197     /* Statistics. */
198     double stat[FRQ_ST_count];
199
200     /* Variable attributes. */
201     int width;
202   };
203
204 struct frq_proc
205   {
206     struct pool *pool;
207
208     struct var_freqs *vars;
209     size_t n_vars;
210
211     /* Percentiles to calculate and possibly display. */
212     struct percentile *percentiles;
213     int n_percentiles, n_show_percentiles;
214
215     /* Frequency table display. */
216     int max_categories;         /* Maximum categories to show. */
217     int sort;                   /* FRQ_AVALUE or FRQ_DVALUE
218                                    or FRQ_AFREQ or FRQ_DFREQ. */
219
220     /* Statistics; number of statistics. */
221     unsigned long stats;
222     int n_stats;
223
224     /* Histogram and pie chart settings. */
225     struct frq_chart *hist, *pie;
226   };
227
228
229 struct freq_compare_aux
230   {
231     bool by_freq;
232     bool ascending_freq;
233
234     int width;
235     bool ascending_value;
236   };
237
238 static void calc_stats (const struct var_freqs *vf, double d[FRQ_ST_count]);
239
240 static void do_piechart(const struct frq_chart *pie, 
241                         const struct variable *var,
242                         const struct freq_tab *frq_tab);
243
244 static void dump_statistics (const struct frq_proc *frq, 
245                              const struct var_freqs *vf,
246                              const struct variable *wv);
247
248 static int
249 compare_freq (const void *a_, const void *b_, const void *aux_)
250 {
251   const struct freq_compare_aux *aux = aux_;
252   const struct freq *a = a_;
253   const struct freq *b = b_;
254
255   if (aux->by_freq && a->count != b->count)
256     {
257       int cmp = a->count > b->count ? 1 : -1;
258       return aux->ascending_freq ? cmp : -cmp;
259     }
260   else
261     {
262       int cmp = value_compare_3way (&a->value, &b->value, aux->width);
263       return aux->ascending_value ? cmp : -cmp;
264     }
265 }
266
267 /* Create a gsl_histogram from a freq_tab */
268 static struct histogram *
269 freq_tab_to_hist (const struct frq_proc *frq, const struct freq_tab *ft,
270                   const struct variable *var);
271
272
273 /* Displays a full frequency table for variable V. */
274 static void
275 dump_freq_table (const struct var_freqs *vf, const struct variable *wv)
276 {
277   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
278   const struct freq_tab *ft = &vf->tab;
279   int n_categories;
280   struct freq *f;
281   struct tab_table *t;
282   int r, x;
283   double cum_total = 0.0;
284   double cum_freq = 0.0;
285
286   static const char *headings[] = {
287     N_("Value Label"),
288     N_("Value"),
289     N_("Frequency"),
290     N_("Percent"),
291     N_("Valid Percent"),
292     N_("Cum Percent")
293   };
294
295   n_categories = ft->n_valid + ft->n_missing;
296   t = tab_create (6, n_categories + 2);
297   tab_set_format (t, RC_WEIGHT, wfmt);
298   tab_headers (t, 0, 0, 1, 0);
299
300   for (x = 0; x < 6; x++)
301     tab_text (t, x, 0, TAB_CENTER | TAT_TITLE, gettext (headings[x]));
302
303   r = 1;
304   for (f = ft->valid; f < ft->missing; f++)
305     {
306       const char *label;
307       double percent, valid_percent;
308
309       cum_freq += f->count;
310
311       percent = f->count / ft->total_cases * 100.0;
312       valid_percent = f->count / ft->valid_cases * 100.0;
313       cum_total += valid_percent;
314
315       label = var_lookup_value_label (vf->var, &f->value);
316       if (label != NULL)
317         tab_text (t, 0, r, TAB_LEFT, label);
318
319       tab_value (t, 1, r, TAB_NONE, &f->value, vf->var, NULL);
320       tab_double (t, 2, r, TAB_NONE, f->count, NULL, RC_WEIGHT);
321       tab_double (t, 3, r, TAB_NONE, percent, NULL, RC_OTHER);
322       tab_double (t, 4, r, TAB_NONE, valid_percent, NULL, RC_OTHER);
323       tab_double (t, 5, r, TAB_NONE, cum_total, NULL, RC_OTHER);
324       r++;
325     }
326   for (; f < &ft->valid[n_categories]; f++)
327     {
328       const char *label;
329
330       cum_freq += f->count;
331
332       label = var_lookup_value_label (vf->var, &f->value);
333       if (label != NULL)
334         tab_text (t, 0, r, TAB_LEFT, label);
335
336       tab_value (t, 1, r, TAB_NONE, &f->value, vf->var, NULL);
337       tab_double (t, 2, r, TAB_NONE, f->count, NULL, RC_WEIGHT);
338       tab_double (t, 3, r, TAB_NONE,
339                   f->count / ft->total_cases * 100.0, NULL, RC_OTHER);
340       tab_text (t, 4, r, TAB_NONE, _("Missing"));
341       r++;
342     }
343
344   tab_box (t, TAL_1, TAL_1, -1, TAL_1, 0, 0, 5, r);
345   tab_hline (t, TAL_2, 0, 5, 1);
346   tab_hline (t, TAL_2, 0, 5, r);
347   tab_joint_text (t, 0, r, 1, r, TAB_RIGHT | TAT_TITLE, _("Total"));
348   tab_vline (t, TAL_0, 1, r, r);
349   tab_double (t, 2, r, TAB_NONE, cum_freq, NULL, RC_WEIGHT);
350   tab_double (t, 3, r, TAB_NONE, 100.0, &F_5_1, RC_OTHER);
351   tab_double (t, 4, r, TAB_NONE, 100.0, &F_5_1, RC_OTHER);
352
353   tab_title (t, "%s", var_to_string (vf->var));
354   tab_submit (t);
355 }
356 \f
357 /* Statistical display. */
358
359 static double
360 calc_percentile (double p, double valid_cases, double x1, double x2)
361 {
362   double s, dummy;
363
364   s = (settings_get_algorithm () != COMPATIBLE
365        ? modf ((valid_cases - 1) * p, &dummy)
366        : modf ((valid_cases + 1) * p - 1, &dummy));
367
368   return x1 + (x2 - x1) * s;
369 }
370
371 /* Calculates all of the percentiles for VF within FRQ. */
372 static void
373 calc_percentiles (const struct frq_proc *frq, const struct var_freqs *vf)
374 {
375   const struct freq_tab *ft = &vf->tab;
376   double W = ft->valid_cases;
377   const struct freq *f;
378   int percentile_idx = 0;
379   double  rank = 0;
380
381   for (f = ft->valid; f < ft->missing; f++)
382     {
383       rank += f->count;
384       for (; percentile_idx < frq->n_percentiles; percentile_idx++)
385         {
386           struct percentile *pc = &frq->percentiles[percentile_idx];
387           double tp;
388
389           tp = (settings_get_algorithm () == ENHANCED
390                 ? (W - 1) * pc->p
391                 : (W + 1) * pc->p - 1);
392
393           if (rank <= tp)
394             break;
395
396           if (tp + 1 < rank || f + 1 >= ft->missing)
397             pc->value = f->value.f;
398           else
399             pc->value = calc_percentile (pc->p, W, f->value.f, f[1].value.f);
400         }
401     }
402   for (; percentile_idx < frq->n_percentiles; percentile_idx++)
403     {
404       struct percentile *pc = &frq->percentiles[percentile_idx];
405       pc->value = ft->valid[ft->n_valid - 1].value.f;
406     }
407 }
408
409 /* Returns true iff the value in struct freq F is non-missing
410    for variable V. */
411 static bool
412 not_missing (const void *f_, const void *v_)
413 {
414   const struct freq *f = f_;
415   const struct variable *v = v_;
416
417   return !var_is_value_missing (v, &f->value, MV_ANY);
418 }
419
420
421 /* Summarizes the frequency table data for variable V. */
422 static void
423 postprocess_freq_tab (const struct frq_proc *frq, struct var_freqs *vf)
424 {
425   struct freq_tab *ft = &vf->tab;
426   struct freq_compare_aux aux;
427   size_t count;
428   struct freq *freqs, *f;
429   size_t i;
430
431   /* Extract data from hash table. */
432   count = hmap_count (&ft->data);
433   freqs = freq_hmap_extract (&ft->data);
434
435   /* Put data into ft. */
436   ft->valid = freqs;
437   ft->n_valid = partition (freqs, count, sizeof *freqs, not_missing, vf->var);
438   ft->missing = freqs + ft->n_valid;
439   ft->n_missing = count - ft->n_valid;
440
441   /* Sort data. */
442   aux.by_freq = frq->sort == FRQ_AFREQ || frq->sort == FRQ_DFREQ;
443   aux.ascending_freq = frq->sort != FRQ_DFREQ;
444   aux.width = vf->width;
445   aux.ascending_value = frq->sort != FRQ_DVALUE;
446   sort (ft->valid, ft->n_valid, sizeof *ft->valid, compare_freq, &aux);
447   sort (ft->missing, ft->n_missing, sizeof *ft->missing, compare_freq, &aux);
448
449   /* Summary statistics. */
450   ft->valid_cases = 0.0;
451   for(i = 0 ;  i < ft->n_valid ; ++i )
452     {
453       f = &ft->valid[i];
454       ft->valid_cases += f->count;
455
456     }
457
458   ft->total_cases = ft->valid_cases ;
459   for(i = 0 ;  i < ft->n_missing ; ++i )
460     {
461       f = &ft->missing[i];
462       ft->total_cases += f->count;
463     }
464
465 }
466
467 /* Frees the frequency table for variable V. */
468 static void
469 cleanup_freq_tab (struct var_freqs *vf)
470 {
471   free (vf->tab.valid);
472   freq_hmap_destroy (&vf->tab.data, vf->width);
473 }
474
475 /* Add data from case C to the frequency table. */
476 static void
477 calc (struct frq_proc *frq, const struct ccase *c, const struct dataset *ds)
478 {
479   double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
480   size_t i;
481
482   for (i = 0; i < frq->n_vars; i++)
483     {
484       struct var_freqs *vf = &frq->vars[i];
485       const union value *value = case_data (c, vf->var);
486       size_t hash = value_hash (value, vf->width, 0);
487       struct freq *f;
488
489       f = freq_hmap_search (&vf->tab.data, value, vf->width, hash);
490       if (f == NULL)
491         f = freq_hmap_insert (&vf->tab.data, value, vf->width, hash);
492
493       f->count += weight;
494     }
495 }
496
497 /* Prepares each variable that is the target of FREQUENCIES by setting
498    up its hash table. */
499 static void
500 precalc (struct frq_proc *frq, struct casereader *input, struct dataset *ds)
501 {
502   struct ccase *c;
503   size_t i;
504
505   c = casereader_peek (input, 0);
506   if (c != NULL)
507     {
508       output_split_file_values (ds, c);
509       case_unref (c);
510     }
511
512   for (i = 0; i < frq->n_vars; i++)
513     hmap_init (&frq->vars[i].tab.data);
514 }
515
516 /* Finishes up with the variables after frequencies have been
517    calculated.  Displays statistics, percentiles, ... */
518 static void
519 postcalc (struct frq_proc *frq, const struct dataset *ds)
520 {
521   const struct dictionary *dict = dataset_dict (ds);
522   const struct variable *wv = dict_get_weight (dict);
523   size_t i;
524
525   for (i = 0; i < frq->n_vars; i++)
526     {
527       struct var_freqs *vf = &frq->vars[i];
528
529       postprocess_freq_tab (frq, vf);
530
531       /* Frequencies tables. */
532       if (vf->tab.n_valid + vf->tab.n_missing <= frq->max_categories)
533         dump_freq_table (vf, wv);
534
535       calc_percentiles (frq, vf);
536
537       /* Statistics. */
538       if (frq->n_stats)
539         dump_statistics (frq, vf, wv);
540
541       if (frq->hist && var_is_numeric (vf->var) && vf->tab.n_valid > 0)
542         {
543           double d[FRQ_ST_count];
544           struct histogram *histogram;
545
546           calc_stats (vf, d);
547
548           histogram = freq_tab_to_hist (frq, &vf->tab, vf->var);
549
550           if ( histogram)
551             {
552               chart_item_submit (histogram_chart_create (
553                                histogram->gsl_hist, var_to_string(vf->var),
554                                vf->tab.valid_cases,
555                                d[FRQ_ST_MEAN],
556                                d[FRQ_ST_STDDEV],
557                                frq->hist->draw_normal));
558
559               statistic_destroy (&histogram->parent);
560             }
561         }
562
563       if (frq->pie)
564         do_piechart(frq->pie, vf->var, &vf->tab);
565
566       cleanup_freq_tab (vf);
567     }
568 }
569
570 int
571 cmd_frequencies (struct lexer *lexer, struct dataset *ds)
572 {
573   int i;
574   struct frq_proc frq;
575   const struct variable **vars;
576
577   bool sbc_barchart = false;
578   bool sbc_piechart = false;
579   bool sbc_histogram = false;
580
581   double pie_min = -DBL_MAX;
582   double pie_max = DBL_MAX;
583   bool pie_missing = false;
584
585   double hi_min = -DBL_MAX;
586   double hi_max = DBL_MAX;
587   int hi_scale = FRQ_FREQ;
588   int hi_freq = INT_MIN;
589   int hi_pcnt = INT_MIN;
590   int hi_norm = FRQ_NONORMAL;
591
592   frq.pool = pool_create ();
593   frq.sort = FRQ_AVALUE;
594
595   frq.vars = NULL;
596   frq.n_vars = 0;
597   
598   frq.stats = BIT_INDEX (FRQ_ST_MEAN) 
599     | BIT_INDEX (FRQ_ST_STDDEV) 
600     | BIT_INDEX (FRQ_ST_MINIMUM)
601     | BIT_INDEX (FRQ_ST_MAXIMUM);
602
603   frq.n_stats = 4;
604
605   frq.max_categories = INT_MAX;
606
607   frq.percentiles = NULL;
608   frq.n_percentiles = 0;
609   frq.n_show_percentiles = 0;
610
611   frq.hist = NULL;
612   frq.pie = NULL;
613
614
615   /* Accept an optional, completely pointless "/VARIABLES=" */
616   lex_match (lexer, T_SLASH);
617   if (lex_match_id  (lexer, "VARIABLES"))
618     {
619       if (! lex_force_match (lexer, T_EQUALS) )
620         goto error;
621     }
622
623   if (!parse_variables_const (lexer, dataset_dict (ds),
624                               &vars,
625                               &frq.n_vars,
626                               PV_NO_DUPLICATE))
627     goto error;
628
629   frq.vars = xzalloc (frq.n_vars * sizeof (*frq.vars));
630   for (i = 0; i < frq.n_vars; ++i)
631     {
632       frq.vars[i].var = vars[i];
633       frq.vars[i].width = var_get_width (vars[i]);
634     }
635
636   while (lex_token (lexer) != T_ENDCMD)
637     {
638       lex_match (lexer, T_SLASH);
639
640       if (lex_match_id (lexer, "STATISTICS"))
641         {
642           frq.stats = BIT_INDEX (FRQ_ST_MEAN) 
643             | BIT_INDEX (FRQ_ST_STDDEV) 
644             | BIT_INDEX (FRQ_ST_MINIMUM)
645             | BIT_INDEX (FRQ_ST_MAXIMUM);
646           
647           frq.n_stats = 4;
648
649           if (lex_match (lexer, T_EQUALS))
650             {
651               frq.n_stats = 0;
652               frq.stats = 0;
653             }
654
655           while (lex_token (lexer) != T_ENDCMD
656                  && lex_token (lexer) != T_SLASH)
657             {
658               if (lex_match_id (lexer, "DEFAULT"))
659                 {
660                   frq.stats = BIT_INDEX (FRQ_ST_MEAN) 
661                     | BIT_INDEX (FRQ_ST_STDDEV) 
662                     | BIT_INDEX (FRQ_ST_MINIMUM)
663                     | BIT_INDEX (FRQ_ST_MAXIMUM);
664
665                   frq.n_stats = 4;
666                 }
667               else if (lex_match_id (lexer, "MEAN"))
668                 {
669                   frq.stats |= BIT_INDEX (FRQ_ST_MEAN);
670                   frq.n_stats++;
671                 }
672               else if (lex_match_id (lexer, "SEMEAN"))
673                 {
674                   frq.stats |= BIT_INDEX (FRQ_ST_SEMEAN);
675                   frq.n_stats++;
676                 }
677               else if (lex_match_id (lexer, "MEDIAN"))
678                 {
679                   frq.stats |= BIT_INDEX (FRQ_ST_MEDIAN);
680                   frq.n_stats++;
681                 }
682               else if (lex_match_id (lexer, "MODE"))
683                 {
684                   frq.stats |= BIT_INDEX (FRQ_ST_MODE);
685                   frq.n_stats++;
686                 }
687               else if (lex_match_id (lexer, "STDDEV"))
688                 {
689                   frq.stats |= BIT_INDEX (FRQ_ST_STDDEV);
690                   frq.n_stats++;
691                 }
692               else if (lex_match_id (lexer, "VARIANCE"))
693                 {
694                   frq.stats |= BIT_INDEX (FRQ_ST_MEAN);
695                   frq.n_stats++;
696                 }
697               else if (lex_match_id (lexer, "KURTOSIS"))
698                 {
699                   frq.stats |= BIT_INDEX (FRQ_ST_KURTOSIS);
700                   frq.n_stats++;
701                 }
702               else if (lex_match_id (lexer, "SKEWNESS"))
703                 {
704                   frq.stats |= BIT_INDEX (FRQ_ST_SKEWNESS);
705                   frq.n_stats++;
706                 }
707               else if (lex_match_id (lexer, "RANGE"))
708                 {
709                   frq.stats |= BIT_INDEX (FRQ_ST_RANGE);
710                   frq.n_stats++;
711                 }
712               else if (lex_match_id (lexer, "MINIMUM"))
713                 {
714                   frq.stats |= BIT_INDEX (FRQ_ST_MINIMUM);
715                   frq.n_stats++;
716                 }
717               else if (lex_match_id (lexer, "MAXIMUM"))
718                 {
719                   frq.stats |= BIT_INDEX (FRQ_ST_MAXIMUM);
720                   frq.n_stats++;
721                 }
722               else if (lex_match_id (lexer, "SUM"))
723                 {
724                   frq.stats |= BIT_INDEX (FRQ_ST_SUM);
725                   frq.n_stats++;
726                 }
727               else if (lex_match_id (lexer, "SESKEWNESS"))
728                 {
729                   frq.stats |= BIT_INDEX (FRQ_ST_SESKEWNESS);
730                   frq.n_stats++;
731                 }
732               else if (lex_match_id (lexer, "SEKURTOSIS"))
733                 {
734                   frq.stats |= BIT_INDEX (FRQ_ST_SEKURTOSIS);
735                   frq.n_stats++;
736                 }
737               else if (lex_match_id (lexer, "NONE"))
738                 {
739                   frq.stats = 0;
740                   frq.n_stats = 0;
741                 }
742               else if (lex_match (lexer, T_ALL))
743                 {
744                   frq.stats = ~0;
745                   frq.n_stats = FRQ_ST_count;
746                 }
747               else
748                 {
749                   lex_error (lexer, NULL);
750                   goto error;
751                 }
752             }
753         }
754       else if (lex_match_id (lexer, "PERCENTILES"))
755         {
756           lex_match (lexer, T_EQUALS);
757           while (lex_token (lexer) != T_ENDCMD
758                  && lex_token (lexer) != T_SLASH)
759             {
760               if (lex_force_num (lexer))
761                 {
762                   frq.percentiles =
763                     xrealloc (frq.percentiles, 
764                               (frq.n_percentiles + 1)
765                               * sizeof (*frq.percentiles));
766                   frq.percentiles[frq.n_percentiles].p = lex_number (lexer)  / 100.0;
767                   frq.percentiles[frq.n_percentiles].show = true;
768                   lex_get (lexer);
769                   frq.n_percentiles++;
770                   frq.n_show_percentiles++;
771                 }
772               else
773                 {
774                   lex_error (lexer, NULL);
775                   goto error;
776                 }
777               lex_match (lexer, T_COMMA);
778             }
779         }
780       else if (lex_match_id (lexer, "FORMAT"))
781         {
782           lex_match (lexer, T_EQUALS);
783           while (lex_token (lexer) != T_ENDCMD
784                  && lex_token (lexer) != T_SLASH)
785             {
786               if (lex_match_id (lexer, "TABLE"))
787                 {
788                   
789                 }
790               else if (lex_match_id (lexer, "NOTABLE"))
791                 {
792                   frq.max_categories = 0;
793                 }
794               else if (lex_match_id (lexer, "AVALUE"))
795                 {
796                   frq.sort = FRQ_AVALUE;
797                 }
798               else if (lex_match_id (lexer, "DVALUE"))
799                 {
800                   frq.sort = FRQ_DVALUE;
801                 }
802               else if (lex_match_id (lexer, "AFREQ"))
803                 {
804                   frq.sort = FRQ_AFREQ;
805                 }
806               else if (lex_match_id (lexer, "DFREQ"))
807                 {
808                   frq.sort = FRQ_DFREQ;
809                 }
810               else
811                 {
812                   lex_error (lexer, NULL);
813                   goto error;
814                 }
815             }
816         }
817       else if (lex_match_id (lexer, "NTILES"))
818         {
819           lex_match (lexer, T_EQUALS);
820
821           if (lex_force_int (lexer))
822             {
823               int i;
824               int n = lex_integer (lexer);
825               lex_get (lexer);
826               for (i = 0; i < n + 1; ++i)
827                 {
828                   frq.percentiles =
829                     xrealloc (frq.percentiles, 
830                               (frq.n_percentiles + 1)
831                               * sizeof (*frq.percentiles));
832                   frq.percentiles[frq.n_percentiles].p =
833                     i / (double) n ;
834                   frq.percentiles[frq.n_percentiles].show = true;
835
836                   frq.n_percentiles++;
837                   frq.n_show_percentiles++;
838                 }
839             }
840           else
841             {
842               lex_error (lexer, NULL);
843               goto error;
844             }
845         }
846       else if (lex_match_id (lexer, "ALGORITHM"))
847         {
848           lex_match (lexer, T_EQUALS);
849
850           if (lex_match_id (lexer, "COMPATIBLE"))
851             {
852               settings_set_cmd_algorithm (COMPATIBLE);
853             }
854           else if (lex_match_id (lexer, "ENHANCED"))
855             {
856               settings_set_cmd_algorithm (ENHANCED);
857             }
858           else
859             {
860               lex_error (lexer, NULL);
861               goto error;
862             }
863         }
864       else if (lex_match_id (lexer, "HISTOGRAM"))
865         {
866           lex_match (lexer, T_EQUALS);
867           sbc_histogram = true;
868
869           while (lex_token (lexer) != T_ENDCMD
870                  && lex_token (lexer) != T_SLASH)
871             {
872               if (lex_match_id (lexer, "NORMAL"))
873                 {
874                   hi_norm = FRQ_NORMAL;
875                 }
876               else if (lex_match_id (lexer, "NONORMAL"))
877                 {
878                   hi_norm = FRQ_NONORMAL;
879                 }
880               else if (lex_match_id (lexer, "FREQ"))
881                 {
882                   hi_scale = FRQ_FREQ;
883                   if (lex_match (lexer, T_LPAREN))
884                     {
885                       if (lex_force_int (lexer))
886                         {
887                           hi_freq = lex_integer (lexer);
888                           if (hi_freq <= 0)
889                             {
890                               lex_error (lexer, _("Histogram frequency must be greater than zero."));
891                             }
892                           lex_get (lexer);
893                           lex_force_match (lexer, T_RPAREN);
894                         }
895                     }
896                 }
897               else if (lex_match_id (lexer, "PERCENT"))
898                 {
899                   hi_scale = FRQ_PERCENT;
900                   if (lex_match (lexer, T_LPAREN))
901                     {
902                       if (lex_force_int (lexer))
903                         {
904                           hi_pcnt = lex_integer (lexer);
905                           if (hi_pcnt <= 0)
906                             {
907                               lex_error (lexer, _("Histogram percentage must be greater than zero."));
908                             }
909                           lex_get (lexer);
910                           lex_force_match (lexer, T_RPAREN);
911                         }
912                     }
913                 }
914               else if (lex_match_id (lexer, "MINIMUM"))
915                 {
916                   lex_force_match (lexer, T_LPAREN);
917                   if (lex_force_num (lexer))
918                     {
919                       hi_min = lex_number (lexer);
920                       lex_get (lexer);
921                     }
922                   lex_force_match (lexer, T_RPAREN);
923                 }
924               else if (lex_match_id (lexer, "MAXIMUM"))
925                 {
926                   lex_force_match (lexer, T_LPAREN);
927                   if (lex_force_num (lexer))
928                     {
929                       hi_max = lex_number (lexer);
930                       lex_get (lexer);
931                     }
932                   lex_force_match (lexer, T_RPAREN);
933                 }
934               else
935                 {
936                   lex_error (lexer, NULL);
937                   goto error;
938                 }
939             }
940         }
941       else if (lex_match_id (lexer, "PIECHART"))
942         {
943           lex_match (lexer, T_EQUALS);
944           while (lex_token (lexer) != T_ENDCMD
945                  && lex_token (lexer) != T_SLASH)
946             {
947               if (lex_match_id (lexer, "MINIMUM"))
948                 {
949                   lex_force_match (lexer, T_LPAREN);
950                   if (lex_force_num (lexer))
951                     {
952                       pie_min = lex_number (lexer);
953                       lex_get (lexer);
954                     }
955                   lex_force_match (lexer, T_RPAREN);
956                 }
957               else if (lex_match_id (lexer, "MAXIMUM"))
958                 {
959                   lex_force_match (lexer, T_LPAREN);
960                   if (lex_force_num (lexer))
961                     {
962                       pie_max = lex_number (lexer);
963                       lex_get (lexer);
964                     }
965                   lex_force_match (lexer, T_RPAREN);
966                 }
967               else if (lex_match_id (lexer, "MISSING"))
968                 {
969                   pie_missing = true;
970                 }
971               else if (lex_match_id (lexer, "NOMISSING"))
972                 {
973                   pie_missing = false;
974                 }
975               else
976                 {
977                   lex_error (lexer, NULL);
978                   goto error;
979                 }
980             }
981           sbc_piechart = true;
982         }
983       else if (lex_match_id (lexer, "MISSING"))
984         {
985           lex_match (lexer, T_EQUALS);
986
987           while (lex_token (lexer) != T_ENDCMD
988                  && lex_token (lexer) != T_SLASH)
989             {
990               if (lex_match_id (lexer, "EXCLUDE"))
991                 {
992                 }
993               else if (lex_match_id (lexer, "INCLUDE"))
994                 {
995                 }
996               else
997                 {
998                   lex_error (lexer, NULL);
999                   goto error;
1000                 }
1001             }
1002         }
1003       else
1004         {
1005           lex_error (lexer, NULL);
1006           goto error;
1007         }
1008     }
1009
1010   if (frq.stats & BIT_INDEX (FRQ_ST_MEDIAN))
1011     {
1012         frq.percentiles =
1013           xrealloc (frq.percentiles, 
1014                     (frq.n_percentiles + 1)
1015                     * sizeof (*frq.percentiles));
1016         
1017         frq.percentiles[frq.n_percentiles].p = 0.50;
1018         frq.percentiles[frq.n_percentiles].show = true;
1019
1020         frq.n_percentiles++;
1021     }
1022
1023
1024 /* Figure out which charts the user requested.  */
1025
1026   {
1027     if (sbc_barchart)
1028       msg (SW, _("Bar charts are not implemented."));
1029
1030     if (sbc_histogram)
1031       {
1032         struct frq_chart *hist;
1033
1034         hist = frq.hist = xmalloc (sizeof *frq.hist);
1035         hist->x_min = hi_min;
1036         hist->x_max = hi_max;
1037         hist->y_scale = hi_scale;
1038         hist->y_max = hi_scale == FRQ_FREQ ? hi_freq : hi_pcnt;
1039         hist->draw_normal = hi_norm != FRQ_NONORMAL;
1040         hist->include_missing = false;
1041
1042         if (hist->x_min != SYSMIS && hist->x_max != SYSMIS
1043             && hist->x_min >= hist->x_max)
1044           {
1045             msg (SE, _("%s for histogram must be greater than or equal to %s, "
1046                        "but %s was specified as %.15g and %s as %.15g.  "
1047                        "%s and %s will be ignored."),
1048                  "MAX", "MIN", 
1049                  "MIN", hist->x_min, 
1050                  "MAX", hist->x_max,
1051                  "MIN", "MAX");
1052             hist->x_min = hist->x_max = SYSMIS;
1053           }
1054
1055         frq.percentiles =
1056           xrealloc (frq.percentiles, 
1057                     (frq.n_percentiles + 2)
1058                     * sizeof (*frq.percentiles));
1059         
1060         frq.percentiles[frq.n_percentiles].p = 0.25;
1061         frq.percentiles[frq.n_percentiles].show = false;
1062
1063         frq.percentiles[frq.n_percentiles + 1].p = 0.75;
1064         frq.percentiles[frq.n_percentiles + 1].show = false;
1065         
1066         frq.n_percentiles+=2;
1067       }
1068
1069     if (sbc_piechart)
1070       {
1071         struct frq_chart *pie;
1072
1073         pie = frq.pie = xmalloc (sizeof *frq.pie);
1074         pie->x_min = pie_min;
1075         pie->x_max = pie_max;
1076         pie->include_missing = pie_missing;
1077
1078         if (pie->x_min != SYSMIS && pie->x_max != SYSMIS
1079             && pie->x_min >= pie->x_max)
1080           {
1081             msg (SE, _("%s for pie chart must be greater than or equal to %s, "
1082                        "but %s was specified as %.15g and %s as %.15g.  "
1083                        "%s and %s will be ignored."), 
1084                  "MAX", "MIN", 
1085                  "MIN", pie->x_min,
1086                  "MAX", pie->x_max,
1087                  "MIN", "MAX");
1088             pie->x_min = pie->x_max = SYSMIS;
1089           }
1090       }
1091   }
1092
1093   {
1094     int i,o;
1095     double previous_p = -1;
1096     qsort (frq.percentiles, frq.n_percentiles,
1097            sizeof (*frq.percentiles), 
1098            ptile_3way);
1099
1100     frq.n_show_percentiles = 0;
1101     for (i = o = 0; i < frq.n_percentiles; ++i)
1102       {
1103         frq.percentiles[o].p = frq.percentiles[i].p;
1104
1105         if (frq.percentiles[i].show)
1106           frq.percentiles[o].show = true;
1107
1108         if (frq.percentiles[i].p != previous_p)
1109           {
1110             if (frq.percentiles[i].show)
1111               frq.n_show_percentiles++;
1112
1113             o++;
1114           }
1115
1116         previous_p = frq.percentiles[i].p;
1117       }
1118
1119     frq.n_percentiles = o;
1120   }
1121
1122   {
1123     struct casegrouper *grouper;
1124     struct casereader *group;
1125     bool ok;
1126
1127     grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
1128     while (casegrouper_get_next_group (grouper, &group))
1129       {
1130         struct ccase *c;
1131         precalc (&frq, group, ds);
1132         for (; (c = casereader_read (group)) != NULL; case_unref (c))
1133           calc (&frq, c, ds);
1134         postcalc (&frq, ds);
1135       }
1136     ok = casegrouper_destroy (grouper);
1137     ok = proc_commit (ds) && ok;
1138   }
1139
1140
1141   return CMD_SUCCESS;
1142
1143  error:
1144
1145   return CMD_FAILURE;
1146 }
1147
1148 static double
1149 calculate_iqr (const struct frq_proc *frq)
1150 {
1151   double q1 = SYSMIS;
1152   double q3 = SYSMIS;
1153   int i;
1154
1155   /* This cannot work unless the 25th and 75th percentile are calculated */
1156   assert (frq->n_percentiles >= 2);
1157   for (i = 0; i < frq->n_percentiles; i++)
1158     {
1159       struct percentile *pc = &frq->percentiles[i];
1160
1161       if (fabs (0.25 - pc->p) < DBL_EPSILON)
1162         q1 = pc->value;
1163       else if (fabs (0.75 - pc->p) < DBL_EPSILON)
1164         q3 = pc->value;
1165     }
1166
1167   return q1 == SYSMIS || q3 == SYSMIS ? SYSMIS : q3 - q1;
1168 }
1169
1170 static bool
1171 chart_includes_value (const struct frq_chart *chart,
1172                       const struct variable *var,
1173                       const union value *value)
1174 {
1175   if (!chart->include_missing && var_is_value_missing (var, value, MV_ANY))
1176     return false;
1177
1178   if (var_is_numeric (var)
1179       && ((chart->x_min != SYSMIS && value->f < chart->x_min)
1180           || (chart->x_max != SYSMIS && value->f > chart->x_max)))
1181     return false;
1182
1183   return true;
1184 }
1185
1186 /* Create a gsl_histogram from a freq_tab */
1187 static struct histogram *
1188 freq_tab_to_hist (const struct frq_proc *frq, const struct freq_tab *ft,
1189                   const struct variable *var)
1190 {
1191   double x_min, x_max, valid_freq;
1192   int i;
1193   double bin_width;
1194   struct histogram *histogram;
1195   double iqr;
1196
1197   /* Find out the extremes of the x value, within the range to be included in
1198      the histogram, and sum the total frequency of those values. */
1199   x_min = DBL_MAX;
1200   x_max = -DBL_MAX;
1201   valid_freq = 0;
1202   for (i = 0; i < ft->n_valid; i++)
1203     {
1204       const struct freq *f = &ft->valid[i];
1205       if (chart_includes_value (frq->hist, var, &f->value))
1206         {
1207           x_min = MIN (x_min, f->value.f);
1208           x_max = MAX (x_max, f->value.f);
1209           valid_freq += f->count;
1210         }
1211     }
1212
1213
1214   iqr = calculate_iqr (frq);
1215
1216   if (iqr > 0)
1217     /* Freedman-Diaconis' choice of bin width. */
1218     bin_width = 2 * iqr / pow (valid_freq, 1.0 / 3.0);
1219
1220   else
1221     /* Sturges Rule */
1222     bin_width = (x_max - x_min) / (1 + log2 (valid_freq));
1223
1224   histogram = histogram_create (bin_width, x_min, x_max);
1225
1226   if ( histogram == NULL)
1227     return NULL;
1228
1229   for (i = 0; i < ft->n_valid; i++)
1230     {
1231       const struct freq *f = &ft->valid[i];
1232       if (chart_includes_value (frq->hist, var, &f->value))
1233         histogram_add (histogram, f->value.f, f->count);
1234     }
1235
1236   return histogram;
1237 }
1238
1239 static int
1240 add_slice (const struct frq_chart *pie, const struct freq *freq,
1241            const struct variable *var, struct slice *slice)
1242 {
1243   if (chart_includes_value (pie, var, &freq->value))
1244     {
1245       ds_init_empty (&slice->label);
1246       var_append_value_name (var, &freq->value, &slice->label);
1247       slice->magnitude = freq->count;
1248       return 1;
1249     }
1250   else
1251     return 0;
1252 }
1253
1254 /* Allocate an array of slices and fill them from the data in frq_tab
1255    n_slices will contain the number of slices allocated.
1256    The caller is responsible for freeing slices
1257 */
1258 static struct slice *
1259 freq_tab_to_slice_array(const struct frq_chart *pie,
1260                         const struct freq_tab *frq_tab,
1261                         const struct variable *var,
1262                         int *n_slicesp)
1263 {
1264   struct slice *slices;
1265   int n_slices;
1266   int i;
1267
1268   slices = xnmalloc (frq_tab->n_valid + frq_tab->n_missing, sizeof *slices);
1269   n_slices = 0;
1270
1271   for (i = 0; i < frq_tab->n_valid; i++)
1272     n_slices += add_slice (pie, &frq_tab->valid[i], var, &slices[n_slices]);
1273   for (i = 0; i < frq_tab->n_missing; i++)
1274     n_slices += add_slice (pie, &frq_tab->missing[i], var, &slices[n_slices]);
1275
1276   *n_slicesp = n_slices;
1277   return slices;
1278 }
1279
1280
1281 static void
1282 do_piechart(const struct frq_chart *pie, const struct variable *var,
1283             const struct freq_tab *frq_tab)
1284 {
1285   struct slice *slices;
1286   int n_slices, i;
1287
1288   slices = freq_tab_to_slice_array (pie, frq_tab, var, &n_slices);
1289
1290   if (n_slices < 2)
1291     msg (SW, _("Omitting pie chart for %s, which has only %d unique values."),
1292          var_get_name (var), n_slices);
1293   else if (n_slices > 50)
1294     msg (SW, _("Omitting pie chart for %s, which has over 50 unique values."),
1295          var_get_name (var));
1296   else
1297     chart_item_submit (piechart_create (var_to_string(var), slices, n_slices));
1298
1299   for (i = 0; i < n_slices; i++)
1300     ds_destroy (&slices[i].label);
1301   free (slices);
1302 }
1303
1304 /* Calculates all the pertinent statistics for VF, putting them in array
1305    D[]. */
1306 static void
1307 calc_stats (const struct var_freqs *vf, double d[FRQ_ST_count])
1308 {
1309   const struct freq_tab *ft = &vf->tab;
1310   double W = ft->valid_cases;
1311   const struct freq *f;
1312   struct moments *m;
1313   int most_often = -1;
1314   double X_mode = SYSMIS;
1315
1316   /* Calculate the mode. */
1317   for (f = ft->valid; f < ft->missing; f++)
1318     {
1319       if (most_often < f->count)
1320         {
1321           most_often = f->count;
1322           X_mode = f->value.f;
1323         }
1324       else if (most_often == f->count)
1325         {
1326           /* A duplicate mode is undefined.
1327              FIXME: keep track of *all* the modes. */
1328           X_mode = SYSMIS;
1329         }
1330     }
1331
1332   /* Calculate moments. */
1333   m = moments_create (MOMENT_KURTOSIS);
1334   for (f = ft->valid; f < ft->missing; f++)
1335     moments_pass_one (m, f->value.f, f->count);
1336   for (f = ft->valid; f < ft->missing; f++)
1337     moments_pass_two (m, f->value.f, f->count);
1338   moments_calculate (m, NULL, &d[FRQ_ST_MEAN], &d[FRQ_ST_VARIANCE],
1339                      &d[FRQ_ST_SKEWNESS], &d[FRQ_ST_KURTOSIS]);
1340   moments_destroy (m);
1341
1342   /* Formulae below are taken from _SPSS Statistical Algorithms_. */
1343   d[FRQ_ST_MINIMUM] = ft->valid[0].value.f;
1344   d[FRQ_ST_MAXIMUM] = ft->valid[ft->n_valid - 1].value.f;
1345   d[FRQ_ST_MODE] = X_mode;
1346   d[FRQ_ST_RANGE] = d[FRQ_ST_MAXIMUM] - d[FRQ_ST_MINIMUM];
1347   d[FRQ_ST_SUM] = d[FRQ_ST_MEAN] * W;
1348   d[FRQ_ST_STDDEV] = sqrt (d[FRQ_ST_VARIANCE]);
1349   d[FRQ_ST_SEMEAN] = d[FRQ_ST_STDDEV] / sqrt (W);
1350   d[FRQ_ST_SESKEWNESS] = calc_seskew (W);
1351   d[FRQ_ST_SEKURTOSIS] = calc_sekurt (W);
1352 }
1353
1354 /* Displays a table of all the statistics requested for variable V. */
1355 static void
1356 dump_statistics (const struct frq_proc *frq, const struct var_freqs *vf,
1357                  const struct variable *wv)
1358 {
1359   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
1360   const struct freq_tab *ft = &vf->tab;
1361   double stat_value[FRQ_ST_count];
1362   struct tab_table *t;
1363   int i, r = 2; /* N missing and N valid are always dumped */
1364
1365   if (var_is_alpha (vf->var))
1366     return;
1367
1368   calc_stats (vf, stat_value);
1369
1370   t = tab_create (3, ((frq->stats & BIT_INDEX (FRQ_ST_MEDIAN)) ? frq->n_stats - 1 : frq->n_stats)
1371                                 + frq->n_show_percentiles + 2);
1372
1373   tab_set_format (t, RC_WEIGHT, wfmt);
1374   tab_box (t, TAL_1, TAL_1, -1, -1 , 0 , 0 , 2, tab_nr(t) - 1) ;
1375
1376   tab_vline (t, TAL_1 , 2, 0, tab_nr(t) - 1);
1377   tab_vline (t, TAL_GAP , 1, 0, tab_nr(t) - 1 ) ;
1378
1379   for (i = 0; i < FRQ_ST_count; i++)
1380     {
1381       if (FRQ_ST_MEDIAN == i)
1382         continue;
1383
1384       if (frq->stats & BIT_INDEX (i))
1385       {
1386         tab_text (t, 0, r, TAB_LEFT | TAT_TITLE,
1387                       gettext (st_name[i]));
1388
1389         if (vf->tab.n_valid <= 0 && r >= 2)
1390           tab_text (t, 2, r, 0,   ".");
1391         else
1392           tab_double (t, 2, r, TAB_NONE, stat_value[i], NULL, RC_OTHER);
1393         r++;
1394       }
1395     }
1396
1397   tab_text (t, 0, 0, TAB_LEFT | TAT_TITLE, _("N"));
1398   tab_text (t, 1, 0, TAB_LEFT | TAT_TITLE, _("Valid"));
1399   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Missing"));
1400
1401   tab_double (t, 2, 0, TAB_NONE, ft->valid_cases, NULL, RC_WEIGHT);
1402   tab_double (t, 2, 1, TAB_NONE, ft->total_cases - ft->valid_cases, NULL, RC_WEIGHT);
1403
1404   for (i = 0; i < frq->n_percentiles; i++)
1405     {
1406       const struct percentile *pc = &frq->percentiles[i];
1407
1408       if (!pc->show)
1409         continue;
1410
1411       if ( i == 0 )
1412         {
1413           tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Percentiles"));
1414         }
1415
1416       if (vf->tab.n_valid <= 0)
1417         {
1418           tab_text (t, 2, r, 0,   ".");
1419           ++r;
1420           continue;
1421         }
1422
1423       if (pc->p == 0.5)
1424         tab_text (t, 1, r, TAB_LEFT, _("50 (Median)"));
1425       else
1426         tab_double (t, 1, r, TAB_LEFT, pc->p * 100, NULL, RC_INTEGER);
1427       tab_double (t, 2, r, TAB_NONE, pc->value,
1428                   var_get_print_format (vf->var), RC_OTHER);
1429
1430       ++r;
1431     }
1432
1433   tab_title (t, "%s", var_to_string (vf->var));
1434
1435   tab_submit (t);
1436 }
1437