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