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