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