pintos-os.org Git - pspp/blob - 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 }