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