HISTOGRAMS: Fix bin width problems on large numbers of bins
[pspp] / src / language / stats / examine.c
1 /*
2   PSPP - a program for statistical analysis.
3   Copyright (C) 2012 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/tab.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 histogram;
133   bool boxplot;
134   bool spreadlevel;
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
200 /* Returns an array of (iact->n_vars) pointers to union value initialised to NULL.
201    The caller must free this array when no longer required. */
202 static const union value **
203 previous_value_alloc (const struct interaction *iact)
204 {
205   int ivar_idx;
206
207   const union value **prev_val = xcalloc (iact->n_vars, sizeof (*prev_val));
208
209   for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
210     prev_val[ivar_idx] = NULL;
211
212   return prev_val;
213 }
214
215 /* Set the contents of PREV_VAL to the values of C indexed by the variables of IACT */
216 static int
217 previous_value_record (const struct interaction *iact, const struct ccase *c, const union value **prev_val)
218 {
219   int ivar_idx;
220   int diff_idx = -1;
221
222   for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
223     {
224       const struct variable *ivar = iact->vars[ivar_idx];
225       const int width = var_get_width (ivar);
226       const union value *val = case_data (c, ivar);
227                   
228       if (prev_val[ivar_idx])
229         if (! value_equal (prev_val[ivar_idx], val, width))
230           {
231             diff_idx = ivar_idx;
232             break;
233           }
234     }
235
236   for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
237     {
238       const struct variable *ivar = iact->vars[ivar_idx];
239       const union value *val = case_data (c, ivar);
240       
241       prev_val[ivar_idx] = val;
242     }
243   return diff_idx;
244 }
245
246
247 static void
248 show_boxplot_grouped (const struct examine *cmd, int iact_idx)
249 {
250   int v;
251
252   const struct interaction *iact = cmd->iacts[iact_idx];
253   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
254
255   for (v = 0; v < cmd->n_dep_vars; ++v)
256     {
257       double y_min = DBL_MAX;
258       double y_max = -DBL_MAX;
259       int grp;
260       struct boxplot *boxplot;
261       struct string title;
262       ds_init_empty (&title);
263
264       if (iact->n_vars > 0)
265         {
266           struct string istr;
267           ds_init_empty (&istr);
268           interaction_to_string (iact, &istr);
269           ds_put_format (&title, _("Boxplot of %s vs. %s"),
270                          var_to_string (cmd->dep_vars[v]),
271                          ds_cstr (&istr));
272           ds_destroy (&istr);
273         }
274       else
275         ds_put_format (&title, _("Boxplot of %s"), var_to_string (cmd->dep_vars[v]));
276       
277       for (grp = 0; grp < n_cats; ++grp)
278         {
279           const struct exploratory_stats *es =
280             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
281
282           if ( y_min > es[v].minimum)
283             y_min = es[v].minimum;
284
285           if ( y_max < es[v].maximum)
286             y_max = es[v].maximum;
287         }
288       
289       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
290
291       ds_destroy (&title);
292
293       for (grp = 0; grp < n_cats; ++grp)
294         {
295           int ivar_idx;
296           struct string label;
297
298           const struct ccase *c =
299             categoricals_get_case_by_category_real (cmd->cats,  iact_idx, grp);
300
301           const struct exploratory_stats *es =
302             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
303
304           ds_init_empty (&label);
305           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
306             {
307               const struct variable *ivar = iact->vars[ivar_idx];
308               const union value *val = case_data (c, ivar);
309               
310               ds_put_cstr (&label, var_to_string (ivar));
311               ds_put_cstr (&label, " = ");
312               append_value_name (ivar, val, &label);
313               ds_put_cstr (&label, "; ");
314             }
315
316           boxplot_add_box (boxplot, es[v].box_whisker, ds_cstr (&label));
317
318           ds_destroy (&label);
319         }
320       
321       boxplot_submit (boxplot);
322     }
323 }
324
325 static void
326 show_boxplot_variabled (const struct examine *cmd, int iact_idx)
327 {
328   int grp;
329   const struct interaction *iact = cmd->iacts[iact_idx];
330   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
331
332   for (grp = 0; grp < n_cats; ++grp)
333     {
334       struct boxplot *boxplot;
335       int v;
336       double y_min = DBL_MAX;
337       double y_max = -DBL_MAX;
338
339       const struct ccase *c =
340         categoricals_get_case_by_category_real (cmd->cats,  iact_idx, grp);
341
342       struct string title;
343       ds_init_empty (&title);
344
345       for (v = 0; v < cmd->n_dep_vars; ++v)
346         {
347           const struct exploratory_stats *es =
348             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
349
350           if ( y_min > es[v].minimum)
351             y_min = es[v].minimum;
352
353           if ( y_max < es[v].maximum)
354             y_max = es[v].maximum;
355         }
356
357       if ( iact->n_vars == 0)
358         ds_put_format (&title, _("Boxplot"));
359       else
360         {
361           int ivar_idx;
362           struct string label;
363           ds_init_empty (&label);
364           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
365             {
366               const struct variable *ivar = iact->vars[ivar_idx];
367               const union value *val = case_data (c, ivar);
368               
369               ds_put_cstr (&label, var_to_string (ivar));
370               ds_put_cstr (&label, " = ");
371               append_value_name (ivar, val, &label);
372               ds_put_cstr (&label, "; ");
373             }
374
375           ds_put_format (&title, _("Boxplot of %s"),
376                          ds_cstr (&label));
377
378           ds_destroy (&label);
379         }
380
381       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
382
383       ds_destroy (&title);
384
385       for (v = 0; v < cmd->n_dep_vars; ++v)
386         {
387           const struct exploratory_stats *es =
388             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
389
390           boxplot_add_box (boxplot, es[v].box_whisker, 
391                            var_to_string (cmd->dep_vars[v]));
392         }
393
394       boxplot_submit (boxplot);
395     }
396 }
397
398
399 static void
400 show_npplot (const struct examine *cmd, int iact_idx)
401 {
402   const struct interaction *iact = cmd->iacts[iact_idx];
403   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
404
405   int v;
406
407   for (v = 0; v < cmd->n_dep_vars; ++v)
408     {
409       int grp;
410       for (grp = 0; grp < n_cats; ++grp)
411         {
412           struct chart_item *npp, *dnpp;
413           struct casereader *reader;
414           struct np *np;
415
416           int ivar_idx;
417           const struct ccase *c =
418             categoricals_get_case_by_category_real (cmd->cats,
419                                                     iact_idx, grp);
420
421           const struct exploratory_stats *es =
422             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
423
424           struct string label;
425           ds_init_cstr (&label, 
426                         var_to_string (cmd->dep_vars[v]));
427
428           if ( iact->n_vars > 0)
429             {
430               ds_put_cstr (&label, " (");
431               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
432                 {
433                   const struct variable *ivar = iact->vars[ivar_idx];
434                   const union value *val = case_data (c, ivar);
435                   
436                   ds_put_cstr (&label, var_to_string (ivar));
437                   ds_put_cstr (&label, " = ");
438                   append_value_name (ivar, val, &label);
439                   ds_put_cstr (&label, "; ");
440                   
441                 }
442               ds_put_cstr (&label, ")");
443             }
444           
445           np = es[v].np;
446           reader = casewriter_make_reader (np->writer);
447           np->writer = NULL;
448
449           npp = np_plot_create (np, reader, ds_cstr (&label));
450           dnpp = dnp_plot_create (np, reader, ds_cstr (&label));
451
452           if (npp == NULL || dnpp == NULL)
453             {
454               msg (MW, _("Not creating NP plot because data set is empty."));
455               chart_item_unref (npp);
456               chart_item_unref (dnpp);
457             }
458           else
459             {
460               chart_item_submit (npp);
461               chart_item_submit (dnpp);
462             }
463           casereader_destroy (reader);
464
465           ds_destroy (&label);
466         }
467     }
468 }
469
470 static void
471 show_spreadlevel (const struct examine *cmd, int iact_idx)
472 {
473   const struct interaction *iact = cmd->iacts[iact_idx];
474   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
475
476   int v;
477
478   /* Spreadlevel when there are no levels is not useful */
479   if (iact->n_vars == 0)
480     return;
481
482   for (v = 0; v < cmd->n_dep_vars; ++v)
483     {
484       int grp;
485       struct chart_item *sl;
486
487       struct string label;
488       ds_init_cstr (&label, 
489                     var_to_string (cmd->dep_vars[v]));
490
491       if (iact->n_vars > 0)
492         {
493           ds_put_cstr (&label, " (");
494           interaction_to_string (iact, &label);
495           ds_put_cstr (&label, ")");
496         }
497       
498       sl = spreadlevel_plot_create (ds_cstr (&label), cmd->sl_power);
499
500       for (grp = 0; grp < n_cats; ++grp)
501         {
502           const struct exploratory_stats *es =
503             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
504
505           double median = percentile_calculate (es[v].quartiles[1], cmd->pc_alg);
506
507           double iqr = percentile_calculate (es[v].quartiles[2], cmd->pc_alg) -
508             percentile_calculate (es[v].quartiles[0], cmd->pc_alg);
509
510           spreadlevel_plot_add (sl, iqr, median);
511         }
512
513       if (sl == NULL)
514         msg (MW, _("Not creating spreadlevel chart for %s"), ds_cstr (&label));
515       else 
516         chart_item_submit (sl);
517
518       ds_destroy (&label);
519     }
520 }
521
522
523 static void
524 show_histogram (const struct examine *cmd, int iact_idx)
525 {
526   const struct interaction *iact = cmd->iacts[iact_idx];
527   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
528
529   int v;
530
531   for (v = 0; v < cmd->n_dep_vars; ++v)
532     {
533       int grp;
534       for (grp = 0; grp < n_cats; ++grp)
535         {
536           double n, mean, var;
537           int ivar_idx;
538           const struct ccase *c =
539             categoricals_get_case_by_category_real (cmd->cats,
540                                                     iact_idx, grp);
541
542           const struct exploratory_stats *es =
543             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
544
545           struct string label;
546
547           if (es[v].histogram == NULL)
548             continue;
549
550           ds_init_cstr (&label, 
551                         var_to_string (cmd->dep_vars[v]));
552
553           if ( iact->n_vars > 0)
554             {
555               ds_put_cstr (&label, " (");
556               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
557                 {
558                   const struct variable *ivar = iact->vars[ivar_idx];
559                   const union value *val = case_data (c, ivar);
560                   
561                   ds_put_cstr (&label, var_to_string (ivar));
562                   ds_put_cstr (&label, " = ");
563                   append_value_name (ivar, val, &label);
564                   ds_put_cstr (&label, "; ");
565                   
566                 }
567               ds_put_cstr (&label, ")");
568             }
569
570
571           moments_calculate (es[v].mom, &n, &mean, &var, NULL, NULL);
572
573           chart_item_submit
574             ( histogram_chart_create (es[v].histogram->gsl_hist,
575                                       ds_cstr (&label), n, mean,
576                                       sqrt (var), false));
577
578           
579           ds_destroy (&label);
580         }
581     }
582 }
583
584 static void
585 percentiles_report (const struct examine *cmd, int iact_idx)
586 {
587   const struct interaction *iact = cmd->iacts[iact_idx];
588   int i, v;
589   const int heading_columns = 1 + iact->n_vars + 1;
590   const int heading_rows = 2;
591   struct tab_table *t;
592
593   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
594
595   const int rows_per_cat = 2;
596   const int rows_per_var = n_cats * rows_per_cat;
597
598   const int nr = heading_rows + cmd->n_dep_vars * rows_per_var;
599   const int nc = heading_columns + cmd->n_percentiles;
600
601   t = tab_create (nc, nr);
602   tab_title (t, _("Percentiles"));
603
604   tab_headers (t, heading_columns, 0, heading_rows, 0);
605
606   /* Internal Vertical lines */
607   tab_box (t, -1, -1, -1, TAL_1,
608            heading_columns, 0, nc - 1, nr - 1);
609
610   /* External Frame */
611   tab_box (t, TAL_2, TAL_2, -1, -1,
612            0, 0, nc - 1, nr - 1);
613
614   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
615   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
616
617   tab_joint_text (t, heading_columns, 0,
618                   nc - 1, 0,
619                   TAT_TITLE | TAB_CENTER,
620                   _("Percentiles")
621                   );
622
623   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
624
625
626   for (i = 0; i < cmd->n_percentiles; ++i)
627     {
628       tab_text_format (t, heading_columns + i, 1,
629                        TAT_TITLE | TAB_CENTER,
630                        _("%g"), cmd->ptiles[i]);
631     }
632
633   for (i = 0; i < iact->n_vars; ++i)
634     {
635       tab_text (t,
636                 1 + i, 1,
637                 TAT_TITLE,
638                 var_to_string (iact->vars[i])
639                 );
640     }
641
642
643
644   if (n_cats > 0)
645     {
646       tab_vline (t, TAL_1, heading_columns - 1, heading_rows, nr - 1);
647
648       for (v = 0; v < cmd->n_dep_vars; ++v)
649         {
650           const union value **prev_vals = previous_value_alloc (iact);
651
652           int ivar_idx;
653           if ( v > 0 )
654             tab_hline (t, TAL_1, 0, nc - 1, heading_rows + v * rows_per_var);
655         
656           tab_text (t,
657                     0, heading_rows + v * rows_per_var,
658                     TAT_TITLE | TAB_LEFT,
659                     var_to_string (cmd->dep_vars[v])
660                     );
661
662           for (i = 0; i < n_cats; ++i)
663             {
664               const struct ccase *c =
665                 categoricals_get_case_by_category_real (cmd->cats,
666                                                         iact_idx, i);
667
668               const struct exploratory_stats *ess =
669                 categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, i);
670
671               const struct exploratory_stats *es = ess + v;
672
673               int diff_idx = previous_value_record (iact, c, prev_vals);
674
675               double hinges[3];
676               int p;
677
678               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
679                 {
680                   const struct variable *ivar = iact->vars[ivar_idx];
681                   const union value *val = case_data (c, ivar);
682
683                   if (( diff_idx != -1 && diff_idx <= ivar_idx)
684                       || i == 0)
685                     {              
686                       struct string str;
687                       ds_init_empty (&str);
688                       append_value_name (ivar, val, &str);
689               
690                       tab_text (t,
691                                 1 + ivar_idx,
692                                 heading_rows + v * rows_per_var + i * rows_per_cat,
693                                 TAT_TITLE | TAB_LEFT,
694                                 ds_cstr (&str)
695                                 );
696                   
697                       ds_destroy (&str);
698                     }
699                 }
700
701               if ( diff_idx != -1 && diff_idx < iact->n_vars)
702                 {
703                   tab_hline (t, TAL_1, 1 + diff_idx, nc - 1,
704                              heading_rows + v * rows_per_var + i * rows_per_cat
705                              );
706                 }
707
708               tab_text (t, heading_columns - 1, 
709                         heading_rows + v * rows_per_var + i * rows_per_cat,
710                         TAT_TITLE | TAB_LEFT,
711                         gettext (ptile_alg_desc [cmd->pc_alg]));
712
713               tukey_hinges_calculate (es->hinges, hinges);
714
715               for (p = 0; p < cmd->n_percentiles; ++p)
716                 {
717                   tab_double (t, heading_columns + p, 
718                               heading_rows + v * rows_per_var + i * rows_per_cat,
719                               0,
720                               percentile_calculate (es->percentiles[p], cmd->pc_alg),
721                               0);
722               
723                   if (cmd->ptiles[p] == 25.0)
724                     {
725                       tab_double (t, heading_columns + p, 
726                                   heading_rows + v * rows_per_var + i * rows_per_cat + 1,
727                                   0,
728                                   hinges[0],
729                                   0);
730                     }
731                   else if (cmd->ptiles[p] == 50.0)
732                     {
733                       tab_double (t, heading_columns + p, 
734                                   heading_rows + v * rows_per_var + i * rows_per_cat + 1,
735                                   0,
736                                   hinges[1],
737                                   0);
738                     }
739                   else if (cmd->ptiles[p] == 75.0)
740                     {
741                       tab_double (t, heading_columns + p, 
742                                   heading_rows + v * rows_per_var + i * rows_per_cat + 1,
743                                   0,
744                                   hinges[2],
745                                   0);
746                     }
747                 }
748
749
750               tab_text (t, heading_columns - 1, 
751                         heading_rows + v * rows_per_var + i * rows_per_cat + 1,
752                         TAT_TITLE | TAB_LEFT,
753                         _("Tukey's Hinges"));
754           
755             }
756
757           free (prev_vals);
758         }
759     }
760   tab_submit (t);
761 }
762
763 static void
764 descriptives_report (const struct examine *cmd, int iact_idx)
765 {
766   const struct interaction *iact = cmd->iacts[iact_idx];
767   int i, v;
768   const int heading_columns = 1 + iact->n_vars + 2;
769   const int heading_rows = 1;
770   struct tab_table *t;
771
772   size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
773
774   const int rows_per_cat = 13;
775   const int rows_per_var = n_cats * rows_per_cat;
776
777   const int nr = heading_rows + cmd->n_dep_vars * rows_per_var;
778   const int nc = 2 + heading_columns;
779
780   t = tab_create (nc, nr);
781   tab_title (t, _("Descriptives"));
782
783   tab_headers (t, heading_columns, 0, heading_rows, 0);
784
785   /* Internal Vertical lines */
786   tab_box (t, -1, -1, -1, TAL_1,
787            heading_columns, 0, nc - 1, nr - 1);
788
789   /* External Frame */
790   tab_box (t, TAL_2, TAL_2, -1, -1,
791            0, 0, nc - 1, nr - 1);
792
793   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
794   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
795
796
797   tab_text (t, heading_columns, 0, TAB_CENTER | TAT_TITLE,
798             _("Statistic"));
799
800   tab_text (t, heading_columns + 1, 0, TAB_CENTER | TAT_TITLE,
801             _("Std. Error"));
802
803   for (i = 0; i < iact->n_vars; ++i)
804     {
805       tab_text (t,
806                 1 + i, 0,
807                 TAT_TITLE,
808                 var_to_string (iact->vars[i])
809                 );
810     }
811
812   for (v = 0; v < cmd->n_dep_vars; ++v)
813     {
814       const union value **prev_val = previous_value_alloc (iact);
815
816       int ivar_idx;
817       if ( v > 0 )
818         tab_hline (t, TAL_1, 0, nc - 1, heading_rows + v * rows_per_var);
819         
820       tab_text (t,
821                 0, heading_rows + v * rows_per_var,
822                 TAT_TITLE | TAB_LEFT,
823                 var_to_string (cmd->dep_vars[v])
824                 );
825
826       for (i = 0; i < n_cats; ++i)
827         {
828           const struct ccase *c =
829             categoricals_get_case_by_category_real (cmd->cats,
830                                                     iact_idx, i);
831
832           const struct exploratory_stats *ess =
833             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, i);
834
835           const struct exploratory_stats *es = ess + v;
836
837           const int diff_idx = previous_value_record (iact, c, prev_val);
838
839           double m0, m1, m2, m3, m4;
840           double tval;
841
842           moments_calculate (es->mom, &m0, &m1, &m2, &m3, &m4);
843
844           tval = gsl_cdf_tdist_Qinv ((1.0 - cmd->conf) / 2.0, m0 - 1.0);
845
846           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
847             {
848               const struct variable *ivar = iact->vars[ivar_idx];
849               const union value *val = case_data (c, ivar);
850
851               if (( diff_idx != -1 && diff_idx <= ivar_idx)
852                   || i == 0)
853                 {              
854                   struct string str;
855                   ds_init_empty (&str);
856                   append_value_name (ivar, val, &str);
857               
858                   tab_text (t,
859                             1 + ivar_idx,
860                             heading_rows + v * rows_per_var + i * rows_per_cat,
861                             TAT_TITLE | TAB_LEFT,
862                             ds_cstr (&str)
863                             );
864                   
865                   ds_destroy (&str);
866                 }
867             }
868
869           if ( diff_idx != -1 && diff_idx < iact->n_vars)
870             {
871               tab_hline (t, TAL_1, 1 + diff_idx, nc - 1,
872                          heading_rows + v * rows_per_var + i * rows_per_cat
873                          );
874             }
875
876           tab_text (t,
877                     1 + iact->n_vars,
878                     heading_rows + v * rows_per_var + i * rows_per_cat,
879                     TAB_LEFT,
880                     _("Mean")
881                     );
882
883           tab_double (t,
884                       1 + iact->n_vars + 2,
885                       heading_rows + v * rows_per_var + i * rows_per_cat,
886                       0, m1, 0);
887
888           tab_double (t,
889                       1 + iact->n_vars + 3,
890                       heading_rows + v * rows_per_var + i * rows_per_cat,
891                       0, calc_semean (m2, m0), 0);
892
893           tab_text_format (t,
894                            1 + iact->n_vars,
895                            heading_rows + v * rows_per_var + i * rows_per_cat + 1,
896                            TAB_LEFT,
897                            _("%g%% Confidence Interval for Mean"),
898                            cmd->conf * 100.0
899                            );
900           
901           tab_text (t,
902                     1 + iact->n_vars + 1,
903                     heading_rows + v * rows_per_var + i * rows_per_cat + 1,
904                     TAB_LEFT,
905                     _("Lower Bound")
906                     );
907
908           tab_double (t,
909                       1 + iact->n_vars + 2,
910                       heading_rows + v * rows_per_var + i * rows_per_cat + 1,
911                       0, m1 - tval * calc_semean (m2, m0), 0);
912
913
914           tab_text (t,
915                     1 + iact->n_vars + 1,
916                     heading_rows + v * rows_per_var + i * rows_per_cat + 2,
917                     TAB_LEFT,
918                     _("Upper Bound")
919                     );
920
921           tab_double (t,
922                       1 + iact->n_vars + 2,
923                       heading_rows + v * rows_per_var + i * rows_per_cat + 2,
924                       0, m1 + tval * calc_semean (m2, m0), 0);
925
926
927           tab_text (t,
928                     1 + iact->n_vars,
929                     heading_rows + v * rows_per_var + i * rows_per_cat + 3,
930                     TAB_LEFT,
931                     _("5% Trimmed Mean")
932                     );
933
934           tab_double (t,
935                       1 + iact->n_vars + 2,
936                       heading_rows + v * rows_per_var + i * rows_per_cat + 3,
937                       0,
938                       trimmed_mean_calculate (es->trimmed_mean),
939                       0);
940
941           tab_text (t,
942                     1 + iact->n_vars,
943                     heading_rows + v * rows_per_var + i * rows_per_cat + 4,
944                     TAB_LEFT,
945                     _("Median")
946                     );
947           
948           tab_double (t,
949                       1 + iact->n_vars + 2,
950                       heading_rows + v * rows_per_var + i * rows_per_cat + 4,
951                       0,
952                       percentile_calculate (es->quartiles[1], cmd->pc_alg),
953                       0);
954
955
956           tab_text (t,
957                     1 + iact->n_vars,
958                     heading_rows + v * rows_per_var + i * rows_per_cat + 5,
959                     TAB_LEFT,
960                     _("Variance")
961                     );
962
963           tab_double (t,
964                       1 + iact->n_vars + 2,
965                       heading_rows + v * rows_per_var + i * rows_per_cat + 5,
966                       0, m2, 0);
967
968           tab_text (t,
969                     1 + iact->n_vars,
970                     heading_rows + v * rows_per_var + i * rows_per_cat + 6,
971                     TAB_LEFT,
972                     _("Std. Deviation")
973                     );
974
975           tab_double (t,
976                       1 + iact->n_vars + 2,
977                       heading_rows + v * rows_per_var + i * rows_per_cat + 6,
978                       0, sqrt (m2), 0);
979
980           tab_text (t,
981                     1 + iact->n_vars,
982                     heading_rows + v * rows_per_var + i * rows_per_cat + 7,
983                     TAB_LEFT,
984                     _("Minimum")
985                     );
986
987           tab_double (t,
988                       1 + iact->n_vars + 2,
989                       heading_rows + v * rows_per_var + i * rows_per_cat + 7,
990                       0, 
991                       es->minima[0].val,
992                       0);
993
994           tab_text (t,
995                     1 + iact->n_vars,
996                     heading_rows + v * rows_per_var + i * rows_per_cat + 8,
997                     TAB_LEFT,
998                     _("Maximum")
999                     );
1000
1001           tab_double (t,
1002                       1 + iact->n_vars + 2,
1003                       heading_rows + v * rows_per_var + i * rows_per_cat + 8,
1004                       0, 
1005                       es->maxima[0].val,
1006                       0);
1007
1008           tab_text (t,
1009                     1 + iact->n_vars,
1010                     heading_rows + v * rows_per_var + i * rows_per_cat + 9,
1011                     TAB_LEFT,
1012                     _("Range")
1013                     );
1014
1015           tab_double (t,
1016                       1 + iact->n_vars + 2,
1017                       heading_rows + v * rows_per_var + i * rows_per_cat + 9,
1018                       0, 
1019                       es->maxima[0].val - es->minima[0].val,
1020                       0);
1021
1022           tab_text (t,
1023                     1 + iact->n_vars,
1024                     heading_rows + v * rows_per_var + i * rows_per_cat + 10,
1025                     TAB_LEFT,
1026                     _("Interquartile Range")
1027                     );
1028
1029
1030           tab_double (t,
1031                       1 + iact->n_vars + 2,
1032                       heading_rows + v * rows_per_var + i * rows_per_cat + 10,
1033                       0,
1034                       percentile_calculate (es->quartiles[2], cmd->pc_alg) - 
1035                       percentile_calculate (es->quartiles[0], cmd->pc_alg),
1036                       0);
1037
1038
1039
1040
1041           tab_text (t,
1042                     1 + iact->n_vars,
1043                     heading_rows + v * rows_per_var + i * rows_per_cat + 11,
1044                     TAB_LEFT,
1045                     _("Skewness")
1046                     );
1047
1048           tab_double (t,
1049                       1 + iact->n_vars + 2,
1050                       heading_rows + v * rows_per_var + i * rows_per_cat + 11,
1051                       0, m3, 0);
1052
1053           tab_double (t,
1054                       1 + iact->n_vars + 3,
1055                       heading_rows + v * rows_per_var + i * rows_per_cat + 11,
1056                       0, calc_seskew (m0), 0);
1057
1058           tab_text (t,
1059                     1 + iact->n_vars,
1060                     heading_rows + v * rows_per_var + i * rows_per_cat + 12,
1061                     TAB_LEFT,
1062                     _("Kurtosis")
1063                     );
1064
1065           tab_double (t,
1066                       1 + iact->n_vars + 2,
1067                       heading_rows + v * rows_per_var + i * rows_per_cat + 12,
1068                       0, m4, 0);
1069
1070           tab_double (t,
1071                       1 + iact->n_vars + 3,
1072                       heading_rows + v * rows_per_var + i * rows_per_cat + 12,
1073                       0, calc_sekurt (m0), 0);
1074         }
1075
1076       free (prev_val);
1077     }
1078   tab_submit (t);
1079 }
1080
1081
1082 static void
1083 extremes_report (const struct examine *cmd, int iact_idx)
1084 {
1085   const struct interaction *iact = cmd->iacts[iact_idx];
1086   int i, v;
1087   const int heading_columns = 1 + iact->n_vars + 2;
1088   const int heading_rows = 1;
1089   struct tab_table *t;
1090
1091   size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
1092
1093   const int rows_per_cat = 2 * cmd->disp_extremes;
1094   const int rows_per_var = n_cats * rows_per_cat;
1095
1096   const int nr = heading_rows + cmd->n_dep_vars * rows_per_var;
1097   const int nc = 2 + heading_columns;
1098
1099   t = tab_create (nc, nr);
1100   tab_title (t, _("Extreme Values"));
1101
1102   tab_headers (t, heading_columns, 0, heading_rows, 0);
1103
1104   /* Internal Vertical lines */
1105   tab_box (t, -1, -1, -1, TAL_1,
1106            heading_columns, 0, nc - 1, nr - 1);
1107
1108   /* External Frame */
1109   tab_box (t, TAL_2, TAL_2, -1, -1,
1110            0, 0, nc - 1, nr - 1);
1111
1112   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1113   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1114
1115
1116   if ( cmd->id_var ) 
1117     tab_text (t, heading_columns, 0, TAB_CENTER | TAT_TITLE,
1118               var_to_string (cmd->id_var));
1119   else
1120     tab_text (t, heading_columns, 0, TAB_CENTER | TAT_TITLE,
1121               _("Case Number"));
1122
1123   tab_text (t, heading_columns + 1, 0, TAB_CENTER | TAT_TITLE,
1124             _("Value"));
1125
1126   for (i = 0; i < iact->n_vars; ++i)
1127     {
1128       tab_text (t,
1129                 1 + i, 0,
1130                 TAT_TITLE,
1131                 var_to_string (iact->vars[i])
1132                 );
1133     }
1134
1135   for (v = 0; v < cmd->n_dep_vars; ++v)
1136     {
1137       const union value **prev_val = previous_value_alloc (iact);
1138
1139       int ivar_idx;
1140       if ( v > 0 )
1141         tab_hline (t, TAL_1, 0, nc - 1, heading_rows + v * rows_per_var);
1142         
1143       tab_text (t,
1144                 0, heading_rows + v * rows_per_var,
1145                 TAT_TITLE,
1146                 var_to_string (cmd->dep_vars[v])
1147                 );
1148
1149       for (i = 0; i < n_cats; ++i)
1150         {
1151           int e;
1152           const struct ccase *c =
1153             categoricals_get_case_by_category_real (cmd->cats, iact_idx, i);
1154
1155           const struct exploratory_stats *ess =
1156             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, i);
1157
1158           const struct exploratory_stats *es = ess + v;
1159
1160           int diff_idx = previous_value_record (iact, c, prev_val);
1161
1162           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
1163             {
1164               const struct variable *ivar = iact->vars[ivar_idx];
1165               const union value *val = case_data (c, ivar);
1166
1167               if (( diff_idx != -1 && diff_idx <= ivar_idx)
1168                   || i == 0)
1169                 {              
1170                   struct string str;
1171                   ds_init_empty (&str);
1172                   append_value_name (ivar, val, &str);
1173               
1174                   tab_text (t,
1175                             1 + ivar_idx,
1176                             heading_rows + v * rows_per_var + i * rows_per_cat,
1177                             TAT_TITLE | TAB_LEFT,
1178                             ds_cstr (&str)
1179                             );
1180                   
1181                   ds_destroy (&str);
1182                 }
1183             }
1184
1185           if ( diff_idx != -1 && diff_idx < iact->n_vars)
1186             {
1187               tab_hline (t, TAL_1, 1 + diff_idx, nc - 1,
1188                          heading_rows + v * rows_per_var + i * rows_per_cat
1189                          );
1190             }
1191           
1192           tab_text (t,
1193                     heading_columns - 2,
1194                     heading_rows + v * rows_per_var + i * rows_per_cat,
1195                     TAB_RIGHT,
1196                     _("Highest"));
1197
1198
1199           tab_hline (t, TAL_1, heading_columns - 2, nc - 1,
1200                      heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes
1201                      );
1202
1203           tab_text (t,
1204                     heading_columns - 2,
1205                     heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes,
1206                     TAB_RIGHT,
1207                     _("Lowest"));
1208
1209           for (e = 0 ; e < cmd->disp_extremes; ++e)
1210             {
1211               tab_double (t,
1212                           heading_columns - 1,
1213                           heading_rows + v * rows_per_var + i * rows_per_cat + e,
1214                           TAB_RIGHT,
1215                           e + 1,
1216                           &F_8_0);
1217
1218               /* The casenumber */
1219               if (cmd->id_var)
1220                 tab_value (t,
1221                            heading_columns,
1222                            heading_rows + v * rows_per_var + i * rows_per_cat + e,
1223                            TAB_RIGHT,
1224                            &es->maxima[e].identity,
1225                            cmd->id_var,
1226                            NULL);
1227               else 
1228                 tab_double (t,
1229                           heading_columns,
1230                             heading_rows + v * rows_per_var + i * rows_per_cat + e,
1231                             TAB_RIGHT,
1232                             es->maxima[e].identity.f,
1233                             &F_8_0);
1234
1235               tab_double (t,
1236                          heading_columns + 1,
1237                          heading_rows + v * rows_per_var + i * rows_per_cat + e,
1238                          0,
1239                          es->maxima[e].val,
1240                          var_get_print_format (cmd->dep_vars[v]));
1241                          
1242
1243               tab_double (t,
1244                           heading_columns - 1,
1245                           heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes + e,
1246                           TAB_RIGHT,
1247                           e + 1,
1248                           &F_8_0);
1249
1250               /* The casenumber */
1251               if (cmd->id_var)
1252                 tab_value (t,
1253                            heading_columns,
1254                            heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes + e,
1255                            TAB_RIGHT,
1256                            &es->minima[e].identity,
1257                            cmd->id_var,
1258                            NULL);
1259               else
1260                 tab_double (t,
1261                             heading_columns,
1262                             heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes + e,
1263                             TAB_RIGHT,
1264                             es->minima[e].identity.f,
1265                             &F_8_0);
1266
1267               tab_double (t,
1268                           heading_columns + 1,
1269                           heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes + e,
1270                           0,
1271                           es->minima[e].val,
1272                           var_get_print_format (cmd->dep_vars[v]));
1273             }
1274         }
1275       free (prev_val);
1276     }
1277
1278   tab_submit (t);
1279 }
1280
1281
1282 static void
1283 summary_report (const struct examine *cmd, int iact_idx)
1284 {
1285   const struct interaction *iact = cmd->iacts[iact_idx];
1286   int i, v;
1287   const int heading_columns = 1 + iact->n_vars;
1288   const int heading_rows = 3;
1289   struct tab_table *t;
1290
1291   const struct fmt_spec *wfmt = cmd->wv ? var_get_print_format (cmd->wv) : &F_8_0;
1292
1293   size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
1294
1295   const int nr = heading_rows + n_cats * cmd->n_dep_vars;
1296   const int nc = 6 + heading_columns;
1297
1298   t = tab_create (nc, nr);
1299   tab_title (t, _("Case Processing Summary"));
1300
1301   tab_headers (t, heading_columns, 0, heading_rows, 0);
1302
1303   /* Internal Vertical lines */
1304   tab_box (t, -1, -1, -1, TAL_1,
1305            heading_columns, 0, nc - 1, nr - 1);
1306
1307   /* External Frame */
1308   tab_box (t, TAL_2, TAL_2, -1, -1,
1309            0, 0, nc - 1, nr - 1);
1310
1311   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1312   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1313
1314   tab_joint_text (t, heading_columns, 0,
1315                   nc - 1, 0, TAB_CENTER | TAT_TITLE, _("Cases"));
1316   tab_joint_text (t,
1317                   heading_columns, 1,
1318                   heading_columns + 1, 1,
1319                   TAB_CENTER | TAT_TITLE, _("Valid"));
1320
1321   tab_joint_text (t,
1322                   heading_columns + 2, 1, 
1323                   heading_columns + 3, 1,
1324                   TAB_CENTER | TAT_TITLE, _("Missing"));
1325
1326   tab_joint_text (t,
1327                   heading_columns + 4, 1,
1328                   heading_columns + 5, 1,
1329                   TAB_CENTER | TAT_TITLE, _("Total"));
1330
1331   for (i = 0; i < 3; ++i)
1332     {
1333       tab_text (t, heading_columns + i * 2, 2, TAB_CENTER | TAT_TITLE,
1334                 _("N"));
1335       tab_text (t, heading_columns + i * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1336                 _("Percent"));
1337     }
1338
1339   for (i = 0; i < iact->n_vars; ++i)
1340     {
1341       tab_text (t,
1342                 1 + i, 2,
1343                 TAT_TITLE,
1344                 var_to_string (iact->vars[i])
1345                 );
1346     }
1347
1348   if (n_cats > 0)
1349     for (v = 0; v < cmd->n_dep_vars; ++v)
1350       {
1351         int ivar_idx;
1352         const union value **prev_values = previous_value_alloc (iact);
1353
1354         if ( v > 0 )
1355           tab_hline (t, TAL_1, 0, nc - 1, heading_rows + v * n_cats);
1356
1357         tab_text (t,
1358                   0, heading_rows + n_cats * v,
1359                   TAT_TITLE,
1360                   var_to_string (cmd->dep_vars[v])
1361                   );
1362
1363
1364         for (i = 0; i < n_cats; ++i)
1365           {
1366             double total;
1367             const struct exploratory_stats *es;
1368
1369             const struct ccase *c =
1370               categoricals_get_case_by_category_real (cmd->cats,
1371                                                       iact_idx, i);
1372             if (c)
1373               {
1374                 int diff_idx = previous_value_record (iact, c, prev_values);
1375
1376                 if ( diff_idx != -1 && diff_idx < iact->n_vars - 1)
1377                   tab_hline (t, TAL_1, 1 + diff_idx, nc - 1,
1378                              heading_rows + n_cats * v + i );
1379
1380                 for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
1381                   {
1382                     const struct variable *ivar = iact->vars[ivar_idx];
1383                     const union value *val = case_data (c, ivar);
1384
1385                     if (( diff_idx != -1 && diff_idx <= ivar_idx)
1386                         || i == 0)
1387                       {              
1388                         struct string str;
1389                         ds_init_empty (&str);
1390                         append_value_name (ivar, val, &str);
1391               
1392                         tab_text (t,
1393                                   1 + ivar_idx, heading_rows + n_cats * v + i,
1394                                   TAT_TITLE | TAB_LEFT,
1395                                   ds_cstr (&str)
1396                                   );
1397                   
1398                         ds_destroy (&str);
1399                       }
1400                   }
1401               }
1402
1403
1404             es = categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, i);
1405   
1406           
1407             total = es[v].missing + es[v].non_missing;
1408             tab_double (t, 
1409                         heading_columns + 0,
1410                         heading_rows + n_cats * v + i,
1411                         0,
1412                         es[v].non_missing,
1413                         wfmt);
1414
1415
1416             tab_text_format (t, 
1417                              heading_columns + 1,
1418                              heading_rows + n_cats * v + i,
1419                              0,
1420                              "%g%%",
1421                              100.0 * es[v].non_missing / total
1422                              );
1423
1424
1425             tab_double (t, 
1426                         heading_columns + 2,
1427                         heading_rows + n_cats * v + i,
1428                         0,
1429                         es[v].missing,
1430                         wfmt);
1431
1432             tab_text_format (t, 
1433                              heading_columns + 3,
1434                              heading_rows + n_cats * v + i,
1435                              0,
1436                              "%g%%",
1437                              100.0 * es[v].missing / total
1438                              );
1439             tab_double (t, 
1440                         heading_columns + 4,
1441                         heading_rows + n_cats * v + i,
1442                         0,
1443                         total,
1444                         wfmt);
1445
1446             /* This can only be 100% can't it? */
1447             tab_text_format (t, 
1448                              heading_columns + 5,
1449                              heading_rows + n_cats * v + i,
1450                              0,
1451                              "%g%%",
1452                              100.0 * (es[v].missing + es[v].non_missing)/ total
1453                              );
1454           }
1455         free (prev_values);
1456       }
1457
1458   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1459   tab_hline (t, TAL_1, heading_columns, nc - 1, 2);
1460
1461   tab_submit (t);
1462 }
1463
1464
1465 /* Match a variable.
1466    If the match succeeds, the variable will be placed in VAR.
1467    Returns true if successful */
1468 static bool
1469 lex_match_variable (struct lexer *lexer, 
1470                     const struct dictionary *dict, const struct variable **var)
1471 {
1472   if (lex_token (lexer) !=  T_ID)
1473
1474     return false;
1475
1476   *var = parse_variable_const  (lexer, dict);
1477
1478   if ( *var == NULL)
1479     return false;
1480   return true;
1481 }
1482
1483 /* Attempt to parse an interaction from LEXER */
1484 static struct interaction *
1485 parse_interaction (struct lexer *lexer, struct examine *ex)
1486 {
1487   const struct variable *v = NULL;
1488   struct interaction *iact = NULL;
1489   
1490   if ( lex_match_variable (lexer, ex->dict, &v))
1491     {
1492       iact = interaction_create (v);
1493
1494       while (lex_match (lexer, T_BY))
1495         {
1496           if (!lex_match_variable (lexer, ex->dict, &v))
1497             {
1498               interaction_destroy (iact);
1499               return NULL;
1500             }
1501           interaction_add_variable (iact, v);
1502         }
1503       lex_match (lexer, T_COMMA);
1504     }
1505   
1506   return iact;
1507 }
1508
1509
1510 static void *
1511 create_n (const void *aux1, void *aux2 UNUSED)
1512 {
1513   int v;
1514   
1515   const struct examine *examine = aux1;
1516   struct exploratory_stats *es = pool_calloc (examine->pool, examine->n_dep_vars, sizeof (*es));
1517   struct subcase ordering;
1518   subcase_init (&ordering, 0, 0, SC_ASCEND);
1519
1520   for (v = 0; v < examine->n_dep_vars; v++)
1521     {
1522       es[v].sorted_writer = sort_create_writer (&ordering, examine->ex_proto);
1523       es[v].sorted_reader = NULL;
1524
1525       es[v].mom = moments_create (MOMENT_KURTOSIS);
1526       es[v].cmin = DBL_MAX;
1527
1528       es[v].maximum = -DBL_MAX;
1529       es[v].minimum =  DBL_MAX;
1530     }
1531
1532   subcase_destroy (&ordering);
1533   return es;
1534 }
1535
1536 static void
1537 update_n (const void *aux1, void *aux2 UNUSED, void *user_data,
1538           const struct ccase *c, double weight)
1539 {
1540   int v;
1541   const struct examine *examine = aux1;
1542   struct exploratory_stats *es = user_data;
1543
1544   for (v = 0; v < examine->n_dep_vars; v++)
1545     {
1546       struct ccase *outcase ;
1547       const struct variable *var = examine->dep_vars[v];
1548       const double x = case_data (c, var)->f;
1549       
1550       if (var_is_value_missing (var, case_data (c, var), examine->dep_excl))
1551         {
1552           es[v].missing += weight;
1553           continue;
1554         }
1555
1556       outcase = case_create (examine->ex_proto);
1557
1558       if (x > es[v].maximum)
1559         es[v].maximum = x;
1560
1561       if (x < es[v].minimum)
1562         es[v].minimum =  x;
1563
1564       es[v].non_missing += weight;
1565
1566       moments_pass_one (es[v].mom, x, weight);
1567
1568       /* Save the value and the ID to the writer */
1569       assert (examine->id_idx != -1);
1570       case_data_rw_idx (outcase, EX_VAL)->f = x;
1571       value_copy (case_data_rw_idx (outcase, EX_ID),
1572                   case_data_idx (c, examine->id_idx), examine->id_width);
1573
1574       case_data_rw_idx (outcase, EX_WT)->f = weight;
1575       
1576       es[v].cc += weight;
1577
1578       if (es[v].cmin > weight)
1579         es[v].cmin = weight;
1580
1581       casewriter_write (es[v].sorted_writer, outcase);
1582     }
1583 }
1584
1585 static void
1586 calculate_n (const void *aux1, void *aux2 UNUSED, void *user_data)
1587 {
1588   int v;
1589   const struct examine *examine = aux1;
1590   struct exploratory_stats *es = user_data;
1591
1592   for (v = 0; v < examine->n_dep_vars; v++)
1593     {
1594       int i;
1595       casenumber imin = 0;
1596       double imax = es[v].cc;
1597       struct casereader *reader;
1598       struct ccase *c;
1599       casenumber total_cases;
1600
1601       if (examine->histogram)
1602         {
1603           /* Sturges Rule */
1604           double bin_width = fabs (es[v].minimum - es[v].maximum)
1605             / (1 + log2 (es[v].cc))
1606             ;
1607
1608           es[v].histogram =
1609             histogram_create (bin_width, es[v].minimum, es[v].maximum);
1610         }
1611
1612       es[v].sorted_reader = casewriter_make_reader (es[v].sorted_writer);
1613       total_cases = casereader_count_cases (es[v].sorted_reader);
1614       es[v].sorted_writer = NULL;
1615
1616       es[v].maxima = pool_calloc (examine->pool, examine->calc_extremes, sizeof (*es[v].maxima));
1617       es[v].minima = pool_calloc (examine->pool, examine->calc_extremes, sizeof (*es[v].minima));
1618       for (i = 0; i < examine->calc_extremes; ++i)
1619         {
1620           value_init_pool (examine->pool, &es[v].maxima[i].identity, examine->id_width) ;
1621           value_init_pool (examine->pool, &es[v].minima[i].identity, examine->id_width) ;
1622         }
1623       
1624       for (reader = casereader_clone (es[v].sorted_reader);
1625            (c = casereader_read (reader)) != NULL; case_unref (c))
1626         {
1627           const double val = case_data_idx (c, EX_VAL)->f;
1628           const double wt = case_data_idx (c, EX_WT)->f;   /* FIXME: What about fractional weights ??? */
1629
1630           moments_pass_two (es[v].mom, val, wt);
1631
1632           if (es[v].histogram)
1633             histogram_add (es[v].histogram, val, wt);
1634
1635           if (imin < examine->calc_extremes)
1636             {
1637               int x;
1638               for (x = imin; x < examine->calc_extremes; ++x)
1639                 {
1640                   struct extremity *min = &es[v].minima[x];
1641                   min->val = val;
1642                   value_copy (&min->identity, case_data_idx (c, EX_ID), examine->id_width);
1643                 }
1644               imin += wt;
1645             }
1646
1647           imax -= wt;
1648           if (imax < examine->calc_extremes)
1649             {
1650               int x;
1651
1652               for (x = imax; x < imax + wt; ++x)
1653                 {
1654                   struct extremity *max;
1655
1656                   if (x >= examine->calc_extremes) 
1657                     break;
1658
1659                   max = &es[v].maxima[x];
1660                   max->val = val;
1661                   value_copy (&max->identity, case_data_idx (c, EX_ID), examine->id_width);
1662                 }
1663             }
1664         }
1665       casereader_destroy (reader);
1666
1667       if (examine->calc_extremes > 0)
1668         {
1669           assert (es[v].minima[0].val == es[v].minimum);
1670           assert (es[v].maxima[0].val == es[v].maximum);
1671         }
1672
1673       {
1674         const int n_os = 5 + examine->n_percentiles;
1675         struct order_stats **os ;
1676         es[v].percentiles = pool_calloc (examine->pool, examine->n_percentiles, sizeof (*es[v].percentiles));
1677
1678         es[v].trimmed_mean = trimmed_mean_create (es[v].cc, 0.05);
1679
1680         os = xcalloc (n_os, sizeof *os);
1681         os[0] = &es[v].trimmed_mean->parent;
1682
1683         es[v].quartiles[0] = percentile_create (0.25, es[v].cc);
1684         es[v].quartiles[1] = percentile_create (0.5,  es[v].cc);
1685         es[v].quartiles[2] = percentile_create (0.75, es[v].cc);
1686
1687         os[1] = &es[v].quartiles[0]->parent;
1688         os[2] = &es[v].quartiles[1]->parent;
1689         os[3] = &es[v].quartiles[2]->parent;
1690
1691         es[v].hinges = tukey_hinges_create (es[v].cc, es[v].cmin);
1692         os[4] = &es[v].hinges->parent;
1693
1694         for (i = 0; i < examine->n_percentiles; ++i)
1695           {
1696             es[v].percentiles[i] = percentile_create (examine->ptiles[i] / 100.00, es[v].cc);
1697             os[5 + i] = &es[v].percentiles[i]->parent;
1698           }
1699
1700         order_stats_accumulate_idx (os, n_os,
1701                                     casereader_clone (es[v].sorted_reader),
1702                                     EX_WT, EX_VAL);
1703
1704         free (os);
1705       }
1706
1707       if (examine->boxplot)
1708         {
1709           struct order_stats *os;
1710
1711           es[v].box_whisker = box_whisker_create (es[v].hinges, 
1712                                                   EX_ID, examine->id_var);
1713
1714           os = &es[v].box_whisker->parent;
1715           order_stats_accumulate_idx (&os, 1,
1716                                       casereader_clone (es[v].sorted_reader),
1717                                       EX_WT, EX_VAL);
1718         }
1719
1720       if (examine->npplot)
1721         {
1722           double n, mean, var;
1723           struct order_stats *os;
1724
1725           moments_calculate (es[v].mom, &n, &mean, &var, NULL, NULL);
1726           
1727           es[v].np = np_create (n, mean, var);
1728
1729           os = &es[v].np->parent;
1730
1731           order_stats_accumulate_idx (&os, 1,
1732                                       casereader_clone (es[v].sorted_reader),
1733                                       EX_WT, EX_VAL);
1734         }
1735
1736     }
1737 }
1738
1739 static void
1740 cleanup_exploratory_stats (struct examine *cmd)
1741
1742   int i;
1743   for (i = 0; i < cmd->n_iacts; ++i)
1744     {
1745       int v;
1746       const size_t n_cats =  categoricals_n_count (cmd->cats, i);
1747
1748       for (v = 0; v < cmd->n_dep_vars; ++v)
1749         {
1750           int grp;
1751           for (grp = 0; grp < n_cats; ++grp)
1752             {
1753               int q;
1754               const struct exploratory_stats *es =
1755                 categoricals_get_user_data_by_category_real (cmd->cats, i, grp);
1756
1757               struct order_stats *os = &es[v].hinges->parent;
1758               struct statistic  *stat = &os->parent;
1759               stat->destroy (stat);
1760
1761               for (q = 0; q < 3 ; q++)
1762                 {
1763                   os = &es[v].quartiles[q]->parent;
1764                   stat = &os->parent;
1765                   stat->destroy (stat);
1766                 }
1767
1768               for (q = 0; q < cmd->n_percentiles ; q++)
1769                 {
1770                   os = &es[v].percentiles[q]->parent;
1771                   stat = &os->parent;
1772                   stat->destroy (stat);
1773                 }
1774
1775               os = &es[v].trimmed_mean->parent;
1776               stat = &os->parent;
1777               stat->destroy (stat);
1778
1779               os = &es[v].np->parent;
1780               if (os)
1781                 {
1782                   stat = &os->parent;
1783                   stat->destroy (stat);
1784                 }
1785
1786               statistic_destroy (&es[v].histogram->parent);
1787               moments_destroy (es[v].mom);
1788
1789               casereader_destroy (es[v].sorted_reader);
1790             }
1791         }
1792     }
1793 }
1794
1795
1796 static void
1797 run_examine (struct examine *cmd, struct casereader *input)
1798 {
1799   int i;
1800   struct ccase *c;
1801   struct casereader *reader;
1802
1803   struct payload payload;
1804   payload.create = create_n;
1805   payload.update = update_n;
1806   payload.destroy = calculate_n;
1807   
1808   cmd->wv = dict_get_weight (cmd->dict);
1809
1810   cmd->cats
1811     = categoricals_create (cmd->iacts, cmd->n_iacts,  
1812                            cmd->wv, cmd->dep_excl, cmd->fctr_excl);
1813
1814   categoricals_set_payload (cmd->cats, &payload, cmd, NULL);
1815
1816   if (cmd->id_idx == -1)
1817     {
1818       struct ccase *c = casereader_peek (input,  0);
1819
1820       assert (cmd->id_var == NULL);
1821
1822       cmd->id_idx = case_get_value_cnt (c);
1823       input = casereader_create_arithmetic_sequence (input, 1.0, 1.0);
1824
1825       case_unref (c);
1826     }
1827
1828   /* Remove cases on a listwise basis if requested */
1829   if ( cmd->missing_pw == false)
1830     input = casereader_create_filter_missing (input,
1831                                               cmd->dep_vars,
1832                                               cmd->n_dep_vars,
1833                                               cmd->dep_excl,
1834                                               NULL,
1835                                               NULL);
1836
1837   for (reader = input;
1838        (c = casereader_read (reader)) != NULL; case_unref (c))
1839     {
1840       categoricals_update (cmd->cats, c);
1841     }
1842   casereader_destroy (reader);
1843   categoricals_done (cmd->cats);
1844
1845   for (i = 0; i < cmd->n_iacts; ++i)
1846     {
1847       summary_report (cmd, i);
1848
1849       if (cmd->disp_extremes > 0)
1850         extremes_report (cmd, i);
1851
1852       if (cmd->n_percentiles > 0)
1853         percentiles_report (cmd, i);
1854
1855       if (cmd->boxplot)
1856         {
1857           switch (cmd->boxplot_mode)
1858             {
1859             case BP_GROUPS:
1860               show_boxplot_grouped (cmd, i);
1861               break;
1862             case BP_VARIABLES:
1863               show_boxplot_variabled (cmd, i);
1864               break;
1865             default:
1866               NOT_REACHED ();
1867               break;
1868             }
1869         }
1870
1871       if (cmd->histogram)
1872         show_histogram (cmd, i);
1873
1874       if (cmd->npplot)
1875         show_npplot (cmd, i);
1876
1877       if (cmd->spreadlevel)
1878         show_spreadlevel (cmd, i);
1879
1880       if (cmd->descriptives)
1881         descriptives_report (cmd, i);
1882     }
1883
1884   cleanup_exploratory_stats (cmd);
1885   categoricals_destroy (cmd->cats);
1886 }
1887
1888
1889 int
1890 cmd_examine (struct lexer *lexer, struct dataset *ds)
1891 {
1892   int i;
1893   bool nototals_seen = false;
1894   bool totals_seen = false;
1895
1896   struct interaction **iacts_mem = NULL;
1897   struct examine examine;
1898   bool percentiles_seen = false;
1899
1900   examine.missing_pw = false;
1901   examine.disp_extremes = 0;
1902   examine.calc_extremes = 0;
1903   examine.descriptives = false;
1904   examine.conf = 0.95;
1905   examine.pc_alg = PC_HAVERAGE;
1906   examine.ptiles = NULL;
1907   examine.n_percentiles = 0;
1908   examine.id_idx = -1;
1909   examine.id_width = 0;
1910   examine.id_var = NULL;
1911   examine.boxplot_mode = BP_GROUPS;
1912   
1913   examine.ex_proto = caseproto_create ();
1914
1915   examine.pool = pool_create ();
1916
1917   /* Allocate space for the first interaction.
1918      This is interaction is an empty one (for the totals).
1919      If no totals are requested, we will simply ignore this
1920      interaction.
1921   */
1922   examine.n_iacts = 1;
1923   examine.iacts = iacts_mem = pool_zalloc (examine.pool, sizeof (struct interaction *));
1924   examine.iacts[0] = interaction_create (NULL);
1925
1926   examine.dep_excl = MV_ANY;
1927   examine.fctr_excl = MV_ANY;
1928   examine.histogram = false;
1929   examine.npplot = false;
1930   examine.boxplot = false;
1931   examine.spreadlevel = false;
1932   examine.sl_power = 0;
1933   
1934   examine.dict = dataset_dict (ds);
1935
1936   /* Accept an optional, completely pointless "/VARIABLES=" */
1937   lex_match (lexer, T_SLASH);
1938   if (lex_match_id  (lexer, "VARIABLES"))
1939     {
1940       if (! lex_force_match (lexer, T_EQUALS) )
1941         goto error;
1942     }
1943
1944   if (!parse_variables_const (lexer, examine.dict,
1945                               &examine.dep_vars, &examine.n_dep_vars,
1946                               PV_NO_DUPLICATE | PV_NUMERIC))
1947     goto error;
1948
1949   if (lex_match (lexer, T_BY))
1950     {
1951       struct interaction *iact = NULL;
1952       do
1953         {
1954           iact = parse_interaction (lexer, &examine);
1955           if (iact)
1956             {
1957               examine.n_iacts++;
1958               iacts_mem = 
1959                 pool_nrealloc (examine.pool, iacts_mem,
1960                                examine.n_iacts,
1961                                sizeof (*iacts_mem));
1962               
1963               iacts_mem[examine.n_iacts - 1] = iact;
1964             }
1965         }
1966       while (iact);
1967     }
1968
1969
1970   while (lex_token (lexer) != T_ENDCMD)
1971     {
1972       lex_match (lexer, T_SLASH);
1973
1974       if (lex_match_id (lexer, "STATISTICS"))
1975         {
1976           lex_match (lexer, T_EQUALS);
1977
1978           while (lex_token (lexer) != T_ENDCMD
1979                  && lex_token (lexer) != T_SLASH)
1980             {
1981               if (lex_match_id (lexer, "DESCRIPTIVES"))
1982                 {
1983                   examine.descriptives = true;
1984                 }
1985               else if (lex_match_id (lexer, "EXTREME"))
1986                 {
1987                   int extr = 5;
1988                   if (lex_match (lexer, T_LPAREN))
1989                     {
1990                       extr = lex_integer (lexer);
1991
1992                       if (extr < 0)
1993                         {
1994                           msg (MW, _("%s may not be negative. Using default value (%g)."), "EXTREME", 5.0);
1995                           extr = 5;
1996                         }
1997
1998                       lex_get (lexer);
1999                       if (! lex_force_match (lexer, T_RPAREN))
2000                         goto error;
2001                     }
2002                   examine.disp_extremes  = extr;
2003                 }
2004               else if (lex_match_id (lexer, "NONE"))
2005                 {
2006                 }
2007               else if (lex_match (lexer, T_ALL))
2008                 {
2009                   if (examine.disp_extremes == 0)
2010                     examine.disp_extremes = 5;
2011                 }
2012               else
2013                 {
2014                   lex_error (lexer, NULL);
2015                   goto error;
2016                 }
2017             }
2018         }
2019       else if (lex_match_id (lexer, "PERCENTILES"))
2020         {
2021           percentiles_seen = true;
2022           if (lex_match (lexer, T_LPAREN))
2023             {
2024               while (lex_is_number (lexer))
2025                 {
2026                   double p = lex_number (lexer);
2027                   
2028                   if ( p <= 0 || p >= 100.0)
2029                     {
2030                       lex_error (lexer,
2031                                  _("Percentiles must lie in the range (0, 100)"));
2032                       goto error;
2033                     }
2034
2035                   examine.n_percentiles++;
2036                   examine.ptiles =
2037                     xrealloc (examine.ptiles,
2038                               sizeof (*examine.ptiles) *
2039                               examine.n_percentiles);
2040
2041                   examine.ptiles[examine.n_percentiles - 1] = p;
2042
2043                   lex_get (lexer);
2044                   lex_match (lexer, T_COMMA);
2045                 }
2046               if (!lex_force_match (lexer, T_RPAREN))
2047                 goto error;
2048             }
2049
2050           lex_match (lexer, T_EQUALS);
2051
2052           while (lex_token (lexer) != T_ENDCMD
2053                  && lex_token (lexer) != T_SLASH)
2054             {
2055               if (lex_match_id (lexer, "HAVERAGE"))
2056                 {
2057                   examine.pc_alg = PC_HAVERAGE;
2058                 }
2059               else if (lex_match_id (lexer, "WAVERAGE"))
2060                 {
2061                   examine.pc_alg = PC_WAVERAGE;
2062                 }
2063               else if (lex_match_id (lexer, "ROUND"))
2064                 {
2065                   examine.pc_alg = PC_ROUND;
2066                 }
2067               else if (lex_match_id (lexer, "EMPIRICAL"))
2068                 {
2069                   examine.pc_alg = PC_EMPIRICAL;
2070                 }
2071               else if (lex_match_id (lexer, "AEMPIRICAL"))
2072                 {
2073                   examine.pc_alg = PC_AEMPIRICAL;
2074                 }
2075               else if (lex_match_id (lexer, "NONE"))
2076                 {
2077                   examine.pc_alg = PC_NONE;
2078                 }
2079               else
2080                 {
2081                   lex_error (lexer, NULL);
2082                   goto error;
2083                 }
2084             }
2085         }
2086       else if (lex_match_id (lexer, "TOTAL"))
2087         {
2088           totals_seen = true;
2089         }
2090       else if (lex_match_id (lexer, "NOTOTAL"))
2091         {
2092           nototals_seen = true;
2093         }
2094       else if (lex_match_id (lexer, "MISSING"))
2095         {
2096           lex_match (lexer, T_EQUALS);
2097
2098           while (lex_token (lexer) != T_ENDCMD
2099                  && lex_token (lexer) != T_SLASH)
2100             {
2101               if (lex_match_id (lexer, "LISTWISE"))
2102                 {
2103                   examine.missing_pw = false;
2104                 }
2105               else if (lex_match_id (lexer, "PAIRWISE"))
2106                 {
2107                   examine.missing_pw = true;
2108                 }
2109               else if (lex_match_id (lexer, "EXCLUDE"))
2110                 {
2111                   examine.dep_excl = MV_ANY;
2112                 }
2113               else if (lex_match_id (lexer, "INCLUDE"))
2114                 {
2115                   examine.dep_excl = MV_SYSTEM;
2116                 }
2117               else if (lex_match_id (lexer, "REPORT"))
2118                 {
2119                   examine.fctr_excl = MV_NEVER;
2120                 }
2121               else if (lex_match_id (lexer, "NOREPORT"))
2122                 {
2123                   examine.fctr_excl = MV_ANY;
2124                 }
2125               else
2126                 {
2127                   lex_error (lexer, NULL);
2128                   goto error;
2129                 }
2130             }
2131         }
2132       else if (lex_match_id (lexer, "COMPARE"))
2133         {
2134           lex_match (lexer, T_EQUALS);
2135           if (lex_match_id (lexer, "VARIABLES"))
2136             {
2137               examine.boxplot_mode = BP_VARIABLES;
2138             }
2139           else if (lex_match_id (lexer, "GROUPS"))
2140             {
2141               examine.boxplot_mode = BP_GROUPS;
2142             }
2143           else
2144             {
2145               lex_error (lexer, NULL);
2146               goto error;
2147             }
2148         }
2149       else if (lex_match_id (lexer, "PLOT"))
2150         {
2151           lex_match (lexer, T_EQUALS);
2152
2153           while (lex_token (lexer) != T_ENDCMD
2154                  && lex_token (lexer) != T_SLASH)
2155             {
2156               if (lex_match_id (lexer, "BOXPLOT"))
2157                 {
2158                   examine.boxplot = true;
2159                 }
2160               else if (lex_match_id (lexer, "NPPLOT"))
2161                 {
2162                   examine.npplot = true;
2163                 }
2164               else if (lex_match_id (lexer, "HISTOGRAM"))
2165                 {
2166                   examine.histogram = true;
2167                 }
2168               else if (lex_match_id (lexer, "SPREADLEVEL"))
2169                 {
2170                   examine.spreadlevel = true;
2171                   examine.sl_power = 0;
2172                   if (lex_match (lexer, T_LPAREN))
2173                     {
2174                       examine.sl_power = lex_integer (lexer);
2175
2176                       lex_get (lexer);
2177                       if (! lex_force_match (lexer, T_RPAREN))
2178                         goto error;
2179                     }
2180                 }
2181               else if (lex_match_id (lexer, "NONE"))
2182                 {
2183                   examine.histogram = false;
2184                   examine.npplot = false;
2185                   examine.boxplot = false;
2186                 }
2187               else if (lex_match (lexer, T_ALL))
2188                 {
2189                   examine.histogram = true;
2190                   examine.npplot = true;
2191                   examine.boxplot = true;
2192                 }
2193               else 
2194                 {
2195                   lex_error (lexer, NULL);
2196                   goto error;
2197                 }
2198               lex_match (lexer, T_COMMA);
2199             }          
2200         }
2201       else if (lex_match_id (lexer, "CINTERVAL"))
2202         {
2203           if ( !lex_force_num (lexer))
2204             goto error;
2205         
2206           examine.conf = lex_number (lexer);
2207           lex_get (lexer);
2208         }
2209       else if (lex_match_id (lexer, "ID"))
2210         {
2211           lex_match (lexer, T_EQUALS);
2212
2213           examine.id_var = parse_variable_const (lexer, examine.dict);
2214         }
2215       else
2216         {
2217           lex_error (lexer, NULL);
2218           goto error;
2219         }
2220     }
2221
2222
2223   if ( totals_seen && nototals_seen)
2224     {
2225       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
2226       goto error;
2227     }
2228
2229   /* If totals have been requested or if there are no factors
2230      in this analysis, then the totals need to be included. */
2231   if ( !nototals_seen || examine.n_iacts == 1)
2232     {
2233       examine.iacts = &iacts_mem[0];
2234     }
2235   else
2236     {
2237       examine.n_iacts--;
2238       examine.iacts = &iacts_mem[1];
2239     }
2240
2241
2242   if ( examine.id_var )
2243     {
2244       examine.id_idx = var_get_case_index (examine.id_var);
2245       examine.id_width = var_get_width (examine.id_var);
2246     }
2247
2248   examine.ex_proto = caseproto_add_width (examine.ex_proto, 0); /* value */
2249   examine.ex_proto = caseproto_add_width (examine.ex_proto, examine.id_width);   /* id */
2250   examine.ex_proto = caseproto_add_width (examine.ex_proto, 0); /* weight */
2251
2252
2253   if (examine.disp_extremes > 0)
2254     {
2255       examine.calc_extremes = examine.disp_extremes;
2256     }
2257
2258   if (examine.descriptives && examine.calc_extremes == 0)
2259     {
2260       /* Descriptives always displays the max and min */
2261       examine.calc_extremes = 1;
2262     }
2263
2264   if (percentiles_seen && examine.n_percentiles == 0)
2265     {
2266       examine.n_percentiles = 7;
2267       examine.ptiles = xcalloc (examine.n_percentiles,
2268                                 sizeof (*examine.ptiles));
2269
2270       examine.ptiles[0] = 5;
2271       examine.ptiles[1] = 10;
2272       examine.ptiles[2] = 25;
2273       examine.ptiles[3] = 50;
2274       examine.ptiles[4] = 75;
2275       examine.ptiles[5] = 90;
2276       examine.ptiles[6] = 95;
2277     }
2278
2279   assert (examine.calc_extremes >= examine.disp_extremes);
2280   {
2281     struct casegrouper *grouper;
2282     struct casereader *group;
2283     bool ok;
2284     
2285     grouper = casegrouper_create_splits (proc_open (ds), examine.dict);
2286     while (casegrouper_get_next_group (grouper, &group))
2287       run_examine (&examine, group);
2288     ok = casegrouper_destroy (grouper);
2289     ok = proc_commit (ds) && ok;
2290   }
2291
2292   caseproto_unref (examine.ex_proto);
2293
2294   for (i = 0; i < examine.n_iacts; ++i)
2295     interaction_destroy (examine.iacts[i]);
2296
2297   free (examine.ptiles);
2298   free (examine.dep_vars);
2299   pool_destroy (examine.pool);
2300
2301   return CMD_SUCCESS;
2302
2303  error:
2304   caseproto_unref (examine.ex_proto);
2305   examine.iacts = iacts_mem;
2306   for (i = 0; i < examine.n_iacts; ++i)
2307     interaction_destroy (examine.iacts[i]);
2308   free (examine.dep_vars);
2309   free (examine.ptiles);
2310   pool_destroy (examine.pool);
2311
2312   return CMD_FAILURE;
2313 }