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