pintos-os.org Git - pspp/blob - src/language/stats/means.c

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