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