src/language/stats/crosstabs.q

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