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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2011, 2012 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/tab.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
 506   /* Dependent variable (s) */
 507   if (!parse_variables_const_pool (lexer, cmd->pool, cmd->dict,
 508                               &table->dep_vars, &table->n_dep_vars,
 509                               PV_NO_DUPLICATE | PV_NUMERIC))
 510     return false;
 511
 512   /* Factor variable (s) */
 513   while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 514     {
 515       if (lex_match (lexer, T_BY))
 516         {
 517           table->n_layers++;
 518           table->layers =
 519             pool_realloc (cmd->pool, table->layers,
 520                       sizeof (*table->layers) * table->n_layers);
 521
 522           if (!parse_variables_const_pool
 523               (lexer, cmd->pool, cmd->dict,
 524                &table->layers[table->n_layers - 1].factor_vars,
 525                &table->layers[table->n_layers - 1].n_factor_vars,
 526                PV_NO_DUPLICATE))
 527             return false;
 528
 529         }
 530     }
 531
 532   /* There is always at least one layer.
 533      However the final layer is the total, and not
 534      normally considered by the user as a
 535      layer.
 536   */
 537
 538   table->n_layers++;
 539   table->layers =
 540     pool_realloc (cmd->pool, table->layers,
 541               sizeof (*table->layers) * table->n_layers);
 542   table->layers[table->n_layers - 1].factor_vars = NULL;
 543   table->layers[table->n_layers - 1].n_factor_vars = 0;
 544
 545   return true;
 546 }
 547
 548 /* Match a variable.
 549    If the match succeeds, the variable will be placed in VAR.
 550    Returns true if successful */
 551 static bool
 552 lex_is_variable (struct lexer *lexer, const struct dictionary *dict,
 553                  int n)
 554 {
 555   const char *tstr;
 556   if (lex_next_token (lexer, n) !=  T_ID)
 557     return false;
 558
 559   tstr = lex_next_tokcstr (lexer, n);
 560
 561   if (NULL == dict_lookup_var (dict, tstr) )
 562     return false;
 563
 564   return true;
 565 }
 566
 567
 568 int
 569 cmd_means (struct lexer *lexer, struct dataset *ds)
 570 {
 571   int t;
 572   int i;
 573   int l;
 574   struct means means;
 575   bool more_tables = true;
 576
 577   means.pool = pool_create ();
 578
 579   means.exclude = MV_ANY;
 580   means.dep_exclude = MV_ANY;
 581   means.listwise_exclude = false;
 582   means.table = NULL;
 583   means.n_tables = 0;
 584
 585   means.dict = dataset_dict (ds);
 586
 587   means.n_cells = 3;
 588   means.cells = pool_calloc (means.pool, means.n_cells, sizeof (*means.cells));
 589
 590
 591   /* The first three items (MEAN, COUNT, STDDEV) are the default */
 592   for (i = 0; i < 3; ++i)
 593     means.cells[i] = i;
 594
 595
 596   /*   Optional TABLES =   */
 597   if (lex_match_id (lexer, "TABLES"))
 598     {
 599       lex_force_match (lexer, T_EQUALS);
 600     }
 601
 602
 603   more_tables = true;
 604   /* Parse the "tables" */
 605   while (more_tables)
 606     {
 607       means.n_tables ++;
 608       means.table = pool_realloc (means.pool, means.table, means.n_tables * sizeof (*means.table));
 609
 610       if (! parse_means_table_syntax (lexer, &means,
 611                                       &means.table[means.n_tables - 1]))
 612         {
 613           goto error;
 614         }
 615
 616       /* Look ahead to see if there are more tables to be parsed */
 617       more_tables = false;
 618       if ( T_SLASH == lex_next_token (lexer, 0) )
 619         {
 620           if (lex_is_variable (lexer, means.dict, 1) )
 621             {
 622               more_tables = true;
 623               lex_force_match (lexer, T_SLASH);
 624             }
 625         }
 626     }
 627
 628   /* /MISSING subcommand */
 629   while (lex_token (lexer) != T_ENDCMD)
 630     {
 631       lex_match (lexer, T_SLASH);
 632
 633       if (lex_match_id (lexer, "MISSING"))
 634         {
 635           /*
 636              If no MISSING subcommand is specified, each combination of
 637              a dependent variable and categorical variables is handled
 638              separately.
 639            */
 640           lex_match (lexer, T_EQUALS);
 641           if (lex_match_id (lexer, "INCLUDE"))
 642             {
 643               /*
 644                 Use the subcommand  "/MISSING=INCLUDE" to include user-missing
 645                 values in the analysis.
 646               */
 647
 648               means.exclude = MV_SYSTEM;
 649               means.dep_exclude = MV_SYSTEM;
 650             }
 651           else if (lex_match_id (lexer, "TABLE"))
 652             /*
 653               This is the default. (I think).
 654               Every case containing a complete set of variables for a given
 655               table. If any variable, categorical or dependent for in a table
 656               is missing (as defined by what?), then that variable will
 657               be dropped FOR THAT TABLE ONLY.
 658             */
 659             {
 660               means.listwise_exclude = true;
 661             }
 662           else if (lex_match_id (lexer, "DEPENDENT"))
 663             /*
 664              Use the command "/MISSING=DEPENDENT" to
 665              include user-missing values for the categorical variables,
 666              while excluding them for the dependent variables.
 667
 668              Cases are dropped only when user-missing values
 669              appear in dependent  variables.  User-missing
 670              values for categorical variables are treated according to
 671              their face value.
 672
 673              Cases are ALWAYS dropped when System Missing values appear
 674              in the categorical variables.
 675             */
 676             {
 677               means.dep_exclude = MV_ANY;
 678               means.exclude = MV_SYSTEM;
 679             }
 680           else
 681             {
 682               lex_error (lexer, NULL);
 683               goto error;
 684             }
 685         }
 686       else if (lex_match_id (lexer, "CELLS"))
 687         {
 688           lex_match (lexer, T_EQUALS);
 689
 690           /* The default values become overwritten */
 691           means.n_cells = 0;
 692           while (lex_token (lexer) != T_ENDCMD
 693                  && lex_token (lexer) != T_SLASH)
 694             {
 695               int k = 0;
 696               if (lex_match (lexer, T_ALL))
 697                 {
 698                   int x;
 699                   means.cells =
 700                     pool_realloc (means.pool, means.cells,
 701                               (means.n_cells += n_C) * sizeof (*means.cells));
 702
 703                   for (x = 0; x < n_C; ++x)
 704                     means.cells[means.n_cells - (n_C - 1 - x) - 1] = x;
 705                 }
 706               else if (lex_match_id (lexer, "NONE"))
 707                 {
 708                   /* Do nothing */
 709                 }
 710               else if (lex_match_id (lexer, "DEFAULT"))
 711                 {
 712                   means.cells =
 713                     pool_realloc (means.pool, means.cells,
 714                               (means.n_cells += 3) * sizeof (*means.cells));
 715
 716                   means.cells[means.n_cells - 2 - 1] = MEANS_MEAN;
 717                   means.cells[means.n_cells - 1 - 1] = MEANS_N;
 718                   means.cells[means.n_cells - 0 - 1] = MEANS_STDDEV;
 719                 }
 720               else
 721                 {
 722                   for (; k < n_C; ++k)
 723                     {
 724                       if (lex_match_id (lexer, cell_spec[k].keyword))
 725                         {
 726                           means.cells =
 727                             pool_realloc (means.pool, means.cells,
 728                                       ++means.n_cells * sizeof (*means.cells));
 729
 730                           means.cells[means.n_cells - 1] = k;
 731                           break;
 732                         }
 733                     }
 734                 }
 735               if (k >= n_C)
 736                 {
 737                   lex_error (lexer, NULL);
 738                   goto error;
 739                 }
 740             }
 741         }
 742       else
 743         {
 744           lex_error (lexer, NULL);
 745           goto error;
 746         }
 747     }
 748
 749
 750
 751   for (t = 0; t < means.n_tables; ++t)
 752   {
 753     struct mtable *table = &means.table[t];
 754
 755     table->interactions =
 756       pool_calloc (means.pool, table->n_layers, sizeof (*table->interactions));
 757
 758     table->summary =
 759       pool_calloc (means.pool, table->n_dep_vars * table->n_layers, sizeof (*table->summary));
 760
 761     for (l = 0; l < table->n_layers; ++l)
 762       {
 763         int v;
 764         const struct layer *lyr = &table->layers[l];
 765         const int n_vars = lyr->n_factor_vars;
 766         table->interactions[l] = interaction_create (NULL);
 767         for (v = 0 ; v < n_vars ; ++v)
 768           {
 769             interaction_add_variable (table->interactions[l],
 770                                       lyr->factor_vars[v]);
 771           }
 772       }
 773   }
 774
 775   {
 776     struct casegrouper *grouper;
 777     struct casereader *group;
 778     bool ok;
 779
 780     grouper = casegrouper_create_splits (proc_open (ds), means.dict);
 781     while (casegrouper_get_next_group (grouper, &group))
 782       {
 783         run_means (&means, group, ds);
 784       }
 785     ok = casegrouper_destroy (grouper);
 786     ok = proc_commit (ds) && ok;
 787   }
 788
 789
 790   pool_destroy (means.pool);
 791   return CMD_SUCCESS;
 792
 793 error:
 794
 795   pool_destroy (means.pool);
 796   return CMD_FAILURE;
 797 }
 798
 799
 800 static bool
 801 is_missing (const struct means *cmd,
 802             const struct variable *dvar,
 803             const struct interaction *iact,
 804             const struct ccase *c)
 805 {
 806   if ( interaction_case_is_missing (iact, c, cmd->exclude) )
 807     return true;
 808
 809
 810   if (var_is_value_missing (dvar,
 811                             case_data (c, dvar),
 812                             cmd->dep_exclude))
 813     return true;
 814
 815   return false;
 816 }
 817
 818 static void output_case_processing_summary (const struct mtable *);
 819
 820 static void output_report (const struct means *, int, const struct mtable *);
 821
 822
 823 struct per_cat_data
 824 {
 825   struct per_var_data *pvd;
 826
 827   bool warn;
 828 };
 829
 830
 831 static void
 832 destroy_n (const void *aux1 UNUSED, void *aux2, void *user_data)
 833 {
 834   struct mtable *table = aux2;
 835   int v;
 836   struct per_cat_data *per_cat_data = user_data;
 837   struct per_var_data *pvd = per_cat_data->pvd;
 838
 839   for (v = 0; v < table->n_dep_vars; ++v)
 840     {
 841       struct per_var_data *pp = &pvd[v];
 842       moments1_destroy (pp->mom);
 843     }
 844 }
 845
 846 static void *
 847 create_n (const void *aux1, void *aux2)
 848 {
 849   int i, v;
 850   const struct means *means = aux1;
 851   struct mtable *table = aux2;
 852   struct per_cat_data *per_cat_data = pool_malloc (means->pool, sizeof *per_cat_data);
 853
 854   struct per_var_data *pvd = pool_calloc (means->pool, table->n_dep_vars, sizeof *pvd);
 855
 856   for (v = 0; v < table->n_dep_vars; ++v)
 857     {
 858       enum moment maxmom = MOMENT_KURTOSIS;
 859       struct per_var_data *pp = &pvd[v];
 860
 861       pp->cell_stats = pool_calloc (means->pool, means->n_cells, sizeof *pp->cell_stats);
 862
 863
 864       for (i = 0; i < means->n_cells; ++i)
 865         {
 866           int csi = means->cells[i];
 867           const struct cell_spec *cs = &cell_spec[csi];
 868           if (cs->sc)
 869             {
 870               pp->cell_stats[i] = cs->sc (means->pool);
 871             }
 872         }
 873       pp->mom = moments1_create (maxmom);
 874     }
 875
 876
 877   per_cat_data->pvd = pvd;
 878   per_cat_data->warn = true;
 879   return per_cat_data;
 880 }
 881
 882 static void
 883 update_n (const void *aux1, void *aux2, void *user_data, const struct ccase *c, double weight)
 884 {
 885   int i;
 886   int v = 0;
 887   const struct means *means = aux1;
 888   struct mtable *table = aux2;
 889   struct per_cat_data *per_cat_data = user_data;
 890
 891   for (v = 0; v < table->n_dep_vars; ++v)
 892     {
 893       struct per_var_data *pvd = &per_cat_data->pvd[v];
 894
 895       const double x = case_data (c, table->dep_vars[v])->f;
 896
 897       for (i = 0; i < table->n_layers; ++i)
 898         {
 899           if ( is_missing (means, table->dep_vars[v],
 900                            table->interactions[i], c))
 901             goto end;
 902         }
 903
 904       for (i = 0; i < means->n_cells; ++i)
 905         {
 906           const int csi = means->cells[i];
 907           const struct cell_spec *cs = &cell_spec[csi];
 908
 909
 910           if (cs->su)
 911             cs->su (pvd->cell_stats[i],
 912                     weight, x);
 913         }
 914
 915       moments1_add (pvd->mom, x, weight);
 916
 917     end:
 918       continue;
 919     }
 920 }
 921
 922 static void
 923 calculate_n (const void *aux1, void *aux2, void *user_data)
 924 {
 925   int i;
 926   int v = 0;
 927   struct per_cat_data *per_cat_data = user_data;
 928   const struct means *means = aux1;
 929   struct mtable *table = aux2;
 930
 931   for (v = 0; v < table->n_dep_vars; ++v)
 932     {
 933       struct per_var_data *pvd = &per_cat_data->pvd[v];
 934       for (i = 0; i < means->n_cells; ++i)
 935         {
 936           int csi = means->cells[i];
 937           const struct cell_spec *cs = &cell_spec[csi];
 938
 939           if (cs->su)
 940             cs->sd (pvd, pvd->cell_stats[i]);
 941         }
 942     }
 943 }
 944
 945 static void
 946 run_means (struct means *cmd, struct casereader *input,
 947            const struct dataset *ds UNUSED)
 948 {
 949   int t;
 950   const struct variable *wv = dict_get_weight (cmd->dict);
 951   struct ccase *c;
 952   struct casereader *reader;
 953
 954   struct payload payload;
 955   payload.create = create_n;
 956   payload.update = update_n;
 957   payload.calculate = calculate_n;
 958   payload.destroy = destroy_n;
 959
 960   for (t = 0; t < cmd->n_tables; ++t)
 961   {
 962     struct mtable *table = &cmd->table[t];
 963     table->cats
 964       = categoricals_create (table->interactions,
 965                              table->n_layers, wv, cmd->dep_exclude, cmd->exclude);
 966
 967     categoricals_set_payload (table->cats, &payload, cmd, table);
 968   }
 969
 970   for (reader = input;
 971        (c = casereader_read (reader)) != NULL; case_unref (c))
 972     {
 973       for (t = 0; t < cmd->n_tables; ++t)
 974         {
 975           bool something_missing = false;
 976           int  v;
 977           struct mtable *table = &cmd->table[t];
 978
 979           for (v = 0; v < table->n_dep_vars; ++v)
 980             {
 981               int i;
 982               for (i = 0; i < table->n_layers; ++i)
 983                 {
 984                   const bool missing =
 985                     is_missing (cmd, table->dep_vars[v],
 986                                 table->interactions[i], c);
 987                   if (missing)
 988                     {
 989                       something_missing = true;
 990                       table->summary[v * table->n_layers + i].missing++;
 991                     }
 992                   else
 993                     table->summary[v * table->n_layers  + i].non_missing++;
 994                 }
 995             }
 996           if ( something_missing && cmd->listwise_exclude)
 997             continue;
 998
 999           categoricals_update (table->cats, c);
1000         }
1001     }
1002   casereader_destroy (reader);
1003
1004   for (t = 0; t < cmd->n_tables; ++t)
1005     {
1006       struct mtable *table = &cmd->table[t];
1007
1008       categoricals_done (table->cats);
1009     }
1010
1011
1012   for (t = 0; t < cmd->n_tables; ++t)
1013     {
1014       int i;
1015       const struct mtable *table = &cmd->table[t];
1016
1017       output_case_processing_summary (table);
1018
1019       for (i = 0; i < table->n_layers; ++i)
1020         {
1021           output_report (cmd, i, table);
1022         }
1023       categoricals_destroy (table->cats);
1024     }
1025
1026 }
1027
1028
1029
1030 static void
1031 output_case_processing_summary (const struct mtable *table)
1032 {
1033   int i, v;
1034   const int heading_columns = 1;
1035   const int heading_rows = 3;
1036   struct tab_table *t;
1037
1038   const int nr = heading_rows + table->n_layers * table->n_dep_vars;
1039   const int nc = 7;
1040
1041   t = tab_create (nc, nr);
1042   tab_title (t, _("Case Processing Summary"));
1043
1044   tab_headers (t, heading_columns, 0, heading_rows, 0);
1045
1046   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1047
1048   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1049   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1050
1051
1052   tab_joint_text (t, heading_columns, 0,
1053                   nc - 1, 0, TAB_CENTER | TAT_TITLE, _("Cases"));
1054
1055   tab_joint_text (t, 1, 1, 2, 1, TAB_CENTER | TAT_TITLE, _("Included"));
1056   tab_joint_text (t, 3, 1, 4, 1, TAB_CENTER | TAT_TITLE, _("Excluded"));
1057   tab_joint_text (t, 5, 1, 6, 1, TAB_CENTER | TAT_TITLE, _("Total"));
1058
1059   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1060   tab_hline (t, TAL_1, heading_columns, nc - 1, 2);
1061
1062
1063   for (i = 0; i < 3; ++i)
1064     {
1065       tab_text (t, heading_columns + i * 2, 2, TAB_CENTER | TAT_TITLE,
1066                 _("N"));
1067       tab_text (t, heading_columns + i * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1068                 _("Percent"));
1069     }
1070
1071   for (v = 0; v < table->n_dep_vars; ++v)
1072     {
1073       const struct variable *var = table->dep_vars[v];
1074       const char *dv_name = var_to_string (var);
1075       for (i = 0; i < table->n_layers; ++i)
1076         {
1077           const int row = v * table->n_layers + i;
1078           const struct interaction *iact = table->interactions[i];
1079           casenumber n_total;
1080
1081           struct string str;
1082           ds_init_cstr (&str, dv_name);
1083           ds_put_cstr (&str, ": ");
1084
1085           interaction_to_string (iact, &str);
1086
1087           tab_text (t, 0, row + heading_rows,
1088                     TAB_LEFT | TAT_TITLE, ds_cstr (&str));
1089
1090
1091           n_total = table->summary[row].missing +
1092             table->summary[row].non_missing;
1093
1094           tab_double (t, 1, row + heading_rows,
1095                       0, table->summary[row].non_missing, &F_8_0);
1096
1097           tab_text_format (t, 2, row + heading_rows,
1098                            0, _("%g%%"),
1099                            table->summary[row].non_missing / (double) n_total * 100.0);
1100
1101
1102           tab_double (t, 3, row + heading_rows,
1103                       0, table->summary[row].missing, &F_8_0);
1104
1105
1106           tab_text_format (t, 4, row + heading_rows,
1107                            0, _("%g%%"),
1108                            table->summary[row].missing / (double) n_total * 100.0);
1109
1110
1111           tab_double (t, 5, row + heading_rows,
1112                       0, table->summary[row].missing +
1113                       table->summary[row].non_missing, &F_8_0);
1114
1115           tab_text_format (t, 6, row + heading_rows,
1116                            0, _("%g%%"),
1117                            n_total / (double) n_total * 100.0);
1118
1119
1120           ds_destroy (&str);
1121         }
1122     }
1123
1124   tab_submit (t);
1125 }
1126
1127
1128 static void
1129 output_report (const struct means *cmd,  int iact_idx,
1130                const struct mtable *table)
1131 {
1132   int grp;
1133   int i;
1134
1135   const struct interaction *iact = table->interactions[iact_idx];
1136
1137   const int heading_columns = 1 + iact->n_vars;
1138   const int heading_rows = 1;
1139   struct tab_table *t;
1140
1141   const int n_cats = categoricals_n_count (table->cats, iact_idx);
1142
1143   const int nr = n_cats * table->n_dep_vars + heading_rows;
1144
1145   const int nc = heading_columns + cmd->n_cells;
1146
1147   t = tab_create (nc, nr);
1148   tab_title (t, _("Report"));
1149
1150   tab_headers (t, heading_columns, 0, heading_rows, 0);
1151
1152   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1153
1154   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1155   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1156
1157   for (i = 0; i < iact->n_vars; ++i)
1158     {
1159       tab_text (t, 1 + i, 0, TAB_CENTER | TAT_TITLE,
1160                 var_to_string (iact->vars[i]));
1161     }
1162
1163   for (i = 0; i < cmd->n_cells; ++i)
1164     {
1165       tab_text (t, heading_columns + i, 0,
1166                 TAB_CENTER | TAT_TITLE,
1167                 gettext (cell_spec[cmd->cells[i]].title));
1168     }
1169
1170
1171   for (i = 0; i < n_cats; ++i)
1172     {
1173       int v, dv;
1174       const struct ccase *c =
1175         categoricals_get_case_by_category_real (table->cats, iact_idx, i);
1176
1177       for (dv = 0; dv < table->n_dep_vars; ++dv)
1178         {
1179           tab_text (t, 0,
1180                     heading_rows + dv * n_cats,
1181                     TAB_RIGHT | TAT_TITLE,
1182                     var_to_string (table->dep_vars[dv])
1183                     );
1184
1185           if ( dv > 0)
1186             tab_hline (t, TAL_1, 0, nc - 1, heading_rows + dv * n_cats);
1187
1188           for (v = 0; v < iact->n_vars; ++v)
1189             {
1190               const struct variable *var = iact->vars[v];
1191               const union value *val = case_data (c, var);
1192               struct string str;
1193               ds_init_empty (&str);
1194               var_append_value_name (var, val, &str);
1195
1196               tab_text (t, 1 + v, heading_rows + dv * n_cats + i,
1197                         TAB_RIGHT | TAT_TITLE, ds_cstr (&str));
1198
1199               ds_destroy (&str);
1200             }
1201         }
1202     }
1203
1204   for (grp = 0; grp < n_cats; ++grp)
1205     {
1206       int dv;
1207       struct per_cat_data *per_cat_data =
1208         categoricals_get_user_data_by_category_real (table->cats, iact_idx, grp);
1209
1210       for (dv = 0; dv < table->n_dep_vars; ++dv)
1211         {
1212           const struct per_var_data *pvd = &per_cat_data->pvd[dv];
1213           for (i = 0; i < cmd->n_cells; ++i)
1214             {
1215               const int csi = cmd->cells[i];
1216               const struct cell_spec *cs = &cell_spec[csi];
1217
1218               double result = cs->sd (pvd, pvd->cell_stats[i]);
1219
1220               tab_double (t, heading_columns + i,
1221                           heading_rows + grp + dv * n_cats,
1222                           0, result, 0);
1223             }
1224         }
1225     }
1226
1227   tab_submit (t);
1228 }