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