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