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