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/message.h"
  39 #include "libpspp/misc.h"
  40 #include "math/correlation.h"
  41 #include "math/covariance.h"
  42 #include "math/moments.h"
  43 #include "output/chart-item.h"
  44 #include "output/charts/scree.h"
  45 #include "output/tab.h"
  46
  47 #include "gettext.h"
  48 #define _(msgid) gettext (msgid)
  49 #define N_(msgid) msgid
  50
  51 enum method
  52   {
  53     METHOD_CORR,
  54     METHOD_COV
  55   };
  56
  57 enum missing_type
  58   {
  59     MISS_LISTWISE,
  60     MISS_PAIRWISE,
  61     MISS_MEANSUB,
  62   };
  63
  64 enum extraction_method
  65   {
  66     EXTRACTION_PC,
  67     EXTRACTION_PAF,
  68   };
  69
  70 enum plot_opts
  71   {
  72     PLOT_SCREE = 0x0001,
  73     PLOT_ROTATION = 0x0002
  74   };
  75
  76 enum print_opts
  77   {
  78     PRINT_UNIVARIATE  = 0x0001,
  79     PRINT_DETERMINANT = 0x0002,
  80     PRINT_INV         = 0x0004,
  81     PRINT_AIC         = 0x0008,
  82     PRINT_SIG         = 0x0010,
  83     PRINT_COVARIANCE  = 0x0020,
  84     PRINT_CORRELATION = 0x0040,
  85     PRINT_ROTATION    = 0x0080,
  86     PRINT_EXTRACTION  = 0x0100,
  87     PRINT_INITIAL     = 0x0200,
  88     PRINT_KMO         = 0x0400,
  89     PRINT_REPR        = 0x0800,
  90     PRINT_FSCORE      = 0x1000
  91   };
  92
  93 enum rotation_type
  94   {
  95     ROT_VARIMAX = 0,
  96     ROT_EQUAMAX,
  97     ROT_QUARTIMAX,
  98     ROT_NONE
  99   };
 100
 101 typedef void (*rotation_coefficients) (double *x, double *y,
 102                                     double a, double b, double c, double d,
 103                                     const gsl_matrix *loadings );
 104
 105
 106 static void
 107 varimax_coefficients (double *x, double *y,
 108                       double a, double b, double c, double d,
 109                       const gsl_matrix *loadings )
 110 {
 111   *x = d - 2 * a * b / loadings->size1;
 112   *y = c - (a * a - b * b) / loadings->size1;
 113 }
 114
 115 static void
 116 equamax_coefficients (double *x, double *y,
 117                       double a, double b, double c, double d,
 118                       const gsl_matrix *loadings )
 119 {
 120   *x = d - loadings->size2 * a * b / loadings->size1;
 121   *y = c - loadings->size2 * (a * a - b * b) / (2 * loadings->size1);
 122 }
 123
 124 static void
 125 quartimax_coefficients (double *x, double *y,
 126                       double a UNUSED, double b UNUSED, double c, double d,
 127                       const gsl_matrix *loadings UNUSED)
 128 {
 129   *x = d ;
 130   *y = c ;
 131 }
 132
 133 static const rotation_coefficients rotation_coeff[3] = {
 134   varimax_coefficients,
 135   equamax_coefficients,
 136   quartimax_coefficients
 137 };
 138
 139
 140 struct cmd_factor
 141 {
 142   size_t n_vars;
 143   const struct variable **vars;
 144
 145   const struct variable *wv;
 146
 147   enum method method;
 148   enum missing_type missing_type;
 149   enum mv_class exclude;
 150   enum print_opts print;
 151   enum extraction_method extraction;
 152   enum plot_opts plot;
 153   enum rotation_type rotation;
 154
 155   /* Extraction Criteria */
 156   int n_factors;
 157   double min_eigen;
 158   double econverge;
 159   int iterations;
 160
 161   double rconverge;
 162
 163   /* Format */
 164   double blank;
 165   bool sort;
 166 };
 167
 168 struct idata
 169 {
 170   /* Intermediate values used in calculation */
 171
 172   const gsl_matrix *corr ;  /* The correlation matrix */
 173   gsl_matrix *cov ;         /* The covariance matrix */
 174   const gsl_matrix *n ;     /* Matrix of number of samples */
 175
 176   gsl_vector *eval ;  /* The eigenvalues */
 177   gsl_matrix *evec ;  /* The eigenvectors */
 178
 179   int n_extractions;
 180
 181   gsl_vector *msr ;  /* Multiple Squared Regressions */
 182
 183   double detR;  /* The determinant of the correlation matrix */
 184 };
 185
 186 static struct idata *
 187 idata_alloc (size_t n_vars)
 188 {
 189   struct idata *id = xzalloc (sizeof (*id));
 190
 191   id->n_extractions = 0;
 192   id->msr = gsl_vector_alloc (n_vars);
 193
 194   id->eval = gsl_vector_alloc (n_vars);
 195   id->evec = gsl_matrix_alloc (n_vars, n_vars);
 196
 197   return id;
 198 }
 199
 200 static void
 201 idata_free (struct idata *id)
 202 {
 203   gsl_vector_free (id->msr);
 204   gsl_vector_free (id->eval);
 205   gsl_matrix_free (id->evec);
 206   if (id->cov != NULL)
 207     gsl_matrix_free (id->cov);
 208   if (id->corr != NULL)
 209     gsl_matrix_free (id->corr);
 210
 211   free (id);
 212 }
 213
 214
 215 static gsl_matrix *
 216 anti_image (const gsl_matrix *m)
 217 {
 218   int i, j;
 219   gsl_matrix *a;
 220   assert (m->size1 == m->size2);
 221
 222   a = gsl_matrix_alloc (m->size1, m->size2);
 223
 224   for (i = 0; i < m->size1; ++i)
 225     {
 226       for (j = 0; j < m->size2; ++j)
 227         {
 228           double *p = gsl_matrix_ptr (a, i, j);
 229           *p = gsl_matrix_get (m, i, j);
 230           *p /= gsl_matrix_get (m, i, i);
 231           *p /= gsl_matrix_get (m, j, j);
 232         }
 233     }
 234
 235   return a;
 236 }
 237
 238
 239 /* Return the sum of all the elements excluding row N */
 240 static double
 241 ssq_od_n (const gsl_matrix *m, int n)
 242 {
 243   int i, j;
 244   double ss = 0;
 245   assert (m->size1 == m->size2);
 246
 247   assert (n < m->size1);
 248
 249   for (i = 0; i < m->size1; ++i)
 250     {
 251       if (i == n ) continue;
 252       for (j = 0; j < m->size2; ++j)
 253         {
 254           ss += pow2 (gsl_matrix_get (m, i, j));
 255         }
 256     }
 257
 258   return ss;
 259 }
 260
 261
 262
 263 #if 0
 264 static void
 265 dump_matrix (const gsl_matrix *m)
 266 {
 267   size_t i, j;
 268
 269   for (i = 0 ; i < m->size1; ++i)
 270     {
 271       for (j = 0 ; j < m->size2; ++j)
 272         printf ("%02f ", gsl_matrix_get (m, i, j));
 273       printf ("\n");
 274     }
 275 }
 276
 277 static void
 278 dump_matrix_permute (const gsl_matrix *m, const gsl_permutation *p)
 279 {
 280   size_t i, j;
 281
 282   for (i = 0 ; i < m->size1; ++i)
 283     {
 284       for (j = 0 ; j < m->size2; ++j)
 285         printf ("%02f ", gsl_matrix_get (m, gsl_permutation_get (p, i), j));
 286       printf ("\n");
 287     }
 288 }
 289
 290
 291 static void
 292 dump_vector (const gsl_vector *v)
 293 {
 294   size_t i;
 295   for (i = 0 ; i < v->size; ++i)
 296     {
 297       printf ("%02f\n", gsl_vector_get (v, i));
 298     }
 299   printf ("\n");
 300 }
 301 #endif
 302
 303
 304 static int
 305 n_extracted_factors (const struct cmd_factor *factor, struct idata *idata)
 306 {
 307   int i;
 308
 309   /* If there is a cached value, then return that. */
 310   if ( idata->n_extractions != 0)
 311     return idata->n_extractions;
 312
 313   /* Otherwise, if the number of factors has been explicitly requested,
 314      use that. */
 315   if (factor->n_factors > 0)
 316     {
 317       idata->n_extractions = factor->n_factors;
 318       goto finish;
 319     }
 320
 321   /* Use the MIN_EIGEN setting. */
 322   for (i = 0 ; i < idata->eval->size; ++i)
 323     {
 324       double evali = fabs (gsl_vector_get (idata->eval, i));
 325
 326       idata->n_extractions = i;
 327
 328       if (evali < factor->min_eigen)
 329         goto finish;
 330     }
 331
 332  finish:
 333   return idata->n_extractions;
 334 }
 335
 336
 337 /* Returns a newly allocated matrix identical to M.
 338    It it the callers responsibility to free the returned value.
 339 */
 340 static gsl_matrix *
 341 matrix_dup (const gsl_matrix *m)
 342 {
 343   gsl_matrix *n =  gsl_matrix_alloc (m->size1, m->size2);
 344
 345   gsl_matrix_memcpy (n, m);
 346
 347   return n;
 348 }
 349
 350
 351 struct smr_workspace
 352 {
 353   /* Copy of the subject */
 354   gsl_matrix *m;
 355
 356   gsl_matrix *inverse;
 357
 358   gsl_permutation *perm;
 359
 360   gsl_matrix *result1;
 361   gsl_matrix *result2;
 362 };
 363
 364
 365 static struct smr_workspace *ws_create (const gsl_matrix *input)
 366 {
 367   struct smr_workspace *ws = xmalloc (sizeof (*ws));
 368
 369   ws->m = gsl_matrix_alloc (input->size1, input->size2);
 370   ws->inverse = gsl_matrix_calloc (input->size1 - 1, input->size2 - 1);
 371   ws->perm = gsl_permutation_alloc (input->size1 - 1);
 372   ws->result1 = gsl_matrix_calloc (input->size1 - 1, 1);
 373   ws->result2 = gsl_matrix_calloc (1, 1);
 374
 375   return ws;
 376 }
 377
 378 static void
 379 ws_destroy (struct smr_workspace *ws)
 380 {
 381   gsl_matrix_free (ws->result2);
 382   gsl_matrix_free (ws->result1);
 383   gsl_permutation_free (ws->perm);
 384   gsl_matrix_free (ws->inverse);
 385   gsl_matrix_free (ws->m);
 386
 387   free (ws);
 388 }
 389
 390
 391 /*
 392    Return the square of the regression coefficient for VAR regressed against all other variables.
 393  */
 394 static double
 395 squared_multiple_correlation (const gsl_matrix *corr, int var, struct smr_workspace *ws)
 396 {
 397   /* For an explanation of what this is doing, see
 398      http://www.visualstatistics.net/Visual%20Statistics%20Multimedia/multiple_regression_analysis.htm
 399   */
 400
 401   int signum = 0;
 402   gsl_matrix_view rxx;
 403
 404   gsl_matrix_memcpy (ws->m, corr);
 405
 406   gsl_matrix_swap_rows (ws->m, 0, var);
 407   gsl_matrix_swap_columns (ws->m, 0, var);
 408
 409   rxx = gsl_matrix_submatrix (ws->m, 1, 1, ws->m->size1 - 1, ws->m->size1 - 1);
 410
 411   gsl_linalg_LU_decomp (&rxx.matrix, ws->perm, &signum);
 412
 413   gsl_linalg_LU_invert (&rxx.matrix, ws->perm, ws->inverse);
 414
 415   {
 416     gsl_matrix_const_view rxy = gsl_matrix_const_submatrix (ws->m, 1, 0, ws->m->size1 - 1, 1);
 417     gsl_matrix_const_view ryx = gsl_matrix_const_submatrix (ws->m, 0, 1, 1, ws->m->size1 - 1);
 418
 419     gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans,
 420                     1.0, ws->inverse, &rxy.matrix, 0.0, ws->result1);
 421
 422     gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans,
 423                     1.0, &ryx.matrix, ws->result1, 0.0, ws->result2);
 424   }
 425
 426   return gsl_matrix_get (ws->result2, 0, 0);
 427 }
 428
 429
 430
 431 static double the_communality (const gsl_matrix *evec, const gsl_vector *eval, int n, int n_factors);
 432
 433
 434 struct factor_matrix_workspace
 435 {
 436   size_t n_factors;
 437   gsl_eigen_symmv_workspace *eigen_ws;
 438
 439   gsl_vector *eval ;
 440   gsl_matrix *evec ;
 441
 442   gsl_matrix *gamma ;
 443
 444   gsl_matrix *r;
 445 };
 446
 447 static struct factor_matrix_workspace *
 448 factor_matrix_workspace_alloc (size_t n, size_t nf)
 449 {
 450   struct factor_matrix_workspace *ws = xmalloc (sizeof (*ws));
 451
 452   ws->n_factors = nf;
 453   ws->gamma = gsl_matrix_calloc (nf, nf);
 454   ws->eigen_ws = gsl_eigen_symmv_alloc (n);
 455   ws->eval = gsl_vector_alloc (n);
 456   ws->evec = gsl_matrix_alloc (n, n);
 457   ws->r  = gsl_matrix_alloc (n, n);
 458
 459   return ws;
 460 }
 461
 462 static void
 463 factor_matrix_workspace_free (struct factor_matrix_workspace *ws)
 464 {
 465   gsl_eigen_symmv_free (ws->eigen_ws);
 466   gsl_vector_free (ws->eval);
 467   gsl_matrix_free (ws->evec);
 468   gsl_matrix_free (ws->gamma);
 469   gsl_matrix_free (ws->r);
 470   free (ws);
 471 }
 472
 473 /*
 474   Shift P left by OFFSET places, and overwrite TARGET
 475   with the shifted result.
 476   Positions in TARGET less than OFFSET are unchanged.
 477 */
 478 static void
 479 perm_shift_apply (gsl_permutation *target, const gsl_permutation *p,
 480                   size_t offset)
 481 {
 482   size_t i;
 483   assert (target->size == p->size);
 484   assert (offset <= target->size);
 485
 486   for (i = 0; i < target->size - offset; ++i)
 487     {
 488       target->data[i] = p->data [i + offset];
 489     }
 490 }
 491
 492
 493 /*
 494    Indirectly sort the rows of matrix INPUT, storing the sort order in PERM.
 495    The sort criteria are as follows:
 496
 497    Rows are sorted on the first column, until the absolute value of an
 498    element in a subsequent column  is greater than that of the first
 499    column.  Thereafter, rows will be sorted on the second column,
 500    until the absolute value of an element in a subsequent column
 501    exceeds that of the second column ...
 502 */
 503 static void
 504 sort_matrix_indirect (const gsl_matrix *input, gsl_permutation *perm)
 505 {
 506   const size_t n = perm->size;
 507   const size_t m = input->size2;
 508   int i, j;
 509   gsl_matrix *mat ;
 510   int column_n = 0;
 511   int row_n = 0;
 512   gsl_permutation *p;
 513
 514   assert (perm->size == input->size1);
 515
 516   p = gsl_permutation_alloc (n);
 517
 518   /* Copy INPUT into MAT, discarding the sign */
 519   mat = gsl_matrix_alloc (n, m);
 520   for (i = 0 ; i < mat->size1; ++i)
 521     {
 522       for (j = 0 ; j < mat->size2; ++j)
 523         {
 524           double x = gsl_matrix_get (input, i, j);
 525           gsl_matrix_set (mat, i, j, fabs (x));
 526         }
 527     }
 528
 529   while (column_n < m && row_n < n)
 530     {
 531       gsl_vector_const_view columni = gsl_matrix_const_column (mat, column_n);
 532       gsl_sort_vector_index (p, &columni.vector);
 533
 534       for (i = 0 ; i < n; ++i)
 535         {
 536           gsl_vector_view row = gsl_matrix_row (mat, p->data[n - 1 - i]);
 537           size_t maxindex = gsl_vector_max_index (&row.vector);
 538
 539           if ( maxindex > column_n )
 540             break;
 541
 542           /* All subsequent elements of this row, are of no interest.
 543              So set them all to a highly negative value */
 544           for (j = column_n + 1; j < row.vector.size ; ++j)
 545             gsl_vector_set (&row.vector, j, -DBL_MAX);
 546         }
 547
 548       perm_shift_apply (perm, p, row_n);
 549       row_n += i;
 550
 551       column_n++;
 552     }
 553
 554   gsl_permutation_free (p);
 555   gsl_matrix_free (mat);
 556
 557   assert ( 0 == gsl_permutation_valid (perm));
 558
 559   /* We want the biggest value to be first */
 560   gsl_permutation_reverse (perm);
 561 }
 562
 563
 564 static void
 565 drot_go (double phi, double *l0, double *l1)
 566 {
 567   double r0 = cos (phi) * *l0 + sin (phi) * *l1;
 568   double r1 = - sin (phi) * *l0 + cos (phi) * *l1;
 569
 570   *l0 = r0;
 571   *l1 = r1;
 572 }
 573
 574
 575 static gsl_matrix *
 576 clone_matrix (const gsl_matrix *m)
 577 {
 578   int j, k;
 579   gsl_matrix *c = gsl_matrix_calloc (m->size1, m->size2);
 580
 581   for (j = 0 ; j < c->size1; ++j)
 582     {
 583       for (k = 0 ; k < c->size2; ++k)
 584         {
 585           const double *v = gsl_matrix_const_ptr (m, j, k);
 586           gsl_matrix_set (c, j, k, *v);
 587         }
 588     }
 589
 590   return c;
 591 }
 592
 593
 594 static double
 595 initial_sv (const gsl_matrix *fm)
 596 {
 597   int j, k;
 598
 599   double sv = 0.0;
 600   for (j = 0 ; j < fm->size2; ++j)
 601     {
 602       double l4s = 0;
 603       double l2s = 0;
 604
 605       for (k = j + 1 ; k < fm->size2; ++k)
 606         {
 607           double lambda = gsl_matrix_get (fm, k, j);
 608           double lambda_sq = lambda * lambda;
 609           double lambda_4 = lambda_sq * lambda_sq;
 610
 611           l4s += lambda_4;
 612           l2s += lambda_sq;
 613         }
 614       sv += ( fm->size1 * l4s - (l2s * l2s) ) / (fm->size1 * fm->size1 );
 615     }
 616   return sv;
 617 }
 618
 619 static void
 620 rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
 621         const gsl_vector *communalities,
 622         gsl_matrix *result,
 623         gsl_vector *rotated_loadings
 624         )
 625 {
 626   int j, k;
 627   int i;
 628   double prev_sv;
 629
 630   /* First get a normalised version of UNROT */
 631   gsl_matrix *normalised = gsl_matrix_calloc (unrot->size1, unrot->size2);
 632   gsl_matrix *h_sqrt = gsl_matrix_calloc (communalities->size, communalities->size);
 633   gsl_matrix *h_sqrt_inv ;
 634
 635   /* H is the diagonal matrix containing the absolute values of the communalities */
 636   for (i = 0 ; i < communalities->size ; ++i)
 637     {
 638       double *ptr = gsl_matrix_ptr (h_sqrt, i, i);
 639       *ptr = fabs (gsl_vector_get (communalities, i));
 640     }
 641
 642   /* Take the square root of the communalities */
 643   gsl_linalg_cholesky_decomp (h_sqrt);
 644
 645
 646   /* Save a copy of h_sqrt and invert it */
 647   h_sqrt_inv = clone_matrix (h_sqrt);
 648   gsl_linalg_cholesky_decomp (h_sqrt_inv);
 649   gsl_linalg_cholesky_invert (h_sqrt_inv);
 650
 651   /* normalised vertion is H^{1/2} x UNROT */
 652   gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans, 1.0, h_sqrt_inv, unrot, 0.0, normalised);
 653
 654   gsl_matrix_free (h_sqrt_inv);
 655
 656
 657   /* Now perform the rotation iterations */
 658
 659   prev_sv = initial_sv (normalised);
 660   for (i = 0 ; i < cf->iterations ; ++i)
 661     {
 662       double sv = 0.0;
 663       for (j = 0 ; j < normalised->size2; ++j)
 664         {
 665           /* These variables relate to the convergence criterium */
 666           double l4s = 0;
 667           double l2s = 0;
 668
 669           for (k = j + 1 ; k < normalised->size2; ++k)
 670             {
 671               int p;
 672               double a = 0.0;
 673               double b = 0.0;
 674               double c = 0.0;
 675               double d = 0.0;
 676               double x, y;
 677               double phi;
 678
 679               for (p = 0; p < normalised->size1; ++p)
 680                 {
 681                   double jv = gsl_matrix_get (normalised, p, j);
 682                   double kv = gsl_matrix_get (normalised, p, k);
 683
 684                   double u = jv * jv - kv * kv;
 685                   double v = 2 * jv * kv;
 686                   a += u;
 687                   b += v;
 688                   c +=  u * u - v * v;
 689                   d += 2 * u * v;
 690                 }
 691
 692               rotation_coeff [cf->rotation] (&x, &y, a, b, c, d, normalised);
 693
 694               phi = atan2 (x,  y) / 4.0 ;
 695
 696               /* Don't bother rotating if the angle is small */
 697               if ( fabs (sin (phi) ) <= pow (10.0, -15.0))
 698                   continue;
 699
 700               for (p = 0; p < normalised->size1; ++p)
 701                 {
 702                   double *lambda0 = gsl_matrix_ptr (normalised, p, j);
 703                   double *lambda1 = gsl_matrix_ptr (normalised, p, k);
 704                   drot_go (phi, lambda0, lambda1);
 705                 }
 706
 707               /* Calculate the convergence criterium */
 708               {
 709                 double lambda = gsl_matrix_get (normalised, k, j);
 710                 double lambda_sq = lambda * lambda;
 711                 double lambda_4 = lambda_sq * lambda_sq;
 712
 713                 l4s += lambda_4;
 714                 l2s += lambda_sq;
 715               }
 716             }
 717           sv += ( normalised->size1 * l4s - (l2s * l2s) ) / (normalised->size1 * normalised->size1 );
 718         }
 719
 720       if ( fabs (sv - prev_sv) <= cf->rconverge)
 721         break;
 722
 723       prev_sv = sv;
 724     }
 725
 726   gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans, 1.0,
 727                   h_sqrt, normalised,  0.0,   result);
 728
 729   gsl_matrix_free (h_sqrt);
 730   gsl_matrix_free (normalised);
 731
 732
 733   /* reflect negative sums and populate the rotated loadings vector*/
 734   for (i = 0 ; i < result->size2; ++i)
 735     {
 736       double ssq = 0.0;
 737       double sum = 0.0;
 738       for (j = 0 ; j < result->size1; ++j)
 739         {
 740           double s = gsl_matrix_get (result, j, i);
 741           ssq += s * s;
 742           sum += gsl_matrix_get (result, j, i);
 743         }
 744
 745       gsl_vector_set (rotated_loadings, i, ssq);
 746
 747       if ( sum < 0 )
 748         for (j = 0 ; j < result->size1; ++j)
 749           {
 750             double *lambda = gsl_matrix_ptr (result, j, i);
 751             *lambda = - *lambda;
 752           }
 753     }
 754 }
 755
 756
 757 /*
 758   Get an approximation for the factor matrix into FACTORS, and the communalities into COMMUNALITIES.
 759   R is the matrix to be analysed.
 760   WS is a pointer to a structure which must have been initialised with factor_matrix_workspace_init.
 761  */
 762 static void
 763 iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matrix *factors,
 764                        struct factor_matrix_workspace *ws)
 765 {
 766   size_t i;
 767   gsl_matrix_view mv ;
 768
 769   assert (r->size1 == r->size2);
 770   assert (r->size1 == communalities->size);
 771
 772   assert (factors->size1 == r->size1);
 773   assert (factors->size2 == ws->n_factors);
 774
 775   gsl_matrix_memcpy (ws->r, r);
 776
 777   /* Apply Communalities to diagonal of correlation matrix */
 778   for (i = 0 ; i < communalities->size ; ++i)
 779     {
 780       double *x = gsl_matrix_ptr (ws->r, i, i);
 781       *x = gsl_vector_get (communalities, i);
 782     }
 783
 784   gsl_eigen_symmv (ws->r, ws->eval, ws->evec, ws->eigen_ws);
 785
 786   mv = gsl_matrix_submatrix (ws->evec, 0, 0, ws->evec->size1, ws->n_factors);
 787
 788   /* Gamma is the diagonal matrix containing the absolute values of the eigenvalues */
 789   for (i = 0 ; i < ws->n_factors ; ++i)
 790     {
 791       double *ptr = gsl_matrix_ptr (ws->gamma, i, i);
 792       *ptr = fabs (gsl_vector_get (ws->eval, i));
 793     }
 794
 795   /* Take the square root of gamma */
 796   gsl_linalg_cholesky_decomp (ws->gamma);
 797
 798   gsl_blas_dgemm (CblasNoTrans,  CblasNoTrans, 1.0, &mv.matrix, ws->gamma, 0.0, factors);
 799
 800   for (i = 0 ; i < r->size1 ; ++i)
 801     {
 802       double h = the_communality (ws->evec, ws->eval, i, ws->n_factors);
 803       gsl_vector_set (communalities, i, h);
 804     }
 805 }
 806
 807
 808
 809 static bool run_factor (struct dataset *ds, const struct cmd_factor *factor);
 810
 811
 812 int
 813 cmd_factor (struct lexer *lexer, struct dataset *ds)
 814 {
 815   bool extraction_seen = false;
 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           extraction_seen = true;
1001           lex_match (lexer, T_EQUALS);
1002           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1003             {
1004               if (lex_match_id (lexer, "PAF"))
1005                 {
1006                   factor.extraction = EXTRACTION_PAF;
1007                 }
1008               else if (lex_match_id (lexer, "PC"))
1009                 {
1010                   factor.extraction = EXTRACTION_PC;
1011                 }
1012               else if (lex_match_id (lexer, "PA1"))
1013                 {
1014                   factor.extraction = EXTRACTION_PC;
1015                 }
1016               else if (lex_match_id (lexer, "DEFAULT"))
1017                 {
1018                   factor.extraction = EXTRACTION_PC;
1019                 }
1020               else
1021                 {
1022                   lex_error (lexer, NULL);
1023                   goto error;
1024                 }
1025             }
1026         }
1027       else if (lex_match_id (lexer, "FORMAT"))
1028         {
1029           lex_match (lexer, T_EQUALS);
1030           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1031             {
1032               if (lex_match_id (lexer, "SORT"))
1033                 {
1034                   factor.sort = true;
1035                 }
1036               else if (lex_match_id (lexer, "BLANK"))
1037                 {
1038                   if ( lex_force_match (lexer, T_LPAREN))
1039                     {
1040                       lex_force_num (lexer);
1041                       factor.blank = lex_number (lexer);
1042                       lex_get (lexer);
1043                       lex_force_match (lexer, T_RPAREN);
1044                     }
1045                 }
1046               else if (lex_match_id (lexer, "DEFAULT"))
1047                 {
1048                   factor.blank = 0;
1049                   factor.sort = false;
1050                 }
1051               else
1052                 {
1053                   lex_error (lexer, NULL);
1054                   goto error;
1055                 }
1056             }
1057         }
1058       else if (lex_match_id (lexer, "PRINT"))
1059         {
1060           factor.print = 0;
1061           lex_match (lexer, T_EQUALS);
1062           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1063             {
1064               if (lex_match_id (lexer, "UNIVARIATE"))
1065                 {
1066                   factor.print |= PRINT_UNIVARIATE;
1067                 }
1068               else if (lex_match_id (lexer, "DET"))
1069                 {
1070                   factor.print |= PRINT_DETERMINANT;
1071                 }
1072 #if FACTOR_FULLY_IMPLEMENTED
1073               else if (lex_match_id (lexer, "INV"))
1074                 {
1075                 }
1076               else if (lex_match_id (lexer, "AIC"))
1077                 {
1078                 }
1079 #endif
1080               else if (lex_match_id (lexer, "SIG"))
1081                 {
1082                   factor.print |= PRINT_SIG;
1083                 }
1084               else if (lex_match_id (lexer, "CORRELATION"))
1085                 {
1086                   factor.print |= PRINT_CORRELATION;
1087                 }
1088 #if FACTOR_FULLY_IMPLEMENTED
1089               else if (lex_match_id (lexer, "COVARIANCE"))
1090                 {
1091                 }
1092 #endif
1093               else if (lex_match_id (lexer, "ROTATION"))
1094                 {
1095                   factor.print |= PRINT_ROTATION;
1096                 }
1097               else if (lex_match_id (lexer, "EXTRACTION"))
1098                 {
1099                   factor.print |= PRINT_EXTRACTION;
1100                 }
1101               else if (lex_match_id (lexer, "INITIAL"))
1102                 {
1103                   factor.print |= PRINT_INITIAL;
1104                 }
1105               else if (lex_match_id (lexer, "KMO"))
1106                 {
1107                   factor.print |= PRINT_KMO;
1108                 }
1109 #if FACTOR_FULLY_IMPLEMENTED
1110               else if (lex_match_id (lexer, "REPR"))
1111                 {
1112                 }
1113               else if (lex_match_id (lexer, "FSCORE"))
1114                 {
1115                 }
1116 #endif
1117               else if (lex_match (lexer, T_ALL))
1118                 {
1119                   factor.print = 0xFFFF;
1120                 }
1121               else if (lex_match_id (lexer, "DEFAULT"))
1122                 {
1123                   factor.print |= PRINT_INITIAL ;
1124                   factor.print |= PRINT_EXTRACTION ;
1125                   factor.print |= PRINT_ROTATION ;
1126                 }
1127               else
1128                 {
1129                   lex_error (lexer, NULL);
1130                   goto error;
1131                 }
1132             }
1133         }
1134       else if (lex_match_id (lexer, "MISSING"))
1135         {
1136           lex_match (lexer, T_EQUALS);
1137           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
1138             {
1139               if (lex_match_id (lexer, "INCLUDE"))
1140                 {
1141                   factor.exclude = MV_SYSTEM;
1142                 }
1143               else if (lex_match_id (lexer, "EXCLUDE"))
1144                 {
1145                   factor.exclude = MV_ANY;
1146                 }
1147               else if (lex_match_id (lexer, "LISTWISE"))
1148                 {
1149                   factor.missing_type = MISS_LISTWISE;
1150                 }
1151               else if (lex_match_id (lexer, "PAIRWISE"))
1152                 {
1153                   factor.missing_type = MISS_PAIRWISE;
1154                 }
1155               else if (lex_match_id (lexer, "MEANSUB"))
1156                 {
1157                   factor.missing_type = MISS_MEANSUB;
1158                 }
1159               else
1160                 {
1161                   lex_error (lexer, NULL);
1162                   goto error;
1163                 }
1164             }
1165         }
1166       else
1167         {
1168           lex_error (lexer, NULL);
1169           goto error;
1170         }
1171     }
1172
1173   if ( factor.rotation == ROT_NONE )
1174     factor.print &= ~PRINT_ROTATION;
1175
1176   if ( ! run_factor (ds, &factor))
1177     goto error;
1178
1179   free (factor.vars);
1180   return CMD_SUCCESS;
1181
1182  error:
1183   free (factor.vars);
1184   return CMD_FAILURE;
1185 }
1186
1187 static void do_factor (const struct cmd_factor *factor, struct casereader *group);
1188
1189
1190 static bool
1191 run_factor (struct dataset *ds, const struct cmd_factor *factor)
1192 {
1193   struct dictionary *dict = dataset_dict (ds);
1194   bool ok;
1195   struct casereader *group;
1196
1197   struct casegrouper *grouper = casegrouper_create_splits (proc_open (ds), dict);
1198
1199   while (casegrouper_get_next_group (grouper, &group))
1200     {
1201       if ( factor->missing_type == MISS_LISTWISE )
1202         group  = casereader_create_filter_missing (group, factor->vars, factor->n_vars,
1203                                                    factor->exclude,
1204                                                    NULL,  NULL);
1205       do_factor (factor, group);
1206     }
1207
1208   ok = casegrouper_destroy (grouper);
1209   ok = proc_commit (ds) && ok;
1210
1211   return ok;
1212 }
1213
1214
1215 /* Return the communality of variable N, calculated to N_FACTORS */
1216 static double
1217 the_communality (const gsl_matrix *evec, const gsl_vector *eval, int n, int n_factors)
1218 {
1219   size_t i;
1220
1221   double comm = 0;
1222
1223   assert (n >= 0);
1224   assert (n < eval->size);
1225   assert (n < evec->size1);
1226   assert (n_factors <= eval->size);
1227
1228   for (i = 0 ; i < n_factors; ++i)
1229     {
1230       double evali = fabs (gsl_vector_get (eval, i));
1231
1232       double eveci = gsl_matrix_get (evec, n, i);
1233
1234       comm += pow2 (eveci) * evali;
1235     }
1236
1237   return comm;
1238 }
1239
1240 /* Return the communality of variable N, calculated to N_FACTORS */
1241 static double
1242 communality (struct idata *idata, int n, int n_factors)
1243 {
1244   return the_communality (idata->evec, idata->eval, n, n_factors);
1245 }
1246
1247
1248 static void
1249 show_scree (const struct cmd_factor *f, struct idata *idata)
1250 {
1251   struct scree *s;
1252   const char *label ;
1253
1254   if ( !(f->plot & PLOT_SCREE) )
1255     return;
1256
1257
1258   label = f->extraction == EXTRACTION_PC ? _("Component Number") : _("Factor Number");
1259
1260   s = scree_create (idata->eval, label);
1261
1262   scree_submit (s);
1263 }
1264
1265 static void
1266 show_communalities (const struct cmd_factor * factor,
1267                     const gsl_vector *initial, const gsl_vector *extracted)
1268 {
1269   int i;
1270   int c = 0;
1271   const int heading_columns = 1;
1272   int nc = heading_columns;
1273   const int heading_rows = 1;
1274   const int nr = heading_rows + factor->n_vars;
1275   struct tab_table *t;
1276
1277   if (factor->print & PRINT_EXTRACTION)
1278     nc++;
1279
1280   if (factor->print & PRINT_INITIAL)
1281     nc++;
1282
1283   /* No point having a table with only headings */
1284   if (nc <= 1)
1285     return;
1286
1287   t = tab_create (nc, nr);
1288
1289   tab_title (t, _("Communalities"));
1290
1291   tab_headers (t, heading_columns, 0, heading_rows, 0);
1292
1293   c = 1;
1294   if (factor->print & PRINT_INITIAL)
1295     tab_text (t, c++, 0, TAB_CENTER | TAT_TITLE, _("Initial"));
1296
1297   if (factor->print & PRINT_EXTRACTION)
1298     tab_text (t, c++, 0, TAB_CENTER | TAT_TITLE, _("Extraction"));
1299
1300   /* Outline the box */
1301   tab_box (t,
1302            TAL_2, TAL_2,
1303            -1, -1,
1304            0, 0,
1305            nc - 1, nr - 1);
1306
1307   /* Vertical lines */
1308   tab_box (t,
1309            -1, -1,
1310            -1, TAL_1,
1311            heading_columns, 0,
1312            nc - 1, nr - 1);
1313
1314   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1315   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1316
1317   for (i = 0 ; i < factor->n_vars; ++i)
1318     {
1319       c = 0;
1320       tab_text (t, c++, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[i]));
1321
1322       if (factor->print & PRINT_INITIAL)
1323         tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL);
1324
1325       if (factor->print & PRINT_EXTRACTION)
1326         tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL);
1327     }
1328
1329   tab_submit (t);
1330 }
1331
1332
1333 static void
1334 show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const char *title, const gsl_matrix *fm)
1335 {
1336   int i;
1337   const int n_factors = idata->n_extractions;
1338
1339   const int heading_columns = 1;
1340   const int heading_rows = 2;
1341   const int nr = heading_rows + factor->n_vars;
1342   const int nc = heading_columns + n_factors;
1343   gsl_permutation *perm;
1344
1345   struct tab_table *t = tab_create (nc, nr);
1346
1347   /*
1348   if ( factor->extraction == EXTRACTION_PC )
1349     tab_title (t, _("Component Matrix"));
1350   else
1351     tab_title (t, _("Factor Matrix"));
1352   */
1353
1354   tab_title (t, "%s", title);
1355
1356   tab_headers (t, heading_columns, 0, heading_rows, 0);
1357
1358   if ( factor->extraction == EXTRACTION_PC )
1359     tab_joint_text (t,
1360                     1, 0,
1361                     nc - 1, 0,
1362                     TAB_CENTER | TAT_TITLE, _("Component"));
1363   else
1364     tab_joint_text (t,
1365                     1, 0,
1366                     nc - 1, 0,
1367                     TAB_CENTER | TAT_TITLE, _("Factor"));
1368
1369
1370   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1371
1372
1373   /* Outline the box */
1374   tab_box (t,
1375            TAL_2, TAL_2,
1376            -1, -1,
1377            0, 0,
1378            nc - 1, nr - 1);
1379
1380   /* Vertical lines */
1381   tab_box (t,
1382            -1, -1,
1383            -1, TAL_1,
1384            heading_columns, 1,
1385            nc - 1, nr - 1);
1386
1387   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1388   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1389
1390
1391   /* Initialise to the identity permutation */
1392   perm = gsl_permutation_calloc (factor->n_vars);
1393
1394   if ( factor->sort)
1395     sort_matrix_indirect (fm, perm);
1396
1397   for (i = 0 ; i < n_factors; ++i)
1398     {
1399       tab_text_format (t, heading_columns + i, 1, TAB_CENTER | TAT_TITLE, _("%d"), i + 1);
1400     }
1401
1402   for (i = 0 ; i < factor->n_vars; ++i)
1403     {
1404       int j;
1405       const int matrix_row = perm->data[i];
1406       tab_text (t, 0, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[matrix_row]));
1407
1408       for (j = 0 ; j < n_factors; ++j)
1409         {
1410           double x = gsl_matrix_get (fm, matrix_row, j);
1411
1412           if ( fabs (x) < factor->blank)
1413             continue;
1414
1415           tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL);
1416         }
1417     }
1418
1419   gsl_permutation_free (perm);
1420
1421   tab_submit (t);
1422 }
1423
1424
1425 static void
1426 show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
1427                          const gsl_vector *initial_eigenvalues,
1428                          const gsl_vector *extracted_eigenvalues,
1429                          const gsl_vector *rotated_loadings)
1430 {
1431   size_t i;
1432   int c = 0;
1433   const int heading_columns = 1;
1434   const int heading_rows = 2;
1435   const int nr = heading_rows + factor->n_vars;
1436
1437   struct tab_table *t ;
1438
1439   double i_total = 0.0;
1440   double i_cum = 0.0;
1441
1442   double e_total = 0.0;
1443   double e_cum = 0.0;
1444
1445   double r_cum = 0.0;
1446
1447   int nc = heading_columns;
1448
1449   if (factor->print & PRINT_EXTRACTION)
1450     nc += 3;
1451
1452   if (factor->print & PRINT_INITIAL)
1453     nc += 3;
1454
1455   if (factor->print & PRINT_ROTATION)
1456     nc += 3;
1457
1458   /* No point having a table with only headings */
1459   if ( nc <= heading_columns)
1460     return;
1461
1462   t = tab_create (nc, nr);
1463
1464   tab_title (t, _("Total Variance Explained"));
1465
1466   tab_headers (t, heading_columns, 0, heading_rows, 0);
1467
1468   /* Outline the box */
1469   tab_box (t,
1470            TAL_2, TAL_2,
1471            -1, -1,
1472            0, 0,
1473            nc - 1, nr - 1);
1474
1475   /* Vertical lines */
1476   tab_box (t,
1477            -1, -1,
1478            -1, TAL_1,
1479            heading_columns, 0,
1480            nc - 1, nr - 1);
1481
1482   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1483   tab_hline (t, TAL_1, 1, nc - 1, 1);
1484
1485   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1486
1487
1488   if ( factor->extraction == EXTRACTION_PC)
1489     tab_text (t, 0, 1, TAB_LEFT | TAT_TITLE, _("Component"));
1490   else
1491     tab_text (t, 0, 1, TAB_LEFT | TAT_TITLE, _("Factor"));
1492
1493   c = 1;
1494   if (factor->print & PRINT_INITIAL)
1495     {
1496       tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Initial Eigenvalues"));
1497       c += 3;
1498     }
1499
1500   if (factor->print & PRINT_EXTRACTION)
1501     {
1502       tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Extraction Sums of Squared Loadings"));
1503       c += 3;
1504     }
1505
1506   if (factor->print & PRINT_ROTATION)
1507     {
1508       tab_joint_text (t, c, 0, c + 2, 0, TAB_CENTER | TAT_TITLE, _("Rotation Sums of Squared Loadings"));
1509       c += 3;
1510     }
1511
1512   for (i = 0; i < (nc - heading_columns) / 3 ; ++i)
1513     {
1514       tab_text (t, i * 3 + 1, 1, TAB_CENTER | TAT_TITLE, _("Total"));
1515       /* xgettext:no-c-format */
1516       tab_text (t, i * 3 + 2, 1, TAB_CENTER | TAT_TITLE, _("% of Variance"));
1517       tab_text (t, i * 3 + 3, 1, TAB_CENTER | TAT_TITLE, _("Cumulative %"));
1518
1519       tab_vline (t, TAL_2, heading_columns + i * 3, 0, nr - 1);
1520     }
1521
1522   for (i = 0 ; i < initial_eigenvalues->size; ++i)
1523     i_total += gsl_vector_get (initial_eigenvalues, i);
1524
1525   if ( factor->extraction == EXTRACTION_PAF)
1526     {
1527       e_total = factor->n_vars;
1528     }
1529   else
1530     {
1531       e_total = i_total;
1532     }
1533
1534   for (i = 0 ; i < factor->n_vars; ++i)
1535     {
1536       const double i_lambda = gsl_vector_get (initial_eigenvalues, i);
1537       double i_percent = 100.0 * i_lambda / i_total ;
1538
1539       const double e_lambda = gsl_vector_get (extracted_eigenvalues, i);
1540       double e_percent = 100.0 * e_lambda / e_total ;
1541
1542       c = 0;
1543
1544       tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%zu"), i + 1);
1545
1546       i_cum += i_percent;
1547       e_cum += e_percent;
1548
1549       /* Initial Eigenvalues */
1550       if (factor->print & PRINT_INITIAL)
1551       {
1552         tab_double (t, c++, i + heading_rows, 0, i_lambda, NULL);
1553         tab_double (t, c++, i + heading_rows, 0, i_percent, NULL);
1554         tab_double (t, c++, i + heading_rows, 0, i_cum, NULL);
1555       }
1556
1557
1558       if (factor->print & PRINT_EXTRACTION)
1559         {
1560           if (i < idata->n_extractions)
1561             {
1562               /* Sums of squared loadings */
1563               tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL);
1564               tab_double (t, c++, i + heading_rows, 0, e_percent, NULL);
1565               tab_double (t, c++, i + heading_rows, 0, e_cum, NULL);
1566             }
1567         }
1568
1569       if (rotated_loadings != NULL)
1570         {
1571           const double r_lambda = gsl_vector_get (rotated_loadings, i);
1572           double r_percent = 100.0 * r_lambda / e_total ;
1573
1574           if (factor->print & PRINT_ROTATION)
1575             {
1576               if (i < idata->n_extractions)
1577                 {
1578                   r_cum += r_percent;
1579                   tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL);
1580                   tab_double (t, c++, i + heading_rows, 0, r_percent, NULL);
1581                   tab_double (t, c++, i + heading_rows, 0, r_cum, NULL);
1582                 }
1583             }
1584         }
1585     }
1586
1587   tab_submit (t);
1588 }
1589
1590
1591 static void
1592 show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata)
1593 {
1594   struct tab_table *t ;
1595   size_t i, j;
1596   int y_pos_corr = -1;
1597   int y_pos_sig = -1;
1598   int suffix_rows = 0;
1599
1600   const int heading_rows = 1;
1601   const int heading_columns = 2;
1602
1603   int nc = heading_columns ;
1604   int nr = heading_rows ;
1605   int n_data_sets = 0;
1606
1607   if (factor->print & PRINT_CORRELATION)
1608     {
1609       y_pos_corr = n_data_sets;
1610       n_data_sets++;
1611       nc = heading_columns + factor->n_vars;
1612     }
1613
1614   if (factor->print & PRINT_SIG)
1615     {
1616       y_pos_sig = n_data_sets;
1617       n_data_sets++;
1618       nc = heading_columns + factor->n_vars;
1619     }
1620
1621   nr += n_data_sets * factor->n_vars;
1622
1623   if (factor->print & PRINT_DETERMINANT)
1624     suffix_rows = 1;
1625
1626   /* If the table would contain only headings, don't bother rendering it */
1627   if (nr <= heading_rows && suffix_rows == 0)
1628     return;
1629
1630   t = tab_create (nc, nr + suffix_rows);
1631
1632   tab_title (t, _("Correlation Matrix"));
1633
1634   tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
1635
1636   if (nr > heading_rows)
1637     {
1638       tab_headers (t, heading_columns, 0, heading_rows, 0);
1639
1640       tab_vline (t, TAL_2, 2, 0, nr - 1);
1641
1642       /* Outline the box */
1643       tab_box (t,
1644                TAL_2, TAL_2,
1645                -1, -1,
1646                0, 0,
1647                nc - 1, nr - 1);
1648
1649       /* Vertical lines */
1650       tab_box (t,
1651                -1, -1,
1652                -1, TAL_1,
1653                heading_columns, 0,
1654                nc - 1, nr - 1);
1655
1656
1657       for (i = 0; i < factor->n_vars; ++i)
1658         tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i]));
1659
1660
1661       for (i = 0 ; i < n_data_sets; ++i)
1662         {
1663           int y = heading_rows + i * factor->n_vars;
1664           size_t v;
1665           for (v = 0; v < factor->n_vars; ++v)
1666             tab_text (t, 1, y + v, TAT_TITLE, var_to_string (factor->vars[v]));
1667
1668           tab_hline (t, TAL_1, 0, nc - 1, y);
1669         }
1670
1671       if (factor->print & PRINT_CORRELATION)
1672         {
1673           const double y = heading_rows + y_pos_corr;
1674           tab_text (t, 0, y, TAT_TITLE, _("Correlations"));
1675
1676           for (i = 0; i < factor->n_vars; ++i)
1677             {
1678               for (j = 0; j < factor->n_vars; ++j)
1679                 tab_double (t, heading_columns + i,  y + j, 0, gsl_matrix_get (idata->corr, i, j), NULL);
1680             }
1681         }
1682
1683       if (factor->print & PRINT_SIG)
1684         {
1685           const double y = heading_rows + y_pos_sig * factor->n_vars;
1686           tab_text (t, 0, y, TAT_TITLE, _("Sig. (1-tailed)"));
1687
1688           for (i = 0; i < factor->n_vars; ++i)
1689             {
1690               for (j = 0; j < factor->n_vars; ++j)
1691                 {
1692                   double rho = gsl_matrix_get (idata->corr, i, j);
1693                   double w = gsl_matrix_get (idata->n, i, j);
1694
1695                   if (i == j)
1696                     continue;
1697
1698                   tab_double (t, heading_columns + i,  y + j, 0, significance_of_correlation (rho, w), NULL);
1699                 }
1700             }
1701         }
1702     }
1703
1704   if (factor->print & PRINT_DETERMINANT)
1705     {
1706       tab_text (t, 0, nr, TAB_LEFT | TAT_TITLE, _("Determinant"));
1707
1708       tab_double (t, 1, nr, 0, idata->detR, NULL);
1709     }
1710
1711   tab_submit (t);
1712 }
1713
1714
1715
1716 static void
1717 do_factor (const struct cmd_factor *factor, struct casereader *r)
1718 {
1719   struct ccase *c;
1720   const gsl_matrix *var_matrix;
1721   const gsl_matrix *mean_matrix;
1722
1723   const gsl_matrix *analysis_matrix;
1724   struct idata *idata = idata_alloc (factor->n_vars);
1725
1726   struct covariance *cov = covariance_1pass_create (factor->n_vars, factor->vars,
1727                                               factor->wv, factor->exclude);
1728
1729   for ( ; (c = casereader_read (r) ); case_unref (c))
1730     {
1731       covariance_accumulate (cov, c);
1732     }
1733
1734   idata->cov = covariance_calculate (cov);
1735
1736   if (idata->cov == NULL)
1737     {
1738       msg (MW, _("The dataset contains no complete observations. No analysis will be performed."));
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