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 tmp = hsh_hash_bytes (val_max, var_get_width (v_max));
 354       tmp ^= hsh_hash_bytes (val_min, var_get_width (v_min));
 355       tmp += *n_vars * (*n_vars + 1 + idx_max) + idx_min;
 356       return (size_t) tmp;
 357     }
 358   return -1u;
 359 }
 360
 361 /*
 362   Make a hash table consisting of struct covariance_accumulators.
 363   This allows the accumulation of the elements of a covariance matrix
 364   in a single data pass. Call covariance_accumulate () for each case
 365   in the data.
 366  */
 367 static struct hsh_table *
 368 covariance_hsh_create (size_t *n_vars)
 369 {
 370   return hsh_create (*n_vars * *n_vars, covariance_accumulator_compare,
 371                      covariance_accumulator_hash, covariance_accumulator_free,
 372                      n_vars);
 373 }
 374
 375 static void
 376 covariance_accumulator_free (void *c_, const void *aux UNUSED)
 377 {
 378   struct covariance_accumulator *c = c_;
 379   assert (c != NULL);
 380   free (c);
 381 }
 382
 383 /*
 384   Hash comparison. Returns 0 for a match, or a non-zero int
 385   otherwise. The sign of a non-zero return value *should* indicate the
 386   position of C relative to the covariance_accumulator described by
 387   the other arguments. But for now, it just returns 1 for any
 388   non-match.  This should be changed when someone figures out how to
 389   compute a sensible sign for the return value.
 390  */
 391 static int
 392 match_nodes (const struct covariance_accumulator *c,
 393              const struct variable *v1, const struct variable *v2,
 394              const union value *val1, const union value *val2)
 395 {
 396   if (var_get_dict_index (v1) == var_get_dict_index (c->v1))
 397     if (var_get_dict_index (v2) == var_get_dict_index (c->v2))
 398       {
 399         if (var_is_numeric (v1) && var_is_numeric (v2))
 400           {
 401             return 0;
 402           }
 403         if (var_is_numeric (v1) && var_is_alpha (v2))
 404           {
 405             if (!compare_values_short (val2, c->val2, v2))
 406               {
 407                 return 0;
 408               }
 409           }
 410         if (var_is_alpha (v1) && var_is_numeric (v2))
 411           {
 412             if (!compare_values_short (val1, c->val1, v1))
 413               {
 414                 return 0;
 415               }
 416           }
 417         if (var_is_alpha (v1) && var_is_alpha (v2))
 418           {
 419             if (!compare_values_short (val1, c->val1, v1))
 420               {
 421                 if (!compare_values_short (val2, c->val2, v2))
 422                   {
 423                     return 0;
 424                   }
 425               }
 426           }
 427       }
 428   return 1;
 429 }
 430
 431 /*
 432   This function is meant to be used as a comparison function for
 433   a struct hsh_table in src/libpspp/hash.c.
 434 */
 435 static int
 436 covariance_accumulator_compare (const void *a1_, const void *a2_,
 437                                 const void *aux UNUSED)
 438 {
 439   const struct covariance_accumulator *a1 = a1_;
 440   const struct covariance_accumulator *a2 = a2_;
 441
 442   if (a1 == NULL && a2 == NULL)
 443     return 0;
 444
 445   if (a1 == NULL || a2 == NULL)
 446     return 1;
 447
 448   return match_nodes (a1, a2->v1, a2->v2, a2->val1, a2->val2);
 449 }
 450
 451 static unsigned int
 452 hash_numeric_alpha (const struct variable *v1, const struct variable *v2,
 453                     const union value *val, size_t n_vars)
 454 {
 455   unsigned int result = -1u;
 456   if (var_is_numeric (v1) && var_is_alpha (v2))
 457     {
 458       result = n_vars * ((n_vars + 1) + var_get_dict_index (v1))
 459         + var_get_dict_index (v2) + hsh_hash_string (val->s);
 460     }
 461   else if (var_is_alpha (v1) && var_is_numeric (v2))
 462     {
 463       result = hash_numeric_alpha (v2, v1, val, n_vars);
 464     }
 465   return result;
 466 }
 467
 468
 469 static double
 470 update_product (const struct variable *v1, const struct variable *v2,
 471                 const union value *val1, const union value *val2)
 472 {
 473   assert (v1 != NULL);
 474   assert (v2 != NULL);
 475   assert (val1 != NULL);
 476   assert (val2 != NULL);
 477   if (var_is_alpha (v1) && var_is_alpha (v2))
 478     {
 479       return 1.0;
 480     }
 481   if (var_is_numeric (v1) && var_is_numeric (v2))
 482     {
 483       return (val1->f * val2->f);
 484     }
 485   if (var_is_numeric (v1) && var_is_alpha (v2))
 486     {
 487       return (val1->f);
 488     }
 489   if (var_is_numeric (v2) && var_is_alpha (v1))
 490     {
 491       update_product (v2, v1, val2, val1);
 492     }
 493   return 0.0;
 494 }
 495 static double
 496 update_sum (const struct variable *var, const union value *val, double weight)
 497 {
 498   assert (var != NULL);
 499   assert (val != NULL);
 500   if (var_is_alpha (var))
 501     {
 502       return weight;
 503     }
 504   return val->f;
 505 }
 506 static struct covariance_accumulator *
 507 get_new_covariance_accumulator (const struct variable *v1,
 508                                 const struct variable *v2,
 509                                 const union value *val1,
 510                                 const union value *val2)
 511 {
 512   if ((v1 != NULL) && (v2 != NULL) && (val1 != NULL) && (val2 != NULL))
 513     {
 514       struct covariance_accumulator *ca;
 515       ca = xmalloc (sizeof (*ca));
 516       ca->v1 = v1;
 517       ca->v2 = v2;
 518       ca->val1 = val1;
 519       ca->val2 = val2;
 520       return ca;
 521     }
 522   return NULL;
 523 }
 524
 525 static const struct variable **
 526 get_covariance_variables (const struct covariance_matrix *cov)
 527 {
 528   return cov->v_variables;
 529 }
 530
 531 static void
 532 update_hash_entry (struct hsh_table *c,
 533                    const struct variable *v1,
 534                    const struct variable *v2,
 535                    const union value *val1, const union value *val2,
 536                    const struct interaction_value *i_val1,
 537                    const struct interaction_value *i_val2)
 538 {
 539   struct covariance_accumulator *ca;
 540   struct covariance_accumulator *new_entry;
 541   double iv_f1;
 542   double iv_f2;
 543
 544   iv_f1 = interaction_value_get_nonzero_entry (i_val1);
 545   iv_f2 = interaction_value_get_nonzero_entry (i_val2);
 546   ca = get_new_covariance_accumulator (v1, v2, val1, val2);
 547   ca->dot_product = update_product (ca->v1, ca->v2, ca->val1, ca->val2);
 548   ca->dot_product *= iv_f1 * iv_f2;
 549   ca->sum1 = update_sum (ca->v1, ca->val1, iv_f1);
 550   ca->sum2 = update_sum (ca->v2, ca->val2, iv_f2);
 551   ca->ssize = 1.0;
 552   new_entry = hsh_insert (c, ca);
 553
 554   if (new_entry != NULL)
 555     {
 556       new_entry->dot_product += ca->dot_product;
 557       new_entry->ssize += 1.0;
 558       new_entry->sum1 += ca->sum1;
 559       new_entry->sum2 += ca->sum2;
 560       /*
 561         If DOT_PRODUCT is null, CA was not already in the hash
 562         hable, so we don't free it because it was just inserted.
 563         If DOT_PRODUCT was not null, CA is already in the hash table.
 564         Unnecessary now, it must be freed here.
 565       */
 566       free (ca);
 567     }
 568 }
 569
 570 /*
 571   Compute the covariance matrix in a single data-pass. Cases with
 572   missing values are dropped pairwise, in other words, only if one of
 573   the two values necessary to accumulate the inner product is missing.
 574
 575   Do not call this function directly. Call it through the struct
 576   covariance_matrix ACCUMULATE member function, for example,
 577   cov->accumulate (cov, ccase).
 578  */
 579 static void
 580 covariance_accumulate_pairwise (struct covariance_matrix *cov,
 581                                 const struct ccase *ccase,
 582                                 const struct interaction_variable **i_var,
 583                                 size_t n_intr)
 584 {
 585   size_t i;
 586   size_t j;
 587   const union value *val1;
 588   const union value *val2;
 589   const struct variable **v_variables;
 590   struct interaction_value *i_val1 = NULL;
 591   struct interaction_value *i_val2 = NULL;
 592
 593   assert (cov != NULL);
 594   assert (ccase != NULL);
 595
 596   v_variables = get_covariance_variables (cov);
 597   assert (v_variables != NULL);
 598
 599   for (i = 0; i < cov->n_variables; ++i)
 600     {
 601       if (is_interaction (v_variables[i], i_var, n_intr))
 602         {
 603           i_val1 = interaction_case_data (ccase, v_variables[i], i_var, n_intr);
 604           val1 = interaction_value_get (i_val1);
 605         }
 606       else
 607         {
 608           val1 = case_data (ccase, v_variables[i]);
 609         }
 610       if (!var_is_value_missing (v_variables[i], val1, cov->missing_value))
 611         {
 612           cat_value_update (v_variables[i], val1);
 613           if (var_is_numeric (v_variables[i]))
 614             cov->update_moments (cov, i, val1->f);
 615
 616           for (j = i; j < cov->n_variables; j++)
 617             {
 618               if (is_interaction (v_variables[j], i_var, n_intr))
 619                 {
 620                   i_val2 = interaction_case_data (ccase, v_variables[j], i_var, n_intr);
 621                   val2 = interaction_value_get (i_val2);
 622                 }
 623               else
 624                 {
 625                   val2 = case_data (ccase, v_variables[j]);
 626                 }
 627               if (!var_is_value_missing
 628                   (v_variables[j], val2, cov->missing_value))
 629                 {
 630                   update_hash_entry (cov->ca, v_variables[i], v_variables[j],
 631                                      val1, val2, i_val1, i_val2);
 632                   if (j != i)
 633                     update_hash_entry (cov->ca, v_variables[j],
 634                                        v_variables[i], val2, val1, i_val2, i_val1);
 635                 }
 636             }
 637         }
 638     }
 639 }
 640
 641 /*
 642   Compute the covariance matrix in a single data-pass. Cases with
 643   missing values are dropped listwise. In other words, if one of the
 644   values for any variable in a case is missing, the entire case is
 645   skipped.
 646
 647   The caller must use a casefilter to remove the cases with missing
 648   values before calling covariance_accumulate_listwise. This function
 649   assumes that CCASE has already passed through this filter, and
 650   contains no missing values.
 651
 652   Do not call this function directly. Call it through the struct
 653   covariance_matrix ACCUMULATE member function, for example,
 654   cov->accumulate (cov, ccase).
 655  */
 656 static void
 657 covariance_accumulate_listwise (struct covariance_matrix *cov,
 658                                 const struct ccase *ccase,
 659                                 const struct interaction_variable **i_var,
 660                                 size_t n_intr)
 661 {
 662   size_t i;
 663   size_t j;
 664   const union value *val1;
 665   const union value *val2;
 666   const struct variable **v_variables;
 667   struct interaction_value *i_val1 = NULL;
 668   struct interaction_value *i_val2 = NULL;
 669
 670   assert (cov != NULL);
 671   assert (ccase != NULL);
 672
 673   v_variables = get_covariance_variables (cov);
 674   assert (v_variables != NULL);
 675
 676   for (i = 0; i < cov->n_variables; ++i)
 677     {
 678       if (is_interaction (v_variables[i], i_var, n_intr))
 679         {
 680           i_val1 = interaction_case_data (ccase, v_variables[i], i_var, n_intr);
 681           val1 = interaction_value_get (i_val1);
 682         }
 683       else
 684         {
 685           val1 = case_data (ccase, v_variables[i]);
 686         }
 687       cat_value_update (v_variables[i], val1);
 688       if (var_is_numeric (v_variables[i]))
 689         cov->update_moments (cov, i, val1->f);
 690
 691       for (j = i; j < cov->n_variables; j++)
 692         {
 693           if (is_interaction (v_variables[j], i_var, n_intr))
 694             {
 695               i_val2 = interaction_case_data (ccase, v_variables[j], i_var, n_intr);
 696               val2 = interaction_value_get (i_val2);
 697             }
 698           else
 699             {
 700               val2 = case_data (ccase, v_variables[j]);
 701             }
 702           update_hash_entry (cov->ca, v_variables[i], v_variables[j],
 703                              val1, val2, i_val1, i_val2);
 704           if (j != i)
 705             update_hash_entry (cov->ca, v_variables[j], v_variables[i],
 706                                val2, val1, i_val2, i_val1);
 707         }
 708     }
 709 }
 710
 711 /*
 712   Call this function during the data pass. Each case will be added to
 713   a hash containing all values of the covariance matrix. After the
 714   data have been passed, call covariance_matrix_compute to put the
 715   values in the struct covariance_matrix.
 716  */
 717 void
 718 covariance_matrix_accumulate (struct covariance_matrix *cov,
 719                               const struct ccase *ccase, void **aux, size_t n_intr)
 720 {
 721   cov->accumulate (cov, ccase, (const struct interaction_variable **) aux, n_intr);
 722 }
 723 /*
 724   If VAR is categorical with d categories, its first category should
 725   correspond to the origin in d-dimensional Euclidean space.
 726  */
 727 static bool
 728 is_origin (const struct variable *var, const union value *val)
 729 {
 730   if (cat_value_find (var, val) == 0)
 731     {
 732       return true;
 733     }
 734   return false;
 735 }
 736
 737 /*
 738   Return the subscript of the column of the design matrix
 739   corresponding to VAL. If VAR is categorical with d categories, its
 740   first category should correspond to the origin in d-dimensional
 741   Euclidean space, so there is no subscript for this value.
 742  */
 743 static size_t
 744 get_exact_subscript (const struct design_matrix *dm, const struct variable *var,
 745                      const union value *val)
 746 {
 747   size_t result;
 748
 749   if (is_origin (var, val))
 750     {
 751       return -1u;
 752     }
 753
 754   result = design_matrix_var_to_column (dm, var);
 755   if (var_is_alpha (var))
 756     {
 757       result += cat_value_find (var, val) - 1;
 758     }
 759   return result;
 760 }
 761
 762 static void
 763 covariance_matrix_insert (struct design_matrix *cov,
 764                           const struct variable *v1,
 765                           const struct variable *v2, const union value *val1,
 766                           const union value *val2, double product)
 767 {
 768   size_t row;
 769   size_t col;
 770   size_t i;
 771   const union value *tmp_val;
 772
 773   assert (cov != NULL);
 774
 775   row = get_exact_subscript (cov, v1, val1);
 776   col = get_exact_subscript (cov, v2, val2);
 777   if (row != -1u && col != -1u)
 778     {
 779       gsl_matrix_set (cov->m, row, col, product);
 780     }
 781 }
 782
 783
 784 static bool
 785 is_covariance_contributor (const struct covariance_accumulator *ca, const struct design_matrix *dm,
 786                            size_t i, size_t j)
 787 {
 788   size_t k;
 789   const struct variable *v1;
 790   const struct variable *v2;
 791
 792   assert (dm != NULL);
 793   v1 = design_matrix_col_to_var (dm, i);
 794   if (var_get_dict_index (v1) == var_get_dict_index(ca->v1))
 795     {
 796       v2 = design_matrix_col_to_var (dm, j);
 797       if (var_get_dict_index (v2) == var_get_dict_index (ca->v2))
 798         {
 799           k = get_exact_subscript (dm, v1, ca->val1);
 800           if (k == i)
 801             {
 802               k = get_exact_subscript (dm, v2, ca->val2);
 803               if (k == j)
 804                 {
 805                   return true;
 806                 }
 807             }
 808         }
 809     }
 810   return false;
 811 }
 812
 813 static struct design_matrix *
 814 covariance_accumulator_to_matrix (struct covariance_matrix *cov)
 815 {
 816   size_t i;
 817   size_t j;
 818   double tmp;
 819   struct covariance_accumulator *entry;
 820   struct design_matrix *result = NULL;
 821   struct hsh_iterator iter;
 822
 823   result = covariance_matrix_create (cov->n_variables, cov->v_variables);
 824
 825   for (i = 0; i < design_matrix_get_n_cols (result); i++)
 826     {
 827       for (j = i; j < design_matrix_get_n_cols (result); j++)
 828         {
 829           entry = hsh_first (cov->ca, &iter);
 830
 831           while (entry != NULL)
 832             {
 833               /*
 834                 We compute the centered, un-normalized covariance matrix.
 835               */
 836               if (is_covariance_contributor (entry, result, i, j))
 837                 {
 838                   tmp = entry->dot_product - entry->sum1 * entry->sum2 / entry->ssize;
 839                   covariance_matrix_insert (result, entry->v1, entry->v2, entry->val1,
 840                                             entry->val2, tmp);
 841                 }
 842               entry = hsh_next (cov->ca, &iter);
 843             }
 844         }
 845     }
 846   return result;
 847 }
 848
 849
 850 /*
 851   Call this function after passing the data.
 852  */
 853 void
 854 covariance_matrix_compute (struct covariance_matrix *cov)
 855 {
 856   if (cov->n_pass == ONE_PASS)
 857     {
 858       cov->cov = covariance_accumulator_to_matrix (cov);
 859     }
 860 }
 861
 862 struct design_matrix *
 863 covariance_to_design (const struct covariance_matrix *c)
 864 {
 865   if (c != NULL)
 866     {
 867       return c->cov;
 868     }
 869   return NULL;
 870 }