categoricals.c (struct payload): Rename destroy to calculate
[pspp] / src / math / covariance.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2009, 2010, 2011 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 "math/covariance.h"
20
21 #include <gsl/gsl_matrix.h>
22
23 #include "data/case.h"
24 #include "data/variable.h"
25 #include "libpspp/assertion.h"
26 #include "libpspp/misc.h"
27 #include "math/categoricals.h"
28 #include "math/interaction.h"
29 #include "math/moments.h"
30
31 #include "gl/xalloc.h"
32
33 #define n_MOMENTS (MOMENT_VARIANCE + 1)
34
35
36 /* Create a new matrix of NEW_SIZE x NEW_SIZE and copy the elements of
37    matrix IN into it.  IN must be a square matrix, and in normal usage
38    it will be smaller than NEW_SIZE.
39    IN is destroyed by this function.  The return value must be destroyed
40    when no longer required.
41 */
42 static gsl_matrix *
43 resize_matrix (gsl_matrix *in, size_t new_size)
44 {
45   size_t i, j;
46
47   gsl_matrix *out = NULL;
48
49   assert (in->size1 == in->size2);
50
51   if (new_size <= in->size1)
52     return in;
53
54   out = gsl_matrix_calloc (new_size, new_size);
55
56   for (i = 0; i < in->size1; ++i)
57     {
58       for (j = 0; j < in->size2; ++j)
59         {
60           double x = gsl_matrix_get (in, i, j);
61
62           gsl_matrix_set (out, i, j, x);
63         }
64     }
65     
66   gsl_matrix_free (in);
67
68   return out;
69 }
70
71 struct covariance
72 {
73   /* The variables for which the covariance matrix is to be calculated. */
74   size_t n_vars;
75   const struct variable *const *vars;
76
77   /* Categorical variables. */
78   struct categoricals *categoricals;
79
80   /* Array containing number of categories per categorical variable. */
81   size_t *n_categories;
82
83   /* Dimension of the covariance matrix. */
84   size_t dim;
85
86   /* The weight variable (or NULL if none) */
87   const struct variable *wv;
88
89   /* A set of matrices containing the 0th, 1st and 2nd moments */
90   gsl_matrix **moments;
91
92   /* The class of missing values to exclude */
93   enum mv_class exclude;
94
95   /* An array of doubles representing the covariance matrix.
96      Only the top triangle is included, and no diagonals */
97   double *cm;
98   int n_cm;
99
100   /* 1 for single pass algorithm; 
101      2 for double pass algorithm
102   */
103   short passes;
104
105   /*
106     0 : No pass has  been made
107     1 : First pass has been started
108     2 : Second pass has been 
109     
110     IE: How many passes have been (partially) made. */
111   short state;
112
113   /* Flags indicating that the first case has been seen */
114   bool pass_one_first_case_seen;
115   bool pass_two_first_case_seen;
116 };
117
118
119
120 /* Return a matrix containing the M th moments.
121    The matrix is of size  NxN where N is the number of variables.
122    Each row represents the moments of a variable.
123    In the absence of missing values, the columns of this matrix will
124    be identical.  If missing values are involved, then element (i,j)
125    is the moment of the i th variable, when paired with the j th variable.
126  */
127 const gsl_matrix *
128 covariance_moments (const struct covariance *cov, int m)
129 {
130   return cov->moments[m];
131 }
132
133
134
135 /* Create a covariance struct.
136  */
137 struct covariance *
138 covariance_1pass_create (size_t n_vars, const struct variable *const *vars,
139                          const struct variable *weight, enum mv_class exclude)
140 {
141   size_t i;
142   struct covariance *cov = xzalloc (sizeof *cov);
143
144   cov->passes = 1;
145   cov->state = 0;
146   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
147   
148   cov->vars = vars;
149
150   cov->wv = weight;
151   cov->n_vars = n_vars;
152   cov->dim = n_vars;
153
154   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
155   
156   for (i = 0; i < n_MOMENTS; ++i)
157     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
158
159   cov->exclude = exclude;
160
161   cov->n_cm = (n_vars * (n_vars - 1)  ) / 2;
162
163
164   cov->cm = xcalloc (cov->n_cm, sizeof *cov->cm);
165   cov->categoricals = NULL;
166
167   return cov;
168 }
169
170 /*
171   Create a covariance struct for a two-pass algorithm. If categorical
172   variables are involed, the dimension cannot be know until after the
173   first data pass, so the actual covariances will not be allocated
174   until then.
175  */
176 struct covariance *
177 covariance_2pass_create (size_t n_vars, const struct variable *const *vars,
178                          struct categoricals *cats,
179                          const struct variable *wv, enum mv_class exclude)
180 {
181   size_t i;
182   struct covariance *cov = xmalloc (sizeof *cov);
183
184   cov->passes = 2;
185   cov->state = 0;
186   cov->pass_one_first_case_seen = cov->pass_two_first_case_seen = false;
187   
188   cov->vars = vars;
189
190   cov->wv = wv;
191   cov->n_vars = n_vars;
192   cov->dim = n_vars;
193
194   cov->moments = xmalloc (sizeof *cov->moments * n_MOMENTS);
195   
196   for (i = 0; i < n_MOMENTS; ++i)
197     cov->moments[i] = gsl_matrix_calloc (n_vars, n_vars);
198
199   cov->exclude = exclude;
200
201   cov->n_cm = -1;
202   cov->cm = NULL;
203
204   cov->categoricals = cats;
205
206   return cov;
207 }
208
209 /* Return an integer, which can be used to index 
210    into COV->cm, to obtain the I, J th element
211    of the covariance matrix.  If COV->cm does not
212    contain that element, then a negative value
213    will be returned.
214 */
215 static int
216 cm_idx (const struct covariance *cov, int i, int j)
217 {
218   int as;
219   const int n2j = cov->dim - 2 - j;
220   const int nj = cov->dim - 2 ;
221   
222   assert (i >= 0);
223   assert (j < cov->dim);
224
225   if ( i == 0)
226     return -1;
227
228   if (j >= cov->dim - 1)
229     return -1;
230
231   if ( i <= j) 
232     return -1 ;
233
234   as = nj * (nj + 1) ;
235   as -= n2j * (n2j + 1) ; 
236   as /= 2;
237
238   return i - 1 + as;
239 }
240
241
242 /*
243   Returns true iff the variable corresponding to the Ith element of the covariance matrix 
244    has a missing value for case C
245 */
246 static bool
247 is_missing (const struct covariance *cov, int i, const struct ccase *c)
248 {
249   const struct variable *var = i < cov->n_vars ?
250     cov->vars[i] : 
251     categoricals_get_interaction_by_subscript (cov->categoricals, i - cov->n_vars)->vars[0];
252
253   const union value *val = case_data (c, var);
254
255   return var_is_value_missing (var, val, cov->exclude);
256 }
257
258
259 static double
260 get_val (const struct covariance *cov, int i, const struct ccase *c)
261 {
262   if ( i < cov->n_vars)
263     {
264       const struct variable *var = cov->vars[i];
265
266       const union value *val = case_data (c, var);
267
268       return val->f;
269     }
270
271   return categoricals_get_code_for_case (cov->categoricals, i - cov->n_vars, c);
272 }
273
274 #if 0
275 void
276 dump_matrix (const gsl_matrix *m)
277 {
278   size_t i, j;
279
280   for (i = 0 ; i < m->size1; ++i)
281     {
282       for (j = 0 ; j < m->size2; ++j)
283         printf ("%02f ", gsl_matrix_get (m, i, j));
284       printf ("\n");
285     }
286 }
287 #endif
288
289 /* Call this function for every case in the data set */
290 void
291 covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
292 {
293   size_t i, j, m;
294   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
295
296   assert (cov->passes == 2);
297   if (!cov->pass_one_first_case_seen)
298     {
299       assert (cov->state == 0);
300       cov->state = 1;
301     }
302
303   if (cov->categoricals)
304     categoricals_update (cov->categoricals, c);
305
306   for (i = 0 ; i < cov->dim; ++i)
307     {
308       double v1 = get_val (cov, i, c);
309
310       if ( is_missing (cov, i, c))
311         continue;
312
313       for (j = 0 ; j < cov->dim; ++j)
314         {
315           double pwr = 1.0;
316
317           if ( is_missing (cov, j, c))
318             continue;
319
320           for (m = 0 ; m <= MOMENT_MEAN; ++m)
321             {
322               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
323
324               *x += pwr * weight;
325               pwr *= v1;
326             }
327         }
328     }
329
330   cov->pass_one_first_case_seen = true;
331 }
332
333
334 /* Call this function for every case in the data set */
335 void
336 covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
337 {
338   size_t i, j;
339   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
340
341   assert (cov->passes == 2);
342   assert (cov->state >= 1);
343
344   if (! cov->pass_two_first_case_seen)
345     {
346       size_t m;
347       assert (cov->state == 1);
348       cov->state = 2;
349
350       if (cov->categoricals)
351         categoricals_done (cov->categoricals);
352
353       cov->dim = cov->n_vars;
354       
355       if (cov->categoricals)
356         cov->dim += categoricals_df_total (cov->categoricals);
357
358       cov->n_cm = (cov->dim * (cov->dim - 1)  ) / 2;
359       cov->cm = xcalloc (cov->n_cm, sizeof *cov->cm);
360
361       /* Grow the moment matrices so that they're large enough to accommodate the
362          categorical elements */
363       for (i = 0; i < n_MOMENTS; ++i)
364         {
365           cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
366         }
367
368       /* Populate the moments matrices with the categorical value elements */
369       for (i = cov->n_vars; i < cov->dim; ++i)
370         {
371           for (j = 0 ; j < cov->dim ; ++j) /* FIXME: This is WRONG !!! */
372             {
373               double w = categoricals_get_weight_by_subscript (cov->categoricals, i - cov->n_vars);
374
375               gsl_matrix_set (cov->moments[MOMENT_NONE], i, j, w);
376
377               w = categoricals_get_sum_by_subscript (cov->categoricals, i - cov->n_vars);
378
379               gsl_matrix_set (cov->moments[MOMENT_MEAN], i, j, w);
380             }
381         }
382
383       /* FIXME: This is WRONG!!  It must be fixed to properly handle missing values.  For
384        now it assumes there are none */
385       for (m = 0 ; m < n_MOMENTS; ++m)
386         {
387           for (i = 0 ; i < cov->dim ; ++i)
388             {
389               double x = gsl_matrix_get (cov->moments[m], i, cov->n_vars -1);
390               for (j = cov->n_vars; j < cov->dim; ++j)
391                 {
392                   gsl_matrix_set (cov->moments[m], i, j, x);
393                 }
394             }
395         }
396
397       /* Divide the means by the number of samples */
398       for (i = 0; i < cov->dim; ++i)
399         {
400           for (j = 0; j < cov->dim; ++j)
401             {
402               double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
403               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
404             }
405         }
406     }
407
408   for (i = 0 ; i < cov->dim; ++i)
409     {
410       double v1 = get_val (cov, i, c);
411
412       if ( is_missing (cov, i, c))
413         continue;
414
415       for (j = 0 ; j < cov->dim; ++j)
416         {
417           int idx;
418           double ss ;
419           double v2 = get_val (cov, j, c);
420
421           const double s = pow2 (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
422
423           if ( is_missing (cov, j, c))
424             continue;
425
426           {
427             double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
428             *x += s;
429           }
430
431           ss = 
432             (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
433             * 
434             (v2 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
435             * weight
436             ;
437
438           idx = cm_idx (cov, i, j);
439           if (idx >= 0)
440             {
441               cov->cm [idx] += ss;
442             }
443         }
444     }
445
446   cov->pass_two_first_case_seen = true;
447 }
448
449
450 /* Call this function for every case in the data set.
451    After all cases have been passed, call covariance_calculate
452  */
453 void
454 covariance_accumulate (struct covariance *cov, const struct ccase *c)
455 {
456   size_t i, j, m;
457   const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
458
459   assert (cov->passes == 1);
460
461   if ( !cov->pass_one_first_case_seen)
462     {
463       assert ( cov->state == 0);
464       cov->state = 1;
465     }
466
467   for (i = 0 ; i < cov->dim; ++i)
468     {
469       const union value *val1 = case_data (c, cov->vars[i]);
470
471       if ( is_missing (cov, i, c))
472         continue;
473
474       for (j = 0 ; j < cov->dim; ++j)
475         {
476           double pwr = 1.0;
477           int idx;
478           const union value *val2 = case_data (c, cov->vars[j]);
479
480           if ( is_missing (cov, j, c))
481             continue;
482
483           idx = cm_idx (cov, i, j);
484           if (idx >= 0)
485             {
486               cov->cm [idx] += val1->f * val2->f * weight;
487             }
488
489           for (m = 0 ; m < n_MOMENTS; ++m)
490             {
491               double *x = gsl_matrix_ptr (cov->moments[m], i, j);
492
493               *x += pwr * weight;
494               pwr *= val1->f;
495             }
496         }
497     }
498
499   cov->pass_one_first_case_seen = true;
500 }
501
502
503 /* 
504    Allocate and return a gsl_matrix containing the covariances of the
505    data.
506 */
507 static gsl_matrix *
508 cm_to_gsl (struct covariance *cov)
509 {
510   int i, j;
511   gsl_matrix *m = gsl_matrix_calloc (cov->dim, cov->dim);
512
513   /* Copy the non-diagonal elements from cov->cm */
514   for ( j = 0 ; j < cov->dim - 1; ++j)
515     {
516       for (i = j+1 ; i < cov->dim; ++i)
517         {
518           double x = cov->cm [cm_idx (cov, i, j)];
519           gsl_matrix_set (m, i, j, x);
520           gsl_matrix_set (m, j, i, x);
521         }
522     }
523
524   /* Copy the diagonal elements from cov->moments[2] */
525   for (j = 0 ; j < cov->dim ; ++j)
526     {
527       double sigma = gsl_matrix_get (cov->moments[2], j, j);
528       gsl_matrix_set (m, j, j, sigma);
529     }
530
531   return m;
532 }
533
534
535 static gsl_matrix *
536 covariance_calculate_double_pass (struct covariance *cov)
537 {
538   size_t i, j;
539   for (i = 0 ; i < cov->dim; ++i)
540     {
541       for (j = 0 ; j < cov->dim; ++j)
542         {
543           int idx;
544           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
545           *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
546
547           idx = cm_idx (cov, i, j);
548           if ( idx >= 0)
549             {
550               x = &cov->cm [idx];
551               *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
552             }
553         }
554     }
555
556   return  cm_to_gsl (cov);
557 }
558
559 static gsl_matrix *
560 covariance_calculate_single_pass (struct covariance *cov)
561 {
562   size_t i, j;
563   size_t m;
564
565   for (m = 0; m < n_MOMENTS; ++m)
566     {
567       /* Divide the moments by the number of samples */
568       if ( m > 0)
569         {
570           for (i = 0 ; i < cov->dim; ++i)
571             {
572               for (j = 0 ; j < cov->dim; ++j)
573                 {
574                   double *x = gsl_matrix_ptr (cov->moments[m], i, j);
575                   *x /= gsl_matrix_get (cov->moments[0], i, j);
576
577                   if ( m == MOMENT_VARIANCE)
578                     *x -= pow2 (gsl_matrix_get (cov->moments[1], i, j));
579                 }
580             }
581         }
582     }
583
584   /* Centre the moments */
585   for ( j = 0 ; j < cov->dim - 1; ++j)
586     {
587       for (i = j + 1 ; i < cov->dim; ++i)
588         {
589           double *x = &cov->cm [cm_idx (cov, i, j)];
590           
591           *x /= gsl_matrix_get (cov->moments[0], i, j);
592
593           *x -=
594             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j) 
595             *
596             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i); 
597         }
598     }
599
600   return cm_to_gsl (cov);
601 }
602
603
604 /* Return a pointer to gsl_matrix containing the pairwise covariances.  The
605    caller owns the returned matrix and must free it when it is no longer
606    needed.
607
608    Call this function only after all data have been accumulated.  */
609 gsl_matrix *
610 covariance_calculate (struct covariance *cov)
611 {
612   if ( cov->state <= 0 )
613     return NULL;
614
615   switch (cov->passes)
616     {
617     case 1:
618       return covariance_calculate_single_pass (cov);  
619       break;
620     case 2:
621       return covariance_calculate_double_pass (cov);  
622       break;
623     default:
624       NOT_REACHED ();
625     }
626 }
627
628 /*
629   Covariance computed without dividing by the sample size.
630  */
631 static gsl_matrix *
632 covariance_calculate_double_pass_unnormalized (struct covariance *cov)
633 {
634   return  cm_to_gsl (cov);
635 }
636
637 static gsl_matrix *
638 covariance_calculate_single_pass_unnormalized (struct covariance *cov)
639 {
640   size_t i, j;
641
642   for (i = 0 ; i < cov->dim; ++i)
643     {
644       for (j = 0 ; j < cov->dim; ++j)
645         {
646           double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
647           *x -= pow2 (gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
648             / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
649         }
650     }
651   for ( j = 0 ; j < cov->dim - 1; ++j)
652     {
653       for (i = j + 1 ; i < cov->dim; ++i)
654         {
655           double *x = &cov->cm [cm_idx (cov, i, j)];
656           
657           *x -=
658             gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j) 
659             *
660             gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i) 
661           / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
662         }
663     }
664
665   return cm_to_gsl (cov);
666 }
667
668
669 /* Return a pointer to gsl_matrix containing the pairwise covariances.  The
670    caller owns the returned matrix and must free it when it is no longer
671    needed.
672
673    Call this function only after all data have been accumulated.  */
674 gsl_matrix *
675 covariance_calculate_unnormalized (struct covariance *cov)
676 {
677   if ( cov->state <= 0 )
678     return NULL;
679
680   switch (cov->passes)
681     {
682     case 1:
683       return covariance_calculate_single_pass_unnormalized (cov);  
684       break;
685     case 2:
686       return covariance_calculate_double_pass_unnormalized (cov);  
687       break;
688     default:
689       NOT_REACHED ();
690     }
691 }
692
693 /* Function to access the categoricals used by COV
694    The return value is owned by the COV
695 */
696 const struct categoricals *
697 covariance_get_categoricals (const struct covariance *cov)
698 {
699   return cov->categoricals;
700 }
701
702
703 /* Destroy the COV object */
704 void
705 covariance_destroy (struct covariance *cov)
706 {
707   size_t i;
708
709   categoricals_destroy (cov->categoricals);
710
711   for (i = 0; i < n_MOMENTS; ++i)
712     gsl_matrix_free (cov->moments[i]);
713
714   free (cov->moments);
715   free (cov->cm);
716   free (cov);
717 }
718
719 size_t
720 covariance_dim (const struct covariance * cov)
721 {
722   return (cov->dim);
723 }
724
725 \f
726
727 /*
728   Routines to assist debugging.
729   The following are not thoroughly tested and in certain respects
730   unreliable.  They should only be
731   used for aids to development. Not as user accessible code.
732 */
733
734 #include "libpspp/str.h"
735 #include "output/tab.h"
736 #include "data/format.h"
737
738
739 /* Create a table which can be populated with the encodings for
740    the covariance matrix COV */
741 struct tab_table *
742 covariance_dump_enc_header (const struct covariance *cov, int length)
743 {
744   struct tab_table *t = tab_create (cov->dim,  length);
745   int n;
746   int i;
747
748   tab_title (t, "Covariance Encoding");
749
750   tab_box (t, 
751            TAL_2, TAL_2, 0, 0,
752            0, 0,   tab_nc (t) - 1,   tab_nr (t) - 1);
753
754   tab_hline (t, TAL_2, 0, tab_nc (t) - 1, 1);
755
756
757   for (i = 0 ; i < cov->n_vars; ++i)
758     {
759       tab_text (t, i, 0, TAT_TITLE, var_get_name (cov->vars[i]));
760       tab_vline (t, TAL_1, i + 1, 0, tab_nr (t) - 1);
761     }
762
763   n = 0;
764   while (i < cov->dim)
765     {
766       struct string str;
767       int idx = i - cov->n_vars;
768       const struct interaction *iact =
769         categoricals_get_interaction_by_subscript (cov->categoricals, idx);
770       int df;
771
772       ds_init_empty (&str);
773       interaction_to_string (iact, &str);
774
775       df = categoricals_df (cov->categoricals, n);
776
777       tab_joint_text (t,
778                       i, 0,
779                       i + df - 1, 0,
780                       TAT_TITLE, ds_cstr (&str));
781
782       if (i + df < tab_nr (t) - 1)
783         tab_vline (t, TAL_1, i + df, 0, tab_nr (t) - 1);
784
785       i += df;
786       n++;
787       ds_destroy (&str);
788     }
789
790   return t;
791 }
792
793
794 /*
795   Append table T, which should have been returned by covariance_dump_enc_header
796   with an entry corresponding to case C for the covariance matrix COV
797  */
798 void
799 covariance_dump_enc (const struct covariance *cov, const struct ccase *c,
800                      struct tab_table *t)
801 {
802   static int row = 0;
803   int i;
804   ++row;
805   for (i = 0 ; i < cov->dim; ++i)
806     {
807       double v = get_val (cov, i, c);
808       tab_double (t, i, row, 0, v, i < cov->n_vars ? NULL : &F_8_0);
809     }
810 }