src/math/covariance-matrix.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2008, 2009 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   size_t idx_max;
 387   size_t idx_min;
 388   const struct variable *v_min;
 389   const struct variable *v_max;
 390   const union value *val_min;
 391   const union value *val_max;
 392
 393   /*
 394      Order everything by the variables' indices. This ensures we get the
 395      same key regardless of the order in which the variables are stored
 396      and passed around.
 397    */
 398   v_min =
 399     (var_get_dict_index (ca->v1) <
 400      var_get_dict_index (ca->v2)) ? ca->v1 : ca->v2;
 401   v_max = (ca->v1 == v_min) ? ca->v2 : ca->v1;
 402
 403   val_min = (v_min == ca->v1) ? ca->val1 : ca->val2;
 404   val_max = (ca->val1 == val_min) ? ca->val2 : ca->val1;
 405
 406   idx_min = var_get_dict_index (v_min);
 407   idx_max = var_get_dict_index (v_max);
 408
 409   if (var_is_numeric (v_max) && var_is_numeric (v_min))
 410     {
 411       return (*n_vars * idx_max + idx_min);
 412     }
 413   if (var_is_numeric (v_max) && var_is_alpha (v_min))
 414     {
 415       return hash_numeric_alpha (v_max, v_min, val_min, *n_vars);
 416     }
 417   if (var_is_alpha (v_max) && var_is_numeric (v_min))
 418     {
 419       return (hash_numeric_alpha (v_min, v_max, val_max, *n_vars));
 420     }
 421   if (var_is_alpha (v_max) && var_is_alpha (v_min))
 422     {
 423       unsigned hash = value_hash (val_max, var_get_width (v_max), 0);
 424       hash = value_hash (val_min, var_get_width (v_min), hash);
 425       return hash_int (*n_vars * (*n_vars + 1 + idx_max) + idx_min, hash);
 426     }
 427   return -1u;
 428 }
 429
 430 /*
 431   Make a hash table consisting of struct covariance_accumulators.
 432   This allows the accumulation of the elements of a covariance matrix
 433   in a single data pass. Call covariance_accumulate () for each case
 434   in the data.
 435  */
 436 static struct hsh_table *
 437 covariance_hsh_create (size_t *n_vars)
 438 {
 439   return hsh_create (*n_vars * *n_vars, covariance_accumulator_compare,
 440                      covariance_accumulator_hash, covariance_accumulator_free,
 441                      n_vars);
 442 }
 443
 444 static void
 445 covariance_accumulator_free (void *c_, const void *aux UNUSED)
 446 {
 447   struct covariance_accumulator *c = c_;
 448   assert (c != NULL);
 449   free (c);
 450 }
 451
 452 static int
 453 ordered_match_nodes (const struct covariance_accumulator *c, const struct variable *v1,
 454                      const struct variable *v2, const union value *val1, const union value *val2)
 455 {
 456   size_t result;
 457   size_t m;
 458
 459   result = var_get_dict_index (v1) ^ var_get_dict_index (c->v1);
 460   m = var_get_dict_index (v2) ^ var_get_dict_index (c->v2);
 461   result = result|m;
 462   if (var_is_alpha (v1))
 463     {
 464       result |= value_compare_3way (val1, c->val1, var_get_width (v1));
 465       if (var_is_alpha (v2))
 466         {
 467           result |= value_compare_3way (val2, c->val2, var_get_width (v2));
 468         }
 469     }
 470   else if (var_is_alpha (v2))
 471     {
 472       result |= value_compare_3way (val2, c->val2, var_get_width (v2));
 473     }
 474   return result;
 475 }
 476
 477 /*
 478   Hash comparison. Returns 0 for a match, or a non-zero int
 479   otherwise. The sign of a non-zero return value *should* indicate the
 480   position of C relative to the covariance_accumulator described by
 481   the other arguments. But for now, it just returns 1 for any
 482   non-match.  This should be changed when someone figures out how to
 483   compute a sensible sign for the return value.
 484  */
 485 static int
 486 match_nodes (const struct covariance_accumulator *c,
 487              const struct variable *v1, const struct variable *v2,
 488              const union value *val1, const union value *val2)
 489 {
 490   size_t n;
 491   size_t m;
 492
 493   n = ordered_match_nodes (c, v1, v2, val1, val2);
 494   m = ordered_match_nodes (c, v2, v1, val2, val1);
 495   return (n & m);
 496 }
 497
 498 /*
 499   This function is meant to be used as a comparison function for
 500   a struct hsh_table in src/libpspp/hash.c.
 501 */
 502 static int
 503 covariance_accumulator_compare (const void *a1_, const void *a2_,
 504                                 const void *aux UNUSED)
 505 {
 506   const struct covariance_accumulator *a1 = a1_;
 507   const struct covariance_accumulator *a2 = a2_;
 508
 509   if (a1 == NULL && a2 == NULL)
 510     return 0;
 511
 512   if (a1 == NULL || a2 == NULL)
 513     return 1;
 514
 515   return match_nodes (a1, a2->v1, a2->v2, a2->val1, a2->val2);
 516 }
 517
 518 static unsigned int
 519 hash_numeric_alpha (const struct variable *v1, const struct variable *v2,
 520                     const union value *val, size_t n_vars)
 521 {
 522   unsigned int result = -1u;
 523   if (var_is_numeric (v1) && var_is_alpha (v2))
 524     {
 525       result = n_vars * ((n_vars + 1) + var_get_dict_index (v1))
 526         + var_get_dict_index (v2) + value_hash (val, var_get_width (v2), 0);
 527     }
 528   else if (var_is_alpha (v1) && var_is_numeric (v2))
 529     {
 530       result = hash_numeric_alpha (v2, v1, val, n_vars);
 531     }
 532   return result;
 533 }
 534
 535
 536 static double
 537 update_product (const struct variable *v1, const struct variable *v2,
 538                 const union value *val1, const union value *val2)
 539 {
 540   assert (v1 != NULL);
 541   assert (v2 != NULL);
 542   assert (val1 != NULL);
 543   assert (val2 != NULL);
 544   if (var_is_alpha (v1) && var_is_alpha (v2))
 545     {
 546       return 1.0;
 547     }
 548   if (var_is_numeric (v1) && var_is_numeric (v2))
 549     {
 550       return (val1->f * val2->f);
 551     }
 552   if (var_is_numeric (v1) && var_is_alpha (v2))
 553     {
 554       return val1->f;
 555     }
 556   if (var_is_numeric (v2) && var_is_alpha (v1))
 557     {
 558       return val2->f;
 559     }
 560   else
 561     {
 562       return 0.0;
 563     }
 564 }
 565 static double
 566 update_sum (const struct variable *var, const union value *val, double weight)
 567 {
 568   assert (var != NULL);
 569   assert (val != NULL);
 570   if (var_is_alpha (var))
 571     {
 572       return weight;
 573     }
 574   return val->f;
 575 }
 576 static struct covariance_accumulator *
 577 get_new_covariance_accumulator (const struct variable *v1,
 578                                 const struct variable *v2,
 579                                 const union value *val1,
 580                                 const union value *val2)
 581 {
 582   if ((v1 != NULL) && (v2 != NULL) && (val1 != NULL) && (val2 != NULL))
 583     {
 584       struct covariance_accumulator *ca;
 585       ca = xmalloc (sizeof (*ca));
 586       ca->v1 = v1;
 587       ca->v2 = v2;
 588       ca->val1 = val1;
 589       ca->val2 = val2;
 590       return ca;
 591     }
 592   return NULL;
 593 }
 594
 595 static const struct variable **
 596 get_covariance_variables (const struct covariance_matrix *cov)
 597 {
 598   return cov->v_variables;
 599 }
 600
 601 static void
 602 update_hash_entry_intr (struct hsh_table *c,
 603                         const struct variable *v1,
 604                         const struct variable *v2,
 605                         const union value *val1, const union value *val2,
 606                         const struct interaction_value *i_val1,
 607                         const struct interaction_value *i_val2)
 608 {
 609   struct covariance_accumulator *ca;
 610   struct covariance_accumulator *new_entry;
 611   double iv_f1;
 612   double iv_f2;
 613
 614   iv_f1 = interaction_value_get_nonzero_entry (i_val1);
 615   iv_f2 = interaction_value_get_nonzero_entry (i_val2);
 616   ca = get_new_covariance_accumulator (v1, v2, val1, val2);
 617   ca->dot_product = update_product (ca->v1, ca->v2, ca->val1, ca->val2);
 618   ca->dot_product *= iv_f1 * iv_f2;
 619   ca->sum1 = update_sum (ca->v1, ca->val1, iv_f1);
 620   ca->sum2 = update_sum (ca->v2, ca->val2, iv_f2);
 621   ca->ssize = 1.0;
 622   new_entry = hsh_insert (c, ca);
 623
 624   if (new_entry != NULL)
 625     {
 626       new_entry->dot_product += ca->dot_product;
 627       new_entry->ssize += 1.0;
 628       new_entry->sum1 += ca->sum1;
 629       new_entry->sum2 += ca->sum2;
 630       /*
 631         If DOT_PRODUCT is null, CA was not already in the hash
 632         hable, so we don't free it because it was just inserted.
 633         If DOT_PRODUCT was not null, CA is already in the hash table.
 634         Unnecessary now, it must be freed here.
 635       */
 636       free (ca);
 637     }
 638 }
 639
 640 static void
 641 update_hash_entry (struct hsh_table *c,
 642                    const struct variable *v1,
 643                    const struct variable *v2,
 644                    const union value *val1, const union value *val2)
 645 {
 646   struct covariance_accumulator *ca;
 647   struct covariance_accumulator *new_entry;
 648
 649   ca = get_new_covariance_accumulator (v1, v2, val1, val2);
 650   ca->dot_product = update_product (ca->v1, ca->v2, ca->val1, ca->val2);
 651   ca->sum1 = update_sum (ca->v1, ca->val1, 1.0);
 652   ca->sum2 = update_sum (ca->v2, ca->val2, 1.0);
 653   ca->ssize = 1.0;
 654   new_entry = hsh_insert (c, ca);
 655
 656   if (new_entry != NULL)
 657     {
 658       new_entry->dot_product += ca->dot_product;
 659       new_entry->ssize += 1.0;
 660       new_entry->sum1 += ca->sum1;
 661       new_entry->sum2 += ca->sum2;
 662       /*
 663         If DOT_PRODUCT is null, CA was not already in the hash
 664         hable, so we don't free it because it was just inserted.
 665         If DOT_PRODUCT was not null, CA is already in the hash table.
 666         Unnecessary now, it must be freed here.
 667       */
 668       free (ca);
 669     }
 670 }
 671
 672 static void
 673 inner_intr_loop (struct covariance_matrix *cov, const struct ccase  *ccase, const struct variable *var1,
 674                  const union value *val1, const struct interaction_variable **i_var,
 675                  const struct interaction_value *i_val1, size_t j)
 676 {
 677   struct variable *var2;
 678   union value *val2;
 679   struct interaction_value *i_val2;
 680
 681   var2 = interaction_get_variable (i_var[j]);
 682   i_val2 = interaction_case_data (ccase, i_var[j]);
 683   val2 = interaction_value_get (i_val2);
 684
 685   if (!var_is_value_missing (var2, val2, cov->missing_value))
 686     {
 687       update_hash_entry_intr (cov->ca, var1, var2, val1, val2, i_val1, i_val2);
 688     }
 689 }
 690 /*
 691   Compute the covariance matrix in a single data-pass. Cases with
 692   missing values are dropped pairwise, in other words, only if one of
 693   the two values necessary to accumulate the inner product is missing.
 694
 695   Do not call this function directly. Call it through the struct
 696   covariance_matrix ACCUMULATE member function, for example,
 697   cov->accumulate (cov, ccase).
 698  */
 699 static void
 700 covariance_accumulate_pairwise (struct covariance_matrix *cov,
 701                                 const struct ccase *ccase,
 702                                 const struct interaction_variable **i_var,
 703                                 size_t n_intr)
 704 {
 705   size_t i;
 706   size_t j;
 707   const union value *val1;
 708   const union value *val2;
 709   const struct variable **v_variables;
 710   const struct variable *var1;
 711   const struct variable *var2;
 712   struct interaction_value *i_val1 = NULL;
 713   struct interaction_value *i_val2 = NULL;
 714
 715   assert (cov != NULL);
 716   assert (ccase != NULL);
 717
 718   v_variables = get_covariance_variables (cov);
 719   assert (v_variables != NULL);
 720
 721   for (i = 0; i < cov->n_variables; ++i)
 722     {
 723       var1 = v_variables[i];
 724       val1 = case_data (ccase, var1);
 725       if (!var_is_value_missing (var1, val1, cov->missing_value))
 726         {
 727           cat_value_update (var1, val1);
 728           if (var_is_numeric (var1))
 729             cov->update_moments (cov, i, val1->f);
 730
 731           for (j = i; j < cov->n_variables; j++)
 732             {
 733               var2 = v_variables[j];
 734               val2 = case_data (ccase, var2);
 735               if (!var_is_value_missing
 736                   (var2, val2, cov->missing_value))
 737                 {
 738                   update_hash_entry (cov->ca, var1, var2, val1, val2);
 739                 }
 740             }
 741           for (j = 0; j < cov->n_intr; j++)
 742             {
 743               inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
 744             }
 745         }
 746     }
 747   for (i = 0; i < cov->n_intr; i++)
 748     {
 749       var1 = interaction_get_variable (i_var[i]);
 750       i_val1 = interaction_case_data (ccase, i_var[i]);
 751       val1 = interaction_value_get (i_val1);
 752       cat_value_update (var1, val1);
 753       if (!var_is_value_missing (var1, val1, cov->missing_value))
 754         {
 755           for (j = i; j < cov->n_intr; j++)
 756             {
 757               inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
 758             }
 759         }
 760     }
 761 }
 762
 763 /*
 764   Compute the covariance matrix in a single data-pass. Cases with
 765   missing values are dropped listwise. In other words, if one of the
 766   values for any variable in a case is missing, the entire case is
 767   skipped.
 768
 769   The caller must use a casefilter to remove the cases with missing
 770   values before calling covariance_accumulate_listwise. This function
 771   assumes that CCASE has already passed through this filter, and
 772   contains no missing values.
 773
 774   Do not call this function directly. Call it through the struct
 775   covariance_matrix ACCUMULATE member function, for example,
 776   cov->accumulate (cov, ccase).
 777  */
 778 static void
 779 covariance_accumulate_listwise (struct covariance_matrix *cov,
 780                                 const struct ccase *ccase,
 781                                 const struct interaction_variable **i_var,
 782                                 size_t n_intr)
 783 {
 784   size_t i;
 785   size_t j;
 786   const union value *val1;
 787   const union value *val2;
 788   const struct variable **v_variables;
 789   struct interaction_value *i_val1 = NULL;
 790   struct interaction_value *i_val2 = NULL;
 791
 792   assert (cov != NULL);
 793   assert (ccase != NULL);
 794
 795   v_variables = get_covariance_variables (cov);
 796   assert (v_variables != NULL);
 797
 798   for (i = 0; i < cov->n_variables; ++i)
 799     {
 800       val1 = case_data (ccase, v_variables[i]);
 801       cat_value_update (v_variables[i], val1);
 802       if (var_is_numeric (v_variables[i]))
 803         cov->update_moments (cov, i, val1->f);
 804
 805       for (j = i; j < cov->n_variables; j++)
 806         {
 807           update_hash_entry (cov->ca, v_variables[i], v_variables[j],
 808                              val1, val2);
 809         }
 810     }
 811 }
 812
 813 /*
 814   Call this function during the data pass. Each case will be added to
 815   a hash containing all values of the covariance matrix. After the
 816   data have been passed, call covariance_matrix_compute to put the
 817   values in the struct covariance_matrix.
 818  */
 819 void
 820 covariance_matrix_accumulate (struct covariance_matrix *cov,
 821                               const struct ccase *ccase, void **aux, size_t n_intr)
 822 {
 823   cov->accumulate (cov, ccase, (const struct interaction_variable **) aux, n_intr);
 824 }
 825
 826 /*
 827   Return the value corresponding to subscript TARGET. If that value corresponds
 828   to the origin, return NULL.
 829  */
 830 static const union value *
 831 get_value_from_subscript (const struct design_matrix *dm, size_t target)
 832 {
 833   const union value *result = NULL;
 834   const struct variable *var;
 835   size_t i;
 836
 837   var = design_matrix_col_to_var (dm, target);
 838   if (var_is_numeric (var))
 839     {
 840       return NULL;
 841     }
 842   for (i = 0; i < cat_get_n_categories (var); i++)
 843     {
 844       result = cat_subscript_to_value (i, var);
 845       if (dm_get_exact_subscript (dm, var, result) == target)
 846         {
 847           return result;
 848         }
 849     }
 850   return NULL;
 851 }
 852
 853 static bool
 854 is_covariance_contributor (const struct covariance_accumulator *ca, const struct design_matrix *dm,
 855                            size_t i, size_t j)
 856 {
 857   size_t k;
 858   const struct variable *v1;
 859   const struct variable *v2;
 860
 861   assert (dm != NULL);
 862   v1 = design_matrix_col_to_var (dm, i);
 863   v2 = design_matrix_col_to_var (dm, j);
 864   if (var_get_dict_index (v1) == var_get_dict_index(ca->v1))
 865     {
 866       if (var_get_dict_index (v2) == var_get_dict_index (ca->v2))
 867         {
 868           k = dm_get_exact_subscript (dm, v1, ca->val1);
 869           if (k == i)
 870             {
 871               k = dm_get_exact_subscript (dm, v2, ca->val2);
 872               if (k == j)
 873                 {
 874                   return true;
 875                 }
 876             }
 877         }
 878     }
 879   else if (var_get_dict_index (v1) == var_get_dict_index (ca->v2))
 880     {
 881       if (var_get_dict_index (v2) == var_get_dict_index (ca->v1))
 882         {
 883           k = dm_get_exact_subscript (dm, v1, ca->val2);
 884           if (k == i)
 885             {
 886               k = dm_get_exact_subscript (dm, v2, ca->val1);
 887               if (k == j)
 888                 {
 889                   return true;
 890                 }
 891             }
 892         }
 893     }
 894
 895   return false;
 896 }
 897 static double
 898 get_sum (const struct covariance_matrix *cov, size_t i)
 899 {
 900   size_t k;
 901   double mean;
 902   double n;
 903   const struct variable *var;
 904   const union value *val = NULL;
 905
 906   assert ( cov != NULL);
 907   var = design_matrix_col_to_var (cov->cov, i);
 908   if (var != NULL)
 909     {
 910       if (var_is_alpha (var))
 911         {
 912           val = get_value_from_subscript (cov->cov, i);
 913           k = cat_value_find (var, val);
 914           return cat_get_category_count (k, var);
 915         }
 916       else
 917         {
 918           k = 0;
 919           while (cov->v_variables[k] != var && k  < cov->n_variables)
 920             {
 921               k++;
 922             }
 923           if (k < cov->n_variables)
 924             {
 925               moments1_calculate (cov->m1[k], &n, &mean, NULL, NULL, NULL);
 926               return mean * n;
 927             }
 928         }
 929     }
 930
 931   return 0.0;
 932 }
 933 static void
 934 update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_accumulator *ca)
 935 {
 936   const struct variable *var;
 937   double tmp;
 938   var = design_matrix_col_to_var (dm, i);
 939   if (var_get_dict_index (ca->v1) == var_get_dict_index (var))
 940     {
 941       var = design_matrix_col_to_var (dm, j);
 942       if (var_get_dict_index (ca->v2) == var_get_dict_index (var))
 943         {
 944           tmp = design_matrix_get_element (dm, i, j);
 945           tmp += ca->ssize;
 946           design_matrix_set_element (dm, i, j, tmp);
 947         }
 948     }
 949 }
 950 static void
 951 covariance_accumulator_to_matrix (struct covariance_matrix *cov)
 952 {
 953   size_t i;
 954   size_t j;
 955   double sum_i = 0.0;
 956   double sum_j = 0.0;
 957   double tmp = 0.0;
 958   struct covariance_accumulator *entry;
 959   struct hsh_iterator iter;
 960
 961   cov->cov = covariance_matrix_create_s (cov);
 962   cov->ssize = covariance_matrix_create_s (cov);
 963   entry = hsh_first (cov->ca, &iter);
 964   while (entry != NULL)
 965     {
 966       entry = hsh_next (cov->ca, &iter);
 967     }
 968
 969   for (i = 0; i < design_matrix_get_n_cols (cov->cov); i++)
 970     {
 971       sum_i = get_sum (cov, i);
 972       for (j = i; j < design_matrix_get_n_cols (cov->cov); j++)
 973         {
 974           sum_j = get_sum (cov, j);
 975           entry = hsh_first (cov->ca, &iter);
 976           while (entry != NULL)
 977             {
 978               update_ssize (cov->ssize, i, j, entry);
 979               /*
 980                 We compute the centered, un-normalized covariance matrix.
 981               */
 982               if (is_covariance_contributor (entry, cov->cov, i, j))
 983                 {
 984                   design_matrix_set_element (cov->cov, i, j, entry->dot_product);
 985                 }
 986               entry = hsh_next (cov->ca, &iter);
 987             }
 988           tmp = design_matrix_get_element (cov->cov, i, j);
 989           tmp -= sum_i * sum_j / design_matrix_get_element (cov->ssize, i, j);
 990           design_matrix_set_element (cov->cov, i, j, tmp);
 991           design_matrix_set_element (cov->cov, j, i, tmp);
 992         }
 993     }
 994 }
 995
 996
 997 /*
 998   Call this function after passing the data.
 999  */
1000 void
1001 covariance_matrix_compute (struct covariance_matrix *cov)
1002 {
1003   if (cov->n_pass == ONE_PASS)
1004     {
1005       covariance_accumulator_to_matrix (cov);
1006     }
1007 }
1008
1009 struct design_matrix *
1010 covariance_to_design (const struct covariance_matrix *c)
1011 {
1012   if (c != NULL)
1013     {
1014       return c->cov;
1015     }
1016   return NULL;
1017 }
1018 size_t
1019 covariance_matrix_get_n_rows (const struct covariance_matrix *c)
1020 {
1021   return design_matrix_get_n_rows (c->cov);
1022 }
1023
1024 double
1025 covariance_matrix_get_element (const struct covariance_matrix *c, size_t row, size_t col)
1026 {
1027   return (design_matrix_get_element (c->cov, row, col));
1028 }
1029