Populate the moment matrices with the categorical entries
[pspp] / src / math / covariance.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2009 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 #include <config.h>
18
19 #include <libpspp/assertion.h>
20 #include "covariance.h"
21 #include <gl/xalloc.h>
22 #include "moments.h"
23 #include <gsl/gsl_matrix.h>
24 #include <data/case.h>
25 #include <data/variable.h>
26 #include <libpspp/misc.h>
27 #include "categoricals.h"
28
29 #define n_MOMENTS (MOMENT_VARIANCE + 1)
30
31
32 /* Create a new matrix of NEW_SIZE x NEW_SIZE and copy the elements of
33    matrix IN into it.  IN must be a square matrix, and in normal usage
34    it will be smaller than NEW_SIZE.
35    IN is destroyed by this function.  The return value must be destroyed
36    when no longer required.
37 */
38 static gsl_matrix *
39 resize_matrix (gsl_matrix *in, size_t new_size)
40 {
41   size_t i, j;
42
43   gsl_matrix *out = NULL;
44
45   assert (in->size1 == in->size2);
46
47   if (new_size <= in->size1)
48     return in;
49
50   out = gsl_matrix_calloc (new_size, new_size);
51
52   for (i = 0; i < in->size1; ++i)
53     {
54       for (j = 0; j < in->size2; ++j)
55         {
56           double x = gsl_matrix_get (in, i, j);
57
58           gsl_matrix_set (out, i, j, x);
59         }
60     }
61     
62   gsl_matrix_free (in);
63
64   return out;
65 }
66
67 struct covariance
68 {
69   /* The variables for which the covariance matrix is to be calculated. */
70   size_t n_vars;
71   const struct variable **vars;
72
73   /* Categorical variables. */
74   struct categoricals *categoricals;
75
76   /* Array containing number of categories per categorical variable. */
77   size_t *n_categories;
78
79   /* Dimension of the covariance matrix. */
80   size_t dim;
81
82   /* The weight variable (or NULL if none) */
83   const struct variable *wv;
84
85   /* A set of matrices containing the 0th, 1st and 2nd moments */
86   gsl_matrix **moments;
87
88   /* The class of missing values to exclude */
89   enum mv_class exclude;
90
91   /* An array of doubles representing the covariance matrix.
92      Only the top triangle is included, and no diagonals */
93   double *cm;
94   int n_cm;
95
96   /* 1 for single pass algorithm; 
97      2 for double pass algorithm
98   */
99   short passes;
100
101   /*
102     0 : No pass has  been made
103     1 : First pass has been started
104     2 : Second pass has been 
105     
106     IE: How many passes have been (partially) made. */
107   short state;
108
109   /* Flags indicating that the first case has been seen */
110   bool pass_one_first_case_seen;
111   bool pass_two_first_case_seen;
112 };
113
114
115
116 /* Return a matrix containing the M th moments.
117    The matrix is of size  NxN where N is the number of variables.
118    Each row represents the moments of a variable.
119    In the absence of missing values, the columns of this matrix will
120    be identical.  If missing values are involved, then element (i,j)
121    is the moment of the i th variable, when paired with the j th variable.
122  */
123 const gsl_matrix *
124 covariance_moments (const struct covariance *cov, int m)
125 {
126   return cov->moments[m];
127 }
128
129
130
131 /* Create a covariance struct.
132  */
133 struct covariance *
134 covariance_1pass_create (size_t n_vars, const struct variable **vars,
135                          const struct variable *weight, enum mv_class exclude)
136 {
137   size_t i;
138   struct covariance *cov = xmalloc (sizeof *cov);
139
140   cov->passes = 1;
141   cov->state = 0;
142   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
143   
144   cov->vars = vars;
145
146   cov->wv = weight;
147   cov->n_vars = n_vars;
148   cov->dim = n_vars;
149
150   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
151   
152   for (i = 0; i < n_MOMENTS; ++i)
153     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
154
155   cov->exclude = exclude;
156
157   cov->n_cm = (n_vars * (n_vars - 1)  ) / 2;
158
159   cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
160
161   return cov;
162 }
163
164 /*
165   Create a covariance struct for a two-pass algorithm. If categorical
166   variables are involed, the dimension cannot be know until after the
167   first data pass, so the actual covariances will not be allocated
168   until then.
169  */
170 struct covariance *
171 covariance_2pass_create (size_t n_vars, const struct variable **vars,
172                          size_t n_catvars, const struct variable **catvars, 
173                          const struct variable *wv, enum mv_class exclude)
174 {
175   size_t i;
176   struct covariance *cov = xmalloc (sizeof *cov);
177
178   cov->passes = 2;
179   cov->state = 0;
180   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
181   
182   cov->vars = vars;
183
184   cov->wv = wv;
185   cov->n_vars = n_vars;
186   cov->dim = n_vars;
187
188   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
189   
190   for (i = 0; i < n_MOMENTS; ++i)
191     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
192
193   cov->exclude = exclude;
194
195   cov->n_cm = -1;
196   cov->cm = NULL;
197
198   cov->categoricals = categoricals_create (catvars, n_catvars, wv, exclude);
199
200   return cov;
201 }
202
203 /* Return an integer, which can be used to index 
204    into COV->cm, to obtain the I, J th element
205    of the covariance matrix.  If COV->cm does not
206    contain that element, then a negative value
207    will be returned.
208 */
209 static int
210 cm_idx (const struct covariance *cov, int i, int j)
211 {
212   int as;
213   const int n2j = cov->dim - 2 - j;
214   const int nj = cov->dim - 2 ;
215   
216   assert (i >= 0);
217   assert (j < cov->dim);
218
219   if ( i == 0)
220     return -1;
221
222   if (j >= cov->dim - 1)
223     return -1;
224
225   if ( i <= j) 
226     return -1 ;
227
228   as = nj * (nj + 1) ;
229   as -= n2j * (n2j + 1) ; 
230   as /= 2;
231
232   return i - 1 + as;
233 }
234
235
236 /*
237   Returns true iff the variable corresponding to the Ith element of the covariance matrix 
238    has a missing value for case C
239 */
240 static bool
241 is_missing (const struct covariance *cov, int i, const struct ccase *c)
242 {
243   const struct variable *var = i < cov->n_vars ?
244     cov->vars[i] : 
245     categoricals_get_variable_by_subscript (cov->categoricals, i - cov->n_vars);
246
247   const union value *val = case_data (c, var);
248
249   return var_is_value_missing (var, val, cov->exclude);
250 }
251
252
253 static double
254 get_val (const struct covariance *cov, int i, const struct ccase *c)
255 {
256   if ( i < cov->n_vars)
257     {
258       const struct variable *var = cov->vars[i];
259
260       const union value *val = case_data (c, var);
261
262       return val->f;
263     }
264
265   return categoricals_get_binary_by_subscript (cov->categoricals, i - cov->n_vars, c);
266 }
267
268 void
269 dump_matrix (const gsl_matrix *m)
270 {
271   size_t i, j;
272
273   for (i = 0 ; i < m->size1; ++i)
274     {
275       for (j = 0 ; j < m->size2; ++j)
276         printf ("%02f ", gsl_matrix_get (m, i, j));
277       printf ("\n");
278     }
279 }
280
281 /* Call this function for every case in the data set */
282 void
283 covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
284 {
285   size_t i, j, m;
286   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
287
288   assert (cov->passes == 2);
289   if (!cov->pass_one_first_case_seen)
290     {
291       assert (cov->state == 0);
292       cov->state = 1;
293     }
294
295   categoricals_update (cov->categoricals, c);
296
297   for (i = 0 ; i < cov->dim; ++i)
298     {
299       double v1 = get_val (cov, i, c);
300
301       if ( is_missing (cov, i, c))
302         continue;
303
304       for (j = 0 ; j < cov->dim; ++j)
305         {
306           double pwr = 1.0;
307
308           if ( is_missing (cov, j, c))
309             continue;
310
311           for (m = 0 ; m <= MOMENT_MEAN; ++m)
312             {
313               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
314
315               *x += pwr * weight;
316               pwr *= v1;
317             }
318         }
319     }
320
321   cov->pass_one_first_case_seen = true;
322 }
323
324
325 /* Call this function for every case in the data set */
326 void
327 covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
328 {
329   size_t i, j;
330   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
331
332   assert (cov->passes == 2);
333   assert (cov->state >= 1);
334
335   if (! cov->pass_two_first_case_seen)
336     {
337       size_t m;
338       assert (cov->state == 1);
339       cov->state = 2;
340
341       cov->dim = cov->n_vars + categoricals_total (cov->categoricals);
342       cov->n_cm = (cov->dim * (cov->dim - 1)  ) / 2;
343       cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
344
345       /* Grow the moment matrices so that they're large enough to accommodate the
346          categorical elements */
347       for (i = 0; i < n_MOMENTS; ++i)
348         {
349           cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
350         }
351
352       categoricals_done (cov->categoricals);
353
354       /* Populate the moments matrices with the categorical value elements */
355       for (i = cov->n_vars; i < cov->dim; ++i)
356         {
357           for (j = 0 ; j < cov->dim ; ++j) /* FIXME: This is WRONG !!! */
358             {
359               double w = categoricals_get_weight_by_subscript (cov->categoricals, i - cov->n_vars);
360
361               gsl_matrix_set (cov->moments[MOMENT_NONE], i, j, w);
362
363               w = categoricals_get_sum_by_subscript (cov->categoricals, i - cov->n_vars);
364
365               gsl_matrix_set (cov->moments[MOMENT_MEAN], i, j, w);
366             }
367         }
368
369       /* FIXME: This is WRONG!!  It must be fixed to properly handle missing values.  For
370        now it assumes there are none */
371       for (m = 0 ; m < n_MOMENTS; ++m)
372         {
373           for (i = 0 ; i < cov->dim ; ++i)
374             {
375               double x = gsl_matrix_get (cov->moments[m], i, cov->n_vars -1);
376               for (j = cov->n_vars; j < cov->dim; ++j)
377                 {
378                   gsl_matrix_set (cov->moments[m], i, j, x);
379                 }
380             }
381         }
382
383       /* Divide the means by the number of samples */
384       for (i = 0; i < cov->dim; ++i)
385         {
386           for (j = 0; j < cov->dim; ++j)
387             {
388               double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
389               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
390             }
391         }
392     }
393
394   for (i = 0 ; i < cov->dim; ++i)
395     {
396       double v1 = get_val (cov, i, c);
397
398       if ( is_missing (cov, i, c))
399         continue;
400
401       for (j = 0 ; j < cov->dim; ++j)
402         {
403           int idx;
404           double ss ;
405           double v2 = get_val (cov, j, c);
406
407           const double s = pow2 (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
408
409           if ( is_missing (cov, j, c))
410             continue;
411
412           {
413             double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
414             *x += s;
415           }
416
417           ss = 
418             (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
419             * 
420             (v2 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
421             * weight
422             ;
423
424           idx = cm_idx (cov, i, j);
425           if (idx >= 0)
426             {
427               cov->cm [idx] += ss;
428             }
429         }
430     }
431
432   cov->pass_two_first_case_seen = true;
433 }
434
435
436 /* Call this function for every case in the data set.
437    After all cases have been passed, call covariance_calculate
438  */
439 void
440 covariance_accumulate (struct covariance *cov, const struct ccase *c)
441 {
442   size_t i, j, m;
443   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
444
445   assert (cov->passes == 1);
446
447   if ( !cov->pass_one_first_case_seen)
448     {
449       assert ( cov->state == 0);
450       cov->state = 1;
451     }
452
453   for (i = 0 ; i < cov->dim; ++i)
454     {
455       const union value *val1 = case_data (c, cov->vars[i]);
456
457       if ( is_missing (cov, i, c))
458         continue;
459
460       for (j = 0 ; j < cov->dim; ++j)
461         {
462           double pwr = 1.0;
463           int idx;
464           const union value *val2 = case_data (c, cov->vars[j]);
465
466           if ( is_missing (cov, j, c))
467             continue;
468
469           idx = cm_idx (cov, i, j);
470           if (idx >= 0)
471             {
472               cov->cm [idx] += val1->f * val2->f * weight;
473             }
474
475           for (m = 0 ; m < n_MOMENTS; ++m)
476             {
477               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
478
479               *x += pwr * weight;
480               pwr *= val1->f;
481             }
482         }
483     }
484
485   cov->pass_one_first_case_seen = true;
486 }
487
488
489 /* 
490    Allocate and return a gsl_matrix containing the covariances of the
491    data.
492 */
493 static gsl_matrix *
494 cm_to_gsl (struct covariance *cov)
495 {
496   int i, j;
497   gsl_matrix *m = gsl_matrix_calloc (cov->dim, cov->dim);
498
499   /* Copy the non-diagonal elements from cov->cm */
500   for ( j = 0 ; j < cov->dim - 1; ++j)
501     {
502       for (i = j+1 ; i < cov->dim; ++i)
503         {
504           double x = cov->cm [cm_idx (cov, i, j)];
505           gsl_matrix_set (m, i, j, x);
506           gsl_matrix_set (m, j, i, x);
507         }
508     }
509
510   /* Copy the diagonal elements from cov->moments[2] */
511   for (j = 0 ; j < cov->dim ; ++j)
512     {
513       double sigma = gsl_matrix_get (cov->moments[2], j, j);
514       gsl_matrix_set (m, j, j, sigma);
515     }
516
517   return m;
518 }
519
520
521 static const gsl_matrix *
522 covariance_calculate_double_pass (struct covariance *cov)
523 {
524   size_t i, j;
525   for (i = 0 ; i < cov->dim; ++i)
526     {
527       for (j = 0 ; j < cov->dim; ++j)
528         {
529           int idx;
530           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
531           *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
532
533           idx = cm_idx (cov, i, j);
534           if ( idx >= 0)
535             {
536               x = &cov->cm [idx];
537               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
538             }
539         }
540     }
541
542   return  cm_to_gsl (cov);
543 }
544
545 static const gsl_matrix *
546 covariance_calculate_single_pass (struct covariance *cov)
547 {
548   size_t i, j;
549   size_t m;
550
551   for (m = 0; m < n_MOMENTS; ++m)
552     {
553       /* Divide the moments by the number of samples */
554       if ( m > 0)
555         {
556           for (i = 0 ; i < cov->dim; ++i)
557             {
558               for (j = 0 ; j < cov->dim; ++j)
559                 {
560                   double *x = gsl_matrix_ptr (cov->moments[m], i, j);
561                   *x /= gsl_matrix_get (cov->moments[0], i, j);
562
563                   if ( m == MOMENT_VARIANCE)
564                     *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
565                 }
566             }
567         }
568     }
569
570   /* Centre the moments */
571   for ( j = 0 ; j < cov->dim - 1; ++j)
572     {
573       for (i = j + 1 ; i < cov->dim; ++i)
574         {
575           double *x = &cov->cm [cm_idx (cov, i, j)];
576           
577           *x /= gsl_matrix_get (cov->moments[0], i, j);
578
579           *x -=
580             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j) 
581             *
582             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i); 
583         }
584     }
585
586   return cm_to_gsl (cov);
587 }
588
589
590
591 /* 
592    Return a pointer to gsl_matrix containing the pairwise covariances.
593    The matrix remains owned by the COV object, and must not be freed.
594    Call this function only after all data have been accumulated.
595 */
596 const gsl_matrix *
597 covariance_calculate (struct covariance *cov)
598 {
599   assert ( cov->state > 0 );
600
601   switch (cov->passes)
602     {
603     case 1:
604       return covariance_calculate_single_pass (cov);  
605       break;
606     case 2:
607       return covariance_calculate_double_pass (cov);  
608       break;
609     default:
610       NOT_REACHED ();
611     }
612 }
613
614
615
616
617 /* Destroy the COV object */
618 void
619 covariance_destroy (struct covariance *cov)
620 {
621   size_t i;
622   free (cov->vars);
623   categoricals_destroy (cov->categoricals);
624
625   for (i = 0; i < n_MOMENTS; ++i)
626     gsl_matrix_free (cov->moments[i]);
627
628   free (cov->moments);
629   free (cov->cm);
630   free (cov);
631 }