ONEWAY: Fix printf format specifier.
[pspp] / src / language / stats / factor.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 <gsl/gsl_vector.h>
20 #include <gsl/gsl_linalg.h>
21 #include <gsl/gsl_matrix.h>
22 #include <gsl/gsl_eigen.h> 
23 #include <gsl/gsl_blas.h> 
24 #include <gsl/gsl_sort_vector.h>
25 #include <gsl/gsl_cdf.h>
26
27 #include "data/casegrouper.h"
28 #include "data/casereader.h"
29 #include "data/casewriter.h"
30 #include "data/dataset.h"
31 #include "data/dictionary.h"
32 #include "data/format.h"
33 #include "data/subcase.h"
34 #include "language/command.h"
35 #include "language/lexer/lexer.h"
36 #include "language/lexer/value-parser.h"
37 #include "language/lexer/variable-parser.h"
38 #include "libpspp/message.h"
39 #include "libpspp/misc.h"
40 #include "math/correlation.h"
41 #include "math/covariance.h"
42 #include "math/moments.h"
43 #include "output/chart-item.h"
44 #include "output/charts/scree.h"
45 #include "output/tab.h"
46
47 #include "gettext.h"
48 #define _(msgid) gettext (msgid)
49 #define N_(msgid) msgid
50
51 enum method
52   {
53     METHOD_CORR,
54     METHOD_COV
55   };
56
57 enum missing_type
58   {
59     MISS_LISTWISE,
60     MISS_PAIRWISE,
61     MISS_MEANSUB,
62   };
63
64 enum extraction_method
65   {
66     EXTRACTION_PC,
67     EXTRACTION_PAF,
68   };
69
70 enum plot_opts
71   {
72     PLOT_SCREE = 0x0001,
73     PLOT_ROTATION = 0x0002
74   };
75
76 enum print_opts
77   {
78     PRINT_UNIVARIATE  = 0x0001,
79     PRINT_DETERMINANT = 0x0002,
80     PRINT_INV         = 0x0004,
81     PRINT_AIC         = 0x0008,
82     PRINT_SIG         = 0x0010,
83     PRINT_COVARIANCE  = 0x0020,
84     PRINT_CORRELATION = 0x0040,
85     PRINT_ROTATION    = 0x0080,
86     PRINT_EXTRACTION  = 0x0100,
87     PRINT_INITIAL     = 0x0200,
88     PRINT_KMO         = 0x0400,
89     PRINT_REPR        = 0x0800, 
90     PRINT_FSCORE      = 0x1000
91   };
92
93 enum rotation_type
94   {
95     ROT_VARIMAX = 0,
96     ROT_EQUAMAX,
97     ROT_QUARTIMAX,
98     ROT_NONE
99   };
100
101 typedef void (*rotation_coefficients) (double *x, double *y,
102                                     double a, double b, double c, double d,
103                                     const gsl_matrix *loadings );
104
105
106 static void
107 varimax_coefficients (double *x, double *y,
108                       double a, double b, double c, double d,
109                       const gsl_matrix *loadings )
110 {
111   *x = d - 2 * a * b / loadings->size1;
112   *y = c - (a * a - b * b) / loadings->size1;
113 }
114
115 static void
116 equamax_coefficients (double *x, double *y,
117                       double a, double b, double c, double d,
118                       const gsl_matrix *loadings )
119 {
120   *x = d - loadings->size2 * a * b / loadings->size1;
121   *y = c - loadings->size2 * (a * a - b * b) / (2 * loadings->size1);
122 }
123
124 static void
125 quartimax_coefficients (double *x, double *y,
126                       double a UNUSED, double b UNUSED, double c, double d,
127                       const gsl_matrix *loadings UNUSED)
128 {
129   *x = d ;
130   *y = c ;
131 }
132
133 static const rotation_coefficients rotation_coeff[3] = {
134   varimax_coefficients,
135   equamax_coefficients,
136   quartimax_coefficients
137 };
138
139
140 struct cmd_factor 
141 {
142   size_t n_vars;
143   const struct variable **vars;
144
145   const struct variable *wv;
146
147   enum method method;
148   enum missing_type missing_type;
149   enum mv_class exclude;
150   enum print_opts print;
151   enum extraction_method extraction;
152   enum plot_opts plot;
153   enum rotation_type rotation;
154
155   /* Extraction Criteria */
156   int n_factors;
157   double min_eigen;
158   double econverge;
159   int iterations;
160
161   double rconverge;
162
163   /* Format */
164   double blank;
165   bool sort;
166 };
167
168 struct idata
169 {
170   /* Intermediate values used in calculation */
171
172   const gsl_matrix *corr ;  /* The correlation matrix */
173   gsl_matrix *cov ;         /* The covariance matrix */
174   const gsl_matrix *n ;     /* Matrix of number of samples */
175
176   gsl_vector *eval ;  /* The eigenvalues */
177   gsl_matrix *evec ;  /* The eigenvectors */
178
179   int n_extractions;
180
181   gsl_vector *msr ;  /* Multiple Squared Regressions */
182
183   double detR;  /* The determinant of the correlation matrix */
184 };
185
186 static struct idata *
187 idata_alloc (size_t n_vars)
188 {
189   struct idata *id = xzalloc (sizeof (*id));
190
191   id->n_extractions = 0;
192   id->msr = gsl_vector_alloc (n_vars);
193
194   id->eval = gsl_vector_alloc (n_vars);
195   id->evec = gsl_matrix_alloc (n_vars, n_vars);
196
197   return id;
198 }
199
200 static void
201 idata_free (struct idata *id)
202 {
203   gsl_vector_free (id->msr);
204   gsl_vector_free (id->eval);
205   gsl_matrix_free (id->evec);
206   if (id->cov != NULL)
207     gsl_matrix_free (id->cov);
208
209   free (id);
210 }
211
212
213 static gsl_matrix *
214 anti_image (const gsl_matrix *m)
215 {
216   int i, j;
217   gsl_matrix *a;
218   assert (m->size1 == m->size2);
219
220   a = gsl_matrix_alloc (m->size1, m->size2);
221   
222   for (i = 0; i < m->size1; ++i)
223     {
224       for (j = 0; j < m->size2; ++j)
225         {
226           double *p = gsl_matrix_ptr (a, i, j);
227           *p = gsl_matrix_get (m, i, j);
228           *p /= gsl_matrix_get (m, i, i);
229           *p /= gsl_matrix_get (m, j, j);
230         }
231     }
232
233   return a;
234 }
235
236
237 /* Return the sum of all the elements excluding row N */
238 static double
239 ssq_od_n (const gsl_matrix *m, int n)
240 {
241   int i, j;
242   double ss = 0;
243   assert (m->size1 == m->size2);
244
245   assert (n < m->size1);
246   
247   for (i = 0; i < m->size1; ++i)
248     {
249       if (i == n ) continue;
250       for (j = 0; j < m->size2; ++j)
251         {
252           ss += pow2 (gsl_matrix_get (m, i, j));
253         }
254     }
255
256   return ss;
257 }
258
259
260
261 #if 0
262 static void
263 dump_matrix (const gsl_matrix *m)
264 {
265   size_t i, j;
266
267   for (i = 0 ; i < m->size1; ++i)
268     {
269       for (j = 0 ; j < m->size2; ++j)
270         printf ("%02f ", gsl_matrix_get (m, i, j));
271       printf ("\n");
272     }
273 }
274
275 static void
276 dump_matrix_permute (const gsl_matrix *m, const gsl_permutation *p)
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, gsl_permutation_get (p, i), j));
284       printf ("\n");
285     }
286 }
287
288
289 static void
290 dump_vector (const gsl_vector *v)
291 {
292   size_t i;
293   for (i = 0 ; i < v->size; ++i)
294     {
295       printf ("%02f\n", gsl_vector_get (v, i));
296     }
297   printf ("\n");
298 }
299 #endif
300
301
302 static int 
303 n_extracted_factors (const struct cmd_factor *factor, struct idata *idata)
304 {
305   int i;
306   
307   /* If there is a cached value, then return that. */
308   if ( idata->n_extractions != 0)
309     return idata->n_extractions;
310
311   /* Otherwise, if the number of factors has been explicitly requested,
312      use that. */
313   if (factor->n_factors > 0)
314     {
315       idata->n_extractions = factor->n_factors;
316       goto finish;
317     }
318   
319   /* Use the MIN_EIGEN setting. */
320   for (i = 0 ; i < idata->eval->size; ++i)
321     {
322       double evali = fabs (gsl_vector_get (idata->eval, i));
323
324       idata->n_extractions = i;
325
326       if (evali < factor->min_eigen)
327         goto finish;
328     }
329
330  finish:
331   return idata->n_extractions;
332 }
333
334
335 /* Returns a newly allocated matrix identical to M.
336    It it the callers responsibility to free the returned value.
337 */
338 static gsl_matrix *
339 matrix_dup (const gsl_matrix *m)
340 {
341   gsl_matrix *n =  gsl_matrix_alloc (m->size1, m->size2);
342
343   gsl_matrix_memcpy (n, m);
344
345   return n;
346 }
347
348
349 struct smr_workspace
350 {
351   /* Copy of the subject */
352   gsl_matrix *m;
353   
354   gsl_matrix *inverse;
355
356   gsl_permutation *perm;
357
358   gsl_matrix *result1;
359   gsl_matrix *result2;
360 };
361
362
363 static struct smr_workspace *ws_create (const gsl_matrix *input)
364 {
365   struct smr_workspace *ws = xmalloc (sizeof (*ws));
366   
367   ws->m = gsl_matrix_alloc (input->size1, input->size2);
368   ws->inverse = gsl_matrix_calloc (input->size1 - 1, input->size2 - 1);
369   ws->perm = gsl_permutation_alloc (input->size1 - 1);
370   ws->result1 = gsl_matrix_calloc (input->size1 - 1, 1);
371   ws->result2 = gsl_matrix_calloc (1, 1);
372
373   return ws;
374 }
375
376 static void
377 ws_destroy (struct smr_workspace *ws)
378 {
379   gsl_matrix_free (ws->result2);
380   gsl_matrix_free (ws->result1);
381   gsl_permutation_free (ws->perm);
382   gsl_matrix_free (ws->inverse);
383   gsl_matrix_free (ws->m);
384
385   free (ws);
386 }
387
388
389 /* 
390    Return the square of the regression coefficient for VAR regressed against all other variables.
391  */
392 static double
393 squared_multiple_correlation (const gsl_matrix *corr, int var, struct smr_workspace *ws)
394 {
395   /* For an explanation of what this is doing, see 
396      http://www.visualstatistics.net/Visual%20Statistics%20Multimedia/multiple_regression_analysis.htm
397   */
398
399   int signum = 0;
400   gsl_matrix_view rxx;
401
402   gsl_matrix_memcpy (ws->m, corr);
403
404   gsl_matrix_swap_rows (ws->m, 0, var);
405   gsl_matrix_swap_columns (ws->m, 0, var);
406
407   rxx = gsl_matrix_submatrix (ws->m, 1, 1, ws->m->size1 - 1, ws->m->size1 - 1); 
408
409   gsl_linalg_LU_decomp (&rxx.matrix, ws->perm, &signum);
410
411   gsl_linalg_LU_invert (&rxx.matrix, ws->perm, ws->inverse);
412
413   {
414     gsl_matrix_const_view rxy = gsl_matrix_const_submatrix (ws->m, 1, 0, ws->m->size1 - 1, 1);
415     gsl_matrix_const_view ryx = gsl_matrix_const_submatrix (ws->m, 0, 1, 1, ws->m->size1 - 1);
416
417     gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans,
418                     1.0, ws->inverse, &rxy.matrix, 0.0, ws->result1);
419
420     gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans,
421                     1.0, &ryx.matrix, ws->result1, 0.0, ws->result2);
422   }
423
424   return gsl_matrix_get (ws->result2, 0, 0);
425 }
426
427
428
429 static double the_communality (const gsl_matrix *evec, const gsl_vector *eval, int n, int n_factors);
430
431
432 struct factor_matrix_workspace
433 {
434   size_t n_factors;
435   gsl_eigen_symmv_workspace *eigen_ws;
436
437   gsl_vector *eval ;
438   gsl_matrix *evec ;
439
440   gsl_matrix *gamma ;
441
442   gsl_matrix *r;
443 };
444
445 static struct factor_matrix_workspace *
446 factor_matrix_workspace_alloc (size_t n, size_t nf)
447 {
448   struct factor_matrix_workspace *ws = xmalloc (sizeof (*ws));
449
450   ws->n_factors = nf;
451   ws->gamma = gsl_matrix_calloc (nf, nf);
452   ws->eigen_ws = gsl_eigen_symmv_alloc (n);
453   ws->eval = gsl_vector_alloc (n);
454   ws->evec = gsl_matrix_alloc (n, n);
455   ws->r  = gsl_matrix_alloc (n, n);
456   
457   return ws;
458 }
459
460 static void
461 factor_matrix_workspace_free (struct factor_matrix_workspace *ws)
462 {
463   gsl_eigen_symmv_free (ws->eigen_ws);
464   gsl_vector_free (ws->eval);
465   gsl_matrix_free (ws->evec);
466   gsl_matrix_free (ws->gamma);
467   gsl_matrix_free (ws->r);
468   free (ws);
469 }
470
471 /*
472   Shift P left by OFFSET places, and overwrite TARGET
473   with the shifted result.
474   Positions in TARGET less than OFFSET are unchanged.
475 */
476 static void
477 perm_shift_apply (gsl_permutation *target, const gsl_permutation *p,
478                   size_t offset)
479 {
480   size_t i;
481   assert (target->size == p->size);
482   assert (offset <= target->size);
483
484   for (i = 0; i < target->size - offset; ++i)
485     {
486       target->data[i] = p->data [i + offset];
487     }
488 }
489
490
491 /* 
492    Indirectly sort the rows of matrix INPUT, storing the sort order in PERM.
493    The sort criteria are as follows:
494    
495    Rows are sorted on the first column, until the absolute value of an
496    element in a subsequent column  is greater than that of the first
497    column.  Thereafter, rows will be sorted on the second column,
498    until the absolute value of an element in a subsequent column
499    exceeds that of the second column ...
500 */
501 static void
502 sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm)
503 {
504   const size_t n = perm->size;
505   const size_t m = input->size2;
506   int i, j;
507   gsl_matrix *mat ;
508   int column_n = 0;
509   int row_n = 0;
510   gsl_permutation *p;
511
512   assert (perm->size == input->size1);
513
514   p = gsl_permutation_alloc (n);
515
516   /* Copy INPUT into MAT, discarding the sign */
517   mat = gsl_matrix_alloc (n, m);
518   for (i = 0 ; i < mat->size1; ++i)
519     {
520       for (j = 0 ; j < mat->size2; ++j)
521         {
522           double x = gsl_matrix_get (input, i, j);
523           gsl_matrix_set (mat, i, j, fabs (x));
524         }
525     }
526
527   while (column_n < m && row_n < n) 
528     {
529       gsl_vector_const_view columni = gsl_matrix_const_column (mat, column_n);
530       gsl_sort_vector_index (p, &columni.vector);
531
532       for (i = 0 ; i < n; ++i)
533         {
534           gsl_vector_view row = gsl_matrix_row (mat, p->data[n - 1 - i]);
535           size_t maxindex = gsl_vector_max_index (&row.vector);
536           
537           if ( maxindex > column_n )
538             break;
539
540           /* All subsequent elements of this row, are of no interest.
541              So set them all to a highly negative value */
542           for (j = column_n + 1; j < row.vector.size ; ++j)
543             gsl_vector_set (&row.vector, j, -DBL_MAX);
544         }
545
546       perm_shift_apply (perm, p, row_n);
547       row_n += i;
548
549       column_n++;
550     }
551
552   gsl_permutation_free (p);
553   gsl_matrix_free (mat);
554   
555   assert ( 0 == gsl_permutation_valid (perm));
556
557   /* We want the biggest value to be first */
558   gsl_permutation_reverse (perm);    
559 }
560
561
562 static void
563 drot_go (double phi, double *l0, double *l1)
564 {
565   double r0 = cos (phi) * *l0 + sin (phi) * *l1;
566   double r1 = - sin (phi) * *l0 + cos (phi) * *l1;
567
568   *l0 = r0;
569   *l1 = r1;
570 }
571
572
573 static gsl_matrix *
574 clone_matrix (const gsl_matrix *m)
575 {
576   int j, k;
577   gsl_matrix *c = gsl_matrix_calloc (m->size1, m->size2);
578
579   for (j = 0 ; j < c->size1; ++j)
580     {
581       for (k = 0 ; k < c->size2; ++k)
582         {
583           const double *v = gsl_matrix_const_ptr (m, j, k);
584           gsl_matrix_set (c, j, k, *v);
585         }
586     }
587
588   return c;
589 }
590
591
592 static double 
593 initial_sv (const gsl_matrix *fm)
594 {
595   int j, k;
596
597   double sv = 0.0;
598   for (j = 0 ; j < fm->size2; ++j)
599     {
600       double l4s = 0;
601       double l2s = 0;
602
603       for (k = j + 1 ; k < fm->size2; ++k)
604         {
605           double lambda = gsl_matrix_get (fm, k, j);
606           double lambda_sq = lambda * lambda;
607           double lambda_4 = lambda_sq * lambda_sq;
608
609           l4s += lambda_4;
610           l2s += lambda_sq;
611         }
612       sv += ( fm->size1 * l4s - (l2s * l2s) ) / (fm->size1 * fm->size1 );
613     }
614   return sv;
615 }
616
617 static void
618 rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
619         const gsl_vector *communalities,
620         gsl_matrix *result,
621         gsl_vector *rotated_loadings
622         )
623 {
624   int j, k;
625   int i;
626   double prev_sv;
627
628   /* First get a normalised version of UNROT */
629   gsl_matrix *normalised = gsl_matrix_calloc (unrot->size1, unrot->size2);
630   gsl_matrix *h_sqrt = gsl_matrix_calloc (communalities->size, communalities->size);
631   gsl_matrix *h_sqrt_inv ;
632
633   /* H is the diagonal matrix containing the absolute values of the communalities */
634   for (i = 0 ; i < communalities->size ; ++i)
635     {
636       double *ptr = gsl_matrix_ptr (h_sqrt, i, i);
637       *ptr = fabs (gsl_vector_get (communalities, i));
638     }
639
640   /* Take the square root of the communalities */
641   gsl_linalg_cholesky_decomp (h_sqrt);
642
643
644   /* Save a copy of h_sqrt and invert it */
645   h_sqrt_inv = clone_matrix (h_sqrt);
646   gsl_linalg_cholesky_decomp (h_sqrt_inv);
647   gsl_linalg_cholesky_invert (h_sqrt_inv);
648
649   /* normalised vertion is H^{1/2} x UNROT */
650   gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans, 1.0, h_sqrt_inv, unrot, 0.0, normalised);
651
652   gsl_matrix_free (h_sqrt_inv);
653
654
655   /* Now perform the rotation iterations */
656
657   prev_sv = initial_sv (normalised);
658   for (i = 0 ; i < cf->iterations ; ++i)
659     {
660       double sv = 0.0;
661       for (j = 0 ; j < normalised->size2; ++j)
662         {
663           /* These variables relate to the convergence criterium */
664           double l4s = 0;
665           double l2s = 0;
666
667           for (k = j + 1 ; k < normalised->size2; ++k)
668             {
669               int p;
670               double a = 0.0;
671               double b = 0.0;
672               double c = 0.0;
673               double d = 0.0;
674               double x, y;
675               double phi;
676
677               for (p = 0; p < normalised->size1; ++p)
678                 {
679                   double jv = gsl_matrix_get (normalised, p, j);
680                   double kv = gsl_matrix_get (normalised, p, k);
681               
682                   double u = jv * jv - kv * kv;
683                   double v = 2 * jv * kv;
684                   a += u;
685                   b += v;
686                   c +=  u * u - v * v;
687                   d += 2 * u * v;
688                 }
689
690               rotation_coeff [cf->rotation] (&x, &y, a, b, c, d, normalised);
691
692               phi = atan2 (x,  y) / 4.0 ;
693
694               /* Don't bother rotating if the angle is small */
695               if ( fabs (sin (phi) ) <= pow (10.0, -15.0))
696                   continue;
697
698               for (p = 0; p < normalised->size1; ++p)
699                 {
700                   double *lambda0 = gsl_matrix_ptr (normalised, p, j);
701                   double *lambda1 = gsl_matrix_ptr (normalised, p, k);
702                   drot_go (phi, lambda0, lambda1);
703                 }
704
705               /* Calculate the convergence criterium */
706               {
707                 double lambda = gsl_matrix_get (normalised, k, j);
708                 double lambda_sq = lambda * lambda;
709                 double lambda_4 = lambda_sq * lambda_sq;
710
711                 l4s += lambda_4;
712                 l2s += lambda_sq;
713               }
714             }
715           sv += ( normalised->size1 * l4s - (l2s * l2s) ) / (normalised->size1 * normalised->size1 );
716         }
717
718       if ( fabs (sv - prev_sv) <= cf->rconverge)
719         break;
720
721       prev_sv = sv;
722     }
723
724   gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans, 1.0,
725                   h_sqrt, normalised,  0.0,   result);
726
727   gsl_matrix_free (h_sqrt);
728
729
730   /* reflect negative sums and populate the rotated loadings vector*/
731   for (i = 0 ; i < result->size2; ++i)
732     {
733       double ssq = 0.0;
734       double sum = 0.0;
735       for (j = 0 ; j < result->size1; ++j)
736         {
737           double s = gsl_matrix_get (result, j, i);
738           ssq += s * s;
739           sum += gsl_matrix_get (result, j, i);
740         }
741
742       gsl_vector_set (rotated_loadings, i, ssq);
743
744       if ( sum < 0 )
745         for (j = 0 ; j < result->size1; ++j)
746           {
747             double *lambda = gsl_matrix_ptr (result, j, i);
748             *lambda = - *lambda;
749           }
750     }
751 }
752
753
754 /*
755   Get an approximation for the factor matrix into FACTORS, and the communalities into COMMUNALITIES.
756   R is the matrix to be analysed.
757   WS is a pointer to a structure which must have been initialised with factor_matrix_workspace_init.
758  */
759 static void
760 iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matrix *factors, 
761                        struct factor_matrix_workspace *ws)
762 {
763   size_t i;
764   gsl_matrix_view mv ;
765
766   assert (r->size1 == r->size2);
767   assert (r->size1 == communalities->size);
768
769   assert (factors->size1 == r->size1);
770   assert (factors->size2 == ws->n_factors);
771
772   gsl_matrix_memcpy (ws->r, r);
773
774   /* Apply Communalities to diagonal of correlation matrix */
775   for (i = 0 ; i < communalities->size ; ++i)
776     {
777       double *x = gsl_matrix_ptr (ws->r, i, i);
778       *x = gsl_vector_get (communalities, i);
779     }
780
781   gsl_eigen_symmv (ws->r, ws->eval, ws->evec, ws->eigen_ws);
782
783   mv = gsl_matrix_submatrix (ws->evec, 0, 0, ws->evec->size1, ws->n_factors);
784
785   /* Gamma is the diagonal matrix containing the absolute values of the eigenvalues */
786   for (i = 0 ; i < ws->n_factors ; ++i)
787     {
788       double *ptr = gsl_matrix_ptr (ws->gamma, i, i);
789       *ptr = fabs (gsl_vector_get (ws->eval, i));
790     }
791
792   /* Take the square root of gamma */
793   gsl_linalg_cholesky_decomp (ws->gamma);
794
795   gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans, 1.0, &mv.matrix, ws->gamma, 0.0, factors);
796
797   for (i = 0 ; i < r->size1 ; ++i)
798     {
799       double h = the_communality (ws->evec, ws->eval, i, ws->n_factors);
800       gsl_vector_set (communalities, i, h);
801     }
802 }
803
804
805
806 static bool run_factor (struct dataset *ds, const struct cmd_factor *factor);
807
808
809 int
810 cmd_factor (struct lexer *lexer, struct dataset *ds)
811 {
812   bool extraction_seen = false;
813   const struct dictionary *dict = dataset_dict (ds);
814
815   struct cmd_factor factor;
816   factor.n_vars = 0;
817   factor.vars = NULL;
818   factor.method = METHOD_CORR;
819   factor.missing_type = MISS_LISTWISE;
820   factor.exclude = MV_ANY;
821   factor.print = PRINT_INITIAL | PRINT_EXTRACTION | PRINT_ROTATION;
822   factor.extraction = EXTRACTION_PC;
823   factor.n_factors = 0;
824   factor.min_eigen = SYSMIS;
825   factor.iterations = 25;
826   factor.econverge = 0.001;
827
828   factor.blank = 0;
829   factor.sort = false;
830   factor.plot = 0;
831   factor.rotation = ROT_VARIMAX;
832
833   factor.rconverge = 0.0001;
834
835   factor.wv = dict_get_weight (dict);
836
837   lex_match (lexer, T_SLASH);
838
839   if (!lex_force_match_id (lexer, "VARIABLES"))
840     {
841       goto error;
842     }
843
844   lex_match (lexer, T_EQUALS);
845
846   if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars,
847                               PV_NO_DUPLICATE | PV_NUMERIC))
848     goto error;
849
850   if (factor.n_vars < 2)
851     msg (MW, _("Factor analysis on a single variable is not useful."));
852
853   while (lex_token (lexer) != T_ENDCMD)
854     {
855       lex_match (lexer, T_SLASH);
856
857       if (lex_match_id (lexer, "PLOT"))
858         {
859           lex_match (lexer, T_EQUALS);
860           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
861             {
862               if (lex_match_id (lexer, "EIGEN"))
863                 {
864                   factor.plot |= PLOT_SCREE;
865                 }
866 #if FACTOR_FULLY_IMPLEMENTED
867               else if (lex_match_id (lexer, "ROTATION"))
868                 {
869                 }
870 #endif
871               else
872                 {
873                   lex_error (lexer, NULL);
874                   goto error;
875                 }
876             }
877         }
878       else if (lex_match_id (lexer, "METHOD"))
879         {
880           lex_match (lexer, T_EQUALS);
881           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
882             {
883               if (lex_match_id (lexer, "COVARIANCE"))
884                 {
885                   factor.method = METHOD_COV;
886                 }
887               else if (lex_match_id (lexer, "CORRELATION"))
888                 {
889                   factor.method = METHOD_CORR;
890                 }
891               else
892                 {
893                   lex_error (lexer, NULL);
894                   goto error;
895                 }
896             }
897         }
898       else if (lex_match_id (lexer, "ROTATION"))
899         {
900           lex_match (lexer, T_EQUALS);
901           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
902             {
903               /* VARIMAX and DEFAULT are defaults */
904               if (lex_match_id (lexer, "VARIMAX") || lex_match_id (lexer, "DEFAULT"))
905                 {
906                   factor.rotation = ROT_VARIMAX;
907                 }
908               else if (lex_match_id (lexer, "EQUAMAX"))
909                 {
910                   factor.rotation = ROT_EQUAMAX;
911                 }
912               else if (lex_match_id (lexer, "QUARTIMAX"))
913                 {
914                   factor.rotation = ROT_QUARTIMAX;
915                 }
916               else if (lex_match_id (lexer, "NOROTATE"))
917                 {
918                   factor.rotation = ROT_NONE;
919                 }
920               else
921                 {
922                   lex_error (lexer, NULL);
923                   goto error;
924                 }
925             }
926         }
927       else if (lex_match_id (lexer, "CRITERIA"))
928         {
929           lex_match (lexer, T_EQUALS);
930           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
931             {
932               if (lex_match_id (lexer, "FACTORS"))
933                 {
934                   if ( lex_force_match (lexer, T_LPAREN))
935                     {
936                       lex_force_int (lexer);
937                       factor.n_factors = lex_integer (lexer);
938                       lex_get (lexer);
939                       lex_force_match (lexer, T_RPAREN);
940                     }
941                 }
942               else if (lex_match_id (lexer, "MINEIGEN"))
943                 {
944                   if ( lex_force_match (lexer, T_LPAREN))
945                     {
946                       lex_force_num (lexer);
947                       factor.min_eigen = lex_number (lexer);
948                       lex_get (lexer);
949                       lex_force_match (lexer, T_RPAREN);
950                     }
951                 }
952               else if (lex_match_id (lexer, "ECONVERGE"))
953                 {
954                   if ( lex_force_match (lexer, T_LPAREN))
955                     {
956                       lex_force_num (lexer);
957                       factor.econverge = lex_number (lexer);
958                       lex_get (lexer);
959                       lex_force_match (lexer, T_RPAREN);
960                     }
961                 }
962               else if (lex_match_id (lexer, "RCONVERGE"))
963                 {
964                   if ( lex_force_match (lexer, T_LPAREN))
965                     {
966                       lex_force_num (lexer);
967                       factor.rconverge = lex_number (lexer);
968                       lex_get (lexer);
969                       lex_force_match (lexer, T_RPAREN);
970                     }
971                 }
972               else if (lex_match_id (lexer, "ITERATE"))
973                 {
974                   if ( lex_force_match (lexer, T_LPAREN))
975                     {
976                       lex_force_int (lexer);
977                       factor.iterations = lex_integer (lexer);
978                       lex_get (lexer);
979                       lex_force_match (lexer, T_RPAREN);
980                     }
981                 }
982               else if (lex_match_id (lexer, "DEFAULT"))
983                 {
984                   factor.n_factors = 0;
985                   factor.min_eigen = 1;
986                   factor.iterations = 25;
987                 }
988               else
989                 {
990                   lex_error (lexer, NULL);
991                   goto error;
992                 }
993             }
994         }
995       else if (lex_match_id (lexer, "EXTRACTION"))
996         {
997           extraction_seen = true;
998           lex_match (lexer, T_EQUALS);
999           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1000             {
1001               if (lex_match_id (lexer, "PAF"))
1002                 {
1003                   factor.extraction = EXTRACTION_PAF;
1004                 }
1005               else if (lex_match_id (lexer, "PC"))
1006                 {
1007                   factor.extraction = EXTRACTION_PC;
1008                 }
1009               else if (lex_match_id (lexer, "PA1"))
1010                 {
1011                   factor.extraction = EXTRACTION_PC;
1012                 }
1013               else if (lex_match_id (lexer, "DEFAULT"))
1014                 {
1015                   factor.extraction = EXTRACTION_PC;
1016                 }
1017               else
1018                 {
1019                   lex_error (lexer, NULL);
1020                   goto error;
1021                 }
1022             }
1023         }
1024       else if (lex_match_id (lexer, "FORMAT"))
1025         {
1026           lex_match (lexer, T_EQUALS);
1027           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1028             {
1029               if (lex_match_id (lexer, "SORT"))
1030                 {
1031                   factor.sort = true;
1032                 }
1033               else if (lex_match_id (lexer, "BLANK"))
1034                 {
1035                   if ( lex_force_match (lexer, T_LPAREN))
1036                     {
1037                       lex_force_num (lexer);
1038                       factor.blank = lex_number (lexer);
1039                       lex_get (lexer);
1040                       lex_force_match (lexer, T_RPAREN);
1041                     }
1042                 }
1043               else if (lex_match_id (lexer, "DEFAULT"))
1044                 {
1045                   factor.blank = 0;
1046                   factor.sort = false;
1047                 }
1048               else
1049                 {
1050                   lex_error (lexer, NULL);
1051                   goto error;
1052                 }
1053             }
1054         }
1055       else if (lex_match_id (lexer, "PRINT"))
1056         {
1057           factor.print = 0;
1058           lex_match (lexer, T_EQUALS);
1059           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1060             {
1061               if (lex_match_id (lexer, "UNIVARIATE"))
1062                 {
1063                   factor.print |= PRINT_UNIVARIATE;
1064                 }
1065               else if (lex_match_id (lexer, "DET"))
1066                 {
1067                   factor.print |= PRINT_DETERMINANT;
1068                 }
1069 #if FACTOR_FULLY_IMPLEMENTED
1070               else if (lex_match_id (lexer, "INV"))
1071                 {
1072                 }
1073               else if (lex_match_id (lexer, "AIC"))
1074                 {
1075                 }
1076 #endif
1077               else if (lex_match_id (lexer, "SIG"))
1078                 {
1079                   factor.print |= PRINT_SIG;
1080                 }
1081               else if (lex_match_id (lexer, "CORRELATION"))
1082                 {
1083                   factor.print |= PRINT_CORRELATION;
1084                 }
1085 #if FACTOR_FULLY_IMPLEMENTED
1086               else if (lex_match_id (lexer, "COVARIANCE"))
1087                 {
1088                 }
1089 #endif
1090               else if (lex_match_id (lexer, "ROTATION"))
1091                 {
1092                   factor.print |= PRINT_ROTATION;
1093                 }
1094               else if (lex_match_id (lexer, "EXTRACTION"))
1095                 {
1096                   factor.print |= PRINT_EXTRACTION;
1097                 }
1098               else if (lex_match_id (lexer, "INITIAL"))
1099                 {
1100                   factor.print |= PRINT_INITIAL;
1101                 }
1102               else if (lex_match_id (lexer, "KMO"))
1103                 {
1104                   factor.print |= PRINT_KMO;
1105                 }
1106 #if FACTOR_FULLY_IMPLEMENTED
1107               else if (lex_match_id (lexer, "REPR"))
1108                 {
1109                 }
1110               else if (lex_match_id (lexer, "FSCORE"))
1111                 {
1112                 }
1113 #endif
1114               else if (lex_match (lexer, T_ALL))
1115                 {
1116                   factor.print = 0xFFFF;
1117                 }
1118               else if (lex_match_id (lexer, "DEFAULT"))
1119                 {
1120                   factor.print |= PRINT_INITIAL ;
1121                   factor.print |= PRINT_EXTRACTION ;
1122                   factor.print |= PRINT_ROTATION ;
1123                 }
1124               else
1125                 {
1126                   lex_error (lexer, NULL);
1127                   goto error;
1128                 }
1129             }
1130         }
1131       else if (lex_match_id (lexer, "MISSING"))
1132         {
1133           lex_match (lexer, T_EQUALS);
1134           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1135             {
1136               if (lex_match_id (lexer, "INCLUDE"))
1137                 {
1138                   factor.exclude = MV_SYSTEM;
1139                 }
1140               else if (lex_match_id (lexer, "EXCLUDE"))
1141                 {
1142                   factor.exclude = MV_ANY;
1143                 }
1144               else if (lex_match_id (lexer, "LISTWISE"))
1145                 {
1146                   factor.missing_type = MISS_LISTWISE;
1147                 }
1148               else if (lex_match_id (lexer, "PAIRWISE"))
1149                 {
1150                   factor.missing_type = MISS_PAIRWISE;
1151                 }
1152               else if (lex_match_id (lexer, "MEANSUB"))
1153                 {
1154                   factor.missing_type = MISS_MEANSUB;
1155                 }
1156               else
1157                 {
1158                   lex_error (lexer, NULL);
1159                   goto error;
1160                 }
1161             }
1162         }
1163       else
1164         {
1165           lex_error (lexer, NULL);
1166           goto error;
1167         }
1168     }
1169
1170   if ( factor.rotation == ROT_NONE )
1171     factor.print &= ~PRINT_ROTATION;
1172
1173   if ( ! run_factor (ds, &factor)) 
1174     goto error;
1175
1176   free (factor.vars);
1177   return CMD_SUCCESS;
1178
1179  error:
1180   free (factor.vars);
1181   return CMD_FAILURE;
1182 }
1183
1184 static void do_factor (const struct cmd_factor *factor, struct casereader *group);
1185
1186
1187 static bool
1188 run_factor (struct dataset *ds, const struct cmd_factor *factor)
1189 {
1190   struct dictionary *dict = dataset_dict (ds);
1191   bool ok;
1192   struct casereader *group;
1193
1194   struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
1195
1196   while (casegrouper_get_next_group (grouper, &group))
1197     {
1198       if ( factor->missing_type == MISS_LISTWISE )
1199         group  = casereader_create_filter_missing (group, factor->vars, factor->n_vars,
1200                                                    factor->exclude,
1201                                                    NULL,  NULL);
1202       do_factor (factor, group);
1203     }
1204
1205   ok = casegrouper_destroy (grouper);
1206   ok = proc_commit (ds) && ok;
1207
1208   return ok;
1209 }
1210
1211
1212 /* Return the communality of variable N, calculated to N_FACTORS */
1213 static double
1214 the_communality (const gsl_matrix *evec, const gsl_vector *eval, int n, int n_factors)
1215 {
1216   size_t i;
1217
1218   double comm = 0;
1219
1220   assert (n >= 0);
1221   assert (n < eval->size);
1222   assert (n < evec->size1);
1223   assert (n_factors <= eval->size);
1224
1225   for (i = 0 ; i < n_factors; ++i)
1226     {
1227       double evali = fabs (gsl_vector_get (eval, i));
1228
1229       double eveci = gsl_matrix_get (evec, n, i);
1230
1231       comm += pow2 (eveci) * evali;
1232     }
1233
1234   return comm;
1235 }
1236
1237 /* Return the communality of variable N, calculated to N_FACTORS */
1238 static double
1239 communality (struct idata *idata, int n, int n_factors)
1240 {
1241   return the_communality (idata->evec, idata->eval, n, n_factors);
1242 }
1243
1244
1245 static void
1246 show_scree (const struct cmd_factor *f, struct idata *idata)
1247 {
1248   struct scree *s;
1249   const char *label ;
1250
1251   if ( !(f->plot & PLOT_SCREE) )
1252     return;
1253
1254
1255   label = f->extraction == EXTRACTION_PC ? _("Component Number") : _("Factor Number");
1256
1257   s = scree_create (idata->eval, label);
1258
1259   scree_submit (s);
1260 }
1261
1262 static void
1263 show_communalities (const struct cmd_factor * factor,
1264                     const gsl_vector *initial, const gsl_vector *extracted)
1265 {
1266   int i;
1267   int c = 0;
1268   const int heading_columns = 1;
1269   int nc = heading_columns;
1270   const int heading_rows = 1;
1271   const int nr = heading_rows + factor->n_vars;
1272   struct tab_table *t;
1273
1274   if (factor->print & PRINT_EXTRACTION)
1275     nc++;
1276
1277   if (factor->print & PRINT_INITIAL)
1278     nc++;
1279
1280   /* No point having a table with only headings */
1281   if (nc <= 1)
1282     return;
1283
1284   t = tab_create (nc, nr);
1285
1286   tab_title (t, _("Communalities"));
1287
1288   tab_headers (t, heading_columns, 0, heading_rows, 0);
1289
1290   c = 1;
1291   if (factor->print & PRINT_INITIAL)
1292     tab_text (t, c++, 0, TAB_CENTER | TAT_TITLE, _("Initial"));
1293
1294   if (factor->print & PRINT_EXTRACTION)
1295     tab_text (t, c++, 0, TAB_CENTER | TAT_TITLE, _("Extraction"));
1296
1297   /* Outline the box */
1298   tab_box (t,
1299            TAL_2, TAL_2,
1300            -1, -1,
1301            0, 0,
1302            nc - 1, nr - 1);
1303
1304   /* Vertical lines */
1305   tab_box (t,
1306            -1, -1,
1307            -1, TAL_1,
1308            heading_columns, 0,
1309            nc - 1, nr - 1);
1310
1311   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1312   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1313
1314   for (i = 0 ; i < factor->n_vars; ++i)
1315     {
1316       c = 0;
1317       tab_text (t, c++, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[i]));
1318
1319       if (factor->print & PRINT_INITIAL)
1320         tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL);
1321
1322       if (factor->print & PRINT_EXTRACTION)
1323         tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL);
1324     }
1325
1326   tab_submit (t);
1327 }
1328
1329
1330 static void
1331 show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const char *title, const gsl_matrix *fm)
1332 {
1333   int i;
1334   const int n_factors = idata->n_extractions;
1335
1336   const int heading_columns = 1;
1337   const int heading_rows = 2;
1338   const int nr = heading_rows + factor->n_vars;
1339   const int nc = heading_columns + n_factors;
1340   gsl_permutation *perm;
1341
1342   struct tab_table *t = tab_create (nc, nr);
1343
1344   /* 
1345   if ( factor->extraction == EXTRACTION_PC )
1346     tab_title (t, _("Component Matrix"));
1347   else 
1348     tab_title (t, _("Factor Matrix"));
1349   */
1350
1351   tab_title (t, "%s", title);
1352
1353   tab_headers (t, heading_columns, 0, heading_rows, 0);
1354
1355   if ( factor->extraction == EXTRACTION_PC )
1356     tab_joint_text (t,
1357                     1, 0,
1358                     nc - 1, 0,
1359                     TAB_CENTER | TAT_TITLE, _("Component"));
1360   else
1361     tab_joint_text (t,
1362                     1, 0,
1363                     nc - 1, 0,
1364                     TAB_CENTER | TAT_TITLE, _("Factor"));
1365
1366
1367   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1368
1369
1370   /* Outline the box */
1371   tab_box (t,
1372            TAL_2, TAL_2,
1373            -1, -1,
1374            0, 0,
1375            nc - 1, nr - 1);
1376
1377   /* Vertical lines */
1378   tab_box (t,
1379            -1, -1,
1380            -1, TAL_1,
1381            heading_columns, 1,
1382            nc - 1, nr - 1);
1383
1384   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1385   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1386
1387
1388   /* Initialise to the identity permutation */
1389   perm = gsl_permutation_calloc (factor->n_vars);
1390
1391   if ( factor->sort)
1392     sort_matrix_indirect (fm, perm);
1393
1394   for (i = 0 ; i < n_factors; ++i)
1395     {
1396       tab_text_format (t, heading_columns + i, 1, TAB_CENTER | TAT_TITLE, _("%d"), i + 1);
1397     }
1398
1399   for (i = 0 ; i < factor->n_vars; ++i)
1400     {
1401       int j;
1402       const int matrix_row = perm->data[i];
1403       tab_text (t, 0, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[matrix_row]));
1404
1405       for (j = 0 ; j < n_factors; ++j)
1406         {
1407           double x = gsl_matrix_get (fm, matrix_row, j);
1408
1409           if ( fabs (x) < factor->blank)
1410             continue;
1411
1412           tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL);
1413         }
1414     }
1415
1416   gsl_permutation_free (perm);
1417
1418   tab_submit (t);
1419 }
1420
1421
1422 static void
1423 show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
1424                          const gsl_vector *initial_eigenvalues,
1425                          const gsl_vector *extracted_eigenvalues,
1426                          const gsl_vector *rotated_loadings)
1427 {
1428   size_t i;
1429   int c = 0;
1430   const int heading_columns = 1;
1431   const int heading_rows = 2;
1432   const int nr = heading_rows + factor->n_vars;
1433
1434   struct tab_table *t ;
1435
1436   double i_total = 0.0;
1437   double i_cum = 0.0;
1438
1439   double e_total = 0.0;
1440   double e_cum = 0.0;
1441
1442   double r_cum = 0.0;
1443
1444   int nc = heading_columns;
1445
1446   if (factor->print & PRINT_EXTRACTION)
1447     nc += 3;
1448
1449   if (factor->print & PRINT_INITIAL)
1450     nc += 3;
1451
1452   if (factor->print & PRINT_ROTATION)
1453     nc += 3;
1454
1455   /* No point having a table with only headings */
1456   if ( nc <= heading_columns)
1457     return;
1458
1459   t = tab_create (nc, nr);
1460
1461   tab_title (t, _("Total Variance Explained"));
1462
1463   tab_headers (t, heading_columns, 0, heading_rows, 0);
1464
1465   /* Outline the box */
1466   tab_box (t,
1467            TAL_2, TAL_2,
1468            -1, -1,
1469            0, 0,
1470            nc - 1, nr - 1);
1471
1472   /* Vertical lines */
1473   tab_box (t,
1474            -1, -1,
1475            -1, TAL_1,
1476            heading_columns, 0,
1477            nc - 1, nr - 1);
1478
1479   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1480   tab_hline (t, TAL_1, 1, nc - 1, 1);
1481
1482   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1483
1484
1485   if ( factor->extraction == EXTRACTION_PC)
1486     tab_text (t, 0, 1, TAB_LEFT | TAT_TITLE, _("Component"));
1487   else
1488     tab_text (t, 0, 1, TAB_LEFT | TAT_TITLE, _("Factor"));
1489
1490   c = 1;
1491   if (factor->print & PRINT_INITIAL)
1492     {
1493       tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Initial Eigenvalues"));
1494       c += 3;
1495     }
1496
1497   if (factor->print & PRINT_EXTRACTION)
1498     {
1499       tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Extraction Sums of Squared Loadings"));
1500       c += 3;
1501     }
1502
1503   if (factor->print & PRINT_ROTATION)
1504     {
1505       tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Rotation Sums of Squared Loadings"));
1506       c += 3;
1507     }
1508
1509   for (i = 0; i < (nc - heading_columns) / 3 ; ++i)
1510     {
1511       tab_text (t, i * 3 + 1, 1, TAB_CENTER | TAT_TITLE, _("Total"));
1512       /* xgettext:no-c-format */
1513       tab_text (t, i * 3 + 2, 1, TAB_CENTER | TAT_TITLE, _("% of Variance"));
1514       tab_text (t, i * 3 + 3, 1, TAB_CENTER | TAT_TITLE, _("Cumulative %"));
1515
1516       tab_vline (t, TAL_2, heading_columns + i * 3, 0, nr - 1);
1517     }
1518
1519   for (i = 0 ; i < initial_eigenvalues->size; ++i)
1520     i_total += gsl_vector_get (initial_eigenvalues, i);
1521
1522   if ( factor->extraction == EXTRACTION_PAF)
1523     {
1524       e_total = factor->n_vars;
1525     }
1526   else
1527     {
1528       e_total = i_total;
1529     }
1530
1531   for (i = 0 ; i < factor->n_vars; ++i)
1532     {
1533       const double i_lambda = gsl_vector_get (initial_eigenvalues, i);
1534       double i_percent = 100.0 * i_lambda / i_total ;
1535
1536       const double e_lambda = gsl_vector_get (extracted_eigenvalues, i);
1537       double e_percent = 100.0 * e_lambda / e_total ;
1538
1539       const double r_lambda = gsl_vector_get (rotated_loadings, i);
1540       double r_percent = 100.0 * r_lambda / e_total ;
1541
1542       c = 0;
1543
1544       tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%zu"), i + 1);
1545
1546       i_cum += i_percent;
1547       e_cum += e_percent;
1548       r_cum += r_percent;
1549
1550       /* Initial Eigenvalues */
1551       if (factor->print & PRINT_INITIAL)
1552       {
1553         tab_double (t, c++, i + heading_rows, 0, i_lambda, NULL);
1554         tab_double (t, c++, i + heading_rows, 0, i_percent, NULL);
1555         tab_double (t, c++, i + heading_rows, 0, i_cum, NULL);
1556       }
1557
1558
1559       if (factor->print & PRINT_EXTRACTION)
1560         {
1561           if (i < idata->n_extractions)
1562             {
1563               /* Sums of squared loadings */
1564               tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL);
1565               tab_double (t, c++, i + heading_rows, 0, e_percent, NULL);
1566               tab_double (t, c++, i + heading_rows, 0, e_cum, NULL);
1567             }
1568         }
1569
1570       if (factor->print & PRINT_ROTATION)
1571         {
1572           if (i < idata->n_extractions)
1573             {
1574               tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL);
1575               tab_double (t, c++, i + heading_rows, 0, r_percent, NULL);
1576               tab_double (t, c++, i + heading_rows, 0, r_cum, NULL);
1577             }
1578       }
1579
1580     }
1581
1582   tab_submit (t);
1583 }
1584
1585
1586 static void
1587 show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata)
1588 {
1589   struct tab_table *t ;
1590   size_t i, j;
1591   int y_pos_corr = -1;
1592   int y_pos_sig = -1;
1593   int suffix_rows = 0;
1594
1595   const int heading_rows = 1;
1596   const int heading_columns = 2;
1597
1598   int nc = heading_columns ;
1599   int nr = heading_rows ;
1600   int n_data_sets = 0;
1601
1602   if (factor->print & PRINT_CORRELATION)
1603     {
1604       y_pos_corr = n_data_sets;
1605       n_data_sets++;
1606       nc = heading_columns + factor->n_vars;
1607     }
1608
1609   if (factor->print & PRINT_SIG)
1610     {
1611       y_pos_sig = n_data_sets;
1612       n_data_sets++;
1613       nc = heading_columns + factor->n_vars;
1614     }
1615
1616   nr += n_data_sets * factor->n_vars;
1617
1618   if (factor->print & PRINT_DETERMINANT)
1619     suffix_rows = 1;
1620
1621   /* If the table would contain only headings, don't bother rendering it */
1622   if (nr <= heading_rows && suffix_rows == 0)
1623     return;
1624
1625   t = tab_create (nc, nr + suffix_rows);
1626
1627   tab_title (t, _("Correlation Matrix"));
1628
1629   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1630
1631   if (nr > heading_rows)
1632     {
1633       tab_headers (t, heading_columns, 0, heading_rows, 0);
1634
1635       tab_vline (t, TAL_2, 2, 0, nr - 1);
1636
1637       /* Outline the box */
1638       tab_box (t,
1639                TAL_2, TAL_2,
1640                -1, -1,
1641                0, 0,
1642                nc - 1, nr - 1);
1643
1644       /* Vertical lines */
1645       tab_box (t,
1646                -1, -1,
1647                -1, TAL_1,
1648                heading_columns, 0,
1649                nc - 1, nr - 1);
1650
1651
1652       for (i = 0; i < factor->n_vars; ++i)
1653         tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i]));
1654
1655
1656       for (i = 0 ; i < n_data_sets; ++i)
1657         {
1658           int y = heading_rows + i * factor->n_vars;
1659           size_t v;
1660           for (v = 0; v < factor->n_vars; ++v)
1661             tab_text (t, 1, y + v, TAT_TITLE, var_to_string (factor->vars[v]));
1662
1663           tab_hline (t, TAL_1, 0, nc - 1, y);
1664         }
1665
1666       if (factor->print & PRINT_CORRELATION)
1667         {
1668           const double y = heading_rows + y_pos_corr;
1669           tab_text (t, 0, y, TAT_TITLE, _("Correlations"));
1670
1671           for (i = 0; i < factor->n_vars; ++i)
1672             {
1673               for (j = 0; j < factor->n_vars; ++j)
1674                 tab_double (t, heading_columns + i,  y + j, 0, gsl_matrix_get (idata->corr, i, j), NULL);
1675             }
1676         }
1677
1678       if (factor->print & PRINT_SIG)
1679         {
1680           const double y = heading_rows + y_pos_sig * factor->n_vars;
1681           tab_text (t, 0, y, TAT_TITLE, _("Sig. (1-tailed)"));
1682
1683           for (i = 0; i < factor->n_vars; ++i)
1684             {
1685               for (j = 0; j < factor->n_vars; ++j)
1686                 {
1687                   double rho = gsl_matrix_get (idata->corr, i, j);
1688                   double w = gsl_matrix_get (idata->n, i, j);
1689
1690                   if (i == j)
1691                     continue;
1692
1693                   tab_double (t, heading_columns + i,  y + j, 0, significance_of_correlation (rho, w), NULL);
1694                 }
1695             }
1696         }
1697     }
1698
1699   if (factor->print & PRINT_DETERMINANT)
1700     {
1701       tab_text (t, 0, nr, TAB_LEFT | TAT_TITLE, _("Determinant"));
1702
1703       tab_double (t, 1, nr, 0, idata->detR, NULL);
1704     }
1705
1706   tab_submit (t);
1707 }
1708
1709
1710
1711 static void
1712 do_factor (const struct cmd_factor *factor, struct casereader *r)
1713 {
1714   struct ccase *c;
1715   const gsl_matrix *var_matrix;
1716   const gsl_matrix *mean_matrix;
1717
1718   const gsl_matrix *analysis_matrix;
1719   struct idata *idata = idata_alloc (factor->n_vars);
1720
1721   struct covariance *cov = covariance_1pass_create (factor->n_vars, factor->vars,
1722                                               factor->wv, factor->exclude);
1723
1724   for ( ; (c = casereader_read (r) ); case_unref (c))
1725     {
1726       covariance_accumulate (cov, c);
1727     }
1728
1729   idata->cov = covariance_calculate (cov);
1730
1731   if (idata->cov == NULL)
1732     {
1733       msg (MW, _("The dataset contains no complete observations. No analysis will be performed."));
1734       goto finish;
1735     }
1736
1737   var_matrix = covariance_moments (cov, MOMENT_VARIANCE);
1738   mean_matrix = covariance_moments (cov, MOMENT_MEAN);
1739   idata->n = covariance_moments (cov, MOMENT_NONE);
1740
1741   if ( factor->method == METHOD_CORR)
1742     {
1743       idata->corr = correlation_from_covariance (idata->cov, var_matrix);
1744       
1745       analysis_matrix = idata->corr;
1746     }
1747   else
1748     analysis_matrix = idata->cov;
1749
1750   if (factor->print & PRINT_DETERMINANT
1751       || factor->print & PRINT_KMO)
1752     {
1753       int sign = 0;
1754
1755       const int size = idata->corr->size1;
1756       gsl_permutation *p = gsl_permutation_calloc (size);
1757       gsl_matrix *tmp = gsl_matrix_calloc (size, size);
1758       gsl_matrix_memcpy (tmp, idata->corr);
1759
1760       gsl_linalg_LU_decomp (tmp, p, &sign);
1761       idata->detR = gsl_linalg_LU_det (tmp, sign);
1762       gsl_permutation_free (p);
1763       gsl_matrix_free (tmp);
1764     }
1765
1766   if ( factor->print & PRINT_UNIVARIATE)
1767     {
1768       const struct fmt_spec *wfmt = factor->wv ? var_get_print_format (factor->wv) : & F_8_0;
1769       const int nc = 4;
1770       int i;
1771
1772       const int heading_columns = 1;
1773       const int heading_rows = 1;
1774
1775       const int nr = heading_rows + factor->n_vars;
1776
1777       struct tab_table *t = tab_create (nc, nr);
1778       tab_title (t, _("Descriptive Statistics"));
1779
1780       tab_headers (t, heading_columns, 0, heading_rows, 0);
1781
1782       /* Outline the box */
1783       tab_box (t,
1784                TAL_2, TAL_2,
1785                -1, -1,
1786                0, 0,
1787                nc - 1, nr - 1);
1788
1789       /* Vertical lines */
1790       tab_box (t,
1791                -1, -1,
1792                -1, TAL_1,
1793                heading_columns, 0,
1794                nc - 1, nr - 1);
1795
1796       tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1797       tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1798
1799       tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
1800       tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
1801       tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Analysis N"));
1802
1803       for (i = 0 ; i < factor->n_vars; ++i)
1804         {
1805           const struct variable *v = factor->vars[i];
1806           tab_text (t, 0, i + heading_rows, TAB_LEFT | TAT_TITLE, var_to_string (v));
1807
1808           tab_double (t, 1, i + heading_rows, 0, gsl_matrix_get (mean_matrix, i, i), NULL);
1809           tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (var_matrix, i, i)), NULL);
1810           tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->n, i, i), wfmt);
1811         }
1812
1813       tab_submit (t);
1814     }
1815
1816   if (factor->print & PRINT_KMO)
1817     {
1818       int i;
1819       double sum_ssq_r = 0;
1820       double sum_ssq_a = 0;
1821
1822       double df = factor->n_vars * ( factor->n_vars - 1) / 2;
1823
1824       double w = 0;
1825
1826
1827       double xsq;
1828
1829       const int heading_columns = 2;
1830       const int heading_rows = 0;
1831
1832       const int nr = heading_rows + 4;
1833       const int nc = heading_columns + 1;
1834
1835       gsl_matrix *a, *x;
1836
1837       struct tab_table *t = tab_create (nc, nr);
1838       tab_title (t, _("KMO and Bartlett's Test"));
1839
1840       x  = clone_matrix (idata->corr);
1841       gsl_linalg_cholesky_decomp (x);
1842       gsl_linalg_cholesky_invert (x);
1843
1844       a = anti_image (x);
1845
1846       for (i = 0; i < x->size1; ++i)
1847         {
1848           sum_ssq_r += ssq_od_n (x, i);
1849           sum_ssq_a += ssq_od_n (a, i);
1850         }
1851
1852       gsl_matrix_free (a);
1853       gsl_matrix_free (x);
1854
1855       tab_headers (t, heading_columns, 0, heading_rows, 0);
1856
1857       /* Outline the box */
1858       tab_box (t,
1859                TAL_2, TAL_2,
1860                -1, -1,
1861                0, 0,
1862                nc - 1, nr - 1);
1863
1864       tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1865
1866       tab_text (t, 0, 0, TAT_TITLE | TAB_LEFT, _("Kaiser-Meyer-Olkin Measure of Sampling Adequacy"));
1867
1868       tab_double (t, 2, 0, 0, sum_ssq_r /  (sum_ssq_r + sum_ssq_a), NULL);
1869
1870       tab_text (t, 0, 1, TAT_TITLE | TAB_LEFT, _("Bartlett's Test of Sphericity"));
1871
1872       tab_text (t, 1, 1, TAT_TITLE, _("Approx. Chi-Square"));
1873       tab_text (t, 1, 2, TAT_TITLE, _("df"));
1874       tab_text (t, 1, 3, TAT_TITLE, _("Sig."));
1875
1876
1877       /* The literature doesn't say what to do for the value of W when 
1878          missing values are involved.  The best thing I can think of
1879          is to take the mean average. */
1880       w = 0;
1881       for (i = 0; i < idata->n->size1; ++i)
1882         w += gsl_matrix_get (idata->n, i, i);
1883       w /= idata->n->size1;
1884
1885       xsq = w - 1 - (2 * factor->n_vars + 5) / 6.0;
1886       xsq *= -log (idata->detR);
1887
1888       tab_double (t, 2, 1, 0, xsq, NULL);
1889       tab_double (t, 2, 2, 0, df, &F_8_0);
1890       tab_double (t, 2, 3, 0, gsl_cdf_chisq_Q (xsq, df), NULL);
1891       
1892
1893       tab_submit (t);
1894     }
1895
1896   show_correlation_matrix (factor, idata);
1897
1898   {
1899     gsl_eigen_symmv_workspace *workspace = gsl_eigen_symmv_alloc (factor->n_vars);
1900     
1901     gsl_eigen_symmv (matrix_dup (analysis_matrix), idata->eval, idata->evec, workspace);
1902
1903     gsl_eigen_symmv_free (workspace);
1904   }
1905
1906   gsl_eigen_symmv_sort (idata->eval, idata->evec, GSL_EIGEN_SORT_ABS_DESC);
1907
1908   idata->n_extractions = n_extracted_factors (factor, idata);
1909
1910   if (idata->n_extractions == 0)
1911     {
1912       msg (MW, _("The FACTOR criteria result in zero factors extracted. Therefore no analysis will be performed."));
1913       goto finish;
1914     }
1915
1916   if (idata->n_extractions > factor->n_vars)
1917     {
1918       msg (MW, _("The FACTOR criteria result in more factors than variables, which is not meaningful. No analysis will be performed."));
1919       goto finish;
1920     }
1921     
1922   {
1923     gsl_matrix *rotated_factors = NULL;
1924     gsl_vector *rotated_loadings = NULL;
1925
1926     const gsl_vector *extracted_eigenvalues = NULL;
1927     gsl_vector *initial_communalities = gsl_vector_alloc (factor->n_vars);
1928     gsl_vector *extracted_communalities = gsl_vector_alloc (factor->n_vars);
1929     size_t i;
1930     struct factor_matrix_workspace *fmw = factor_matrix_workspace_alloc (idata->msr->size, idata->n_extractions);
1931     gsl_matrix *factor_matrix = gsl_matrix_calloc (factor->n_vars, fmw->n_factors);
1932
1933     if ( factor->extraction == EXTRACTION_PAF)
1934       {
1935         gsl_vector *diff = gsl_vector_alloc (idata->msr->size);
1936         struct smr_workspace *ws = ws_create (analysis_matrix);
1937
1938         for (i = 0 ; i < factor->n_vars ; ++i)
1939           {
1940             double r2 = squared_multiple_correlation (analysis_matrix, i, ws);
1941
1942             gsl_vector_set (idata->msr, i, r2);
1943           }
1944         ws_destroy (ws);
1945
1946         gsl_vector_memcpy (initial_communalities, idata->msr);
1947
1948         for (i = 0; i < factor->iterations; ++i)
1949           {
1950             double min, max;
1951             gsl_vector_memcpy (diff, idata->msr);
1952
1953             iterate_factor_matrix (analysis_matrix, idata->msr, factor_matrix, fmw);
1954       
1955             gsl_vector_sub (diff, idata->msr);
1956
1957             gsl_vector_minmax (diff, &min, &max);
1958       
1959             if ( fabs (min) < factor->econverge && fabs (max) < factor->econverge)
1960               break;
1961           }
1962         gsl_vector_free (diff);
1963
1964
1965
1966         gsl_vector_memcpy (extracted_communalities, idata->msr);
1967         extracted_eigenvalues = fmw->eval;
1968       }
1969     else if (factor->extraction == EXTRACTION_PC)
1970       {
1971         for (i = 0; i < factor->n_vars; ++i)
1972           gsl_vector_set (initial_communalities, i, communality (idata, i, factor->n_vars));
1973
1974         gsl_vector_memcpy (extracted_communalities, initial_communalities);
1975
1976         iterate_factor_matrix (analysis_matrix, extracted_communalities, factor_matrix, fmw);
1977
1978
1979         extracted_eigenvalues = idata->eval;
1980       }
1981
1982
1983     show_communalities (factor, initial_communalities, extracted_communalities);
1984
1985
1986     if ( factor->rotation != ROT_NONE)
1987       {
1988         rotated_factors = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2);
1989         rotated_loadings = gsl_vector_calloc (factor_matrix->size2);
1990
1991         rotate (factor, factor_matrix, extracted_communalities, rotated_factors, rotated_loadings);
1992       }
1993
1994     show_explained_variance (factor, idata, idata->eval, extracted_eigenvalues, rotated_loadings);
1995
1996     factor_matrix_workspace_free (fmw);
1997
1998     show_scree (factor, idata);
1999
2000     show_factor_matrix (factor, idata,
2001                         factor->extraction == EXTRACTION_PC ? _("Component Matrix") : _("Factor Matrix"),
2002                         factor_matrix);
2003
2004     if ( factor->rotation != ROT_NONE)
2005       {
2006         show_factor_matrix (factor, idata,
2007                             factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") : _("Rotated Factor Matrix"),
2008                             rotated_factors);
2009
2010         gsl_matrix_free (rotated_factors);
2011       }
2012
2013
2014
2015     gsl_vector_free (initial_communalities);
2016     gsl_vector_free (extracted_communalities);
2017   }
2018
2019  finish:
2020
2021   idata_free (idata);
2022
2023   casereader_destroy (r);
2024 }
2025
2026
2027