src/math/covariance-matrix.c

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