src/math/covariance-matrix.c

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