dictionary: Limit split file variables to 8, for compatibility.
[pspp] / src / language / stats / means.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2019 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 "libpspp/hmap.h"
28 #include "libpspp/bt.h"
29 #include "libpspp/hash-functions.h"
30 #include "libpspp/misc.h"
31 #include "libpspp/pool.h"
32
33 #include "language/command.h"
34
35 #include "count-one-bits.h"
36 #include "count-leading-zeros.h"
37
38 #include "output/pivot-table.h"
39
40 #include "means.h"
41
42
43 #include "gettext.h"
44 #define _(msgid) gettext (msgid)
45 #define N_(msgid) (msgid)
46
47
48 /* A "cell" in this procedure represents a distinct value of the
49    procedure's categorical variables,  and a set of summary statistics
50    of all cases which whose categorical variables have that set of
51    values.   For example,  the dataset
52
53    v1    v2    cat1     cat2
54    100   202      0     1
55    100   202      0     2
56    100   202      1     0
57    100   202      0     1
58
59
60    has three cells in layer 0 and two cells in layer 1  in addition
61    to a "grand summary" cell to which all (non-missing) cases
62    contribute.
63
64    The cells form a n-ary tree structure with the "grand summary"
65    cell at the root.
66 */
67 struct cell
68 {
69   struct hmap_node hmap_node; /* Element in hash table. */
70   struct bt_node  bt_node;    /* Element in binary tree */
71
72   int n_children;
73   struct cell_container *children;
74
75   /* The statistics to be calculated for the cell.  */
76   struct statistic **stat;
77
78   /* The parent of this cell, or NULL if this is the root cell.  */
79   const struct cell *parent_cell;
80
81   /* A bit-field variable which indicates which control variables
82      are allocated a fixed value (for this cell),  and which are
83      "wildcards".
84
85      A one indicates a fixed value.  A zero indicates a wildcard.
86      Wildcard values are used to calculate totals and sub-totals.
87   */
88   unsigned int not_wild;
89
90   /* The value(s). */
91   union value *values;
92
93   /* The variables corresponding to the above values.  */
94   const struct variable **vars;
95 };
96
97 /*  A structure used to find the union of all values used
98     within a layer, and to sort those values.  */
99 struct instance
100 {
101   struct hmap_node hmap_node; /* Element in hash table. */
102   struct bt_node  bt_node;    /* Element in binary tree */
103
104   /* A unique, consecutive, zero based index identifying this
105      instance.  */
106   int index;
107
108   /* The top level value of this instance.  */
109   union value value;
110   const struct variable *var;
111 };
112
113
114 static void
115 destroy_workspace (const struct mtable *mt, struct workspace *ws)
116 {
117   for (int l = 0; l < mt->n_layers; ++l)
118     {
119       struct cell_container *instances = ws->instances + l;
120       struct instance *inst;
121       struct instance *next;
122       HMAP_FOR_EACH_SAFE (inst, next, struct instance, hmap_node,
123                           &instances->map)
124         {
125           int width = var_get_width (inst->var);
126           value_destroy (&inst->value, width);
127           free (inst);
128         }
129       hmap_destroy (&instances->map);
130     }
131   free (ws->control_idx);
132   free (ws->instances);
133 }
134
135 /* Destroy CELL.  */
136 static void
137 destroy_cell (const struct means *means,
138               const struct mtable *mt, struct cell *cell)
139 {
140   int idx = 0;
141   for (int i = 0; i < mt->n_layers; ++i)
142     {
143       if (0 == ((cell->not_wild >> i) & 0x1))
144         continue;
145
146       const struct layer *layer = mt->layers[i];
147       for (int cmb = 0; cmb < mt->n_combinations; ++cmb)
148       {
149         struct workspace *ws = mt->ws + cmb;
150         const struct variable *var
151           = layer->factor_vars[ws->control_idx[i]];
152
153         int width = var_get_width (var);
154         value_destroy (&cell->values[idx++], width);
155       }
156     }
157   for (int i = 0; i < cell->n_children; ++i)
158     {
159       struct cell_container *container = cell->children + i;
160       hmap_destroy (&container->map);
161     }
162
163   for (int v = 0; v < mt->n_dep_vars; ++v)
164     {
165       for (int s = 0; s < means->n_statistics; ++s)
166         {
167           stat_destroy *des = cell_spec[means->statistics[s]].sf;
168           des (cell->stat[s + v * means->n_statistics]);
169         }
170     }
171   free (cell->stat);
172
173   free (cell->children);
174   free (cell->values);
175   free (cell->vars);
176   free (cell);
177 }
178
179
180 /* Walk the tree in postorder starting from CELL and destroy all the
181    cells.  */
182 static void
183 means_destroy_cells (const struct means *means, struct cell *cell,
184                      const struct mtable *table)
185 {
186   for (int i = 0; i < cell->n_children; ++i)
187     {
188       struct cell_container *container = cell->children + i;
189       struct cell *sub_cell;
190       struct cell *next;
191       HMAP_FOR_EACH_SAFE (sub_cell,  next, struct cell, hmap_node,
192                           &container->map)
193         {
194           means_destroy_cells (means, sub_cell, table);
195         }
196     }
197
198   destroy_cell (means, table, cell);
199 }
200
201 #if 0
202
203 static void
204 dump_cell (const struct cell *cell, const struct mtable *mt, int level)
205 {
206   for (int l = 0; l < level; ++l)
207     putchar (' ');
208   printf ("%p: ", cell);
209   for (int i = 0; i < mt->n_layers; ++i)
210     {
211       putchar (((cell->not_wild >> i) & 0x1) ? 'w' : '.');
212     }
213   printf (" - ");
214   int x = 0;
215   for (int i = 0; i < mt->n_layers; ++i)
216     {
217       if ((cell->not_wild >> i) & 0x1)
218         {
219           printf ("%s: ", var_get_name (cell->vars[x]));
220           printf ("%g ", cell->values[x++].f);
221         }
222       else
223         printf ("x ");
224     }
225   stat_get *sg = cell_spec[MEANS_N].sd;
226   printf ("--- S1: %g", sg (cell->stat[0]));
227
228   printf ("--- N Children: %d", cell->n_children);
229   //  printf ("--- Level: %d", level);
230   printf ("--- Parent: %p", cell->parent_cell);
231   printf ("\n");
232 }
233
234 static void
235 dump_indeces (const size_t *indexes, int n)
236 {
237   for (int i = 0 ; i < n; ++i)
238     {
239       printf ("%ld; ", indexes[i]);
240     }
241   printf ("\n");
242 }
243
244 /* Dump the tree in pre-order.  */
245 static void
246 dump_tree (const struct cell *cell, const struct mtable *table,
247            int level, const struct cell *parent)
248 {
249   assert (cell->parent_cell == parent);
250   dump_cell (cell, table, level);
251
252   for (int i = 0; i < cell->n_children; ++i)
253     {
254       struct cell_container *container = cell->children + i;
255       struct cell *sub_cell;
256       BT_FOR_EACH (sub_cell, struct cell, bt_node, &container->bt)
257         {
258           dump_tree (sub_cell, table, level + 1, cell);
259         }
260     }
261 }
262
263 #endif
264
265 /* Generate a hash based on the values of the N variables in
266    the array VARS which are taken from the case C.  */
267 static unsigned int
268 generate_hash (const struct mtable *mt,
269                const struct ccase *c,
270                unsigned int not_wild,
271                const struct workspace *ws)
272 {
273   unsigned int hash = 0;
274   for (int i = 0; i < mt->n_layers; ++i)
275     {
276       if (0 == ((not_wild >> i) & 0x1))
277         continue;
278
279       const struct layer *layer = mt->layers[i];
280       const struct variable *var = layer->factor_vars[ws->control_idx[i]];
281       const union value *vv = case_data (c, var);
282       int width = var_get_width (var);
283       hash = hash_int (i, hash);
284       hash = value_hash (vv, width, hash);
285     }
286
287   return hash;
288 }
289
290 /* Create a cell based on the N variables in the array VARS,
291    which are indeces into the case C.
292    The caller is responsible for destroying this cell when
293    no longer needed. */
294 static struct cell *
295 generate_cell (const struct means *means,
296                const struct mtable *mt,
297                const struct ccase *c,
298                unsigned int not_wild,
299                const struct cell *pcell,
300                const struct workspace *ws)
301 {
302   int n_vars = count_one_bits (not_wild);
303   struct cell *cell = XZALLOC (struct cell);
304   cell->values = xcalloc (n_vars, sizeof *cell->values);
305   cell->vars = xcalloc (n_vars, sizeof *cell->vars);
306   cell->not_wild = not_wild;
307
308   cell->parent_cell = pcell;
309   cell->n_children = mt->n_layers -
310     (sizeof (cell->not_wild) * CHAR_BIT) +
311     count_leading_zeros (cell->not_wild);
312
313   int idx = 0;
314   for (int i = 0; i < mt->n_layers; ++i)
315     {
316       if (0 == ((not_wild >> i) & 0x1))
317         continue;
318
319       const struct layer *layer = mt->layers[i];
320       const struct variable *var = layer->factor_vars[ws->control_idx[i]];
321       const union value *vv = case_data (c, var);
322       int width = var_get_width (var);
323       cell->vars[idx] = var;
324       value_clone (&cell->values[idx++], vv, width);
325     }
326   assert (idx == n_vars);
327
328   cell->children = xcalloc (cell->n_children, sizeof *cell->children);
329   for (int i = 0; i < cell->n_children; ++i)
330     {
331       struct cell_container *container = cell->children + i;
332       hmap_init (&container->map);
333     }
334
335   cell->stat = xcalloc (means->n_statistics * mt->n_dep_vars, sizeof *cell->stat);
336   for (int v = 0; v < mt->n_dep_vars; ++v)
337     {
338       for (int stat = 0; stat < means->n_statistics; ++stat)
339         {
340           stat_create *sc = cell_spec[means->statistics[stat]].sc;
341
342           cell->stat[stat + v * means->n_statistics] = sc (means->pool);
343         }
344     }
345   return cell;
346 }
347
348
349 /* If a  cell based on the N variables in the array VARS,
350    which are indeces into the case C and whose hash is HASH,
351    exists in HMAP, then return that cell.
352    Otherwise, return NULL.  */
353 static struct cell *
354 lookup_cell (const struct mtable *mt,
355              struct hmap *hmap,  unsigned int hash,
356              const struct ccase *c,
357              unsigned int not_wild,
358              const struct workspace *ws)
359 {
360   struct cell *cell = NULL;
361   HMAP_FOR_EACH_WITH_HASH (cell, struct cell, hmap_node, hash, hmap)
362     {
363       bool match = true;
364       int idx = 0;
365       if (cell->not_wild != not_wild)
366         continue;
367       for (int i = 0; i < mt->n_layers; ++i)
368         {
369           if (0 == ((cell->not_wild >> i) & 0x1))
370             continue;
371
372           const struct layer *layer = mt->layers[i];
373           const struct variable *var = layer->factor_vars[ws->control_idx[i]];
374           const union value *vv = case_data (c, var);
375           int width = var_get_width (var);
376           assert (var == cell->vars[idx]);
377           if (!value_equal (vv, &cell->values[idx++], width))
378             {
379               match = false;
380               break;
381             }
382         }
383       if (match)
384         return cell;
385     }
386   return NULL;
387 }
388
389
390 /*  A comparison function used to sort cells in a binary tree.
391     Only the innermost value needs to be compared, because no
392     two cells with similar outer values will appear in the same
393     tree/map.   */
394 static int
395 cell_compare_3way (const struct bt_node *a,
396                    const struct bt_node *b,
397                    const void *aux UNUSED)
398 {
399   const struct cell *fa = BT_DATA (a, struct cell, bt_node);
400   const struct cell *fb = BT_DATA (b, struct cell, bt_node);
401
402   assert (fa->not_wild == fb->not_wild);
403   int vidx = count_one_bits (fa->not_wild) - 1;
404   assert (fa->vars[vidx] == fb->vars[vidx]);
405
406   return value_compare_3way (&fa->values[vidx],
407                              &fb->values[vidx],
408                              var_get_width (fa->vars[vidx]));
409 }
410
411 /*  A comparison function used to sort cells in a binary tree.  */
412 static int
413 compare_instance_3way (const struct bt_node *a,
414                        const struct bt_node *b,
415                        const void *aux UNUSED)
416 {
417   const struct instance *fa = BT_DATA (a, struct instance, bt_node);
418   const struct instance *fb = BT_DATA (b, struct instance, bt_node);
419
420   assert (fa->var == fb->var);
421
422   return  value_compare_3way (&fa->value,
423                               &fb->value,
424                               var_get_width (fa->var));
425 }
426
427
428 static void arrange_cells (struct workspace *ws,
429                            struct cell *cell, const struct mtable *table);
430
431
432 /* Iterate CONTAINER's map inserting a copy of its elements into
433    CONTAINER's binary tree.    Also, for each layer in TABLE, create
434    an instance container, containing the union of all elements in
435    CONTAINER.  */
436 static void
437 arrange_cell (struct workspace *ws, struct cell_container *container,
438               const struct mtable *mt)
439 {
440   struct bt *bt = &container->bt;
441   struct hmap *map = &container->map;
442   bt_init (bt, cell_compare_3way, NULL);
443
444   struct cell *cell;
445   HMAP_FOR_EACH (cell, struct cell, hmap_node, map)
446     {
447       bt_insert (bt, &cell->bt_node);
448
449       int idx = 0;
450       for (int i = 0; i < mt->n_layers; ++i)
451         {
452           if (0 == ((cell->not_wild >> i) & 0x1))
453             continue;
454
455           struct cell_container *instances = ws->instances + i;
456           const struct variable *var = cell->vars[idx];
457           int width = var_get_width (var);
458           unsigned int hash
459             = value_hash (&cell->values[idx], width, 0);
460
461           struct instance *inst = NULL;
462           struct instance *next = NULL;
463           HMAP_FOR_EACH_WITH_HASH_SAFE (inst, next, struct instance,
464                                         hmap_node,
465                                         hash, &instances->map)
466             {
467               assert (cell->vars[idx] == var);
468               if (value_equal (&inst->value,
469                                &cell->values[idx],
470                                width))
471                 {
472                   break;
473                 }
474             }
475
476           if (!inst)
477             {
478               inst = xzalloc (sizeof *inst);
479               inst->index = -1;
480               inst->var = var;
481               value_clone (&inst->value, &cell->values[idx],
482                            width);
483               hmap_insert (&instances->map, &inst->hmap_node, hash);
484             }
485
486           idx++;
487         }
488
489       arrange_cells (ws, cell, mt);
490     }
491 }
492
493 /* Arrange the children and then all the subtotals.  */
494 static void
495 arrange_cells (struct workspace *ws, struct cell *cell,
496                const struct mtable *table)
497 {
498   for (int i = 0; i < cell->n_children; ++i)
499     {
500       struct cell_container *container = cell->children + i;
501       arrange_cell (ws, container, table);
502     }
503 }
504
505
506 \f
507
508 /*  If the top level value in CELL, has an instance in the L_IDX'th layer,
509     then return that instance.  Otherwise return NULL.  */
510 static const struct instance *
511 lookup_instance (const struct mtable *mt, const struct workspace *ws,
512                  int l_idx, const struct cell *cell)
513 {
514   const struct layer *layer = mt->layers[l_idx];
515   int n_vals = count_one_bits (cell->not_wild);
516   const struct variable *var = layer->factor_vars[ws->control_idx[l_idx]];
517   const union value *val = cell->values + n_vals - 1;
518   int width = var_get_width (var);
519   unsigned int hash = value_hash (val, width, 0);
520   const struct cell_container *instances = ws->instances + l_idx;
521   struct instance *inst = NULL;
522   struct instance *next;
523   HMAP_FOR_EACH_WITH_HASH_SAFE (inst, next,
524                                 struct instance, hmap_node,
525                                 hash, &instances->map)
526     {
527       if (value_equal (val, &inst->value, width))
528         break;
529     }
530   return inst;
531 }
532
533 /* Enter the values into PT.  */
534 static void
535 populate_table (const struct means *means, const struct mtable *mt,
536                 const struct workspace *ws,
537                 const struct cell *cell,
538                 struct pivot_table *pt)
539 {
540   size_t *indexes = XCALLOC (pt->n_dimensions, size_t);
541   for (int v = 0; v < mt->n_dep_vars; ++v)
542     {
543       for (int s = 0; s < means->n_statistics; ++s)
544         {
545           int i = 0;
546           if (mt->n_dep_vars > 1)
547             indexes[i++] = v;
548           indexes[i++] = s;
549           int stat = means->statistics[s];
550           stat_get *sg = cell_spec[stat].sd;
551           {
552             const struct cell *pc = cell;
553             for (; i < pt->n_dimensions; ++i)
554               {
555                 int l_idx = pt->n_dimensions - i - 1;
556                 const struct cell_container *instances = ws->instances + l_idx;
557                 if (0 == (cell->not_wild >> l_idx & 0x1U))
558                   {
559                     indexes [i] = hmap_count (&instances->map);
560                   }
561                 else
562                   {
563                     assert (pc);
564                     const struct instance *inst
565                       = lookup_instance (mt, ws, l_idx, pc);
566                     assert (inst);
567                     indexes [i] = inst->index;
568                     pc = pc->parent_cell;
569                   }
570               }
571           }
572
573           int idx = s + v * means->n_statistics;
574           struct pivot_value *pv
575             = pivot_value_new_number (sg (cell->stat[idx]));
576           if (NULL == cell_spec[stat].rc)
577             {
578               const struct variable *dv = mt->dep_vars[v];
579               pv->numeric.format = * var_get_print_format (dv);
580             }
581           pivot_table_put (pt, indexes, pt->n_dimensions, pv);
582         }
583     }
584   free (indexes);
585
586   for (int i = 0; i < cell->n_children; ++i)
587     {
588       struct cell_container *container = cell->children + i;
589       struct cell *child = NULL;
590       BT_FOR_EACH (child, struct cell, bt_node, &container->bt)
591         {
592           populate_table (means, mt, ws, child, pt);
593         }
594     }
595 }
596
597 static void
598 create_table_structure (const struct mtable *mt, struct pivot_table *pt,
599                         const struct workspace *ws)
600 {
601   int * lindexes = ws->control_idx;
602   /* The inner layers are situated rightmost in the table.
603      So this iteration is in reverse order.  */
604   for (int l = mt->n_layers -1; l >=0 ; --l)
605     {
606       const struct layer *layer = mt->layers[l];
607       const struct cell_container *instances = ws->instances + l;
608       const struct variable *var = layer->factor_vars[lindexes[l]];
609       struct pivot_dimension *dim_layer
610         = pivot_dimension_create (pt, PIVOT_AXIS_ROW,
611                                   var_to_string (var));
612       dim_layer->root->show_label = true;
613
614       /* Place the values of the control variables as table headings.  */
615       {
616         struct instance *inst = NULL;
617         BT_FOR_EACH (inst, struct instance, bt_node, &instances->bt)
618           {
619             struct substring space = SS_LITERAL_INITIALIZER ("\t ");
620             struct string str;
621             ds_init_empty (&str);
622             var_append_value_name (var,
623                                    &inst->value,
624                                    &str);
625
626             ds_ltrim (&str, space);
627
628             pivot_category_create_leaf (dim_layer->root,
629                                         pivot_value_new_text (ds_cstr (&str)));
630
631             ds_destroy (&str);
632           }
633       }
634
635       pivot_category_create_leaf (dim_layer->root,
636                                   pivot_value_new_text ("Total"));
637     }
638 }
639
640 /* Initialise C_DES with a string describing the control variable
641    relating to MT, LINDEXES.  */
642 static void
643 layers_to_string (const struct mtable *mt, const int *lindexes,
644                   struct string *c_des)
645 {
646   for (int l = 0; l < mt->n_layers; ++l)
647     {
648       const struct layer *layer = mt->layers[l];
649       const struct variable *ctrl_var = layer->factor_vars[lindexes[l]];
650       if (l > 0)
651         ds_put_cstr (c_des, " * ");
652       ds_put_cstr (c_des, var_get_name (ctrl_var));
653     }
654 }
655
656 static void
657 populate_case_processing_summary (struct pivot_category *pc,
658                                   const struct mtable *mt,
659                                   const int *lindexes)
660 {
661   struct string ds;
662   ds_init_empty (&ds);
663   int l = 0;
664   for (l = 0; l < mt->n_layers; ++l)
665     {
666       const struct layer *layer = mt->layers[l];
667       const struct variable *ctrl_var = layer->factor_vars[lindexes[l]];
668       if (l > 0)
669         ds_put_cstr (&ds, " * ");
670       ds_put_cstr (&ds, var_get_name (ctrl_var));
671     }
672   for (int dv = 0; dv < mt->n_dep_vars; ++dv)
673     {
674       struct string dss;
675       ds_init_empty (&dss);
676       ds_put_cstr (&dss, var_get_name (mt->dep_vars[dv]));
677       if (mt->n_layers > 0)
678         {
679           ds_put_cstr (&dss, " * ");
680           ds_put_substring (&dss, ds.ss);
681         }
682       pivot_category_create_leaf (pc,
683                                   pivot_value_new_text (ds_cstr (&dss)));
684       ds_destroy (&dss);
685     }
686
687   ds_destroy (&ds);
688 }
689
690 /* Create the "Case Processing Summary" table.  */
691 static void
692 means_case_processing_summary (const struct mtable *mt)
693 {
694   struct pivot_table *pt = pivot_table_create (N_("Case Processing Summary"));
695
696   struct pivot_dimension *dim_cases =
697     pivot_dimension_create (pt, PIVOT_AXIS_COLUMN, N_("Cases"));
698   dim_cases->root->show_label = true;
699
700   struct pivot_category *cats[3];
701   cats[0] = pivot_category_create_group (dim_cases->root,
702                                          N_("Included"), NULL);
703   cats[1] = pivot_category_create_group (dim_cases->root,
704                                          N_("Excluded"), NULL);
705   cats[2] = pivot_category_create_group (dim_cases->root,
706                                          N_("Total"), NULL);
707   for (int i = 0; i < 3; ++i)
708     {
709       pivot_category_create_leaf_rc (cats[i],
710                                      pivot_value_new_text (N_("N")),
711                                      PIVOT_RC_COUNT);
712       pivot_category_create_leaf_rc (cats[i],
713                                      pivot_value_new_text (N_("Percent")),
714                                      PIVOT_RC_PERCENT);
715     }
716
717   struct pivot_dimension *rows =
718     pivot_dimension_create (pt, PIVOT_AXIS_ROW, N_("Variables"));
719
720   for (int cmb = 0; cmb < mt->n_combinations; ++cmb)
721     {
722       const struct workspace *ws = mt->ws + cmb;
723       populate_case_processing_summary (rows->root, mt, ws->control_idx);
724       for (int dv = 0; dv < mt->n_dep_vars; ++dv)
725         {
726           int idx = cmb * mt->n_dep_vars + dv;
727           const struct summary *summ = mt->summ + idx;
728           double n_included = summ->n_total - summ->n_missing;
729           pivot_table_put2 (pt, 5, idx,
730                             pivot_value_new_number (100.0 * summ->n_total / summ->n_total));
731           pivot_table_put2 (pt, 4, idx,
732                             pivot_value_new_number (summ->n_total));
733
734           pivot_table_put2 (pt, 3, idx,
735                             pivot_value_new_number (100.0 * summ->n_missing / summ->n_total));
736           pivot_table_put2 (pt, 2, idx,
737                             pivot_value_new_number (summ->n_missing));
738
739           pivot_table_put2 (pt, 1, idx,
740                             pivot_value_new_number (100.0 * n_included / summ->n_total));
741           pivot_table_put2 (pt, 0, idx,
742                             pivot_value_new_number (n_included));
743         }
744     }
745
746   pivot_table_submit (pt);
747 }
748
749 static void
750 means_shipout_single (const struct mtable *mt, const struct means *means,
751                       const struct workspace *ws)
752 {
753   struct pivot_table *pt = pivot_table_create (N_("Report"));
754
755   struct pivot_dimension *dim_cells =
756     pivot_dimension_create (pt, PIVOT_AXIS_COLUMN, N_("Statistics"));
757
758   /* Set the statistics headings, eg "Mean", "Std. Dev" etc.  */
759   for (int i = 0; i < means->n_statistics; ++i)
760     {
761       const struct cell_spec *cs = cell_spec + means->statistics[i];
762       pivot_category_create_leaf_rc
763         (dim_cells->root,
764          pivot_value_new_text (gettext (cs->title)), cs->rc);
765     }
766
767   create_table_structure (mt, pt, ws);
768   populate_table (means, mt, ws, ws->root_cell, pt);
769   pivot_table_submit (pt);
770 }
771
772
773 static void
774 means_shipout_multivar (const struct mtable *mt, const struct means *means,
775                         const struct workspace *ws)
776 {
777   struct string dss;
778   ds_init_empty (&dss);
779   for (int dv = 0; dv < mt->n_dep_vars; ++dv)
780     {
781       if (dv > 0)
782         ds_put_cstr (&dss, " * ");
783       ds_put_cstr (&dss, var_get_name (mt->dep_vars[dv]));
784     }
785
786   for (int l = 0; l < mt->n_layers; ++l)
787     {
788       ds_put_cstr (&dss, " * ");
789       const struct layer *layer = mt->layers[l];
790       const struct variable *var = layer->factor_vars[ws->control_idx[l]];
791       ds_put_cstr (&dss, var_get_name (var));
792     }
793
794   struct pivot_table *pt = pivot_table_create (ds_cstr (&dss));
795   ds_destroy (&dss);
796
797   struct pivot_dimension *dim_cells =
798     pivot_dimension_create (pt, PIVOT_AXIS_COLUMN, N_("Variables"));
799
800   for (int i = 0; i < mt->n_dep_vars; ++i)
801     {
802       pivot_category_create_leaf
803         (dim_cells->root,
804          pivot_value_new_variable (mt->dep_vars[i]));
805     }
806
807   struct pivot_dimension *dim_stats
808     = pivot_dimension_create (pt, PIVOT_AXIS_ROW,
809                               N_ ("Statistics"));
810   dim_stats->root->show_label = false;
811
812   for (int i = 0; i < means->n_statistics; ++i)
813     {
814       const struct cell_spec *cs = cell_spec + means->statistics[i];
815       pivot_category_create_leaf_rc
816         (dim_stats->root,
817          pivot_value_new_text (gettext (cs->title)), cs->rc);
818     }
819
820   create_table_structure (mt, pt, ws);
821   populate_table (means, mt, ws, ws->root_cell, pt);
822   pivot_table_submit (pt);
823 }
824
825 static void
826 means_shipout (const struct mtable *mt, const struct means *means)
827 {
828   for (int cmb = 0; cmb < mt->n_combinations; ++cmb)
829     {
830       const struct workspace *ws = mt->ws + cmb;
831       if (ws->root_cell == NULL)
832         {
833           struct string des;
834           ds_init_empty (&des);
835           layers_to_string (mt, ws->control_idx, &des);
836           msg (MW, _("The table \"%s\" has no non-empty control variables."
837                      "  No result for this table will be displayed."),
838                ds_cstr (&des));
839           ds_destroy (&des);
840           continue;
841         }
842       if (mt->n_dep_vars > 1)
843         means_shipout_multivar (mt, means, ws);
844       else
845         means_shipout_single (mt, means, ws);
846     }
847 }
848
849
850 \f
851
852 static bool
853 control_var_missing (const struct means *means,
854                      const struct mtable *mt,
855                      unsigned int not_wild UNUSED,
856                      const struct ccase *c,
857                      const struct workspace *ws)
858 {
859   bool miss = false;
860   for (int l = 0; l < mt->n_layers; ++l)
861     {
862       /* if (0 == ((not_wild >> l) & 0x1)) */
863       /* { */
864       /*   continue; */
865       /* } */
866
867       const struct layer *layer = mt->layers[l];
868       const struct variable *var = layer->factor_vars[ws->control_idx[l]];
869       const union value *vv = case_data (c, var);
870
871       miss = (var_is_value_missing (var, vv) & means->ctrl_exclude) != 0;
872       if (miss)
873         break;
874     }
875
876   return miss;
877 }
878
879 /* Lookup the set of control variables described by MT, C and NOT_WILD,
880    in the hash table MAP.  If there is no such entry, then create a
881    cell with these paremeters and add is to MAP.
882    If the generated cell has childen, repeat for all the children.
883    Returns the root cell.
884 */
885 static struct cell *
886 service_cell_map (const struct means *means, const struct mtable *mt,
887                  const struct ccase *c,
888                  unsigned int not_wild,
889                  struct hmap *map,
890                  const struct cell *pcell,
891                  int level,
892                  const struct workspace *ws)
893 {
894   struct cell *cell = NULL;
895   if (map)
896     {
897       if (!control_var_missing (means, mt, not_wild, c, ws))
898         {
899           /* Lookup this set of values in the cell's hash table.  */
900           unsigned int hash = generate_hash (mt, c, not_wild, ws);
901           cell = lookup_cell (mt, map, hash, c, not_wild, ws);
902
903           /* If it has not been seen before, then create a new
904              subcell, with this set of values, and insert it
905              into the table.  */
906           if (cell == NULL)
907             {
908               cell = generate_cell (means, mt, c, not_wild, pcell, ws);
909               hmap_insert (map, &cell->hmap_node, hash);
910             }
911         }
912     }
913   else
914     {
915       /* This condition should only happen in the root node case. */
916       cell = ws->root_cell;
917       if (cell == NULL &&
918           !control_var_missing (means, mt, not_wild, c, ws))
919         cell = generate_cell (means, mt, c, not_wild, pcell, ws);
920     }
921
922   if (cell)
923     {
924       /* Here is where the business really happens!   After
925          testing for missing values, the cell's statistics
926          are accumulated.  */
927       if (!control_var_missing (means, mt, not_wild, c, ws))
928         {
929           for (int v = 0; v < mt->n_dep_vars; ++v)
930             {
931               const struct variable *dep_var = mt->dep_vars[v];
932               const union value *vv = case_data (c, dep_var);
933               if (var_is_value_missing (dep_var, vv) & means->dep_exclude)
934                 continue;
935
936               for (int stat = 0; stat < means->n_statistics; ++stat)
937                 {
938                   const double weight = dict_get_case_weight (means->dict, c,
939                                                               NULL);
940                   stat_update *su = cell_spec[means->statistics[stat]].su;
941                   su (cell->stat[stat + v * means->n_statistics], weight,
942                       case_num (c, dep_var));
943                 }
944             }
945         }
946
947       /* Recurse into all the children (if there are any).  */
948       for (int i = 0; i < cell->n_children; ++i)
949         {
950           struct cell_container *cc = cell->children + i;
951           service_cell_map (means, mt, c,
952                            not_wild | (0x1U << (i + level)),
953                            &cc->map, cell, level + i + 1, ws);
954         }
955     }
956
957   return cell;
958 }
959
960 /*  Do all the necessary preparation and pre-calculation that
961     needs to be done before iterating the data.  */
962 static void
963 prepare_means (struct means *cmd)
964 {
965   for (int t = 0; t < cmd->n_tables; ++t)
966     {
967       struct mtable *mt = cmd->table + t;
968
969       for (int i = 0; i < mt->n_combinations; ++i)
970         {
971           struct workspace *ws = mt->ws + i;
972           ws->root_cell = NULL;
973           ws->control_idx = xcalloc (mt->n_layers, sizeof *ws->control_idx);
974           ws->instances = xcalloc (mt->n_layers, sizeof *ws->instances);
975           int cmb = i;
976           for (int l = mt->n_layers - 1; l >= 0; --l)
977             {
978               struct cell_container *instances = ws->instances + l;
979               const struct layer *layer = mt->layers[l];
980               ws->control_idx[l] = cmb % layer->n_factor_vars;
981               cmb /= layer->n_factor_vars;
982               hmap_init (&instances->map);
983             }
984         }
985     }
986 }
987
988
989 /* Do all the necessary calculations that occur AFTER iterating
990    the data.  */
991 static void
992 post_means (struct means *cmd)
993 {
994   for (int t = 0; t < cmd->n_tables; ++t)
995     {
996       struct mtable *mt = cmd->table + t;
997       for (int cmb = 0; cmb < mt->n_combinations; ++cmb)
998         {
999           struct workspace *ws = mt->ws + cmb;
1000           if (ws->root_cell == NULL)
1001             continue;
1002           arrange_cells (ws, ws->root_cell, mt);
1003           /*  The root cell should have no parent.  */
1004           assert (ws->root_cell->parent_cell == 0);
1005
1006           for (int l = 0; l < mt->n_layers; ++l)
1007             {
1008               struct cell_container *instances = ws->instances + l;
1009               bt_init (&instances->bt, compare_instance_3way, NULL);
1010
1011               /* Iterate the instance hash table, and insert each instance
1012                  into the binary tree BT.  */
1013               struct instance *inst;
1014               HMAP_FOR_EACH (inst, struct instance, hmap_node,
1015                              &instances->map)
1016                 {
1017                   bt_insert (&instances->bt, &inst->bt_node);
1018                 }
1019
1020               /* Iterate the binary tree (in order) and assign the index
1021                  member accordingly.  */
1022               int index = 0;
1023               BT_FOR_EACH (inst, struct instance, bt_node, &instances->bt)
1024                 {
1025                   inst->index = index++;
1026                 }
1027             }
1028         }
1029     }
1030 }
1031
1032
1033 /* Update the summary information (the missings and the totals).  */
1034 static void
1035 update_summaries (const struct means *means, struct mtable *mt,
1036                   const struct ccase *c, double weight)
1037 {
1038   for (int dv = 0; dv < mt->n_dep_vars; ++dv)
1039     {
1040       for (int cmb = 0; cmb < mt->n_combinations; ++cmb)
1041         {
1042           struct workspace *ws = mt->ws + cmb;
1043           struct summary *summ = mt->summ
1044             + cmb * mt->n_dep_vars + dv;
1045
1046           summ->n_total += weight;
1047           const struct variable *var = mt->dep_vars[dv];
1048           const union value *vv = case_data (c, var);
1049           /* First check if the dependent variable is missing.  */
1050           if (var_is_value_missing (var, vv) & means->dep_exclude)
1051             summ->n_missing += weight;
1052           /* If the dep var is not missing, then check each
1053              control variable.  */
1054           else
1055             for (int l = 0; l < mt->n_layers; ++l)
1056               {
1057                 const struct layer *layer = mt->layers [l];
1058                 const struct variable *var
1059                   = layer->factor_vars[ws->control_idx[l]];
1060                 const union value *vv = case_data (c, var);
1061                 if (var_is_value_missing (var, vv) & means->ctrl_exclude)
1062                   {
1063                     summ->n_missing += weight;
1064                     break;
1065                   }
1066               }
1067         }
1068     }
1069 }
1070
1071
1072 void
1073 run_means (struct means *cmd, struct casereader *input,
1074            const struct dataset *ds UNUSED)
1075 {
1076   struct ccase *c = NULL;
1077   struct casereader *reader;
1078
1079   prepare_means (cmd);
1080
1081   for (reader = input;
1082        (c = casereader_read (reader)) != NULL; case_unref (c))
1083     {
1084       const double weight
1085         = dict_get_case_weight (cmd->dict, c, NULL);
1086       for (int t = 0; t < cmd->n_tables; ++t)
1087         {
1088           struct mtable *mt = cmd->table + t;
1089           update_summaries (cmd, mt, c, weight);
1090
1091           for (int cmb = 0; cmb < mt->n_combinations; ++cmb)
1092             {
1093               struct workspace *ws = mt->ws + cmb;
1094
1095               ws->root_cell = service_cell_map (cmd, mt, c,
1096                                                 0U, NULL, NULL, 0, ws);
1097             }
1098         }
1099     }
1100   casereader_destroy (reader);
1101
1102   post_means (cmd);
1103 }
1104
1105 struct lexer;
1106
1107 int
1108 cmd_means (struct lexer *lexer, struct dataset *ds)
1109 {
1110   struct means means;
1111   means.pool = pool_create ();
1112
1113   means.ctrl_exclude = MV_ANY;
1114   means.dep_exclude = MV_ANY;
1115   means.table = NULL;
1116   means.n_tables = 0;
1117
1118   means.dict = dataset_dict (ds);
1119
1120   means.n_statistics = 3;
1121   means.statistics = pool_calloc (means.pool, 3, sizeof *means.statistics);
1122   means.statistics[0] = MEANS_MEAN;
1123   means.statistics[1] = MEANS_N;
1124   means.statistics[2] = MEANS_STDDEV;
1125
1126   if (! means_parse (lexer, &means))
1127     goto error;
1128
1129   /* Calculate some constant data for each table.  */
1130   for (int t = 0; t < means.n_tables; ++t)
1131     {
1132       struct mtable *mt = means.table + t;
1133       mt->n_combinations = 1;
1134       for (int l = 0; l < mt->n_layers; ++l)
1135         mt->n_combinations *= mt->layers[l]->n_factor_vars;
1136     }
1137
1138   {
1139     struct casegrouper *grouper;
1140     struct casereader *group;
1141     bool ok;
1142
1143     grouper = casegrouper_create_splits (proc_open (ds), means.dict);
1144     while (casegrouper_get_next_group (grouper, &group))
1145       {
1146         /* Allocate the workspaces.  */
1147         for (int t = 0; t < means.n_tables; ++t)
1148         {
1149           struct mtable *mt = means.table + t;
1150           mt->summ = xcalloc (mt->n_combinations * mt->n_dep_vars,
1151                               sizeof (*mt->summ));
1152           mt->ws = xcalloc (mt->n_combinations, sizeof (*mt->ws));
1153         }
1154         run_means (&means, group, ds);
1155         for (int t = 0; t < means.n_tables; ++t)
1156           {
1157             const struct mtable *mt = means.table + t;
1158
1159             means_case_processing_summary (mt);
1160             means_shipout (mt, &means);
1161
1162             for (int i = 0; i < mt->n_combinations; ++i)
1163               {
1164                 struct workspace *ws = mt->ws + i;
1165                 if (ws->root_cell == NULL)
1166                   continue;
1167
1168                 means_destroy_cells (&means, ws->root_cell, mt);
1169               }
1170           }
1171
1172         /* Destroy the workspaces.  */
1173         for (int t = 0; t < means.n_tables; ++t)
1174           {
1175             struct mtable *mt = means.table + t;
1176             free (mt->summ);
1177             for (int i = 0; i < mt->n_combinations; ++i)
1178               {
1179                 struct workspace *ws = mt->ws + i;
1180                 destroy_workspace (mt, ws);
1181               }
1182             free (mt->ws);
1183           }
1184       }
1185     ok = casegrouper_destroy (grouper);
1186     ok = proc_commit (ds) && ok;
1187   }
1188
1189   pool_destroy (means.pool);
1190   return CMD_SUCCESS;
1191
1192  error:
1193
1194   pool_destroy (means.pool);
1195   return CMD_FAILURE;
1196 }