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