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

   1
   2 /* PSPP - a program for statistical analysis.
   3    Copyright (C) 2011, 2012 Free Software Foundation, Inc.
   4
   5    This program is free software: you can redistribute it and/or modify
   6    it under the terms of the GNU General Public License as published by
   7    the Free Software Foundation, either version 3 of the License, or
   8    (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful,
  11    but WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13    GNU General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  17
  18 #include <config.h>
  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 (struct pool *pool);
  58 typedef void stat_update  (void *stat, double w, double x);
  59 typedef double stat_get   (const struct per_var_data *, void *aux);
  60
  61 struct cell_spec
  62 {
  63   /* Printable title for output */
  64   const char *title;
  65
  66   /* Keyword for syntax */
  67   const char *keyword;
  68
  69   stat_create *sc;
  70   stat_update *su;
  71   stat_get *sd;
  72 };
  73
  74 struct harmonic_mean
  75 {
  76   double rsum;
  77   double n;
  78 };
  79
  80 static void *
  81 harmonic_create (struct pool *pool)
  82 {
  83   struct harmonic_mean *hm = pool_alloc (pool, sizeof *hm);
  84
  85   hm->rsum = 0;
  86   hm->n = 0;
  87
  88   return hm;
  89 }
  90
  91
  92 static void
  93 harmonic_update (void *stat, double w, double x)
  94 {
  95   struct harmonic_mean *hm = stat;
  96   hm->rsum  += w / x;
  97   hm->n += w;
  98 }
  99
 100
 101 static double
 102 harmonic_get (const struct per_var_data *pvd UNUSED, void *stat)
 103 {
 104   struct harmonic_mean *hm = stat;
 105
 106   return hm->n / hm->rsum;
 107 }
 108
 109 \f
 110
 111 struct geometric_mean
 112 {
 113   double prod;
 114   double n;
 115 };
 116
 117
 118 static void *
 119 geometric_create (struct pool *pool)
 120 {
 121   struct geometric_mean *gm = pool_alloc (pool, sizeof *gm);
 122
 123   gm->prod = 1.0;
 124   gm->n = 0;
 125
 126   return gm;
 127 }
 128
 129
 130 static void
 131 geometric_update (void *stat, double w, double x)
 132 {
 133   struct geometric_mean *gm = stat;
 134   gm->prod  *=  pow (x, w);
 135   gm->n += w;
 136 }
 137
 138
 139 static double
 140 geometric_get (const struct per_var_data *pvd UNUSED, void *stat)
 141 {
 142   struct geometric_mean *gm = stat;
 143
 144   return pow (gm->prod, 1.0 / gm->n);
 145 }
 146
 147 \f
 148
 149 static double
 150 sum_get (const struct per_var_data *pvd, void *stat UNUSED)
 151 {
 152   double n, mean;
 153
 154   moments1_calculate (pvd->mom, &n, &mean, 0, 0, 0);
 155
 156   return mean * n;
 157 }
 158
 159
 160 static double
 161 n_get (const struct per_var_data *pvd, void *stat UNUSED)
 162 {
 163   double n;
 164
 165   moments1_calculate (pvd->mom, &n, 0, 0, 0, 0);
 166
 167   return n;
 168 }
 169
 170 static double
 171 arithmean_get (const struct per_var_data *pvd, void *stat UNUSED)
 172 {
 173   double n, mean;
 174
 175   moments1_calculate (pvd->mom, &n, &mean, 0, 0, 0);
 176
 177   return mean;
 178 }
 179
 180 static double
 181 variance_get (const struct per_var_data *pvd, void *stat UNUSED)
 182 {
 183   double n, mean, variance;
 184
 185   moments1_calculate (pvd->mom, &n, &mean, &variance, 0, 0);
 186
 187   return variance;
 188 }
 189
 190
 191 static double
 192 stddev_get (const struct per_var_data *pvd, void *stat)
 193 {
 194   return sqrt (variance_get (pvd, stat));
 195 }
 196
 197
 198 \f
 199
 200 static double
 201 skew_get (const struct per_var_data *pvd, void *stat UNUSED)
 202 {
 203   double skew;
 204
 205   moments1_calculate (pvd->mom, NULL, NULL, NULL, &skew, 0);
 206
 207   return skew;
 208 }
 209
 210 static double
 211 sekurt_get (const struct per_var_data *pvd, void *stat UNUSED)
 212 {
 213   double n;
 214
 215   moments1_calculate (pvd->mom, &n, NULL, NULL, NULL, NULL);
 216
 217   return calc_sekurt (n);
 218 }
 219
 220 static double
 221 seskew_get (const struct per_var_data *pvd, void *stat UNUSED)
 222 {
 223   double n;
 224
 225   moments1_calculate (pvd->mom, &n, NULL, NULL, NULL, NULL);
 226
 227   return calc_seskew (n);
 228 }
 229
 230 static double
 231 kurt_get (const struct per_var_data *pvd, void *stat UNUSED)
 232 {
 233   double kurt;
 234
 235   moments1_calculate (pvd->mom, NULL, NULL, NULL, NULL, &kurt);
 236
 237   return kurt;
 238 }
 239
 240 static double
 241 semean_get (const struct per_var_data *pvd, void *stat UNUSED)
 242 {
 243   double n, var;
 244
 245   moments1_calculate (pvd->mom, &n, NULL, &var, NULL, NULL);
 246
 247   return sqrt (var / n);
 248 }
 249
 250 \f
 251
 252 static void *
 253 min_create (struct pool *pool)
 254 {
 255   double *r = pool_alloc (pool, sizeof *r);
 256
 257   *r = DBL_MAX;
 258
 259   return r;
 260 }
 261
 262 static void
 263 min_update (void *stat, double w UNUSED, double x)
 264 {
 265   double *r = stat;
 266
 267   if (x < *r)
 268     *r = x;
 269 }
 270
 271 static double
 272 min_get (const struct per_var_data *pvd UNUSED, void *stat)
 273 {
 274   double *r = stat;
 275
 276   return *r;
 277 }
 278
 279 static void *
 280 max_create (struct pool *pool)
 281 {
 282   double *r = pool_alloc (pool, sizeof *r);
 283
 284   *r = -DBL_MAX;
 285
 286   return r;
 287 }
 288
 289 static void
 290 max_update (void *stat, double w UNUSED, double x)
 291 {
 292   double *r = stat;
 293
 294   if (x > *r)
 295     *r = x;
 296 }
 297
 298 static double
 299 max_get (const struct per_var_data *pvd UNUSED, void *stat)
 300 {
 301   double *r = stat;
 302
 303   return *r;
 304 }
 305
 306 \f
 307
 308 struct range
 309 {
 310   double min;
 311   double max;
 312 };
 313
 314 static void *
 315 range_create (struct pool *pool)
 316 {
 317   struct range *r = pool_alloc (pool, sizeof *r);
 318
 319   r->min = DBL_MAX;
 320   r->max = -DBL_MAX;
 321
 322   return r;
 323 }
 324
 325 static void
 326 range_update (void *stat, double w UNUSED, double x)
 327 {
 328   struct range *r = stat;
 329
 330   if (x > r->max)
 331     r->max = x;
 332
 333   if (x < r->min)
 334     r->min = x;
 335 }
 336
 337 static double
 338 range_get (const struct per_var_data *pvd UNUSED, void *stat)
 339 {
 340   struct range *r = stat;
 341
 342   return r->max - r->min;
 343 }
 344
 345 \f
 346
 347 static void *
 348 last_create (struct pool *pool)
 349 {
 350   double *l = pool_alloc (pool, sizeof *l);
 351
 352   return l;
 353 }
 354
 355 static void
 356 last_update (void *stat, double w UNUSED, double x)
 357 {
 358   double *l = stat;
 359
 360   *l = x;
 361 }
 362
 363 static double
 364 last_get (const struct per_var_data *pvd UNUSED, void *stat)
 365 {
 366   double *l = stat;
 367
 368   return *l;
 369 }
 370
 371
 372 static void *
 373 first_create (struct pool *pool)
 374 {
 375   double *f = pool_alloc (pool, sizeof *f);
 376
 377   *f = SYSMIS;
 378
 379   return f;
 380 }
 381
 382 static void
 383 first_update (void *stat, double w UNUSED, double x)
 384 {
 385   double *f = stat;
 386
 387   if (*f == SYSMIS)
 388      *f = x;
 389 }
 390
 391 static double
 392 first_get (const struct per_var_data *pvd UNUSED,  void *stat)
 393 {
 394   double *f = stat;
 395
 396   return *f;
 397 }
 398
 399 enum
 400   {
 401     MEANS_MEAN = 0,
 402     MEANS_N,
 403     MEANS_STDDEV
 404   };
 405
 406 /* Table of cell_specs */
 407 static const struct cell_spec cell_spec[] = {
 408   {N_("Mean"), "MEAN", NULL, NULL, arithmean_get},
 409   {N_("N"), "COUNT", NULL, NULL, n_get},
 410   {N_("Std. Deviation"), "STDDEV", NULL, NULL, stddev_get},
 411 #if 0
 412   {N_("Median"), "MEDIAN", NULL, NULL, NULL},
 413   {N_("Group Median"), "GMEDIAN", NULL, NULL, NULL},
 414 #endif
 415   {N_("S.E. Mean"), "SEMEAN", NULL, NULL, semean_get},
 416   {N_("Sum"), "SUM", NULL, NULL, sum_get},
 417   {N_("Min"), "MIN", min_create, min_update, min_get},
 418   {N_("Max"), "MAX", max_create, max_update, max_get},
 419   {N_("Range"), "RANGE", range_create, range_update, range_get},
 420   {N_("Variance"), "VARIANCE", NULL, NULL, variance_get},
 421   {N_("Kurtosis"), "KURT", NULL, NULL, kurt_get},
 422   {N_("S.E. Kurt"), "SEKURT", NULL, NULL, sekurt_get},
 423   {N_("Skewness"), "SKEW", NULL, NULL, skew_get},
 424   {N_("S.E. Skew"), "SESKEW", NULL, NULL, seskew_get},
 425   {N_("First"), "FIRST", first_create, first_update, first_get},
 426   {N_("Last"), "LAST", last_create, last_update, last_get},
 427 #if 0
 428   {N_("Percent N"), "NPCT", NULL, NULL, NULL},
 429   {N_("Percent Sum"), "SPCT", NULL, NULL, NULL},
 430 #endif
 431   {N_("Harmonic Mean"), "HARMONIC", harmonic_create, harmonic_update, harmonic_get},
 432   {N_("Geom. Mean"), "GEOMETRIC", geometric_create, geometric_update, geometric_get}
 433 };
 434
 435 #define n_C (sizeof (cell_spec) / sizeof (struct cell_spec))
 436
 437
 438 struct summary
 439 {
 440   casenumber missing;
 441   casenumber non_missing;
 442 };
 443
 444
 445 struct layer
 446 {
 447   size_t n_factor_vars;
 448   const struct variable **factor_vars;
 449 };
 450
 451 /* The thing parsed after TABLES= */
 452 struct mtable
 453 {
 454   size_t n_dep_vars;
 455   const struct variable **dep_vars;
 456
 457   int n_layers;
 458   struct layer *layers;
 459
 460   struct interaction **interactions;
 461   struct summary *summary;
 462
 463   int ii;
 464
 465   struct categoricals *cats;
 466 };
 467
 468 struct means
 469 {
 470   const struct dictionary *dict;
 471
 472   struct mtable *table;
 473   size_t n_tables;
 474
 475   /* Missing value class for categorical variables */
 476   enum mv_class exclude;
 477
 478   /* Missing value class for dependent variables */
 479   enum mv_class dep_exclude;
 480
 481   bool listwise_exclude;
 482
 483   /* an array  indicating which statistics are to be calculated */
 484   int *cells;
 485
 486   /* Size of cells */
 487   int n_cells;
 488
 489   /* Pool on which cell functions may allocate data */
 490   struct pool *pool;
 491 };
 492
 493
 494 static void
 495 run_means (struct means *cmd, struct casereader *input,
 496            const struct dataset *ds);
 497
 498
 499
 500 static bool
 501 parse_means_table_syntax (struct lexer *lexer, const struct means *cmd, struct mtable *table)
 502 {
 503   table->ii = 0;
 504   table->n_layers = 0;
 505   table->layers = 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
 791   pool_destroy (means.pool);
 792   return CMD_SUCCESS;
 793
 794 error:
 795
 796   pool_destroy (means.pool);
 797   return CMD_FAILURE;
 798 }
 799
 800
 801 static bool
 802 is_missing (const struct means *cmd,
 803             const struct variable *dvar,
 804             const struct interaction *iact,
 805             const struct ccase *c)
 806 {
 807   if ( interaction_case_is_missing (iact, c, cmd->exclude) )
 808     return true;
 809
 810
 811   if (var_is_value_missing (dvar,
 812                             case_data (c, dvar),
 813                             cmd->dep_exclude))
 814     return true;
 815
 816   return false;
 817 }
 818
 819 static void output_case_processing_summary (const struct mtable *);
 820
 821 static void output_report (const struct means *, int, const struct mtable *);
 822
 823
 824 struct per_cat_data
 825 {
 826   struct per_var_data *pvd;
 827
 828   bool warn;
 829 };
 830
 831 static void *
 832 create_n (const void *aux1, void *aux2)
 833 {
 834   int i, v;
 835   const struct means *means = aux1;
 836   struct mtable *table = aux2;
 837   struct per_cat_data *per_cat_data = pool_malloc (means->pool, sizeof *per_cat_data);
 838
 839   struct per_var_data *pvd = pool_calloc (means->pool, table->n_dep_vars, sizeof *pvd);
 840
 841   for (v = 0; v < table->n_dep_vars; ++v)
 842     {
 843       enum moment maxmom = MOMENT_KURTOSIS;
 844       struct per_var_data *pp = &pvd[v];
 845
 846       pp->cell_stats = pool_calloc (means->pool, means->n_cells, sizeof *pp->cell_stats);
 847
 848
 849       for (i = 0; i < means->n_cells; ++i)
 850         {
 851           int csi = means->cells[i];
 852           const struct cell_spec *cs = &cell_spec[csi];
 853           if (cs->sc)
 854             {
 855               pp->cell_stats[i] = cs->sc (means->pool);
 856             }
 857         }
 858       pp->mom = moments1_create (maxmom);
 859     }
 860
 861
 862   per_cat_data->pvd = pvd;
 863   per_cat_data->warn = true;
 864   return per_cat_data;
 865 }
 866
 867 static void
 868 update_n (const void *aux1, void *aux2, void *user_data, const struct ccase *c, double weight)
 869 {
 870   int i;
 871   int v = 0;
 872   const struct means *means = aux1;
 873   struct mtable *table = aux2;
 874   struct per_cat_data *per_cat_data = user_data;
 875
 876   for (v = 0; v < table->n_dep_vars; ++v)
 877     {
 878       struct per_var_data *pvd = &per_cat_data->pvd[v];
 879
 880       const double x = case_data (c, table->dep_vars[v])->f;
 881
 882       for (i = 0; i < table->n_layers; ++i)
 883         {
 884           if ( is_missing (means, table->dep_vars[v],
 885                            table->interactions[i], c))
 886             goto end;
 887         }
 888
 889       for (i = 0; i < means->n_cells; ++i)
 890         {
 891           const int csi = means->cells[i];
 892           const struct cell_spec *cs = &cell_spec[csi];
 893
 894
 895           if (cs->su)
 896             cs->su (pvd->cell_stats[i],
 897                     weight, x);
 898         }
 899
 900       moments1_add (pvd->mom, x, weight);
 901
 902     end:
 903       continue;
 904     }
 905 }
 906
 907 static void
 908 calculate_n (const void *aux1, void *aux2, void *user_data)
 909 {
 910   int i;
 911   int v = 0;
 912   struct per_cat_data *per_cat_data = user_data;
 913   const struct means *means = aux1;
 914   struct mtable *table = aux2;
 915
 916   for (v = 0; v < table->n_dep_vars; ++v)
 917     {
 918       struct per_var_data *pvd = &per_cat_data->pvd[v];
 919       for (i = 0; i < means->n_cells; ++i)
 920         {
 921           int csi = means->cells[i];
 922           const struct cell_spec *cs = &cell_spec[csi];
 923
 924           if (cs->su)
 925             cs->sd (pvd, pvd->cell_stats[i]);
 926         }
 927     }
 928 }
 929
 930 static void
 931 run_means (struct means *cmd, struct casereader *input,
 932            const struct dataset *ds UNUSED)
 933 {
 934   int t;
 935   const struct variable *wv = dict_get_weight (cmd->dict);
 936   struct ccase *c;
 937   struct casereader *reader;
 938
 939   struct payload payload;
 940   payload.create = create_n;
 941   payload.update = update_n;
 942   payload.calculate = calculate_n;
 943   payload.destroy = NULL;
 944
 945   for (t = 0; t < cmd->n_tables; ++t)
 946   {
 947     struct mtable *table = &cmd->table[t];
 948     table->cats
 949       = categoricals_create (table->interactions,
 950                              table->n_layers, wv, cmd->dep_exclude, cmd->exclude);
 951
 952     categoricals_set_payload (table->cats, &payload, cmd, table);
 953   }
 954
 955   for (reader = input;
 956        (c = casereader_read (reader)) != NULL; case_unref (c))
 957     {
 958       for (t = 0; t < cmd->n_tables; ++t)
 959         {
 960           bool something_missing = false;
 961           int  v;
 962           struct mtable *table = &cmd->table[t];
 963
 964           for (v = 0; v < table->n_dep_vars; ++v)
 965             {
 966               int i;
 967               for (i = 0; i < table->n_layers; ++i)
 968                 {
 969                   const bool missing =
 970                     is_missing (cmd, table->dep_vars[v],
 971                                 table->interactions[i], c);
 972                   if (missing)
 973                     {
 974                       something_missing = true;
 975                       table->summary[v * table->n_layers + i].missing++;
 976                     }
 977                   else
 978                     table->summary[v * table->n_layers  + i].non_missing++;
 979                 }
 980             }
 981           if ( something_missing && cmd->listwise_exclude)
 982             continue;
 983
 984           categoricals_update (table->cats, c);
 985         }
 986     }
 987   casereader_destroy (reader);
 988
 989   for (t = 0; t < cmd->n_tables; ++t)
 990     {
 991       struct mtable *table = &cmd->table[t];
 992
 993       categoricals_done (table->cats);
 994     }
 995
 996
 997   for (t = 0; t < cmd->n_tables; ++t)
 998     {
 999       int i;
1000       const struct mtable *table = &cmd->table[t];
1001
1002       output_case_processing_summary (table);
1003
1004       for (i = 0; i < table->n_layers; ++i)
1005         {
1006           output_report (cmd, i, table);
1007         }
1008       categoricals_destroy (table->cats);
1009     }
1010
1011 }
1012
1013
1014
1015 static void
1016 output_case_processing_summary (const struct mtable *table)
1017 {
1018   int i, v;
1019   const int heading_columns = 1;
1020   const int heading_rows = 3;
1021   struct tab_table *t;
1022
1023   const int nr = heading_rows + table->n_layers * table->n_dep_vars;
1024   const int nc = 7;
1025
1026   t = tab_create (nc, nr);
1027   tab_title (t, _("Case Processing Summary"));
1028
1029   tab_headers (t, heading_columns, 0, heading_rows, 0);
1030
1031   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1032
1033   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1034   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1035
1036
1037   tab_joint_text (t, heading_columns, 0,
1038                   nc - 1, 0, TAB_CENTER | TAT_TITLE, _("Cases"));
1039
1040   tab_joint_text (t, 1, 1, 2, 1, TAB_CENTER | TAT_TITLE, _("Included"));
1041   tab_joint_text (t, 3, 1, 4, 1, TAB_CENTER | TAT_TITLE, _("Excluded"));
1042   tab_joint_text (t, 5, 1, 6, 1, TAB_CENTER | TAT_TITLE, _("Total"));
1043
1044   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1045   tab_hline (t, TAL_1, heading_columns, nc - 1, 2);
1046
1047
1048   for (i = 0; i < 3; ++i)
1049     {
1050       tab_text (t, heading_columns + i * 2, 2, TAB_CENTER | TAT_TITLE,
1051                 _("N"));
1052       tab_text (t, heading_columns + i * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1053                 _("Percent"));
1054     }
1055
1056   for (v = 0; v < table->n_dep_vars; ++v)
1057     {
1058       const struct variable *var = table->dep_vars[v];
1059       const char *dv_name = var_to_string (var);
1060       for (i = 0; i < table->n_layers; ++i)
1061         {
1062           const int row = v * table->n_layers + i;
1063           const struct interaction *iact = table->interactions[i];
1064           casenumber n_total;
1065
1066           struct string str;
1067           ds_init_cstr (&str, dv_name);
1068           ds_put_cstr (&str, ": ");
1069
1070           interaction_to_string (iact, &str);
1071
1072           tab_text (t, 0, row + heading_rows,
1073                     TAB_LEFT | TAT_TITLE, ds_cstr (&str));
1074
1075
1076           n_total = table->summary[row].missing +
1077             table->summary[row].non_missing;
1078
1079           tab_double (t, 1, row + heading_rows,
1080                       0, table->summary[row].non_missing, &F_8_0);
1081
1082           tab_text_format (t, 2, row + heading_rows,
1083                            0, _("%g%%"),
1084                            table->summary[row].non_missing / (double) n_total * 100.0);
1085
1086
1087           tab_double (t, 3, row + heading_rows,
1088                       0, table->summary[row].missing, &F_8_0);
1089
1090
1091           tab_text_format (t, 4, row + heading_rows,
1092                            0, _("%g%%"),
1093                            table->summary[row].missing / (double) n_total * 100.0);
1094
1095
1096           tab_double (t, 5, row + heading_rows,
1097                       0, table->summary[row].missing +
1098                       table->summary[row].non_missing, &F_8_0);
1099
1100           tab_text_format (t, 6, row + heading_rows,
1101                            0, _("%g%%"),
1102                            n_total / (double) n_total * 100.0);
1103
1104
1105           ds_destroy (&str);
1106         }
1107     }
1108
1109   tab_submit (t);
1110 }
1111
1112
1113 static void
1114 output_report (const struct means *cmd,  int iact_idx,
1115                const struct mtable *table)
1116 {
1117   int grp;
1118   int i;
1119
1120   const struct interaction *iact = table->interactions[iact_idx];
1121
1122   const int heading_columns = 1 + iact->n_vars;
1123   const int heading_rows = 1;
1124   struct tab_table *t;
1125
1126   const int n_cats = categoricals_n_count (table->cats, iact_idx);
1127
1128   const int nr = n_cats * table->n_dep_vars + heading_rows;
1129
1130   const int nc = heading_columns + cmd->n_cells;
1131
1132   t = tab_create (nc, nr);
1133   tab_title (t, _("Report"));
1134
1135   tab_headers (t, heading_columns, 0, heading_rows, 0);
1136
1137   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, nc - 1, nr - 1);
1138
1139   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1140   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1141
1142   for (i = 0; i < iact->n_vars; ++i)
1143     {
1144       tab_text (t, 1 + i, 0, TAB_CENTER | TAT_TITLE,
1145                 var_to_string (iact->vars[i]));
1146     }
1147
1148   for (i = 0; i < cmd->n_cells; ++i)
1149     {
1150       tab_text (t, heading_columns + i, 0,
1151                 TAB_CENTER | TAT_TITLE,
1152                 gettext (cell_spec[cmd->cells[i]].title));
1153     }
1154
1155
1156   for (i = 0; i < n_cats; ++i)
1157     {
1158       int v, dv;
1159       const struct ccase *c =
1160         categoricals_get_case_by_category_real (table->cats, iact_idx, i);
1161
1162       for (dv = 0; dv < table->n_dep_vars; ++dv)
1163         {
1164           tab_text (t, 0,
1165                     heading_rows + dv * n_cats,
1166                     TAB_RIGHT | TAT_TITLE,
1167                     var_to_string (table->dep_vars[dv])
1168                     );
1169
1170           if ( dv > 0)
1171             tab_hline (t, TAL_1, 0, nc - 1, heading_rows + dv * n_cats);
1172
1173           for (v = 0; v < iact->n_vars; ++v)
1174             {
1175               const struct variable *var = iact->vars[v];
1176               const union value *val = case_data (c, var);
1177               struct string str;
1178               ds_init_empty (&str);
1179               var_append_value_name (var, val, &str);
1180
1181               tab_text (t, 1 + v, heading_rows + dv * n_cats + i,
1182                         TAB_RIGHT | TAT_TITLE, ds_cstr (&str));
1183
1184               ds_destroy (&str);
1185             }
1186         }
1187     }
1188
1189   for (grp = 0; grp < n_cats; ++grp)
1190     {
1191       int dv;
1192       struct per_cat_data *per_cat_data =
1193         categoricals_get_user_data_by_category_real (table->cats, iact_idx, grp);
1194
1195       for (dv = 0; dv < table->n_dep_vars; ++dv)
1196         {
1197           const struct per_var_data *pvd = &per_cat_data->pvd[dv];
1198           for (i = 0; i < cmd->n_cells; ++i)
1199             {
1200               const int csi = cmd->cells[i];
1201               const struct cell_spec *cs = &cell_spec[csi];
1202
1203               double result = cs->sd (pvd, pvd->cell_stats[i]);
1204
1205               tab_double (t, heading_columns + i,
1206                           heading_rows + grp + dv * n_cats,
1207                           0, result, 0);
1208             }
1209         }
1210     }
1211
1212   tab_submit (t);
1213 }