e7a4f584b8cef35f6653f6cf511b39741b45ece7
[pspp-builds.git] / src / language / stats / means.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2011, 2012 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 "data/case.h"
20 #include "data/casegrouper.h"
21 #include "data/casereader.h"
22 #include "data/dataset.h"
23 #include "data/dictionary.h"
24 #include "data/format.h"
25 #include "data/variable.h"
26
27 #include "language/command.h"
28 #include "language/lexer/lexer.h"
29 #include "language/lexer/variable-parser.h"
30
31 #include "libpspp/misc.h"
32 #include "libpspp/pool.h"
33
34 #include "math/categoricals.h"
35 #include "math/interaction.h"
36 #include "math/moments.h"
37
38 #include "output/tab.h"
39
40 #include <math.h>
41
42 #include "gettext.h"
43 #define _(msgid) gettext (msgid)
44 #define N_(msgid) (msgid)
45
46
47 struct means;
48
49 struct per_var_data
50 {
51   void **cell_stats;
52   struct moments1 *mom;
53 };
54
55
56 typedef void *stat_create (const struct means *, struct pool *pool);
57
58 typedef void stat_update (const struct means *, void *stat, double w, double x);
59
60 typedef double stat_get (const struct means *, struct per_var_data *, void *aux);
61
62
63 struct cell_spec
64 {
65   /* Printable title for output */
66   const char *title;
67
68   /* Keyword for syntax */
69   const char *keyword;
70
71   stat_create *sc;
72   stat_update *su;
73   stat_get *sd;
74 };
75
76
77 struct harmonic_mean
78 {
79   double rsum;
80   double n;
81 };
82
83
84 static void *
85 harmonic_create (const struct means *means UNUSED, struct pool *pool)
86 {
87   struct harmonic_mean *hm = pool_alloc (pool, sizeof *hm);
88
89   hm->rsum = 0;
90   hm->n = 0;
91
92   return hm;
93 }
94
95
96 static void
97 harmonic_update (const struct means *means UNUSED, void *stat, double w, double x)
98 {
99   struct harmonic_mean *hm = stat;
100   hm->rsum  += w / x;
101   hm->n += w;
102 }
103
104
105 static double
106 harmonic_get (const struct means *means UNUSED, struct per_var_data *pvd UNUSED, void *stat)
107 {
108   struct harmonic_mean *hm = stat;
109
110   return hm->n / hm->rsum;
111 }
112
113 \f
114
115 struct geometric_mean
116 {
117   double prod;
118   double n;
119 };
120
121
122 static void *
123 geometric_create (const struct means *means UNUSED, struct pool *pool)
124 {
125   struct geometric_mean *gm = pool_alloc (pool, sizeof *gm);
126
127   gm->prod = 1.0;
128   gm->n = 0;
129
130   return gm;
131 }
132
133
134 static void
135 geometric_update (const struct means *means UNUSED, void *stat, double w, double x)
136 {
137   struct geometric_mean *gm = stat;
138   gm->prod  *=  pow (x, w);
139   gm->n += w;
140 }
141
142
143 static double
144 geometric_get (const struct means *means UNUSED, struct per_var_data *pvd UNUSED, void *stat)
145 {
146   struct geometric_mean *gm = stat;
147
148   return pow (gm->prod, 1.0 / gm->n);
149 }
150
151 \f
152
153 static double
154 sum_get (const struct means *means UNUSED, struct per_var_data *pvd, void *stat UNUSED)
155 {
156   double n, mean;
157
158   moments1_calculate (pvd->mom, &n, &mean, 0, 0, 0);
159
160   return mean * n;
161 }
162
163
164 static double
165 n_get (const struct means *means UNUSED, struct per_var_data *pvd, void *stat UNUSED)
166 {
167   double n;
168
169   moments1_calculate (pvd->mom, &n, 0, 0, 0, 0);
170
171   return n;
172 }
173
174 static double
175 arithmean_get (const const struct means *means UNUSED, struct per_var_data *pvd, void *stat UNUSED)
176 {
177   double n, mean;
178
179   moments1_calculate (pvd->mom, &n, &mean, 0, 0, 0);
180
181   return mean;
182 }
183
184 static double
185 variance_get (const const struct means *means UNUSED, struct per_var_data *pvd, void *stat UNUSED)
186 {
187   double n, mean, variance;
188
189   moments1_calculate (pvd->mom, &n, &mean, &variance, 0, 0);
190
191   return variance;
192 }
193
194
195 static double
196 stddev_get (const const struct means *means UNUSED, struct per_var_data *pvd, void *stat)
197 {
198   return sqrt (variance_get (means, pvd, stat));
199 }
200
201
202 \f
203
204 static double
205 skew_get (const const struct means *means UNUSED, struct per_var_data *pvd, void *stat UNUSED)
206 {
207   double skew;
208
209   moments1_calculate (pvd->mom, NULL, NULL, NULL, &skew, 0);
210
211   return skew;
212 }
213
214 static double
215 sekurt_get (const const struct means *means UNUSED, struct per_var_data *pvd, void *stat UNUSED)
216 {
217   double n;
218
219   moments1_calculate (pvd->mom, &n, NULL, NULL, NULL, NULL);
220
221   return calc_sekurt (n);
222 }
223
224
225
226 static double
227 seskew_get (const const struct means *means UNUSED, struct per_var_data *pvd, void *stat UNUSED)
228 {
229   double n;
230
231   moments1_calculate (pvd->mom, &n, NULL, NULL, NULL, NULL);
232
233   return calc_seskew (n);
234 }
235
236
237
238 static double
239 kurt_get (const const struct means *means UNUSED, struct per_var_data *pvd, void *stat UNUSED)
240 {
241   double kurt;
242
243   moments1_calculate (pvd->mom, NULL, NULL, NULL, NULL, &kurt);
244
245   return kurt;
246 }
247
248
249 static double
250 semean_get (const struct means *means UNUSED, struct per_var_data *pvd, void *stat UNUSED)
251 {
252   double n, var;
253
254   moments1_calculate (pvd->mom, &n, NULL, &var, NULL, NULL);
255
256   return sqrt (var / n);
257 }
258
259
260
261 \f
262
263 static void *
264 min_create (const struct means *means UNUSED, struct pool *pool)
265 {
266   double *r = pool_alloc (pool, sizeof *r);
267
268   *r = DBL_MAX;
269
270   return r;
271 }
272
273 static void
274 min_update (const struct means *means UNUSED, void *stat, double w UNUSED, double x)
275 {
276   double *r = stat;
277
278   if (x < *r)
279     *r = x;
280 }
281
282 static double
283 min_get (const const struct means *means UNUSED, struct per_var_data *pvd UNUSED, void *stat)
284 {
285   double *r = stat;
286
287   return *r;
288 }
289
290 static void *
291 max_create (const struct means *means UNUSED, struct pool *pool)
292 {
293   double *r = pool_alloc (pool, sizeof *r);
294
295   *r = -DBL_MAX;
296
297   return r;
298 }
299
300 static void
301 max_update (const struct means *means UNUSED, void *stat, double w UNUSED, double x)
302 {
303   double *r = stat;
304
305   if (x > *r)
306     *r = x;
307 }
308
309 static double
310 max_get (const const struct means *means UNUSED, struct per_var_data *pvd UNUSED, void *stat)
311 {
312   double *r = stat;
313
314   return *r;
315 }
316
317 \f
318
319 struct range
320 {
321   double min;
322   double max;
323 };
324
325 static void *
326 range_create (const struct means *means UNUSED, struct pool *pool)
327 {
328   struct range *r = pool_alloc (pool, sizeof *r);
329
330   r->min = DBL_MAX;
331   r->max = -DBL_MAX;
332
333   return r;
334 }
335
336 static void
337 range_update (const struct means *means UNUSED, void *stat, double w UNUSED, double x)
338 {
339   struct range *r = stat;
340
341   if (x > r->max)
342     r->max = x;
343
344   if (x < r->min)
345     r->min = x;
346 }
347
348 static double
349 range_get (const const struct means *means UNUSED, struct per_var_data *pvd UNUSED, void *stat)
350 {
351   struct range *r = stat;
352
353   return r->max - r->min;
354 }
355
356 \f
357
358 static void *
359 last_create (const struct means *means UNUSED, struct pool *pool)
360 {
361   double *l = pool_alloc (pool, sizeof *l);
362
363   return l;
364 }
365
366 static void
367 last_update (const struct means *means UNUSED, void *stat, double w UNUSED, double x)
368 {
369   double *l = stat;
370
371   *l = x;
372 }
373
374 static double
375 last_get (const const struct means *means UNUSED, struct per_var_data *pvd UNUSED, void *stat)
376 {
377   double *l = stat;
378
379   return *l;
380 }
381
382
383 static void *
384 first_create (const struct means *means UNUSED, struct pool *pool)
385 {
386   double *f = pool_alloc (pool, sizeof *f);
387
388    *f = SYSMIS;
389
390   return f;
391 }
392
393 static void
394 first_update (const struct means *means UNUSED, void *stat, double w UNUSED, double x)
395 {
396   double *f = stat;
397
398   if (*f == SYSMIS)
399      *f = x;
400 }
401
402 static double
403 first_get (const const struct means *means UNUSED, struct per_var_data *pvd UNUSED,  void *stat)
404 {
405   double *f = stat;
406
407   return *f;
408 }
409
410
411
412 /* Table of cell_specs */
413 static const struct cell_spec cell_spec[] = {
414   {N_("Mean"), "MEAN", NULL, NULL, arithmean_get},
415   {N_("N"), "COUNT", NULL, NULL, n_get},
416   {N_("Std. Deviation"), "STDDEV", NULL, NULL, stddev_get},
417 #if 0
418   {N_("Median"), "MEDIAN", NULL, NULL, NULL},
419   {N_("Group Median"), "GMEDIAN", NULL, NULL, NULL},
420 #endif
421   {N_("S.E. Mean"), "SEMEAN", NULL, NULL, semean_get},
422   {N_("Sum"), "SUM", NULL, NULL, sum_get},
423   {N_("Min"), "MIN", min_create, min_update, min_get},
424   {N_("Max"), "MAX", max_create, max_update, max_get},
425   {N_("Range"), "RANGE", range_create, range_update, range_get},
426   {N_("Variance"), "VARIANCE", NULL, NULL, variance_get},
427   {N_("Kurtosis"), "KURT", NULL, NULL, kurt_get},
428   {N_("S.E. Kurt"), "SEKURT", NULL, NULL, sekurt_get},
429   {N_("Skewness"), "SKEW", NULL, NULL, skew_get},
430   {N_("S.E. Skew"), "SESKEW", NULL, NULL, seskew_get},
431   {N_("First"), "FIRST", first_create, first_update, first_get},
432   {N_("Last"), "LAST", last_create, last_update, last_get},
433 #if 0
434   {N_("Percent N"), "NPCT", NULL, NULL, NULL},
435   {N_("Percent Sum"), "SPCT", NULL, NULL, NULL},
436 #endif
437   {N_("Harmonic Mean"), "HARMONIC", harmonic_create, harmonic_update, harmonic_get},
438   {N_("Geom. Mean"), "GEOMETRIC", geometric_create, geometric_update, geometric_get}
439 };
440
441 #define n_C (sizeof (cell_spec) / sizeof (struct cell_spec))
442
443
444 struct summary
445 {
446   casenumber missing;
447   casenumber non_missing;
448 };
449
450
451 /* The thing parsed after TABLES= */
452 struct mtable
453 {
454   size_t n_dep_vars;
455   const struct variable **dep_vars;
456
457   size_t n_interactions;
458   struct interaction **interactions;
459   struct summary *summary;
460
461   size_t *n_factor_vars;
462   const struct variable ***factor_vars;
463
464   int ii;
465
466   int n_layers;
467
468   struct categoricals *cats;
469 };
470
471 struct means
472 {
473   const struct dictionary *dict;
474
475   struct mtable *table;
476   size_t n_tables;
477
478   /* Missing value class for categorical variables */
479   enum mv_class exclude;
480
481   /* Missing value class for dependent variables */
482   enum mv_class dep_exclude;
483
484   /* an array  indicating which statistics are to be calculated */
485   int *cells;
486
487   /* Size of cells */
488   int n_cells;
489
490   /* Pool on which cell functions may allocate data */
491   struct pool *pool;
492 };
493
494
495 static void
496 run_means (struct means *cmd, struct casereader *input,
497            const struct dataset *ds);
498
499 /* Append all the variables belonging to layer and all subsequent layers
500    to iact. And then append iact to the means->interaction.
501    This is a recursive function.
502  */
503 static void
504 iact_append_factor (struct mtable *means, int layer,
505                     const struct interaction *iact)
506 {
507   int v;
508   const struct variable **fv;
509
510   if (layer >= means->n_layers)
511     return;
512
513   fv = means->factor_vars[layer];
514
515   for (v = 0; v < means->n_factor_vars[layer]; ++v)
516     {
517       struct interaction *nexti = interaction_clone (iact);
518
519       interaction_add_variable (nexti, fv[v]);
520
521       iact_append_factor (means, layer + 1, nexti);
522
523       if (layer == means->n_layers - 1)
524         {
525           means->interactions[means->ii++] = nexti;
526         }
527     }
528 }
529
530 static bool
531 parse_means_table_syntax (struct lexer *lexer, const struct means *cmd, struct mtable *table)
532 {
533   table->ii = 0;
534   table->n_layers = 0;
535   table->factor_vars = NULL;
536   table->n_factor_vars = NULL;
537
538   /* Dependent variable (s) */
539   if (!parse_variables_const (lexer, cmd->dict,
540                               &table->dep_vars, &table->n_dep_vars,
541                               PV_NO_DUPLICATE | PV_NUMERIC))
542     return false;
543
544   /* Factor variable (s) */
545   while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
546     {
547       if (lex_match (lexer, T_BY))
548         {
549           table->n_layers++;
550           table->factor_vars =
551             xrealloc (table->factor_vars,
552                       sizeof (*table->factor_vars) * table->n_layers);
553
554           table->n_factor_vars =
555             xrealloc (table->n_factor_vars,
556                       sizeof (*table->n_factor_vars) * table->n_layers);
557
558           if (!parse_variables_const (lexer, cmd->dict,
559                                       &table->factor_vars[table->n_layers - 1],
560                                       &table->n_factor_vars[table->n_layers -
561                                                             1],
562                                       PV_NO_DUPLICATE))
563             return false;
564
565         }
566     }
567
568   return true;
569 }
570
571 /* Match a variable.
572    If the match succeeds, the variable will be placed in VAR.
573    Returns true if successful */
574 static bool
575 lex_is_variable (struct lexer *lexer, const struct dictionary *dict,
576                  int n)
577 {
578   const char *tstr;
579   if (lex_next_token (lexer, n) !=  T_ID)
580     return false;
581
582   tstr = lex_next_tokcstr (lexer, n);
583
584   if (NULL == dict_lookup_var (dict, tstr) )
585     return false;
586
587   return true;
588 }
589
590
591 int
592 cmd_means (struct lexer *lexer, struct dataset *ds)
593 {
594   int t;
595   int i;
596   int l;
597   struct means means;
598   bool more_tables = true;
599
600   means.exclude = MV_ANY;
601   means.dep_exclude = MV_ANY;
602   means.table = NULL;
603   means.n_tables = 0;
604
605   means.dict = dataset_dict (ds);
606
607   means.n_cells = 3;
608   means.cells = xcalloc (means.n_cells, sizeof (*means.cells));
609
610
611   /* The first three items (MEAN, COUNT, STDDEV) are the default */
612   for (i = 0; i < 3; ++i)
613     means.cells[i] = i;
614
615
616   /*   Optional TABLES =   */
617   if (lex_match_id (lexer, "TABLES"))
618     {
619       lex_force_match (lexer, T_EQUALS);
620     }
621
622
623   more_tables = true;
624   /* Parse the "tables" */
625   while (more_tables)
626     {
627       means.n_tables ++;
628       means.table = xrealloc (means.table, means.n_tables * sizeof (*means.table));
629
630       if (! parse_means_table_syntax (lexer, &means, 
631                                       &means.table[means.n_tables - 1]))
632         {
633           goto error;
634         }
635
636       /* Look ahead to see if there are more tables to be parsed */
637       more_tables = false;
638       if ( T_SLASH == lex_next_token (lexer, 0) )
639         {
640           if (lex_is_variable (lexer, means.dict, 1) )
641             {
642               more_tables = true;
643               lex_force_match (lexer, T_SLASH);
644             }
645         }
646     }
647
648   /* /MISSING subcommand */
649   while (lex_token (lexer) != T_ENDCMD)
650     {
651       lex_match (lexer, T_SLASH);
652
653       if (lex_match_id (lexer, "MISSING"))
654         {
655           /*
656              If no MISSING subcommand is specified, each combination of
657              a dependent variable and categorical variables is handled
658              separately.
659            */
660           lex_match (lexer, T_EQUALS);
661           if (lex_match_id (lexer, "INCLUDE"))
662             {
663               /*
664                 Use the subcommand  "/MISSING=INCLUDE" to include user-missing
665                 values in the analysis.
666               */
667
668               means.exclude = MV_SYSTEM;
669               means.dep_exclude = MV_SYSTEM;
670             }
671           else if (lex_match_id (lexer, "TABLE"))
672             /*
673               This is the default. (I think).
674               Every case containing a complete set of variables for a given
675               table. If any variable, categorical or dependent for in a table
676               is missing (as defined by what?), then that variable will
677               be dropped FOR THAT TABLE ONLY.
678             */
679             {
680               means.exclude = MV_ANY;
681               means.dep_exclude = MV_ANY;
682             }
683           else if (lex_match_id (lexer, "DEPENDENT"))
684             /*
685              Use the command "/MISSING=DEPENDENT" to
686              include user-missing values for the categorical variables, 
687              while excluding them for the dependent variables.
688
689              Cases are dropped only when user-missing values
690              appear in dependent  variables.  User-missing
691              values for categorical variables are treated according to
692              their face value.
693
694              Cases are ALWAYS dropped when System Missing values appear 
695              in the categorical variables.
696             */
697             {
698               means.dep_exclude = MV_ANY;
699               means.exclude = MV_SYSTEM;
700             }
701           else
702             {
703               lex_error (lexer, NULL);
704               goto error;
705             }
706         }
707       else if (lex_match_id (lexer, "CELLS"))
708         {
709           lex_match (lexer, T_EQUALS);
710
711           /* The default values become overwritten */
712           means.n_cells = 0;
713           while (lex_token (lexer) != T_ENDCMD
714                  && lex_token (lexer) != T_SLASH)
715             {
716               int k;
717               for (k = 0; k < n_C; ++k)
718                 {
719                   if (lex_match_id (lexer, cell_spec[k].keyword))
720                     {
721                       means.cells =
722                         xrealloc (means.cells,
723                                   ++means.n_cells * sizeof (*means.cells));
724
725                       means.cells[means.n_cells - 1] = k;
726                       break;
727                     }
728                 }
729               if (k >= n_C)
730                 {
731                   lex_error (lexer, NULL);
732                   goto error;
733                 }
734             }
735         }
736       else
737         {
738           lex_error (lexer, NULL);
739           goto error;
740         }
741     }
742
743   means.pool = pool_create ();
744
745
746   for (t = 0; t < means.n_tables; ++t)
747   {
748     struct mtable *table = &means.table[t];
749     table->n_interactions = 1;
750     for (l = 0; l < table->n_layers; ++l)
751       {
752         const int n_vars = table->n_factor_vars[l];
753         table->n_interactions *= n_vars;
754       }
755
756     table->interactions =
757       xcalloc (table->n_interactions, sizeof (*table->interactions));
758
759     table->summary =
760       xcalloc (table->n_dep_vars * table->n_interactions, sizeof (*table->summary));
761
762
763     if (table->n_layers > 0)
764       iact_append_factor (table, 0, interaction_create (NULL));
765     else
766       table->interactions[0] = interaction_create (NULL);
767
768   }
769
770
771   {
772     struct casegrouper *grouper;
773     struct casereader *group;
774     bool ok;
775
776     grouper = casegrouper_create_splits (proc_open (ds), means.dict);
777     while (casegrouper_get_next_group (grouper, &group))
778       {
779         run_means (&means, group, ds);
780       }
781     ok = casegrouper_destroy (grouper);
782     ok = proc_commit (ds) && ok;
783   }
784
785
786   return CMD_SUCCESS;
787
788 error:
789
790   return CMD_FAILURE;
791 }
792
793
794 static bool
795 is_missing (const struct means *cmd,
796             const struct variable *dvar,
797             const struct interaction *iact,
798             const struct ccase *c)
799 {
800   if ( interaction_case_is_missing (iact, c, cmd->exclude) )
801     return true;
802
803
804   if (var_is_value_missing (dvar,
805                             case_data (c, dvar),
806                             cmd->dep_exclude))
807     return true;
808
809   return false;
810 }
811
812 static void output_case_processing_summary (const struct mtable *);
813
814 static void output_report (const struct means *, int, const struct mtable *);
815
816
817 struct per_cat_data
818 {
819   struct per_var_data *pvd;
820
821   bool warn;
822 };
823
824 static void *
825 create_n (const void *aux1, void *aux2)
826 {
827   int i, v;
828   const struct means *means = aux1;
829   struct mtable *table = aux2;
830   struct per_cat_data *per_cat_data = xmalloc (sizeof *per_cat_data);
831
832   struct per_var_data *pvd = xcalloc (table->n_dep_vars, sizeof *pvd);
833
834   for (v = 0; v < table->n_dep_vars; ++v)
835     {
836       enum moment maxmom = MOMENT_KURTOSIS;
837       struct per_var_data *pp = &pvd[v];
838
839       pp->cell_stats = xcalloc (means->n_cells, sizeof *pp->cell_stats);
840       
841
842       for (i = 0; i < means->n_cells; ++i)
843         {
844           int csi = means->cells[i];
845           const struct cell_spec *cs = &cell_spec[csi];
846           if (cs->sc)
847             {
848               pp->cell_stats[i] = cs->sc (means, means->pool);
849             }
850         }
851       pp->mom = moments1_create (maxmom);
852     }
853
854
855   per_cat_data->pvd = pvd;
856   per_cat_data->warn = true;
857   return per_cat_data;
858 }
859
860 static void
861 update_n (const void *aux1, void *aux2, void *user_data, const struct ccase *c)
862 {
863   int i;
864   int v = 0;
865   const struct means *means = aux1;
866   struct mtable *table = aux2;
867   struct per_cat_data *per_cat_data = user_data;
868
869   double weight = dict_get_case_weight (means->dict, c, &per_cat_data->warn);
870
871   for (v = 0; v < table->n_dep_vars; ++v)
872     {
873       struct per_var_data *pvd = &per_cat_data->pvd[v];
874
875       const double x = case_data (c, table->dep_vars[v])->f;
876
877       for (i = 0; i < table->n_interactions; ++i)
878         {
879           if ( is_missing (means, table->dep_vars[v], table->interactions[i], c))
880             goto end;
881         }
882
883       for (i = 0; i < means->n_cells; ++i)
884         {
885           const int csi = means->cells[i];
886           const struct cell_spec *cs = &cell_spec[csi];
887
888
889           if (cs->su)
890             cs->su (means,
891                     pvd->cell_stats[i],
892                     weight, x);
893         }
894
895       moments1_add (pvd->mom, x, weight);
896
897     end:
898       continue;
899     }
900 }
901
902 static void
903 calculate_n (const void *aux1, void *aux2, void *user_data)
904 {
905   int i;
906   int v = 0;
907   struct per_cat_data *per_cat_data = user_data;
908   const struct means *means = aux1;
909   struct mtable *table = aux2;
910
911   for (v = 0; v < table->n_dep_vars; ++v)
912     {
913       struct per_var_data *pvd = &per_cat_data->pvd[v];
914       for (i = 0; i < means->n_cells; ++i)
915         {
916           int csi = means->cells[i];
917           const struct cell_spec *cs = &cell_spec[csi];
918
919           if (cs->su)
920             cs->sd (means, pvd, pvd->cell_stats[i]);
921         }
922     }
923 }
924
925
926 static void
927 run_means (struct means *cmd, struct casereader *input,
928            const struct dataset *ds UNUSED)
929 {
930   int i,t;
931   const struct variable *wv = dict_get_weight (cmd->dict);
932   struct ccase *c;
933   struct casereader *reader;
934
935   struct payload payload;
936   payload.create = create_n;
937   payload.update = update_n;
938   payload.destroy = calculate_n;
939   
940   for (t = 0; t < cmd->n_tables; ++t)
941   {
942     struct mtable *table = &cmd->table[t];
943     table->cats
944       = categoricals_create (table->interactions,
945                              table->n_interactions, wv, cmd->exclude);
946
947     categoricals_set_payload (table->cats, &payload, cmd, table);
948   }
949
950   for (reader = casereader_clone (input);
951        (c = casereader_read (reader)) != NULL; case_unref (c))
952     {
953       for (t = 0; t < cmd->n_tables; ++t)
954         {
955           int  v;
956           struct mtable *table = &cmd->table[t];
957
958           for (v = 0; v < table->n_dep_vars; ++v)
959             {
960               int i;
961               for (i = 0; i < table->n_interactions; ++i)
962                 {
963                   const bool missing =
964                     is_missing (cmd, table->dep_vars[v],
965                                 table->interactions[i], c);
966                   if (missing)
967                     table->summary[v * table->n_interactions + i].missing++;
968                   else
969                     table->summary[v * table->n_interactions + i].non_missing++;
970                 }
971             }
972           categoricals_update (table->cats, c);
973         }
974     }
975   casereader_destroy (reader);
976
977   for (t = 0; t < cmd->n_tables; ++t)
978     {
979       struct mtable *table = &cmd->table[t];
980
981       categoricals_done (table->cats);
982     }
983
984
985   for (t = 0; t < cmd->n_tables; ++t)
986     {
987       const struct mtable *table = &cmd->table[t];
988
989       output_case_processing_summary (table);
990
991       for (i = 0; i < table->n_interactions; ++i)
992         {
993           output_report (cmd, i, table);
994         }
995
996       categoricals_destroy (table->cats);
997     }
998 }
999
1000
1001 static void
1002 output_case_processing_summary (const struct mtable *table)
1003 {
1004   int i, v;
1005   const int heading_columns = 1;
1006   const int heading_rows = 3;
1007   struct tab_table *t;
1008
1009   const int nr = heading_rows + table->n_interactions * table->n_dep_vars;
1010   const int nc = 7;
1011
1012   t = tab_create (nc, nr);
1013   tab_title (t, _("Case Processing Summary"));
1014
1015   tab_headers (t, heading_columns, 0, heading_rows, 0);
1016
1017   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1018
1019   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1020   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1021
1022
1023   tab_joint_text (t, heading_columns, 0,
1024                   nc - 1, 0, TAB_CENTER | TAT_TITLE, _("Cases"));
1025
1026   tab_joint_text (t, 1, 1, 2, 1, TAB_CENTER | TAT_TITLE, _("Included"));
1027   tab_joint_text (t, 3, 1, 4, 1, TAB_CENTER | TAT_TITLE, _("Excluded"));
1028   tab_joint_text (t, 5, 1, 6, 1, TAB_CENTER | TAT_TITLE, _("Total"));
1029
1030   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1031   tab_hline (t, TAL_1, heading_columns, nc - 1, 2);
1032
1033
1034   for (i = 0; i < 3; ++i)
1035     {
1036       tab_text (t, heading_columns + i * 2, 2, TAB_CENTER | TAT_TITLE,
1037                 _("N"));
1038       tab_text (t, heading_columns + i * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1039                 _("Percent"));
1040     }
1041
1042   for (v = 0; v < table->n_dep_vars; ++v)
1043     {
1044       const struct variable *var = table->dep_vars[v];
1045       const char *dv_name = var_to_string (var);
1046       for (i = 0; i < table->n_interactions; ++i)
1047         {
1048           const int row = v * table->n_interactions + i;
1049           const struct interaction *iact = table->interactions[i];
1050           casenumber n_total;
1051
1052           struct string str;
1053           ds_init_cstr (&str, dv_name);
1054           ds_put_cstr (&str, ": ");
1055
1056           interaction_to_string (iact, &str);
1057
1058           tab_text (t, 0, row + heading_rows,
1059                     TAB_LEFT | TAT_TITLE, ds_cstr (&str));
1060
1061
1062           n_total = table->summary[row].missing + 
1063             table->summary[row].non_missing;
1064
1065           tab_double (t, 1, row + heading_rows,
1066                       0, table->summary[row].non_missing, &F_8_0);
1067
1068           tab_text_format (t, 2, row + heading_rows,
1069                            0, _("%g%%"), 
1070                            table->summary[row].non_missing / (double) n_total * 100.0);
1071
1072
1073           tab_double (t, 3, row + heading_rows,
1074                       0, table->summary[row].missing, &F_8_0);
1075
1076
1077           tab_text_format (t, 4, row + heading_rows,
1078                            0, _("%g%%"), 
1079                            table->summary[row].missing / (double) n_total * 100.0);
1080
1081
1082           tab_double (t, 5, row + heading_rows,
1083                       0, table->summary[row].missing + 
1084                       table->summary[row].non_missing, &F_8_0);
1085
1086           tab_text_format (t, 6, row + heading_rows,
1087                            0, _("%g%%"), 
1088                            n_total / (double) n_total * 100.0);
1089
1090
1091           ds_destroy (&str);
1092         }
1093     }
1094
1095   tab_submit (t);
1096 }
1097
1098
1099
1100 static void
1101 output_report (const struct means *cmd,  int iact_idx,
1102                const struct mtable *table)
1103 {
1104   int grp;
1105   int i;
1106
1107   const struct interaction *iact = table->interactions[iact_idx];
1108
1109   const int heading_columns = 1 + iact->n_vars;
1110   const int heading_rows = 1;
1111   struct tab_table *t;
1112
1113   const int n_cats = categoricals_n_count (table->cats, iact_idx);
1114
1115   const int nr = n_cats * table->n_dep_vars + heading_rows;
1116
1117   const int nc = heading_columns + cmd->n_cells;
1118
1119   t = tab_create (nc, nr);
1120   tab_title (t, _("Report"));
1121
1122   tab_headers (t, heading_columns, 0, heading_rows, 0);
1123
1124   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1125
1126   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1127   tab_vline (t, TAL_2, iact->n_vars, 0, nr - 1);
1128
1129   for (i = 0; i < iact->n_vars; ++i)
1130     {
1131       tab_text (t, 1 + i, 0, TAB_CENTER | TAT_TITLE,
1132                 var_to_string (iact->vars[i]));
1133     }
1134
1135   for (i = 0; i < cmd->n_cells; ++i)
1136     {
1137       tab_text (t, heading_columns + i, 0,
1138                 TAB_CENTER | TAT_TITLE,
1139                 gettext (cell_spec[cmd->cells[i]].title));
1140     }
1141
1142
1143   for (i = 0; i < n_cats; ++i)
1144     {
1145       int v, dv;
1146       const struct ccase *c =
1147         categoricals_get_case_by_category_real (table->cats, iact_idx, i);
1148
1149       for (dv = 0; dv < table->n_dep_vars; ++dv)
1150         {
1151           tab_text (t, 0,
1152                     heading_rows + dv * n_cats,
1153                     TAB_RIGHT | TAT_TITLE,
1154                     var_get_name (table->dep_vars[dv])
1155                     );
1156
1157           if ( dv > 0)
1158             tab_hline (t, TAL_1, 0, nc - 1, heading_rows + dv * n_cats);
1159
1160           for (v = 0; v < iact->n_vars; ++v)
1161             {
1162               const struct variable *var = iact->vars[v];
1163               const union value *val = case_data (c, var);
1164               struct string str;
1165               ds_init_empty (&str);
1166               var_append_value_name (var, val, &str);
1167
1168               tab_text (t, 1 + v, heading_rows + dv * n_cats + i,
1169                         TAB_RIGHT | TAT_TITLE, ds_cstr (&str));
1170
1171               ds_destroy (&str);
1172             }
1173         }
1174     }
1175
1176   for (grp = 0; grp < n_cats; ++grp)
1177     {
1178       int dv;
1179       struct per_cat_data *per_cat_data =
1180         categoricals_get_user_data_by_category_real (table->cats, iact_idx, grp);
1181
1182       for (dv = 0; dv < table->n_dep_vars; ++dv)
1183         {
1184           const struct per_var_data *pvd = &per_cat_data->pvd[dv];
1185           for (i = 0; i < cmd->n_cells; ++i)
1186             {
1187               const int csi = cmd->cells[i];
1188               const struct cell_spec *cs = &cell_spec[csi];
1189
1190               double result = cs->sd (cmd, pvd, pvd->cell_stats[i]);
1191
1192               tab_double (t, heading_columns + i,
1193                           heading_rows + grp + dv * n_cats,
1194                           0, result, 0);
1195             }
1196         }
1197     }
1198
1199   tab_submit (t);
1200 }