Fixed non-portable shell syntax in tests. Added a dist-hook target to check
[pspp-builds.git] / src / language / stats / examine.q
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2004 Free Software Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation, either version 3 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
16
17 #include <config.h>
18
19 #include <gsl/gsl_cdf.h>
20 #include <libpspp/message.h>
21 #include <math.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #include <data/case.h>
26 #include <data/casegrouper.h>
27 #include <data/casereader.h>
28 #include <data/dictionary.h>
29 #include <data/procedure.h>
30 #include <data/value-labels.h>
31 #include <data/variable.h>
32 #include <language/command.h>
33 #include <language/dictionary/split-file.h>
34 #include <language/lexer/lexer.h>
35 #include <libpspp/compiler.h>
36 #include <libpspp/hash.h>
37 #include <libpspp/message.h>
38 #include <libpspp/misc.h>
39 #include <libpspp/str.h>
40 #include <math/factor-stats.h>
41 #include <math/moments.h>
42 #include <math/percentiles.h>
43 #include <output/charts/box-whisker.h>
44 #include <output/charts/cartesian.h>
45 #include <output/manager.h>
46 #include <output/table.h>
47
48 #include "minmax.h"
49 #include "xalloc.h"
50
51 #include "gettext.h"
52 #define _(msgid) gettext (msgid)
53 #define N_(msgid) msgid
54
55 /* (headers) */
56 #include <output/chart.h>
57 #include <output/charts/plot-hist.h>
58 #include <output/charts/plot-chart.h>
59
60 /* (specification)
61    "EXAMINE" (xmn_):
62    *^variables=custom;
63    +total=custom;
64    +nototal=custom;
65    missing=miss:pairwise/!listwise,
66            rep:report/!noreport,
67            incl:include/!exclude;
68    +compare=cmp:variables/!groups;
69    +percentiles=custom;
70    +id=var;
71    +plot[plt_]=stemleaf,boxplot,npplot,:spreadlevel(*d:n),histogram,all,none;
72    +cinterval=double;
73    +statistics[st_]=descriptives,:extreme(*d:n),all,none.
74 */
75
76 /* (declarations) */
77
78 /* (functions) */
79
80
81
82 static struct cmd_examine cmd;
83
84 static const struct variable **dependent_vars;
85
86 static size_t n_dependent_vars;
87
88
89 struct factor
90 {
91   /* The independent variable */
92   struct variable *indep_var[2];
93
94
95   /* Hash table of factor stats indexed by 2 values */
96   struct hsh_table *fstats;
97
98   /* The hash table after it has been crunched */
99   struct factor_statistics **fs;
100
101   struct factor *next;
102
103 };
104
105 /* Linked list of factors */
106 static struct factor *factors = 0;
107
108 static struct metrics *totals = 0;
109
110 /* Parse the clause specifying the factors */
111 static int examine_parse_independent_vars (struct lexer *lexer, const struct dictionary *dict, struct cmd_examine *cmd);
112
113
114
115 /* Output functions */
116 static void show_summary (const struct variable **dependent_var, int n_dep_var,
117                          const struct factor *f);
118
119 static void show_extremes (const struct variable **dependent_var,
120                           int n_dep_var,
121                           const struct factor *factor,
122                           int n_extremities);
123
124 static void show_descriptives (const struct variable **dependent_var,
125                               int n_dep_var,
126                               struct factor *factor);
127
128 static void show_percentiles (const struct variable **dependent_var,
129                              int n_dep_var,
130                              struct factor *factor);
131
132
133
134
135 void np_plot (const struct metrics *m, const char *factorname);
136
137
138 void box_plot_group (const struct factor *fctr,
139                     const struct variable **vars, int n_vars,
140                     const struct variable *id
141                     ) ;
142
143
144 void box_plot_variables (const struct factor *fctr,
145                         const struct variable **vars, int n_vars,
146                         const struct variable *id
147                         );
148
149
150
151 /* Per Split function */
152 static void run_examine (struct cmd_examine *, struct casereader *,
153                          struct dataset *);
154
155 static void output_examine (void);
156
157
158 void factor_calc (const struct ccase *c, int case_no,
159                   double weight, bool case_missing);
160
161
162 /* Represent a factor as a string, so it can be
163    printed in a human readable fashion */
164 static void factor_to_string (const struct factor *fctr,
165                                const struct factor_statistics *fs,
166                               const struct variable *var,
167                               struct string *str
168                               );
169
170 /* Represent a factor as a string, so it can be
171    printed in a human readable fashion,
172    but sacrificing some readablility for the sake of brevity */
173 static void factor_to_string_concise (const struct factor *fctr,
174                                       const struct factor_statistics *fs,
175                                       struct string *);
176
177
178
179
180 /* Categories of missing values to exclude. */
181 static enum mv_class exclude_values;
182
183 /* PERCENTILES */
184
185 static subc_list_double percentile_list;
186
187 static enum pc_alg percentile_algorithm;
188
189 static short sbc_percentile;
190
191
192 int
193 cmd_examine (struct lexer *lexer, struct dataset *ds)
194 {
195   struct casegrouper *grouper;
196   struct casereader *group;
197   bool ok;
198
199   subc_list_double_create (&percentile_list);
200   percentile_algorithm = PC_HAVERAGE;
201
202   if ( !parse_examine (lexer, ds, &cmd, NULL) )
203     {
204       subc_list_double_destroy (&percentile_list);
205       return CMD_FAILURE;
206     }
207
208   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
209   exclude_values = cmd.incl == XMN_INCLUDE ? MV_SYSTEM : MV_ANY;
210
211   if ( cmd.st_n == SYSMIS )
212     cmd.st_n = 5;
213
214   if ( ! cmd.sbc_cinterval)
215     cmd.n_cinterval[0] = 95.0;
216
217   /* If descriptives have been requested, make sure the
218      quartiles are calculated */
219   if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
220     {
221       subc_list_double_push (&percentile_list, 25);
222       subc_list_double_push (&percentile_list, 50);
223       subc_list_double_push (&percentile_list, 75);
224     }
225
226   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
227   while (casegrouper_get_next_group (grouper, &group))
228     run_examine (&cmd, group, ds);
229   ok = casegrouper_destroy (grouper);
230   ok = proc_commit (ds) && ok;
231
232   if ( totals )
233     {
234       free ( totals );
235     }
236
237   if ( dependent_vars )
238     free (dependent_vars);
239
240   {
241     struct factor *f = factors ;
242     while ( f )
243       {
244         struct factor *ff = f;
245
246         f = f->next;
247         free ( ff->fs );
248         hsh_destroy ( ff->fstats ) ;
249         free ( ff ) ;
250       }
251     factors = 0;
252   }
253
254   subc_list_double_destroy (&percentile_list);
255
256   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
257 };
258
259
260
261 /* Show all the appropriate tables */
262 static void
263 output_examine (void)
264 {
265   struct factor *fctr;
266
267   /* Show totals if appropriate */
268   if ( ! cmd.sbc_nototal || factors == 0 )
269     {
270       show_summary (dependent_vars, n_dependent_vars, 0);
271
272       if ( cmd.sbc_statistics )
273         {
274           if ( cmd.a_statistics[XMN_ST_EXTREME])
275             show_extremes (dependent_vars, n_dependent_vars, 0, cmd.st_n);
276
277           if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES])
278             show_descriptives (dependent_vars, n_dependent_vars, 0);
279
280         }
281       if ( sbc_percentile )
282         show_percentiles (dependent_vars, n_dependent_vars, 0);
283
284       if ( cmd.sbc_plot)
285         {
286           int v;
287           if ( cmd.a_plot[XMN_PLT_STEMLEAF] )
288             msg (SW, _ ("%s is not currently supported."), "STEMLEAF");
289
290           if ( cmd.a_plot[XMN_PLT_SPREADLEVEL] )
291             msg (SW, _ ("%s is not currently supported."), "SPREADLEVEL");
292
293           if ( cmd.a_plot[XMN_PLT_NPPLOT] )
294             {
295               for ( v = 0 ; v < n_dependent_vars; ++v )
296                 np_plot (&totals[v], var_to_string (dependent_vars[v]));
297             }
298
299           if ( cmd.a_plot[XMN_PLT_BOXPLOT] )
300             {
301               if ( cmd.cmp == XMN_GROUPS )
302                 {
303                   box_plot_group (0, (const struct variable **) dependent_vars,
304                                   n_dependent_vars, cmd.v_id);
305                 }
306               else
307                 box_plot_variables (0,
308                                     (const struct variable **) dependent_vars,
309                                     n_dependent_vars, cmd.v_id);
310             }
311
312           if ( cmd.a_plot[XMN_PLT_HISTOGRAM] )
313             {
314               for ( v = 0 ; v < n_dependent_vars; ++v )
315                 {
316                   struct normal_curve normal;
317
318                   normal.N      = totals[v].n;
319                   normal.mean   = totals[v].mean;
320                   normal.stddev = totals[v].stddev;
321
322                   histogram_plot (totals[v].histogram,
323                                  var_to_string (dependent_vars[v]),
324                                  &normal, 0);
325                 }
326             }
327
328         }
329
330     }
331
332
333   /* Show grouped statistics  as appropriate */
334   fctr = factors;
335   while ( fctr )
336     {
337       show_summary (dependent_vars, n_dependent_vars, fctr);
338
339       if ( cmd.sbc_statistics )
340         {
341           if ( cmd.a_statistics[XMN_ST_EXTREME])
342             show_extremes (dependent_vars, n_dependent_vars, fctr, cmd.st_n);
343
344           if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES])
345             show_descriptives (dependent_vars, n_dependent_vars, fctr);
346         }
347
348       if ( sbc_percentile )
349         show_percentiles (dependent_vars, n_dependent_vars, fctr);
350
351
352       if ( cmd.sbc_plot)
353         {
354           size_t v;
355
356           struct factor_statistics **fs = fctr->fs ;
357
358           if ( cmd.a_plot[XMN_PLT_BOXPLOT] )
359             {
360               if ( cmd.cmp == XMN_VARIABLES )
361                 box_plot_variables (fctr,
362                                     (const struct variable **) dependent_vars,
363                                     n_dependent_vars, cmd.v_id);
364               else
365                 box_plot_group (fctr,
366                                 (const struct variable **) dependent_vars,
367                                 n_dependent_vars, cmd.v_id);
368             }
369
370           for ( v = 0 ; v < n_dependent_vars; ++v )
371             {
372
373               for ( fs = fctr->fs ; *fs ; ++fs )
374                 {
375                   struct string str;
376                   ds_init_empty (&str);
377                   factor_to_string (fctr, *fs, dependent_vars[v], &str);
378
379                   if ( cmd.a_plot[XMN_PLT_NPPLOT] )
380                     np_plot (& (*fs)->m[v], ds_cstr (&str));
381
382                   if ( cmd.a_plot[XMN_PLT_HISTOGRAM] )
383                     {
384                       struct normal_curve normal;
385
386                       normal.N      = (*fs)->m[v].n;
387                       normal.mean   = (*fs)->m[v].mean;
388                       normal.stddev = (*fs)->m[v].stddev;
389
390                       histogram_plot ((*fs)->m[v].histogram,
391                                      ds_cstr (&str) ,  &normal, 0);
392                     }
393
394                   ds_destroy (&str);
395
396                 } /* for ( fs .... */
397
398             } /* for ( v = 0 ..... */
399
400         }
401
402       fctr = fctr->next;
403     }
404
405 }
406
407
408 /* Create a hash table of percentiles and their values from the list of
409    percentiles */
410 static struct hsh_table *
411 list_to_ptile_hash (const subc_list_double *l)
412 {
413   int i;
414
415   struct hsh_table *h ;
416
417   h = hsh_create (subc_list_double_count (l),
418                  (hsh_compare_func *) ptile_compare,
419                  (hsh_hash_func *) ptile_hash,
420                  (hsh_free_func *) free,
421                  0);
422
423
424   for ( i = 0 ; i < subc_list_double_count (l) ; ++i )
425     {
426       struct percentile *p = xmalloc (sizeof *p);
427
428       p->p = subc_list_double_at (l,i);
429       p->v = SYSMIS;
430
431       hsh_insert (h, p);
432
433     }
434
435   return h;
436
437 }
438
439 /* Parse the PERCENTILES subcommand */
440 static int
441 xmn_custom_percentiles (struct lexer *lexer, struct dataset *ds UNUSED,
442                        struct cmd_examine *p UNUSED, void *aux UNUSED)
443 {
444   sbc_percentile = 1;
445
446   lex_match (lexer, '=');
447
448   lex_match (lexer, '(');
449
450   while ( lex_is_number (lexer) )
451     {
452       subc_list_double_push (&percentile_list, lex_number (lexer));
453
454       lex_get (lexer);
455
456       lex_match (lexer, ',') ;
457     }
458   lex_match (lexer, ')');
459
460   lex_match (lexer, '=');
461
462   if ( lex_match_id (lexer, "HAVERAGE"))
463     percentile_algorithm = PC_HAVERAGE;
464
465   else if ( lex_match_id (lexer, "WAVERAGE"))
466     percentile_algorithm = PC_WAVERAGE;
467
468   else if ( lex_match_id (lexer, "ROUND"))
469     percentile_algorithm = PC_ROUND;
470
471   else if ( lex_match_id (lexer, "EMPIRICAL"))
472     percentile_algorithm = PC_EMPIRICAL;
473
474   else if ( lex_match_id (lexer, "AEMPIRICAL"))
475     percentile_algorithm = PC_AEMPIRICAL;
476
477   else if ( lex_match_id (lexer, "NONE"))
478     percentile_algorithm = PC_NONE;
479
480
481   if ( 0 == subc_list_double_count (&percentile_list))
482     {
483       subc_list_double_push (&percentile_list, 5);
484       subc_list_double_push (&percentile_list, 10);
485       subc_list_double_push (&percentile_list, 25);
486       subc_list_double_push (&percentile_list, 50);
487       subc_list_double_push (&percentile_list, 75);
488       subc_list_double_push (&percentile_list, 90);
489       subc_list_double_push (&percentile_list, 95);
490     }
491
492   return 1;
493 }
494
495 /* TOTAL and NOTOTAL are simple, mutually exclusive flags */
496 static int
497 xmn_custom_total (struct lexer *lexer UNUSED, struct dataset *ds UNUSED, struct cmd_examine *p, void *aux UNUSED)
498 {
499   if ( p->sbc_nototal )
500     {
501       msg (SE, _ ("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
502       return 0;
503     }
504
505   return 1;
506 }
507
508 static int
509 xmn_custom_nototal (struct lexer *lexer UNUSED, struct dataset *ds UNUSED,
510                     struct cmd_examine *p, void *aux UNUSED)
511 {
512   if ( p->sbc_total )
513     {
514       msg (SE, _ ("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
515       return 0;
516     }
517
518   return 1;
519 }
520
521
522
523 /* Parser for the variables sub command
524    Returns 1 on success */
525 static int
526 xmn_custom_variables (struct lexer *lexer, struct dataset *ds, struct cmd_examine *cmd, void *aux UNUSED)
527 {
528   const struct dictionary *dict = dataset_dict (ds);
529   lex_match (lexer, '=');
530
531   if ( (lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
532       && lex_token (lexer) != T_ALL)
533     {
534       return 2;
535     }
536
537   if (!parse_variables_const (lexer, dict, &dependent_vars, &n_dependent_vars,
538                         PV_NO_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH) )
539     {
540       free (dependent_vars);
541       return 0;
542     }
543
544   assert (n_dependent_vars);
545
546   totals = xnmalloc (n_dependent_vars, sizeof *totals);
547
548   if ( lex_match (lexer, T_BY))
549     {
550       int success ;
551       success =  examine_parse_independent_vars (lexer, dict, cmd);
552       if ( success != 1 ) {
553         free (dependent_vars);
554         free (totals) ;
555       }
556       return success;
557     }
558
559   return 1;
560 }
561
562
563
564 /* Parse the clause specifying the factors */
565 static int
566 examine_parse_independent_vars (struct lexer *lexer, const struct dictionary *dict, struct cmd_examine *cmd)
567 {
568   int success;
569   struct factor *sf = xmalloc (sizeof *sf);
570
571   if ( (lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
572       && lex_token (lexer) != T_ALL)
573     {
574       free ( sf ) ;
575       return 2;
576     }
577
578
579   sf->indep_var[0] = parse_variable (lexer, dict);
580   sf->indep_var[1] = 0;
581
582   if ( lex_token (lexer) == T_BY )
583     {
584
585       lex_match (lexer, T_BY);
586
587       if ( (lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
588           && lex_token (lexer) != T_ALL)
589         {
590           free ( sf ) ;
591           return 2;
592         }
593
594       sf->indep_var[1] = parse_variable (lexer, dict);
595
596     }
597
598
599   sf->fstats = hsh_create (4,
600                           (hsh_compare_func *) factor_statistics_compare,
601                           (hsh_hash_func *) factor_statistics_hash,
602                           (hsh_free_func *) factor_statistics_free,
603                           0);
604
605   sf->next = factors;
606   factors = sf;
607
608   lex_match (lexer, ',');
609
610   if ( lex_token (lexer) == '.' || lex_token (lexer) == '/' )
611     return 1;
612
613   success =  examine_parse_independent_vars (lexer, dict, cmd);
614
615   if ( success != 1 )
616     free ( sf ) ;
617
618   return success;
619 }
620
621
622
623
624 void populate_percentiles (struct tab_table *tbl, int col, int row,
625                           const struct metrics *m);
626
627 void populate_descriptives (struct tab_table *t, int col, int row,
628                            const struct metrics *fs);
629
630 void populate_extremes (struct tab_table *t, int col, int row, int n,
631                        const struct metrics *m);
632
633 void populate_summary (struct tab_table *t, int col, int row,
634                       const struct metrics *m);
635
636
637
638
639 /* Perform calculations for the sub factors */
640 void
641 factor_calc (const struct ccase *c, int case_no, double weight,
642              bool case_missing)
643 {
644   size_t v;
645   struct factor *fctr = factors;
646
647   while ( fctr)
648     {
649       struct factor_statistics **foo ;
650       union value *indep_vals[2] ;
651
652       indep_vals[0] = value_dup (
653                                  case_data (c, fctr->indep_var[0]),
654                                  var_get_width (fctr->indep_var[0])
655                                  );
656
657       if ( fctr->indep_var[1] )
658         indep_vals[1] = value_dup (
659                                    case_data (c, fctr->indep_var[1]),
660                                    var_get_width (fctr->indep_var[1])
661                                    );
662       else
663         {
664           const union value sm = {SYSMIS};
665           indep_vals[1] = value_dup (&sm, 0);
666         }
667
668       assert (fctr->fstats);
669
670       foo = ( struct factor_statistics ** )
671         hsh_probe (fctr->fstats, (void *) &indep_vals);
672
673       if ( !*foo )
674         {
675
676           *foo = create_factor_statistics (n_dependent_vars,
677                                           indep_vals[0],
678                                           indep_vals[1]);
679
680           for ( v =  0 ; v  < n_dependent_vars ; ++v )
681             {
682               metrics_precalc ( & (*foo)->m[v] );
683             }
684
685         }
686       else
687         {
688           free (indep_vals[0]);
689           free (indep_vals[1]);
690         }
691
692       for ( v =  0 ; v  < n_dependent_vars ; ++v )
693         {
694           const struct variable *var = dependent_vars[v];
695           union value *val = value_dup (
696                                         case_data (c, var),
697                                         var_get_width (var)
698                                         );
699
700           if (case_missing || var_is_value_missing (var, val, exclude_values))
701             {
702               free (val);
703               val = NULL;
704             }
705
706           metrics_calc ( & (*foo)->m[v], val, weight, case_no);
707
708           free (val);
709         }
710
711       fctr = fctr->next;
712     }
713 }
714
715 static void
716 run_examine (struct cmd_examine *cmd, struct casereader *input,
717              struct dataset *ds)
718 {
719   struct dictionary *dict = dataset_dict (ds);
720   casenumber case_no;
721   struct ccase c;
722   int v;
723   bool ok;
724
725   struct factor *fctr;
726
727   if (!casereader_peek (input, 0, &c))
728     {
729       casereader_destroy (input);
730       return;
731     }
732   output_split_file_values (ds, &c);
733   case_destroy (&c);
734
735   input = casereader_create_filter_weight (input, dict, NULL, NULL);
736   input = casereader_create_counter (input, &case_no, 0);
737
738   /* Make sure we haven't got rubbish left over from a
739      previous split. */
740   fctr = factors;
741   while (fctr)
742     {
743       struct factor *next = fctr->next;
744
745       hsh_clear (fctr->fstats);
746
747       fctr->fs = 0;
748
749       fctr = next;
750     }
751
752   for ( v = 0 ; v < n_dependent_vars ; ++v )
753     metrics_precalc (&totals[v]);
754
755   for (; casereader_read (input, &c); case_destroy (&c))
756     {
757       bool case_missing = false;
758       const double weight = dict_get_case_weight (dict, &c, NULL);
759
760       if ( cmd->miss == XMN_LISTWISE )
761         {
762           for ( v = 0 ; v < n_dependent_vars ; ++v )
763             {
764               const struct variable *var = dependent_vars[v];
765               union value *val = value_dup (
766                                                   case_data (&c, var),
767                                                   var_get_width (var)
768                                                   );
769
770               if ( var_is_value_missing (var, val, exclude_values))
771                 case_missing = true;
772
773               free (val);
774             }
775         }
776
777       for ( v = 0 ; v < n_dependent_vars ; ++v )
778         {
779           const struct variable *var = dependent_vars[v];
780           union value *val = value_dup (
781                                         case_data (&c, var),
782                                         var_get_width (var)
783                                         );
784
785           if ( var_is_value_missing (var, val, exclude_values)
786                || case_missing )
787             {
788               free (val) ;
789               val = NULL;
790             }
791
792           metrics_calc (&totals[v], val, weight, case_no);
793
794           free (val);
795         }
796
797       factor_calc (&c, case_no, weight, case_missing);
798     }
799   ok = casereader_destroy (input);
800
801   for ( v = 0 ; v < n_dependent_vars ; ++v)
802     {
803       fctr = factors;
804       while ( fctr )
805         {
806           struct hsh_iterator hi;
807           struct factor_statistics *fs;
808
809           for ( fs = hsh_first (fctr->fstats, &hi);
810                 fs != 0 ;
811                 fs = hsh_next (fctr->fstats, &hi))
812             {
813
814               fs->m[v].ptile_hash = list_to_ptile_hash (&percentile_list);
815               fs->m[v].ptile_alg = percentile_algorithm;
816               metrics_postcalc (&fs->m[v]);
817             }
818
819           fctr = fctr->next;
820         }
821
822       totals[v].ptile_hash = list_to_ptile_hash (&percentile_list);
823       totals[v].ptile_alg = percentile_algorithm;
824       metrics_postcalc (&totals[v]);
825     }
826
827
828   /* Make sure that the combination of factors are complete */
829
830   fctr = factors;
831   while ( fctr )
832     {
833       struct hsh_iterator hi;
834       struct hsh_iterator hi0;
835       struct hsh_iterator hi1;
836       struct factor_statistics *fs;
837
838       struct hsh_table *idh0=0;
839       struct hsh_table *idh1=0;
840       union value *val0;
841       union value *val1;
842
843       idh0 = hsh_create (4, (hsh_compare_func *) compare_values,
844                          (hsh_hash_func *) hash_value,
845                         0,0);
846
847       idh1 = hsh_create (4, (hsh_compare_func *) compare_values,
848                          (hsh_hash_func *) hash_value,
849                         0,0);
850
851
852       for ( fs = hsh_first (fctr->fstats, &hi);
853             fs != 0 ;
854             fs = hsh_next (fctr->fstats, &hi))
855         {
856           hsh_insert (idh0, (void *) &fs->id[0]);
857           hsh_insert (idh1, (void *) &fs->id[1]);
858         }
859
860       /* Ensure that the factors combination is complete */
861       for ( val0 = hsh_first (idh0, &hi0);
862             val0 != 0 ;
863             val0 = hsh_next (idh0, &hi0))
864         {
865           for ( val1 = hsh_first (idh1, &hi1);
866                 val1 != 0 ;
867                 val1 = hsh_next (idh1, &hi1))
868             {
869               struct factor_statistics **ffs;
870               union value key[2];
871               key[0] = *val0;
872               key[1] = *val1;
873
874               ffs = (struct factor_statistics **)
875                 hsh_probe (fctr->fstats, (void *) &key );
876
877               if ( !*ffs ) {
878                 size_t i;
879                  (*ffs) = create_factor_statistics (n_dependent_vars,
880                                                    &key[0], &key[1]);
881                 for ( i = 0 ; i < n_dependent_vars ; ++i )
882                   metrics_precalc ( & (*ffs)->m[i]);
883               }
884             }
885         }
886
887       hsh_destroy (idh0);
888       hsh_destroy (idh1);
889
890       fctr->fs = (struct factor_statistics **) hsh_sort_copy (fctr->fstats);
891
892       fctr = fctr->next;
893     }
894
895   if (ok)
896     output_examine ();
897
898
899   if ( totals )
900     {
901       size_t i;
902       for ( i = 0 ; i < n_dependent_vars ; ++i )
903         {
904           metrics_destroy (&totals[i]);
905         }
906     }
907 }
908
909
910 static void
911 show_summary (const struct variable **dependent_var, int n_dep_var,
912              const struct factor *fctr)
913 {
914   static const char *subtitle[]=
915     {
916       N_ ("Valid"),
917       N_ ("Missing"),
918       N_ ("Total")
919     };
920
921   int i;
922   int heading_columns ;
923   int n_cols;
924   const int heading_rows = 3;
925   struct tab_table *tbl;
926
927   int n_rows ;
928   int n_factors = 1;
929
930   if ( fctr )
931     {
932       heading_columns = 2;
933       n_factors = hsh_count (fctr->fstats);
934       n_rows = n_dep_var * n_factors ;
935
936       if ( fctr->indep_var[1] )
937         heading_columns = 3;
938     }
939   else
940     {
941       heading_columns = 1;
942       n_rows = n_dep_var;
943     }
944
945   n_rows += heading_rows;
946
947   n_cols = heading_columns + 6;
948
949   tbl = tab_create (n_cols,n_rows,0);
950   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
951
952   tab_dim (tbl, tab_natural_dimensions);
953
954   /* Outline the box */
955   tab_box (tbl,
956            TAL_2, TAL_2,
957            -1, -1,
958            0, 0,
959            n_cols - 1, n_rows - 1);
960
961   /* Vertical lines for the data only */
962   tab_box (tbl,
963            -1, -1,
964            -1, TAL_1,
965            heading_columns, 0,
966            n_cols - 1, n_rows - 1);
967
968
969   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
970   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, 1 );
971   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, heading_rows -1 );
972
973   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
974
975
976   tab_title (tbl, _ ("Case Processing Summary"));
977
978   tab_joint_text (tbl, heading_columns, 0,
979                  n_cols -1, 0,
980                  TAB_CENTER | TAT_TITLE,
981                  _ ("Cases"));
982
983   /* Remove lines ... */
984   tab_box (tbl,
985            -1, -1,
986            TAL_0, TAL_0,
987            heading_columns, 0,
988            n_cols - 1, 0);
989
990   for ( i = 0 ; i < 3 ; ++i )
991     {
992       tab_text (tbl, heading_columns + i * 2 , 2, TAB_CENTER | TAT_TITLE,
993                 _ ("N"));
994
995       tab_text (tbl, heading_columns + i * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
996                 _ ("Percent"));
997
998       tab_joint_text (tbl, heading_columns + i*2 , 1,
999                      heading_columns + i * 2 + 1, 1,
1000                      TAB_CENTER | TAT_TITLE,
1001                      subtitle[i]);
1002
1003       tab_box (tbl, -1, -1,
1004                TAL_0, TAL_0,
1005                heading_columns + i * 2, 1,
1006                heading_columns + i * 2 + 1, 1);
1007     }
1008
1009
1010   /* Titles for the independent variables */
1011   if ( fctr )
1012     {
1013       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1014                 var_to_string (fctr->indep_var[0]));
1015
1016       if ( fctr->indep_var[1] )
1017         {
1018           tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1019                     var_to_string (fctr->indep_var[1]));
1020         }
1021     }
1022
1023
1024   for ( i = 0 ; i < n_dep_var ; ++i )
1025     {
1026       int n_factors = 1;
1027       if ( fctr )
1028         n_factors = hsh_count (fctr->fstats);
1029
1030       if ( i > 0 )
1031         tab_hline (tbl, TAL_1, 0, n_cols -1 , i * n_factors + heading_rows);
1032
1033       tab_text (tbl,
1034                 0, i * n_factors + heading_rows,
1035                 TAB_LEFT | TAT_TITLE,
1036                 var_to_string (dependent_var[i])
1037                 );
1038
1039       if ( !fctr )
1040         populate_summary (tbl, heading_columns,
1041                          (i * n_factors) + heading_rows,
1042                          &totals[i]);
1043       else
1044         {
1045           struct factor_statistics **fs = fctr->fs;
1046           int count = 0 ;
1047           const union value *prev = NULL;
1048
1049           while (*fs)
1050             {
1051               if ( !prev ||
1052                    0 != compare_values (prev, (*fs)->id[0],
1053                                    var_get_width (fctr->indep_var[0])))
1054                 {
1055                   struct string vstr;
1056                   ds_init_empty (&vstr);
1057                   var_append_value_name (fctr->indep_var[0],
1058                                       (*fs)->id[0], &vstr);
1059
1060                   tab_text (tbl,
1061                             1,
1062                             (i * n_factors ) + count +
1063                             heading_rows,
1064                             TAB_LEFT | TAT_TITLE,
1065                             ds_cstr (&vstr)
1066                             );
1067
1068                   ds_destroy (&vstr);
1069
1070                   if (fctr->indep_var[1] && count > 0 )
1071                     tab_hline (tbl, TAL_1, 1, n_cols - 1,
1072                               (i * n_factors ) + count + heading_rows);
1073                 }
1074
1075               prev = (*fs)->id[0];
1076
1077               if ( fctr->indep_var[1])
1078                 {
1079                   struct string vstr;
1080                   ds_init_empty (&vstr);
1081                   var_append_value_name (fctr->indep_var[1],
1082                                          (*fs)->id[1], &vstr);
1083                   tab_text (tbl,
1084                             2,
1085                             (i * n_factors ) + count +
1086                             heading_rows,
1087                             TAB_LEFT | TAT_TITLE,
1088                             ds_cstr (&vstr)
1089                             );
1090                   ds_destroy (&vstr);
1091                 }
1092
1093               populate_summary (tbl, heading_columns,
1094                                 (i * n_factors) + count
1095                                 + heading_rows,
1096                                 & (*fs)->m[i]);
1097
1098               count++ ;
1099               fs++;
1100             }
1101         }
1102     }
1103
1104   tab_submit (tbl);
1105 }
1106
1107
1108 void
1109 populate_summary (struct tab_table *t, int col, int row,
1110                  const struct metrics *m)
1111
1112 {
1113   const double total = m->n + m->n_missing ;
1114
1115   tab_float (t, col + 0, row + 0, TAB_RIGHT, m->n, 8, 0);
1116   tab_float (t, col + 2, row + 0, TAB_RIGHT, m->n_missing, 8, 0);
1117   tab_float (t, col + 4, row + 0, TAB_RIGHT, total, 8, 0);
1118
1119
1120   if ( total > 0 ) {
1121     tab_text (t, col + 1, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
1122               100.0 * m->n / total );
1123
1124     tab_text (t, col + 3, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
1125               100.0 * m->n_missing / total );
1126
1127     /* This seems a bit pointless !!! */
1128     tab_text (t, col + 5, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
1129               100.0 * total / total );
1130   }
1131 }
1132
1133
1134
1135 static void
1136 show_extremes (const struct variable **dependent_var, int n_dep_var,
1137               const struct factor *fctr, int n_extremities)
1138 {
1139   int i;
1140   int heading_columns ;
1141   int n_cols;
1142   const int heading_rows = 1;
1143   struct tab_table *tbl;
1144
1145   int n_factors = 1;
1146   int n_rows ;
1147
1148   if ( fctr )
1149     {
1150       heading_columns = 2;
1151       n_factors = hsh_count (fctr->fstats);
1152
1153       n_rows = n_dep_var * 2 * n_extremities * n_factors;
1154
1155       if ( fctr->indep_var[1] )
1156         heading_columns = 3;
1157     }
1158   else
1159     {
1160       heading_columns = 1;
1161       n_rows = n_dep_var * 2 * n_extremities;
1162     }
1163
1164   n_rows += heading_rows;
1165
1166   heading_columns += 2;
1167   n_cols = heading_columns + 2;
1168
1169   tbl = tab_create (n_cols,n_rows,0);
1170   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1171
1172   tab_dim (tbl, tab_natural_dimensions);
1173
1174   /* Outline the box, No internal lines*/
1175   tab_box (tbl,
1176            TAL_2, TAL_2,
1177            -1, -1,
1178            0, 0,
1179            n_cols - 1, n_rows - 1);
1180
1181   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1182
1183   tab_title (tbl, _ ("Extreme Values"));
1184
1185   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows -1);
1186   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows -1);
1187
1188   if ( fctr )
1189     {
1190       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1191                 var_to_string (fctr->indep_var[0]));
1192
1193       if ( fctr->indep_var[1] )
1194         tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1195                   var_to_string (fctr->indep_var[1]));
1196     }
1197
1198   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _ ("Value"));
1199   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _ ("Case Number"));
1200
1201   for ( i = 0 ; i < n_dep_var ; ++i )
1202     {
1203
1204       if ( i > 0 )
1205         tab_hline (tbl, TAL_1, 0, n_cols -1 ,
1206                   i * 2 * n_extremities * n_factors + heading_rows);
1207
1208       tab_text (tbl, 0,
1209                 i * 2 * n_extremities * n_factors  + heading_rows,
1210                 TAB_LEFT | TAT_TITLE,
1211                 var_to_string (dependent_var[i])
1212                 );
1213
1214
1215       if ( !fctr )
1216         populate_extremes (tbl, heading_columns - 2,
1217                           i * 2 * n_extremities * n_factors  + heading_rows,
1218                           n_extremities, &totals[i]);
1219
1220       else
1221         {
1222           struct factor_statistics **fs = fctr->fs;
1223           int count = 0 ;
1224           const union value *prev  = NULL;
1225
1226           while (*fs)
1227             {
1228               const int row = heading_rows + ( 2 * n_extremities )  *
1229                  ( ( i  * n_factors  ) +  count );
1230
1231
1232               if ( !prev || 0 != compare_values (prev, (*fs)->id[0],
1233                                         var_get_width (fctr->indep_var[0])))
1234                 {
1235                   struct string vstr;
1236                   ds_init_empty (&vstr);
1237                   var_append_value_name (fctr->indep_var[0],
1238                                       (*fs)->id[0], &vstr);
1239
1240                   if ( count > 0 )
1241                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
1242
1243                   tab_text (tbl,
1244                             1, row,
1245                             TAB_LEFT | TAT_TITLE,
1246                             ds_cstr (&vstr)
1247                             );
1248
1249                   ds_destroy (&vstr);
1250                 }
1251
1252               prev = (*fs)->id[0];
1253
1254               if (fctr->indep_var[1] && count > 0 )
1255                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
1256
1257               if ( fctr->indep_var[1])
1258                 {
1259                   struct string vstr;
1260                   ds_init_empty (&vstr);
1261                   var_append_value_name (fctr->indep_var[1], (*fs)->id[1], &vstr);
1262
1263                 tab_text (tbl, 2, row,
1264                           TAB_LEFT | TAT_TITLE,
1265                             ds_cstr (&vstr)
1266                           );
1267
1268                   ds_destroy (&vstr);
1269                 }
1270
1271               populate_extremes (tbl, heading_columns - 2,
1272                                 row, n_extremities,
1273                                 & (*fs)->m[i]);
1274
1275               count++ ;
1276               fs++;
1277             }
1278         }
1279     }
1280
1281   tab_submit (tbl);
1282 }
1283
1284
1285
1286 /* Fill in the extremities table */
1287 void
1288 populate_extremes (struct tab_table *t,
1289                   int col, int row, int n, const struct metrics *m)
1290 {
1291   int extremity;
1292   int idx=0;
1293
1294
1295   tab_text (t, col, row,
1296            TAB_RIGHT | TAT_TITLE ,
1297            _ ("Highest")
1298            );
1299
1300   tab_text (t, col, row + n ,
1301            TAB_RIGHT | TAT_TITLE ,
1302            _ ("Lowest")
1303            );
1304
1305
1306   tab_hline (t, TAL_1, col, col + 3, row + n );
1307
1308   for (extremity = 0; extremity < n ; ++extremity )
1309     {
1310       /* Highest */
1311       tab_float (t, col + 1, row + extremity,
1312                 TAB_RIGHT,
1313                 extremity + 1, 8, 0);
1314
1315
1316       /* Lowest */
1317       tab_float (t, col + 1, row + extremity + n,
1318                 TAB_RIGHT,
1319                 extremity + 1, 8, 0);
1320
1321     }
1322
1323
1324   /* Lowest */
1325   for (idx = 0, extremity = 0; extremity < n && idx < m->n_data ; ++idx )
1326     {
1327       int j;
1328       const struct weighted_value *wv = m->wvp[idx];
1329       struct case_node *cn = wv->case_nos;
1330
1331
1332       for (j = 0 ; j < wv->w ; ++j  )
1333         {
1334           if ( extremity + j >= n )
1335             break ;
1336
1337           tab_float (t, col + 3, row + extremity + j  + n,
1338                     TAB_RIGHT,
1339                     wv->v.f, 8, 2);
1340
1341           tab_float (t, col + 2, row + extremity + j  + n,
1342                     TAB_RIGHT,
1343                     cn->num, 8, 0);
1344
1345           if ( cn->next )
1346             cn = cn->next;
1347
1348         }
1349
1350       extremity +=  wv->w ;
1351     }
1352
1353
1354   /* Highest */
1355   for (idx = m->n_data - 1, extremity = 0; extremity < n && idx >= 0; --idx )
1356     {
1357       int j;
1358       const struct weighted_value *wv = m->wvp[idx];
1359       struct case_node *cn = wv->case_nos;
1360
1361       for (j = 0 ; j < wv->w ; ++j  )
1362         {
1363           if ( extremity + j >= n )
1364             break ;
1365
1366           tab_float (t, col + 3, row + extremity + j,
1367                     TAB_RIGHT,
1368                     wv->v.f, 8, 2);
1369
1370           tab_float (t, col + 2, row + extremity + j,
1371                     TAB_RIGHT,
1372                     cn->num, 8, 0);
1373
1374           if ( cn->next )
1375             cn = cn->next;
1376
1377         }
1378
1379       extremity +=  wv->w ;
1380     }
1381 }
1382
1383
1384 /* Show the descriptives table */
1385 void
1386 show_descriptives (const struct variable **dependent_var,
1387                   int n_dep_var,
1388                   struct factor *fctr)
1389 {
1390   int i;
1391   int heading_columns ;
1392   int n_cols;
1393   const int n_stat_rows = 13;
1394
1395   const int heading_rows = 1;
1396
1397   struct tab_table *tbl;
1398
1399   int n_factors = 1;
1400   int n_rows ;
1401
1402   if ( fctr )
1403     {
1404       heading_columns = 4;
1405       n_factors = hsh_count (fctr->fstats);
1406
1407       n_rows = n_dep_var * n_stat_rows * n_factors;
1408
1409       if ( fctr->indep_var[1] )
1410         heading_columns = 5;
1411     }
1412   else
1413     {
1414       heading_columns = 3;
1415       n_rows = n_dep_var * n_stat_rows;
1416     }
1417
1418   n_rows += heading_rows;
1419
1420   n_cols = heading_columns + 2;
1421
1422
1423   tbl = tab_create (n_cols, n_rows, 0);
1424
1425   tab_headers (tbl, heading_columns + 1, 0, heading_rows, 0);
1426
1427   tab_dim (tbl, tab_natural_dimensions);
1428
1429   /* Outline the box and have no internal lines*/
1430   tab_box (tbl,
1431            TAL_2, TAL_2,
1432            -1, -1,
1433            0, 0,
1434            n_cols - 1, n_rows - 1);
1435
1436   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1437
1438   tab_vline (tbl, TAL_1, heading_columns, 0, n_rows - 1);
1439   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows - 1);
1440   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1441
1442   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _ ("Statistic"));
1443   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _ ("Std. Error"));
1444
1445   tab_title (tbl, _ ("Descriptives"));
1446
1447
1448   for ( i = 0 ; i < n_dep_var ; ++i )
1449     {
1450       const int row = heading_rows + i * n_stat_rows * n_factors ;
1451
1452       if ( i > 0 )
1453         tab_hline (tbl, TAL_1, 0, n_cols - 1, row );
1454
1455       tab_text (tbl, 0,
1456                 i * n_stat_rows * n_factors  + heading_rows,
1457                 TAB_LEFT | TAT_TITLE,
1458                 var_to_string (dependent_var[i])
1459                 );
1460
1461
1462       if ( fctr  )
1463         {
1464           const union value *prev = NULL;
1465
1466           struct factor_statistics **fs = fctr->fs;
1467           int count = 0;
1468
1469           tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1470                     var_to_string (fctr->indep_var[0]));
1471
1472
1473           if ( fctr->indep_var[1])
1474             tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1475                       var_to_string (fctr->indep_var[1]));
1476
1477           while ( *fs )
1478             {
1479               const int row = heading_rows + n_stat_rows  *
1480                  ( ( i  * n_factors  ) +  count );
1481
1482
1483               if ( !prev || 0 != compare_values (prev, (*fs)->id[0],
1484                                         var_get_width (fctr->indep_var[0])))
1485                 {
1486                   struct string vstr;
1487                   ds_init_empty (&vstr);
1488                   var_append_value_name (fctr->indep_var[0],
1489                                       (*fs)->id[0], &vstr);
1490
1491                   if ( count > 0 )
1492                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
1493
1494                   tab_text (tbl,
1495                             1, row,
1496                             TAB_LEFT | TAT_TITLE,
1497                             ds_cstr (&vstr)
1498                             );
1499
1500                   ds_destroy (&vstr);
1501                 }
1502
1503               prev = (*fs)->id[0];
1504
1505               if (fctr->indep_var[1] && count > 0 )
1506                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
1507
1508               if ( fctr->indep_var[1])
1509                 {
1510                   struct string vstr;
1511                   ds_init_empty (&vstr);
1512                   var_append_value_name (fctr->indep_var[1], (*fs)->id[1], &vstr);
1513
1514                 tab_text (tbl, 2, row,
1515                           TAB_LEFT | TAT_TITLE,
1516                             ds_cstr (&vstr)
1517                           );
1518
1519                   ds_destroy (&vstr);
1520                 }
1521
1522               populate_descriptives (tbl, heading_columns - 2,
1523                                     row, & (*fs)->m[i]);
1524
1525               count++ ;
1526               fs++;
1527             }
1528
1529         }
1530
1531       else
1532         {
1533
1534           populate_descriptives (tbl, heading_columns - 2,
1535                                 i * n_stat_rows * n_factors  + heading_rows,
1536                                 &totals[i]);
1537         }
1538     }
1539
1540   tab_submit (tbl);
1541
1542 }
1543
1544
1545 /* Fill in the descriptives data */
1546 void
1547 populate_descriptives (struct tab_table *tbl, int col, int row,
1548                       const struct metrics *m)
1549 {
1550   const double t = gsl_cdf_tdist_Qinv ((1 - cmd.n_cinterval[0] / 100.0)/2.0,
1551                                       m->n -1);
1552
1553   tab_text (tbl, col,
1554             row,
1555             TAB_LEFT | TAT_TITLE,
1556             _ ("Mean"));
1557
1558   tab_float (tbl, col + 2,
1559              row,
1560              TAB_CENTER,
1561              m->mean,
1562              8,2);
1563
1564   tab_float (tbl, col + 3,
1565              row,
1566              TAB_CENTER,
1567              m->se_mean,
1568              8,3);
1569
1570
1571   tab_text (tbl, col,
1572             row + 1,
1573             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1574             _ ("%g%% Confidence Interval for Mean"), cmd.n_cinterval[0]);
1575
1576
1577   tab_text (tbl, col + 1,
1578             row  + 1,
1579             TAB_LEFT | TAT_TITLE,
1580             _ ("Lower Bound"));
1581
1582   tab_float (tbl, col + 2,
1583              row + 1,
1584              TAB_CENTER,
1585              m->mean - t * m->se_mean,
1586              8,3);
1587
1588   tab_text (tbl, col + 1,
1589             row + 2,
1590             TAB_LEFT | TAT_TITLE,
1591             _ ("Upper Bound"));
1592
1593
1594   tab_float (tbl, col + 2,
1595              row + 2,
1596              TAB_CENTER,
1597              m->mean + t * m->se_mean,
1598              8,3);
1599
1600   tab_text (tbl, col,
1601             row + 3,
1602             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1603             _ ("5%% Trimmed Mean"));
1604
1605   tab_float (tbl, col + 2,
1606              row + 3,
1607              TAB_CENTER,
1608              m->trimmed_mean,
1609              8,2);
1610
1611   tab_text (tbl, col,
1612             row + 4,
1613             TAB_LEFT | TAT_TITLE,
1614             _ ("Median"));
1615
1616   {
1617     struct percentile *p;
1618     double d = 50;
1619
1620     p = hsh_find (m->ptile_hash, &d);
1621
1622     assert (p);
1623
1624
1625     tab_float (tbl, col + 2,
1626                row + 4,
1627                TAB_CENTER,
1628                p->v,
1629                8, 2);
1630   }
1631
1632
1633   tab_text (tbl, col,
1634             row + 5,
1635             TAB_LEFT | TAT_TITLE,
1636             _ ("Variance"));
1637
1638   tab_float (tbl, col + 2,
1639              row + 5,
1640              TAB_CENTER,
1641              m->var,
1642              8,3);
1643
1644
1645   tab_text (tbl, col,
1646             row + 6,
1647             TAB_LEFT | TAT_TITLE,
1648             _ ("Std. Deviation"));
1649
1650
1651   tab_float (tbl, col + 2,
1652              row + 6,
1653              TAB_CENTER,
1654              m->stddev,
1655              8,3);
1656
1657
1658   tab_text (tbl, col,
1659             row + 7,
1660             TAB_LEFT | TAT_TITLE,
1661             _ ("Minimum"));
1662
1663   tab_float (tbl, col + 2,
1664              row + 7,
1665              TAB_CENTER,
1666              m->min,
1667              8,3);
1668
1669   tab_text (tbl, col,
1670             row + 8,
1671             TAB_LEFT | TAT_TITLE,
1672             _ ("Maximum"));
1673
1674   tab_float (tbl, col + 2,
1675              row + 8,
1676              TAB_CENTER,
1677              m->max,
1678              8,3);
1679
1680
1681   tab_text (tbl, col,
1682             row + 9,
1683             TAB_LEFT | TAT_TITLE,
1684             _ ("Range"));
1685
1686
1687   tab_float (tbl, col + 2,
1688              row + 9,
1689              TAB_CENTER,
1690              m->max - m->min,
1691              8,3);
1692
1693   tab_text (tbl, col,
1694             row + 10,
1695             TAB_LEFT | TAT_TITLE,
1696             _ ("Interquartile Range"));
1697
1698   {
1699     struct percentile *p1;
1700     struct percentile *p2;
1701
1702     double d = 75;
1703     p1 = hsh_find (m->ptile_hash, &d);
1704
1705     d = 25;
1706     p2 = hsh_find (m->ptile_hash, &d);
1707
1708     assert (p1);
1709     assert (p2);
1710
1711     tab_float (tbl, col + 2,
1712                row + 10,
1713                TAB_CENTER,
1714                p1->v - p2->v,
1715                8, 2);
1716   }
1717
1718
1719
1720   tab_text (tbl, col,
1721             row + 11,
1722             TAB_LEFT | TAT_TITLE,
1723             _ ("Skewness"));
1724
1725
1726   tab_float (tbl, col + 2,
1727              row + 11,
1728              TAB_CENTER,
1729              m->skewness,
1730              8,3);
1731
1732   /* stderr of skewness */
1733   tab_float (tbl, col + 3,
1734              row + 11,
1735              TAB_CENTER,
1736              calc_seskew (m->n),
1737              8,3);
1738
1739
1740   tab_text (tbl, col,
1741             row + 12,
1742             TAB_LEFT | TAT_TITLE,
1743             _ ("Kurtosis"));
1744
1745
1746   tab_float (tbl, col + 2,
1747              row + 12,
1748              TAB_CENTER,
1749              m->kurtosis,
1750              8,3);
1751
1752   /* stderr of kurtosis */
1753   tab_float (tbl, col + 3,
1754              row + 12,
1755              TAB_CENTER,
1756              calc_sekurt (m->n),
1757              8,3);
1758
1759
1760 }
1761
1762
1763
1764 void
1765 box_plot_variables (const struct factor *fctr,
1766                    const struct variable **vars, int n_vars,
1767                    const struct variable *id)
1768 {
1769
1770   int i;
1771   struct factor_statistics **fs ;
1772
1773   if ( ! fctr )
1774     {
1775       box_plot_group (fctr, vars, n_vars, id);
1776       return;
1777     }
1778
1779   for ( fs = fctr->fs ; *fs ; ++fs )
1780     {
1781       struct string str;
1782       double y_min = DBL_MAX;
1783       double y_max = -DBL_MAX;
1784       struct chart *ch = chart_create ();
1785       ds_init_empty (&str);
1786       factor_to_string (fctr, *fs, 0, &str );
1787
1788       chart_write_title (ch, ds_cstr (&str));
1789
1790       for ( i = 0 ; i < n_vars ; ++i )
1791         {
1792           y_max = MAX (y_max, (*fs)->m[i].max);
1793           y_min = MIN (y_min, (*fs)->m[i].min);
1794         }
1795
1796       boxplot_draw_yscale (ch, y_max, y_min);
1797
1798       for ( i = 0 ; i < n_vars ; ++i )
1799         {
1800
1801           const double box_width = (ch->data_right - ch->data_left)
1802             / (n_vars * 2.0 ) ;
1803
1804           const double box_centre = ( i * 2 + 1) * box_width
1805             + ch->data_left;
1806
1807           boxplot_draw_boxplot (ch,
1808                                box_centre, box_width,
1809                                & (*fs)->m[i],
1810                                var_to_string (vars[i]));
1811
1812
1813         }
1814
1815       chart_submit (ch);
1816       ds_destroy (&str);
1817     }
1818 }
1819
1820
1821
1822 /* Do a box plot, grouping all factors into one plot ;
1823    each dependent variable has its own plot.
1824 */
1825 void
1826 box_plot_group (const struct factor *fctr,
1827                const struct variable **vars,
1828                int n_vars,
1829                const struct variable *id UNUSED)
1830 {
1831
1832   int i;
1833
1834   for ( i = 0 ; i < n_vars ; ++i )
1835     {
1836       struct factor_statistics **fs ;
1837       struct chart *ch;
1838
1839       ch = chart_create ();
1840
1841       boxplot_draw_yscale (ch, totals[i].max, totals[i].min);
1842
1843       if ( fctr )
1844         {
1845           int n_factors = 0;
1846           int f=0;
1847           for ( fs = fctr->fs ; *fs ; ++fs )
1848             ++n_factors;
1849
1850           chart_write_title (ch, _ ("Boxplot of %s vs. %s"),
1851                             var_to_string (vars[i]), var_to_string (fctr->indep_var[0]) );
1852
1853           for ( fs = fctr->fs ; *fs ; ++fs )
1854             {
1855               struct string str;
1856               const double box_width = (ch->data_right - ch->data_left)
1857                 / (n_factors * 2.0 ) ;
1858
1859               const double box_centre = ( f++ * 2 + 1) * box_width
1860                 + ch->data_left;
1861
1862               ds_init_empty (&str);
1863               factor_to_string_concise (fctr, *fs, &str);
1864
1865               boxplot_draw_boxplot (ch,
1866                                    box_centre, box_width,
1867                                    & (*fs)->m[i],
1868                                    ds_cstr (&str));
1869               ds_destroy (&str);
1870             }
1871         }
1872       else if ( ch )
1873         {
1874           const double box_width = (ch->data_right - ch->data_left) / 3.0;
1875           const double box_centre = (ch->data_right + ch->data_left) / 2.0;
1876
1877           chart_write_title (ch, _ ("Boxplot"));
1878
1879           boxplot_draw_boxplot (ch,
1880                                box_centre,    box_width,
1881                                &totals[i],
1882                                var_to_string (vars[i]) );
1883
1884         }
1885
1886       chart_submit (ch);
1887     }
1888 }
1889
1890
1891 /* Plot the normal and detrended normal plots for m
1892    Label the plots with factorname */
1893 void
1894 np_plot (const struct metrics *m, const char *factorname)
1895 {
1896   int i;
1897   double yfirst=0, ylast=0;
1898
1899   /* Normal Plot */
1900   struct chart *np_chart;
1901
1902   /* Detrended Normal Plot */
1903   struct chart *dnp_chart;
1904
1905   /* The slope and intercept of the ideal normal probability line */
1906   const double slope = 1.0 / m->stddev;
1907   const double intercept = - m->mean / m->stddev;
1908
1909   /* Cowardly refuse to plot an empty data set */
1910   if ( m->n_data == 0 )
1911     return ;
1912
1913   np_chart = chart_create ();
1914   dnp_chart = chart_create ();
1915
1916   if ( !np_chart || ! dnp_chart )
1917     return ;
1918
1919   chart_write_title (np_chart, _ ("Normal Q-Q Plot of %s"), factorname);
1920   chart_write_xlabel (np_chart, _ ("Observed Value"));
1921   chart_write_ylabel (np_chart, _ ("Expected Normal"));
1922
1923
1924   chart_write_title (dnp_chart, _ ("Detrended Normal Q-Q Plot of %s"),
1925                     factorname);
1926   chart_write_xlabel (dnp_chart, _ ("Observed Value"));
1927   chart_write_ylabel (dnp_chart, _ ("Dev from Normal"));
1928
1929   yfirst = gsl_cdf_ugaussian_Pinv (m->wvp[0]->rank / ( m->n + 1));
1930   ylast =  gsl_cdf_ugaussian_Pinv (m->wvp[m->n_data-1]->rank / ( m->n + 1));
1931
1932
1933   {
1934     /* Need to make sure that both the scatter plot and the ideal fit into the
1935        plot */
1936     double x_lower = MIN (m->min, (yfirst - intercept) / slope) ;
1937     double x_upper = MAX (m->max, (ylast  - intercept) / slope) ;
1938     double slack = (x_upper - x_lower)  * 0.05 ;
1939
1940     chart_write_xscale (np_chart, x_lower - slack, x_upper + slack, 5);
1941
1942     chart_write_xscale (dnp_chart, m->min, m->max, 5);
1943
1944   }
1945
1946   chart_write_yscale (np_chart, yfirst, ylast, 5);
1947
1948   {
1949     /* We have to cache the detrended data, beacause we need to
1950        find its limits before we can plot it */
1951     double *d_data = xnmalloc (m->n_data, sizeof *d_data);
1952     double d_max = -DBL_MAX;
1953     double d_min = DBL_MAX;
1954     for ( i = 0 ; i < m->n_data; ++i )
1955       {
1956         const double ns = gsl_cdf_ugaussian_Pinv (m->wvp[i]->rank / ( m->n + 1));
1957
1958         chart_datum (np_chart, 0, m->wvp[i]->v.f, ns);
1959
1960         d_data[i] = (m->wvp[i]->v.f - m->mean) / m->stddev  - ns;
1961
1962         if ( d_data[i] < d_min ) d_min = d_data[i];
1963         if ( d_data[i] > d_max ) d_max = d_data[i];
1964       }
1965     chart_write_yscale (dnp_chart, d_min, d_max, 5);
1966
1967     for ( i = 0 ; i < m->n_data; ++i )
1968       chart_datum (dnp_chart, 0, m->wvp[i]->v.f, d_data[i]);
1969
1970     free (d_data);
1971   }
1972
1973   chart_line (np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
1974   chart_line (dnp_chart, 0, 0, m->min, m->max , CHART_DIM_X);
1975
1976   chart_submit (np_chart);
1977   chart_submit (dnp_chart);
1978 }
1979
1980
1981
1982
1983 /* Show the percentiles */
1984 void
1985 show_percentiles (const struct variable **dependent_var,
1986                  int n_dep_var,
1987                  struct factor *fctr)
1988 {
1989   struct tab_table *tbl;
1990   int i;
1991
1992   int n_cols, n_rows;
1993   int n_factors;
1994
1995   struct hsh_table *ptiles ;
1996
1997   int n_heading_columns;
1998   const int n_heading_rows = 2;
1999   const int n_stat_rows = 2;
2000
2001   int n_ptiles ;
2002
2003   if ( fctr )
2004     {
2005       struct factor_statistics **fs = fctr->fs ;
2006       n_heading_columns = 3;
2007       n_factors = hsh_count (fctr->fstats);
2008
2009       ptiles = (*fs)->m[0].ptile_hash;
2010
2011       if ( fctr->indep_var[1] )
2012         n_heading_columns = 4;
2013     }
2014   else
2015     {
2016       n_factors = 1;
2017       n_heading_columns = 2;
2018
2019       ptiles = totals[0].ptile_hash;
2020     }
2021
2022   n_ptiles = hsh_count (ptiles);
2023
2024   n_rows = n_heading_rows + n_dep_var * n_stat_rows * n_factors;
2025
2026   n_cols = n_heading_columns + n_ptiles ;
2027
2028   tbl = tab_create (n_cols, n_rows, 0);
2029
2030   tab_headers (tbl, n_heading_columns + 1, 0, n_heading_rows, 0);
2031
2032   tab_dim (tbl, tab_natural_dimensions);
2033
2034   /* Outline the box and have no internal lines*/
2035   tab_box (tbl,
2036            TAL_2, TAL_2,
2037            -1, -1,
2038            0, 0,
2039            n_cols - 1, n_rows - 1);
2040
2041   tab_hline (tbl, TAL_2, 0, n_cols - 1, n_heading_rows );
2042
2043   tab_vline (tbl, TAL_2, n_heading_columns, 0, n_rows - 1);
2044
2045
2046   tab_title (tbl, _ ("Percentiles"));
2047
2048
2049   tab_hline (tbl, TAL_1, n_heading_columns, n_cols - 1, 1 );
2050
2051
2052   tab_box (tbl,
2053            -1, -1,
2054            -1, TAL_1,
2055            0, n_heading_rows,
2056            n_heading_columns - 1, n_rows - 1);
2057
2058
2059   tab_box (tbl,
2060            -1, -1,
2061            -1, TAL_1,
2062            n_heading_columns, n_heading_rows - 1,
2063            n_cols - 1, n_rows - 1);
2064
2065   tab_joint_text (tbl, n_heading_columns + 1, 0,
2066                  n_cols - 1 , 0,
2067                  TAB_CENTER | TAT_TITLE ,
2068                  _ ("Percentiles"));
2069
2070
2071   {
2072     /* Put in the percentile break points as headings */
2073
2074     struct percentile **p = (struct percentile **) hsh_sort (ptiles);
2075
2076     i = 0;
2077     while ( (*p)  )
2078       {
2079         tab_float (tbl, n_heading_columns + i++ , 1,
2080                   TAB_CENTER,
2081                   (*p)->p, 8, 0);
2082
2083         p++;
2084       }
2085
2086   }
2087
2088   for ( i = 0 ; i < n_dep_var ; ++i )
2089     {
2090       const int n_stat_rows = 2;
2091       const int row = n_heading_rows + i * n_stat_rows * n_factors ;
2092
2093       if ( i > 0 )
2094         tab_hline (tbl, TAL_1, 0, n_cols - 1, row );
2095
2096       tab_text (tbl, 0,
2097                 i * n_stat_rows * n_factors  + n_heading_rows,
2098                 TAB_LEFT | TAT_TITLE,
2099                 var_to_string (dependent_var[i])
2100                 );
2101
2102       if ( fctr  )
2103         {
2104           const union value *prev  = NULL ;
2105           struct factor_statistics **fs = fctr->fs;
2106           int count = 0;
2107
2108           tab_text (tbl, 1, n_heading_rows - 1,
2109                     TAB_CENTER | TAT_TITLE,
2110                     var_to_string (fctr->indep_var[0]));
2111
2112
2113           if ( fctr->indep_var[1])
2114             tab_text (tbl, 2, n_heading_rows - 1, TAB_CENTER | TAT_TITLE,
2115                       var_to_string (fctr->indep_var[1]));
2116
2117           while ( *fs )
2118             {
2119               const int row = n_heading_rows + n_stat_rows  *
2120                  ( ( i  * n_factors  ) +  count );
2121
2122
2123               if ( !prev || 0 != compare_values (prev, (*fs)->id[0],
2124                                         var_get_width (fctr->indep_var[0])))
2125                 {
2126                   struct string vstr;
2127                   ds_init_empty (&vstr);
2128                   var_append_value_name (fctr->indep_var[0],
2129                                       (*fs)->id[0], &vstr);
2130
2131
2132                   if ( count > 0 )
2133                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
2134
2135                   tab_text (tbl,
2136                             1, row,
2137                             TAB_LEFT | TAT_TITLE,
2138                             ds_cstr (&vstr)
2139                             );
2140
2141                   ds_destroy (&vstr);
2142                 }
2143
2144               prev = (*fs)->id[0];
2145
2146               if (fctr->indep_var[1] && count > 0 )
2147                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
2148
2149               if ( fctr->indep_var[1])
2150                 {
2151                   struct string vstr;
2152                   ds_init_empty (&vstr);
2153                   var_append_value_name (fctr->indep_var[1], (*fs)->id[1], &vstr);
2154
2155                 tab_text (tbl, 2, row,
2156                           TAB_LEFT | TAT_TITLE,
2157                             ds_cstr (&vstr)
2158                           );
2159
2160                   ds_destroy (&vstr);
2161                 }
2162
2163
2164               populate_percentiles (tbl, n_heading_columns - 1,
2165                                    row, & (*fs)->m[i]);
2166
2167
2168               count++ ;
2169               fs++;
2170             }
2171
2172
2173         }
2174       else
2175         {
2176           populate_percentiles (tbl, n_heading_columns - 1,
2177                                i * n_stat_rows * n_factors  + n_heading_rows,
2178                                &totals[i]);
2179         }
2180
2181
2182     }
2183
2184
2185   tab_submit (tbl);
2186
2187
2188 }
2189
2190
2191
2192
2193 void
2194 populate_percentiles (struct tab_table *tbl, int col, int row,
2195                      const struct metrics *m)
2196 {
2197   int i;
2198
2199   struct percentile **p = (struct percentile **) hsh_sort (m->ptile_hash);
2200
2201   tab_text (tbl,
2202             col, row + 1,
2203             TAB_LEFT | TAT_TITLE,
2204             _ ("Tukey\'s Hinges")
2205             );
2206
2207   tab_text (tbl,
2208             col, row,
2209             TAB_LEFT | TAT_TITLE,
2210             ptile_alg_desc[m->ptile_alg]
2211             );
2212
2213
2214   i = 0;
2215   while ( (*p)  )
2216     {
2217       tab_float (tbl, col + i + 1 , row,
2218                 TAB_CENTER,
2219                  (*p)->v, 8, 2);
2220       if ( (*p)->p == 25 )
2221         tab_float (tbl, col + i + 1 , row + 1,
2222                   TAB_CENTER,
2223                   m->hinge[0], 8, 2);
2224
2225       if ( (*p)->p == 50 )
2226         tab_float (tbl, col + i + 1 , row + 1,
2227                   TAB_CENTER,
2228                   m->hinge[1], 8, 2);
2229
2230       if ( (*p)->p == 75 )
2231         tab_float (tbl, col + i + 1 , row + 1,
2232                   TAB_CENTER,
2233                   m->hinge[2], 8, 2);
2234
2235
2236       i++;
2237
2238       p++;
2239     }
2240
2241 }
2242
2243 static void
2244 factor_to_string (const struct factor *fctr,
2245                   const struct factor_statistics *fs,
2246                   const struct variable *var,
2247                   struct string *str
2248                   )
2249 {
2250   if (var)
2251     ds_put_format (str, "%s (",var_to_string (var) );
2252
2253
2254   ds_put_format (str,  "%s = ",
2255                  var_to_string (fctr->indep_var[0]));
2256
2257   var_append_value_name (fctr->indep_var[0], fs->id[0], str);
2258
2259   if ( fctr->indep_var[1] )
2260     {
2261       ds_put_format (str, "; %s = )",
2262                      var_to_string (fctr->indep_var[1]));
2263
2264       var_append_value_name (fctr->indep_var[1], fs->id[1], str);
2265     }
2266   else
2267     {
2268       if ( var )
2269         ds_put_cstr (str, ")");
2270     }
2271 }
2272
2273
2274 static void
2275 factor_to_string_concise (const struct factor *fctr,
2276                           const struct factor_statistics *fs,
2277                           struct string *str
2278                           )
2279
2280 {
2281   var_append_value_name (fctr->indep_var[0], fs->id[0], str);
2282
2283   if ( fctr->indep_var[1] )
2284     {
2285       ds_put_cstr (str, ",");
2286
2287       var_append_value_name (fctr->indep_var[1],fs->id[1], str);
2288
2289       ds_put_cstr (str, ")");
2290     }
2291 }
2292
2293 /*
2294   Local Variables:
2295   mode: c
2296   End:
2297 */