pintos-os.org Git - pspp/blob - src/language/stats/crosstabs.q

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