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