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