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