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