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