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