pintos-os.org Git - pspp/blob - src/language/stats/frequencies.c

   1 /*
   2   PSPP - a program for statistical analysis.
   3   Copyright (C) 1997-9, 2000, 2007, 2009, 2010, 2011, 2014, 2015 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
  49 #include "math/histogram.h"
  50 #include "math/moments.h"
  51 #include "math/chart-geometry.h"
  52
  53
  54 #include "output/charts/barchart.h"
  55 #include "output/charts/piechart.h"
  56 #include "output/charts/plot-hist.h"
  57 #include "output/pivot-table.h"
  58
  59 #include "gl/minmax.h"
  60 #include "gl/xalloc.h"
  61
  62 #include "gettext.h"
  63 #define _(msgid) gettext (msgid)
  64 #define N_(msgid) msgid
  65
  66 /* Percentiles to calculate. */
  67
  68 struct percentile
  69 {
  70   double p;        /* the %ile to be calculated */
  71   bool show;       /* True to show this percentile in the statistics box. */
  72 };
  73
  74 static int
  75 ptile_3way (const void *_p1, const void *_p2)
  76 {
  77   const struct percentile *p1 = _p1;
  78   const struct percentile *p2 = _p2;
  79
  80   if (p1->p < p2->p)
  81     return -1;
  82
  83   if (p1->p == p2->p)
  84     {
  85       if (p1->show > p2->show)
  86         return -1;
  87
  88       return (p1->show < p2->show);
  89     }
  90
  91   return (p1->p > p2->p);
  92 }
  93
  94
  95 enum
  96   {
  97     FRQ_NONORMAL,
  98     FRQ_NORMAL
  99   };
 100
 101 enum
 102   {
 103     FRQ_FREQ,
 104     FRQ_PERCENT
 105   };
 106
 107 enum sortprops
 108   {
 109     FRQ_AFREQ,
 110     FRQ_DFREQ,
 111     FRQ_AVALUE,
 112     FRQ_DVALUE
 113   };
 114
 115 /* Array indices for STATISTICS subcommand. */
 116 enum
 117   {
 118     FRQ_ST_MEAN,
 119     FRQ_ST_SEMEAN,
 120     FRQ_ST_MEDIAN,
 121     FRQ_ST_MODE,
 122     FRQ_ST_STDDEV,
 123     FRQ_ST_VARIANCE,
 124     FRQ_ST_KURTOSIS,
 125     FRQ_ST_SEKURTOSIS,
 126     FRQ_ST_SKEWNESS,
 127     FRQ_ST_SESKEWNESS,
 128     FRQ_ST_RANGE,
 129     FRQ_ST_MINIMUM,
 130     FRQ_ST_MAXIMUM,
 131     FRQ_ST_SUM,
 132     FRQ_ST_count
 133   };
 134
 135 /* Description of statistics. */
 136 static const char *st_name[FRQ_ST_count] =
 137 {
 138    N_("Mean"),
 139    N_("S.E. Mean"),
 140    N_("Median"),
 141    N_("Mode"),
 142    N_("Std Dev"),
 143    N_("Variance"),
 144    N_("Kurtosis"),
 145    N_("S.E. Kurt"),
 146    N_("Skewness"),
 147    N_("S.E. Skew"),
 148    N_("Range"),
 149    N_("Minimum"),
 150    N_("Maximum"),
 151    N_("Sum")
 152 };
 153
 154 struct freq_tab
 155   {
 156     struct hmap data;           /* Hash table for accumulating counts. */
 157     struct freq *valid;         /* Valid freqs. */
 158     int n_valid;                /* Number of total freqs. */
 159     const struct dictionary *dict; /* Source of entries in the table. */
 160
 161     struct freq *missing;       /* Missing freqs. */
 162     int n_missing;              /* Number of missing freqs. */
 163
 164     /* Statistics. */
 165     double total_cases;         /* Sum of weights of all cases. */
 166     double valid_cases;         /* Sum of weights of valid cases. */
 167   };
 168
 169 struct frq_chart
 170   {
 171     double x_min;               /* X axis minimum value. */
 172     double x_max;               /* X axis maximum value. */
 173     int y_scale;                /* Y axis scale: FRQ_FREQ or FRQ_PERCENT. */
 174
 175     /* Histograms only. */
 176     double y_max;               /* Y axis maximum value. */
 177     bool draw_normal;           /* Whether to draw normal curve. */
 178
 179     /* Pie charts only. */
 180     bool include_missing;       /* Whether to include missing values. */
 181   };
 182
 183 /* Per-variable frequency data. */
 184 struct var_freqs
 185   {
 186     const struct variable *var;
 187
 188     /* Freqency table. */
 189     struct freq_tab tab;        /* Frequencies table to use. */
 190
 191     /* Statistics. */
 192     double stat[FRQ_ST_count];
 193     double *percentiles;
 194
 195     /* Variable attributes. */
 196     int width;
 197   };
 198
 199 struct frq_proc
 200   {
 201     struct var_freqs *vars;
 202     size_t n_vars;
 203
 204     /* Percentiles to calculate and possibly display. */
 205     struct percentile *percentiles;
 206     size_t median_idx;
 207     size_t n_percentiles;
 208
 209     /* Frequency table display. */
 210     long int max_categories;         /* Maximum categories to show. */
 211     int sort;                   /* FRQ_AVALUE or FRQ_DVALUE
 212                                    or FRQ_AFREQ or FRQ_DFREQ. */
 213
 214     /* Statistics; number of statistics. */
 215     unsigned long stats;
 216     int n_stats;
 217
 218     /* Histogram and pie chart settings. */
 219     struct frq_chart *hist, *pie, *bar;
 220
 221     bool warn;
 222   };
 223
 224
 225 struct freq_compare_aux
 226   {
 227     bool by_freq;
 228     bool ascending_freq;
 229
 230     int width;
 231     bool ascending_value;
 232   };
 233
 234 static void calc_stats (const struct frq_proc *,
 235                         const struct var_freqs *, double d[FRQ_ST_count]);
 236
 237 static void do_piechart(const struct frq_chart *pie,
 238                         const struct variable *var,
 239                         const struct freq_tab *frq_tab);
 240
 241 static void do_barchart(const struct frq_chart *bar,
 242                         const struct variable **var,
 243                         const struct freq_tab *frq_tab);
 244
 245 static struct frq_stats_table *frq_stats_table_submit (
 246   struct frq_stats_table *, const struct frq_proc *,
 247   const struct dictionary *, const struct variable *wv);
 248 static void frq_stats_table_destroy (struct frq_stats_table *);
 249
 250 static int
 251 compare_freq (const void *a_, const void *b_, const void *aux_)
 252 {
 253   const struct freq_compare_aux *aux = aux_;
 254   const struct freq *a = a_;
 255   const struct freq *b = b_;
 256
 257   if (aux->by_freq && a->count != b->count)
 258     {
 259       int cmp = a->count > b->count ? 1 : -1;
 260       return aux->ascending_freq ? cmp : -cmp;
 261     }
 262   else
 263     {
 264       int cmp = value_compare_3way (a->values, b->values, aux->width);
 265       return aux->ascending_value ? cmp : -cmp;
 266     }
 267 }
 268
 269 /* Create a gsl_histogram from a freq_tab */
 270 static struct histogram *freq_tab_to_hist (const struct frq_proc *,
 271                                            const struct var_freqs *);
 272
 273 static void
 274 put_freq_row (struct pivot_table *table, int var_idx,
 275               double frequency, double percent,
 276               double valid_percent, double cum_percent)
 277 {
 278   double entries[] = { frequency, percent, valid_percent, cum_percent };
 279   for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
 280     if (entries[i] != SYSMIS)
 281       pivot_table_put2 (table, i, var_idx,
 282                         pivot_value_new_number (entries[i]));
 283 }
 284
 285 /* Displays a full frequency table for variable V. */
 286 static void
 287 dump_freq_table (const struct var_freqs *vf, const struct variable *wv)
 288 {
 289   const struct freq_tab *ft = &vf->tab;
 290
 291   struct pivot_table *table = pivot_table_create__ (pivot_value_new_variable (
 292                                                       vf->var), "Frequencies");
 293   pivot_table_set_weight_var (table, wv);
 294
 295   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 296                           N_("Frequency"), PIVOT_RC_COUNT,
 297                           N_("Percent"), PIVOT_RC_PERCENT,
 298                           N_("Valid Percent"), PIVOT_RC_PERCENT,
 299                           N_("Cumulative Percent"), PIVOT_RC_PERCENT);
 300
 301   struct pivot_dimension *variable = pivot_dimension_create__ (
 302     table, PIVOT_AXIS_ROW, pivot_value_new_variable (vf->var));
 303
 304   double cum_freq = 0.0;
 305   double cum_percent = 0.0;
 306   struct pivot_category *valid = NULL;
 307   for (const struct freq *f = ft->valid; f < ft->missing; f++)
 308     {
 309       cum_freq += f->count;
 310       double valid_percent = f->count / ft->valid_cases * 100.0;
 311       cum_percent += valid_percent;
 312
 313       if (!valid)
 314         valid = pivot_category_create_group (variable->root, N_("Valid"));
 315       int var_idx = pivot_category_create_leaf (
 316         valid, pivot_value_new_var_value (vf->var, &f->values[0]));
 317       put_freq_row (table, var_idx, f->count,
 318                     f->count / ft->total_cases * 100.0,
 319                     valid_percent, cum_percent);
 320     }
 321
 322   struct pivot_category *missing = NULL;
 323   size_t n_categories = ft->n_valid + ft->n_missing;
 324   for (const struct freq *f = ft->missing; f < &ft->valid[n_categories]; f++)
 325     {
 326       cum_freq += f->count;
 327
 328       if (!missing)
 329         missing = pivot_category_create_group (variable->root, N_("Missing"));
 330       int var_idx = pivot_category_create_leaf (
 331         missing, pivot_value_new_var_value (vf->var, &f->values[0]));
 332       put_freq_row (table, var_idx, f->count,
 333                     f->count / ft->total_cases * 100.0, SYSMIS, SYSMIS);
 334     }
 335
 336   int var_idx = pivot_category_create_leaf (
 337     variable->root, pivot_value_new_text (N_("Total")));
 338   put_freq_row (table, var_idx, cum_freq, cum_percent, SYSMIS, SYSMIS);
 339
 340   pivot_table_submit (table);
 341 }
 342 \f
 343 /* Statistical display. */
 344
 345 static double
 346 calc_percentile (double p, double valid_cases, double x1, double x2)
 347 {
 348   double s, dummy;
 349
 350   s = (settings_get_algorithm () != COMPATIBLE
 351        ? modf ((valid_cases - 1) * p, &dummy)
 352        : modf ((valid_cases + 1) * p - 1, &dummy));
 353
 354   return x1 + (x2 - x1) * s;
 355 }
 356
 357 /* Calculates all of the percentiles for VF within FRQ. */
 358 static void
 359 calc_percentiles (const struct frq_proc *frq, struct var_freqs *vf)
 360 {
 361   if (!frq->n_percentiles)
 362     return;
 363
 364   if (!vf->percentiles)
 365     vf->percentiles = xnmalloc (frq->n_percentiles, sizeof *vf->percentiles);
 366
 367   const struct freq_tab *ft = &vf->tab;
 368   const double W = ft->valid_cases;
 369   size_t idx = 0;
 370
 371   double rank = 0;
 372   for (const struct freq *f = ft->valid; f < ft->missing; f++)
 373     {
 374       rank += f->count;
 375       for (; idx < frq->n_percentiles; idx++)
 376         {
 377           struct percentile *pc = &frq->percentiles[idx];
 378           double tp;
 379
 380           tp = (settings_get_algorithm () == ENHANCED
 381                 ? (W - 1) * pc->p
 382                 : (W + 1) * pc->p - 1);
 383
 384           if (rank <= tp)
 385             break;
 386
 387           if (tp + 1 < rank || f + 1 >= ft->missing)
 388             vf->percentiles[idx] = f->values[0].f;
 389           else
 390             vf->percentiles[idx] = calc_percentile (pc->p, W, f->values[0].f,
 391                                                     f[1].values[0].f);
 392         }
 393     }
 394   for (; idx < frq->n_percentiles; idx++)
 395     vf->percentiles[idx] = (ft->n_valid > 0
 396                             ? ft->valid[ft->n_valid - 1].values[0].f
 397                             : SYSMIS);
 398 }
 399
 400 /* Returns true iff the value in struct freq F is non-missing
 401    for variable V. */
 402 static bool
 403 not_missing (const void *f_, const void *v_)
 404 {
 405   const struct freq *f = f_;
 406   const struct variable *v = v_;
 407
 408   return !var_is_value_missing (v, f->values);
 409 }
 410
 411
 412 /* Summarizes the frequency table data for variable V. */
 413 static void
 414 postprocess_freq_tab (const struct frq_proc *frq, struct var_freqs *vf)
 415 {
 416   struct freq_tab *ft = &vf->tab;
 417   struct freq_compare_aux aux;
 418   size_t count;
 419   struct freq *freqs, *f;
 420   size_t i;
 421
 422   /* Extract data from hash table. */
 423   count = hmap_count (&ft->data);
 424   freqs = freq_hmap_extract (&ft->data);
 425
 426   /* Put data into ft. */
 427   ft->valid = freqs;
 428   ft->n_valid = partition (freqs, count, sizeof *freqs, not_missing, vf->var);
 429   ft->missing = freqs + ft->n_valid;
 430   ft->n_missing = count - ft->n_valid;
 431
 432   /* Sort data. */
 433   aux.by_freq = frq->sort == FRQ_AFREQ || frq->sort == FRQ_DFREQ;
 434   aux.ascending_freq = frq->sort != FRQ_DFREQ;
 435   aux.width = vf->width;
 436   aux.ascending_value = frq->sort != FRQ_DVALUE;
 437   sort (ft->valid, ft->n_valid, sizeof *ft->valid, compare_freq, &aux);
 438   sort (ft->missing, ft->n_missing, sizeof *ft->missing, compare_freq, &aux);
 439
 440   /* Summary statistics. */
 441   ft->valid_cases = 0.0;
 442   for(i = 0 ;  i < ft->n_valid ; ++i)
 443     {
 444       f = &ft->valid[i];
 445       ft->valid_cases += f->count;
 446
 447     }
 448
 449   ft->total_cases = ft->valid_cases ;
 450   for(i = 0 ;  i < ft->n_missing ; ++i)
 451     {
 452       f = &ft->missing[i];
 453       ft->total_cases += f->count;
 454     }
 455
 456 }
 457
 458 /* Add data from case C to the frequency table. */
 459 static void
 460 calc (struct frq_proc *frq, const struct ccase *c, const struct dataset *ds)
 461 {
 462   double weight = dict_get_case_weight (dataset_dict (ds), c, &frq->warn);
 463   size_t i;
 464
 465   for (i = 0; i < frq->n_vars; i++)
 466     {
 467       struct var_freqs *vf = &frq->vars[i];
 468       const union value *value = case_data (c, vf->var);
 469       size_t hash = value_hash (value, vf->width, 0);
 470       struct freq *f;
 471
 472       f = freq_hmap_search (&vf->tab.data, value, vf->width, hash);
 473       if (f == NULL)
 474         f = freq_hmap_insert (&vf->tab.data, value, vf->width, hash);
 475
 476       f->count += weight;
 477     }
 478 }
 479
 480 static void
 481 output_splits_once (bool *need_splits, const struct dataset *ds,
 482                     const struct ccase *c)
 483 {
 484   if (*need_splits)
 485     {
 486       output_split_file_values (ds, c);
 487       *need_splits = false;
 488     }
 489 }
 490
 491 /* Finishes up with the variables after frequencies have been
 492    calculated.  Displays statistics, percentiles, ... */
 493 static struct frq_stats_table *
 494 postcalc (struct frq_proc *frq, const struct dataset *ds,
 495           struct ccase *first, struct frq_stats_table *fst)
 496 {
 497   const struct dictionary *dict = dataset_dict (ds);
 498   const struct variable *wv = dict_get_weight (dict);
 499
 500   for (size_t i = 0; i < frq->n_vars; i++)
 501     {
 502       struct var_freqs *vf = &frq->vars[i];
 503       postprocess_freq_tab (frq, vf);
 504       calc_percentiles (frq, vf);
 505     }
 506
 507   enum split_type st = dict_get_split_type (dict);
 508   bool need_splits = true;
 509   if (frq->n_stats)
 510     {
 511       if (st != SPLIT_LAYERED)
 512         output_splits_once (&need_splits, ds, first);
 513       fst = frq_stats_table_submit (fst, frq, dict, wv);
 514     }
 515
 516   for (size_t i = 0; i < frq->n_vars; i++)
 517     {
 518       struct var_freqs *vf = &frq->vars[i];
 519
 520       /* Frequencies tables. */
 521       if (vf->tab.n_valid + vf->tab.n_missing <= frq->max_categories)
 522         {
 523           output_splits_once (&need_splits, ds, first);
 524           dump_freq_table (vf, wv);
 525         }
 526
 527       if (frq->hist && var_is_numeric (vf->var) && vf->tab.n_valid > 0)
 528         {
 529           double d[FRQ_ST_count];
 530           struct histogram *histogram;
 531
 532           calc_stats (frq, vf, d);
 533
 534           histogram = freq_tab_to_hist (frq, vf);
 535
 536           if (histogram)
 537             {
 538               output_splits_once (&need_splits, ds, first);
 539               chart_submit (histogram_chart_create (
 540                               histogram->gsl_hist, var_to_string(vf->var),
 541                               vf->tab.valid_cases,
 542                               d[FRQ_ST_MEAN],
 543                               d[FRQ_ST_STDDEV],
 544                               frq->hist->draw_normal));
 545
 546               statistic_destroy (&histogram->parent);
 547             }
 548         }
 549
 550       if (frq->pie)
 551         {
 552           output_splits_once (&need_splits, ds, first);
 553           do_piechart(frq->pie, vf->var, &vf->tab);
 554         }
 555
 556       if (frq->bar)
 557         {
 558           output_splits_once (&need_splits, ds, first);
 559           do_barchart(frq->bar, &vf->var, &vf->tab);
 560         }
 561
 562       free (vf->tab.valid);
 563       freq_hmap_destroy (&vf->tab.data, vf->width);
 564     }
 565
 566   return fst;
 567 }
 568
 569 static void
 570 frq_run (struct frq_proc *frq, struct dataset *ds)
 571 {
 572   struct frq_stats_table *fst = NULL;
 573   struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds),
 574                                                            dataset_dict (ds));
 575   struct casereader *group;
 576   while (casegrouper_get_next_group (grouper, &group))
 577     {
 578       for (size_t i = 0; i < frq->n_vars; i++)
 579         hmap_init (&frq->vars[i].tab.data);
 580
 581       struct ccase *first = casereader_peek (group, 0);
 582
 583       struct ccase *c;
 584       for (; (c = casereader_read (group)) != NULL; case_unref (c))
 585         calc (frq, c, ds);
 586       fst = postcalc (frq, ds, first, fst);
 587       casereader_destroy (group);
 588
 589       case_unref (first);
 590     }
 591   frq_stats_table_destroy (fst);
 592   casegrouper_destroy (grouper);
 593   proc_commit (ds);
 594 }
 595
 596 int
 597 cmd_frequencies (struct lexer *lexer, struct dataset *ds)
 598 {
 599   int i;
 600   struct frq_proc frq;
 601   const struct variable **vars = NULL;
 602
 603   bool sbc_barchart = false;
 604   bool sbc_piechart = false;
 605   bool sbc_histogram = false;
 606
 607   double pie_min = -DBL_MAX;
 608   double pie_max = DBL_MAX;
 609   bool pie_missing = true;
 610
 611   double bar_min = -DBL_MAX;
 612   double bar_max = DBL_MAX;
 613   bool bar_freq = true;
 614
 615   double hi_min = -DBL_MAX;
 616   double hi_max = DBL_MAX;
 617   int hi_scale = FRQ_FREQ;
 618   int hi_freq = INT_MIN;
 619   int hi_pcnt = INT_MIN;
 620   int hi_norm = FRQ_NONORMAL;
 621
 622   frq.sort = FRQ_AVALUE;
 623
 624   frq.vars = NULL;
 625   frq.n_vars = 0;
 626
 627   frq.stats = BIT_INDEX (FRQ_ST_MEAN)
 628     | BIT_INDEX (FRQ_ST_STDDEV)
 629     | BIT_INDEX (FRQ_ST_MINIMUM)
 630     | BIT_INDEX (FRQ_ST_MAXIMUM);
 631
 632   frq.n_stats = 4;
 633
 634   frq.max_categories = LONG_MAX;
 635
 636   frq.percentiles = NULL;
 637   frq.n_percentiles = 0;
 638
 639   frq.hist = NULL;
 640   frq.pie = NULL;
 641   frq.bar = NULL;
 642   frq.warn = true;
 643
 644
 645   /* Accept an optional, completely pointless "/VARIABLES=" */
 646   lex_match (lexer, T_SLASH);
 647   if (lex_match_id  (lexer, "VARIABLES"))
 648     {
 649       if (! lex_force_match (lexer, T_EQUALS))
 650         goto error;
 651     }
 652
 653   if (!parse_variables_const (lexer, dataset_dict (ds),
 654                               &vars,
 655                               &frq.n_vars,
 656                               PV_NO_DUPLICATE))
 657     goto error;
 658
 659   frq.vars = xcalloc (frq.n_vars, sizeof (*frq.vars));
 660   for (i = 0; i < frq.n_vars; ++i)
 661     {
 662       frq.vars[i].var = vars[i];
 663       frq.vars[i].width = var_get_width (vars[i]);
 664     }
 665
 666   while (lex_token (lexer) != T_ENDCMD)
 667     {
 668       lex_match (lexer, T_SLASH);
 669
 670       if (lex_match_id (lexer, "STATISTICS"))
 671         {
 672           frq.stats = BIT_INDEX (FRQ_ST_MEAN)
 673             | BIT_INDEX (FRQ_ST_STDDEV)
 674             | BIT_INDEX (FRQ_ST_MINIMUM)
 675             | BIT_INDEX (FRQ_ST_MAXIMUM);
 676
 677           frq.n_stats = 4;
 678
 679           if (lex_match (lexer, T_EQUALS))
 680             {
 681               frq.n_stats = 0;
 682               frq.stats = 0;
 683             }
 684
 685           while (lex_token (lexer) != T_ENDCMD
 686                  && lex_token (lexer) != T_SLASH)
 687             {
 688               if (lex_match_id (lexer, "DEFAULT"))
 689                 {
 690                   frq.stats = BIT_INDEX (FRQ_ST_MEAN)
 691                     | BIT_INDEX (FRQ_ST_STDDEV)
 692                     | BIT_INDEX (FRQ_ST_MINIMUM)
 693                     | BIT_INDEX (FRQ_ST_MAXIMUM);
 694
 695                   frq.n_stats = 4;
 696                 }
 697               else if (lex_match_id (lexer, "MEAN"))
 698                 {
 699                   frq.stats |= BIT_INDEX (FRQ_ST_MEAN);
 700                   frq.n_stats++;
 701                 }
 702               else if (lex_match_id (lexer, "SEMEAN"))
 703                 {
 704                   frq.stats |= BIT_INDEX (FRQ_ST_SEMEAN);
 705                   frq.n_stats++;
 706                 }
 707               else if (lex_match_id (lexer, "MEDIAN"))
 708                 {
 709                   frq.stats |= BIT_INDEX (FRQ_ST_MEDIAN);
 710                   frq.n_stats++;
 711                 }
 712               else if (lex_match_id (lexer, "MODE"))
 713                 {
 714                   frq.stats |= BIT_INDEX (FRQ_ST_MODE);
 715                   frq.n_stats++;
 716                 }
 717               else if (lex_match_id (lexer, "STDDEV"))
 718                 {
 719                   frq.stats |= BIT_INDEX (FRQ_ST_STDDEV);
 720                   frq.n_stats++;
 721                 }
 722               else if (lex_match_id (lexer, "VARIANCE"))
 723                 {
 724                   frq.stats |= BIT_INDEX (FRQ_ST_VARIANCE);
 725                   frq.n_stats++;
 726                 }
 727               else if (lex_match_id (lexer, "KURTOSIS"))
 728                 {
 729                   frq.stats |= BIT_INDEX (FRQ_ST_KURTOSIS);
 730                   frq.n_stats++;
 731                 }
 732               else if (lex_match_id (lexer, "SKEWNESS"))
 733                 {
 734                   frq.stats |= BIT_INDEX (FRQ_ST_SKEWNESS);
 735                   frq.n_stats++;
 736                 }
 737               else if (lex_match_id (lexer, "RANGE"))
 738                 {
 739                   frq.stats |= BIT_INDEX (FRQ_ST_RANGE);
 740                   frq.n_stats++;
 741                 }
 742               else if (lex_match_id (lexer, "MINIMUM"))
 743                 {
 744                   frq.stats |= BIT_INDEX (FRQ_ST_MINIMUM);
 745                   frq.n_stats++;
 746                 }
 747               else if (lex_match_id (lexer, "MAXIMUM"))
 748                 {
 749                   frq.stats |= BIT_INDEX (FRQ_ST_MAXIMUM);
 750                   frq.n_stats++;
 751                 }
 752               else if (lex_match_id (lexer, "SUM"))
 753                 {
 754                   frq.stats |= BIT_INDEX (FRQ_ST_SUM);
 755                   frq.n_stats++;
 756                 }
 757               else if (lex_match_id (lexer, "SESKEWNESS"))
 758                 {
 759                   frq.stats |= BIT_INDEX (FRQ_ST_SESKEWNESS);
 760                   frq.n_stats++;
 761                 }
 762               else if (lex_match_id (lexer, "SEKURTOSIS"))
 763                 {
 764                   frq.stats |= BIT_INDEX (FRQ_ST_SEKURTOSIS);
 765                   frq.n_stats++;
 766                 }
 767               else if (lex_match_id (lexer, "NONE"))
 768                 {
 769                   frq.stats = 0;
 770                   frq.n_stats = 0;
 771                 }
 772               else if (lex_match (lexer, T_ALL))
 773                 {
 774                   frq.stats = ~0;
 775                   frq.n_stats = FRQ_ST_count;
 776                 }
 777               else
 778                 {
 779                   lex_error (lexer, NULL);
 780                   goto error;
 781                 }
 782             }
 783         }
 784       else if (lex_match_id (lexer, "PERCENTILES"))
 785         {
 786           lex_match (lexer, T_EQUALS);
 787           while (lex_token (lexer) != T_ENDCMD
 788                  && lex_token (lexer) != T_SLASH)
 789             {
 790               if (lex_force_num (lexer))
 791                 {
 792                   frq.percentiles =
 793                     xrealloc (frq.percentiles,
 794                               (frq.n_percentiles + 1)
 795                               * sizeof (*frq.percentiles));
 796                   frq.percentiles[frq.n_percentiles].p = lex_number (lexer)  / 100.0;
 797                   frq.percentiles[frq.n_percentiles].show = true;
 798                   lex_get (lexer);
 799                   frq.n_percentiles++;
 800                 }
 801               else
 802                 {
 803                   lex_error (lexer, NULL);
 804                   goto error;
 805                 }
 806               lex_match (lexer, T_COMMA);
 807             }
 808         }
 809       else if (lex_match_id (lexer, "FORMAT"))
 810         {
 811           lex_match (lexer, T_EQUALS);
 812           while (lex_token (lexer) != T_ENDCMD
 813                  && lex_token (lexer) != T_SLASH)
 814             {
 815               if (lex_match_id (lexer, "TABLE"))
 816                 {
 817                 }
 818               else if (lex_match_id (lexer, "NOTABLE"))
 819                 {
 820                   frq.max_categories = 0;
 821                 }
 822               else if (lex_match_id (lexer, "LIMIT"))
 823                 {
 824                   if (!lex_force_match (lexer, T_LPAREN)
 825                       || !lex_force_int_range (lexer, "LIMIT", 0, INT_MAX))
 826                     goto error;
 827
 828                   frq.max_categories = lex_integer (lexer);
 829                   lex_get (lexer);
 830
 831                   if (!lex_force_match (lexer, T_RPAREN))
 832                     goto error;
 833                 }
 834               else if (lex_match_id (lexer, "AVALUE"))
 835                 {
 836                   frq.sort = FRQ_AVALUE;
 837                 }
 838               else if (lex_match_id (lexer, "DVALUE"))
 839                 {
 840                   frq.sort = FRQ_DVALUE;
 841                 }
 842               else if (lex_match_id (lexer, "AFREQ"))
 843                 {
 844                   frq.sort = FRQ_AFREQ;
 845                 }
 846               else if (lex_match_id (lexer, "DFREQ"))
 847                 {
 848                   frq.sort = FRQ_DFREQ;
 849                 }
 850               else
 851                 {
 852                   lex_error (lexer, NULL);
 853                   goto error;
 854                 }
 855             }
 856         }
 857       else if (lex_match_id (lexer, "NTILES"))
 858         {
 859           lex_match (lexer, T_EQUALS);
 860
 861           if (lex_force_int_range (lexer, "NTILES", 0, INT_MAX))
 862             {
 863               int n = lex_integer (lexer);
 864               lex_get (lexer);
 865               for (int i = 0; i < n + 1; ++i)
 866                 {
 867                   frq.percentiles =
 868                     xrealloc (frq.percentiles,
 869                               (frq.n_percentiles + 1)
 870                               * sizeof (*frq.percentiles));
 871                   frq.percentiles[frq.n_percentiles].p =
 872                     i / (double) n ;
 873                   frq.percentiles[frq.n_percentiles].show = true;
 874
 875                   frq.n_percentiles++;
 876                 }
 877             }
 878           else
 879             {
 880               lex_error (lexer, NULL);
 881               goto error;
 882             }
 883         }
 884       else if (lex_match_id (lexer, "ALGORITHM"))
 885         {
 886           lex_match (lexer, T_EQUALS);
 887
 888           if (lex_match_id (lexer, "COMPATIBLE"))
 889             {
 890               settings_set_cmd_algorithm (COMPATIBLE);
 891             }
 892           else if (lex_match_id (lexer, "ENHANCED"))
 893             {
 894               settings_set_cmd_algorithm (ENHANCED);
 895             }
 896           else
 897             {
 898               lex_error (lexer, NULL);
 899               goto error;
 900             }
 901         }
 902       else if (lex_match_id (lexer, "HISTOGRAM"))
 903         {
 904           lex_match (lexer, T_EQUALS);
 905           sbc_histogram = true;
 906
 907           while (lex_token (lexer) != T_ENDCMD
 908                  && lex_token (lexer) != T_SLASH)
 909             {
 910               if (lex_match_id (lexer, "NORMAL"))
 911                 {
 912                   hi_norm = FRQ_NORMAL;
 913                 }
 914               else if (lex_match_id (lexer, "NONORMAL"))
 915                 {
 916                   hi_norm = FRQ_NONORMAL;
 917                 }
 918               else if (lex_match_id (lexer, "FREQ"))
 919                 {
 920                   hi_scale = FRQ_FREQ;
 921                   if (lex_match (lexer, T_LPAREN))
 922                     {
 923                       if (lex_force_int_range (lexer, "FREQ", 1, INT_MAX))
 924                         {
 925                           hi_freq = lex_integer (lexer);
 926                           lex_get (lexer);
 927                           if (! lex_force_match (lexer, T_RPAREN))
 928                             goto error;
 929                         }
 930                     }
 931                 }
 932               else if (lex_match_id (lexer, "PERCENT"))
 933                 {
 934                   hi_scale = FRQ_PERCENT;
 935                   if (lex_match (lexer, T_LPAREN))
 936                     {
 937                       if (lex_force_int_range (lexer, "PERCENT", 1, INT_MAX))
 938                         {
 939                           hi_pcnt = lex_integer (lexer);
 940                           lex_get (lexer);
 941                           if (! lex_force_match (lexer, T_RPAREN))
 942                             goto error;
 943                         }
 944                     }
 945                 }
 946               else if (lex_match_id (lexer, "MINIMUM"))
 947                 {
 948                   if (! lex_force_match (lexer, T_LPAREN))
 949                     goto error;
 950                   if (lex_force_num (lexer))
 951                     {
 952                       hi_min = lex_number (lexer);
 953                       lex_get (lexer);
 954                     }
 955                   if (! lex_force_match (lexer, T_RPAREN))
 956                     goto error;
 957                 }
 958               else if (lex_match_id (lexer, "MAXIMUM"))
 959                 {
 960                   if (! lex_force_match (lexer, T_LPAREN))
 961                     goto error;
 962                   if (lex_force_num (lexer))
 963                     {
 964                       hi_max = lex_number (lexer);
 965                       lex_get (lexer);
 966                     }
 967                   if (! lex_force_match (lexer, T_RPAREN))
 968                     goto error;
 969                 }
 970               else
 971                 {
 972                   lex_error (lexer, NULL);
 973                   goto error;
 974                 }
 975             }
 976         }
 977       else if (lex_match_id (lexer, "PIECHART"))
 978         {
 979           lex_match (lexer, T_EQUALS);
 980           while (lex_token (lexer) != T_ENDCMD
 981                  && lex_token (lexer) != T_SLASH)
 982             {
 983               if (lex_match_id (lexer, "MINIMUM"))
 984                 {
 985                   if (! lex_force_match (lexer, T_LPAREN))
 986                     goto error;
 987                   if (lex_force_num (lexer))
 988                     {
 989                       pie_min = lex_number (lexer);
 990                       lex_get (lexer);
 991                     }
 992                   if (! lex_force_match (lexer, T_RPAREN))
 993                     goto error;
 994                 }
 995               else if (lex_match_id (lexer, "MAXIMUM"))
 996                 {
 997                   if (! lex_force_match (lexer, T_LPAREN))
 998                     goto error;
 999                   if (lex_force_num (lexer))
1000                     {
1001                       pie_max = lex_number (lexer);
1002                       lex_get (lexer);
1003                     }
1004                   if (! lex_force_match (lexer, T_RPAREN))
1005                     goto error;
1006                 }
1007               else if (lex_match_id (lexer, "MISSING"))
1008                 {
1009                   pie_missing = true;
1010                 }
1011               else if (lex_match_id (lexer, "NOMISSING"))
1012                 {
1013                   pie_missing = false;
1014                 }
1015               else
1016                 {
1017                   lex_error (lexer, NULL);
1018                   goto error;
1019                 }
1020             }
1021           sbc_piechart = true;
1022         }
1023       else if (lex_match_id (lexer, "BARCHART"))
1024         {
1025           lex_match (lexer, T_EQUALS);
1026           while (lex_token (lexer) != T_ENDCMD
1027                  && lex_token (lexer) != T_SLASH)
1028             {
1029               if (lex_match_id (lexer, "MINIMUM"))
1030                 {
1031                   if (! lex_force_match (lexer, T_LPAREN))
1032                     goto error;
1033                   if (lex_force_num (lexer))
1034                     {
1035                       bar_min = lex_number (lexer);
1036                       lex_get (lexer);
1037                     }
1038                   if (! lex_force_match (lexer, T_RPAREN))
1039                     goto error;
1040                 }
1041               else if (lex_match_id (lexer, "MAXIMUM"))
1042                 {
1043                   if (! lex_force_match (lexer, T_LPAREN))
1044                     goto error;
1045                   if (lex_force_num (lexer))
1046                     {
1047                       bar_max = lex_number (lexer);
1048                       lex_get (lexer);
1049                     }
1050                   if (! lex_force_match (lexer, T_RPAREN))
1051                     goto error;
1052                 }
1053               else if (lex_match_id (lexer, "FREQ"))
1054                 {
1055                   if (lex_match (lexer, T_LPAREN))
1056                     {
1057                       if (lex_force_num (lexer))
1058                         {
1059                           lex_number (lexer);
1060                           lex_get (lexer);
1061                         }
1062                       if (! lex_force_match (lexer, T_RPAREN))
1063                         goto error;
1064                     }
1065                   bar_freq = true;
1066                 }
1067               else if (lex_match_id (lexer, "PERCENT"))
1068                 {
1069                   if (lex_match (lexer, T_LPAREN))
1070                     {
1071                       if (lex_force_num (lexer))
1072                         {
1073                           lex_number (lexer);
1074                           lex_get (lexer);
1075                         }
1076                       if (! lex_force_match (lexer, T_RPAREN))
1077                         goto error;
1078                     }
1079                   bar_freq = false;
1080                 }
1081               else
1082                 {
1083                   lex_error (lexer, NULL);
1084                   goto error;
1085                 }
1086             }
1087           sbc_barchart = true;
1088         }
1089       else if (lex_match_id (lexer, "MISSING"))
1090         {
1091           lex_match (lexer, T_EQUALS);
1092
1093           while (lex_token (lexer) != T_ENDCMD
1094                  && lex_token (lexer) != T_SLASH)
1095             {
1096               if (lex_match_id (lexer, "EXCLUDE"))
1097                 {
1098                 }
1099               else if (lex_match_id (lexer, "INCLUDE"))
1100                 {
1101                 }
1102               else
1103                 {
1104                   lex_error (lexer, NULL);
1105                   goto error;
1106                 }
1107             }
1108         }
1109       else if (lex_match_id (lexer, "ORDER"))
1110         {
1111           lex_match (lexer, T_EQUALS);
1112           if (!lex_match_id (lexer, "ANALYSIS"))
1113             lex_match_id (lexer, "VARIABLE");
1114         }
1115       else
1116         {
1117           lex_error (lexer, NULL);
1118           goto error;
1119         }
1120     }
1121
1122   if (frq.stats & BIT_INDEX (FRQ_ST_MEDIAN))
1123     {
1124         frq.percentiles =
1125           xrealloc (frq.percentiles,
1126                     (frq.n_percentiles + 1)
1127                     * sizeof (*frq.percentiles));
1128
1129         frq.percentiles[frq.n_percentiles].p = 0.50;
1130         frq.percentiles[frq.n_percentiles].show = false;
1131
1132         frq.n_percentiles++;
1133     }
1134
1135
1136 /* Figure out which charts the user requested.  */
1137
1138   {
1139     if (sbc_histogram)
1140       {
1141         struct frq_chart *hist;
1142
1143         hist = frq.hist = xmalloc (sizeof *frq.hist);
1144         hist->x_min = hi_min;
1145         hist->x_max = hi_max;
1146         hist->y_scale = hi_scale;
1147         hist->y_max = hi_scale == FRQ_FREQ ? hi_freq : hi_pcnt;
1148         hist->draw_normal = hi_norm != FRQ_NONORMAL;
1149         hist->include_missing = false;
1150
1151         if (hist->x_min != SYSMIS && hist->x_max != SYSMIS
1152             && hist->x_min >= hist->x_max)
1153           {
1154             msg (SE, _("%s for histogram must be greater than or equal to %s, "
1155                        "but %s was specified as %.15g and %s as %.15g.  "
1156                        "%s and %s will be ignored."),
1157                  "MAX", "MIN",
1158                  "MIN", hist->x_min,
1159                  "MAX", hist->x_max,
1160                  "MIN", "MAX");
1161             hist->x_min = hist->x_max = SYSMIS;
1162           }
1163
1164         frq.percentiles =
1165           xrealloc (frq.percentiles,
1166                     (frq.n_percentiles + 2)
1167                     * sizeof (*frq.percentiles));
1168
1169         frq.percentiles[frq.n_percentiles].p = 0.25;
1170         frq.percentiles[frq.n_percentiles].show = false;
1171
1172         frq.percentiles[frq.n_percentiles + 1].p = 0.75;
1173         frq.percentiles[frq.n_percentiles + 1].show = false;
1174
1175         frq.n_percentiles+=2;
1176       }
1177
1178     if (sbc_barchart)
1179       {
1180         frq.bar = xmalloc (sizeof *frq.bar);
1181         frq.bar->x_min = bar_min;
1182         frq.bar->x_max = bar_max;
1183         frq.bar->include_missing = false;
1184         frq.bar->y_scale = bar_freq ? FRQ_FREQ : FRQ_PERCENT;
1185       }
1186
1187     if (sbc_piechart)
1188       {
1189         struct frq_chart *pie;
1190
1191         pie = frq.pie = xmalloc (sizeof *frq.pie);
1192         pie->x_min = pie_min;
1193         pie->x_max = pie_max;
1194         pie->include_missing = pie_missing;
1195
1196         if (pie->x_min != SYSMIS && pie->x_max != SYSMIS
1197             && pie->x_min >= pie->x_max)
1198           {
1199             msg (SE, _("%s for pie chart must be greater than or equal to %s, "
1200                        "but %s was specified as %.15g and %s as %.15g.  "
1201                        "%s and %s will be ignored."),
1202                  "MAX", "MIN",
1203                  "MIN", pie->x_min,
1204                  "MAX", pie->x_max,
1205                  "MIN", "MAX");
1206             pie->x_min = pie->x_max = SYSMIS;
1207           }
1208       }
1209   }
1210
1211   {
1212     int i,o;
1213     double previous_p = -1;
1214     qsort (frq.percentiles, frq.n_percentiles,
1215            sizeof (*frq.percentiles),
1216            ptile_3way);
1217
1218     for (i = o = 0; i < frq.n_percentiles; ++i)
1219       {
1220         if (frq.percentiles[i].p != previous_p)
1221           {
1222             frq.percentiles[o].p = frq.percentiles[i].p;
1223             frq.percentiles[o].show = frq.percentiles[i].show;
1224             o++;
1225           }
1226         else if (frq.percentiles[i].show &&
1227                  !frq.percentiles[o].show)
1228           {
1229             frq.percentiles[o].show = true;
1230           }
1231         previous_p = frq.percentiles[i].p;
1232       }
1233
1234     frq.n_percentiles = o;
1235
1236     frq.median_idx = SIZE_MAX;
1237     for (i = 0; i < frq.n_percentiles; i++)
1238       if (frq.percentiles[i].p == 0.5)
1239         {
1240           frq.median_idx = i;
1241           break;
1242         }
1243   }
1244
1245   frq_run (&frq, ds);
1246
1247   free (vars);
1248   for (size_t i = 0; i < frq.n_vars; i++)
1249     free (frq.vars[i].percentiles);
1250   free (frq.vars);
1251   free (frq.bar);
1252   free (frq.pie);
1253   free (frq.hist);
1254   free (frq.percentiles);
1255
1256   return CMD_SUCCESS;
1257
1258  error:
1259
1260   free (vars);
1261   for (size_t i = 0; i < frq.n_vars; i++)
1262     free (frq.vars[i].percentiles);
1263   free (frq.vars);
1264   free (frq.bar);
1265   free (frq.pie);
1266   free (frq.hist);
1267   free (frq.percentiles);
1268
1269   return CMD_FAILURE;
1270 }
1271
1272 static double
1273 calculate_iqr (const struct frq_proc *frq, const struct var_freqs *vf)
1274 {
1275   double q1 = SYSMIS;
1276   double q3 = SYSMIS;
1277   int i;
1278
1279   /* This cannot work unless the 25th and 75th percentile are calculated */
1280   assert (frq->n_percentiles >= 2);
1281   for (i = 0; i < frq->n_percentiles; i++)
1282     {
1283       struct percentile *pc = &frq->percentiles[i];
1284
1285       if (fabs (0.25 - pc->p) < DBL_EPSILON)
1286         q1 = vf->percentiles[i];
1287       else if (fabs (0.75 - pc->p) < DBL_EPSILON)
1288         q3 = vf->percentiles[i];
1289     }
1290
1291   return q1 == SYSMIS || q3 == SYSMIS ? SYSMIS : q3 - q1;
1292 }
1293
1294 static bool
1295 chart_includes_value (const struct frq_chart *chart,
1296                       const struct variable *var,
1297                       const union value *value)
1298 {
1299   if (!chart->include_missing && var_is_value_missing (var, value))
1300     return false;
1301
1302   if (var_is_numeric (var)
1303       && ((chart->x_min != SYSMIS && value->f < chart->x_min)
1304           || (chart->x_max != SYSMIS && value->f > chart->x_max)))
1305     return false;
1306
1307   return true;
1308 }
1309
1310 /* Create a gsl_histogram from a freq_tab */
1311 static struct histogram *
1312 freq_tab_to_hist (const struct frq_proc *frq, const struct var_freqs *vf)
1313 {
1314   /* Find out the extremes of the x value, within the range to be included in
1315      the histogram, and sum the total frequency of those values. */
1316   double x_min = DBL_MAX;
1317   double x_max = -DBL_MAX;
1318   double valid_freq = 0;
1319   for (int i = 0; i < vf->tab.n_valid; i++)
1320     {
1321       const struct freq *f = &vf->tab.valid[i];
1322       if (chart_includes_value (frq->hist, vf->var, f->values))
1323         {
1324           x_min = MIN (x_min, f->values[0].f);
1325           x_max = MAX (x_max, f->values[0].f);
1326           valid_freq += f->count;
1327         }
1328     }
1329
1330   if (valid_freq <= 0)
1331     return NULL;
1332
1333   double iqr = calculate_iqr (frq, vf);
1334
1335   double bin_width =
1336     (iqr > 0
1337      ? 2 * iqr / pow (valid_freq, 1.0 / 3.0)       /* Freedman-Diaconis. */
1338      : (x_max - x_min) / (1 + log2 (valid_freq))); /* Sturges */
1339
1340   struct histogram *histogram = histogram_create (bin_width, x_min, x_max);
1341   if (histogram == NULL)
1342     return NULL;
1343
1344   for (int i = 0; i < vf->tab.n_valid; i++)
1345     {
1346       const struct freq *f = &vf->tab.valid[i];
1347       if (chart_includes_value (frq->hist, vf->var, f->values))
1348         histogram_add (histogram, f->values[0].f, f->count);
1349     }
1350
1351   return histogram;
1352 }
1353
1354
1355 /* Allocate an array of struct freqs and fill them from the data in FRQ_TAB,
1356    according to the parameters of CATCHART
1357    N_SLICES will contain the number of slices allocated.
1358    The caller is responsible for freeing slices
1359 */
1360 static struct freq *
1361 pick_cat_counts (const struct frq_chart *catchart,
1362                  const struct freq_tab *frq_tab,
1363                  int *n_slicesp)
1364 {
1365   int n_slices = 0;
1366   int i;
1367   struct freq *slices = xnmalloc (frq_tab->n_valid + frq_tab->n_missing, sizeof *slices);
1368
1369   for (i = 0; i < frq_tab->n_valid; i++)
1370     {
1371       const struct freq *f = &frq_tab->valid[i];
1372       if (f->count > catchart->x_max)
1373         continue;
1374
1375       if (f->count < catchart->x_min)
1376         continue;
1377
1378       slices[n_slices] = *f;
1379
1380       n_slices++;
1381     }
1382
1383   if (catchart->include_missing)
1384     {
1385       for (i = 0; i < frq_tab->n_missing; i++)
1386         {
1387           const struct freq *f = &frq_tab->missing[i];
1388           slices[n_slices].count += f->count;
1389
1390           if (i == 0)
1391             slices[n_slices].values[0] = f->values[0];
1392         }
1393
1394       if (frq_tab->n_missing > 0)
1395         n_slices++;
1396     }
1397
1398   *n_slicesp = n_slices;
1399   return slices;
1400 }
1401
1402
1403 /* Allocate an array of struct freqs and fill them from the data in FRQ_TAB,
1404    according to the parameters of CATCHART
1405    N_SLICES will contain the number of slices allocated.
1406    The caller is responsible for freeing slices
1407 */
1408 static struct freq **
1409 pick_cat_counts_ptr (const struct frq_chart *catchart,
1410                      const struct freq_tab *frq_tab,
1411                      int *n_slicesp)
1412 {
1413   int n_slices = 0;
1414   int i;
1415   struct freq **slices = xnmalloc (frq_tab->n_valid + frq_tab->n_missing, sizeof *slices);
1416
1417   for (i = 0; i < frq_tab->n_valid; i++)
1418     {
1419       struct freq *f = &frq_tab->valid[i];
1420       if (f->count > catchart->x_max)
1421         continue;
1422
1423       if (f->count < catchart->x_min)
1424         continue;
1425
1426       slices[n_slices] = f;
1427
1428       n_slices++;
1429     }
1430
1431   if (catchart->include_missing)
1432     {
1433       for (i = 0; i < frq_tab->n_missing; i++)
1434         {
1435           const struct freq *f = &frq_tab->missing[i];
1436           if (i == 0)
1437             {
1438               slices[n_slices] = xmalloc (sizeof (struct freq));
1439               slices[n_slices]->values[0] = f->values[0];
1440             }
1441
1442           slices[n_slices]->count += f->count;
1443
1444         }
1445     }
1446
1447   *n_slicesp = n_slices;
1448   return slices;
1449 }
1450
1451
1452
1453 static void
1454 do_piechart(const struct frq_chart *pie, const struct variable *var,
1455             const struct freq_tab *frq_tab)
1456 {
1457   int n_slices;
1458   struct freq *slices = pick_cat_counts (pie, frq_tab, &n_slices);
1459
1460   if (n_slices < 2)
1461     msg (SW, _("Omitting pie chart for %s, which has only %d unique values."),
1462          var_get_name (var), n_slices);
1463   else if (n_slices > 50)
1464     msg (SW, _("Omitting pie chart for %s, which has over 50 unique values."),
1465          var_get_name (var));
1466   else
1467     chart_submit (piechart_create (var, slices, n_slices));
1468
1469   free (slices);
1470 }
1471
1472
1473 static void
1474 do_barchart(const struct frq_chart *bar, const struct variable **var,
1475             const struct freq_tab *frq_tab)
1476 {
1477   int n_slices;
1478   struct freq **slices = pick_cat_counts_ptr (bar, frq_tab, &n_slices);
1479
1480   if (n_slices < 1)
1481     msg (SW, _("Omitting bar chart, which has no values."));
1482   else
1483     chart_submit (barchart_create (var, 1,
1484                                    (bar->y_scale == FRQ_FREQ) ? _("Count") : _("Percent"),
1485                                    (bar->y_scale == FRQ_PERCENT),
1486                                    slices, n_slices));
1487   free (slices);
1488 }
1489
1490
1491 /* Calculates all the pertinent statistics for VF, putting them in array
1492    D[]. */
1493 static void
1494 calc_stats (const struct frq_proc *frq, const struct var_freqs *vf,
1495             double d[FRQ_ST_count])
1496 {
1497   const struct freq_tab *ft = &vf->tab;
1498   double W = ft->valid_cases;
1499   const struct freq *f;
1500   struct moments *m;
1501   int most_often = -1;
1502   double X_mode = SYSMIS;
1503
1504   /* Calculate the mode. */
1505   for (f = ft->valid; f < ft->missing; f++)
1506     {
1507       if (most_often < f->count)
1508         {
1509           most_often = f->count;
1510           X_mode = f->values[0].f;
1511         }
1512       else if (most_often == f->count)
1513         {
1514           /* A duplicate mode is undefined.
1515              FIXME: keep track of *all* the modes. */
1516           X_mode = SYSMIS;
1517         }
1518     }
1519
1520   /* Calculate moments. */
1521   m = moments_create (MOMENT_KURTOSIS);
1522   for (f = ft->valid; f < ft->missing; f++)
1523     moments_pass_one (m, f->values[0].f, f->count);
1524   for (f = ft->valid; f < ft->missing; f++)
1525     moments_pass_two (m, f->values[0].f, f->count);
1526   moments_calculate (m, NULL, &d[FRQ_ST_MEAN], &d[FRQ_ST_VARIANCE],
1527                      &d[FRQ_ST_SKEWNESS], &d[FRQ_ST_KURTOSIS]);
1528   moments_destroy (m);
1529
1530   /* Formulae below are taken from _SPSS Statistical Algorithms_. */
1531   if (ft->n_valid > 0)
1532     {
1533       d[FRQ_ST_MINIMUM] = ft->valid[0].values[0].f;
1534       d[FRQ_ST_MAXIMUM] = ft->valid[ft->n_valid - 1].values[0].f;
1535       d[FRQ_ST_RANGE] = d[FRQ_ST_MAXIMUM] - d[FRQ_ST_MINIMUM];
1536     }
1537   else
1538     {
1539       d[FRQ_ST_MINIMUM] = SYSMIS;
1540       d[FRQ_ST_MAXIMUM] = SYSMIS;
1541       d[FRQ_ST_RANGE] = SYSMIS;
1542     }
1543   d[FRQ_ST_MODE] = X_mode;
1544   d[FRQ_ST_SUM] = d[FRQ_ST_MEAN] * W;
1545   d[FRQ_ST_STDDEV] = sqrt (d[FRQ_ST_VARIANCE]);
1546   d[FRQ_ST_SEMEAN] = d[FRQ_ST_STDDEV] / sqrt (W);
1547   d[FRQ_ST_SESKEWNESS] = calc_seskew (W);
1548   d[FRQ_ST_SEKURTOSIS] = calc_sekurt (W);
1549   d[FRQ_ST_MEDIAN] = (frq->median_idx != SIZE_MAX
1550                       ? vf->percentiles[frq->median_idx]
1551                       : SYSMIS);
1552 }
1553
1554 static bool
1555 all_string_variables (const struct frq_proc *frq)
1556 {
1557   for (size_t i = 0; i < frq->n_vars; i++)
1558     if (var_is_numeric (frq->vars[i].var))
1559       return false;
1560
1561   return true;
1562 }
1563 \f
1564 struct frq_stats_table
1565   {
1566     struct pivot_table *table;
1567     struct pivot_splits *splits;
1568   };
1569
1570 /* Displays a table of all the statistics requested. */
1571 static struct frq_stats_table *
1572 frq_stats_table_create (const struct frq_proc *frq,
1573                         const struct dictionary *dict,
1574                         const struct variable *wv)
1575 {
1576   if (all_string_variables (frq))
1577     return NULL;
1578
1579   struct pivot_table *table = pivot_table_create (N_("Statistics"));
1580   pivot_table_set_weight_var (table, wv);
1581
1582   struct pivot_splits *splits = pivot_splits_create (table, dict);
1583
1584   struct pivot_dimension *variables
1585     = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Variables"));
1586   for (size_t i = 0; i < frq->n_vars; i++)
1587     if (!var_is_alpha (frq->vars[i].var))
1588       pivot_category_create_leaf (variables->root,
1589                                   pivot_value_new_variable (frq->vars[i].var));
1590
1591   struct pivot_dimension *statistics = pivot_dimension_create (
1592     table, PIVOT_AXIS_ROW, N_("Statistics"));
1593   struct pivot_category *n = pivot_category_create_group (
1594     statistics->root, N_("N"));
1595   pivot_category_create_leaves (n,
1596                                 N_("Valid"), PIVOT_RC_COUNT,
1597                                 N_("Missing"), PIVOT_RC_COUNT);
1598   for (int i = 0; i < FRQ_ST_count; i++)
1599     if (frq->stats & BIT_INDEX (i))
1600       pivot_category_create_leaf (statistics->root,
1601                                   pivot_value_new_text (st_name[i]));
1602   struct pivot_category *percentiles = NULL;
1603   for (size_t i = 0; i < frq->n_percentiles; i++)
1604     {
1605       const struct percentile *pc = &frq->percentiles[i];
1606
1607       if (!pc->show)
1608         continue;
1609
1610       if (!percentiles)
1611         percentiles = pivot_category_create_group (
1612           statistics->root, N_("Percentiles"));
1613       pivot_category_create_leaf (percentiles, pivot_value_new_integer (
1614                                     pc->p * 100.0));
1615     }
1616
1617   struct frq_stats_table *fst = xmalloc (sizeof *fst);
1618   *fst = (struct frq_stats_table) { .table = table, .splits = splits };
1619   return fst;
1620 }
1621
1622 static struct frq_stats_table *
1623 frq_stats_table_submit (struct frq_stats_table *fst,
1624                         const struct frq_proc *frq,
1625                         const struct dictionary *dict,
1626                         const struct variable *wv)
1627 {
1628   if (!fst)
1629     {
1630       fst = frq_stats_table_create (frq, dict, wv);
1631       if (!fst)
1632         return NULL;
1633     }
1634
1635   int var_idx = 0;
1636   for (size_t i = 0; i < frq->n_vars; i++)
1637     {
1638       struct var_freqs *vf = &frq->vars[i];
1639       if (var_is_alpha (vf->var))
1640         continue;
1641
1642       const struct freq_tab *ft = &vf->tab;
1643
1644       int row = 0;
1645       pivot_table_put2 (fst->table, var_idx, row++,
1646                         pivot_value_new_number (ft->valid_cases));
1647       pivot_table_put2 (fst->table, var_idx, row++,
1648                         pivot_value_new_number (
1649                           ft->total_cases - ft->valid_cases));
1650
1651       double stat_values[FRQ_ST_count];
1652       calc_stats (frq, vf, stat_values);
1653       for (int j = 0; j < FRQ_ST_count; j++)
1654         {
1655           if (!(frq->stats & BIT_INDEX (j)))
1656             continue;
1657
1658           union value v = { .f = vf->tab.n_valid ? stat_values[j] : SYSMIS };
1659           struct pivot_value *pv
1660             = (j == FRQ_ST_MODE || j == FRQ_ST_MINIMUM || j == FRQ_ST_MAXIMUM
1661                ? pivot_value_new_var_value (vf->var, &v)
1662                : pivot_value_new_number (v.f));
1663           pivot_table_put2 (fst->table, var_idx, row++, pv);
1664         }
1665
1666       for (size_t j = 0; j < frq->n_percentiles; j++)
1667         {
1668           const struct percentile *pc = &frq->percentiles[j];
1669           if (!pc->show)
1670             continue;
1671
1672           union value v = {
1673             .f = vf->tab.n_valid ? vf->percentiles[j] : SYSMIS
1674           };
1675           pivot_table_put2 (fst->table, var_idx, row++,
1676                             pivot_value_new_var_value (vf->var, &v));
1677         }
1678
1679       var_idx++;
1680     }
1681
1682   if (!fst->splits)
1683     {
1684       frq_stats_table_destroy (fst);
1685       return NULL;
1686     }
1687   return fst;
1688 }
1689
1690 static void
1691 frq_stats_table_destroy (struct frq_stats_table *fst)
1692 {
1693   if (!fst)
1694     return;
1695
1696   pivot_table_submit (fst->table);
1697   pivot_splits_destroy (fst->splits);
1698   free (fst);
1699 }