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 (lexer, 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             xrealloc (table->layers,
 520                       sizeof (*table->layers) * table->n_layers);
 521
 522           if (!parse_variables_const
 523               (lexer, 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     xrealloc (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.exclude = MV_ANY;
 578   means.dep_exclude = MV_ANY;
 579   means.listwise_exclude = false;
 580   means.table = NULL;
 581   means.n_tables = 0;
 582
 583   means.dict = dataset_dict (ds);
 584
 585   means.n_cells = 3;
 586   means.cells = xcalloc (means.n_cells, sizeof (*means.cells));
 587
 588
 589   /* The first three items (MEAN, COUNT, STDDEV) are the default */
 590   for (i = 0; i < 3; ++i)
 591     means.cells[i] = i;
 592
 593
 594   /*   Optional TABLES =   */
 595   if (lex_match_id (lexer, "TABLES"))
 596     {
 597       lex_force_match (lexer, T_EQUALS);
 598     }
 599
 600
 601   more_tables = true;
 602   /* Parse the "tables" */
 603   while (more_tables)
 604     {
 605       means.n_tables ++;
 606       means.table = xrealloc (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_force_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                     xrealloc (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                     xrealloc (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                             xrealloc (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   means.pool = pool_create ();
 748
 749
 750   for (t = 0; t < means.n_tables; ++t)
 751   {
 752     struct mtable *table = &means.table[t];
 753
 754     table->interactions =
 755       xcalloc (table->n_layers, sizeof (*table->interactions));
 756
 757     table->summary =
 758       xcalloc (table->n_dep_vars * table->n_layers, sizeof (*table->summary));
 759
 760     for (l = 0; l < table->n_layers; ++l)
 761       {
 762         int v;
 763         const struct layer *lyr = &table->layers[l];
 764         const int n_vars = lyr->n_factor_vars;
 765         table->interactions[l] = interaction_create (NULL);
 766         for (v = 0 ; v < n_vars ; ++v)
 767           {
 768             interaction_add_variable (table->interactions[l],
 769                                       lyr->factor_vars[v]);
 770           }
 771       }
 772   }
 773
 774   {
 775     struct casegrouper *grouper;
 776     struct casereader *group;
 777     bool ok;
 778
 779     grouper = casegrouper_create_splits (proc_open (ds), means.dict);
 780     while (casegrouper_get_next_group (grouper, &group))
 781       {
 782         run_means (&means, group, ds);
 783       }
 784     ok = casegrouper_destroy (grouper);
 785     ok = proc_commit (ds) && ok;
 786   }
 787
 788
 789   return CMD_SUCCESS;
 790
 791 error:
 792
 793   return CMD_FAILURE;
 794 }
 795
 796
 797 static bool
 798 is_missing (const struct means *cmd,
 799             const struct variable *dvar,
 800             const struct interaction *iact,
 801             const struct ccase *c)
 802 {
 803   if ( interaction_case_is_missing (iact, c, cmd->exclude) )
 804     return true;
 805
 806
 807   if (var_is_value_missing (dvar,
 808                             case_data (c, dvar),
 809                             cmd->dep_exclude))
 810     return true;
 811
 812   return false;
 813 }
 814
 815 static void output_case_processing_summary (const struct mtable *);
 816
 817 static void output_report (const struct means *, int, const struct mtable *);
 818
 819
 820 struct per_cat_data
 821 {
 822   struct per_var_data *pvd;
 823
 824   bool warn;
 825 };
 826
 827 static void *
 828 create_n (const void *aux1, void *aux2)
 829 {
 830   int i, v;
 831   const struct means *means = aux1;
 832   struct mtable *table = aux2;
 833   struct per_cat_data *per_cat_data = xmalloc (sizeof *per_cat_data);
 834
 835   struct per_var_data *pvd = xcalloc (table->n_dep_vars, sizeof *pvd);
 836
 837   for (v = 0; v < table->n_dep_vars; ++v)
 838     {
 839       enum moment maxmom = MOMENT_KURTOSIS;
 840       struct per_var_data *pp = &pvd[v];
 841
 842       pp->cell_stats = xcalloc (means->n_cells, sizeof *pp->cell_stats);
 843
 844
 845       for (i = 0; i < means->n_cells; ++i)
 846         {
 847           int csi = means->cells[i];
 848           const struct cell_spec *cs = &cell_spec[csi];
 849           if (cs->sc)
 850             {
 851               pp->cell_stats[i] = cs->sc (means->pool);
 852             }
 853         }
 854       pp->mom = moments1_create (maxmom);
 855     }
 856
 857
 858   per_cat_data->pvd = pvd;
 859   per_cat_data->warn = true;
 860   return per_cat_data;
 861 }
 862
 863 static void
 864 update_n (const void *aux1, void *aux2, void *user_data, const struct ccase *c, double weight)
 865 {
 866   int i;
 867   int v = 0;
 868   const struct means *means = aux1;
 869   struct mtable *table = aux2;
 870   struct per_cat_data *per_cat_data = user_data;
 871
 872   for (v = 0; v < table->n_dep_vars; ++v)
 873     {
 874       struct per_var_data *pvd = &per_cat_data->pvd[v];
 875
 876       const double x = case_data (c, table->dep_vars[v])->f;
 877
 878       for (i = 0; i < table->n_layers; ++i)
 879         {
 880           if ( is_missing (means, table->dep_vars[v],
 881                            table->interactions[i], c))
 882             goto end;
 883         }
 884
 885       for (i = 0; i < means->n_cells; ++i)
 886         {
 887           const int csi = means->cells[i];
 888           const struct cell_spec *cs = &cell_spec[csi];
 889
 890
 891           if (cs->su)
 892             cs->su (pvd->cell_stats[i],
 893                     weight, x);
 894         }
 895
 896       moments1_add (pvd->mom, x, weight);
 897
 898     end:
 899       continue;
 900     }
 901 }
 902
 903 static void
 904 calculate_n (const void *aux1, void *aux2, void *user_data)
 905 {
 906   int i;
 907   int v = 0;
 908   struct per_cat_data *per_cat_data = user_data;
 909   const struct means *means = aux1;
 910   struct mtable *table = aux2;
 911
 912   for (v = 0; v < table->n_dep_vars; ++v)
 913     {
 914       struct per_var_data *pvd = &per_cat_data->pvd[v];
 915       for (i = 0; i < means->n_cells; ++i)
 916         {
 917           int csi = means->cells[i];
 918           const struct cell_spec *cs = &cell_spec[csi];
 919
 920           if (cs->su)
 921             cs->sd (pvd, pvd->cell_stats[i]);
 922         }
 923     }
 924 }
 925
 926 static void
 927 run_means (struct means *cmd, struct casereader *input,
 928            const struct dataset *ds UNUSED)
 929 {
 930   int t;
 931   const struct variable *wv = dict_get_weight (cmd->dict);
 932   struct ccase *c;
 933   struct casereader *reader;
 934
 935   struct payload payload;
 936   payload.create = create_n;
 937   payload.update = update_n;
 938   payload.destroy = calculate_n;
 939
 940   for (t = 0; t < cmd->n_tables; ++t)
 941   {
 942     struct mtable *table = &cmd->table[t];
 943     table->cats
 944       = categoricals_create (table->interactions,
 945                              table->n_layers, wv, cmd->exclude);
 946
 947     categoricals_set_payload (table->cats, &payload, cmd, table);
 948   }
 949
 950   for (reader = casereader_clone (input);
 951        (c = casereader_read (reader)) != NULL; case_unref (c))
 952     {
 953       for (t = 0; t < cmd->n_tables; ++t)
 954         {
 955           bool something_missing = false;
 956           int  v;
 957           struct mtable *table = &cmd->table[t];
 958
 959           for (v = 0; v < table->n_dep_vars; ++v)
 960             {
 961               int i;
 962               for (i = 0; i < table->n_layers; ++i)
 963                 {
 964                   const bool missing =
 965                     is_missing (cmd, table->dep_vars[v],
 966                                 table->interactions[i], c);
 967                   if (missing)
 968                     {
 969                       something_missing = true;
 970                       table->summary[v * table->n_layers + i].missing++;
 971                     }
 972                   else
 973                     table->summary[v * table->n_layers  + i].non_missing++;
 974                 }
 975             }
 976           if ( something_missing && cmd->listwise_exclude)
 977             continue;
 978
 979           categoricals_update (table->cats, c);
 980         }
 981     }
 982   casereader_destroy (reader);
 983
 984   for (t = 0; t < cmd->n_tables; ++t)
 985     {
 986       struct mtable *table = &cmd->table[t];
 987
 988       categoricals_done (table->cats);
 989     }
 990
 991
 992   for (t = 0; t < cmd->n_tables; ++t)
 993     {
 994       int i;
 995       const struct mtable *table = &cmd->table[t];
 996
 997       output_case_processing_summary (table);
 998
 999       for (i = 0; i < table->n_layers; ++i)
1000         {
1001           output_report (cmd, i, table);
1002         }
1003       categoricals_destroy (table->cats);
1004     }
1005
1006 }
1007
1008
1009
1010 static void
1011 output_case_processing_summary (const struct mtable *table)
1012 {
1013   int i, v;
1014   const int heading_columns = 1;
1015   const int heading_rows = 3;
1016   struct tab_table *t;
1017
1018   const int nr = heading_rows + table->n_layers * table->n_dep_vars;
1019   const int nc = 7;
1020
1021   t = tab_create (nc, nr);
1022   tab_title (t, _("Case Processing Summary"));
1023
1024   tab_headers (t, heading_columns, 0, heading_rows, 0);
1025
1026   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1027
1028   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1029   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1030
1031
1032   tab_joint_text (t, heading_columns, 0,
1033                   nc - 1, 0, TAB_CENTER | TAT_TITLE, _("Cases"));
1034
1035   tab_joint_text (t, 1, 1, 2, 1, TAB_CENTER | TAT_TITLE, _("Included"));
1036   tab_joint_text (t, 3, 1, 4, 1, TAB_CENTER | TAT_TITLE, _("Excluded"));
1037   tab_joint_text (t, 5, 1, 6, 1, TAB_CENTER | TAT_TITLE, _("Total"));
1038
1039   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1040   tab_hline (t, TAL_1, heading_columns, nc - 1, 2);
1041
1042
1043   for (i = 0; i < 3; ++i)
1044     {
1045       tab_text (t, heading_columns + i * 2, 2, TAB_CENTER | TAT_TITLE,
1046                 _("N"));
1047       tab_text (t, heading_columns + i * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1048                 _("Percent"));
1049     }
1050
1051   for (v = 0; v < table->n_dep_vars; ++v)
1052     {
1053       const struct variable *var = table->dep_vars[v];
1054       const char *dv_name = var_to_string (var);
1055       for (i = 0; i < table->n_layers; ++i)
1056         {
1057           const int row = v * table->n_layers + i;
1058           const struct interaction *iact = table->interactions[i];
1059           casenumber n_total;
1060
1061           struct string str;
1062           ds_init_cstr (&str, dv_name);
1063           ds_put_cstr (&str, ": ");
1064
1065           interaction_to_string (iact, &str);
1066
1067           tab_text (t, 0, row + heading_rows,
1068                     TAB_LEFT | TAT_TITLE, ds_cstr (&str));
1069
1070
1071           n_total = table->summary[row].missing +
1072             table->summary[row].non_missing;
1073
1074           tab_double (t, 1, row + heading_rows,
1075                       0, table->summary[row].non_missing, &F_8_0);
1076
1077           tab_text_format (t, 2, row + heading_rows,
1078                            0, _("%g%%"),
1079                            table->summary[row].non_missing / (double) n_total * 100.0);
1080
1081
1082           tab_double (t, 3, row + heading_rows,
1083                       0, table->summary[row].missing, &F_8_0);
1084
1085
1086           tab_text_format (t, 4, row + heading_rows,
1087                            0, _("%g%%"),
1088                            table->summary[row].missing / (double) n_total * 100.0);
1089
1090
1091           tab_double (t, 5, row + heading_rows,
1092                       0, table->summary[row].missing +
1093                       table->summary[row].non_missing, &F_8_0);
1094
1095           tab_text_format (t, 6, row + heading_rows,
1096                            0, _("%g%%"),
1097                            n_total / (double) n_total * 100.0);
1098
1099
1100           ds_destroy (&str);
1101         }
1102     }
1103
1104   tab_submit (t);
1105 }
1106
1107
1108 static void
1109 output_report (const struct means *cmd,  int iact_idx,
1110                const struct mtable *table)
1111 {
1112   int grp;
1113   int i;
1114
1115   const struct interaction *iact = table->interactions[iact_idx];
1116
1117   const int heading_columns = 1 + iact->n_vars;
1118   const int heading_rows = 1;
1119   struct tab_table *t;
1120
1121   const int n_cats = categoricals_n_count (table->cats, iact_idx);
1122
1123   const int nr = n_cats * table->n_dep_vars + heading_rows;
1124
1125   const int nc = heading_columns + cmd->n_cells;
1126
1127   t = tab_create (nc, nr);
1128   tab_title (t, _("Report"));
1129
1130   tab_headers (t, heading_columns, 0, heading_rows, 0);
1131
1132   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1133
1134   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1135   tab_vline (t, TAL_2, iact->n_vars, 0, nr - 1);
1136
1137   for (i = 0; i < iact->n_vars; ++i)
1138     {
1139       tab_text (t, 1 + i, 0, TAB_CENTER | TAT_TITLE,
1140                 var_to_string (iact->vars[i]));
1141     }
1142
1143   for (i = 0; i < cmd->n_cells; ++i)
1144     {
1145       tab_text (t, heading_columns + i, 0,
1146                 TAB_CENTER | TAT_TITLE,
1147                 gettext (cell_spec[cmd->cells[i]].title));
1148     }
1149
1150
1151   for (i = 0; i < n_cats; ++i)
1152     {
1153       int v, dv;
1154       const struct ccase *c =
1155         categoricals_get_case_by_category_real (table->cats, iact_idx, i);
1156
1157       for (dv = 0; dv < table->n_dep_vars; ++dv)
1158         {
1159           tab_text (t, 0,
1160                     heading_rows + dv * n_cats,
1161                     TAB_RIGHT | TAT_TITLE,
1162                     var_get_name (table->dep_vars[dv])
1163                     );
1164
1165           if ( dv > 0)
1166             tab_hline (t, TAL_1, 0, nc - 1, heading_rows + dv * n_cats);
1167
1168           for (v = 0; v < iact->n_vars; ++v)
1169             {
1170               const struct variable *var = iact->vars[v];
1171               const union value *val = case_data (c, var);
1172               struct string str;
1173               ds_init_empty (&str);
1174               var_append_value_name (var, val, &str);
1175
1176               tab_text (t, 1 + v, heading_rows + dv * n_cats + i,
1177                         TAB_RIGHT | TAT_TITLE, ds_cstr (&str));
1178
1179               ds_destroy (&str);
1180             }
1181         }
1182     }
1183
1184   for (grp = 0; grp < n_cats; ++grp)
1185     {
1186       int dv;
1187       struct per_cat_data *per_cat_data =
1188         categoricals_get_user_data_by_category_real (table->cats, iact_idx, grp);
1189
1190       for (dv = 0; dv < table->n_dep_vars; ++dv)
1191         {
1192           const struct per_var_data *pvd = &per_cat_data->pvd[dv];
1193           for (i = 0; i < cmd->n_cells; ++i)
1194             {
1195               const int csi = cmd->cells[i];
1196               const struct cell_spec *cs = &cell_spec[csi];
1197
1198               double result = cs->sd (pvd, pvd->cell_stats[i]);
1199
1200               tab_double (t, heading_columns + i,
1201                           heading_rows + grp + dv * n_cats,
1202                           0, result, 0);
1203             }
1204         }
1205     }
1206
1207   tab_submit (t);
1208 }