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