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