pintos-os.org Git - pspp/blob - src/language/stats/means.c

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