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