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

   1 /*
   2   PSPP - a program for statistical analysis.
   3   Copyright (C) 2012, 2013, 2016, 2019  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
  21 #include <math.h>
  22 #include <gsl/gsl_cdf.h>
  23
  24 #include "libpspp/assertion.h"
  25 #include "libpspp/message.h"
  26 #include "libpspp/pool.h"
  27
  28
  29 #include "data/dataset.h"
  30 #include "data/dictionary.h"
  31 #include "data/casegrouper.h"
  32 #include "data/casereader.h"
  33 #include "data/casewriter.h"
  34 #include "data/caseproto.h"
  35 #include "data/subcase.h"
  36
  37
  38 #include "data/format.h"
  39
  40 #include "math/interaction.h"
  41 #include "math/box-whisker.h"
  42 #include "math/categoricals.h"
  43 #include "math/chart-geometry.h"
  44 #include "math/histogram.h"
  45 #include "math/moments.h"
  46 #include "math/np.h"
  47 #include "math/sort.h"
  48 #include "math/order-stats.h"
  49 #include "math/percentiles.h"
  50 #include "math/shapiro-wilk.h"
  51 #include "math/tukey-hinges.h"
  52 #include "math/trimmed-mean.h"
  53
  54 #include "output/charts/boxplot.h"
  55 #include "output/charts/np-plot.h"
  56 #include "output/charts/spreadlevel-plot.h"
  57 #include "output/charts/plot-hist.h"
  58
  59 #include "language/command.h"
  60 #include "language/lexer/lexer.h"
  61 #include "language/lexer/value-parser.h"
  62 #include "language/lexer/variable-parser.h"
  63
  64 #include "output/pivot-table.h"
  65
  66 #include "gettext.h"
  67 #define _(msgid) gettext (msgid)
  68 #define N_(msgid) msgid
  69
  70 static void
  71 append_value_name (const struct variable *var, const union value *val, struct string *str)
  72 {
  73   var_append_value_name (var, val, str);
  74   if (var_is_value_missing (var, val))
  75     ds_put_cstr (str, _(" (missing)"));
  76 }
  77
  78 enum bp_mode
  79   {
  80     BP_GROUPS,
  81     BP_VARIABLES
  82   };
  83
  84
  85 /* Indices for the ex_proto member (below) */
  86 enum
  87   {
  88     EX_VAL,  /* value */
  89     EX_ID,   /* identity */
  90     EX_WT    /* weight */
  91   };
  92
  93
  94 #define PLOT_HISTOGRAM      0x1
  95 #define PLOT_BOXPLOT        0x2
  96 #define PLOT_NPPLOT         0x4
  97 #define PLOT_SPREADLEVEL    0x8
  98
  99 struct examine
 100 {
 101   struct pool *pool;
 102
 103   /* A caseproto used to contain the data subsets under examination,
 104      see (enum above)   */
 105   struct caseproto *ex_proto;
 106
 107   size_t n_dep_vars;
 108   const struct variable **dep_vars;
 109
 110   size_t n_iacts;
 111   struct interaction **iacts;
 112
 113   enum mv_class dep_excl;
 114   enum mv_class fctr_excl;
 115
 116   const struct dictionary *dict;
 117
 118   struct categoricals *cats;
 119
 120   /* how many extremities to display */
 121   int disp_extremes;
 122   int calc_extremes;
 123   bool descriptives;
 124
 125   double conf;
 126
 127   bool missing_pw;
 128
 129   /* The case index of the ID value (or -1) if not applicable */
 130   size_t id_idx;
 131   int id_width;
 132
 133   enum pc_alg pc_alg;
 134   double *ptiles;
 135   size_t n_percentiles;
 136
 137   unsigned int plot;
 138   float sl_power;
 139
 140   enum bp_mode boxplot_mode;
 141
 142   const struct variable *id_var;
 143
 144   const struct variable *wv;
 145 };
 146
 147 struct extremity
 148 {
 149   /* The value of this extremity */
 150   double val;
 151
 152   /* Either the casenumber or the value of the variable specified
 153      by the /ID subcommand which corresponds to this extremity */
 154   union value identity;
 155 };
 156
 157 struct exploratory_stats
 158 {
 159   double missing;
 160   double non_missing;
 161
 162   struct moments *mom;
 163
 164   /* Most operations need a sorted reader/writer */
 165   struct casewriter *sorted_writer;
 166   struct casereader *sorted_reader;
 167
 168   struct extremity *minima;
 169   struct extremity *maxima;
 170
 171   /*
 172      Minimum should alway equal mimima[0].val.
 173      Likewise, maximum should alway equal maxima[0].val.
 174      This redundancy exists as an optimisation effort.
 175      Some statistics (eg histogram) require early calculation
 176      of the min and max
 177   */
 178   double minimum;
 179   double maximum;
 180
 181   struct trimmed_mean *trimmed_mean;
 182   struct percentile *quartiles[3];
 183   struct percentile **percentiles;
 184   struct shapiro_wilk *shapiro_wilk;
 185
 186   struct tukey_hinges *hinges;
 187
 188   /* The data for the NP Plots */
 189   struct np *np;
 190
 191   struct histogram *histogram;
 192
 193   /* The data for the box plots */
 194   struct box_whisker *box_whisker;
 195
 196   /* Total weight */
 197   double cc;
 198
 199   /* The minimum weight */
 200   double cmin;
 201 };
 202
 203 static void
 204 show_boxplot_grouped (const struct examine *cmd, int iact_idx)
 205 {
 206   int v;
 207
 208   const struct interaction *iact = cmd->iacts[iact_idx];
 209   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 210
 211   for (v = 0; v < cmd->n_dep_vars; ++v)
 212     {
 213       double y_min = DBL_MAX;
 214       double y_max = -DBL_MAX;
 215       int grp;
 216       struct boxplot *boxplot;
 217       struct string title;
 218       ds_init_empty (&title);
 219
 220       if (iact->n_vars > 0)
 221         {
 222           struct string istr;
 223           ds_init_empty (&istr);
 224           interaction_to_string (iact, &istr);
 225           ds_put_format (&title, _("Boxplot of %s vs. %s"),
 226                          var_to_string (cmd->dep_vars[v]),
 227                          ds_cstr (&istr));
 228           ds_destroy (&istr);
 229         }
 230       else
 231         ds_put_format (&title, _("Boxplot of %s"), var_to_string (cmd->dep_vars[v]));
 232
 233       for (grp = 0; grp < n_cats; ++grp)
 234         {
 235           const struct exploratory_stats *es =
 236             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 237
 238           if (y_min > es[v].minimum)
 239             y_min = es[v].minimum;
 240
 241           if (y_max < es[v].maximum)
 242             y_max = es[v].maximum;
 243         }
 244
 245       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
 246
 247       ds_destroy (&title);
 248
 249       for (grp = 0; grp < n_cats; ++grp)
 250         {
 251           int ivar_idx;
 252           struct string label;
 253
 254           const struct ccase *c =
 255             categoricals_get_case_by_category_real (cmd->cats,  iact_idx, grp);
 256
 257           struct exploratory_stats *es =
 258             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 259
 260           ds_init_empty (&label);
 261           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 262             {
 263               struct string l;
 264               const struct variable *ivar = iact->vars[ivar_idx];
 265               const union value *val = case_data (c, ivar);
 266               ds_init_empty (&l);
 267
 268               append_value_name (ivar, val, &l);
 269               ds_ltrim (&l, ss_cstr (" "));
 270
 271               ds_put_substring (&label, l.ss);
 272               if (ivar_idx < iact->n_vars - 1)
 273                 ds_put_cstr (&label, "; ");
 274
 275               ds_destroy (&l);
 276             }
 277
 278           boxplot_add_box (boxplot, es[v].box_whisker, ds_cstr (&label));
 279           es[v].box_whisker = NULL;
 280
 281           ds_destroy (&label);
 282         }
 283
 284       boxplot_submit (boxplot);
 285     }
 286 }
 287
 288 static void
 289 show_boxplot_variabled (const struct examine *cmd, int iact_idx)
 290 {
 291   int grp;
 292   const struct interaction *iact = cmd->iacts[iact_idx];
 293   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 294
 295   for (grp = 0; grp < n_cats; ++grp)
 296     {
 297       struct boxplot *boxplot;
 298       int v;
 299       double y_min = DBL_MAX;
 300       double y_max = -DBL_MAX;
 301
 302       const struct ccase *c =
 303         categoricals_get_case_by_category_real (cmd->cats,  iact_idx, grp);
 304
 305       struct string title;
 306       ds_init_empty (&title);
 307
 308       for (v = 0; v < cmd->n_dep_vars; ++v)
 309         {
 310           const struct exploratory_stats *es =
 311             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 312
 313           if (y_min > es[v].minimum)
 314             y_min = es[v].minimum;
 315
 316           if (y_max < es[v].maximum)
 317             y_max = es[v].maximum;
 318         }
 319
 320       if (iact->n_vars == 0)
 321         ds_put_format (&title, _("Boxplot"));
 322       else
 323         {
 324           int ivar_idx;
 325           struct string label;
 326           ds_init_empty (&label);
 327           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 328             {
 329               const struct variable *ivar = iact->vars[ivar_idx];
 330               const union value *val = case_data (c, ivar);
 331
 332               ds_put_cstr (&label, var_to_string (ivar));
 333               ds_put_cstr (&label, " = ");
 334               append_value_name (ivar, val, &label);
 335               ds_put_cstr (&label, "; ");
 336             }
 337
 338           ds_put_format (&title, _("Boxplot of %s"),
 339                          ds_cstr (&label));
 340
 341           ds_destroy (&label);
 342         }
 343
 344       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
 345
 346       ds_destroy (&title);
 347
 348       for (v = 0; v < cmd->n_dep_vars; ++v)
 349         {
 350           struct exploratory_stats *es =
 351             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 352
 353           boxplot_add_box (boxplot, es[v].box_whisker,
 354                            var_to_string (cmd->dep_vars[v]));
 355           es[v].box_whisker = NULL;
 356         }
 357
 358       boxplot_submit (boxplot);
 359     }
 360 }
 361
 362
 363 static void
 364 show_npplot (const struct examine *cmd, int iact_idx)
 365 {
 366   const struct interaction *iact = cmd->iacts[iact_idx];
 367   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 368
 369   int v;
 370
 371   for (v = 0; v < cmd->n_dep_vars; ++v)
 372     {
 373       int grp;
 374       for (grp = 0; grp < n_cats; ++grp)
 375         {
 376           struct chart *npp, *dnpp;
 377           struct casereader *reader;
 378           struct np *np;
 379
 380           int ivar_idx;
 381           const struct ccase *c =
 382             categoricals_get_case_by_category_real (cmd->cats,
 383                                                     iact_idx, grp);
 384
 385           const struct exploratory_stats *es =
 386             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 387
 388           struct string label;
 389           ds_init_cstr (&label,
 390                         var_to_string (cmd->dep_vars[v]));
 391
 392           if (iact->n_vars > 0)
 393             {
 394               ds_put_cstr (&label, " (");
 395               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 396                 {
 397                   const struct variable *ivar = iact->vars[ivar_idx];
 398                   const union value *val = case_data (c, ivar);
 399
 400                   ds_put_cstr (&label, var_to_string (ivar));
 401                   ds_put_cstr (&label, " = ");
 402                   append_value_name (ivar, val, &label);
 403                   ds_put_cstr (&label, "; ");
 404
 405                 }
 406               ds_put_cstr (&label, ")");
 407             }
 408
 409           np = es[v].np;
 410           reader = casewriter_make_reader (np->writer);
 411           np->writer = NULL;
 412
 413           npp = np_plot_create (np, reader, ds_cstr (&label));
 414           dnpp = dnp_plot_create (np, reader, ds_cstr (&label));
 415
 416           if (npp == NULL || dnpp == NULL)
 417             {
 418               msg (MW, _("Not creating NP plot because data set is empty."));
 419               chart_unref (npp);
 420               chart_unref (dnpp);
 421             }
 422           else
 423             {
 424               chart_submit (npp);
 425               chart_submit (dnpp);
 426             }
 427           casereader_destroy (reader);
 428
 429           ds_destroy (&label);
 430         }
 431     }
 432 }
 433
 434 static void
 435 show_spreadlevel (const struct examine *cmd, int iact_idx)
 436 {
 437   const struct interaction *iact = cmd->iacts[iact_idx];
 438   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 439
 440   int v;
 441
 442   /* Spreadlevel when there are no levels is not useful */
 443   if (iact->n_vars == 0)
 444     return;
 445
 446   for (v = 0; v < cmd->n_dep_vars; ++v)
 447     {
 448       int grp;
 449       struct chart *sl;
 450
 451       struct string label;
 452       ds_init_cstr (&label,
 453                     var_to_string (cmd->dep_vars[v]));
 454
 455       if (iact->n_vars > 0)
 456         {
 457           ds_put_cstr (&label, " (");
 458           interaction_to_string (iact, &label);
 459           ds_put_cstr (&label, ")");
 460         }
 461
 462       sl = spreadlevel_plot_create (ds_cstr (&label), cmd->sl_power);
 463
 464       for (grp = 0; grp < n_cats; ++grp)
 465         {
 466           const struct exploratory_stats *es =
 467             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 468
 469           double median = percentile_calculate (es[v].quartiles[1], cmd->pc_alg);
 470
 471           double iqr = percentile_calculate (es[v].quartiles[2], cmd->pc_alg) -
 472             percentile_calculate (es[v].quartiles[0], cmd->pc_alg);
 473
 474           spreadlevel_plot_add (sl, iqr, median);
 475         }
 476
 477       if (sl == NULL)
 478         msg (MW, _("Not creating spreadlevel chart for %s"), ds_cstr (&label));
 479       else
 480         chart_submit (sl);
 481
 482       ds_destroy (&label);
 483     }
 484 }
 485
 486
 487 static void
 488 show_histogram (const struct examine *cmd, int iact_idx)
 489 {
 490   const struct interaction *iact = cmd->iacts[iact_idx];
 491   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 492
 493   int v;
 494
 495   for (v = 0; v < cmd->n_dep_vars; ++v)
 496     {
 497       int grp;
 498       for (grp = 0; grp < n_cats; ++grp)
 499         {
 500           double n, mean, var;
 501           int ivar_idx;
 502           const struct ccase *c =
 503             categoricals_get_case_by_category_real (cmd->cats,
 504                                                     iact_idx, grp);
 505
 506           const struct exploratory_stats *es =
 507             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 508
 509           struct string label;
 510
 511           if (es[v].histogram == NULL)
 512             continue;
 513
 514           ds_init_cstr (&label,
 515                         var_to_string (cmd->dep_vars[v]));
 516
 517           if (iact->n_vars > 0)
 518             {
 519               ds_put_cstr (&label, " (");
 520               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 521                 {
 522                   const struct variable *ivar = iact->vars[ivar_idx];
 523                   const union value *val = case_data (c, ivar);
 524
 525                   ds_put_cstr (&label, var_to_string (ivar));
 526                   ds_put_cstr (&label, " = ");
 527                   append_value_name (ivar, val, &label);
 528                   ds_put_cstr (&label, "; ");
 529
 530                 }
 531               ds_put_cstr (&label, ")");
 532             }
 533
 534
 535           moments_calculate (es[v].mom, &n, &mean, &var, NULL, NULL);
 536
 537           chart_submit
 538             (histogram_chart_create (es[v].histogram->gsl_hist,
 539                                       ds_cstr (&label), n, mean,
 540                                       sqrt (var), false));
 541
 542
 543           ds_destroy (&label);
 544         }
 545     }
 546 }
 547
 548 static struct pivot_value *
 549 new_value_with_missing_footnote (const struct variable *var,
 550                                  const union value *value,
 551                                  struct pivot_footnote *missing_footnote)
 552 {
 553   struct pivot_value *pv = pivot_value_new_var_value (var, value);
 554   if (var_is_value_missing (var, value) == MV_USER)
 555     pivot_value_add_footnote (pv, missing_footnote);
 556   return pv;
 557 }
 558
 559 static void
 560 create_interaction_dimensions (struct pivot_table *table,
 561                                const struct categoricals *cats,
 562                                const struct interaction *iact,
 563                                struct pivot_footnote *missing_footnote)
 564 {
 565   for (size_t i = iact->n_vars; i-- > 0;)
 566     {
 567       const struct variable *var = iact->vars[i];
 568       struct pivot_dimension *d = pivot_dimension_create__ (
 569         table, PIVOT_AXIS_ROW, pivot_value_new_variable (var));
 570       d->root->show_label = true;
 571
 572       size_t n;
 573       union value *values = categoricals_get_var_values (cats, var, &n);
 574       for (size_t j = 0; j < n; j++)
 575         pivot_category_create_leaf (
 576           d->root, new_value_with_missing_footnote (var, &values[j],
 577                                                     missing_footnote));
 578     }
 579 }
 580
 581 static struct pivot_footnote *
 582 create_missing_footnote (struct pivot_table *table)
 583 {
 584   return pivot_table_create_footnote (
 585     table, pivot_value_new_text (N_("User-missing value.")));
 586 }
 587
 588 static void
 589 percentiles_report (const struct examine *cmd, int iact_idx)
 590 {
 591   struct pivot_table *table = pivot_table_create (N_("Percentiles"));
 592
 593   struct pivot_dimension *percentiles = pivot_dimension_create (
 594     table, PIVOT_AXIS_COLUMN, N_("Percentiles"));
 595   percentiles->root->show_label = true;
 596   for (int i = 0; i < cmd->n_percentiles; ++i)
 597     pivot_category_create_leaf (
 598       percentiles->root,
 599       pivot_value_new_user_text_nocopy (xasprintf ("%g", cmd->ptiles[i])));
 600
 601   pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Statistics"),
 602                           N_("Weighted Average"), N_("Tukey's Hinges"));
 603
 604   const struct interaction *iact = cmd->iacts[iact_idx];
 605   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 606   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 607
 608   struct pivot_dimension *dep_dim = pivot_dimension_create (
 609     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 610
 611   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 612
 613   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 614   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 615     {
 616       indexes[table->n_dimensions - 1] = pivot_category_create_leaf (
 617         dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 618
 619       for (size_t i = 0; i < n_cats; ++i)
 620         {
 621           for (size_t j = 0; j < iact->n_vars; j++)
 622             {
 623               int idx = categoricals_get_value_index_by_category_real (
 624                 cmd->cats, iact_idx, i, j);
 625               indexes[table->n_dimensions - 2 - j] = idx;
 626             }
 627
 628           const struct exploratory_stats *ess
 629             = categoricals_get_user_data_by_category_real (
 630               cmd->cats, iact_idx, i);
 631           const struct exploratory_stats *es = ess + v;
 632
 633           double hinges[3];
 634           tukey_hinges_calculate (es->hinges, hinges);
 635
 636           for (size_t pc_idx = 0; pc_idx < cmd->n_percentiles; ++pc_idx)
 637             {
 638               indexes[0] = pc_idx;
 639
 640               indexes[1] = 0;
 641               double value = percentile_calculate (es->percentiles[pc_idx],
 642                                                    cmd->pc_alg);
 643               pivot_table_put (table, indexes, table->n_dimensions,
 644                                pivot_value_new_number (value));
 645
 646               double hinge = (cmd->ptiles[pc_idx] == 25.0 ? hinges[0]
 647                               : cmd->ptiles[pc_idx] == 50.0 ? hinges[1]
 648                               : cmd->ptiles[pc_idx] == 75.0 ? hinges[2]
 649                               : SYSMIS);
 650               if (hinge != SYSMIS)
 651                 {
 652                   indexes[1] = 1;
 653                   pivot_table_put (table, indexes, table->n_dimensions,
 654                                    pivot_value_new_number (hinge));
 655                 }
 656             }
 657         }
 658
 659     }
 660   free (indexes);
 661
 662   pivot_table_submit (table);
 663 }
 664
 665 static void
 666 normality_report (const struct examine *cmd, int iact_idx)
 667 {
 668   struct pivot_table *table = pivot_table_create (N_("Tests of Normality"));
 669
 670   struct pivot_dimension *test =
 671     pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Shapiro-Wilk"),
 672                             N_("Statistic"),
 673                             N_("df"), PIVOT_RC_COUNT,
 674                             N_("Sig."));
 675
 676   test->root->show_label = true;
 677
 678   const struct interaction *iact = cmd->iacts[iact_idx];
 679   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 680   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 681
 682   struct pivot_dimension *dep_dim = pivot_dimension_create (
 683     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 684
 685   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 686
 687   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 688   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 689     {
 690       indexes[table->n_dimensions - 1] =
 691         pivot_category_create_leaf (dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 692
 693       for (size_t i = 0; i < n_cats; ++i)
 694         {
 695           indexes[1] = i;
 696
 697           const struct exploratory_stats *es
 698             = categoricals_get_user_data_by_category_real (
 699               cmd->cats, iact_idx, i);
 700
 701           struct shapiro_wilk *sw =  es[v].shapiro_wilk;
 702
 703           if (sw == NULL)
 704             continue;
 705
 706           double w = shapiro_wilk_calculate (sw);
 707
 708           int j = 0;
 709           indexes[0] = j;
 710
 711           pivot_table_put (table, indexes, table->n_dimensions,
 712                            pivot_value_new_number (w));
 713
 714           indexes[0] = ++j;
 715           pivot_table_put (table, indexes, table->n_dimensions,
 716                            pivot_value_new_number (sw->n));
 717
 718           indexes[0] = ++j;
 719           pivot_table_put (table, indexes, table->n_dimensions,
 720                            pivot_value_new_number (shapiro_wilk_significance (sw->n, w)));
 721         }
 722     }
 723
 724   free (indexes);
 725
 726   pivot_table_submit (table);
 727 }
 728
 729
 730 static void
 731 descriptives_report (const struct examine *cmd, int iact_idx)
 732 {
 733   struct pivot_table *table = pivot_table_create (N_("Descriptives"));
 734
 735   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Aspect"),
 736                           N_("Statistic"), N_("Std. Error"));
 737
 738   struct pivot_dimension *statistics = pivot_dimension_create (
 739     table, PIVOT_AXIS_ROW, N_("Statistics"), N_("Mean"));
 740   struct pivot_category *interval = pivot_category_create_group__ (
 741     statistics->root,
 742     pivot_value_new_text_format (N_("%g%% Confidence Interval for Mean"),
 743                                  cmd->conf * 100.0));
 744   pivot_category_create_leaves (interval, N_("Lower Bound"),
 745                                 N_("Upper Bound"));
 746   pivot_category_create_leaves (
 747     statistics->root, N_("5% Trimmed Mean"), N_("Median"), N_("Variance"),
 748     N_("Std. Deviation"), N_("Minimum"), N_("Maximum"), N_("Range"),
 749     N_("Interquartile Range"), N_("Skewness"), N_("Kurtosis"));
 750
 751   const struct interaction *iact = cmd->iacts[iact_idx];
 752   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 753   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 754
 755   struct pivot_dimension *dep_dim = pivot_dimension_create (
 756     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 757
 758   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 759
 760   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 761   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 762     {
 763       indexes[table->n_dimensions - 1] = pivot_category_create_leaf (
 764         dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 765
 766       for (size_t i = 0; i < n_cats; ++i)
 767         {
 768           for (size_t j = 0; j < iact->n_vars; j++)
 769             {
 770               int idx = categoricals_get_value_index_by_category_real (
 771                 cmd->cats, iact_idx, i, j);
 772               indexes[table->n_dimensions - 2 - j] = idx;
 773             }
 774
 775           const struct exploratory_stats *ess
 776             = categoricals_get_user_data_by_category_real (cmd->cats,
 777                                                            iact_idx, i);
 778           const struct exploratory_stats *es = ess + v;
 779
 780           double m0, m1, m2, m3, m4;
 781           moments_calculate (es->mom, &m0, &m1, &m2, &m3, &m4);
 782           double tval = gsl_cdf_tdist_Qinv ((1.0 - cmd->conf) / 2.0, m0 - 1.0);
 783
 784           struct entry
 785             {
 786               int stat_idx;
 787               int aspect_idx;
 788               double x;
 789             }
 790           entries[] = {
 791             { 0, 0, m1 },
 792             { 0, 1, calc_semean (m2, m0) },
 793             { 1, 0, m1 - tval * calc_semean (m2, m0) },
 794             { 2, 0, m1 + tval * calc_semean (m2, m0) },
 795             { 3, 0, trimmed_mean_calculate (es->trimmed_mean) },
 796             { 4, 0, percentile_calculate (es->quartiles[1], cmd->pc_alg) },
 797             { 5, 0, m2 },
 798             { 6, 0, sqrt (m2) },
 799             { 7, 0, es->minima[0].val },
 800             { 8, 0, es->maxima[0].val },
 801             { 9, 0, es->maxima[0].val - es->minima[0].val },
 802             { 10, 0, (percentile_calculate (es->quartiles[2], cmd->pc_alg) -
 803                       percentile_calculate (es->quartiles[0], cmd->pc_alg)) },
 804             { 11, 0, m3 },
 805             { 11, 1, calc_seskew (m0) },
 806             { 12, 0, m4 },
 807             { 12, 1, calc_sekurt (m0) },
 808           };
 809           for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 810             {
 811               const struct entry *e = &entries[j];
 812               indexes[0] = e->aspect_idx;
 813               indexes[1] = e->stat_idx;
 814               pivot_table_put (table, indexes, table->n_dimensions,
 815                                pivot_value_new_number (e->x));
 816             }
 817         }
 818     }
 819
 820   free (indexes);
 821
 822   pivot_table_submit (table);
 823 }
 824
 825
 826 static void
 827 extremes_report (const struct examine *cmd, int iact_idx)
 828 {
 829   struct pivot_table *table = pivot_table_create (N_("Extreme Values"));
 830
 831   struct pivot_dimension *statistics = pivot_dimension_create (
 832     table, PIVOT_AXIS_COLUMN, N_("Statistics"));
 833   pivot_category_create_leaf (statistics->root,
 834                               (cmd->id_var
 835                                ? pivot_value_new_variable (cmd->id_var)
 836                                : pivot_value_new_text (N_("Case Number"))));
 837   pivot_category_create_leaves (statistics->root, N_("Value"));
 838
 839   struct pivot_dimension *order = pivot_dimension_create (
 840     table, PIVOT_AXIS_ROW, N_("Order"));
 841   for (size_t i = 0; i < cmd->disp_extremes; i++)
 842     pivot_category_create_leaf (order->root, pivot_value_new_integer (i + 1));
 843
 844   pivot_dimension_create (table, PIVOT_AXIS_ROW,
 845                           /* TRANSLATORS: This is a noun, not an adjective.  */
 846                           N_("Extreme"),
 847                           N_("Highest"), N_("Lowest"));
 848
 849   const struct interaction *iact = cmd->iacts[iact_idx];
 850   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 851   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 852
 853   struct pivot_dimension *dep_dim = pivot_dimension_create (
 854     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 855
 856   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 857
 858   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 859   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 860     {
 861       indexes[table->n_dimensions - 1] = pivot_category_create_leaf (
 862         dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 863
 864       for (size_t i = 0; i < n_cats; ++i)
 865         {
 866           for (size_t j = 0; j < iact->n_vars; j++)
 867             {
 868               int idx = categoricals_get_value_index_by_category_real (
 869                 cmd->cats, iact_idx, i, j);
 870               indexes[table->n_dimensions - 2 - j] = idx;
 871             }
 872
 873           const struct exploratory_stats *ess
 874             = categoricals_get_user_data_by_category_real (cmd->cats,
 875                                                            iact_idx, i);
 876           const struct exploratory_stats *es = ess + v;
 877
 878           for (int e = 0 ; e < cmd->disp_extremes; ++e)
 879             {
 880               indexes[1] = e;
 881
 882               for (size_t j = 0; j < 2; j++)
 883                 {
 884                   const struct extremity *extremity
 885                     = j ? &es->minima[e] : &es->maxima[e];
 886                   indexes[2] = j;
 887
 888                   indexes[0] = 0;
 889                   pivot_table_put (
 890                     table, indexes, table->n_dimensions,
 891                     (cmd->id_var
 892                      ? new_value_with_missing_footnote (cmd->id_var,
 893                                                         &extremity->identity,
 894                                                         missing_footnote)
 895                      : pivot_value_new_integer (extremity->identity.f)));
 896
 897                   indexes[0] = 1;
 898                   union value val = { .f = extremity->val };
 899                   pivot_table_put (
 900                     table, indexes, table->n_dimensions,
 901                     new_value_with_missing_footnote (cmd->dep_vars[v], &val,
 902                                                      missing_footnote));
 903                 }
 904             }
 905         }
 906     }
 907   free (indexes);
 908
 909   pivot_table_submit (table);
 910 }
 911
 912
 913 static void
 914 summary_report (const struct examine *cmd, int iact_idx)
 915 {
 916   struct pivot_table *table = pivot_table_create (
 917     N_("Case Processing Summary"));
 918   pivot_table_set_weight_var (table, dict_get_weight (cmd->dict));
 919
 920   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 921                           N_("N"), PIVOT_RC_COUNT,
 922                           N_("Percent"), PIVOT_RC_PERCENT);
 923   struct pivot_dimension *cases = pivot_dimension_create (
 924     table, PIVOT_AXIS_COLUMN, N_("Cases"), N_("Valid"), N_("Missing"),
 925     N_("Total"));
 926   cases->root->show_label = true;
 927
 928   const struct interaction *iact = cmd->iacts[iact_idx];
 929   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 930   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 931
 932   struct pivot_dimension *dep_dim = pivot_dimension_create (
 933     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 934
 935   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 936
 937   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 938   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 939     {
 940       indexes[table->n_dimensions - 1] = pivot_category_create_leaf (
 941         dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 942
 943       for (size_t i = 0; i < n_cats; ++i)
 944         {
 945           for (size_t j = 0; j < iact->n_vars; j++)
 946             {
 947               int idx = categoricals_get_value_index_by_category_real (
 948                 cmd->cats, iact_idx, i, j);
 949               indexes[table->n_dimensions - 2 - j] = idx;
 950             }
 951
 952           const struct exploratory_stats *es
 953             = categoricals_get_user_data_by_category_real (
 954               cmd->cats, iact_idx, i);
 955
 956           double total = es[v].missing + es[v].non_missing;
 957           struct entry
 958             {
 959               int stat_idx;
 960               int case_idx;
 961               double x;
 962             }
 963           entries[] = {
 964             { 0, 0, es[v].non_missing },
 965             { 1, 0, 100.0 * es[v].non_missing / total },
 966             { 0, 1, es[v].missing },
 967             { 1, 1, 100.0 * es[v].missing / total },
 968             { 0, 2, total },
 969             { 1, 2, 100.0 },
 970           };
 971           for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 972             {
 973               const struct entry *e = &entries[j];
 974               indexes[0] = e->stat_idx;
 975               indexes[1] = e->case_idx;
 976               pivot_table_put (table, indexes, table->n_dimensions,
 977                                pivot_value_new_number (e->x));
 978             }
 979         }
 980     }
 981
 982   free (indexes);
 983
 984   pivot_table_submit (table);
 985 }
 986
 987 /* Attempt to parse an interaction from LEXER */
 988 static struct interaction *
 989 parse_interaction (struct lexer *lexer, struct examine *ex)
 990 {
 991   const struct variable *v = NULL;
 992   struct interaction *iact = NULL;
 993
 994   if (lex_match_variable (lexer, ex->dict, &v))
 995     {
 996       iact = interaction_create (v);
 997
 998       while (lex_match (lexer, T_BY))
 999         {
1000           if (!lex_match_variable (lexer, ex->dict, &v))
1001             {
1002               interaction_destroy (iact);
1003               return NULL;
1004             }
1005           interaction_add_variable (iact, v);
1006         }
1007       lex_match (lexer, T_COMMA);
1008     }
1009
1010   return iact;
1011 }
1012
1013
1014 static void *
1015 create_n (const void *aux1, void *aux2 UNUSED)
1016 {
1017   int v;
1018
1019   const struct examine *examine = aux1;
1020   struct exploratory_stats *es = pool_calloc (examine->pool, examine->n_dep_vars, sizeof (*es));
1021   struct subcase ordering;
1022   subcase_init (&ordering, 0, 0, SC_ASCEND);
1023
1024   for (v = 0; v < examine->n_dep_vars; v++)
1025     {
1026       es[v].sorted_writer = sort_create_writer (&ordering, examine->ex_proto);
1027       es[v].sorted_reader = NULL;
1028
1029       es[v].mom = moments_create (MOMENT_KURTOSIS);
1030       es[v].cmin = DBL_MAX;
1031
1032       es[v].maximum = -DBL_MAX;
1033       es[v].minimum =  DBL_MAX;
1034     }
1035
1036   subcase_destroy (&ordering);
1037   return es;
1038 }
1039
1040 static void
1041 update_n (const void *aux1, void *aux2 UNUSED, void *user_data,
1042           const struct ccase *c, double weight)
1043 {
1044   int v;
1045   const struct examine *examine = aux1;
1046   struct exploratory_stats *es = user_data;
1047
1048   bool this_case_is_missing = false;
1049   /* LISTWISE missing must be dealt with here */
1050   if (!examine->missing_pw)
1051     {
1052       for (v = 0; v < examine->n_dep_vars; v++)
1053         {
1054           const struct variable *var = examine->dep_vars[v];
1055
1056           if (var_is_value_missing (var, case_data (c, var))
1057               & examine->dep_excl)
1058             {
1059               es[v].missing += weight;
1060               this_case_is_missing = true;
1061             }
1062         }
1063     }
1064
1065   if (this_case_is_missing)
1066     return;
1067
1068   for (v = 0; v < examine->n_dep_vars; v++)
1069     {
1070       struct ccase *outcase ;
1071       const struct variable *var = examine->dep_vars[v];
1072       const double x = case_num (c, var);
1073
1074       if (var_is_value_missing (var, case_data (c, var)) & examine->dep_excl)
1075         {
1076           es[v].missing += weight;
1077           continue;
1078         }
1079
1080       outcase = case_create (examine->ex_proto);
1081
1082       if (x > es[v].maximum)
1083         es[v].maximum = x;
1084
1085       if (x < es[v].minimum)
1086         es[v].minimum =  x;
1087
1088       es[v].non_missing += weight;
1089
1090       moments_pass_one (es[v].mom, x, weight);
1091
1092       /* Save the value and the ID to the writer */
1093       assert (examine->id_idx != -1);
1094       *case_num_rw_idx (outcase, EX_VAL) = x;
1095       value_copy (case_data_rw_idx (outcase, EX_ID),
1096                   case_data_idx (c, examine->id_idx), examine->id_width);
1097
1098       *case_num_rw_idx (outcase, EX_WT) = weight;
1099
1100       es[v].cc += weight;
1101
1102       if (es[v].cmin > weight)
1103         es[v].cmin = weight;
1104
1105       casewriter_write (es[v].sorted_writer, outcase);
1106     }
1107 }
1108
1109 static void
1110 calculate_n (const void *aux1, void *aux2 UNUSED, void *user_data)
1111 {
1112   int v;
1113   const struct examine *examine = aux1;
1114   struct exploratory_stats *es = user_data;
1115
1116   for (v = 0; v < examine->n_dep_vars; v++)
1117     {
1118       int i;
1119       casenumber imin = 0;
1120       casenumber imax;
1121       struct casereader *reader;
1122       struct ccase *c;
1123
1124       if (examine->plot & PLOT_HISTOGRAM && es[v].non_missing > 0)
1125         {
1126           /* Sturges Rule */
1127           double bin_width = fabs (es[v].minimum - es[v].maximum)
1128             / (1 + log2 (es[v].cc))
1129             ;
1130
1131           es[v].histogram =
1132             histogram_create (bin_width, es[v].minimum, es[v].maximum);
1133         }
1134
1135       es[v].sorted_reader = casewriter_make_reader (es[v].sorted_writer);
1136       es[v].sorted_writer = NULL;
1137
1138       imax = casereader_get_n_cases (es[v].sorted_reader);
1139
1140       es[v].maxima = pool_calloc (examine->pool, examine->calc_extremes, sizeof (*es[v].maxima));
1141       es[v].minima = pool_calloc (examine->pool, examine->calc_extremes, sizeof (*es[v].minima));
1142       for (i = 0; i < examine->calc_extremes; ++i)
1143         {
1144           value_init_pool (examine->pool, &es[v].maxima[i].identity, examine->id_width) ;
1145           value_init_pool (examine->pool, &es[v].minima[i].identity, examine->id_width) ;
1146         }
1147
1148       bool warn = true;
1149       for (reader = casereader_clone (es[v].sorted_reader);
1150            (c = casereader_read (reader)) != NULL; case_unref (c))
1151         {
1152           const double val = case_num_idx (c, EX_VAL);
1153           double wt = case_num_idx (c, EX_WT);
1154           wt = var_force_valid_weight (examine->wv, wt, &warn);
1155
1156           moments_pass_two (es[v].mom, val, wt);
1157
1158           if (es[v].histogram)
1159             histogram_add (es[v].histogram, val, wt);
1160
1161           if (imin < examine->calc_extremes)
1162             {
1163               int x;
1164               for (x = imin; x < examine->calc_extremes; ++x)
1165                 {
1166                   struct extremity *min = &es[v].minima[x];
1167                   min->val = val;
1168                   value_copy (&min->identity, case_data_idx (c, EX_ID), examine->id_width);
1169                 }
1170               imin ++;
1171             }
1172
1173           imax --;
1174           if (imax < examine->calc_extremes)
1175             {
1176               int x;
1177
1178               for (x = imax; x < imax + 1; ++x)
1179                 {
1180                   struct extremity *max;
1181
1182                   if (x >= examine->calc_extremes)
1183                     break;
1184
1185                   max = &es[v].maxima[x];
1186                   max->val = val;
1187                   value_copy (&max->identity, case_data_idx (c, EX_ID), examine->id_width);
1188                 }
1189             }
1190         }
1191       casereader_destroy (reader);
1192
1193       if (examine->calc_extremes > 0 && es[v].non_missing > 0)
1194         {
1195           assert (es[v].minima[0].val == es[v].minimum);
1196           assert (es[v].maxima[0].val == es[v].maximum);
1197         }
1198
1199       {
1200         const int n_os = 5 + examine->n_percentiles;
1201         es[v].percentiles = pool_calloc (examine->pool, examine->n_percentiles, sizeof (*es[v].percentiles));
1202
1203         es[v].trimmed_mean = trimmed_mean_create (es[v].cc, 0.05);
1204         es[v].shapiro_wilk = NULL;
1205
1206         struct order_stats **os = XCALLOC (n_os, struct order_stats *);
1207         os[0] = &es[v].trimmed_mean->parent;
1208
1209         es[v].quartiles[0] = percentile_create (0.25, es[v].cc);
1210         es[v].quartiles[1] = percentile_create (0.5,  es[v].cc);
1211         es[v].quartiles[2] = percentile_create (0.75, es[v].cc);
1212
1213         os[1] = &es[v].quartiles[0]->parent;
1214         os[2] = &es[v].quartiles[1]->parent;
1215         os[3] = &es[v].quartiles[2]->parent;
1216
1217         es[v].hinges = tukey_hinges_create (es[v].cc, es[v].cmin);
1218         os[4] = &es[v].hinges->parent;
1219
1220         for (i = 0; i < examine->n_percentiles; ++i)
1221           {
1222             es[v].percentiles[i] = percentile_create (examine->ptiles[i] / 100.00, es[v].cc);
1223             os[5 + i] = &es[v].percentiles[i]->parent;
1224           }
1225
1226         order_stats_accumulate_idx (os, n_os,
1227                                     casereader_clone (es[v].sorted_reader),
1228                                     EX_WT, EX_VAL);
1229
1230         free (os);
1231       }
1232
1233       if (examine->plot & PLOT_BOXPLOT)
1234         {
1235           struct order_stats *os;
1236
1237           es[v].box_whisker = box_whisker_create (es[v].hinges,
1238                                                   EX_ID, examine->id_var);
1239
1240           os = &es[v].box_whisker->parent;
1241           order_stats_accumulate_idx (&os, 1,
1242                                       casereader_clone (es[v].sorted_reader),
1243                                       EX_WT, EX_VAL);
1244         }
1245
1246       if (examine->plot)
1247         {
1248           double mean;
1249
1250           moments_calculate (es[v].mom, NULL, &mean, NULL, NULL, NULL);
1251
1252           es[v].shapiro_wilk = shapiro_wilk_create (es[v].non_missing, mean);
1253
1254           if (es[v].shapiro_wilk)
1255             {
1256               struct order_stats *os = &es[v].shapiro_wilk->parent;
1257               order_stats_accumulate_idx (&os, 1,
1258                                           casereader_clone (es[v].sorted_reader),
1259                                           EX_WT, EX_VAL);
1260             }
1261         }
1262
1263       if (examine->plot & PLOT_NPPLOT)
1264         {
1265           double n, mean, var;
1266           struct order_stats *os;
1267
1268           moments_calculate (es[v].mom, &n, &mean, &var, NULL, NULL);
1269
1270           es[v].np = np_create (n, mean, var);
1271
1272           os = &es[v].np->parent;
1273
1274           order_stats_accumulate_idx (&os, 1,
1275                                       casereader_clone (es[v].sorted_reader),
1276                                       EX_WT, EX_VAL);
1277         }
1278
1279     }
1280 }
1281
1282 static void
1283 cleanup_exploratory_stats (struct examine *cmd)
1284 {
1285   int i;
1286   for (i = 0; i < cmd->n_iacts; ++i)
1287     {
1288       int v;
1289       const size_t n_cats =  categoricals_n_count (cmd->cats, i);
1290
1291       for (v = 0; v < cmd->n_dep_vars; ++v)
1292         {
1293           int grp;
1294           for (grp = 0; grp < n_cats; ++grp)
1295             {
1296               int q;
1297               const struct exploratory_stats *es =
1298                 categoricals_get_user_data_by_category_real (cmd->cats, i, grp);
1299
1300               struct order_stats *os = &es[v].hinges->parent;
1301               struct statistic  *stat = &os->parent;
1302               stat->destroy (stat);
1303
1304               for (q = 0; q < 3 ; q++)
1305                 {
1306                   os = &es[v].quartiles[q]->parent;
1307                   stat = &os->parent;
1308                   stat->destroy (stat);
1309                 }
1310
1311               for (q = 0; q < cmd->n_percentiles ; q++)
1312                 {
1313                   os = &es[v].percentiles[q]->parent;
1314                   stat = &os->parent;
1315                   stat->destroy (stat);
1316                 }
1317
1318               if (es[v].shapiro_wilk)
1319                 {
1320                   stat = &es[v].shapiro_wilk->parent.parent;
1321                   stat->destroy (stat);
1322                 }
1323
1324               os = &es[v].trimmed_mean->parent;
1325               stat = &os->parent;
1326               stat->destroy (stat);
1327
1328               os = &es[v].np->parent;
1329               if (os)
1330                 {
1331                   stat = &os->parent;
1332                   stat->destroy (stat);
1333                 }
1334
1335               statistic_destroy (&es[v].histogram->parent);
1336               moments_destroy (es[v].mom);
1337
1338               if (es[v].box_whisker)
1339                 {
1340                   stat = &es[v].box_whisker->parent.parent;
1341                   stat->destroy (stat);
1342                 }
1343
1344               casereader_destroy (es[v].sorted_reader);
1345             }
1346         }
1347     }
1348 }
1349
1350
1351 static void
1352 run_examine (struct examine *cmd, struct casereader *input)
1353 {
1354   int i;
1355   struct ccase *c;
1356   struct casereader *reader;
1357
1358   struct payload payload;
1359   payload.create = create_n;
1360   payload.update = update_n;
1361   payload.calculate = calculate_n;
1362   payload.destroy = NULL;
1363
1364   cmd->wv = dict_get_weight (cmd->dict);
1365
1366   cmd->cats
1367     = categoricals_create (cmd->iacts, cmd->n_iacts, cmd->wv, cmd->fctr_excl);
1368
1369   categoricals_set_payload (cmd->cats, &payload, cmd, NULL);
1370
1371   if (cmd->id_var == NULL)
1372     {
1373       struct ccase *c = casereader_peek (input,  0);
1374
1375       cmd->id_idx = case_get_n_values (c);
1376       input = casereader_create_arithmetic_sequence (input, 1.0, 1.0);
1377
1378       case_unref (c);
1379     }
1380
1381   for (reader = input;
1382        (c = casereader_read (reader)) != NULL; case_unref (c))
1383     {
1384       categoricals_update (cmd->cats, c);
1385     }
1386   casereader_destroy (reader);
1387   categoricals_done (cmd->cats);
1388
1389   for (i = 0; i < cmd->n_iacts; ++i)
1390     {
1391       summary_report (cmd, i);
1392
1393       const size_t n_cats =  categoricals_n_count (cmd->cats, i);
1394       if (n_cats == 0)
1395         continue;
1396
1397       if (cmd->disp_extremes > 0)
1398         extremes_report (cmd, i);
1399
1400       if (cmd->n_percentiles > 0)
1401         percentiles_report (cmd, i);
1402
1403       if (cmd->plot & PLOT_BOXPLOT)
1404         {
1405           switch (cmd->boxplot_mode)
1406             {
1407             case BP_GROUPS:
1408               show_boxplot_grouped (cmd, i);
1409               break;
1410             case BP_VARIABLES:
1411               show_boxplot_variabled (cmd, i);
1412               break;
1413             default:
1414               NOT_REACHED ();
1415               break;
1416             }
1417         }
1418
1419       if (cmd->plot & PLOT_HISTOGRAM)
1420         show_histogram (cmd, i);
1421
1422       if (cmd->plot & PLOT_NPPLOT)
1423         show_npplot (cmd, i);
1424
1425       if (cmd->plot & PLOT_SPREADLEVEL)
1426         show_spreadlevel (cmd, i);
1427
1428       if (cmd->descriptives)
1429         descriptives_report (cmd, i);
1430
1431       if (cmd->plot)
1432         normality_report (cmd, i);
1433     }
1434
1435   cleanup_exploratory_stats (cmd);
1436   categoricals_destroy (cmd->cats);
1437 }
1438
1439
1440 int
1441 cmd_examine (struct lexer *lexer, struct dataset *ds)
1442 {
1443   int i;
1444   bool nototals_seen = false;
1445   bool totals_seen = false;
1446
1447   struct interaction **iacts_mem = NULL;
1448   struct examine examine;
1449   bool percentiles_seen = false;
1450
1451   examine.missing_pw = false;
1452   examine.disp_extremes = 0;
1453   examine.calc_extremes = 0;
1454   examine.descriptives = false;
1455   examine.conf = 0.95;
1456   examine.pc_alg = PC_HAVERAGE;
1457   examine.ptiles = NULL;
1458   examine.n_percentiles = 0;
1459   examine.id_idx = -1;
1460   examine.id_width = 0;
1461   examine.id_var = NULL;
1462   examine.boxplot_mode = BP_GROUPS;
1463
1464   examine.ex_proto = caseproto_create ();
1465
1466   examine.pool = pool_create ();
1467
1468   /* Allocate space for the first interaction.
1469      This is interaction is an empty one (for the totals).
1470      If no totals are requested, we will simply ignore this
1471      interaction.
1472   */
1473   examine.n_iacts = 1;
1474   examine.iacts = iacts_mem = pool_zalloc (examine.pool, sizeof (struct interaction *));
1475   examine.iacts[0] = interaction_create (NULL);
1476
1477   examine.dep_excl = MV_ANY;
1478   examine.fctr_excl = MV_ANY;
1479   examine.plot = 0;
1480   examine.sl_power = 0;
1481   examine.dep_vars = NULL;
1482   examine.n_dep_vars = 0;
1483   examine.dict = dataset_dict (ds);
1484
1485   /* Accept an optional, completely pointless "/VARIABLES=" */
1486   lex_match (lexer, T_SLASH);
1487   if (lex_match_id  (lexer, "VARIABLES"))
1488     {
1489       if (! lex_force_match (lexer, T_EQUALS))
1490         goto error;
1491     }
1492
1493   if (!parse_variables_const (lexer, examine.dict,
1494                               &examine.dep_vars, &examine.n_dep_vars,
1495                               PV_NO_DUPLICATE | PV_NUMERIC))
1496     goto error;
1497
1498   if (lex_match (lexer, T_BY))
1499     {
1500       struct interaction *iact = NULL;
1501       do
1502         {
1503           iact = parse_interaction (lexer, &examine);
1504           if (iact)
1505             {
1506               examine.n_iacts++;
1507               iacts_mem =
1508                 pool_nrealloc (examine.pool, iacts_mem,
1509                                examine.n_iacts,
1510                                sizeof (*iacts_mem));
1511
1512               iacts_mem[examine.n_iacts - 1] = iact;
1513             }
1514         }
1515       while (iact);
1516     }
1517
1518
1519   while (lex_token (lexer) != T_ENDCMD)
1520     {
1521       lex_match (lexer, T_SLASH);
1522
1523       if (lex_match_id (lexer, "STATISTICS"))
1524         {
1525           lex_match (lexer, T_EQUALS);
1526
1527           while (lex_token (lexer) != T_ENDCMD
1528                  && lex_token (lexer) != T_SLASH)
1529             {
1530               if (lex_match_id (lexer, "DESCRIPTIVES"))
1531                 {
1532                   examine.descriptives = true;
1533                 }
1534               else if (lex_match_id (lexer, "EXTREME"))
1535                 {
1536                   int extr = 5;
1537                   if (lex_match (lexer, T_LPAREN))
1538                     {
1539                       if (!lex_force_int_range (lexer, "EXTREME", 0, INT_MAX))
1540                         goto error;
1541                       extr = lex_integer (lexer);
1542
1543                       lex_get (lexer);
1544                       if (! lex_force_match (lexer, T_RPAREN))
1545                         goto error;
1546                     }
1547                   examine.disp_extremes  = extr;
1548                 }
1549               else if (lex_match_id (lexer, "NONE"))
1550                 {
1551                 }
1552               else if (lex_match (lexer, T_ALL))
1553                 {
1554                   if (examine.disp_extremes == 0)
1555                     examine.disp_extremes = 5;
1556                 }
1557               else
1558                 {
1559                   lex_error (lexer, NULL);
1560                   goto error;
1561                 }
1562             }
1563         }
1564       else if (lex_match_id (lexer, "PERCENTILES"))
1565         {
1566           percentiles_seen = true;
1567           if (lex_match (lexer, T_LPAREN))
1568             {
1569               while (lex_is_number (lexer))
1570                 {
1571                   if (!lex_force_num_range_open (lexer, "PERCENTILES", 0, 100))
1572                     goto error;
1573                   double p = lex_number (lexer);
1574
1575                   examine.n_percentiles++;
1576                   examine.ptiles =
1577                     xrealloc (examine.ptiles,
1578                               sizeof (*examine.ptiles) *
1579                               examine.n_percentiles);
1580
1581                   examine.ptiles[examine.n_percentiles - 1] = p;
1582
1583                   lex_get (lexer);
1584                   lex_match (lexer, T_COMMA);
1585                 }
1586               if (!lex_force_match (lexer, T_RPAREN))
1587                 goto error;
1588             }
1589
1590           lex_match (lexer, T_EQUALS);
1591
1592           while (lex_token (lexer) != T_ENDCMD
1593                  && lex_token (lexer) != T_SLASH)
1594             {
1595               if (lex_match_id (lexer, "HAVERAGE"))
1596                 {
1597                   examine.pc_alg = PC_HAVERAGE;
1598                 }
1599               else if (lex_match_id (lexer, "WAVERAGE"))
1600                 {
1601                   examine.pc_alg = PC_WAVERAGE;
1602                 }
1603               else if (lex_match_id (lexer, "ROUND"))
1604                 {
1605                   examine.pc_alg = PC_ROUND;
1606                 }
1607               else if (lex_match_id (lexer, "EMPIRICAL"))
1608                 {
1609                   examine.pc_alg = PC_EMPIRICAL;
1610                 }
1611               else if (lex_match_id (lexer, "AEMPIRICAL"))
1612                 {
1613                   examine.pc_alg = PC_AEMPIRICAL;
1614                 }
1615               else if (lex_match_id (lexer, "NONE"))
1616                 {
1617                   examine.pc_alg = PC_NONE;
1618                 }
1619               else
1620                 {
1621                   lex_error (lexer, NULL);
1622                   goto error;
1623                 }
1624             }
1625         }
1626       else if (lex_match_id (lexer, "TOTAL"))
1627         {
1628           totals_seen = true;
1629         }
1630       else if (lex_match_id (lexer, "NOTOTAL"))
1631         {
1632           nototals_seen = true;
1633         }
1634       else if (lex_match_id (lexer, "MISSING"))
1635         {
1636           lex_match (lexer, T_EQUALS);
1637
1638           while (lex_token (lexer) != T_ENDCMD
1639                  && lex_token (lexer) != T_SLASH)
1640             {
1641               if (lex_match_id (lexer, "LISTWISE"))
1642                 {
1643                   examine.missing_pw = false;
1644                 }
1645               else if (lex_match_id (lexer, "PAIRWISE"))
1646                 {
1647                   examine.missing_pw = true;
1648                 }
1649               else if (lex_match_id (lexer, "EXCLUDE"))
1650                 {
1651                   examine.dep_excl = MV_ANY;
1652                 }
1653               else if (lex_match_id (lexer, "INCLUDE"))
1654                 {
1655                   examine.dep_excl = MV_SYSTEM;
1656                 }
1657               else if (lex_match_id (lexer, "REPORT"))
1658                 {
1659                   examine.fctr_excl = 0;
1660                 }
1661               else if (lex_match_id (lexer, "NOREPORT"))
1662                 {
1663                   examine.fctr_excl = MV_ANY;
1664                 }
1665               else
1666                 {
1667                   lex_error (lexer, NULL);
1668                   goto error;
1669                 }
1670             }
1671         }
1672       else if (lex_match_id (lexer, "COMPARE"))
1673         {
1674           lex_match (lexer, T_EQUALS);
1675           if (lex_match_id (lexer, "VARIABLES"))
1676             {
1677               examine.boxplot_mode = BP_VARIABLES;
1678             }
1679           else if (lex_match_id (lexer, "GROUPS"))
1680             {
1681               examine.boxplot_mode = BP_GROUPS;
1682             }
1683           else
1684             {
1685               lex_error (lexer, NULL);
1686               goto error;
1687             }
1688         }
1689       else if (lex_match_id (lexer, "PLOT"))
1690         {
1691           lex_match (lexer, T_EQUALS);
1692
1693           while (lex_token (lexer) != T_ENDCMD
1694                  && lex_token (lexer) != T_SLASH)
1695             {
1696               if (lex_match_id (lexer, "BOXPLOT"))
1697                 {
1698                   examine.plot |= PLOT_BOXPLOT;
1699                 }
1700               else if (lex_match_id (lexer, "NPPLOT"))
1701                 {
1702                   examine.plot |= PLOT_NPPLOT;
1703                 }
1704               else if (lex_match_id (lexer, "HISTOGRAM"))
1705                 {
1706                   examine.plot |= PLOT_HISTOGRAM;
1707                 }
1708               else if (lex_match_id (lexer, "SPREADLEVEL"))
1709                 {
1710                   examine.plot |= PLOT_SPREADLEVEL;
1711                   examine.sl_power = 0;
1712                   if (lex_match (lexer, T_LPAREN) && lex_force_num (lexer))
1713                     {
1714                       examine.sl_power = lex_number (lexer);
1715
1716                       lex_get (lexer);
1717                       if (! lex_force_match (lexer, T_RPAREN))
1718                         goto error;
1719                     }
1720                 }
1721               else if (lex_match_id (lexer, "NONE"))
1722                 {
1723                   examine.plot = 0;
1724                 }
1725               else if (lex_match (lexer, T_ALL))
1726                 {
1727                   examine.plot = ~0;
1728                 }
1729               else
1730                 {
1731                   lex_error (lexer, NULL);
1732                   goto error;
1733                 }
1734               lex_match (lexer, T_COMMA);
1735             }
1736         }
1737       else if (lex_match_id (lexer, "CINTERVAL"))
1738         {
1739           if (!lex_force_num (lexer))
1740             goto error;
1741
1742           examine.conf = lex_number (lexer);
1743           lex_get (lexer);
1744         }
1745       else if (lex_match_id (lexer, "ID"))
1746         {
1747           lex_match (lexer, T_EQUALS);
1748
1749           examine.id_var = parse_variable_const (lexer, examine.dict);
1750         }
1751       else
1752         {
1753           lex_error (lexer, NULL);
1754           goto error;
1755         }
1756     }
1757
1758
1759   if (totals_seen && nototals_seen)
1760     {
1761       msg (SE, _("%s and %s are mutually exclusive"), "TOTAL", "NOTOTAL");
1762       goto error;
1763     }
1764
1765   /* If totals have been requested or if there are no factors
1766      in this analysis, then the totals need to be included. */
1767   if (!nototals_seen || examine.n_iacts == 1)
1768     {
1769       examine.iacts = &iacts_mem[0];
1770     }
1771   else
1772     {
1773       examine.n_iacts--;
1774       examine.iacts = &iacts_mem[1];
1775       interaction_destroy (iacts_mem[0]);
1776     }
1777
1778
1779   if (examine.id_var)
1780     {
1781       examine.id_idx = var_get_case_index (examine.id_var);
1782       examine.id_width = var_get_width (examine.id_var);
1783     }
1784
1785   examine.ex_proto = caseproto_add_width (examine.ex_proto, 0); /* value */
1786   examine.ex_proto = caseproto_add_width (examine.ex_proto, examine.id_width);   /* id */
1787   examine.ex_proto = caseproto_add_width (examine.ex_proto, 0); /* weight */
1788
1789
1790   if (examine.disp_extremes > 0)
1791     {
1792       examine.calc_extremes = examine.disp_extremes;
1793     }
1794
1795   if (examine.descriptives && examine.calc_extremes == 0)
1796     {
1797       /* Descriptives always displays the max and min */
1798       examine.calc_extremes = 1;
1799     }
1800
1801   if (percentiles_seen && examine.n_percentiles == 0)
1802     {
1803       examine.n_percentiles = 7;
1804       examine.ptiles = xcalloc (examine.n_percentiles, sizeof (*examine.ptiles));
1805
1806       examine.ptiles[0] = 5;
1807       examine.ptiles[1] = 10;
1808       examine.ptiles[2] = 25;
1809       examine.ptiles[3] = 50;
1810       examine.ptiles[4] = 75;
1811       examine.ptiles[5] = 90;
1812       examine.ptiles[6] = 95;
1813     }
1814
1815   assert (examine.calc_extremes >= examine.disp_extremes);
1816   {
1817     struct casegrouper *grouper;
1818     struct casereader *group;
1819     bool ok;
1820
1821     grouper = casegrouper_create_splits (proc_open (ds), examine.dict);
1822     while (casegrouper_get_next_group (grouper, &group))
1823       run_examine (&examine, group);
1824     ok = casegrouper_destroy (grouper);
1825     ok = proc_commit (ds) && ok;
1826   }
1827
1828   caseproto_unref (examine.ex_proto);
1829
1830   for (i = 0; i < examine.n_iacts; ++i)
1831     interaction_destroy (examine.iacts[i]);
1832   free (examine.ptiles);
1833   free (examine.dep_vars);
1834   pool_destroy (examine.pool);
1835
1836   return CMD_SUCCESS;
1837
1838  error:
1839   caseproto_unref (examine.ex_proto);
1840   examine.iacts = iacts_mem;
1841   for (i = 0; i < examine.n_iacts; ++i)
1842     interaction_destroy (examine.iacts[i]);
1843   free (examine.dep_vars);
1844   free (examine.ptiles);
1845   pool_destroy (examine.pool);
1846
1847   return CMD_FAILURE;
1848 }