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 /* Table of cell_specs */
 399 static const struct cell_spec cell_spec[] = {
 400   {N_("Mean"), "MEAN", NULL, NULL, arithmean_get},
 401   {N_("N"), "COUNT", NULL, NULL, n_get},
 402   {N_("Std. Deviation"), "STDDEV", NULL, NULL, stddev_get},
 403 #if 0
 404   {N_("Median"), "MEDIAN", NULL, NULL, NULL},
 405   {N_("Group Median"), "GMEDIAN", NULL, NULL, NULL},
 406 #endif
 407   {N_("S.E. Mean"), "SEMEAN", NULL, NULL, semean_get},
 408   {N_("Sum"), "SUM", NULL, NULL, sum_get},
 409   {N_("Min"), "MIN", min_create, min_update, min_get},
 410   {N_("Max"), "MAX", max_create, max_update, max_get},
 411   {N_("Range"), "RANGE", range_create, range_update, range_get},
 412   {N_("Variance"), "VARIANCE", NULL, NULL, variance_get},
 413   {N_("Kurtosis"), "KURT", NULL, NULL, kurt_get},
 414   {N_("S.E. Kurt"), "SEKURT", NULL, NULL, sekurt_get},
 415   {N_("Skewness"), "SKEW", NULL, NULL, skew_get},
 416   {N_("S.E. Skew"), "SESKEW", NULL, NULL, seskew_get},
 417   {N_("First"), "FIRST", first_create, first_update, first_get},
 418   {N_("Last"), "LAST", last_create, last_update, last_get},
 419 #if 0
 420   {N_("Percent N"), "NPCT", NULL, NULL, NULL},
 421   {N_("Percent Sum"), "SPCT", NULL, NULL, NULL},
 422 #endif
 423   {N_("Harmonic Mean"), "HARMONIC", harmonic_create, harmonic_update, harmonic_get},
 424   {N_("Geom. Mean"), "GEOMETRIC", geometric_create, geometric_update, geometric_get}
 425 };
 426
 427 #define n_C (sizeof (cell_spec) / sizeof (struct cell_spec))
 428
 429
 430 struct summary
 431 {
 432   casenumber missing;
 433   casenumber non_missing;
 434 };
 435
 436
 437 /* The thing parsed after TABLES= */
 438 struct mtable
 439 {
 440   size_t n_dep_vars;
 441   const struct variable **dep_vars;
 442
 443   size_t n_interactions;
 444   struct interaction **interactions;
 445   struct summary *summary;
 446
 447   size_t *n_factor_vars;
 448   const struct variable ***factor_vars;
 449
 450   int ii;
 451
 452   int n_layers;
 453
 454   struct categoricals *cats;
 455 };
 456
 457 struct means
 458 {
 459   const struct dictionary *dict;
 460
 461   struct mtable *table;
 462   size_t n_tables;
 463
 464   /* Missing value class for categorical variables */
 465   enum mv_class exclude;
 466
 467   /* Missing value class for dependent variables */
 468   enum mv_class dep_exclude;
 469
 470   /* an array  indicating which statistics are to be calculated */
 471   int *cells;
 472
 473   /* Size of cells */
 474   int n_cells;
 475
 476   /* Pool on which cell functions may allocate data */
 477   struct pool *pool;
 478 };
 479
 480
 481 static void
 482 run_means (struct means *cmd, struct casereader *input,
 483            const struct dataset *ds);
 484
 485 /* Append all the variables belonging to layer and all subsequent layers
 486    to iact. And then append iact to the means->interaction.
 487    This is a recursive function.
 488  */
 489 static void
 490 iact_append_factor (struct mtable *means, int layer,
 491                     const struct interaction *iact)
 492 {
 493   int v;
 494   const struct variable **fv;
 495
 496   if (layer >= means->n_layers)
 497     return;
 498
 499   fv = means->factor_vars[layer];
 500
 501   for (v = 0; v < means->n_factor_vars[layer]; ++v)
 502     {
 503       struct interaction *nexti = interaction_clone (iact);
 504
 505       interaction_add_variable (nexti, fv[v]);
 506
 507       iact_append_factor (means, layer + 1, nexti);
 508
 509       if (layer == means->n_layers - 1)
 510         {
 511           means->interactions[means->ii++] = nexti;
 512         }
 513     }
 514 }
 515
 516 static bool
 517 parse_means_table_syntax (struct lexer *lexer, const struct means *cmd, struct mtable *table)
 518 {
 519   table->ii = 0;
 520   table->n_layers = 0;
 521   table->factor_vars = NULL;
 522   table->n_factor_vars = NULL;
 523
 524   /* Dependent variable (s) */
 525   if (!parse_variables_const (lexer, cmd->dict,
 526                               &table->dep_vars, &table->n_dep_vars,
 527                               PV_NO_DUPLICATE | PV_NUMERIC))
 528     return false;
 529
 530   /* Factor variable (s) */
 531   while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 532     {
 533       if (lex_match (lexer, T_BY))
 534         {
 535           table->n_layers++;
 536           table->factor_vars =
 537             xrealloc (table->factor_vars,
 538                       sizeof (*table->factor_vars) * table->n_layers);
 539
 540           table->n_factor_vars =
 541             xrealloc (table->n_factor_vars,
 542                       sizeof (*table->n_factor_vars) * table->n_layers);
 543
 544           if (!parse_variables_const (lexer, cmd->dict,
 545                                       &table->factor_vars[table->n_layers - 1],
 546                                       &table->n_factor_vars[table->n_layers -
 547                                                             1],
 548                                       PV_NO_DUPLICATE))
 549             return false;
 550
 551         }
 552     }
 553
 554   return true;
 555 }
 556
 557 /* Match a variable.
 558    If the match succeeds, the variable will be placed in VAR.
 559    Returns true if successful */
 560 static bool
 561 lex_is_variable (struct lexer *lexer, const struct dictionary *dict,
 562                  int n)
 563 {
 564   const char *tstr;
 565   if (lex_next_token (lexer, n) !=  T_ID)
 566     return false;
 567
 568   tstr = lex_next_tokcstr (lexer, n);
 569
 570   if (NULL == dict_lookup_var (dict, tstr) )
 571     return false;
 572
 573   return true;
 574 }
 575
 576
 577 int
 578 cmd_means (struct lexer *lexer, struct dataset *ds)
 579 {
 580   int t;
 581   int i;
 582   int l;
 583   struct means means;
 584   bool more_tables = true;
 585
 586   means.exclude = MV_ANY;
 587   means.dep_exclude = MV_ANY;
 588   means.table = NULL;
 589   means.n_tables = 0;
 590
 591   means.dict = dataset_dict (ds);
 592
 593   means.n_cells = 3;
 594   means.cells = xcalloc (means.n_cells, sizeof (*means.cells));
 595
 596
 597   /* The first three items (MEAN, COUNT, STDDEV) are the default */
 598   for (i = 0; i < 3; ++i)
 599     means.cells[i] = i;
 600
 601
 602   /*   Optional TABLES =   */
 603   if (lex_match_id (lexer, "TABLES"))
 604     {
 605       lex_force_match (lexer, T_EQUALS);
 606     }
 607
 608
 609   more_tables = true;
 610   /* Parse the "tables" */
 611   while (more_tables)
 612     {
 613       means.n_tables ++;
 614       means.table = xrealloc (means.table, means.n_tables * sizeof (*means.table));
 615
 616       if (! parse_means_table_syntax (lexer, &means,
 617                                       &means.table[means.n_tables - 1]))
 618         {
 619           goto error;
 620         }
 621
 622       /* Look ahead to see if there are more tables to be parsed */
 623       more_tables = false;
 624       if ( T_SLASH == lex_next_token (lexer, 0) )
 625         {
 626           if (lex_is_variable (lexer, means.dict, 1) )
 627             {
 628               more_tables = true;
 629               lex_force_match (lexer, T_SLASH);
 630             }
 631         }
 632     }
 633
 634   /* /MISSING subcommand */
 635   while (lex_token (lexer) != T_ENDCMD)
 636     {
 637       lex_match (lexer, T_SLASH);
 638
 639       if (lex_match_id (lexer, "MISSING"))
 640         {
 641           /*
 642              If no MISSING subcommand is specified, each combination of
 643              a dependent variable and categorical variables is handled
 644              separately.
 645            */
 646           lex_match (lexer, T_EQUALS);
 647           if (lex_match_id (lexer, "INCLUDE"))
 648             {
 649               /*
 650                 Use the subcommand  "/MISSING=INCLUDE" to include user-missing
 651                 values in the analysis.
 652               */
 653
 654               means.exclude = MV_SYSTEM;
 655               means.dep_exclude = MV_SYSTEM;
 656             }
 657           else if (lex_match_id (lexer, "TABLE"))
 658             /*
 659               This is the default. (I think).
 660               Every case containing a complete set of variables for a given
 661               table. If any variable, categorical or dependent for in a table
 662               is missing (as defined by what?), then that variable will
 663               be dropped FOR THAT TABLE ONLY.
 664             */
 665             {
 666               means.exclude = MV_ANY;
 667               means.dep_exclude = MV_ANY;
 668             }
 669           else if (lex_match_id (lexer, "DEPENDENT"))
 670             /*
 671              Use the command "/MISSING=DEPENDENT" to
 672              include user-missing values for the categorical variables,
 673              while excluding them for the dependent variables.
 674
 675              Cases are dropped only when user-missing values
 676              appear in dependent  variables.  User-missing
 677              values for categorical variables are treated according to
 678              their face value.
 679
 680              Cases are ALWAYS dropped when System Missing values appear
 681              in the categorical variables.
 682             */
 683             {
 684               means.dep_exclude = MV_ANY;
 685               means.exclude = MV_SYSTEM;
 686             }
 687           else
 688             {
 689               lex_error (lexer, NULL);
 690               goto error;
 691             }
 692         }
 693       else if (lex_match_id (lexer, "CELLS"))
 694         {
 695           lex_match (lexer, T_EQUALS);
 696
 697           /* The default values become overwritten */
 698           means.n_cells = 0;
 699           while (lex_token (lexer) != T_ENDCMD
 700                  && lex_token (lexer) != T_SLASH)
 701             {
 702               int k;
 703               for (k = 0; k < n_C; ++k)
 704                 {
 705                   if (lex_match_id (lexer, cell_spec[k].keyword))
 706                     {
 707                       means.cells =
 708                         xrealloc (means.cells,
 709                                   ++means.n_cells * sizeof (*means.cells));
 710
 711                       means.cells[means.n_cells - 1] = k;
 712                       break;
 713                     }
 714                 }
 715               if (k >= n_C)
 716                 {
 717                   lex_error (lexer, NULL);
 718                   goto error;
 719                 }
 720             }
 721         }
 722       else
 723         {
 724           lex_error (lexer, NULL);
 725           goto error;
 726         }
 727     }
 728
 729   means.pool = pool_create ();
 730
 731
 732   for (t = 0; t < means.n_tables; ++t)
 733   {
 734     struct mtable *table = &means.table[t];
 735     table->n_interactions = 1;
 736     for (l = 0; l < table->n_layers; ++l)
 737       {
 738         const int n_vars = table->n_factor_vars[l];
 739         table->n_interactions *= n_vars;
 740       }
 741
 742     table->interactions =
 743       xcalloc (table->n_interactions, sizeof (*table->interactions));
 744
 745     table->summary =
 746       xcalloc (table->n_dep_vars * table->n_interactions, sizeof (*table->summary));
 747
 748
 749     if (table->n_layers > 0)
 750       iact_append_factor (table, 0, interaction_create (NULL));
 751     else
 752       table->interactions[0] = interaction_create (NULL);
 753
 754   }
 755
 756
 757   {
 758     struct casegrouper *grouper;
 759     struct casereader *group;
 760     bool ok;
 761
 762     grouper = casegrouper_create_splits (proc_open (ds), means.dict);
 763     while (casegrouper_get_next_group (grouper, &group))
 764       {
 765         run_means (&means, group, ds);
 766       }
 767     ok = casegrouper_destroy (grouper);
 768     ok = proc_commit (ds) && ok;
 769   }
 770
 771
 772   return CMD_SUCCESS;
 773
 774 error:
 775
 776   return CMD_FAILURE;
 777 }
 778
 779
 780 static bool
 781 is_missing (const struct means *cmd,
 782             const struct variable *dvar,
 783             const struct interaction *iact,
 784             const struct ccase *c)
 785 {
 786   if ( interaction_case_is_missing (iact, c, cmd->exclude) )
 787     return true;
 788
 789
 790   if (var_is_value_missing (dvar,
 791                             case_data (c, dvar),
 792                             cmd->dep_exclude))
 793     return true;
 794
 795   return false;
 796 }
 797
 798 static void output_case_processing_summary (const struct mtable *);
 799
 800 static void output_report (const struct means *, int, const struct mtable *);
 801
 802
 803 struct per_cat_data
 804 {
 805   struct per_var_data *pvd;
 806
 807   bool warn;
 808 };
 809
 810 static void *
 811 create_n (const void *aux1, void *aux2)
 812 {
 813   int i, v;
 814   const struct means *means = aux1;
 815   struct mtable *table = aux2;
 816   struct per_cat_data *per_cat_data = xmalloc (sizeof *per_cat_data);
 817
 818   struct per_var_data *pvd = xcalloc (table->n_dep_vars, sizeof *pvd);
 819
 820   for (v = 0; v < table->n_dep_vars; ++v)
 821     {
 822       enum moment maxmom = MOMENT_KURTOSIS;
 823       struct per_var_data *pp = &pvd[v];
 824
 825       pp->cell_stats = xcalloc (means->n_cells, sizeof *pp->cell_stats);
 826
 827
 828       for (i = 0; i < means->n_cells; ++i)
 829         {
 830           int csi = means->cells[i];
 831           const struct cell_spec *cs = &cell_spec[csi];
 832           if (cs->sc)
 833             {
 834               pp->cell_stats[i] = cs->sc (means->pool);
 835             }
 836         }
 837       pp->mom = moments1_create (maxmom);
 838     }
 839
 840
 841   per_cat_data->pvd = pvd;
 842   per_cat_data->warn = true;
 843   return per_cat_data;
 844 }
 845
 846 static void
 847 update_n (const void *aux1, void *aux2, void *user_data, const struct ccase *c, double weight)
 848 {
 849   int i;
 850   int v = 0;
 851   const struct means *means = aux1;
 852   struct mtable *table = aux2;
 853   struct per_cat_data *per_cat_data = user_data;
 854
 855   for (v = 0; v < table->n_dep_vars; ++v)
 856     {
 857       struct per_var_data *pvd = &per_cat_data->pvd[v];
 858
 859       const double x = case_data (c, table->dep_vars[v])->f;
 860
 861       for (i = 0; i < table->n_interactions; ++i)
 862         {
 863           if ( is_missing (means, table->dep_vars[v], table->interactions[i], c))
 864             goto end;
 865         }
 866
 867       for (i = 0; i < means->n_cells; ++i)
 868         {
 869           const int csi = means->cells[i];
 870           const struct cell_spec *cs = &cell_spec[csi];
 871
 872
 873           if (cs->su)
 874             cs->su (pvd->cell_stats[i],
 875                     weight, x);
 876         }
 877
 878       moments1_add (pvd->mom, x, weight);
 879
 880     end:
 881       continue;
 882     }
 883 }
 884
 885 static void
 886 calculate_n (const void *aux1, void *aux2, void *user_data)
 887 {
 888   int i;
 889   int v = 0;
 890   struct per_cat_data *per_cat_data = user_data;
 891   const struct means *means = aux1;
 892   struct mtable *table = aux2;
 893
 894   for (v = 0; v < table->n_dep_vars; ++v)
 895     {
 896       struct per_var_data *pvd = &per_cat_data->pvd[v];
 897       for (i = 0; i < means->n_cells; ++i)
 898         {
 899           int csi = means->cells[i];
 900           const struct cell_spec *cs = &cell_spec[csi];
 901
 902           if (cs->su)
 903             cs->sd (pvd, pvd->cell_stats[i]);
 904         }
 905     }
 906 }
 907
 908
 909 static void
 910 run_means (struct means *cmd, struct casereader *input,
 911            const struct dataset *ds UNUSED)
 912 {
 913   int i,t;
 914   const struct variable *wv = dict_get_weight (cmd->dict);
 915   struct ccase *c;
 916   struct casereader *reader;
 917
 918   struct payload payload;
 919   payload.create = create_n;
 920   payload.update = update_n;
 921   payload.destroy = calculate_n;
 922
 923   for (t = 0; t < cmd->n_tables; ++t)
 924   {
 925     struct mtable *table = &cmd->table[t];
 926     table->cats
 927       = categoricals_create (table->interactions,
 928                              table->n_interactions, wv, cmd->exclude);
 929
 930     categoricals_set_payload (table->cats, &payload, cmd, table);
 931   }
 932
 933   for (reader = casereader_clone (input);
 934        (c = casereader_read (reader)) != NULL; case_unref (c))
 935     {
 936       for (t = 0; t < cmd->n_tables; ++t)
 937         {
 938           int  v;
 939           struct mtable *table = &cmd->table[t];
 940
 941           for (v = 0; v < table->n_dep_vars; ++v)
 942             {
 943               int i;
 944               for (i = 0; i < table->n_interactions; ++i)
 945                 {
 946                   const bool missing =
 947                     is_missing (cmd, table->dep_vars[v],
 948                                 table->interactions[i], c);
 949                   if (missing)
 950                     table->summary[v * table->n_interactions + i].missing++;
 951                   else
 952                     table->summary[v * table->n_interactions + i].non_missing++;
 953                 }
 954             }
 955           categoricals_update (table->cats, c);
 956         }
 957     }
 958   casereader_destroy (reader);
 959
 960   for (t = 0; t < cmd->n_tables; ++t)
 961     {
 962       struct mtable *table = &cmd->table[t];
 963
 964       categoricals_done (table->cats);
 965     }
 966
 967
 968   for (t = 0; t < cmd->n_tables; ++t)
 969     {
 970       const struct mtable *table = &cmd->table[t];
 971
 972       output_case_processing_summary (table);
 973
 974       for (i = 0; i < table->n_interactions; ++i)
 975         {
 976           output_report (cmd, i, table);
 977         }
 978
 979       categoricals_destroy (table->cats);
 980     }
 981 }
 982
 983
 984 static void
 985 output_case_processing_summary (const struct mtable *table)
 986 {
 987   int i, v;
 988   const int heading_columns = 1;
 989   const int heading_rows = 3;
 990   struct tab_table *t;
 991
 992   const int nr = heading_rows + table->n_interactions * table->n_dep_vars;
 993   const int nc = 7;
 994
 995   t = tab_create (nc, nr);
 996   tab_title (t, _("Case Processing Summary"));
 997
 998   tab_headers (t, heading_columns, 0, heading_rows, 0);
 999
1000   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1001
1002   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1003   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1004
1005
1006   tab_joint_text (t, heading_columns, 0,
1007                   nc - 1, 0, TAB_CENTER | TAT_TITLE, _("Cases"));
1008
1009   tab_joint_text (t, 1, 1, 2, 1, TAB_CENTER | TAT_TITLE, _("Included"));
1010   tab_joint_text (t, 3, 1, 4, 1, TAB_CENTER | TAT_TITLE, _("Excluded"));
1011   tab_joint_text (t, 5, 1, 6, 1, TAB_CENTER | TAT_TITLE, _("Total"));
1012
1013   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1014   tab_hline (t, TAL_1, heading_columns, nc - 1, 2);
1015
1016
1017   for (i = 0; i < 3; ++i)
1018     {
1019       tab_text (t, heading_columns + i * 2, 2, TAB_CENTER | TAT_TITLE,
1020                 _("N"));
1021       tab_text (t, heading_columns + i * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1022                 _("Percent"));
1023     }
1024
1025   for (v = 0; v < table->n_dep_vars; ++v)
1026     {
1027       const struct variable *var = table->dep_vars[v];
1028       const char *dv_name = var_to_string (var);
1029       for (i = 0; i < table->n_interactions; ++i)
1030         {
1031           const int row = v * table->n_interactions + i;
1032           const struct interaction *iact = table->interactions[i];
1033           casenumber n_total;
1034
1035           struct string str;
1036           ds_init_cstr (&str, dv_name);
1037           ds_put_cstr (&str, ": ");
1038
1039           interaction_to_string (iact, &str);
1040
1041           tab_text (t, 0, row + heading_rows,
1042                     TAB_LEFT | TAT_TITLE, ds_cstr (&str));
1043
1044
1045           n_total = table->summary[row].missing +
1046             table->summary[row].non_missing;
1047
1048           tab_double (t, 1, row + heading_rows,
1049                       0, table->summary[row].non_missing, &F_8_0);
1050
1051           tab_text_format (t, 2, row + heading_rows,
1052                            0, _("%g%%"),
1053                            table->summary[row].non_missing / (double) n_total * 100.0);
1054
1055
1056           tab_double (t, 3, row + heading_rows,
1057                       0, table->summary[row].missing, &F_8_0);
1058
1059
1060           tab_text_format (t, 4, row + heading_rows,
1061                            0, _("%g%%"),
1062                            table->summary[row].missing / (double) n_total * 100.0);
1063
1064
1065           tab_double (t, 5, row + heading_rows,
1066                       0, table->summary[row].missing +
1067                       table->summary[row].non_missing, &F_8_0);
1068
1069           tab_text_format (t, 6, row + heading_rows,
1070                            0, _("%g%%"),
1071                            n_total / (double) n_total * 100.0);
1072
1073
1074           ds_destroy (&str);
1075         }
1076     }
1077
1078   tab_submit (t);
1079 }
1080
1081
1082
1083 static void
1084 output_report (const struct means *cmd,  int iact_idx,
1085                const struct mtable *table)
1086 {
1087   int grp;
1088   int i;
1089
1090   const struct interaction *iact = table->interactions[iact_idx];
1091
1092   const int heading_columns = 1 + iact->n_vars;
1093   const int heading_rows = 1;
1094   struct tab_table *t;
1095
1096   const int n_cats = categoricals_n_count (table->cats, iact_idx);
1097
1098   const int nr = n_cats * table->n_dep_vars + heading_rows;
1099
1100   const int nc = heading_columns + cmd->n_cells;
1101
1102   t = tab_create (nc, nr);
1103   tab_title (t, _("Report"));
1104
1105   tab_headers (t, heading_columns, 0, heading_rows, 0);
1106
1107   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1108
1109   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1110   tab_vline (t, TAL_2, iact->n_vars, 0, nr - 1);
1111
1112   for (i = 0; i < iact->n_vars; ++i)
1113     {
1114       tab_text (t, 1 + i, 0, TAB_CENTER | TAT_TITLE,
1115                 var_to_string (iact->vars[i]));
1116     }
1117
1118   for (i = 0; i < cmd->n_cells; ++i)
1119     {
1120       tab_text (t, heading_columns + i, 0,
1121                 TAB_CENTER | TAT_TITLE,
1122                 gettext (cell_spec[cmd->cells[i]].title));
1123     }
1124
1125
1126   for (i = 0; i < n_cats; ++i)
1127     {
1128       int v, dv;
1129       const struct ccase *c =
1130         categoricals_get_case_by_category_real (table->cats, iact_idx, i);
1131
1132       for (dv = 0; dv < table->n_dep_vars; ++dv)
1133         {
1134           tab_text (t, 0,
1135                     heading_rows + dv * n_cats,
1136                     TAB_RIGHT | TAT_TITLE,
1137                     var_get_name (table->dep_vars[dv])
1138                     );
1139
1140           if ( dv > 0)
1141             tab_hline (t, TAL_1, 0, nc - 1, heading_rows + dv * n_cats);
1142
1143           for (v = 0; v < iact->n_vars; ++v)
1144             {
1145               const struct variable *var = iact->vars[v];
1146               const union value *val = case_data (c, var);
1147               struct string str;
1148               ds_init_empty (&str);
1149               var_append_value_name (var, val, &str);
1150
1151               tab_text (t, 1 + v, heading_rows + dv * n_cats + i,
1152                         TAB_RIGHT | TAT_TITLE, ds_cstr (&str));
1153
1154               ds_destroy (&str);
1155             }
1156         }
1157     }
1158
1159   for (grp = 0; grp < n_cats; ++grp)
1160     {
1161       int dv;
1162       struct per_cat_data *per_cat_data =
1163         categoricals_get_user_data_by_category_real (table->cats, iact_idx, grp);
1164
1165       for (dv = 0; dv < table->n_dep_vars; ++dv)
1166         {
1167           const struct per_var_data *pvd = &per_cat_data->pvd[dv];
1168           for (i = 0; i < cmd->n_cells; ++i)
1169             {
1170               const int csi = cmd->cells[i];
1171               const struct cell_spec *cs = &cell_spec[csi];
1172
1173               double result = cs->sd (pvd, pvd->cell_stats[i]);
1174
1175               tab_double (t, heading_columns + i,
1176                           heading_rows + grp + dv * n_cats,
1177                           0, result, 0);
1178             }
1179         }
1180     }
1181
1182   tab_submit (t);
1183 }