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