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