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