pintos-os.org Git - pspp/blob - 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