pintos-os.org Git - pspp-builds.git/blob

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2008, 2009, 2010 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 /*
  18   Create and update the values in the covariance matrix.
  19 */
  20 #include <assert.h>
  21 #include <config.h>
  22 #include <data/case.h>
  23 #include <data/category.h>
  24 #include <data/variable.h>
  25 #include <data/value.h>
  26 #include <libpspp/hash.h>
  27 #include <libpspp/hash-functions.h>
  28 #include <math/covariance-matrix.h>
  29 #include <math/moments.h>
  30 #include <string.h>
  31 #include <xalloc.h>
  32
  33 /*
  34   Structure used to accumulate the covariance matrix in a single data
  35   pass.  Before passing the data, we do not know how many categories
  36   there are in each categorical variable. Therefore we do not know the
  37   size of the covariance matrix. To get around this problem, we
  38   accumulate the elements of the covariance matrix in pointers to
  39   COVARIANC_ACCUMULATOR. These values are then used to populate
  40   the covariance matrix.
  41  */
  42 struct covariance_accumulator
  43 {
  44   const struct variable *v1;
  45   const struct variable *v2;
  46   const union value *val1;
  47   const union value *val2;
  48   double dot_product;
  49   double sum1;
  50   double sum2;
  51   double ssize;
  52 };
  53
  54
  55
  56 struct covariance_matrix
  57 {
  58   struct design_matrix *cov;
  59   struct design_matrix *ssize;
  60   struct hsh_table *ca;
  61   struct moments1 **m1;
  62   struct moments **m;
  63   const struct variable **v_variables;
  64   const struct interaction_variable **interactions;
  65   size_t n_variables;
  66   size_t n_intr;
  67   int n_pass;
  68   int missing_handling;
  69   enum mv_class missing_value;
  70   void (*accumulate) (struct covariance_matrix *, const struct ccase *,
  71                       const struct interaction_variable **, size_t);
  72   void (*update_moments) (struct covariance_matrix *, size_t, double);
  73 };
  74
  75
  76
  77 static struct hsh_table *covariance_hsh_create (size_t *);
  78 static hsh_hash_func covariance_accumulator_hash;
  79 static unsigned int hash_numeric_alpha (const struct variable *,
  80                                         const struct variable *,
  81                                         const union value *, size_t);
  82 static hsh_compare_func covariance_accumulator_compare;
  83 static hsh_free_func covariance_accumulator_free;
  84 static void update_moments1 (struct covariance_matrix *, size_t, double);
  85 static void update_moments2 (struct covariance_matrix *, size_t, double);
  86 static struct covariance_accumulator *get_new_covariance_accumulator (const
  87                                                                       struct
  88                                                                       variable
  89                                                                       *,
  90                                                                       const
  91                                                                       struct
  92                                                                       variable
  93                                                                       *,
  94                                                                       const
  95                                                                       union
  96                                                                       value *,
  97                                                                       const
  98                                                                       union
  99                                                                       value
 100                                                                       *);
 101 static void covariance_accumulate_listwise (struct covariance_matrix *,
 102                                             const struct ccase *,
 103                                             const struct interaction_variable **,
 104                                             size_t);
 105 static void covariance_accumulate_pairwise (struct covariance_matrix *,
 106                                             const struct ccase *,
 107                                             const struct interaction_variable **,
 108                                             size_t);
 109
 110 struct covariance_matrix *
 111 covariance_matrix_init (size_t n_variables,
 112                         const struct variable *v_variables[], int n_pass,
 113                         int missing_handling, enum mv_class missing_value)
 114 {
 115   size_t i;
 116   struct covariance_matrix *result = NULL;
 117
 118   result = xmalloc (sizeof (*result));
 119   result->cov = NULL;
 120   result->n_variables = n_variables;
 121   result->ca = covariance_hsh_create (&result->n_variables);
 122   result->m = NULL;
 123   result->m1 = NULL;
 124   result->n_intr = 0;
 125   result->missing_handling = missing_handling;
 126   result->missing_value = missing_value;
 127   result->accumulate = (result->missing_handling == LISTWISE) ?
 128     covariance_accumulate_listwise : covariance_accumulate_pairwise;
 129   if (n_pass == ONE_PASS)
 130     {
 131       result->update_moments = update_moments1;
 132       result->m1 = xnmalloc (n_variables, sizeof (*result->m1));
 133       for (i = 0; i < n_variables; i++)
 134         {
 135           result->m1[i] = moments1_create (MOMENT_MEAN);
 136         }
 137     }
 138   else
 139     {
 140       result->update_moments = update_moments2;
 141       result->m = xnmalloc (n_variables, sizeof (*result->m));
 142       for (i = 0; i < n_variables; i++)
 143         {
 144           result->m[i] = moments_create (MOMENT_MEAN);
 145         }
 146     }
 147   result->v_variables = v_variables;
 148
 149   result->n_pass = n_pass;
 150
 151   return result;
 152 }
 153 void
 154 covariance_interaction_set (struct covariance_matrix *cov,
 155                             const struct interaction_variable **intr, size_t n_intr)
 156 {
 157   cov->interactions = intr;
 158   cov->n_intr = n_intr;
 159 }
 160
 161 static size_t
 162 get_n_rows (size_t n_variables, const struct variable *v_variables[])
 163 {
 164   size_t i;
 165   size_t result = 0;
 166   for (i = 0; i < n_variables; i++)
 167     {
 168       if (var_is_numeric (v_variables[i]))
 169         {
 170           result++;
 171         }
 172       else if (var_is_alpha (v_variables[i]))
 173         {
 174           size_t n_categories = cat_get_n_categories (v_variables[i]);
 175           result += n_categories - 1;
 176         }
 177     }
 178   return result;
 179 }
 180 /*
 181   The covariances are stored in a DESIGN_MATRIX structure.
 182  */
 183 struct design_matrix *
 184 covariance_matrix_create (size_t n_variables,
 185                           const struct variable *v_variables[])
 186 {
 187   size_t n_rows = get_n_rows (n_variables, v_variables);
 188   return design_matrix_create (n_variables, v_variables, n_rows);
 189 }
 190
 191 static size_t
 192 get_n_rows_s (const struct variable *var)
 193 {
 194   size_t result = 0;
 195   if (var_is_numeric (var))
 196     {
 197       result++;
 198     }
 199   else
 200     {
 201       result += cat_get_n_categories (var) - 1;
 202     }
 203   return result;
 204 }
 205 static struct design_matrix *
 206 covariance_matrix_create_s (struct covariance_matrix *cov)
 207 {
 208   struct variable **v_variables;
 209   size_t n_variables;
 210   size_t n_rows = 0;
 211   size_t i;
 212   size_t j;
 213
 214   n_variables = cov->n_variables + cov->n_intr;
 215   v_variables = xnmalloc (n_variables, sizeof (*v_variables));
 216   for (i = 0; i < cov->n_variables; i++)
 217     {
 218       v_variables[i] = cov->v_variables[i];
 219       n_rows += get_n_rows_s (v_variables[i]);
 220     }
 221   for (j = 0; j < cov->n_intr; j++)
 222     {
 223       v_variables[i + j] = interaction_get_variable (cov->interactions[j]);
 224       n_rows += get_n_rows_s (v_variables[i]);
 225     }
 226   return design_matrix_create (n_variables, v_variables, n_rows);
 227 }
 228
 229 static void
 230 update_moments1 (struct covariance_matrix *cov, size_t i, double x)
 231 {
 232   assert (cov->m1 != NULL);
 233   moments1_add (cov->m1[i], x, 1.0);
 234 }
 235
 236 static void
 237 update_moments2 (struct covariance_matrix *cov, size_t i, double x)
 238 {
 239   assert (cov->m != NULL);
 240   moments_pass_one (cov->m[i], x, 1.0);
 241 }
 242
 243 void
 244 covariance_matrix_destroy (struct covariance_matrix *cov)
 245 {
 246   size_t i;
 247
 248   assert (cov != NULL);
 249   design_matrix_destroy (cov->cov);
 250   design_matrix_destroy (cov->ssize);
 251   hsh_destroy (cov->ca);
 252   if (cov->n_pass == ONE_PASS)
 253     {
 254       for (i = 0; i < cov->n_variables; i++)
 255         {
 256           moments1_destroy (cov->m1[i]);
 257         }
 258       free (cov->m1);
 259     }
 260   else
 261     {
 262       for (i = 0; i < cov->n_variables; i++)
 263         {
 264           moments_destroy (cov->m[i]);
 265         }
 266       free (cov->m);
 267     }
 268 }
 269
 270 /*
 271   Update the covariance matrix with the new entries, assuming that ROW
 272   corresponds to a categorical variable and V2 is numeric.
 273  */
 274 static void
 275 covariance_update_categorical_numeric (struct design_matrix *cov, double mean,
 276                                        size_t row,
 277                                        const struct variable *v2, double x,
 278                                        const union value *val2)
 279 {
 280   size_t col;
 281   double tmp;
 282
 283   assert (var_is_numeric (v2));
 284
 285   col = design_matrix_var_to_column (cov, v2);
 286   assert (val2 != NULL);
 287   tmp = design_matrix_get_element (cov, row, col);
 288   design_matrix_set_element (cov, row, col, (val2->f - mean) * x + tmp);
 289   design_matrix_set_element (cov, col, row, (val2->f - mean) * x + tmp);
 290 }
 291 static void
 292 column_iterate (struct design_matrix *cov, const struct variable *v,
 293                 double ssize, double x, const union value *val1, size_t row)
 294 {
 295   int width = var_get_width (v);
 296   size_t col;
 297   size_t i;
 298   double y;
 299   double tmp;
 300   const union value *tmp_val;
 301
 302   col = design_matrix_var_to_column (cov, v);
 303   for (i = 0; i < cat_get_n_categories (v) - 1; i++)
 304     {
 305       col += i;
 306       y = -1.0 * cat_get_category_count (i, v) / ssize;
 307       tmp_val = cat_subscript_to_value (i, v);
 308       if (!value_equal (tmp_val, val1, width))
 309         {
 310           y += -1.0;
 311         }
 312       tmp = design_matrix_get_element (cov, row, col);
 313       design_matrix_set_element (cov, row, col, x * y + tmp);
 314       design_matrix_set_element (cov, col, row, x * y + tmp);
 315     }
 316 }
 317
 318 /*
 319   Call this function in the second data pass. The central moments are
 320   MEAN1 and MEAN2. Any categorical variables should already have their
 321   values summarized in in its OBS_VALS element.
 322  */
 323 void
 324 covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
 325                      double ssize, const struct variable *v1,
 326                      const struct variable *v2, const union value *val1,
 327                      const union value *val2)
 328 {
 329   size_t row;
 330   size_t col;
 331   size_t i;
 332   double x;
 333   const union value *tmp_val;
 334
 335   if (var_is_alpha (v1))
 336     {
 337       row = design_matrix_var_to_column (cov, v1);
 338       for (i = 0; i < cat_get_n_categories (v1) - 1; i++)
 339         {
 340           row += i;
 341           x = -1.0 * cat_get_category_count (i, v1) / ssize;
 342           tmp_val = cat_subscript_to_value (i, v1);
 343           if (!value_equal (tmp_val, val1, var_get_width (v1)))
 344             {
 345               x += 1.0;
 346             }
 347           if (var_is_numeric (v2))
 348             {
 349               covariance_update_categorical_numeric (cov, mean2, row,
 350                                                      v2, x, val2);
 351             }
 352           else
 353             {
 354               column_iterate (cov, v1, ssize, x, val1, row);
 355               column_iterate (cov, v2, ssize, x, val2, row);
 356             }
 357         }
 358     }
 359   else if (var_is_alpha (v2))
 360     {
 361       /*
 362          Reverse the orders of V1, V2, etc. and put ourselves back
 363          in the previous IF scope.
 364        */
 365       covariance_pass_two (cov, mean2, mean1, ssize, v2, v1, val2, val1);
 366     }
 367   else
 368     {
 369       /*
 370          Both variables are numeric.
 371        */
 372       row = design_matrix_var_to_column (cov, v1);
 373       col = design_matrix_var_to_column (cov, v2);
 374       x = (val1->f - mean1) * (val2->f - mean2);
 375       x += design_matrix_get_element (cov, col, row);
 376       design_matrix_set_element (cov, row, col, x);
 377       design_matrix_set_element (cov, col, row, x);
 378     }
 379 }
 380
 381 static unsigned int
 382 covariance_accumulator_hash (const void *h, const void *aux)
 383 {
 384   struct covariance_accumulator *ca = (struct covariance_accumulator *) h;
 385   size_t *n_vars = (size_t *) aux;
 386   const struct variable *v_min;
 387   const struct variable *v_max;
 388   const union value *val_min;
 389   const union value *val_max;
 390
 391   /*
 392      Order everything by the variables' addresses. This ensures we get the
 393      same key regardless of the order in which the variables are stored
 394      and passed around.
 395    */
 396   if (ca->v1 < ca->v2)
 397     {
 398       v_min = ca->v1;
 399       v_max = ca->v2;
 400       val_min = ca->val1;
 401       val_max = ca->val2;
 402     }
 403   else
 404     {
 405       v_min = ca->v2;
 406       v_max = ca->v2;
 407       val_min = ca->val2;
 408       val_max = ca->val1;
 409     }
 410
 411   if (var_is_numeric (v_max) && var_is_numeric (v_min))
 412     {
 413       return hash_pointer (v_min, hash_pointer (v_max, 0));
 414     }
 415   if (var_is_numeric (v_max) && var_is_alpha (v_min))
 416     {
 417       return hash_numeric_alpha (v_max, v_min, val_min, *n_vars);
 418     }
 419   if (var_is_alpha (v_max) && var_is_numeric (v_min))
 420     {
 421       return (hash_numeric_alpha (v_min, v_max, val_max, *n_vars));
 422     }
 423   if (var_is_alpha (v_max) && var_is_alpha (v_min))
 424     {
 425       unsigned hash = value_hash (val_max, var_get_width (v_max), 0);
 426       hash = value_hash (val_min, var_get_width (v_min), hash);
 427       hash = hash_pointer (v_min, hash);
 428       return hash_pointer (v_max, hash);
 429     }
 430   return -1u;
 431 }
 432
 433 /*
 434   Make a hash table consisting of struct covariance_accumulators.
 435   This allows the accumulation of the elements of a covariance matrix
 436   in a single data pass. Call covariance_accumulate () for each case
 437   in the data.
 438  */
 439 static struct hsh_table *
 440 covariance_hsh_create (size_t *n_vars)
 441 {
 442   return hsh_create (*n_vars * *n_vars, covariance_accumulator_compare,
 443                      covariance_accumulator_hash, covariance_accumulator_free,
 444                      n_vars);
 445 }
 446
 447 static void
 448 covariance_accumulator_free (void *c_, const void *aux UNUSED)
 449 {
 450   struct covariance_accumulator *c = c_;
 451   assert (c != NULL);
 452   free (c);
 453 }
 454
 455 static int
 456 ordered_match_nodes (const struct covariance_accumulator *c, const struct variable *v1,
 457                      const struct variable *v2, const union value *val1, const union value *val2)
 458 {
 459   return (v1 != c->v1
 460           || v2 != c->v2
 461           || (var_is_alpha (v1)
 462               && !value_equal (val1, c->val1, var_get_width (v1)))
 463           || (var_is_alpha (v2)
 464               && !value_equal (val2, c->val2, var_get_width (v2))));
 465 }
 466
 467 /*
 468   Hash comparison. Returns 0 for a match, or a non-zero int
 469   otherwise. The sign of a non-zero return value *should* indicate the
 470   position of C relative to the covariance_accumulator described by
 471   the other arguments. But for now, it just returns 1 for any
 472   non-match.  This should be changed when someone figures out how to
 473   compute a sensible sign for the return value.
 474  */
 475 static int
 476 match_nodes (const struct covariance_accumulator *c,
 477              const struct variable *v1, const struct variable *v2,
 478              const union value *val1, const union value *val2)
 479 {
 480   return (ordered_match_nodes (c, v1, v2, val1, val2)
 481           && ordered_match_nodes (c, v2, v1, val2, val1));
 482 }
 483
 484 /*
 485   This function is meant to be used as a comparison function for
 486   a struct hsh_table in src/libpspp/hash.c.
 487 */
 488 static int
 489 covariance_accumulator_compare (const void *a1_, const void *a2_,
 490                                 const void *aux UNUSED)
 491 {
 492   const struct covariance_accumulator *a1 = a1_;
 493   const struct covariance_accumulator *a2 = a2_;
 494
 495   if (a1 == NULL && a2 == NULL)
 496     return 0;
 497
 498   if (a1 == NULL || a2 == NULL)
 499     return 1;
 500
 501   return match_nodes (a1, a2->v1, a2->v2, a2->val1, a2->val2);
 502 }
 503
 504 static unsigned int
 505 hash_numeric_alpha (const struct variable *v1, const struct variable *v2,
 506                     const union value *val, size_t n_vars)
 507 {
 508   unsigned int result = -1u;
 509   if (var_is_numeric (v1) && var_is_alpha (v2))
 510     {
 511       result = hash_pointer (v1, 0);
 512       result = hash_pointer (v2, result);
 513       result = value_hash (val, var_get_width (v2), result);
 514     }
 515   else if (var_is_alpha (v1) && var_is_numeric (v2))
 516     {
 517       result = hash_numeric_alpha (v2, v1, val, n_vars);
 518     }
 519   return result;
 520 }
 521
 522
 523 static double
 524 update_product (const struct variable *v1, const struct variable *v2,
 525                 const union value *val1, const union value *val2)
 526 {
 527   assert (v1 != NULL);
 528   assert (v2 != NULL);
 529   assert (val1 != NULL);
 530   assert (val2 != NULL);
 531   if (var_is_alpha (v1) && var_is_alpha (v2))
 532     {
 533       return 1.0;
 534     }
 535   if (var_is_numeric (v1) && var_is_numeric (v2))
 536     {
 537       return (val1->f * val2->f);
 538     }
 539   if (var_is_numeric (v1) && var_is_alpha (v2))
 540     {
 541       return val1->f;
 542     }
 543   if (var_is_numeric (v2) && var_is_alpha (v1))
 544     {
 545       return val2->f;
 546     }
 547   else
 548     {
 549       return 0.0;
 550     }
 551 }
 552 static double
 553 update_sum (const struct variable *var, const union value *val, double weight)
 554 {
 555   assert (var != NULL);
 556   assert (val != NULL);
 557   if (var_is_alpha (var))
 558     {
 559       return weight;
 560     }
 561   return val->f;
 562 }
 563 static struct covariance_accumulator *
 564 get_new_covariance_accumulator (const struct variable *v1,
 565                                 const struct variable *v2,
 566                                 const union value *val1,
 567                                 const union value *val2)
 568 {
 569   if ((v1 != NULL) && (v2 != NULL) && (val1 != NULL) && (val2 != NULL))
 570     {
 571       struct covariance_accumulator *ca;
 572       ca = xmalloc (sizeof (*ca));
 573       ca->v1 = v1;
 574       ca->v2 = v2;
 575       ca->val1 = val1;
 576       ca->val2 = val2;
 577       return ca;
 578     }
 579   return NULL;
 580 }
 581
 582 static const struct variable **
 583 get_covariance_variables (const struct covariance_matrix *cov)
 584 {
 585   return cov->v_variables;
 586 }
 587
 588 static void
 589 update_hash_entry_intr (struct hsh_table *c,
 590                         const struct variable *v1,
 591                         const struct variable *v2,
 592                         const union value *val1, const union value *val2,
 593                         const struct interaction_value *i_val1,
 594                         const struct interaction_value *i_val2)
 595 {
 596   struct covariance_accumulator *ca;
 597   struct covariance_accumulator *new_entry;
 598   double iv_f1;
 599   double iv_f2;
 600
 601   iv_f1 = interaction_value_get_nonzero_entry (i_val1);
 602   iv_f2 = interaction_value_get_nonzero_entry (i_val2);
 603   ca = get_new_covariance_accumulator (v1, v2, val1, val2);
 604   ca->dot_product = update_product (ca->v1, ca->v2, ca->val1, ca->val2);
 605   ca->dot_product *= iv_f1 * iv_f2;
 606   ca->sum1 = update_sum (ca->v1, ca->val1, iv_f1);
 607   ca->sum2 = update_sum (ca->v2, ca->val2, iv_f2);
 608   ca->ssize = 1.0;
 609   new_entry = hsh_insert (c, ca);
 610
 611   if (new_entry != NULL)
 612     {
 613       new_entry->dot_product += ca->dot_product;
 614       new_entry->ssize += 1.0;
 615       new_entry->sum1 += ca->sum1;
 616       new_entry->sum2 += ca->sum2;
 617       /*
 618         If DOT_PRODUCT is null, CA was not already in the hash
 619         hable, so we don't free it because it was just inserted.
 620         If DOT_PRODUCT was not null, CA is already in the hash table.
 621         Unnecessary now, it must be freed here.
 622       */
 623       free (ca);
 624     }
 625 }
 626
 627 static void
 628 update_hash_entry (struct hsh_table *c,
 629                    const struct variable *v1,
 630                    const struct variable *v2,
 631                    const union value *val1, const union value *val2)
 632 {
 633   struct covariance_accumulator *ca;
 634   struct covariance_accumulator *new_entry;
 635
 636   ca = get_new_covariance_accumulator (v1, v2, val1, val2);
 637   ca->dot_product = update_product (ca->v1, ca->v2, ca->val1, ca->val2);
 638   ca->sum1 = update_sum (ca->v1, ca->val1, 1.0);
 639   ca->sum2 = update_sum (ca->v2, ca->val2, 1.0);
 640   ca->ssize = 1.0;
 641   new_entry = hsh_insert (c, ca);
 642
 643   if (new_entry != NULL)
 644     {
 645       new_entry->dot_product += ca->dot_product;
 646       new_entry->ssize += 1.0;
 647       new_entry->sum1 += ca->sum1;
 648       new_entry->sum2 += ca->sum2;
 649       /*
 650         If DOT_PRODUCT is null, CA was not already in the hash
 651         hable, so we don't free it because it was just inserted.
 652         If DOT_PRODUCT was not null, CA is already in the hash table.
 653         Unnecessary now, it must be freed here.
 654       */
 655       free (ca);
 656     }
 657 }
 658
 659 static void
 660 inner_intr_loop (struct covariance_matrix *cov, const struct ccase  *ccase, const struct variable *var1,
 661                  const union value *val1, const struct interaction_variable **i_var,
 662                  const struct interaction_value *i_val1, size_t j)
 663 {
 664   struct variable *var2;
 665   union value *val2;
 666   struct interaction_value *i_val2;
 667
 668   var2 = interaction_get_variable (i_var[j]);
 669   i_val2 = interaction_case_data (ccase, i_var[j]);
 670   val2 = interaction_value_get (i_val2);
 671
 672   if (!var_is_value_missing (var2, val2, cov->missing_value))
 673     {
 674       update_hash_entry_intr (cov->ca, var1, var2, val1, val2, i_val1, i_val2);
 675     }
 676 }
 677 /*
 678   Compute the covariance matrix in a single data-pass. Cases with
 679   missing values are dropped pairwise, in other words, only if one of
 680   the two values necessary to accumulate the inner product is missing.
 681
 682   Do not call this function directly. Call it through the struct
 683   covariance_matrix ACCUMULATE member function, for example,
 684   cov->accumulate (cov, ccase).
 685  */
 686 static void
 687 covariance_accumulate_pairwise (struct covariance_matrix *cov,
 688                                 const struct ccase *ccase,
 689                                 const struct interaction_variable **i_var,
 690                                 size_t n_intr)
 691 {
 692   size_t i;
 693   size_t j;
 694   const union value *val1;
 695   const union value *val2;
 696   const struct variable **v_variables;
 697   const struct variable *var1;
 698   const struct variable *var2;
 699   struct interaction_value *i_val1 = NULL;
 700   struct interaction_value *i_val2 = NULL;
 701
 702   assert (cov != NULL);
 703   assert (ccase != NULL);
 704
 705   v_variables = get_covariance_variables (cov);
 706   assert (v_variables != NULL);
 707
 708   for (i = 0; i < cov->n_variables; ++i)
 709     {
 710       var1 = v_variables[i];
 711       val1 = case_data (ccase, var1);
 712       if (!var_is_value_missing (var1, val1, cov->missing_value))
 713         {
 714           cat_value_update (var1, val1);
 715           if (var_is_numeric (var1))
 716             cov->update_moments (cov, i, val1->f);
 717
 718           for (j = i; j < cov->n_variables; j++)
 719             {
 720               var2 = v_variables[j];
 721               val2 = case_data (ccase, var2);
 722               if (!var_is_value_missing
 723                   (var2, val2, cov->missing_value))
 724                 {
 725                   update_hash_entry (cov->ca, var1, var2, val1, val2);
 726                 }
 727             }
 728           for (j = 0; j < cov->n_intr; j++)
 729             {
 730               inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
 731             }
 732         }
 733     }
 734   for (i = 0; i < cov->n_intr; i++)
 735     {
 736       var1 = interaction_get_variable (i_var[i]);
 737       i_val1 = interaction_case_data (ccase, i_var[i]);
 738       val1 = interaction_value_get (i_val1);
 739       cat_value_update (var1, val1);
 740       if (!var_is_value_missing (var1, val1, cov->missing_value))
 741         {
 742           for (j = i; j < cov->n_intr; j++)
 743             {
 744               inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
 745             }
 746         }
 747     }
 748 }
 749
 750 /*
 751   Compute the covariance matrix in a single data-pass. Cases with
 752   missing values are dropped listwise. In other words, if one of the
 753   values for any variable in a case is missing, the entire case is
 754   skipped.
 755
 756   The caller must use a casefilter to remove the cases with missing
 757   values before calling covariance_accumulate_listwise. This function
 758   assumes that CCASE has already passed through this filter, and
 759   contains no missing values.
 760
 761   Do not call this function directly. Call it through the struct
 762   covariance_matrix ACCUMULATE member function, for example,
 763   cov->accumulate (cov, ccase).
 764  */
 765 static void
 766 covariance_accumulate_listwise (struct covariance_matrix *cov,
 767                                 const struct ccase *ccase,
 768                                 const struct interaction_variable **i_var,
 769                                 size_t n_intr)
 770 {
 771   size_t i;
 772   size_t j;
 773   const union value *val1;
 774   const union value *val2;
 775   const struct variable **v_variables;
 776   struct interaction_value *i_val1 = NULL;
 777   struct interaction_value *i_val2 = NULL;
 778
 779   assert (cov != NULL);
 780   assert (ccase != NULL);
 781
 782   v_variables = get_covariance_variables (cov);
 783   assert (v_variables != NULL);
 784
 785   for (i = 0; i < cov->n_variables; ++i)
 786     {
 787       val1 = case_data (ccase, v_variables[i]);
 788       cat_value_update (v_variables[i], val1);
 789       if (var_is_numeric (v_variables[i]))
 790         cov->update_moments (cov, i, val1->f);
 791
 792       for (j = i; j < cov->n_variables; j++)
 793         {
 794           update_hash_entry (cov->ca, v_variables[i], v_variables[j],
 795                              val1, val2);
 796         }
 797     }
 798 }
 799
 800 /*
 801   Call this function during the data pass. Each case will be added to
 802   a hash containing all values of the covariance matrix. After the
 803   data have been passed, call covariance_matrix_compute to put the
 804   values in the struct covariance_matrix.
 805  */
 806 void
 807 covariance_matrix_accumulate (struct covariance_matrix *cov,
 808                               const struct ccase *ccase, void **aux, size_t n_intr)
 809 {
 810   cov->accumulate (cov, ccase, (const struct interaction_variable **) aux, n_intr);
 811 }
 812
 813 /*
 814   Return the value corresponding to subscript TARGET. If that value corresponds
 815   to the origin, return NULL.
 816  */
 817 static const union value *
 818 get_value_from_subscript (const struct design_matrix *dm, size_t target)
 819 {
 820   const union value *result = NULL;
 821   const struct variable *var;
 822   size_t i;
 823
 824   var = design_matrix_col_to_var (dm, target);
 825   if (var_is_numeric (var))
 826     {
 827       return NULL;
 828     }
 829   for (i = 0; i < cat_get_n_categories (var); i++)
 830     {
 831       result = cat_subscript_to_value (i, var);
 832       if (dm_get_exact_subscript (dm, var, result) == target)
 833         {
 834           return result;
 835         }
 836     }
 837   return NULL;
 838 }
 839
 840 static bool
 841 is_covariance_contributor (const struct covariance_accumulator *ca, const struct design_matrix *dm,
 842                            size_t i, size_t j)
 843 {
 844   size_t k;
 845   const struct variable *v1;
 846   const struct variable *v2;
 847
 848   assert (dm != NULL);
 849   v1 = design_matrix_col_to_var (dm, i);
 850   v2 = design_matrix_col_to_var (dm, j);
 851   if (v1 == ca->v1)
 852     {
 853       if (v2 == ca->v2)
 854         {
 855           k = dm_get_exact_subscript (dm, v1, ca->val1);
 856           if (k == i)
 857             {
 858               k = dm_get_exact_subscript (dm, v2, ca->val2);
 859               if (k == j)
 860                 {
 861                   return true;
 862                 }
 863             }
 864         }
 865     }
 866   else if (v1 == ca->v2)
 867     {
 868       if (v2 == ca->v1)
 869         {
 870           k = dm_get_exact_subscript (dm, v1, ca->val2);
 871           if (k == i)
 872             {
 873               k = dm_get_exact_subscript (dm, v2, ca->val1);
 874               if (k == j)
 875                 {
 876                   return true;
 877                 }
 878             }
 879         }
 880     }
 881
 882   return false;
 883 }
 884 static double
 885 get_sum (const struct covariance_matrix *cov, size_t i)
 886 {
 887   size_t k;
 888   double mean;
 889   double n;
 890   const struct variable *var;
 891   const union value *val = NULL;
 892
 893   assert ( cov != NULL);
 894   var = design_matrix_col_to_var (cov->cov, i);
 895   if (var != NULL)
 896     {
 897       if (var_is_alpha (var))
 898         {
 899           val = get_value_from_subscript (cov->cov, i);
 900           k = cat_value_find (var, val);
 901           return cat_get_category_count (k, var);
 902         }
 903       else
 904         {
 905           k = 0;
 906           while (cov->v_variables[k] != var && k  < cov->n_variables)
 907             {
 908               k++;
 909             }
 910           if (k < cov->n_variables)
 911             {
 912               moments1_calculate (cov->m1[k], &n, &mean, NULL, NULL, NULL);
 913               return mean * n;
 914             }
 915         }
 916     }
 917
 918   return 0.0;
 919 }
 920 static void
 921 update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_accumulator *ca)
 922 {
 923   const struct variable *var;
 924   double tmp;
 925   var = design_matrix_col_to_var (dm, i);
 926   if (ca->v1 == var)
 927     {
 928       var = design_matrix_col_to_var (dm, j);
 929       if (ca->v2 == var)
 930         {
 931           tmp = design_matrix_get_element (dm, i, j);
 932           tmp += ca->ssize;
 933           design_matrix_set_element (dm, i, j, tmp);
 934         }
 935     }
 936 }
 937 static void
 938 covariance_accumulator_to_matrix (struct covariance_matrix *cov)
 939 {
 940   size_t i;
 941   size_t j;
 942   double sum_i = 0.0;
 943   double sum_j = 0.0;
 944   double tmp = 0.0;
 945   struct covariance_accumulator *entry;
 946   struct hsh_iterator iter;
 947
 948   cov->cov = covariance_matrix_create_s (cov);
 949   cov->ssize = covariance_matrix_create_s (cov);
 950   entry = hsh_first (cov->ca, &iter);
 951   while (entry != NULL)
 952     {
 953       entry = hsh_next (cov->ca, &iter);
 954     }
 955
 956   for (i = 0; i < design_matrix_get_n_cols (cov->cov); i++)
 957     {
 958       sum_i = get_sum (cov, i);
 959       for (j = i; j < design_matrix_get_n_cols (cov->cov); j++)
 960         {
 961           sum_j = get_sum (cov, j);
 962           entry = hsh_first (cov->ca, &iter);
 963           while (entry != NULL)
 964             {
 965               update_ssize (cov->ssize, i, j, entry);
 966               /*
 967                 We compute the centered, un-normalized covariance matrix.
 968               */
 969               if (is_covariance_contributor (entry, cov->cov, i, j))
 970                 {
 971                   design_matrix_set_element (cov->cov, i, j, entry->dot_product);
 972                 }
 973               entry = hsh_next (cov->ca, &iter);
 974             }
 975           tmp = design_matrix_get_element (cov->cov, i, j);
 976           tmp -= sum_i * sum_j / design_matrix_get_element (cov->ssize, i, j);
 977           design_matrix_set_element (cov->cov, i, j, tmp);
 978           design_matrix_set_element (cov->cov, j, i, tmp);
 979         }
 980     }
 981 }
 982
 983
 984 /*
 985   Call this function after passing the data.
 986  */
 987 void
 988 covariance_matrix_compute (struct covariance_matrix *cov)
 989 {
 990   if (cov->n_pass == ONE_PASS)
 991     {
 992       covariance_accumulator_to_matrix (cov);
 993     }
 994 }
 995
 996 struct design_matrix *
 997 covariance_to_design (const struct covariance_matrix *c)
 998 {
 999   if (c != NULL)
1000     {
1001       return c->cov;
1002     }
1003   return NULL;
1004 }
1005 size_t
1006 covariance_matrix_get_n_rows (const struct covariance_matrix *c)
1007 {
1008   return design_matrix_get_n_rows (c->cov);
1009 }
1010
1011 double
1012 covariance_matrix_get_element (const struct covariance_matrix *c, size_t row, size_t col)
1013 {
1014   return (design_matrix_get_element (c->cov, row, col));
1015 }
1016