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           struct pivot_value *pv
 571             = pivot_value_new_number (sg (cell->stat[idx]));
 572           if (NULL == cell_spec[stat].rc)
 573             {
 574               const struct variable *dv = mt->dep_vars[v];
 575               pv->numeric.format = * var_get_print_format (dv);
 576             }
 577           pivot_table_put (pt, indexes, pt->n_dimensions, pv);
 578         }
 579     }
 580   free (indexes);
 581
 582   for (int i = 0; i < cell->n_children; ++i)
 583     {
 584       struct cell_container *container = cell->children + i;
 585       struct cell *child = NULL;
 586       BT_FOR_EACH (child, struct cell, bt_node, &container->bt)
 587         {
 588           populate_table (means, mt, ws, child, pt);
 589         }
 590     }
 591 }
 592
 593 static void
 594 create_table_structure (const struct mtable *mt, struct pivot_table *pt,
 595                         const struct workspace *ws)
 596 {
 597   int * lindexes = ws->control_idx;
 598   /* The inner layers are situated rightmost in the table.
 599      So this iteration is in reverse order.  */
 600   for (int l = mt->n_layers -1; l >=0 ; --l)
 601     {
 602       const struct layer *layer = mt->layers[l];
 603       const struct cell_container *instances = ws->instances + l;
 604       const struct variable *var = layer->factor_vars[lindexes[l]];
 605       struct pivot_dimension *dim_layer
 606         = pivot_dimension_create (pt, PIVOT_AXIS_ROW,
 607                                   var_to_string (var));
 608       dim_layer->root->show_label = true;
 609
 610       /* Place the values of the control variables as table headings.  */
 611       {
 612         struct instance *inst = NULL;
 613         BT_FOR_EACH (inst, struct instance, bt_node, &instances->bt)
 614           {
 615             struct substring space = SS_LITERAL_INITIALIZER ("\t ");
 616             struct string str;
 617             ds_init_empty (&str);
 618             var_append_value_name (var,
 619                                    &inst->value,
 620                                    &str);
 621
 622             ds_ltrim (&str, space);
 623
 624             pivot_category_create_leaf (dim_layer->root,
 625                                         pivot_value_new_text (ds_cstr (&str)));
 626
 627             ds_destroy (&str);
 628           }
 629       }
 630
 631       pivot_category_create_leaf (dim_layer->root,
 632                                   pivot_value_new_text ("Total"));
 633     }
 634 }
 635
 636 /* Initialise C_DES with a string describing the control variable
 637    relating to MT, LINDEXES.  */
 638 static void
 639 layers_to_string (const struct mtable *mt, const int *lindexes,
 640                   struct string *c_des)
 641 {
 642   for (int l = 0; l < mt->n_layers; ++l)
 643     {
 644       const struct layer *layer = mt->layers[l];
 645       const struct variable *ctrl_var = layer->factor_vars[lindexes[l]];
 646       if (l > 0)
 647         ds_put_cstr (c_des, " * ");
 648       ds_put_cstr (c_des, var_get_name (ctrl_var));
 649     }
 650 }
 651
 652 static void
 653 populate_case_processing_summary (struct pivot_category *pc,
 654                                   const struct mtable *mt,
 655                                   const int *lindexes)
 656 {
 657   struct string ds;
 658   ds_init_empty (&ds);
 659   int l = 0;
 660   for (l = 0; l < mt->n_layers; ++l)
 661     {
 662       const struct layer *layer = mt->layers[l];
 663       const struct variable *ctrl_var = layer->factor_vars[lindexes[l]];
 664       if (l > 0)
 665         ds_put_cstr (&ds, " * ");
 666       ds_put_cstr (&ds, var_get_name (ctrl_var));
 667     }
 668   for (int dv = 0; dv < mt->n_dep_vars; ++dv)
 669     {
 670       struct string dss;
 671       ds_init_empty (&dss);
 672       ds_put_cstr (&dss, var_get_name (mt->dep_vars[dv]));
 673       if (mt->n_layers > 0)
 674         {
 675           ds_put_cstr (&dss, " * ");
 676           ds_put_substring (&dss, ds.ss);
 677         }
 678       pivot_category_create_leaf (pc,
 679                                   pivot_value_new_text (ds_cstr (&dss)));
 680       ds_destroy (&dss);
 681     }
 682
 683   ds_destroy (&ds);
 684 }
 685
 686 /* Create the "Case Processing Summary" table.  */
 687 static void
 688 means_case_processing_summary (const struct mtable *mt)
 689 {
 690   struct pivot_table *pt = pivot_table_create (N_("Case Processing Summary"));
 691
 692   struct pivot_dimension *dim_cases =
 693     pivot_dimension_create (pt, PIVOT_AXIS_COLUMN, N_("Cases"));
 694   dim_cases->root->show_label = true;
 695
 696   struct pivot_category *cats[3];
 697   cats[0] = pivot_category_create_group (dim_cases->root,
 698                                          N_("Included"), NULL);
 699   cats[1] = pivot_category_create_group (dim_cases->root,
 700                                          N_("Excluded"), NULL);
 701   cats[2] = pivot_category_create_group (dim_cases->root,
 702                                          N_("Total"), NULL);
 703   for (int i = 0; i < 3; ++i)
 704     {
 705       pivot_category_create_leaf_rc (cats[i],
 706                                      pivot_value_new_text (N_("N")),
 707                                      PIVOT_RC_COUNT);
 708       pivot_category_create_leaf_rc (cats[i],
 709                                      pivot_value_new_text (N_("Percent")),
 710                                      PIVOT_RC_PERCENT);
 711     }
 712
 713   struct pivot_dimension *rows =
 714     pivot_dimension_create (pt, PIVOT_AXIS_ROW, N_("Variables"));
 715
 716   for (int cmb = 0; cmb < mt->n_combinations; ++cmb)
 717     {
 718       const struct workspace *ws = mt->ws + cmb;
 719       populate_case_processing_summary (rows->root, mt, ws->control_idx);
 720       for (int dv = 0; dv < mt->n_dep_vars; ++dv)
 721         {
 722           int idx = cmb * mt->n_dep_vars + dv;
 723           const struct summary *summ = mt->summ + idx;
 724           double n_included = summ->n_total - summ->n_missing;
 725           pivot_table_put2 (pt, 5, idx,
 726                             pivot_value_new_number (100.0 * summ->n_total / summ->n_total));
 727           pivot_table_put2 (pt, 4, idx,
 728                             pivot_value_new_number (summ->n_total));
 729
 730           pivot_table_put2 (pt, 3, idx,
 731                             pivot_value_new_number (100.0 * summ->n_missing / summ->n_total));
 732           pivot_table_put2 (pt, 2, idx,
 733                             pivot_value_new_number (summ->n_missing));
 734
 735           pivot_table_put2 (pt, 1, idx,
 736                             pivot_value_new_number (100.0 * n_included / summ->n_total));
 737           pivot_table_put2 (pt, 0, idx,
 738                             pivot_value_new_number (n_included));
 739         }
 740     }
 741
 742   pivot_table_submit (pt);
 743 }
 744
 745 static void
 746 means_shipout_single (const struct mtable *mt, const struct means *means,
 747                       const struct workspace *ws)
 748 {
 749   struct pivot_table *pt = pivot_table_create (N_("Report"));
 750   pt->omit_empty = true;
 751
 752   struct pivot_dimension *dim_cells =
 753     pivot_dimension_create (pt, PIVOT_AXIS_COLUMN, N_("Statistics"));
 754
 755   /* Set the statistics headings, eg "Mean", "Std. Dev" etc.  */
 756   for (int i = 0; i < means->n_statistics; ++i)
 757     {
 758       const struct cell_spec *cs = cell_spec + means->statistics[i];
 759       pivot_category_create_leaf_rc
 760         (dim_cells->root,
 761          pivot_value_new_text (gettext (cs->title)), cs->rc);
 762     }
 763
 764   create_table_structure (mt, pt, ws);
 765   populate_table (means, mt, ws, ws->root_cell, pt);
 766   pivot_table_submit (pt);
 767 }
 768
 769
 770 static void
 771 means_shipout_multivar (const struct mtable *mt, const struct means *means,
 772                         const struct workspace *ws)
 773 {
 774   struct string dss;
 775   ds_init_empty (&dss);
 776   for (int dv = 0; dv < mt->n_dep_vars; ++dv)
 777     {
 778       if (dv > 0)
 779         ds_put_cstr (&dss, " * ");
 780       ds_put_cstr (&dss, var_get_name (mt->dep_vars[dv]));
 781     }
 782
 783   for (int l = 0; l < mt->n_layers; ++l)
 784     {
 785       ds_put_cstr (&dss, " * ");
 786       const struct layer *layer = mt->layers[l];
 787       const struct variable *var = layer->factor_vars[ws->control_idx[l]];
 788       ds_put_cstr (&dss, var_get_name (var));
 789     }
 790
 791   struct pivot_table *pt = pivot_table_create (ds_cstr (&dss));
 792   pt->omit_empty = true;
 793   ds_destroy (&dss);
 794
 795   struct pivot_dimension *dim_cells =
 796     pivot_dimension_create (pt, PIVOT_AXIS_COLUMN, N_("Variables"));
 797
 798   for (int i = 0; i < mt->n_dep_vars; ++i)
 799     {
 800       pivot_category_create_leaf
 801         (dim_cells->root,
 802          pivot_value_new_variable (mt->dep_vars[i]));
 803     }
 804
 805   struct pivot_dimension *dim_stats
 806     = pivot_dimension_create (pt, PIVOT_AXIS_ROW,
 807                               N_ ("Statistics"));
 808   dim_stats->root->show_label = false;
 809
 810   for (int i = 0; i < means->n_statistics; ++i)
 811     {
 812       const struct cell_spec *cs = cell_spec + means->statistics[i];
 813       pivot_category_create_leaf_rc
 814         (dim_stats->root,
 815          pivot_value_new_text (gettext (cs->title)), cs->rc);
 816     }
 817
 818   create_table_structure (mt, pt, ws);
 819   populate_table (means, mt, ws, ws->root_cell, pt);
 820   pivot_table_submit (pt);
 821 }
 822
 823 static void
 824 means_shipout (const struct mtable *mt, const struct means *means)
 825 {
 826   for (int cmb = 0; cmb < mt->n_combinations; ++cmb)
 827     {
 828       const struct workspace *ws = mt->ws + cmb;
 829       if (ws->root_cell == NULL)
 830         {
 831           struct string des;
 832           ds_init_empty (&des);
 833           layers_to_string (mt, ws->control_idx, &des);
 834           msg (MW, _("The table \"%s\" has no non-empty control variables."
 835                      "  No result for this table will be displayed."),
 836                ds_cstr (&des));
 837           ds_destroy (&des);
 838           continue;
 839         }
 840       if (mt->n_dep_vars > 1)
 841         means_shipout_multivar (mt, means, ws);
 842       else
 843         means_shipout_single (mt, means, ws);
 844     }
 845 }
 846
 847
 848 \f
 849
 850 static bool
 851 control_var_missing (const struct means *means,
 852                      const struct mtable *mt,
 853                      unsigned int not_wild UNUSED,
 854                      const struct ccase *c,
 855                      const struct workspace *ws)
 856 {
 857   bool miss = false;
 858   for (int l = 0; l < mt->n_layers; ++l)
 859     {
 860       /* if (0 == ((not_wild >> l) & 0x1)) */
 861       /* { */
 862       /*   continue; */
 863       /* } */
 864
 865       const struct layer *layer = mt->layers[l];
 866       const struct variable *var = layer->factor_vars[ws->control_idx[l]];
 867       const union value *vv = case_data (c, var);
 868
 869       miss = var_is_value_missing (var, vv, means->ctrl_exclude);
 870       if (miss)
 871         break;
 872     }
 873
 874   return miss;
 875 }
 876
 877 /* Lookup the set of control variables described by MT, C and NOT_WILD,
 878    in the hash table MAP.  If there is no such entry, then create a
 879    cell with these paremeters and add is to MAP.
 880    If the generated cell has childen, repeat for all the children.
 881    Returns the root cell.
 882 */
 883 static struct cell *
 884 service_cell_map (const struct means *means, const struct mtable *mt,
 885                  const struct ccase *c,
 886                  unsigned int not_wild,
 887                  struct hmap *map,
 888                  const struct cell *pcell,
 889                  int level,
 890                  const struct workspace *ws)
 891 {
 892   struct cell *cell = NULL;
 893   if (map)
 894     {
 895       if (!control_var_missing (means, mt, not_wild, c, ws))
 896         {
 897           /* Lookup this set of values in the cell's hash table.  */
 898           unsigned int hash = generate_hash (mt, c, not_wild, ws);
 899           cell = lookup_cell (mt, map, hash, c, not_wild, ws);
 900
 901           /* If it has not been seen before, then create a new
 902              subcell, with this set of values, and insert it
 903              into the table.  */
 904           if (cell == NULL)
 905             {
 906               cell = generate_cell (means, mt, c, not_wild, pcell, ws);
 907               hmap_insert (map, &cell->hmap_node, hash);
 908             }
 909         }
 910     }
 911   else
 912     {
 913       /* This condition should only happen in the root node case. */
 914       cell = ws->root_cell;
 915       if (cell == NULL &&
 916           !control_var_missing (means, mt, not_wild, c, ws))
 917         cell = generate_cell (means, mt, c, not_wild, pcell, ws);
 918     }
 919
 920   if (cell)
 921     {
 922       /* Here is where the business really happens!   After
 923          testing for missing values, the cell's statistics
 924          are accumulated.  */
 925       if (!control_var_missing (means, mt, not_wild, c, ws))
 926         {
 927           for (int v = 0; v < mt->n_dep_vars; ++v)
 928             {
 929               const struct variable *dep_var = mt->dep_vars[v];
 930               const union value *vv = case_data (c, dep_var);
 931               if (var_is_value_missing (dep_var, vv, means->dep_exclude))
 932                 continue;
 933
 934               for (int stat = 0; stat < means->n_statistics; ++stat)
 935                 {
 936                   const double weight = dict_get_case_weight (means->dict, c,
 937                                                               NULL);
 938                   stat_update *su = cell_spec[means->statistics[stat]].su;
 939                   su (cell->stat[stat + v * means->n_statistics], weight,
 940                       case_data (c, dep_var)->f);
 941                 }
 942             }
 943         }
 944
 945       /* Recurse into all the children (if there are any).  */
 946       for (int i = 0; i < cell->n_children; ++i)
 947         {
 948           struct cell_container *cc = cell->children + i;
 949           service_cell_map (means, mt, c,
 950                            not_wild | (0x1U << (i + level)),
 951                            &cc->map, cell, level + i + 1, ws);
 952         }
 953     }
 954
 955   return cell;
 956 }
 957
 958 /*  Do all the necessary preparation and pre-calculation that
 959     needs to be done before iterating the data.  */
 960 static void
 961 prepare_means (struct means *cmd)
 962 {
 963   for (int t = 0; t < cmd->n_tables; ++t)
 964     {
 965       struct mtable *mt = cmd->table + t;
 966
 967       for (int i = 0; i < mt->n_combinations; ++i)
 968         {
 969           struct workspace *ws = mt->ws + i;
 970           ws->root_cell = NULL;
 971           ws->control_idx = xzalloc (mt->n_layers
 972                                          * sizeof *ws->control_idx);
 973           ws->instances = xzalloc (mt->n_layers
 974                                          * 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 = xzalloc (mt->n_combinations * mt->n_dep_vars
1151                               * sizeof (*mt->summ));
1152           mt->ws = xzalloc (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 }