src/language/stats/crosstabs.q

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