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 void dump_statistics (const struct frq_proc *frq,
 246                              const struct variable *wv);
 247
 248 static int
 249 compare_freq (const void *a_, const void *b_, const void *aux_)
 250 {
 251   const struct freq_compare_aux *aux = aux_;
 252   const struct freq *a = a_;
 253   const struct freq *b = b_;
 254
 255   if (aux->by_freq && a->count != b->count)
 256     {
 257       int cmp = a->count > b->count ? 1 : -1;
 258       return aux->ascending_freq ? cmp : -cmp;
 259     }
 260   else
 261     {
 262       int cmp = value_compare_3way (a->values, b->values, aux->width);
 263       return aux->ascending_value ? cmp : -cmp;
 264     }
 265 }
 266
 267 /* Create a gsl_histogram from a freq_tab */
 268 static struct histogram *freq_tab_to_hist (const struct frq_proc *,
 269                                            const struct var_freqs *);
 270
 271 static void
 272 put_freq_row (struct pivot_table *table, int var_idx,
 273               double frequency, double percent,
 274               double valid_percent, double cum_percent)
 275 {
 276   double entries[] = { frequency, percent, valid_percent, cum_percent };
 277   for (size_t i = 0; i < sizeof entries / sizeof *entries; i++)
 278     if (entries[i] != SYSMIS)
 279       pivot_table_put2 (table, i, var_idx,
 280                         pivot_value_new_number (entries[i]));
 281 }
 282
 283 /* Displays a full frequency table for variable V. */
 284 static void
 285 dump_freq_table (const struct var_freqs *vf, const struct variable *wv)
 286 {
 287   const struct freq_tab *ft = &vf->tab;
 288
 289   struct pivot_table *table = pivot_table_create__ (pivot_value_new_variable (
 290                                                       vf->var), "Frequencies");
 291   pivot_table_set_weight_var (table, wv);
 292
 293   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 294                           N_("Frequency"), PIVOT_RC_COUNT,
 295                           N_("Percent"), PIVOT_RC_PERCENT,
 296                           N_("Valid Percent"), PIVOT_RC_PERCENT,
 297                           N_("Cumulative Percent"), PIVOT_RC_PERCENT);
 298
 299   struct pivot_dimension *variable = pivot_dimension_create__ (
 300     table, PIVOT_AXIS_ROW, pivot_value_new_variable (vf->var));
 301
 302   double cum_freq = 0.0;
 303   double cum_percent = 0.0;
 304   struct pivot_category *valid = NULL;
 305   for (const struct freq *f = ft->valid; f < ft->missing; f++)
 306     {
 307       cum_freq += f->count;
 308       double valid_percent = f->count / ft->valid_cases * 100.0;
 309       cum_percent += valid_percent;
 310
 311       if (!valid)
 312         valid = pivot_category_create_group (variable->root, N_("Valid"));
 313       int var_idx = pivot_category_create_leaf (
 314         valid, pivot_value_new_var_value (vf->var, &f->values[0]));
 315       put_freq_row (table, var_idx, f->count,
 316                     f->count / ft->total_cases * 100.0,
 317                     valid_percent, cum_percent);
 318     }
 319
 320   struct pivot_category *missing = NULL;
 321   size_t n_categories = ft->n_valid + ft->n_missing;
 322   for (const struct freq *f = ft->missing; f < &ft->valid[n_categories]; f++)
 323     {
 324       cum_freq += f->count;
 325
 326       if (!missing)
 327         missing = pivot_category_create_group (variable->root, N_("Missing"));
 328       int var_idx = pivot_category_create_leaf (
 329         missing, pivot_value_new_var_value (vf->var, &f->values[0]));
 330       put_freq_row (table, var_idx, f->count,
 331                     f->count / ft->total_cases * 100.0, SYSMIS, SYSMIS);
 332     }
 333
 334   int var_idx = pivot_category_create_leaf (
 335     variable->root, pivot_value_new_text (N_("Total")));
 336   put_freq_row (table, var_idx, cum_freq, cum_percent, SYSMIS, SYSMIS);
 337
 338   pivot_table_submit (table);
 339 }
 340 \f
 341 /* Statistical display. */
 342
 343 static double
 344 calc_percentile (double p, double valid_cases, double x1, double x2)
 345 {
 346   double s, dummy;
 347
 348   s = (settings_get_algorithm () != COMPATIBLE
 349        ? modf ((valid_cases - 1) * p, &dummy)
 350        : modf ((valid_cases + 1) * p - 1, &dummy));
 351
 352   return x1 + (x2 - x1) * s;
 353 }
 354
 355 /* Calculates all of the percentiles for VF within FRQ. */
 356 static void
 357 calc_percentiles (const struct frq_proc *frq, struct var_freqs *vf)
 358 {
 359   if (!frq->n_percentiles)
 360     return;
 361
 362   if (!vf->percentiles)
 363     vf->percentiles = xnmalloc (frq->n_percentiles, sizeof *vf->percentiles);
 364
 365   const struct freq_tab *ft = &vf->tab;
 366   const double W = ft->valid_cases;
 367   size_t idx = 0;
 368
 369   double rank = 0;
 370   for (const struct freq *f = ft->valid; f < ft->missing; f++)
 371     {
 372       rank += f->count;
 373       for (; idx < frq->n_percentiles; idx++)
 374         {
 375           struct percentile *pc = &frq->percentiles[idx];
 376           double tp;
 377
 378           tp = (settings_get_algorithm () == ENHANCED
 379                 ? (W - 1) * pc->p
 380                 : (W + 1) * pc->p - 1);
 381
 382           if (rank <= tp)
 383             break;
 384
 385           if (tp + 1 < rank || f + 1 >= ft->missing)
 386             vf->percentiles[idx] = f->values[0].f;
 387           else
 388             vf->percentiles[idx] = calc_percentile (pc->p, W, f->values[0].f,
 389                                                     f[1].values[0].f);
 390         }
 391     }
 392   for (; idx < frq->n_percentiles; idx++)
 393     vf->percentiles[idx] = (ft->n_valid > 0
 394                             ? ft->valid[ft->n_valid - 1].values[0].f
 395                             : SYSMIS);
 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);
 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_idx = SIZE_MAX;
1206     for (i = 0; i < frq.n_percentiles; i++)
1207       if (frq.percentiles[i].p == 0.5)
1208         {
1209           frq.median_idx = 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   free (vars);
1235   for (size_t i = 0; i < frq.n_vars; i++)
1236     free (frq.vars[i].percentiles);
1237   free (frq.vars);
1238   free (frq.bar);
1239   free (frq.pie);
1240   free (frq.hist);
1241   free (frq.percentiles);
1242
1243   return CMD_SUCCESS;
1244
1245  error:
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_FAILURE;
1257 }
1258
1259 static double
1260 calculate_iqr (const struct frq_proc *frq, const struct var_freqs *vf)
1261 {
1262   double q1 = SYSMIS;
1263   double q3 = SYSMIS;
1264   int i;
1265
1266   /* This cannot work unless the 25th and 75th percentile are calculated */
1267   assert (frq->n_percentiles >= 2);
1268   for (i = 0; i < frq->n_percentiles; i++)
1269     {
1270       struct percentile *pc = &frq->percentiles[i];
1271
1272       if (fabs (0.25 - pc->p) < DBL_EPSILON)
1273         q1 = vf->percentiles[i];
1274       else if (fabs (0.75 - pc->p) < DBL_EPSILON)
1275         q3 = vf->percentiles[i];
1276     }
1277
1278   return q1 == SYSMIS || q3 == SYSMIS ? SYSMIS : q3 - q1;
1279 }
1280
1281 static bool
1282 chart_includes_value (const struct frq_chart *chart,
1283                       const struct variable *var,
1284                       const union value *value)
1285 {
1286   if (!chart->include_missing && var_is_value_missing (var, value))
1287     return false;
1288
1289   if (var_is_numeric (var)
1290       && ((chart->x_min != SYSMIS && value->f < chart->x_min)
1291           || (chart->x_max != SYSMIS && value->f > chart->x_max)))
1292     return false;
1293
1294   return true;
1295 }
1296
1297 /* Create a gsl_histogram from a freq_tab */
1298 static struct histogram *
1299 freq_tab_to_hist (const struct frq_proc *frq, const struct var_freqs *vf)
1300 {
1301   /* Find out the extremes of the x value, within the range to be included in
1302      the histogram, and sum the total frequency of those values. */
1303   double x_min = DBL_MAX;
1304   double x_max = -DBL_MAX;
1305   double valid_freq = 0;
1306   for (int i = 0; i < vf->tab.n_valid; i++)
1307     {
1308       const struct freq *f = &vf->tab.valid[i];
1309       if (chart_includes_value (frq->hist, vf->var, f->values))
1310         {
1311           x_min = MIN (x_min, f->values[0].f);
1312           x_max = MAX (x_max, f->values[0].f);
1313           valid_freq += f->count;
1314         }
1315     }
1316
1317   if (valid_freq <= 0)
1318     return NULL;
1319
1320   double iqr = calculate_iqr (frq, vf);
1321
1322   double bin_width =
1323     (iqr > 0
1324      ? 2 * iqr / pow (valid_freq, 1.0 / 3.0)       /* Freedman-Diaconis. */
1325      : (x_max - x_min) / (1 + log2 (valid_freq))); /* Sturges */
1326
1327   struct histogram *histogram = histogram_create (bin_width, x_min, x_max);
1328   if (histogram == NULL)
1329     return NULL;
1330
1331   for (int i = 0; i < vf->tab.n_valid; i++)
1332     {
1333       const struct freq *f = &vf->tab.valid[i];
1334       if (chart_includes_value (frq->hist, vf->var, f->values))
1335         histogram_add (histogram, f->values[0].f, f->count);
1336     }
1337
1338   return histogram;
1339 }
1340
1341
1342 /* Allocate an array of struct freqs and fill them from the data in FRQ_TAB,
1343    according to the parameters of CATCHART
1344    N_SLICES will contain the number of slices allocated.
1345    The caller is responsible for freeing slices
1346 */
1347 static struct freq *
1348 pick_cat_counts (const struct frq_chart *catchart,
1349                  const struct freq_tab *frq_tab,
1350                  int *n_slicesp)
1351 {
1352   int n_slices = 0;
1353   int i;
1354   struct freq *slices = xnmalloc (frq_tab->n_valid + frq_tab->n_missing, sizeof *slices);
1355
1356   for (i = 0; i < frq_tab->n_valid; i++)
1357     {
1358       const struct freq *f = &frq_tab->valid[i];
1359       if (f->count > catchart->x_max)
1360         continue;
1361
1362       if (f->count < catchart->x_min)
1363         continue;
1364
1365       slices[n_slices] = *f;
1366
1367       n_slices++;
1368     }
1369
1370   if (catchart->include_missing)
1371     {
1372       for (i = 0; i < frq_tab->n_missing; i++)
1373         {
1374           const struct freq *f = &frq_tab->missing[i];
1375           slices[n_slices].count += f->count;
1376
1377           if (i == 0)
1378             slices[n_slices].values[0] = f->values[0];
1379         }
1380
1381       if (frq_tab->n_missing > 0)
1382         n_slices++;
1383     }
1384
1385   *n_slicesp = n_slices;
1386   return slices;
1387 }
1388
1389
1390 /* Allocate an array of struct freqs and fill them from the data in FRQ_TAB,
1391    according to the parameters of CATCHART
1392    N_SLICES will contain the number of slices allocated.
1393    The caller is responsible for freeing slices
1394 */
1395 static struct freq **
1396 pick_cat_counts_ptr (const struct frq_chart *catchart,
1397                      const struct freq_tab *frq_tab,
1398                      int *n_slicesp)
1399 {
1400   int n_slices = 0;
1401   int i;
1402   struct freq **slices = xnmalloc (frq_tab->n_valid + frq_tab->n_missing, sizeof *slices);
1403
1404   for (i = 0; i < frq_tab->n_valid; i++)
1405     {
1406       struct freq *f = &frq_tab->valid[i];
1407       if (f->count > catchart->x_max)
1408         continue;
1409
1410       if (f->count < catchart->x_min)
1411         continue;
1412
1413       slices[n_slices] = f;
1414
1415       n_slices++;
1416     }
1417
1418   if (catchart->include_missing)
1419     {
1420       for (i = 0; i < frq_tab->n_missing; i++)
1421         {
1422           const struct freq *f = &frq_tab->missing[i];
1423           if (i == 0)
1424             {
1425               slices[n_slices] = xmalloc (sizeof (struct freq));
1426               slices[n_slices]->values[0] = f->values[0];
1427             }
1428
1429           slices[n_slices]->count += f->count;
1430
1431         }
1432     }
1433
1434   *n_slicesp = n_slices;
1435   return slices;
1436 }
1437
1438
1439
1440 static void
1441 do_piechart(const struct frq_chart *pie, const struct variable *var,
1442             const struct freq_tab *frq_tab)
1443 {
1444   int n_slices;
1445   struct freq *slices = pick_cat_counts (pie, frq_tab, &n_slices);
1446
1447   if (n_slices < 2)
1448     msg (SW, _("Omitting pie chart for %s, which has only %d unique values."),
1449          var_get_name (var), n_slices);
1450   else if (n_slices > 50)
1451     msg (SW, _("Omitting pie chart for %s, which has over 50 unique values."),
1452          var_get_name (var));
1453   else
1454     chart_submit (piechart_create (var, slices, n_slices));
1455
1456   free (slices);
1457 }
1458
1459
1460 static void
1461 do_barchart(const struct frq_chart *bar, const struct variable **var,
1462             const struct freq_tab *frq_tab)
1463 {
1464   int n_slices;
1465   struct freq **slices = pick_cat_counts_ptr (bar, frq_tab, &n_slices);
1466
1467   if (n_slices < 1)
1468     msg (SW, _("Omitting bar chart, which has no values."));
1469   else
1470     chart_submit (barchart_create (var, 1,
1471                                    (bar->y_scale == FRQ_FREQ) ? _("Count") : _("Percent"),
1472                                    (bar->y_scale == FRQ_PERCENT),
1473                                    slices, n_slices));
1474   free (slices);
1475 }
1476
1477
1478 /* Calculates all the pertinent statistics for VF, putting them in array
1479    D[]. */
1480 static void
1481 calc_stats (const struct frq_proc *frq, const struct var_freqs *vf,
1482             double d[FRQ_ST_count])
1483 {
1484   const struct freq_tab *ft = &vf->tab;
1485   double W = ft->valid_cases;
1486   const struct freq *f;
1487   struct moments *m;
1488   int most_often = -1;
1489   double X_mode = SYSMIS;
1490
1491   /* Calculate the mode. */
1492   for (f = ft->valid; f < ft->missing; f++)
1493     {
1494       if (most_often < f->count)
1495         {
1496           most_often = f->count;
1497           X_mode = f->values[0].f;
1498         }
1499       else if (most_often == f->count)
1500         {
1501           /* A duplicate mode is undefined.
1502              FIXME: keep track of *all* the modes. */
1503           X_mode = SYSMIS;
1504         }
1505     }
1506
1507   /* Calculate moments. */
1508   m = moments_create (MOMENT_KURTOSIS);
1509   for (f = ft->valid; f < ft->missing; f++)
1510     moments_pass_one (m, f->values[0].f, f->count);
1511   for (f = ft->valid; f < ft->missing; f++)
1512     moments_pass_two (m, f->values[0].f, f->count);
1513   moments_calculate (m, NULL, &d[FRQ_ST_MEAN], &d[FRQ_ST_VARIANCE],
1514                      &d[FRQ_ST_SKEWNESS], &d[FRQ_ST_KURTOSIS]);
1515   moments_destroy (m);
1516
1517   /* Formulae below are taken from _SPSS Statistical Algorithms_. */
1518   if (ft->n_valid > 0)
1519     {
1520       d[FRQ_ST_MINIMUM] = ft->valid[0].values[0].f;
1521       d[FRQ_ST_MAXIMUM] = ft->valid[ft->n_valid - 1].values[0].f;
1522       d[FRQ_ST_RANGE] = d[FRQ_ST_MAXIMUM] - d[FRQ_ST_MINIMUM];
1523     }
1524   else
1525     {
1526       d[FRQ_ST_MINIMUM] = SYSMIS;
1527       d[FRQ_ST_MAXIMUM] = SYSMIS;
1528       d[FRQ_ST_RANGE] = SYSMIS;
1529     }
1530   d[FRQ_ST_MODE] = X_mode;
1531   d[FRQ_ST_SUM] = d[FRQ_ST_MEAN] * W;
1532   d[FRQ_ST_STDDEV] = sqrt (d[FRQ_ST_VARIANCE]);
1533   d[FRQ_ST_SEMEAN] = d[FRQ_ST_STDDEV] / sqrt (W);
1534   d[FRQ_ST_SESKEWNESS] = calc_seskew (W);
1535   d[FRQ_ST_SEKURTOSIS] = calc_sekurt (W);
1536   d[FRQ_ST_MEDIAN] = (frq->median_idx != SIZE_MAX
1537                       ? vf->percentiles[frq->median_idx]
1538                       : SYSMIS);
1539 }
1540
1541 static bool
1542 all_string_variables (const struct frq_proc *frq)
1543 {
1544   for (size_t i = 0; i < frq->n_vars; i++)
1545     if (var_is_numeric (frq->vars[i].var))
1546       return false;
1547
1548   return true;
1549 }
1550
1551 /* Displays a table of all the statistics requested. */
1552 static void
1553 dump_statistics (const struct frq_proc *frq, const struct variable *wv)
1554 {
1555   if (all_string_variables (frq))
1556     return;
1557
1558   struct pivot_table *table = pivot_table_create (N_("Statistics"));
1559   pivot_table_set_weight_var (table, wv);
1560
1561   struct pivot_dimension *variables
1562     = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Variables"));
1563
1564   struct pivot_dimension *statistics = pivot_dimension_create (
1565     table, PIVOT_AXIS_ROW, N_("Statistics"));
1566   struct pivot_category *n = pivot_category_create_group (
1567     statistics->root, N_("N"));
1568   pivot_category_create_leaves (n,
1569                                 N_("Valid"), PIVOT_RC_COUNT,
1570                                 N_("Missing"), PIVOT_RC_COUNT);
1571   for (int i = 0; i < FRQ_ST_count; i++)
1572     if (frq->stats & BIT_INDEX (i))
1573       pivot_category_create_leaf (statistics->root,
1574                                   pivot_value_new_text (st_name[i]));
1575   struct pivot_category *percentiles = NULL;
1576   for (size_t i = 0; i < frq->n_percentiles; i++)
1577     {
1578       const struct percentile *pc = &frq->percentiles[i];
1579
1580       if (!pc->show)
1581         continue;
1582
1583       if (!percentiles)
1584         percentiles = pivot_category_create_group (
1585           statistics->root, N_("Percentiles"));
1586       pivot_category_create_leaf (percentiles, pivot_value_new_integer (
1587                                     pc->p * 100.0));
1588     }
1589
1590   for (size_t i = 0; i < frq->n_vars; i++)
1591     {
1592       struct var_freqs *vf = &frq->vars[i];
1593       if (var_is_alpha (vf->var))
1594         continue;
1595
1596       const struct freq_tab *ft = &vf->tab;
1597
1598       int var_idx = pivot_category_create_leaf (
1599         variables->root, pivot_value_new_variable (vf->var));
1600
1601       int row = 0;
1602       pivot_table_put2 (table, var_idx, row++,
1603                         pivot_value_new_number (ft->valid_cases));
1604       pivot_table_put2 (table, var_idx, row++,
1605                         pivot_value_new_number (
1606                           ft->total_cases - ft->valid_cases));
1607
1608       double stat_values[FRQ_ST_count];
1609       calc_stats (frq, vf, stat_values);
1610       for (int j = 0; j < FRQ_ST_count; j++)
1611         {
1612           if (!(frq->stats & BIT_INDEX (j)))
1613             continue;
1614
1615           union value v = { .f = vf->tab.n_valid ? stat_values[j] : SYSMIS };
1616           struct pivot_value *pv
1617             = (j == FRQ_ST_MODE || j == FRQ_ST_MINIMUM || j == FRQ_ST_MAXIMUM
1618                ? pivot_value_new_var_value (vf->var, &v)
1619                : pivot_value_new_number (v.f));
1620           pivot_table_put2 (table, var_idx, row++, pv);
1621         }
1622
1623       for (size_t j = 0; j < frq->n_percentiles; j++)
1624         {
1625           const struct percentile *pc = &frq->percentiles[j];
1626           if (!pc->show)
1627             continue;
1628
1629           union value v = {
1630             .f = vf->tab.n_valid ? vf->percentiles[j] : SYSMIS
1631           };
1632           pivot_table_put2 (table, var_idx, row++,
1633                             pivot_value_new_var_value (vf->var, &v));
1634         }
1635     }
1636
1637   pivot_table_submit (table);
1638 }