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