src/crosstabs.q

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