src/language/stats/crosstabs.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 1997-9, 2000, 2006, 2009 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 /* FIXME:
  18
  19    - Pearson's R (but not Spearman!) is off a little.
  20    - T values for Spearman's R and Pearson's R are wrong.
  21    - How to calculate significance of symmetric and directional measures?
  22    - Asymmetric ASEs and T values for lambda are wrong.
  23    - ASE of Goodman and Kruskal's tau is not calculated.
  24    - ASE of symmetric somers' d is wrong.
  25    - Approx. T of uncertainty coefficient is wrong.
  26
  27 */
  28
  29 #include <config.h>
  30
  31 #include <ctype.h>
  32 #include <gsl/gsl_cdf.h>
  33 #include <stdlib.h>
  34 #include <stdio.h>
  35
  36 #include <data/case.h>
  37 #include <data/casegrouper.h>
  38 #include <data/casereader.h>
  39 #include <data/data-out.h>
  40 #include <data/dictionary.h>
  41 #include <data/format.h>
  42 #include <data/procedure.h>
  43 #include <data/value-labels.h>
  44 #include <data/variable.h>
  45 #include <language/command.h>
  46 #include <language/dictionary/split-file.h>
  47 #include <language/lexer/lexer.h>
  48 #include <language/lexer/variable-parser.h>
  49 #include <libpspp/array.h>
  50 #include <libpspp/assertion.h>
  51 #include <libpspp/compiler.h>
  52 #include <libpspp/hash.h>
  53 #include <libpspp/hmap.h>
  54 #include <libpspp/hmapx.h>
  55 #include <libpspp/message.h>
  56 #include <libpspp/misc.h>
  57 #include <libpspp/pool.h>
  58 #include <libpspp/str.h>
  59 #include <output/output.h>
  60 #include <output/table.h>
  61
  62 #include "minmax.h"
  63 #include "xalloc.h"
  64 #include "xsize.h"
  65
  66 #include "gettext.h"
  67 #define _(msgid) gettext (msgid)
  68 #define N_(msgid) msgid
  69
  70 /* (headers) */
  71
  72 /* (specification)
  73    crosstabs (crs_):
  74      *^tables=custom;
  75      +variables=custom;
  76      missing=miss:!table/include/report;
  77      +write[wr_]=none,cells,all;
  78      +format=fmt:!labels/nolabels/novallabs,
  79              val:!avalue/dvalue,
  80              indx:!noindex/index,
  81              tabl:!tables/notables,
  82              box:!box/nobox,
  83              pivot:!pivot/nopivot;
  84      +cells[cl_]=count,expected,row,column,total,residual,sresidual,
  85                  asresidual,all,none;
  86      +statistics[st_]=chisq,phi,cc,lambda,uc,none,btau,ctau,risk,gamma,d,
  87                       kappa,eta,corr,all.
  88 */
  89 /* (declarations) */
  90 /* (functions) */
  91
  92 /* Number of chi-square statistics. */
  93 #define N_CHISQ 5
  94
  95 /* Number of symmetric statistics. */
  96 #define N_SYMMETRIC 9
  97
  98 /* Number of directional statistics. */
  99 #define N_DIRECTIONAL 13
 100
 101 /* A single table entry for general mode. */
 102 struct table_entry
 103   {
 104     struct hmap_node node;      /* Entry in hash table. */
 105     double freq;                /* Frequency count. */
 106     union value values[1];      /* Values. */
 107   };
 108
 109 static size_t
 110 table_entry_size (size_t n_values)
 111 {
 112   return (offsetof (struct table_entry, values)
 113           + n_values * sizeof (union value));
 114 }
 115
 116 /* Indexes into the 'vars' member of struct pivot_table and
 117    struct crosstab member. */
 118 enum
 119   {
 120     ROW_VAR = 0,                /* Row variable. */
 121     COL_VAR = 1                 /* Column variable. */
 122     /* Higher indexes cause multiple tables to be output. */
 123   };
 124
 125 /* A crosstabulation of 2 or more variables. */
 126 struct pivot_table
 127   {
 128     struct fmt_spec weight_format; /* Format for weight variable. */
 129     double missing;             /* Weight of missing cases. */
 130
 131     /* Variables (2 or more). */
 132     int n_vars;
 133     const struct variable **vars;
 134
 135     /* Constants (0 or more). */
 136     int n_consts;
 137     const struct variable **const_vars;
 138     union value *const_values;
 139
 140     /* Data. */
 141     struct hmap data;
 142     struct table_entry **entries;
 143     size_t n_entries;
 144
 145     /* Column values, number of columns. */
 146     union value *cols;
 147     int n_cols;
 148
 149     /* Row values, number of rows. */
 150     union value *rows;
 151     int n_rows;
 152
 153     /* Number of statistically interesting columns/rows
 154        (columns/rows with data in them). */
 155     int ns_cols, ns_rows;
 156
 157     /* Matrix contents. */
 158     double *mat;                /* Matrix proper. */
 159     double *row_tot;            /* Row totals. */
 160     double *col_tot;            /* Column totals. */
 161     double total;               /* Grand total. */
 162   };
 163
 164 /* A crosstabulation of exactly 2 variables, conditional on zero
 165    or more other variables having given values. */
 166 struct crosstab
 167   {
 168     /* Case counts. */
 169     double missing;
 170
 171     /* Variables. */
 172     int n_vars;                 /* Number of variables (at least 2). */
 173     const struct variable **vars;
 174     union value *values;       /* Values of variables beyond 2. */
 175
 176     /* Data. */
 177     struct table_entry **entries;
 178     size_t n_entries;
 179
 180     /* Column values, number of columns. */
 181     union value *cols;
 182     int n_cols;
 183
 184     /* Row values, number of rows. */
 185     union value *rows;
 186     int n_rows;
 187
 188     /* Number of statistically interesting columns/rows
 189        (columns/rows with data in them). */
 190     int ns_cols, ns_rows;
 191
 192     /* Matrix contents. */
 193     double *mat;                /* Matrix proper. */
 194     double *row_tot;            /* Row totals. */
 195     double *col_tot;            /* Column totals. */
 196     double total;               /* Grand total. */
 197   };
 198
 199 /* Integer mode variable info. */
 200 struct var_range
 201   {
 202     int min;                    /* Minimum value. */
 203     int max;                    /* Maximum value + 1. */
 204     int count;                  /* max - min. */
 205   };
 206
 207 static inline struct var_range *
 208 get_var_range (const struct variable *v)
 209 {
 210   return var_get_aux (v);
 211 }
 212
 213 struct crosstabs_proc
 214   {
 215     enum { INTEGER, GENERAL } mode;
 216     enum mv_class exclude;
 217     bool pivot;
 218     bool bad_warn;
 219     struct fmt_spec weight_format;
 220
 221     /* Variables specifies on VARIABLES. */
 222     const struct variable **variables;
 223     size_t n_variables;
 224
 225     /* TABLES. */
 226     struct pivot_table *pivots;
 227     int n_pivots;
 228
 229     /* CELLS. */
 230     int n_cells;                /* Number of cells requested. */
 231     unsigned int cells;         /* Bit k is 1 if cell k is requested. */
 232     int a_cells[CRS_CL_count];  /* 0...n_cells-1 are the requested cells. */
 233
 234     /* STATISTICS. */
 235     unsigned int statistics;    /* Bit k is 1 if statistic k is requested. */
 236   };
 237
 238 static void
 239 init_proc (struct crosstabs_proc *proc, struct dataset *ds)
 240 {
 241   const struct variable *wv = dict_get_weight (dataset_dict (ds));
 242   proc->bad_warn = true;
 243   proc->variables = NULL;
 244   proc->n_variables = 0;
 245   proc->pivots = NULL;
 246   proc->n_pivots = 0;
 247   proc->weight_format = wv ? *var_get_print_format (wv) : F_8_0;
 248 }
 249
 250 static void
 251 free_proc (struct crosstabs_proc *proc UNUSED)
 252 {
 253   /* XXX */
 254 }
 255
 256 static int internal_cmd_crosstabs (struct lexer *lexer, struct dataset *ds,
 257                                    struct crosstabs_proc *);
 258 static bool should_tabulate_case (const struct pivot_table *,
 259                                   const struct ccase *, enum mv_class exclude);
 260 static void tabulate_general_case (struct pivot_table *, const struct ccase *,
 261                                    double weight);
 262 static void tabulate_integer_case (struct pivot_table *, const struct ccase *,
 263                                    double weight);
 264 static void postcalc (struct crosstabs_proc *);
 265 static void submit (struct crosstabs_proc *, struct pivot_table *,
 266                     struct tab_table *);
 267
 268 /* Parse and execute CROSSTABS, then clean up. */
 269 int
 270 cmd_crosstabs (struct lexer *lexer, struct dataset *ds)
 271 {
 272   struct crosstabs_proc proc;
 273   int result;
 274
 275   init_proc (&proc, ds);
 276   result = internal_cmd_crosstabs (lexer, ds, &proc);
 277   free_proc (&proc);
 278
 279   return result;
 280 }
 281
 282 /* Parses and executes the CROSSTABS procedure. */
 283 static int
 284 internal_cmd_crosstabs (struct lexer *lexer, struct dataset *ds,
 285                         struct crosstabs_proc *proc)
 286 {
 287   struct casegrouper *grouper;
 288   struct casereader *input, *group;
 289   struct cmd_crosstabs cmd;
 290   struct pivot_table *pt;
 291   bool ok;
 292   int i;
 293
 294   if (!parse_crosstabs (lexer, ds, &cmd, proc))
 295     return CMD_FAILURE;
 296
 297   proc->mode = proc->n_variables ? INTEGER : GENERAL;
 298
 299   /* CELLS. */
 300   if (!cmd.sbc_cells)
 301     proc->cells = 1u << CRS_CL_COUNT;
 302   else if (cmd.a_cells[CRS_CL_ALL])
 303     proc->cells = UINT_MAX;
 304   else
 305     {
 306       proc->cells = 0;
 307       for (i = 0; i < CRS_CL_count; i++)
 308         if (cmd.a_cells[i])
 309           proc->cells |= 1u << i;
 310       if (proc->cells == 0)
 311         proc->cells = ((1u << CRS_CL_COUNT)
 312                        | (1u << CRS_CL_ROW)
 313                        | (1u << CRS_CL_COLUMN)
 314                        | (1u << CRS_CL_TOTAL));
 315     }
 316   proc->cells &= ((1u << CRS_CL_count) - 1);
 317   proc->cells &= ~((1u << CRS_CL_NONE) | (1u << CRS_CL_ALL));
 318   proc->n_cells = 0;
 319   for (i = 0; i < CRS_CL_count; i++)
 320     if (proc->cells & (1u << i))
 321       proc->a_cells[proc->n_cells++] = i;
 322
 323   /* STATISTICS. */
 324   if (cmd.a_statistics[CRS_ST_ALL])
 325     proc->statistics = UINT_MAX;
 326   else if (cmd.sbc_statistics)
 327     {
 328       int i;
 329
 330       proc->statistics = 0;
 331       for (i = 0; i < CRS_ST_count; i++)
 332         if (cmd.a_statistics[i])
 333           proc->statistics |= 1u << i;
 334       if (proc->statistics == 0)
 335         proc->statistics |= 1u << CRS_ST_CHISQ;
 336     }
 337   else
 338     proc->statistics = 0;
 339
 340   /* MISSING. */
 341   proc->exclude = (cmd.miss == CRS_TABLE ? MV_ANY
 342                    : cmd.miss == CRS_INCLUDE ? MV_SYSTEM
 343                    : MV_NEVER);
 344   if (proc->mode == GENERAL && proc->mode == MV_NEVER)
 345     {
 346       msg (SE, _("Missing mode REPORT not allowed in general mode.  "
 347                  "Assuming MISSING=TABLE."));
 348       proc->mode = MV_ANY;
 349     }
 350
 351   /* PIVOT. */
 352   proc->pivot = cmd.pivot == CRS_PIVOT;
 353
 354   input = casereader_create_filter_weight (proc_open (ds), dataset_dict (ds),
 355                                            NULL, NULL);
 356   grouper = casegrouper_create_splits (input, dataset_dict (ds));
 357   while (casegrouper_get_next_group (grouper, &group))
 358     {
 359       struct ccase *c;
 360
 361       /* Output SPLIT FILE variables. */
 362       c = casereader_peek (group, 0);
 363       if (c != NULL)
 364         {
 365           output_split_file_values (ds, c);
 366           case_unref (c);
 367         }
 368
 369       /* Tabulate. */
 370       for (; (c = casereader_read (group)) != NULL; case_unref (c))
 371         for (pt = &proc->pivots[0]; pt < &proc->pivots[proc->n_pivots]; pt++)
 372           {
 373             double weight = dict_get_case_weight (dataset_dict (ds), c,
 374                                                   &proc->bad_warn);
 375             if (should_tabulate_case (pt, c, proc->exclude))
 376               {
 377                 if (proc->mode == GENERAL)
 378                   tabulate_general_case (pt, c, weight);
 379                 else
 380                   tabulate_integer_case (pt, c, weight);
 381               }
 382             else
 383               pt->missing += weight;
 384           }
 385       casereader_destroy (group);
 386
 387       /* Output. */
 388       postcalc (proc);
 389     }
 390   ok = casegrouper_destroy (grouper);
 391   ok = proc_commit (ds) && ok;
 392
 393   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
 394 }
 395
 396 /* Parses the TABLES subcommand. */
 397 static int
 398 crs_custom_tables (struct lexer *lexer, struct dataset *ds,
 399                    struct cmd_crosstabs *cmd UNUSED, void *proc_)
 400 {
 401   struct crosstabs_proc *proc = proc_;
 402   struct const_var_set *var_set;
 403   int n_by;
 404   const struct variable ***by = NULL;
 405   int *by_iter;
 406   size_t *by_nvar = NULL;
 407   size_t nx = 1;
 408   bool ok = false;
 409   int i;
 410
 411   /* Ensure that this is a TABLES subcommand. */
 412   if (!lex_match_id (lexer, "TABLES")
 413       && (lex_token (lexer) != T_ID ||
 414           dict_lookup_var (dataset_dict (ds), lex_tokid (lexer)) == NULL)
 415       && lex_token (lexer) != T_ALL)
 416     return 2;
 417   lex_match (lexer, '=');
 418
 419   if (proc->variables != NULL)
 420     var_set = const_var_set_create_from_array (proc->variables,
 421                                                proc->n_variables);
 422   else
 423     var_set = const_var_set_create_from_dict (dataset_dict (ds));
 424   assert (var_set != NULL);
 425
 426   for (n_by = 0; ;)
 427     {
 428       by = xnrealloc (by, n_by + 1, sizeof *by);
 429       by_nvar = xnrealloc (by_nvar, n_by + 1, sizeof *by_nvar);
 430       if (!parse_const_var_set_vars (lexer, var_set, &by[n_by], &by_nvar[n_by],
 431                                      PV_NO_DUPLICATE | PV_NO_SCRATCH))
 432         goto done;
 433       if (xalloc_oversized (nx, by_nvar[n_by]))
 434         {
 435           msg (SE, _("Too many cross-tabulation variables or dimensions."));
 436           goto done;
 437         }
 438       nx *= by_nvar[n_by];
 439       n_by++;
 440
 441       if (!lex_match (lexer, T_BY))
 442         {
 443           if (n_by < 2)
 444             {
 445               lex_error (lexer, _("expecting BY"));
 446               goto done;
 447             }
 448           else
 449             break;
 450         }
 451     }
 452
 453   by_iter = xcalloc (n_by, sizeof *by_iter);
 454   proc->pivots = xnrealloc (proc->pivots,
 455                             proc->n_pivots + nx, sizeof *proc->pivots);
 456   for (i = 0; i < nx; i++)
 457     {
 458       struct pivot_table *pt = &proc->pivots[proc->n_pivots++];
 459       int j;
 460
 461       pt->weight_format = proc->weight_format;
 462       pt->missing = 0.;
 463       pt->n_vars = n_by;
 464       pt->vars = xmalloc (n_by * sizeof *pt->vars);
 465       pt->n_consts = 0;
 466       pt->const_vars = NULL;
 467       pt->const_values = NULL;
 468       hmap_init (&pt->data);
 469       pt->entries = NULL;
 470       pt->n_entries = 0;
 471
 472       for (j = 0; j < n_by; j++)
 473         pt->vars[j] = by[j][by_iter[j]];
 474
 475       for (j = n_by - 1; j >= 0; j--)
 476         {
 477           if (++by_iter[j] < by_nvar[j])
 478             break;
 479           by_iter[j] = 0;
 480         }
 481     }
 482   free (by_iter);
 483   ok = true;
 484
 485 done:
 486   /* All return paths lead here. */
 487   for (i = 0; i < n_by; i++)
 488     free (by[i]);
 489   free (by);
 490   free (by_nvar);
 491
 492   const_var_set_destroy (var_set);
 493
 494   return ok;
 495 }
 496
 497 /* Parses the VARIABLES subcommand. */
 498 static int
 499 crs_custom_variables (struct lexer *lexer, struct dataset *ds,
 500                       struct cmd_crosstabs *cmd UNUSED, void *proc_)
 501 {
 502   struct crosstabs_proc *proc = proc_;
 503   if (proc->n_pivots)
 504     {
 505       msg (SE, _("VARIABLES must be specified before TABLES."));
 506       return 0;
 507     }
 508
 509   lex_match (lexer, '=');
 510
 511   for (;;)
 512     {
 513       size_t orig_nv = proc->n_variables;
 514       size_t i;
 515
 516       long min, max;
 517
 518       if (!parse_variables_const (lexer, dataset_dict (ds),
 519                                   &proc->variables, &proc->n_variables,
 520                                   (PV_APPEND | PV_NUMERIC
 521                                    | PV_NO_DUPLICATE | PV_NO_SCRATCH)))
 522         return 0;
 523
 524       if (lex_token (lexer) != '(')
 525         {
 526           lex_error (lexer, "expecting `('");
 527           goto lossage;
 528         }
 529       lex_get (lexer);
 530
 531       if (!lex_force_int (lexer))
 532         goto lossage;
 533       min = lex_integer (lexer);
 534       lex_get (lexer);
 535
 536       lex_match (lexer, ',');
 537
 538       if (!lex_force_int (lexer))
 539         goto lossage;
 540       max = lex_integer (lexer);
 541       if (max < min)
 542         {
 543           msg (SE, _("Maximum value (%ld) less than minimum value (%ld)."),
 544                max, min);
 545           goto lossage;
 546         }
 547       lex_get (lexer);
 548
 549       if (lex_token (lexer) != ')')
 550         {
 551           lex_error (lexer, "expecting `)'");
 552           goto lossage;
 553         }
 554       lex_get (lexer);
 555
 556       for (i = orig_nv; i < proc->n_variables; i++)
 557         {
 558           struct var_range *vr = xmalloc (sizeof *vr);
 559           vr->min = min;
 560           vr->max = max + 1.;
 561           vr->count = max - min + 1;
 562           var_attach_aux (proc->variables[i], vr, var_dtor_free);
 563         }
 564
 565       if (lex_token (lexer) == '/')
 566         break;
 567     }
 568
 569   return 1;
 570
 571  lossage:
 572   free (proc->variables);
 573   proc->variables = NULL;
 574   proc->n_variables = 0;
 575   return 0;
 576 }
 577 \f
 578 /* Data file processing. */
 579
 580 static bool
 581 should_tabulate_case (const struct pivot_table *pt, const struct ccase *c,
 582                       enum mv_class exclude)
 583 {
 584   int j;
 585   for (j = 0; j < pt->n_vars; j++)
 586     {
 587       const struct variable *var = pt->vars[j];
 588       struct var_range *range = get_var_range (var);
 589
 590       if (var_is_value_missing (var, case_data (c, var), exclude))
 591         return false;
 592
 593       if (range != NULL)
 594         {
 595           double num = case_num (c, var);
 596           if (num < range->min || num > range->max)
 597             return false;
 598         }
 599     }
 600   return true;
 601 }
 602
 603 static void
 604 tabulate_integer_case (struct pivot_table *pt, const struct ccase *c,
 605                        double weight)
 606 {
 607   struct table_entry *te;
 608   size_t hash;
 609   int j;
 610
 611   hash = 0;
 612   for (j = 0; j < pt->n_vars; j++)
 613     {
 614       /* Throw away fractional parts of values. */
 615       hash = hash_int (case_num (c, pt->vars[j]), hash);
 616     }
 617
 618   HMAP_FOR_EACH_WITH_HASH (te, struct table_entry, node, hash, &pt->data)
 619     {
 620       for (j = 0; j < pt->n_vars; j++)
 621         if ((int) case_num (c, pt->vars[j]) != (int) te->values[j].f)
 622           goto no_match;
 623
 624       /* Found an existing entry. */
 625       te->freq += weight;
 626       return;
 627
 628     no_match: ;
 629     }
 630
 631   /* No existing entry.  Create a new one. */
 632   te = xmalloc (table_entry_size (pt->n_vars));
 633   te->freq = weight;
 634   for (j = 0; j < pt->n_vars; j++)
 635     te->values[j].f = (int) case_num (c, pt->vars[j]);
 636   hmap_insert (&pt->data, &te->node, hash);
 637 }
 638
 639 static void
 640 tabulate_general_case (struct pivot_table *pt, const struct ccase *c,
 641                        double weight)
 642 {
 643   struct table_entry *te;
 644   size_t hash;
 645   int j;
 646
 647   hash = 0;
 648   for (j = 0; j < pt->n_vars; j++)
 649     {
 650       const struct variable *var = pt->vars[j];
 651       hash = value_hash (case_data (c, var), var_get_width (var), hash);
 652     }
 653
 654   HMAP_FOR_EACH_WITH_HASH (te, struct table_entry, node, hash, &pt->data)
 655     {
 656       for (j = 0; j < pt->n_vars; j++)
 657         {
 658           const struct variable *var = pt->vars[j];
 659           if (!value_equal (case_data (c, var), &te->values[j],
 660                             var_get_width (var)))
 661             goto no_match;
 662         }
 663
 664       /* Found an existing entry. */
 665       te->freq += weight;
 666       return;
 667
 668     no_match: ;
 669     }
 670
 671   /* No existing entry.  Create a new one. */
 672   te = xmalloc (table_entry_size (pt->n_vars));
 673   te->freq = weight;
 674   for (j = 0; j < pt->n_vars; j++)
 675     {
 676       const struct variable *var = pt->vars[j];
 677       int width = var_get_width (var);
 678       value_init (&te->values[j], width);
 679       value_copy (&te->values[j], case_data (c, var), width);
 680     }
 681   hmap_insert (&pt->data, &te->node, hash);
 682 }
 683 \f
 684 /* Post-data reading calculations. */
 685
 686 static int compare_table_entry_vars_3way (const struct table_entry *a,
 687                                           const struct table_entry *b,
 688                                           const struct pivot_table *pt,
 689                                           int idx0, int idx1);
 690 static int compare_table_entry_3way (const void *ap_, const void *bp_,
 691                                      const void *pt_);
 692 static void enum_var_values (const struct pivot_table *, int var_idx,
 693                              union value **valuesp, int *n_values);
 694 static void output_pivot_table (struct crosstabs_proc *,
 695                                 struct pivot_table *);
 696 static void make_pivot_table_subset (struct pivot_table *pt,
 697                                      size_t row0, size_t row1,
 698                                      struct pivot_table *subset);
 699 static void make_summary_table (struct crosstabs_proc *);
 700 static bool find_crosstab (struct pivot_table *, size_t *row0p, size_t *row1p);
 701
 702 static void
 703 postcalc (struct crosstabs_proc *proc)
 704 {
 705   struct pivot_table *pt;
 706
 707   /* Convert hash tables into sorted arrays of entries. */
 708   for (pt = &proc->pivots[0]; pt < &proc->pivots[proc->n_pivots]; pt++)
 709     {
 710       struct table_entry *e;
 711       size_t i;
 712
 713       pt->n_entries = hmap_count (&pt->data);
 714       pt->entries = xnmalloc (pt->n_entries, sizeof *pt->entries);
 715       i = 0;
 716       HMAP_FOR_EACH (e, struct table_entry, node, &pt->data)
 717         pt->entries[i++] = e;
 718       hmap_destroy (&pt->data);
 719
 720       sort (pt->entries, pt->n_entries, sizeof *pt->entries,
 721             compare_table_entry_3way, pt);
 722     }
 723
 724   make_summary_table (proc);
 725
 726   /* Output each pivot table. */
 727   for (pt = &proc->pivots[0]; pt < &proc->pivots[proc->n_pivots]; pt++)
 728     {
 729       if (proc->pivot || pt->n_vars == 2)
 730         output_pivot_table (proc, pt);
 731       else
 732         {
 733           size_t row0 = 0, row1 = 0;
 734           while (find_crosstab (pt, &row0, &row1))
 735             {
 736               struct pivot_table subset;
 737               make_pivot_table_subset (pt, row0, row1, &subset);
 738               output_pivot_table (proc, &subset);
 739             }
 740         }
 741     }
 742
 743   /* XXX clear output and prepare for next split file. */
 744 }
 745
 746 static void
 747 make_pivot_table_subset (struct pivot_table *pt, size_t row0, size_t row1,
 748                          struct pivot_table *subset)
 749 {
 750   *subset = *pt;
 751   if (pt->n_vars > 2)
 752     {
 753       assert (pt->n_consts == 0);
 754       subset->missing = pt->missing;
 755       subset->n_vars = 2;
 756       subset->vars = pt->vars;
 757       subset->n_consts = pt->n_vars - 2;
 758       subset->const_vars = pt->vars + 2;
 759       subset->const_values = &pt->entries[row0]->values[2];
 760     }
 761   subset->entries = &pt->entries[row0];
 762   subset->n_entries = row1 - row0;
 763 }
 764
 765 static int
 766 compare_table_entry_var_3way (const struct table_entry *a,
 767                               const struct table_entry *b,
 768                               const struct pivot_table *pt,
 769                               int idx)
 770 {
 771   return value_compare_3way (&a->values[idx], &b->values[idx],
 772                              var_get_width (pt->vars[idx]));
 773 }
 774
 775 static int
 776 compare_table_entry_vars_3way (const struct table_entry *a,
 777                                const struct table_entry *b,
 778                                const struct pivot_table *pt,
 779                                int idx0, int idx1)
 780 {
 781   int i;
 782
 783   for (i = idx1 - 1; i >= idx0; i--)
 784     {
 785       int cmp = compare_table_entry_var_3way (a, b, pt, i);
 786       if (cmp != 0)
 787         return cmp;
 788     }
 789   return 0;
 790 }
 791
 792 /* Compare the struct table_entry at *AP to the one at *BP and
 793    return a strcmp()-type result. */
 794 static int
 795 compare_table_entry_3way (const void *ap_, const void *bp_, const void *pt_)
 796 {
 797   const struct table_entry *const *ap = ap_;
 798   const struct table_entry *const *bp = bp_;
 799   const struct table_entry *a = *ap;
 800   const struct table_entry *b = *bp;
 801   const struct pivot_table *pt = pt_;
 802   int cmp;
 803
 804   cmp = compare_table_entry_vars_3way (a, b, pt, 2, pt->n_vars);
 805   if (cmp != 0)
 806     return cmp;
 807
 808   cmp = compare_table_entry_var_3way (a, b, pt, ROW_VAR);
 809   if (cmp != 0)
 810     return cmp;
 811
 812   return compare_table_entry_var_3way (a, b, pt, COL_VAR);
 813 }
 814
 815 static int
 816 find_first_difference (const struct pivot_table *pt, size_t row)
 817 {
 818   if (row == 0)
 819     return pt->n_vars - 1;
 820   else
 821     {
 822       const struct table_entry *a = pt->entries[row];
 823       const struct table_entry *b = pt->entries[row - 1];
 824       int col;
 825
 826       for (col = pt->n_vars - 1; col >= 0; col--)
 827         if (compare_table_entry_var_3way (a, b, pt, col))
 828           return col;
 829       NOT_REACHED ();
 830     }
 831 }
 832
 833 /* Output a table summarizing the cases processed. */
 834 static void
 835 make_summary_table (struct crosstabs_proc *proc)
 836 {
 837   struct tab_table *summary;
 838   struct pivot_table *pt;
 839   struct string name;
 840   int i;
 841
 842   summary = tab_create (7, 3 + proc->n_pivots, 1);
 843   tab_title (summary, _("Summary."));
 844   tab_headers (summary, 1, 0, 3, 0);
 845   tab_joint_text (summary, 1, 0, 6, 0, TAB_CENTER, _("Cases"));
 846   tab_joint_text (summary, 1, 1, 2, 1, TAB_CENTER, _("Valid"));
 847   tab_joint_text (summary, 3, 1, 4, 1, TAB_CENTER, _("Missing"));
 848   tab_joint_text (summary, 5, 1, 6, 1, TAB_CENTER, _("Total"));
 849   tab_hline (summary, TAL_1, 1, 6, 1);
 850   tab_hline (summary, TAL_1, 1, 6, 2);
 851   tab_vline (summary, TAL_1, 3, 1, 1);
 852   tab_vline (summary, TAL_1, 5, 1, 1);
 853   for (i = 0; i < 3; i++)
 854     {
 855       tab_text (summary, 1 + i * 2, 2, TAB_RIGHT, _("N"));
 856       tab_text (summary, 2 + i * 2, 2, TAB_RIGHT, _("Percent"));
 857     }
 858   tab_offset (summary, 0, 3);
 859
 860   ds_init_empty (&name);
 861   for (pt = &proc->pivots[0]; pt < &proc->pivots[proc->n_pivots]; pt++)
 862     {
 863       double valid;
 864       double n[3];
 865       size_t i;
 866
 867       tab_hline (summary, TAL_1, 0, 6, 0);
 868
 869       ds_clear (&name);
 870       for (i = 0; i < pt->n_vars; i++)
 871         {
 872           if (i > 0)
 873             ds_put_cstr (&name, " * ");
 874           ds_put_cstr (&name, var_to_string (pt->vars[i]));
 875         }
 876       tab_text (summary, 0, 0, TAB_LEFT, ds_cstr (&name));
 877
 878       valid = 0.;
 879       for (i = 0; i < pt->n_entries; i++)
 880         valid += pt->entries[i]->freq;
 881
 882       n[0] = valid;
 883       n[1] = pt->missing;
 884       n[2] = n[0] + n[1];
 885       for (i = 0; i < 3; i++)
 886         {
 887           tab_double (summary, i * 2 + 1, 0, TAB_RIGHT, n[i],
 888                       &proc->weight_format);
 889           tab_text (summary, i * 2 + 2, 0, TAB_RIGHT | TAT_PRINTF, "%.1f%%",
 890                     n[i] / n[2] * 100.);
 891         }
 892
 893       tab_next_row (summary);
 894     }
 895   ds_destroy (&name);
 896
 897   submit (proc, NULL, summary);
 898 }
 899 \f
 900 /* Output. */
 901
 902 static struct tab_table *create_crosstab_table (struct crosstabs_proc *,
 903                                                 struct pivot_table *);
 904 static struct tab_table *create_chisq_table (struct pivot_table *);
 905 static struct tab_table *create_sym_table (struct pivot_table *);
 906 static struct tab_table *create_risk_table (struct pivot_table *);
 907 static struct tab_table *create_direct_table (struct pivot_table *);
 908 static void display_dimensions (struct crosstabs_proc *, struct pivot_table *,
 909                                 struct tab_table *, int first_difference);
 910 static void display_crosstabulation (struct crosstabs_proc *,
 911                                      struct pivot_table *,
 912                                      struct tab_table *);
 913 static void display_chisq (struct pivot_table *, struct tab_table *,
 914                            bool *showed_fisher);
 915 static void display_symmetric (struct crosstabs_proc *, struct pivot_table *,
 916                                struct tab_table *);
 917 static void display_risk (struct pivot_table *, struct tab_table *);
 918 static void display_directional (struct crosstabs_proc *, struct pivot_table *,
 919                                  struct tab_table *);
 920 static void crosstabs_dim (struct tab_table *, struct outp_driver *,
 921                            void *proc);
 922 static void table_value_missing (struct crosstabs_proc *proc,
 923                                  struct tab_table *table, int c, int r,
 924                                  unsigned char opt, const union value *v,
 925                                  const struct variable *var);
 926 static void delete_missing (struct pivot_table *);
 927 static void build_matrix (struct pivot_table *);
 928
 929 /* Output pivot table beginning at PB and continuing until PE,
 930    exclusive.  For efficiency, *MATP is a pointer to a matrix that can
 931    hold *MAXROWS entries. */
 932 static void
 933 output_pivot_table (struct crosstabs_proc *proc, struct pivot_table *pt)
 934 {
 935   struct tab_table *table = NULL; /* Crosstabulation table. */
 936   struct tab_table *chisq = NULL; /* Chi-square table. */
 937   bool showed_fisher = false;
 938   struct tab_table *sym = NULL;   /* Symmetric measures table. */
 939   struct tab_table *risk = NULL;  /* Risk estimate table. */
 940   struct tab_table *direct = NULL; /* Directional measures table. */
 941   size_t row0, row1;
 942
 943   enum_var_values (pt, COL_VAR, &pt->cols, &pt->n_cols);
 944
 945   if (proc->cells)
 946     table = create_crosstab_table (proc, pt);
 947   if (proc->statistics & (1u << CRS_ST_CHISQ))
 948     chisq = create_chisq_table (pt);
 949   if (proc->statistics & ((1u << CRS_ST_PHI) | (1u << CRS_ST_CC)
 950                           | (1u << CRS_ST_BTAU) | (1u << CRS_ST_CTAU)
 951                           | (1u << CRS_ST_GAMMA) | (1u << CRS_ST_CORR)
 952                           | (1u << CRS_ST_KAPPA)))
 953     sym = create_sym_table (pt);
 954   if (proc->statistics & (1u << CRS_ST_RISK))
 955     risk = create_risk_table (pt);
 956   if (proc->statistics & ((1u << CRS_ST_LAMBDA) | (1u << CRS_ST_UC)
 957                           | (1u << CRS_ST_D) | (1u << CRS_ST_ETA)))
 958     direct = create_direct_table (pt);
 959
 960   row0 = row1 = 0;
 961   while (find_crosstab (pt, &row0, &row1))
 962     {
 963       struct pivot_table x;
 964       int first_difference;
 965
 966       make_pivot_table_subset (pt, row0, row1, &x);
 967
 968       /* Find all the row variable values. */
 969       enum_var_values (&x, ROW_VAR, &x.rows, &x.n_rows);
 970
 971       if (size_overflow_p (xtimes (xtimes (x.n_rows, x.n_cols),
 972                                    sizeof (double))))
 973         xalloc_die ();
 974       x.row_tot = xmalloc (x.n_rows * sizeof *x.row_tot);
 975       x.col_tot = xmalloc (x.n_cols * sizeof *x.col_tot);
 976       x.mat = xmalloc (x.n_rows * x.n_cols * sizeof *x.mat);
 977
 978       /* Allocate table space for the matrix. */
 979       if (table
 980           && tab_row (table) + (x.n_rows + 1) * proc->n_cells > tab_nr (table))
 981         tab_realloc (table, -1,
 982                      MAX (tab_nr (table) + (x.n_rows + 1) * proc->n_cells,
 983                           tab_nr (table) * pt->n_entries / x.n_entries));
 984
 985       build_matrix (&x);
 986
 987       /* Find the first variable that differs from the last subtable. */
 988       first_difference = find_first_difference (pt, row0);
 989       if (table)
 990         {
 991           display_dimensions (proc, &x, table, first_difference);
 992           display_crosstabulation (proc, &x, table);
 993         }
 994
 995       if (proc->exclude == MV_NEVER)
 996         delete_missing (&x);
 997
 998       if (chisq)
 999         {
1000           display_dimensions (proc, &x, chisq, first_difference);
1001           display_chisq (pt, chisq, &showed_fisher);
1002         }
1003       if (sym)
1004         {
1005           display_dimensions (proc, &x, sym, first_difference);
1006           display_symmetric (proc, pt, sym);
1007         }
1008       if (risk)
1009         {
1010           display_dimensions (proc, &x, risk, first_difference);
1011           display_risk (pt, risk);
1012         }
1013       if (direct)
1014         {
1015           display_dimensions (proc, &x, direct, first_difference);
1016           display_directional (proc, pt, direct);
1017         }
1018
1019       /* XXX Free data in x. */
1020       free (x.rows);
1021     }
1022
1023   submit (proc, NULL, table);
1024
1025   if (chisq)
1026     {
1027       if (!showed_fisher)
1028         tab_resize (chisq, 4 + (pt->n_vars - 2), -1);
1029       submit (proc, pt, chisq);
1030     }
1031
1032   submit (proc, pt, sym);
1033   submit (proc, pt, risk);
1034   submit (proc, pt, direct);
1035
1036   free (pt->cols);
1037 }
1038
1039 static void
1040 build_matrix (struct pivot_table *x)
1041 {
1042   const int col_var_width = var_get_width (x->vars[COL_VAR]);
1043   const int row_var_width = var_get_width (x->vars[ROW_VAR]);
1044   int col, row;
1045   double *mp;
1046   struct table_entry **p;
1047
1048   mp = x->mat;
1049   col = row = 0;
1050   for (p = x->entries; p < &x->entries[x->n_entries]; p++)
1051     {
1052       const struct table_entry *te = *p;
1053
1054       while (!value_equal (&x->rows[row], &te->values[ROW_VAR], row_var_width))
1055         {
1056           for (; col < x->n_cols; col++)
1057             *mp++ = 0.0;
1058           col = 0;
1059           row++;
1060         }
1061
1062       while (!value_equal (&x->cols[col], &te->values[COL_VAR], col_var_width))
1063         {
1064           *mp++ = 0.0;
1065           col++;
1066         }
1067
1068       *mp++ = te->freq;
1069       if (++col >= x->n_cols)
1070         {
1071           col = 0;
1072           row++;
1073         }
1074     }
1075   while (mp < &x->mat[x->n_cols * x->n_rows])
1076     *mp++ = 0.0;
1077   assert (mp == &x->mat[x->n_cols * x->n_rows]);
1078
1079   /* Column totals, row totals, ns_rows. */
1080   mp = x->mat;
1081   for (col = 0; col < x->n_cols; col++)
1082     x->col_tot[col] = 0.0;
1083   for (row = 0; row < x->n_rows; row++)
1084     x->row_tot[row] = 0.0;
1085   x->ns_rows = 0;
1086   for (row = 0; row < x->n_rows; row++)
1087     {
1088       bool row_is_empty = true;
1089       for (col = 0; col < x->n_cols; col++)
1090         {
1091           if (*mp != 0.0)
1092             {
1093               row_is_empty = false;
1094               x->col_tot[col] += *mp;
1095               x->row_tot[row] += *mp;
1096             }
1097           mp++;
1098         }
1099       if (!row_is_empty)
1100         x->ns_rows++;
1101     }
1102   assert (mp == &x->mat[x->n_cols * x->n_rows]);
1103
1104   /* ns_cols. */
1105   x->ns_cols = 0;
1106   for (col = 0; col < x->n_cols; col++)
1107     for (row = 0; row < x->n_rows; row++)
1108       if (x->mat[col + row * x->n_cols] != 0.0)
1109         {
1110           x->ns_cols++;
1111           break;
1112         }
1113
1114   /* Grand total. */
1115   x->total = 0.0;
1116   for (col = 0; col < x->n_cols; col++)
1117     x->total += x->col_tot[col];
1118 }
1119
1120 static struct tab_table *
1121 create_crosstab_table (struct crosstabs_proc *proc, struct pivot_table *pt)
1122 {
1123   struct tuple
1124     {
1125       int value;
1126       const char *name;
1127     };
1128   static const struct tuple names[] =
1129     {
1130       {CRS_CL_COUNT, N_("count")},
1131       {CRS_CL_ROW, N_("row %")},
1132       {CRS_CL_COLUMN, N_("column %")},
1133       {CRS_CL_TOTAL, N_("total %")},
1134       {CRS_CL_EXPECTED, N_("expected")},
1135       {CRS_CL_RESIDUAL, N_("residual")},
1136       {CRS_CL_SRESIDUAL, N_("std. resid.")},
1137       {CRS_CL_ASRESIDUAL, N_("adj. resid.")},
1138     };
1139   const int n_names = sizeof names / sizeof *names;
1140   const struct tuple *t;
1141
1142   struct tab_table *table;
1143   struct string title;
1144   int i;
1145
1146   table = tab_create (pt->n_consts + 1 + pt->n_cols + 1,
1147                       (pt->n_entries / pt->n_cols) * 3 / 2 * proc->n_cells + 10,
1148                       true);
1149   tab_headers (table, pt->n_consts + 1, 0, 2, 0);
1150
1151   /* First header line. */
1152   tab_joint_text (table, pt->n_consts + 1, 0,
1153                   (pt->n_consts + 1) + (pt->n_cols - 1), 0,
1154                   TAB_CENTER | TAT_TITLE, var_get_name (pt->vars[COL_VAR]));
1155
1156   tab_hline (table, TAL_1, pt->n_consts + 1,
1157              pt->n_consts + 2 + pt->n_cols - 2, 1);
1158
1159   /* Second header line. */
1160   for (i = 2; i < pt->n_consts + 2; i++)
1161     tab_joint_text (table, pt->n_consts + 2 - i - 1, 0,
1162                     pt->n_consts + 2 - i - 1, 1,
1163                     TAB_RIGHT | TAT_TITLE, var_to_string (pt->vars[i]));
1164   tab_text (table, pt->n_consts + 2 - 2, 1, TAB_RIGHT | TAT_TITLE,
1165             var_get_name (pt->vars[ROW_VAR]));
1166   for (i = 0; i < pt->n_cols; i++)
1167     table_value_missing (proc, table, pt->n_consts + 2 + i - 1, 1, TAB_RIGHT,
1168                          &pt->cols[i], pt->vars[COL_VAR]);
1169   tab_text (table, pt->n_consts + 2 + pt->n_cols - 1, 1, TAB_CENTER, _("Total"));
1170
1171   tab_hline (table, TAL_1, 0, pt->n_consts + 2 + pt->n_cols - 1, 2);
1172   tab_vline (table, TAL_1, pt->n_consts + 2 + pt->n_cols - 1, 0, 1);
1173
1174   /* Title. */
1175   ds_init_empty (&title);
1176   for (i = 0; i < pt->n_consts + 2; i++)
1177     {
1178       if (i)
1179         ds_put_cstr (&title, " * ");
1180       ds_put_cstr (&title, var_get_name (pt->vars[i]));
1181     }
1182   for (i = 0; i < pt->n_consts; i++)
1183     {
1184       const struct variable *var = pt->const_vars[i];
1185       ds_put_format (&title, ", %s=", var_get_name (var));
1186       data_out (&pt->const_values[i], var_get_print_format (var),
1187                 ds_put_uninit (&title, var_get_width (var)));
1188       /* XXX remove any leading space in what was just inserted.  */
1189     }
1190
1191   ds_put_cstr (&title, " [");
1192   i = 0;
1193   for (t = names; t < &names[n_names]; t++)
1194     if (proc->cells & (1u << t->value))
1195       {
1196         if (i++)
1197           ds_put_cstr (&title, ", ");
1198         ds_put_cstr (&title, gettext (t->name));
1199       }
1200   ds_put_cstr (&title, "].");
1201
1202   tab_title (table, "%s", ds_cstr (&title));
1203   ds_destroy (&title);
1204
1205   tab_offset (table, 0, 2);
1206   return table;
1207 }
1208
1209 static struct tab_table *
1210 create_chisq_table (struct pivot_table *pt)
1211 {
1212   struct tab_table *chisq;
1213
1214   chisq = tab_create (6 + (pt->n_vars - 2),
1215                       pt->n_entries / pt->n_cols * 3 / 2 * N_CHISQ + 10,
1216                       1);
1217   tab_headers (chisq, 1 + (pt->n_vars - 2), 0, 1, 0);
1218
1219   tab_title (chisq, _("Chi-square tests."));
1220
1221   tab_offset (chisq, pt->n_vars - 2, 0);
1222   tab_text (chisq, 0, 0, TAB_LEFT | TAT_TITLE, _("Statistic"));
1223   tab_text (chisq, 1, 0, TAB_RIGHT | TAT_TITLE, _("Value"));
1224   tab_text (chisq, 2, 0, TAB_RIGHT | TAT_TITLE, _("df"));
1225   tab_text (chisq, 3, 0, TAB_RIGHT | TAT_TITLE,
1226             _("Asymp. Sig. (2-sided)"));
1227   tab_text (chisq, 4, 0, TAB_RIGHT | TAT_TITLE,
1228             _("Exact. Sig. (2-sided)"));
1229   tab_text (chisq, 5, 0, TAB_RIGHT | TAT_TITLE,
1230             _("Exact. Sig. (1-sided)"));
1231   chisq = 0;
1232   tab_offset (chisq, 0, 1);
1233
1234   return chisq;
1235 }
1236
1237 /* Symmetric measures. */
1238 static struct tab_table *
1239 create_sym_table (struct pivot_table *pt)
1240 {
1241   struct tab_table *sym;
1242
1243   sym = tab_create (6 + (pt->n_vars - 2),
1244                     pt->n_entries / pt->n_cols * 7 + 10, 1);
1245   tab_headers (sym, 2 + (pt->n_vars - 2), 0, 1, 0);
1246   tab_title (sym, _("Symmetric measures."));
1247
1248   tab_offset (sym, pt->n_vars - 2, 0);
1249   tab_text (sym, 0, 0, TAB_LEFT | TAT_TITLE, _("Category"));
1250   tab_text (sym, 1, 0, TAB_LEFT | TAT_TITLE, _("Statistic"));
1251   tab_text (sym, 2, 0, TAB_RIGHT | TAT_TITLE, _("Value"));
1252   tab_text (sym, 3, 0, TAB_RIGHT | TAT_TITLE, _("Asymp. Std. Error"));
1253   tab_text (sym, 4, 0, TAB_RIGHT | TAT_TITLE, _("Approx. T"));
1254   tab_text (sym, 5, 0, TAB_RIGHT | TAT_TITLE, _("Approx. Sig."));
1255   tab_offset (sym, 0, 1);
1256
1257   return sym;
1258 }
1259
1260 /* Risk estimate. */
1261 static struct tab_table *
1262 create_risk_table (struct pivot_table *pt)
1263 {
1264   struct tab_table *risk;
1265
1266   risk = tab_create (4 + (pt->n_vars - 2), pt->n_entries / pt->n_cols * 4 + 10,
1267                      1);
1268   tab_headers (risk, 1 + pt->n_vars - 2, 0, 2, 0);
1269   tab_title (risk, _("Risk estimate."));
1270
1271   tab_offset (risk, pt->n_vars - 2, 0);
1272   tab_joint_text (risk, 2, 0, 3, 0, TAB_CENTER | TAT_TITLE | TAT_PRINTF,
1273                   _("95%% Confidence Interval"));
1274   tab_text (risk, 0, 1, TAB_LEFT | TAT_TITLE, _("Statistic"));
1275   tab_text (risk, 1, 1, TAB_RIGHT | TAT_TITLE, _("Value"));
1276   tab_text (risk, 2, 1, TAB_RIGHT | TAT_TITLE, _("Lower"));
1277   tab_text (risk, 3, 1, TAB_RIGHT | TAT_TITLE, _("Upper"));
1278   tab_hline (risk, TAL_1, 2, 3, 1);
1279   tab_vline (risk, TAL_1, 2, 0, 1);
1280   tab_offset (risk, 0, 2);
1281
1282   return risk;
1283 }
1284
1285 /* Directional measures. */
1286 static struct tab_table *
1287 create_direct_table (struct pivot_table *pt)
1288 {
1289   struct tab_table *direct;
1290
1291   direct = tab_create (7 + (pt->n_vars - 2),
1292                        pt->n_entries / pt->n_cols * 7 + 10, 1);
1293   tab_headers (direct, 3 + (pt->n_vars - 2), 0, 1, 0);
1294   tab_title (direct, _("Directional measures."));
1295
1296   tab_offset (direct, pt->n_vars - 2, 0);
1297   tab_text (direct, 0, 0, TAB_LEFT | TAT_TITLE, _("Category"));
1298   tab_text (direct, 1, 0, TAB_LEFT | TAT_TITLE, _("Statistic"));
1299   tab_text (direct, 2, 0, TAB_LEFT | TAT_TITLE, _("Type"));
1300   tab_text (direct, 3, 0, TAB_RIGHT | TAT_TITLE, _("Value"));
1301   tab_text (direct, 4, 0, TAB_RIGHT | TAT_TITLE, _("Asymp. Std. Error"));
1302   tab_text (direct, 5, 0, TAB_RIGHT | TAT_TITLE, _("Approx. T"));
1303   tab_text (direct, 6, 0, TAB_RIGHT | TAT_TITLE, _("Approx. Sig."));
1304   tab_offset (direct, 0, 1);
1305
1306   return direct;
1307 }
1308
1309
1310 /* Delete missing rows and columns for statistical analysis when
1311    /MISSING=REPORT. */
1312 static void
1313 delete_missing (struct pivot_table *pt)
1314 {
1315   int r, c;
1316
1317   for (r = 0; r < pt->n_rows; r++)
1318     if (var_is_num_missing (pt->vars[ROW_VAR], pt->rows[r].f, MV_USER))
1319       {
1320         for (c = 0; c < pt->n_cols; c++)
1321           pt->mat[c + r * pt->n_cols] = 0.;
1322         pt->ns_rows--;
1323       }
1324
1325
1326   for (c = 0; c < pt->n_cols; c++)
1327     if (var_is_num_missing (pt->vars[COL_VAR], pt->cols[c].f, MV_USER))
1328       {
1329         for (r = 0; r < pt->n_rows; r++)
1330           pt->mat[c + r * pt->n_cols] = 0.;
1331         pt->ns_cols--;
1332       }
1333 }
1334
1335 /* Prepare table T for submission, and submit it. */
1336 static void
1337 submit (struct crosstabs_proc *proc, struct pivot_table *pt,
1338         struct tab_table *t)
1339 {
1340   int i;
1341
1342   if (t == NULL)
1343     return;
1344
1345   tab_resize (t, -1, 0);
1346   if (tab_nr (t) == tab_t (t))
1347     {
1348       tab_destroy (t);
1349       return;
1350     }
1351   tab_offset (t, 0, 0);
1352   if (pt != NULL)
1353     for (i = 2; i < pt->n_vars; i++)
1354       tab_text (t, pt->n_vars - i - 1, 0, TAB_RIGHT | TAT_TITLE,
1355                 var_to_string (pt->vars[i]));
1356   tab_box (t, TAL_2, TAL_2, -1, -1, 0, 0, tab_nc (t) - 1, tab_nr (t) - 1);
1357   tab_box (t, -1, -1, -1, TAL_1, tab_l (t), tab_t (t) - 1, tab_nc (t) - 1,
1358            tab_nr (t) - 1);
1359   tab_box (t, -1, -1, -1, TAL_GAP, 0, tab_t (t), tab_l (t) - 1,
1360            tab_nr (t) - 1);
1361   tab_vline (t, TAL_2, tab_l (t), 0, tab_nr (t) - 1);
1362   tab_dim (t, crosstabs_dim, proc);
1363   tab_submit (t);
1364 }
1365
1366 /* Sets the widths of all the columns and heights of all the rows in
1367    table T for driver D. */
1368 static void
1369 crosstabs_dim (struct tab_table *t, struct outp_driver *d, void *proc_)
1370 {
1371   struct crosstabs_proc *proc = proc_;
1372   int i;
1373
1374   /* Width of a numerical column. */
1375   int c = outp_string_width (d, "0.000000", OUTP_PROPORTIONAL);
1376   if (proc->exclude == MV_NEVER)
1377     c += outp_string_width (d, "M", OUTP_PROPORTIONAL);
1378
1379   /* Set width for header columns. */
1380   if (t->l != 0)
1381     {
1382       size_t i;
1383       int w;
1384
1385       w = d->width - c * (t->nc - t->l);
1386       for (i = 0; i <= t->nc; i++)
1387         w -= t->wrv[i];
1388       w /= t->l;
1389
1390       if (w < d->prop_em_width * 8)
1391         w = d->prop_em_width * 8;
1392
1393       if (w > d->prop_em_width * 15)
1394         w = d->prop_em_width * 15;
1395
1396       for (i = 0; i < t->l; i++)
1397         t->w[i] = w;
1398     }
1399
1400   for (i = t->l; i < t->nc; i++)
1401     t->w[i] = c;
1402
1403   for (i = 0; i < t->nr; i++)
1404     t->h[i] = tab_natural_height (t, d, i);
1405 }
1406
1407 static bool
1408 find_crosstab (struct pivot_table *pt, size_t *row0p, size_t *row1p)
1409 {
1410   size_t row0 = *row1p;
1411   size_t row1;
1412
1413   if (row0 >= pt->n_entries)
1414     return false;
1415
1416   for (row1 = row0 + 1; row1 < pt->n_entries; row1++)
1417     {
1418       struct table_entry *a = pt->entries[row0];
1419       struct table_entry *b = pt->entries[row1];
1420       if (compare_table_entry_vars_3way (a, b, pt, 2, pt->n_vars) != 0)
1421         break;
1422     }
1423   *row0p = row0;
1424   *row1p = row1;
1425   return true;
1426 }
1427
1428 /* Compares `union value's A_ and B_ and returns a strcmp()-like
1429    result.  WIDTH_ points to an int which is either 0 for a
1430    numeric value or a string width for a string value. */
1431 static int
1432 compare_value_3way (const void *a_, const void *b_, const void *width_)
1433 {
1434   const union value *a = a_;
1435   const union value *b = b_;
1436   const int *width = width_;
1437
1438   return value_compare_3way (a, b, *width);
1439 }
1440
1441 /* Given an array of ENTRY_CNT table_entry structures starting at
1442    ENTRIES, creates a sorted list of the values that the variable
1443    with index VAR_IDX takes on.  The values are returned as a
1444    malloc()'d array stored in *VALUES, with the number of values
1445    stored in *VALUE_CNT.
1446    */
1447 static void
1448 enum_var_values (const struct pivot_table *pt, int var_idx,
1449                  union value **valuesp, int *n_values)
1450 {
1451   const struct variable *var = pt->vars[var_idx];
1452   struct var_range *range = get_var_range (var);
1453   union value *values;
1454   size_t i;
1455
1456   if (range)
1457     {
1458       values = *valuesp = xnmalloc (range->count, sizeof *values);
1459       *n_values = range->count;
1460       for (i = 0; i < range->count; i++)
1461         values[i].f = range->min + i;
1462     }
1463   else
1464     {
1465       int width = var_get_width (var);
1466       struct hmapx_node *node;
1467       const union value *iter;
1468       struct hmapx set;
1469
1470       hmapx_init (&set);
1471       for (i = 0; i < pt->n_entries; i++)
1472         {
1473           const struct table_entry *te = pt->entries[i];
1474           const union value *value = &te->values[var_idx];
1475           size_t hash = value_hash (value, width, 0);
1476
1477           HMAPX_FOR_EACH_WITH_HASH (iter, node, hash, &set)
1478             if (value_equal (iter, value, width))
1479               goto next_entry;
1480
1481           hmapx_insert (&set, (union value *) value, hash);
1482
1483         next_entry: ;
1484         }
1485
1486       *n_values = hmapx_count (&set);
1487       values = *valuesp = xnmalloc (*n_values, sizeof *values);
1488       i = 0;
1489       HMAPX_FOR_EACH (iter, node, &set)
1490         values[i++] = *iter;
1491       hmapx_destroy (&set);
1492
1493       sort (values, *n_values, sizeof *values, compare_value_3way, &width);
1494     }
1495 }
1496
1497 /* Sets cell (C,R) in TABLE, with options OPT, to have a value taken
1498    from V, displayed with print format spec from variable VAR.  When
1499    in REPORT missing-value mode, missing values have an M appended. */
1500 static void
1501 table_value_missing (struct crosstabs_proc *proc,
1502                      struct tab_table *table, int c, int r, unsigned char opt,
1503                      const union value *v, const struct variable *var)
1504 {
1505   struct substring s;
1506   const struct fmt_spec *print = var_get_print_format (var);
1507
1508   const char *label = var_lookup_value_label (var, v);
1509   if (label)
1510     {
1511       tab_text (table, c, r, TAB_LEFT, label);
1512       return;
1513     }
1514
1515   s.string = tab_alloc (table, print->w);
1516   data_out (v, print, s.string);
1517   s.length = print->w;
1518   if (proc->exclude == MV_NEVER && var_is_num_missing (var, v->f, MV_USER))
1519     s.string[s.length++] = 'M';
1520   while (s.length && *s.string == ' ')
1521     {
1522       s.length--;
1523       s.string++;
1524     }
1525   tab_raw (table, c, r, opt, &s);
1526 }
1527
1528 /* Draws a line across TABLE at the current row to indicate the most
1529    major dimension variable with index FIRST_DIFFERENCE out of N_VARS
1530    that changed, and puts the values that changed into the table.  TB
1531    and PT must be the corresponding table_entry and crosstab,
1532    respectively. */
1533 static void
1534 display_dimensions (struct crosstabs_proc *proc, struct pivot_table *pt,
1535                     struct tab_table *table, int first_difference)
1536 {
1537   tab_hline (table, TAL_1, pt->n_vars - first_difference - 1, tab_nc (table) - 1, 0);
1538
1539   for (; first_difference >= 2; first_difference--)
1540     table_value_missing (proc, table, pt->n_vars - first_difference - 1, 0,
1541                          TAB_RIGHT, &pt->entries[0]->values[first_difference],
1542                          pt->vars[first_difference]);
1543 }
1544
1545 /* Put VALUE into cell (C,R) of TABLE, suffixed with character
1546    SUFFIX if nonzero.  If MARK_MISSING is true the entry is
1547    additionally suffixed with a letter `M'. */
1548 static void
1549 format_cell_entry (struct tab_table *table, int c, int r, double value,
1550                    char suffix, bool mark_missing)
1551 {
1552   const struct fmt_spec f = {FMT_F, 10, 1};
1553   union value v;
1554   struct substring s;
1555
1556   s.length = 10;
1557   s.string = tab_alloc (table, 16);
1558   v.f = value;
1559   data_out (&v, &f, s.string);
1560   while (*s.string == ' ')
1561     {
1562       s.length--;
1563       s.string++;
1564     }
1565   if (suffix != 0)
1566     s.string[s.length++] = suffix;
1567   if (mark_missing)
1568     s.string[s.length++] = 'M';
1569
1570   tab_raw (table, c, r, TAB_RIGHT, &s);
1571 }
1572
1573 /* Displays the crosstabulation table. */
1574 static void
1575 display_crosstabulation (struct crosstabs_proc *proc, struct pivot_table *pt,
1576                          struct tab_table *table)
1577 {
1578   int last_row;
1579   int r, c, i;
1580   double *mp;
1581
1582   for (r = 0; r < pt->n_rows; r++)
1583     table_value_missing (proc, table, pt->n_vars - 2, r * proc->n_cells,
1584                          TAB_RIGHT, &pt->rows[r], pt->vars[ROW_VAR]);
1585
1586   tab_text (table, pt->n_vars - 2, pt->n_rows * proc->n_cells,
1587             TAB_LEFT, _("Total"));
1588
1589   /* Put in the actual cells. */
1590   mp = pt->mat;
1591   tab_offset (table, pt->n_vars - 1, -1);
1592   for (r = 0; r < pt->n_rows; r++)
1593     {
1594       if (proc->n_cells > 1)
1595         tab_hline (table, TAL_1, -1, pt->n_cols, 0);
1596       for (c = 0; c < pt->n_cols; c++)
1597         {
1598           bool mark_missing = false;
1599           double expected_value = pt->row_tot[r] * pt->col_tot[c] / pt->total;
1600           if (proc->exclude == MV_NEVER
1601               && (var_is_num_missing (pt->vars[COL_VAR], pt->cols[c].f, MV_USER)
1602                   || var_is_num_missing (pt->vars[ROW_VAR], pt->rows[r].f,
1603                                          MV_USER)))
1604             mark_missing = true;
1605           for (i = 0; i < proc->n_cells; i++)
1606             {
1607               double v;
1608               int suffix = 0;
1609
1610               switch (proc->a_cells[i])
1611                 {
1612                 case CRS_CL_COUNT:
1613                   v = *mp;
1614                   break;
1615                 case CRS_CL_ROW:
1616                   v = *mp / pt->row_tot[r] * 100.;
1617                   suffix = '%';
1618                   break;
1619                 case CRS_CL_COLUMN:
1620                   v = *mp / pt->col_tot[c] * 100.;
1621                   suffix = '%';
1622                   break;
1623                 case CRS_CL_TOTAL:
1624                   v = *mp / pt->total * 100.;
1625                   suffix = '%';
1626                   break;
1627                 case CRS_CL_EXPECTED:
1628                   v = expected_value;
1629                   break;
1630                 case CRS_CL_RESIDUAL:
1631                   v = *mp - expected_value;
1632                   break;
1633                 case CRS_CL_SRESIDUAL:
1634                   v = (*mp - expected_value) / sqrt (expected_value);
1635                   break;
1636                 case CRS_CL_ASRESIDUAL:
1637                   v = ((*mp - expected_value)
1638                        / sqrt (expected_value
1639                                * (1. - pt->row_tot[r] / pt->total)
1640                                * (1. - pt->col_tot[c] / pt->total)));
1641                   break;
1642                 default:
1643                   NOT_REACHED ();
1644                 }
1645               format_cell_entry (table, c, i, v, suffix, mark_missing);
1646             }
1647
1648           mp++;
1649         }
1650
1651       tab_offset (table, -1, tab_row (table) + proc->n_cells);
1652     }
1653
1654   /* Row totals. */
1655   tab_offset (table, -1, tab_row (table) - proc->n_cells * pt->n_rows);
1656   for (r = 0; r < pt->n_rows; r++)
1657     {
1658       bool mark_missing = false;
1659
1660       if (proc->exclude == MV_NEVER
1661           && var_is_num_missing (pt->vars[ROW_VAR], pt->rows[r].f, MV_USER))
1662         mark_missing = true;
1663
1664       for (i = 0; i < proc->n_cells; i++)
1665         {
1666           char suffix = 0;
1667           double v;
1668
1669           switch (proc->a_cells[i])
1670             {
1671             case CRS_CL_COUNT:
1672               v = pt->row_tot[r];
1673               break;
1674             case CRS_CL_ROW:
1675               v = 100.0;
1676               suffix = '%';
1677               break;
1678             case CRS_CL_COLUMN:
1679               v = pt->row_tot[r] / pt->total * 100.;
1680               suffix = '%';
1681               break;
1682             case CRS_CL_TOTAL:
1683               v = pt->row_tot[r] / pt->total * 100.;
1684               suffix = '%';
1685               break;
1686             case CRS_CL_EXPECTED:
1687             case CRS_CL_RESIDUAL:
1688             case CRS_CL_SRESIDUAL:
1689             case CRS_CL_ASRESIDUAL:
1690               v = 0.;
1691               break;
1692             default:
1693               NOT_REACHED ();
1694             }
1695
1696           format_cell_entry (table, pt->n_cols, 0, v, suffix, mark_missing);
1697           tab_next_row (table);
1698         }
1699     }
1700
1701   /* Column totals, grand total. */
1702   last_row = 0;
1703   if (proc->n_cells > 1)
1704     tab_hline (table, TAL_1, -1, pt->n_cols, 0);
1705   for (c = 0; c <= pt->n_cols; c++)
1706     {
1707       double ct = c < pt->n_cols ? pt->col_tot[c] : pt->total;
1708       bool mark_missing = false;
1709       int i;
1710
1711       if (proc->exclude == MV_NEVER && c < pt->n_cols
1712           && var_is_num_missing (pt->vars[COL_VAR], pt->cols[c].f, MV_USER))
1713         mark_missing = true;
1714
1715       for (i = 0; i < proc->n_cells; i++)
1716         {
1717           char suffix = 0;
1718           double v;
1719
1720           switch (proc->a_cells[i])
1721             {
1722             case CRS_CL_COUNT:
1723               v = ct;
1724               break;
1725             case CRS_CL_ROW:
1726               v = ct / pt->total * 100.;
1727               suffix = '%';
1728               break;
1729             case CRS_CL_COLUMN:
1730               v = 100.;
1731               suffix = '%';
1732               break;
1733             case CRS_CL_TOTAL:
1734               v = ct / pt->total * 100.;
1735               suffix = '%';
1736               break;
1737             case CRS_CL_EXPECTED:
1738             case CRS_CL_RESIDUAL:
1739             case CRS_CL_SRESIDUAL:
1740             case CRS_CL_ASRESIDUAL:
1741               continue;
1742             default:
1743               NOT_REACHED ();
1744             }
1745
1746           format_cell_entry (table, c, i, v, suffix, mark_missing);
1747         }
1748       last_row = i;
1749     }
1750
1751   tab_offset (table, -1, tab_row (table) + last_row);
1752   tab_offset (table, 0, -1);
1753 }
1754
1755 static void calc_r (struct pivot_table *,
1756                     double *PT, double *Y, double *, double *, double *);
1757 static void calc_chisq (struct pivot_table *,
1758                         double[N_CHISQ], int[N_CHISQ], double *, double *);
1759
1760 /* Display chi-square statistics. */
1761 static void
1762 display_chisq (struct pivot_table *pt, struct tab_table *chisq,
1763                bool *showed_fisher)
1764 {
1765   static const char *chisq_stats[N_CHISQ] =
1766     {
1767       N_("Pearson Chi-Square"),
1768       N_("Likelihood Ratio"),
1769       N_("Fisher's Exact Test"),
1770       N_("Continuity Correction"),
1771       N_("Linear-by-Linear Association"),
1772     };
1773   double chisq_v[N_CHISQ];
1774   double fisher1, fisher2;
1775   int df[N_CHISQ];
1776   int s = 0;
1777
1778   int i;
1779
1780   calc_chisq (pt, chisq_v, df, &fisher1, &fisher2);
1781
1782   tab_offset (chisq, pt->n_vars - 2, -1);
1783
1784   for (i = 0; i < N_CHISQ; i++)
1785     {
1786       if ((i != 2 && chisq_v[i] == SYSMIS)
1787           || (i == 2 && fisher1 == SYSMIS))
1788         continue;
1789       s = 1;
1790
1791       tab_text (chisq, 0, 0, TAB_LEFT, gettext (chisq_stats[i]));
1792       if (i != 2)
1793         {
1794           tab_double (chisq, 1, 0, TAB_RIGHT, chisq_v[i], NULL);
1795           tab_double (chisq, 2, 0, TAB_RIGHT, df[i], &pt->weight_format);
1796           tab_double (chisq, 3, 0, TAB_RIGHT,
1797                      gsl_cdf_chisq_Q (chisq_v[i], df[i]), NULL);
1798         }
1799       else
1800         {
1801           *showed_fisher = true;
1802           tab_double (chisq, 4, 0, TAB_RIGHT, fisher2, NULL);
1803           tab_double (chisq, 5, 0, TAB_RIGHT, fisher1, NULL);
1804         }
1805       tab_next_row (chisq);
1806     }
1807
1808   tab_text (chisq, 0, 0, TAB_LEFT, _("N of Valid Cases"));
1809   tab_double (chisq, 1, 0, TAB_RIGHT, pt->total, &pt->weight_format);
1810   tab_next_row (chisq);
1811
1812   tab_offset (chisq, 0, -1);
1813 }
1814
1815 static int calc_symmetric (struct crosstabs_proc *, struct pivot_table *,
1816                            double[N_SYMMETRIC], double[N_SYMMETRIC],
1817                            double[N_SYMMETRIC],
1818                            double[3], double[3], double[3]);
1819
1820 /* Display symmetric measures. */
1821 static void
1822 display_symmetric (struct crosstabs_proc *proc, struct pivot_table *pt,
1823                    struct tab_table *sym)
1824 {
1825   static const char *categories[] =
1826     {
1827       N_("Nominal by Nominal"),
1828       N_("Ordinal by Ordinal"),
1829       N_("Interval by Interval"),
1830       N_("Measure of Agreement"),
1831     };
1832
1833   static const char *stats[N_SYMMETRIC] =
1834     {
1835       N_("Phi"),
1836       N_("Cramer's V"),
1837       N_("Contingency Coefficient"),
1838       N_("Kendall's tau-b"),
1839       N_("Kendall's tau-c"),
1840       N_("Gamma"),
1841       N_("Spearman Correlation"),
1842       N_("Pearson's R"),
1843       N_("Kappa"),
1844     };
1845
1846   static const int stats_categories[N_SYMMETRIC] =
1847     {
1848       0, 0, 0, 1, 1, 1, 1, 2, 3,
1849     };
1850
1851   int last_cat = -1;
1852   double sym_v[N_SYMMETRIC], sym_ase[N_SYMMETRIC], sym_t[N_SYMMETRIC];
1853   double somers_d_v[3], somers_d_ase[3], somers_d_t[3];
1854   int i;
1855
1856   if (!calc_symmetric (proc, pt, sym_v, sym_ase, sym_t,
1857                        somers_d_v, somers_d_ase, somers_d_t))
1858     return;
1859
1860   tab_offset (sym, pt->n_vars - 2, -1);
1861
1862   for (i = 0; i < N_SYMMETRIC; i++)
1863     {
1864       if (sym_v[i] == SYSMIS)
1865         continue;
1866
1867       if (stats_categories[i] != last_cat)
1868         {
1869           last_cat = stats_categories[i];
1870           tab_text (sym, 0, 0, TAB_LEFT, gettext (categories[last_cat]));
1871         }
1872
1873       tab_text (sym, 1, 0, TAB_LEFT, gettext (stats[i]));
1874       tab_double (sym, 2, 0, TAB_RIGHT, sym_v[i], NULL);
1875       if (sym_ase[i] != SYSMIS)
1876         tab_double (sym, 3, 0, TAB_RIGHT, sym_ase[i], NULL);
1877       if (sym_t[i] != SYSMIS)
1878         tab_double (sym, 4, 0, TAB_RIGHT, sym_t[i], NULL);
1879       /*tab_double (sym, 5, 0, TAB_RIGHT, normal_sig (sym_v[i]), NULL);*/
1880       tab_next_row (sym);
1881     }
1882
1883   tab_text (sym, 0, 0, TAB_LEFT, _("N of Valid Cases"));
1884   tab_double (sym, 2, 0, TAB_RIGHT, pt->total, &pt->weight_format);
1885   tab_next_row (sym);
1886
1887   tab_offset (sym, 0, -1);
1888 }
1889
1890 static int calc_risk (struct pivot_table *,
1891                       double[], double[], double[], union value *);
1892
1893 /* Display risk estimate. */
1894 static void
1895 display_risk (struct pivot_table *pt, struct tab_table *risk)
1896 {
1897   char buf[256];
1898   double risk_v[3], lower[3], upper[3];
1899   union value c[2];
1900   int i;
1901
1902   if (!calc_risk (pt, risk_v, upper, lower, c))
1903     return;
1904
1905   tab_offset (risk, pt->n_vars - 2, -1);
1906
1907   for (i = 0; i < 3; i++)
1908     {
1909       const struct variable *cv = pt->vars[COL_VAR];
1910       const struct variable *rv = pt->vars[ROW_VAR];
1911       int cvw = var_get_width (cv);
1912       int rvw = var_get_width (rv);
1913
1914       if (risk_v[i] == SYSMIS)
1915         continue;
1916
1917       switch (i)
1918         {
1919         case 0:
1920           if (var_is_numeric (cv))
1921             sprintf (buf, _("Odds Ratio for %s (%g / %g)"),
1922                      var_get_name (cv), c[0].f, c[1].f);
1923           else
1924             sprintf (buf, _("Odds Ratio for %s (%.*s / %.*s)"),
1925                      var_get_name (cv),
1926                      cvw, value_str (&c[0], cvw),
1927                      cvw, value_str (&c[1], cvw));
1928           break;
1929         case 1:
1930         case 2:
1931           if (var_is_numeric (rv))
1932             sprintf (buf, _("For cohort %s = %g"),
1933                      var_get_name (rv), pt->rows[i - 1].f);
1934           else
1935             sprintf (buf, _("For cohort %s = %.*s"),
1936                      var_get_name (rv),
1937                      rvw, value_str (&pt->rows[i - 1], rvw));
1938           break;
1939         }
1940
1941       tab_text (risk, 0, 0, TAB_LEFT, buf);
1942       tab_double (risk, 1, 0, TAB_RIGHT, risk_v[i], NULL);
1943       tab_double (risk, 2, 0, TAB_RIGHT, lower[i], NULL);
1944       tab_double (risk, 3, 0, TAB_RIGHT, upper[i], NULL);
1945       tab_next_row (risk);
1946     }
1947
1948   tab_text (risk, 0, 0, TAB_LEFT, _("N of Valid Cases"));
1949   tab_double (risk, 1, 0, TAB_RIGHT, pt->total, &pt->weight_format);
1950   tab_next_row (risk);
1951
1952   tab_offset (risk, 0, -1);
1953 }
1954
1955 static int calc_directional (struct crosstabs_proc *, struct pivot_table *,
1956                              double[N_DIRECTIONAL], double[N_DIRECTIONAL],
1957                              double[N_DIRECTIONAL]);
1958
1959 /* Display directional measures. */
1960 static void
1961 display_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
1962                      struct tab_table *direct)
1963 {
1964   static const char *categories[] =
1965     {
1966       N_("Nominal by Nominal"),
1967       N_("Ordinal by Ordinal"),
1968       N_("Nominal by Interval"),
1969     };
1970
1971   static const char *stats[] =
1972     {
1973       N_("Lambda"),
1974       N_("Goodman and Kruskal tau"),
1975       N_("Uncertainty Coefficient"),
1976       N_("Somers' d"),
1977       N_("Eta"),
1978     };
1979
1980   static const char *types[] =
1981     {
1982       N_("Symmetric"),
1983       N_("%s Dependent"),
1984       N_("%s Dependent"),
1985     };
1986
1987   static const int stats_categories[N_DIRECTIONAL] =
1988     {
1989       0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 2,
1990     };
1991
1992   static const int stats_stats[N_DIRECTIONAL] =
1993     {
1994       0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4,
1995     };
1996
1997   static const int stats_types[N_DIRECTIONAL] =
1998     {
1999       0, 1, 2, 1, 2, 0, 1, 2, 0, 1, 2, 1, 2,
2000     };
2001
2002   static const int *stats_lookup[] =
2003     {
2004       stats_categories,
2005       stats_stats,
2006       stats_types,
2007     };
2008
2009   static const char **stats_names[] =
2010     {
2011       categories,
2012       stats,
2013       types,
2014     };
2015
2016   int last[3] =
2017     {
2018       -1, -1, -1,
2019     };
2020
2021   double direct_v[N_DIRECTIONAL];
2022   double direct_ase[N_DIRECTIONAL];
2023   double direct_t[N_DIRECTIONAL];
2024
2025   int i;
2026
2027   if (!calc_directional (proc, pt, direct_v, direct_ase, direct_t))
2028     return;
2029
2030   tab_offset (direct, pt->n_vars - 2, -1);
2031
2032   for (i = 0; i < N_DIRECTIONAL; i++)
2033     {
2034       if (direct_v[i] == SYSMIS)
2035         continue;
2036
2037       {
2038         int j;
2039
2040         for (j = 0; j < 3; j++)
2041           if (last[j] != stats_lookup[j][i])
2042             {
2043               if (j < 2)
2044                 tab_hline (direct, TAL_1, j, 6, 0);
2045
2046               for (; j < 3; j++)
2047                 {
2048                   const char *string;
2049                   int k = last[j] = stats_lookup[j][i];
2050
2051                   if (k == 0)
2052                     string = NULL;
2053                   else if (k == 1)
2054                     string = var_get_name (pt->vars[0]);
2055                   else
2056                     string = var_get_name (pt->vars[1]);
2057
2058                   tab_text (direct, j, 0, TAB_LEFT | TAT_PRINTF,
2059                             gettext (stats_names[j][k]), string);
2060                 }
2061             }
2062       }
2063
2064       tab_double (direct, 3, 0, TAB_RIGHT, direct_v[i], NULL);
2065       if (direct_ase[i] != SYSMIS)
2066         tab_double (direct, 4, 0, TAB_RIGHT, direct_ase[i], NULL);
2067       if (direct_t[i] != SYSMIS)
2068         tab_double (direct, 5, 0, TAB_RIGHT, direct_t[i], NULL);
2069       /*tab_double (direct, 6, 0, TAB_RIGHT, normal_sig (direct_v[i]), NULL);*/
2070       tab_next_row (direct);
2071     }
2072
2073   tab_offset (direct, 0, -1);
2074 }
2075 \f
2076 /* Statistical calculations. */
2077
2078 /* Returns the value of the gamma (factorial) function for an integer
2079    argument PT. */
2080 static double
2081 gamma_int (double pt)
2082 {
2083   double r = 1;
2084   int i;
2085
2086   for (i = 2; i < pt; i++)
2087     r *= i;
2088   return r;
2089 }
2090
2091 /* Calculate P_r as specified in _SPSS Statistical Algorithms_,
2092    Appendix 5. */
2093 static inline double
2094 Pr (int a, int b, int c, int d)
2095 {
2096   return (gamma_int (a + b + 1.) / gamma_int (a + 1.)
2097           * gamma_int (c + d + 1.) / gamma_int (b + 1.)
2098           * gamma_int (a + c + 1.) / gamma_int (c + 1.)
2099           * gamma_int (b + d + 1.) / gamma_int (d + 1.)
2100           / gamma_int (a + b + c + d + 1.));
2101 }
2102
2103 /* Swap the contents of A and B. */
2104 static inline void
2105 swap (int *a, int *b)
2106 {
2107   int t = *a;
2108   *a = *b;
2109   *b = t;
2110 }
2111
2112 /* Calculate significance for Fisher's exact test as specified in
2113    _SPSS Statistical Algorithms_, Appendix 5. */
2114 static void
2115 calc_fisher (int a, int b, int c, int d, double *fisher1, double *fisher2)
2116 {
2117   int pt;
2118
2119   if (MIN (c, d) < MIN (a, b))
2120     swap (&a, &c), swap (&b, &d);
2121   if (MIN (b, d) < MIN (a, c))
2122     swap (&a, &b), swap (&c, &d);
2123   if (b * c < a * d)
2124     {
2125       if (b < c)
2126         swap (&a, &b), swap (&c, &d);
2127       else
2128         swap (&a, &c), swap (&b, &d);
2129     }
2130
2131   *fisher1 = 0.;
2132   for (pt = 0; pt <= a; pt++)
2133     *fisher1 += Pr (a - pt, b + pt, c + pt, d - pt);
2134
2135   *fisher2 = *fisher1;
2136   for (pt = 1; pt <= b; pt++)
2137     *fisher2 += Pr (a + pt, b - pt, c - pt, d + pt);
2138 }
2139
2140 /* Calculates chi-squares into CHISQ.  MAT is a matrix with N_COLS
2141    columns with values COLS and N_ROWS rows with values ROWS.  Values
2142    in the matrix sum to pt->total. */
2143 static void
2144 calc_chisq (struct pivot_table *pt,
2145             double chisq[N_CHISQ], int df[N_CHISQ],
2146             double *fisher1, double *fisher2)
2147 {
2148   int r, c;
2149
2150   chisq[0] = chisq[1] = 0.;
2151   chisq[2] = chisq[3] = chisq[4] = SYSMIS;
2152   *fisher1 = *fisher2 = SYSMIS;
2153
2154   df[0] = df[1] = (pt->ns_cols - 1) * (pt->ns_rows - 1);
2155
2156   if (pt->ns_rows <= 1 || pt->ns_cols <= 1)
2157     {
2158       chisq[0] = chisq[1] = SYSMIS;
2159       return;
2160     }
2161
2162   for (r = 0; r < pt->n_rows; r++)
2163     for (c = 0; c < pt->n_cols; c++)
2164       {
2165         const double expected = pt->row_tot[r] * pt->col_tot[c] / pt->total;
2166         const double freq = pt->mat[pt->n_cols * r + c];
2167         const double residual = freq - expected;
2168
2169         chisq[0] += residual * residual / expected;
2170         if (freq)
2171           chisq[1] += freq * log (expected / freq);
2172       }
2173
2174   if (chisq[0] == 0.)
2175     chisq[0] = SYSMIS;
2176
2177   if (chisq[1] != 0.)
2178     chisq[1] *= -2.;
2179   else
2180     chisq[1] = SYSMIS;
2181
2182   /* Calculate Yates and Fisher exact test. */
2183   if (pt->ns_cols == 2 && pt->ns_rows == 2)
2184     {
2185       double f11, f12, f21, f22;
2186
2187       {
2188         int nz_cols[2];
2189         int i, j;
2190
2191         for (i = j = 0; i < pt->n_cols; i++)
2192           if (pt->col_tot[i] != 0.)
2193             {
2194               nz_cols[j++] = i;
2195               if (j == 2)
2196                 break;
2197             }
2198
2199         assert (j == 2);
2200
2201         f11 = pt->mat[nz_cols[0]];
2202         f12 = pt->mat[nz_cols[1]];
2203         f21 = pt->mat[nz_cols[0] + pt->n_cols];
2204         f22 = pt->mat[nz_cols[1] + pt->n_cols];
2205       }
2206
2207       /* Yates. */
2208       {
2209         const double pt_ = fabs (f11 * f22 - f12 * f21) - 0.5 * pt->total;
2210
2211         if (pt_ > 0.)
2212           chisq[3] = (pt->total * pow2 (pt_)
2213                       / (f11 + f12) / (f21 + f22)
2214                       / (f11 + f21) / (f12 + f22));
2215         else
2216           chisq[3] = 0.;
2217
2218         df[3] = 1.;
2219       }
2220
2221       /* Fisher. */
2222       if (f11 < 5. || f12 < 5. || f21 < 5. || f22 < 5.)
2223         calc_fisher (f11 + .5, f12 + .5, f21 + .5, f22 + .5, fisher1, fisher2);
2224     }
2225
2226   /* Calculate Mantel-Haenszel. */
2227   if (var_is_numeric (pt->vars[ROW_VAR]) && var_is_numeric (pt->vars[COL_VAR]))
2228     {
2229       double r, ase_0, ase_1;
2230       calc_r (pt, (double *) pt->rows, (double *) pt->cols, &r, &ase_0, &ase_1);
2231
2232       chisq[4] = (pt->total - 1.) * r * r;
2233       df[4] = 1;
2234     }
2235 }
2236
2237 /* Calculate the value of Pearson's r.  r is stored into R, ase_1 into
2238    ASE_1, and ase_0 into ASE_0.  The row and column values must be
2239    passed in PT and Y. */
2240 static void
2241 calc_r (struct pivot_table *pt,
2242         double *PT, double *Y, double *r, double *ase_0, double *ase_1)
2243 {
2244   double SX, SY, S, T;
2245   double Xbar, Ybar;
2246   double sum_XYf, sum_X2Y2f;
2247   double sum_Xr, sum_X2r;
2248   double sum_Yc, sum_Y2c;
2249   int i, j;
2250
2251   for (sum_X2Y2f = sum_XYf = 0., i = 0; i < pt->n_rows; i++)
2252     for (j = 0; j < pt->n_cols; j++)
2253       {
2254         double fij = pt->mat[j + i * pt->n_cols];
2255         double product = PT[i] * Y[j];
2256         double temp = fij * product;
2257         sum_XYf += temp;
2258         sum_X2Y2f += temp * product;
2259       }
2260
2261   for (sum_Xr = sum_X2r = 0., i = 0; i < pt->n_rows; i++)
2262     {
2263       sum_Xr += PT[i] * pt->row_tot[i];
2264       sum_X2r += pow2 (PT[i]) * pt->row_tot[i];
2265     }
2266   Xbar = sum_Xr / pt->total;
2267
2268   for (sum_Yc = sum_Y2c = 0., i = 0; i < pt->n_cols; i++)
2269     {
2270       sum_Yc += Y[i] * pt->col_tot[i];
2271       sum_Y2c += Y[i] * Y[i] * pt->col_tot[i];
2272     }
2273   Ybar = sum_Yc / pt->total;
2274
2275   S = sum_XYf - sum_Xr * sum_Yc / pt->total;
2276   SX = sum_X2r - pow2 (sum_Xr) / pt->total;
2277   SY = sum_Y2c - pow2 (sum_Yc) / pt->total;
2278   T = sqrt (SX * SY);
2279   *r = S / T;
2280   *ase_0 = sqrt ((sum_X2Y2f - pow2 (sum_XYf) / pt->total) / (sum_X2r * sum_Y2c));
2281
2282   {
2283     double s, c, y, t;
2284
2285     for (s = c = 0., i = 0; i < pt->n_rows; i++)
2286       for (j = 0; j < pt->n_cols; j++)
2287         {
2288           double Xresid, Yresid;
2289           double temp;
2290
2291           Xresid = PT[i] - Xbar;
2292           Yresid = Y[j] - Ybar;
2293           temp = (T * Xresid * Yresid
2294                   - ((S / (2. * T))
2295                      * (Xresid * Xresid * SY + Yresid * Yresid * SX)));
2296           y = pt->mat[j + i * pt->n_cols] * temp * temp - c;
2297           t = s + y;
2298           c = (t - s) - y;
2299           s = t;
2300         }
2301     *ase_1 = sqrt (s) / (T * T);
2302   }
2303 }
2304
2305 /* Calculate symmetric statistics and their asymptotic standard
2306    errors.  Returns 0 if none could be calculated. */
2307 static int
2308 calc_symmetric (struct crosstabs_proc *proc, struct pivot_table *pt,
2309                 double v[N_SYMMETRIC], double ase[N_SYMMETRIC],
2310                 double t[N_SYMMETRIC],
2311                 double somers_d_v[3], double somers_d_ase[3],
2312                 double somers_d_t[3])
2313 {
2314   int q, i;
2315
2316   q = MIN (pt->ns_rows, pt->ns_cols);
2317   if (q <= 1)
2318     return 0;
2319
2320   for (i = 0; i < N_SYMMETRIC; i++)
2321     v[i] = ase[i] = t[i] = SYSMIS;
2322
2323   /* Phi, Cramer's V, contingency coefficient. */
2324   if (proc->statistics & ((1u << CRS_ST_PHI) | (1u << CRS_ST_CC)))
2325     {
2326       double Xp = 0.;   /* Pearson chi-square. */
2327       int r, c;
2328
2329       for (r = 0; r < pt->n_rows; r++)
2330         for (c = 0; c < pt->n_cols; c++)
2331           {
2332             const double expected = pt->row_tot[r] * pt->col_tot[c] / pt->total;
2333             const double freq = pt->mat[pt->n_cols * r + c];
2334             const double residual = freq - expected;
2335
2336             Xp += residual * residual / expected;
2337           }
2338
2339       if (proc->statistics & (1u << CRS_ST_PHI))
2340         {
2341           v[0] = sqrt (Xp / pt->total);
2342           v[1] = sqrt (Xp / (pt->total * (q - 1)));
2343         }
2344       if (proc->statistics & (1u << CRS_ST_CC))
2345         v[2] = sqrt (Xp / (Xp + pt->total));
2346     }
2347
2348   if (proc->statistics & ((1u << CRS_ST_BTAU) | (1u << CRS_ST_CTAU)
2349                           | (1u << CRS_ST_GAMMA) | (1u << CRS_ST_D)))
2350     {
2351       double *cum;
2352       double Dr, Dc;
2353       double P, Q;
2354       double btau_cum, ctau_cum, gamma_cum, d_yx_cum, d_xy_cum;
2355       double btau_var;
2356       int r, c;
2357
2358       Dr = Dc = pow2 (pt->total);
2359       for (r = 0; r < pt->n_rows; r++)
2360         Dr -= pow2 (pt->row_tot[r]);
2361       for (c = 0; c < pt->n_cols; c++)
2362         Dc -= pow2 (pt->col_tot[c]);
2363
2364       cum = xnmalloc (pt->n_cols * pt->n_rows, sizeof *cum);
2365       for (c = 0; c < pt->n_cols; c++)
2366         {
2367           double ct = 0.;
2368
2369           for (r = 0; r < pt->n_rows; r++)
2370             cum[c + r * pt->n_cols] = ct += pt->mat[c + r * pt->n_cols];
2371         }
2372
2373       /* P and Q. */
2374       {
2375         int i, j;
2376         double Cij, Dij;
2377
2378         P = Q = 0.;
2379         for (i = 0; i < pt->n_rows; i++)
2380           {
2381             Cij = Dij = 0.;
2382
2383             for (j = 1; j < pt->n_cols; j++)
2384               Cij += pt->col_tot[j] - cum[j + i * pt->n_cols];
2385
2386             if (i > 0)
2387               for (j = 1; j < pt->n_cols; j++)
2388                 Dij += cum[j + (i - 1) * pt->n_cols];
2389
2390             for (j = 0;;)
2391               {
2392                 double fij = pt->mat[j + i * pt->n_cols];
2393                 P += fij * Cij;
2394                 Q += fij * Dij;
2395
2396                 if (++j == pt->n_cols)
2397                   break;
2398                 assert (j < pt->n_cols);
2399
2400                 Cij -= pt->col_tot[j] - cum[j + i * pt->n_cols];
2401                 Dij += pt->col_tot[j - 1] - cum[j - 1 + i * pt->n_cols];
2402
2403                 if (i > 0)
2404                   {
2405                     Cij += cum[j - 1 + (i - 1) * pt->n_cols];
2406                     Dij -= cum[j + (i - 1) * pt->n_cols];
2407                   }
2408               }
2409           }
2410       }
2411
2412       if (proc->statistics & (1u << CRS_ST_BTAU))
2413         v[3] = (P - Q) / sqrt (Dr * Dc);
2414       if (proc->statistics & (1u << CRS_ST_CTAU))
2415         v[4] = (q * (P - Q)) / (pow2 (pt->total) * (q - 1));
2416       if (proc->statistics & (1u << CRS_ST_GAMMA))
2417         v[5] = (P - Q) / (P + Q);
2418
2419       /* ASE for tau-b, tau-c, gamma.  Calculations could be
2420          eliminated here, at expense of memory.  */
2421       {
2422         int i, j;
2423         double Cij, Dij;
2424
2425         btau_cum = ctau_cum = gamma_cum = d_yx_cum = d_xy_cum = 0.;
2426         for (i = 0; i < pt->n_rows; i++)
2427           {
2428             Cij = Dij = 0.;
2429
2430             for (j = 1; j < pt->n_cols; j++)
2431               Cij += pt->col_tot[j] - cum[j + i * pt->n_cols];
2432
2433             if (i > 0)
2434               for (j = 1; j < pt->n_cols; j++)
2435                 Dij += cum[j + (i - 1) * pt->n_cols];
2436
2437             for (j = 0;;)
2438               {
2439                 double fij = pt->mat[j + i * pt->n_cols];
2440
2441                 if (proc->statistics & (1u << CRS_ST_BTAU))
2442                   {
2443                     const double temp = (2. * sqrt (Dr * Dc) * (Cij - Dij)
2444                                          + v[3] * (pt->row_tot[i] * Dc
2445                                                    + pt->col_tot[j] * Dr));
2446                     btau_cum += fij * temp * temp;
2447                   }
2448
2449                 {
2450                   const double temp = Cij - Dij;
2451                   ctau_cum += fij * temp * temp;
2452                 }
2453
2454                 if (proc->statistics & (1u << CRS_ST_GAMMA))
2455                   {
2456                     const double temp = Q * Cij - P * Dij;
2457                     gamma_cum += fij * temp * temp;
2458                   }
2459
2460                 if (proc->statistics & (1u << CRS_ST_D))
2461                   {
2462                     d_yx_cum += fij * pow2 (Dr * (Cij - Dij)
2463                                             - (P - Q) * (pt->total - pt->row_tot[i]));
2464                     d_xy_cum += fij * pow2 (Dc * (Dij - Cij)
2465                                             - (Q - P) * (pt->total - pt->col_tot[j]));
2466                   }
2467
2468                 if (++j == pt->n_cols)
2469                   break;
2470                 assert (j < pt->n_cols);
2471
2472                 Cij -= pt->col_tot[j] - cum[j + i * pt->n_cols];
2473                 Dij += pt->col_tot[j - 1] - cum[j - 1 + i * pt->n_cols];
2474
2475                 if (i > 0)
2476                   {
2477                     Cij += cum[j - 1 + (i - 1) * pt->n_cols];
2478                     Dij -= cum[j + (i - 1) * pt->n_cols];
2479                   }
2480               }
2481           }
2482       }
2483
2484       btau_var = ((btau_cum
2485                    - (pt->total * pow2 (pt->total * (P - Q) / sqrt (Dr * Dc) * (Dr + Dc))))
2486                   / pow2 (Dr * Dc));
2487       if (proc->statistics & (1u << CRS_ST_BTAU))
2488         {
2489           ase[3] = sqrt (btau_var);
2490           t[3] = v[3] / (2 * sqrt ((ctau_cum - (P - Q) * (P - Q) / pt->total)
2491                                    / (Dr * Dc)));
2492         }
2493       if (proc->statistics & (1u << CRS_ST_CTAU))
2494         {
2495           ase[4] = ((2 * q / ((q - 1) * pow2 (pt->total)))
2496                     * sqrt (ctau_cum - (P - Q) * (P - Q) / pt->total));
2497           t[4] = v[4] / ase[4];
2498         }
2499       if (proc->statistics & (1u << CRS_ST_GAMMA))
2500         {
2501           ase[5] = ((4. / ((P + Q) * (P + Q))) * sqrt (gamma_cum));
2502           t[5] = v[5] / (2. / (P + Q)
2503                          * sqrt (ctau_cum - (P - Q) * (P - Q) / pt->total));
2504         }
2505       if (proc->statistics & (1u << CRS_ST_D))
2506         {
2507           somers_d_v[0] = (P - Q) / (.5 * (Dc + Dr));
2508           somers_d_ase[0] = 2. * btau_var / (Dr + Dc) * sqrt (Dr * Dc);
2509           somers_d_t[0] = (somers_d_v[0]
2510                            / (4 / (Dc + Dr)
2511                               * sqrt (ctau_cum - pow2 (P - Q) / pt->total)));
2512           somers_d_v[1] = (P - Q) / Dc;
2513           somers_d_ase[1] = 2. / pow2 (Dc) * sqrt (d_xy_cum);
2514           somers_d_t[1] = (somers_d_v[1]
2515                            / (2. / Dc
2516                               * sqrt (ctau_cum - pow2 (P - Q) / pt->total)));
2517           somers_d_v[2] = (P - Q) / Dr;
2518           somers_d_ase[2] = 2. / pow2 (Dr) * sqrt (d_yx_cum);
2519           somers_d_t[2] = (somers_d_v[2]
2520                            / (2. / Dr
2521                               * sqrt (ctau_cum - pow2 (P - Q) / pt->total)));
2522         }
2523
2524       free (cum);
2525     }
2526
2527   /* Spearman correlation, Pearson's r. */
2528   if (proc->statistics & (1u << CRS_ST_CORR))
2529     {
2530       double *R = xmalloc (sizeof *R * pt->n_rows);
2531       double *C = xmalloc (sizeof *C * pt->n_cols);
2532
2533       {
2534         double y, t, c = 0., s = 0.;
2535         int i = 0;
2536
2537         for (;;)
2538           {
2539             R[i] = s + (pt->row_tot[i] + 1.) / 2.;
2540             y = pt->row_tot[i] - c;
2541             t = s + y;
2542             c = (t - s) - y;
2543             s = t;
2544             if (++i == pt->n_rows)
2545               break;
2546             assert (i < pt->n_rows);
2547           }
2548       }
2549
2550       {
2551         double y, t, c = 0., s = 0.;
2552         int j = 0;
2553
2554         for (;;)
2555           {
2556             C[j] = s + (pt->col_tot[j] + 1.) / 2;
2557             y = pt->col_tot[j] - c;
2558             t = s + y;
2559             c = (t - s) - y;
2560             s = t;
2561             if (++j == pt->n_cols)
2562               break;
2563             assert (j < pt->n_cols);
2564           }
2565       }
2566
2567       calc_r (pt, R, C, &v[6], &t[6], &ase[6]);
2568       t[6] = v[6] / t[6];
2569
2570       free (R);
2571       free (C);
2572
2573       calc_r (pt, (double *) pt->rows, (double *) pt->cols, &v[7], &t[7], &ase[7]);
2574       t[7] = v[7] / t[7];
2575     }
2576
2577   /* Cohen's kappa. */
2578   if (proc->statistics & (1u << CRS_ST_KAPPA) && pt->ns_rows == pt->ns_cols)
2579     {
2580       double sum_fii, sum_rici, sum_fiiri_ci, sum_fijri_ci2, sum_riciri_ci;
2581       int i, j;
2582
2583       for (sum_fii = sum_rici = sum_fiiri_ci = sum_riciri_ci = 0., i = j = 0;
2584            i < pt->ns_rows; i++, j++)
2585         {
2586           double prod, sum;
2587
2588           while (pt->col_tot[j] == 0.)
2589             j++;
2590
2591           prod = pt->row_tot[i] * pt->col_tot[j];
2592           sum = pt->row_tot[i] + pt->col_tot[j];
2593
2594           sum_fii += pt->mat[j + i * pt->n_cols];
2595           sum_rici += prod;
2596           sum_fiiri_ci += pt->mat[j + i * pt->n_cols] * sum;
2597           sum_riciri_ci += prod * sum;
2598         }
2599       for (sum_fijri_ci2 = 0., i = 0; i < pt->ns_rows; i++)
2600         for (j = 0; j < pt->ns_cols; j++)
2601           {
2602             double sum = pt->row_tot[i] + pt->col_tot[j];
2603             sum_fijri_ci2 += pt->mat[j + i * pt->n_cols] * sum * sum;
2604           }
2605
2606       v[8] = (pt->total * sum_fii - sum_rici) / (pow2 (pt->total) - sum_rici);
2607
2608       ase[8] = sqrt ((pow2 (pt->total) * sum_rici
2609                       + sum_rici * sum_rici
2610                       - pt->total * sum_riciri_ci)
2611                      / (pt->total * (pow2 (pt->total) - sum_rici) * (pow2 (pt->total) - sum_rici)));
2612 #if 0
2613       t[8] = v[8] / sqrt (pt->total * (((sum_fii * (pt->total - sum_fii))
2614                                 / pow2 (pow2 (pt->total) - sum_rici))
2615                                + ((2. * (pt->total - sum_fii)
2616                                    * (2. * sum_fii * sum_rici
2617                                       - pt->total * sum_fiiri_ci))
2618                                   / cube (pow2 (pt->total) - sum_rici))
2619                                + (pow2 (pt->total - sum_fii)
2620                                   * (pt->total * sum_fijri_ci2 - 4.
2621                                      * sum_rici * sum_rici)
2622                                   / pow4 (pow2 (pt->total) - sum_rici))));
2623 #else
2624       t[8] = v[8] / ase[8];
2625 #endif
2626     }
2627
2628   return 1;
2629 }
2630
2631 /* Calculate risk estimate. */
2632 static int
2633 calc_risk (struct pivot_table *pt,
2634            double *value, double *upper, double *lower, union value *c)
2635 {
2636   double f11, f12, f21, f22;
2637   double v;
2638
2639   {
2640     int i;
2641
2642     for (i = 0; i < 3; i++)
2643       value[i] = upper[i] = lower[i] = SYSMIS;
2644   }
2645
2646   if (pt->ns_rows != 2 || pt->ns_cols != 2)
2647     return 0;
2648
2649   {
2650     int nz_cols[2];
2651     int i, j;
2652
2653     for (i = j = 0; i < pt->n_cols; i++)
2654       if (pt->col_tot[i] != 0.)
2655         {
2656           nz_cols[j++] = i;
2657           if (j == 2)
2658             break;
2659         }
2660
2661     assert (j == 2);
2662
2663     f11 = pt->mat[nz_cols[0]];
2664     f12 = pt->mat[nz_cols[1]];
2665     f21 = pt->mat[nz_cols[0] + pt->n_cols];
2666     f22 = pt->mat[nz_cols[1] + pt->n_cols];
2667
2668     c[0] = pt->cols[nz_cols[0]];
2669     c[1] = pt->cols[nz_cols[1]];
2670   }
2671
2672   value[0] = (f11 * f22) / (f12 * f21);
2673   v = sqrt (1. / f11 + 1. / f12 + 1. / f21 + 1. / f22);
2674   lower[0] = value[0] * exp (-1.960 * v);
2675   upper[0] = value[0] * exp (1.960 * v);
2676
2677   value[1] = (f11 * (f21 + f22)) / (f21 * (f11 + f12));
2678   v = sqrt ((f12 / (f11 * (f11 + f12)))
2679             + (f22 / (f21 * (f21 + f22))));
2680   lower[1] = value[1] * exp (-1.960 * v);
2681   upper[1] = value[1] * exp (1.960 * v);
2682
2683   value[2] = (f12 * (f21 + f22)) / (f22 * (f11 + f12));
2684   v = sqrt ((f11 / (f12 * (f11 + f12)))
2685             + (f21 / (f22 * (f21 + f22))));
2686   lower[2] = value[2] * exp (-1.960 * v);
2687   upper[2] = value[2] * exp (1.960 * v);
2688
2689   return 1;
2690 }
2691
2692 /* Calculate directional measures. */
2693 static int
2694 calc_directional (struct crosstabs_proc *proc, struct pivot_table *pt,
2695                   double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
2696                   double t[N_DIRECTIONAL])
2697 {
2698   {
2699     int i;
2700
2701     for (i = 0; i < N_DIRECTIONAL; i++)
2702       v[i] = ase[i] = t[i] = SYSMIS;
2703   }
2704
2705   /* Lambda. */
2706   if (proc->statistics & (1u << CRS_ST_LAMBDA))
2707     {
2708       double *fim = xnmalloc (pt->n_rows, sizeof *fim);
2709       int *fim_index = xnmalloc (pt->n_rows, sizeof *fim_index);
2710       double *fmj = xnmalloc (pt->n_cols, sizeof *fmj);
2711       int *fmj_index = xnmalloc (pt->n_cols, sizeof *fmj_index);
2712       double sum_fim, sum_fmj;
2713       double rm, cm;
2714       int rm_index, cm_index;
2715       int i, j;
2716
2717       /* Find maximum for each row and their sum. */
2718       for (sum_fim = 0., i = 0; i < pt->n_rows; i++)
2719         {
2720           double max = pt->mat[i * pt->n_cols];
2721           int index = 0;
2722
2723           for (j = 1; j < pt->n_cols; j++)
2724             if (pt->mat[j + i * pt->n_cols] > max)
2725               {
2726                 max = pt->mat[j + i * pt->n_cols];
2727                 index = j;
2728               }
2729
2730           sum_fim += fim[i] = max;
2731           fim_index[i] = index;
2732         }
2733
2734       /* Find maximum for each column. */
2735       for (sum_fmj = 0., j = 0; j < pt->n_cols; j++)
2736         {
2737           double max = pt->mat[j];
2738           int index = 0;
2739
2740           for (i = 1; i < pt->n_rows; i++)
2741             if (pt->mat[j + i * pt->n_cols] > max)
2742               {
2743                 max = pt->mat[j + i * pt->n_cols];
2744                 index = i;
2745               }
2746
2747           sum_fmj += fmj[j] = max;
2748           fmj_index[j] = index;
2749         }
2750
2751       /* Find maximum row total. */
2752       rm = pt->row_tot[0];
2753       rm_index = 0;
2754       for (i = 1; i < pt->n_rows; i++)
2755         if (pt->row_tot[i] > rm)
2756           {
2757             rm = pt->row_tot[i];
2758             rm_index = i;
2759           }
2760
2761       /* Find maximum column total. */
2762       cm = pt->col_tot[0];
2763       cm_index = 0;
2764       for (j = 1; j < pt->n_cols; j++)
2765         if (pt->col_tot[j] > cm)
2766           {
2767             cm = pt->col_tot[j];
2768             cm_index = j;
2769           }
2770
2771       v[0] = (sum_fim + sum_fmj - cm - rm) / (2. * pt->total - rm - cm);
2772       v[1] = (sum_fmj - rm) / (pt->total - rm);
2773       v[2] = (sum_fim - cm) / (pt->total - cm);
2774
2775       /* ASE1 for Y given PT. */
2776       {
2777         double accum;
2778
2779         for (accum = 0., i = 0; i < pt->n_rows; i++)
2780           for (j = 0; j < pt->n_cols; j++)
2781             {
2782               const int deltaj = j == cm_index;
2783               accum += (pt->mat[j + i * pt->n_cols]
2784                         * pow2 ((j == fim_index[i])
2785                                - deltaj
2786                                + v[0] * deltaj));
2787             }
2788
2789         ase[2] = sqrt (accum - pt->total * v[0]) / (pt->total - cm);
2790       }
2791
2792       /* ASE0 for Y given PT. */
2793       {
2794         double accum;
2795
2796         for (accum = 0., i = 0; i < pt->n_rows; i++)
2797           if (cm_index != fim_index[i])
2798             accum += (pt->mat[i * pt->n_cols + fim_index[i]]
2799                       + pt->mat[i * pt->n_cols + cm_index]);
2800         t[2] = v[2] / (sqrt (accum - pow2 (sum_fim - cm) / pt->total) / (pt->total - cm));
2801       }
2802
2803       /* ASE1 for PT given Y. */
2804       {
2805         double accum;
2806
2807         for (accum = 0., i = 0; i < pt->n_rows; i++)
2808           for (j = 0; j < pt->n_cols; j++)
2809             {
2810               const int deltaj = i == rm_index;
2811               accum += (pt->mat[j + i * pt->n_cols]
2812                         * pow2 ((i == fmj_index[j])
2813                                - deltaj
2814                                + v[0] * deltaj));
2815             }
2816
2817         ase[1] = sqrt (accum - pt->total * v[0]) / (pt->total - rm);
2818       }
2819
2820       /* ASE0 for PT given Y. */
2821       {
2822         double accum;
2823
2824         for (accum = 0., j = 0; j < pt->n_cols; j++)
2825           if (rm_index != fmj_index[j])
2826             accum += (pt->mat[j + pt->n_cols * fmj_index[j]]
2827                       + pt->mat[j + pt->n_cols * rm_index]);
2828         t[1] = v[1] / (sqrt (accum - pow2 (sum_fmj - rm) / pt->total) / (pt->total - rm));
2829       }
2830
2831       /* Symmetric ASE0 and ASE1. */
2832       {
2833         double accum0;
2834         double accum1;
2835
2836         for (accum0 = accum1 = 0., i = 0; i < pt->n_rows; i++)
2837           for (j = 0; j < pt->n_cols; j++)
2838             {
2839               int temp0 = (fmj_index[j] == i) + (fim_index[i] == j);
2840               int temp1 = (i == rm_index) + (j == cm_index);
2841               accum0 += pt->mat[j + i * pt->n_cols] * pow2 (temp0 - temp1);
2842               accum1 += (pt->mat[j + i * pt->n_cols]
2843                          * pow2 (temp0 + (v[0] - 1.) * temp1));
2844             }
2845         ase[0] = sqrt (accum1 - 4. * pt->total * v[0] * v[0]) / (2. * pt->total - rm - cm);
2846         t[0] = v[0] / (sqrt (accum0 - pow2 ((sum_fim + sum_fmj - cm - rm) / pt->total))
2847                        / (2. * pt->total - rm - cm));
2848       }
2849
2850       free (fim);
2851       free (fim_index);
2852       free (fmj);
2853       free (fmj_index);
2854
2855       {
2856         double sum_fij2_ri, sum_fij2_ci;
2857         double sum_ri2, sum_cj2;
2858
2859         for (sum_fij2_ri = sum_fij2_ci = 0., i = 0; i < pt->n_rows; i++)
2860           for (j = 0; j < pt->n_cols; j++)
2861             {
2862               double temp = pow2 (pt->mat[j + i * pt->n_cols]);
2863               sum_fij2_ri += temp / pt->row_tot[i];
2864               sum_fij2_ci += temp / pt->col_tot[j];
2865             }
2866
2867         for (sum_ri2 = 0., i = 0; i < pt->n_rows; i++)
2868           sum_ri2 += pow2 (pt->row_tot[i]);
2869
2870         for (sum_cj2 = 0., j = 0; j < pt->n_cols; j++)
2871           sum_cj2 += pow2 (pt->col_tot[j]);
2872
2873         v[3] = (pt->total * sum_fij2_ci - sum_ri2) / (pow2 (pt->total) - sum_ri2);
2874         v[4] = (pt->total * sum_fij2_ri - sum_cj2) / (pow2 (pt->total) - sum_cj2);
2875       }
2876     }
2877
2878   if (proc->statistics & (1u << CRS_ST_UC))
2879     {
2880       double UX, UY, UXY, P;
2881       double ase1_yx, ase1_xy, ase1_sym;
2882       int i, j;
2883
2884       for (UX = 0., i = 0; i < pt->n_rows; i++)
2885         if (pt->row_tot[i] > 0.)
2886           UX -= pt->row_tot[i] / pt->total * log (pt->row_tot[i] / pt->total);
2887
2888       for (UY = 0., j = 0; j < pt->n_cols; j++)
2889         if (pt->col_tot[j] > 0.)
2890           UY -= pt->col_tot[j] / pt->total * log (pt->col_tot[j] / pt->total);
2891
2892       for (UXY = P = 0., i = 0; i < pt->n_rows; i++)
2893         for (j = 0; j < pt->n_cols; j++)
2894           {
2895             double entry = pt->mat[j + i * pt->n_cols];
2896
2897             if (entry <= 0.)
2898               continue;
2899
2900             P += entry * pow2 (log (pt->col_tot[j] * pt->row_tot[i] / (pt->total * entry)));
2901             UXY -= entry / pt->total * log (entry / pt->total);
2902           }
2903
2904       for (ase1_yx = ase1_xy = ase1_sym = 0., i = 0; i < pt->n_rows; i++)
2905         for (j = 0; j < pt->n_cols; j++)
2906           {
2907             double entry = pt->mat[j + i * pt->n_cols];
2908
2909             if (entry <= 0.)
2910               continue;
2911
2912             ase1_yx += entry * pow2 (UY * log (entry / pt->row_tot[i])
2913                                     + (UX - UXY) * log (pt->col_tot[j] / pt->total));
2914             ase1_xy += entry * pow2 (UX * log (entry / pt->col_tot[j])
2915                                     + (UY - UXY) * log (pt->row_tot[i] / pt->total));
2916             ase1_sym += entry * pow2 ((UXY
2917                                       * log (pt->row_tot[i] * pt->col_tot[j] / pow2 (pt->total)))
2918                                      - (UX + UY) * log (entry / pt->total));
2919           }
2920
2921       v[5] = 2. * ((UX + UY - UXY) / (UX + UY));
2922       ase[5] = (2. / (pt->total * pow2 (UX + UY))) * sqrt (ase1_sym);
2923       t[5] = v[5] / ((2. / (pt->total * (UX + UY)))
2924                      * sqrt (P - pow2 (UX + UY - UXY) / pt->total));
2925
2926       v[6] = (UX + UY - UXY) / UX;
2927       ase[6] = sqrt (ase1_xy) / (pt->total * UX * UX);
2928       t[6] = v[6] / (sqrt (P - pt->total * pow2 (UX + UY - UXY)) / (pt->total * UX));
2929
2930       v[7] = (UX + UY - UXY) / UY;
2931       ase[7] = sqrt (ase1_yx) / (pt->total * UY * UY);
2932       t[7] = v[7] / (sqrt (P - pt->total * pow2 (UX + UY - UXY)) / (pt->total * UY));
2933     }
2934
2935   /* Somers' D. */
2936   if (proc->statistics & (1u << CRS_ST_D))
2937     {
2938       double v_dummy[N_SYMMETRIC];
2939       double ase_dummy[N_SYMMETRIC];
2940       double t_dummy[N_SYMMETRIC];
2941       double somers_d_v[3];
2942       double somers_d_ase[3];
2943       double somers_d_t[3];
2944
2945       if (calc_symmetric (proc, pt, v_dummy, ase_dummy, t_dummy,
2946                           somers_d_v, somers_d_ase, somers_d_t))
2947         {
2948           int i;
2949           for (i = 0; i < 3; i++)
2950             {
2951               v[8 + i] = somers_d_v[i];
2952               ase[8 + i] = somers_d_ase[i];
2953               t[8 + i] = somers_d_t[i];
2954             }
2955         }
2956     }
2957
2958   /* Eta. */
2959   if (proc->statistics & (1u << CRS_ST_ETA))
2960     {
2961       {
2962         double sum_Xr, sum_X2r;
2963         double SX, SXW;
2964         int i, j;
2965
2966         for (sum_Xr = sum_X2r = 0., i = 0; i < pt->n_rows; i++)
2967           {
2968             sum_Xr += pt->rows[i].f * pt->row_tot[i];
2969             sum_X2r += pow2 (pt->rows[i].f) * pt->row_tot[i];
2970           }
2971         SX = sum_X2r - pow2 (sum_Xr) / pt->total;
2972
2973         for (SXW = 0., j = 0; j < pt->n_cols; j++)
2974           {
2975             double cum;
2976
2977             for (cum = 0., i = 0; i < pt->n_rows; i++)
2978               {
2979                 SXW += pow2 (pt->rows[i].f) * pt->mat[j + i * pt->n_cols];
2980                 cum += pt->rows[i].f * pt->mat[j + i * pt->n_cols];
2981               }
2982
2983             SXW -= cum * cum / pt->col_tot[j];
2984           }
2985         v[11] = sqrt (1. - SXW / SX);
2986       }
2987
2988       {
2989         double sum_Yc, sum_Y2c;
2990         double SY, SYW;
2991         int i, j;
2992
2993         for (sum_Yc = sum_Y2c = 0., i = 0; i < pt->n_cols; i++)
2994           {
2995             sum_Yc += pt->cols[i].f * pt->col_tot[i];
2996             sum_Y2c += pow2 (pt->cols[i].f) * pt->col_tot[i];
2997           }
2998         SY = sum_Y2c - sum_Yc * sum_Yc / pt->total;
2999
3000         for (SYW = 0., i = 0; i < pt->n_rows; i++)
3001           {
3002             double cum;
3003
3004             for (cum = 0., j = 0; j < pt->n_cols; j++)
3005               {
3006                 SYW += pow2 (pt->cols[j].f) * pt->mat[j + i * pt->n_cols];
3007                 cum += pt->cols[j].f * pt->mat[j + i * pt->n_cols];
3008               }
3009
3010             SYW -= cum * cum / pt->row_tot[i];
3011           }
3012         v[12] = sqrt (1. - SYW / SY);
3013       }
3014     }
3015
3016   return 1;
3017 }
3018
3019 /*
3020    Local Variables:
3021    mode: c
3022    End:
3023 */