Added result_class parameter to tab_double and updated all callers. Removed tab_fixed
[pspp] / src / language / stats / factor.c
index 83dc13de7ed5a28f92e16cd8e24e2a2a02c77b25..4094a22b432467ce25dba36a3b3f7c55f2861013 100644 (file)
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2009, 2010 Free Software Foundation, Inc.
+   Copyright (C) 2009, 2010, 2011, 2012 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
 
 #include <config.h>
 
-
 #include <gsl/gsl_vector.h>
 #include <gsl/gsl_linalg.h>
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_eigen.h> 
 #include <gsl/gsl_blas.h> 
 #include <gsl/gsl_sort_vector.h>
-
-#include <math/covariance.h>
-
-#include <math/correlation.h>
-#include <math/moments.h>
-#include <data/procedure.h>
-#include <language/lexer/variable-parser.h>
-#include <language/lexer/value-parser.h>
-#include <language/command.h>
-#include <language/lexer/lexer.h>
-
-#include <data/casegrouper.h>
-#include <data/casereader.h>
-#include <data/casewriter.h>
-#include <data/dictionary.h>
-#include <data/format.h>
-#include <data/subcase.h>
-
-#include <libpspp/misc.h>
-#include <libpspp/message.h>
-
-#include <output/tab.h>
-
-#include <output/charts/scree.h>
-#include <output/chart-item.h>
+#include <gsl/gsl_cdf.h>
+
+#include "data/casegrouper.h"
+#include "data/casereader.h"
+#include "data/casewriter.h"
+#include "data/dataset.h"
+#include "data/dictionary.h"
+#include "data/format.h"
+#include "data/subcase.h"
+#include "language/command.h"
+#include "language/lexer/lexer.h"
+#include "language/lexer/value-parser.h"
+#include "language/lexer/variable-parser.h"
+#include "libpspp/cast.h"
+#include "libpspp/message.h"
+#include "libpspp/misc.h"
+#include "math/correlation.h"
+#include "math/covariance.h"
+#include "math/moments.h"
+#include "output/chart-item.h"
+#include "output/charts/scree.h"
+#include "output/tab.h"
 
 #include "gettext.h"
 #define _(msgid) gettext (msgid)
@@ -156,12 +152,13 @@ struct cmd_factor
   enum extraction_method extraction;
   enum plot_opts plot;
   enum rotation_type rotation;
+  int rotation_iterations;
 
   /* Extraction Criteria */
   int n_factors;
   double min_eigen;
   double econverge;
-  int iterations;
+  int extraction_iterations;
 
   double rconverge;
 
@@ -175,7 +172,7 @@ struct idata
   /* Intermediate values used in calculation */
 
   const gsl_matrix *corr ;  /* The correlation matrix */
-  const gsl_matrix *cov ;   /* The covariance matrix */
+  gsl_matrix *cov ;         /* The covariance matrix */
   const gsl_matrix *n ;     /* Matrix of number of samples */
 
   gsl_vector *eval ;  /* The eigenvalues */
@@ -184,6 +181,8 @@ struct idata
   int n_extractions;
 
   gsl_vector *msr ;  /* Multiple Squared Regressions */
+
+  double detR;  /* The determinant of the correlation matrix */
 };
 
 static struct idata *
@@ -206,11 +205,64 @@ idata_free (struct idata *id)
   gsl_vector_free (id->msr);
   gsl_vector_free (id->eval);
   gsl_matrix_free (id->evec);
+  if (id->cov != NULL)
+    gsl_matrix_free (id->cov);
+  if (id->corr != NULL)
+    gsl_matrix_free (CONST_CAST (gsl_matrix *, id->corr));
 
   free (id);
 }
 
 
+static gsl_matrix *
+anti_image (const gsl_matrix *m)
+{
+  int i, j;
+  gsl_matrix *a;
+  assert (m->size1 == m->size2);
+
+  a = gsl_matrix_alloc (m->size1, m->size2);
+  
+  for (i = 0; i < m->size1; ++i)
+    {
+      for (j = 0; j < m->size2; ++j)
+       {
+         double *p = gsl_matrix_ptr (a, i, j);
+         *p = gsl_matrix_get (m, i, j);
+         *p /= gsl_matrix_get (m, i, i);
+         *p /= gsl_matrix_get (m, j, j);
+       }
+    }
+
+  return a;
+}
+
+
+/* Return the sum of all the elements excluding row N */
+static double
+ssq_od_n (const gsl_matrix *m, int n)
+{
+  int i, j;
+  double ss = 0;
+  assert (m->size1 == m->size2);
+
+  assert (n < m->size1);
+  
+  for (i = 0; i < m->size1; ++i)
+    {
+      if (i == n ) continue;
+      for (j = 0; j < m->size2; ++j)
+       {
+         ss += pow2 (gsl_matrix_get (m, i, j));
+       }
+    }
+
+  return ss;
+}
+
+
+
+#if 0
 static void
 dump_matrix (const gsl_matrix *m)
 {
@@ -224,7 +276,6 @@ dump_matrix (const gsl_matrix *m)
     }
 }
 
-
 static void
 dump_matrix_permute (const gsl_matrix *m, const gsl_permutation *p)
 {
@@ -249,6 +300,7 @@ dump_vector (const gsl_vector *v)
     }
   printf ("\n");
 }
+#endif
 
 
 static int 
@@ -607,7 +659,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
   /* Now perform the rotation iterations */
 
   prev_sv = initial_sv (normalised);
-  for (i = 0 ; i < cf->iterations ; ++i)
+  for (i = 0 ; i < cf->rotation_iterations ; ++i)
     {
       double sv = 0.0;
       for (j = 0 ; j < normalised->size2; ++j)
@@ -677,6 +729,7 @@ rotate (const struct cmd_factor *cf, const gsl_matrix *unrot,
                  h_sqrt, normalised,  0.0,   result);
 
   gsl_matrix_free (h_sqrt);
+  gsl_matrix_free (normalised);
 
 
   /* reflect negative sums and populate the rotated loadings vector*/
@@ -761,9 +814,8 @@ static bool run_factor (struct dataset *ds, const struct cmd_factor *factor);
 int
 cmd_factor (struct lexer *lexer, struct dataset *ds)
 {
-  bool extraction_seen = false;
   const struct dictionary *dict = dataset_dict (ds);
-
+  int n_iterations = 25;
   struct cmd_factor factor;
   factor.n_vars = 0;
   factor.vars = NULL;
@@ -774,7 +826,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
   factor.extraction = EXTRACTION_PC;
   factor.n_factors = 0;
   factor.min_eigen = SYSMIS;
-  factor.iterations = 25;
+  factor.extraction_iterations = 25;
+  factor.rotation_iterations = 25;
   factor.econverge = 0.001;
 
   factor.blank = 0;
@@ -786,14 +839,14 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
 
   factor.wv = dict_get_weight (dict);
 
-  lex_match (lexer, '/');
+  lex_match (lexer, T_SLASH);
 
   if (!lex_force_match_id (lexer, "VARIABLES"))
     {
       goto error;
     }
 
-  lex_match (lexer, '=');
+  lex_match (lexer, T_EQUALS);
 
   if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars,
                              PV_NO_DUPLICATE | PV_NUMERIC))
@@ -802,14 +855,14 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
   if (factor.n_vars < 2)
     msg (MW, _("Factor analysis on a single variable is not useful."));
 
-  while (lex_token (lexer) != '.')
+  while (lex_token (lexer) != T_ENDCMD)
     {
-      lex_match (lexer, '/');
+      lex_match (lexer, T_SLASH);
 
       if (lex_match_id (lexer, "PLOT"))
        {
-          lex_match (lexer, '=');
-          while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
+          lex_match (lexer, T_EQUALS);
+          while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
            {
              if (lex_match_id (lexer, "EIGEN"))
                {
@@ -829,8 +882,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
        }
       else if (lex_match_id (lexer, "METHOD"))
        {
-          lex_match (lexer, '=');
-          while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
+          lex_match (lexer, T_EQUALS);
+          while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
            {
              if (lex_match_id (lexer, "COVARIANCE"))
                {
@@ -849,8 +902,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
        }
       else if (lex_match_id (lexer, "ROTATION"))
        {
-          lex_match (lexer, '=');
-          while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
+          lex_match (lexer, T_EQUALS);
+          while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
            {
              /* VARIMAX and DEFAULT are defaults */
              if (lex_match_id (lexer, "VARIMAX") || lex_match_id (lexer, "DEFAULT"))
@@ -875,67 +928,68 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                  goto error;
                }
            }
+          factor.rotation_iterations = n_iterations;
        }
       else if (lex_match_id (lexer, "CRITERIA"))
        {
-          lex_match (lexer, '=');
-          while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
+          lex_match (lexer, T_EQUALS);
+          while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
            {
              if (lex_match_id (lexer, "FACTORS"))
                {
-                 if ( lex_force_match (lexer, '('))
+                 if ( lex_force_match (lexer, T_LPAREN))
                    {
                      lex_force_int (lexer);
                      factor.n_factors = lex_integer (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, ')');
+                     lex_force_match (lexer, T_RPAREN);
                    }
                }
              else if (lex_match_id (lexer, "MINEIGEN"))
                {
-                 if ( lex_force_match (lexer, '('))
+                 if ( lex_force_match (lexer, T_LPAREN))
                    {
                      lex_force_num (lexer);
                      factor.min_eigen = lex_number (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, ')');
+                     lex_force_match (lexer, T_RPAREN);
                    }
                }
              else if (lex_match_id (lexer, "ECONVERGE"))
                {
-                 if ( lex_force_match (lexer, '('))
+                 if ( lex_force_match (lexer, T_LPAREN))
                    {
                      lex_force_num (lexer);
                      factor.econverge = lex_number (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, ')');
+                     lex_force_match (lexer, T_RPAREN);
                    }
                }
              else if (lex_match_id (lexer, "RCONVERGE"))
                {
-                 if ( lex_force_match (lexer, '('))
+                 if ( lex_force_match (lexer, T_LPAREN))
                    {
                      lex_force_num (lexer);
                      factor.rconverge = lex_number (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, ')');
+                     lex_force_match (lexer, T_RPAREN);
                    }
                }
              else if (lex_match_id (lexer, "ITERATE"))
                {
-                 if ( lex_force_match (lexer, '('))
+                 if ( lex_force_match (lexer, T_LPAREN))
                    {
                      lex_force_int (lexer);
-                     factor.iterations = lex_integer (lexer);
+                     n_iterations = lex_integer (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, ')');
+                     lex_force_match (lexer, T_RPAREN);
                    }
                }
              else if (lex_match_id (lexer, "DEFAULT"))
                {
                  factor.n_factors = 0;
                  factor.min_eigen = 1;
-                 factor.iterations = 25;
+                 n_iterations = 25;
                }
              else
                {
@@ -946,9 +1000,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
        }
       else if (lex_match_id (lexer, "EXTRACTION"))
        {
-         extraction_seen = true;
-          lex_match (lexer, '=');
-          while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
+          lex_match (lexer, T_EQUALS);
+          while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
            {
              if (lex_match_id (lexer, "PAF"))
                {
@@ -972,11 +1025,12 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                  goto error;
                }
            }
+          factor.extraction_iterations = n_iterations;
        }
       else if (lex_match_id (lexer, "FORMAT"))
        {
-          lex_match (lexer, '=');
-          while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
+          lex_match (lexer, T_EQUALS);
+          while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
            {
              if (lex_match_id (lexer, "SORT"))
                {
@@ -984,12 +1038,12 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                }
              else if (lex_match_id (lexer, "BLANK"))
                {
-                 if ( lex_force_match (lexer, '('))
+                 if ( lex_force_match (lexer, T_LPAREN))
                    {
                      lex_force_num (lexer);
                      factor.blank = lex_number (lexer);
                      lex_get (lexer);
-                     lex_force_match (lexer, ')');
+                     lex_force_match (lexer, T_RPAREN);
                    }
                }
              else if (lex_match_id (lexer, "DEFAULT"))
@@ -1007,8 +1061,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
       else if (lex_match_id (lexer, "PRINT"))
        {
          factor.print = 0;
-          lex_match (lexer, '=');
-          while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
+          lex_match (lexer, T_EQUALS);
+          while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
             {
               if (lex_match_id (lexer, "UNIVARIATE"))
                {
@@ -1051,10 +1105,11 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                {
                  factor.print |= PRINT_INITIAL;
                }
-#if FACTOR_FULLY_IMPLEMENTED
              else if (lex_match_id (lexer, "KMO"))
                {
+                 factor.print |= PRINT_KMO;
                }
+#if FACTOR_FULLY_IMPLEMENTED
              else if (lex_match_id (lexer, "REPR"))
                {
                }
@@ -1081,8 +1136,8 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
        }
       else if (lex_match_id (lexer, "MISSING"))
         {
-          lex_match (lexer, '=');
-          while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
+          lex_match (lexer, T_EQUALS);
+          while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
             {
              if (lex_match_id (lexer, "INCLUDE"))
                {
@@ -1268,10 +1323,10 @@ show_communalities (const struct cmd_factor * factor,
       tab_text (t, c++, i + heading_rows, TAT_TITLE, var_to_string (factor->vars[i]));
 
       if (factor->print & PRINT_INITIAL)
-       tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL);
+       tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (initial, i), NULL, RC_OTHER);
 
       if (factor->print & PRINT_EXTRACTION)
-       tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL);
+       tab_double (t, c++, i + heading_rows, 0, gsl_vector_get (extracted, i), NULL, RC_OTHER);
     }
 
   tab_submit (t);
@@ -1299,7 +1354,7 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const
     tab_title (t, _("Factor Matrix"));
   */
 
-  tab_title (t, title);
+  tab_title (t, "%s", title);
 
   tab_headers (t, heading_columns, 0, heading_rows, 0);
 
@@ -1360,7 +1415,7 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const
          if ( fabs (x) < factor->blank)
            continue;
 
-         tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL);
+         tab_double (t, heading_columns + j, heading_rows + i, 0, x, NULL, RC_OTHER);
        }
     }
 
@@ -1487,23 +1542,19 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
       const double e_lambda = gsl_vector_get (extracted_eigenvalues, i);
       double e_percent = 100.0 * e_lambda / e_total ;
 
-      const double r_lambda = gsl_vector_get (rotated_loadings, i);
-      double r_percent = 100.0 * r_lambda / e_total ;
-
       c = 0;
 
-      tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%d"), i + 1);
+      tab_text_format (t, c++, i + heading_rows, TAB_LEFT | TAT_TITLE, _("%zu"), i + 1);
 
       i_cum += i_percent;
       e_cum += e_percent;
-      r_cum += r_percent;
 
       /* Initial Eigenvalues */
       if (factor->print & PRINT_INITIAL)
       {
-       tab_double (t, c++, i + heading_rows, 0, i_lambda, NULL);
-       tab_double (t, c++, i + heading_rows, 0, i_percent, NULL);
-       tab_double (t, c++, i + heading_rows, 0, i_cum, NULL);
+       tab_double (t, c++, i + heading_rows, 0, i_lambda, NULL, RC_OTHER);
+       tab_double (t, c++, i + heading_rows, 0, i_percent, NULL, RC_OTHER);
+       tab_double (t, c++, i + heading_rows, 0, i_cum, NULL, RC_OTHER);
       }
 
 
@@ -1512,22 +1563,28 @@ show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
          if (i < idata->n_extractions)
            {
              /* Sums of squared loadings */
-             tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL);
-             tab_double (t, c++, i + heading_rows, 0, e_percent, NULL);
-             tab_double (t, c++, i + heading_rows, 0, e_cum, NULL);
+             tab_double (t, c++, i + heading_rows, 0, e_lambda, NULL, RC_OTHER);
+             tab_double (t, c++, i + heading_rows, 0, e_percent, NULL, RC_OTHER);
+             tab_double (t, c++, i + heading_rows, 0, e_cum, NULL, RC_OTHER);
            }
        }
 
-      if (factor->print & PRINT_ROTATION)
-       {
-         if (i < idata->n_extractions)
-           {
-             tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL);
-             tab_double (t, c++, i + heading_rows, 0, r_percent, NULL);
-             tab_double (t, c++, i + heading_rows, 0, r_cum, NULL);
-           }
-      }
+      if (rotated_loadings != NULL)
+        {
+          const double r_lambda = gsl_vector_get (rotated_loadings, i);
+          double r_percent = 100.0 * r_lambda / e_total ;
 
+          if (factor->print & PRINT_ROTATION)
+            {
+              if (i < idata->n_extractions)
+                {
+                  r_cum += r_percent;
+                  tab_double (t, c++, i + heading_rows, 0, r_lambda, NULL, RC_OTHER);
+                  tab_double (t, c++, i + heading_rows, 0, r_percent, NULL, RC_OTHER);
+                  tab_double (t, c++, i + heading_rows, 0, r_cum, NULL, RC_OTHER);
+                }
+            }
+        }
     }
 
   tab_submit (t);
@@ -1622,14 +1679,14 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id
          for (i = 0; i < factor->n_vars; ++i)
            {
              for (j = 0; j < factor->n_vars; ++j)
-               tab_double (t, heading_columns + i,  y + j, 0, gsl_matrix_get (idata->corr, i, j), NULL);
+               tab_double (t, heading_columns + i,  y + j, 0, gsl_matrix_get (idata->corr, i, j), NULL, RC_OTHER);
            }
        }
 
       if (factor->print & PRINT_SIG)
        {
          const double y = heading_rows + y_pos_sig * factor->n_vars;
-         tab_text (t, 0, y, TAT_TITLE, _("Sig. 1-tailed"));
+         tab_text (t, 0, y, TAT_TITLE, _("Sig. (1-tailed)"));
 
          for (i = 0; i < factor->n_vars; ++i)
            {
@@ -1641,7 +1698,7 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id
                  if (i == j)
                    continue;
 
-                 tab_double (t, heading_columns + i,  y + j, 0, significance_of_correlation (rho, w), NULL);
+                 tab_double (t, heading_columns + i,  y + j, 0, significance_of_correlation (rho, w), NULL, RC_PVALUE);
                }
            }
        }
@@ -1649,22 +1706,9 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id
 
   if (factor->print & PRINT_DETERMINANT)
     {
-      int sign = 0;
-      double det = 0.0;
-
-      const int size = idata->corr->size1;
-      gsl_permutation *p = gsl_permutation_calloc (size);
-      gsl_matrix *tmp = gsl_matrix_calloc (size, size);
-      gsl_matrix_memcpy (tmp, idata->corr);
-
-      gsl_linalg_LU_decomp (tmp, p, &sign);
-      det = gsl_linalg_LU_det (tmp, sign);
-      gsl_permutation_free (p);
-      gsl_matrix_free (tmp);
-
-
       tab_text (t, 0, nr, TAB_LEFT | TAT_TITLE, _("Determinant"));
-      tab_double (t, 1, nr, 0, det, NULL);
+
+      tab_double (t, 1, nr, 0, idata->detR, NULL, RC_OTHER);
     }
 
   tab_submit (t);
@@ -1682,7 +1726,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
   const gsl_matrix *analysis_matrix;
   struct idata *idata = idata_alloc (factor->n_vars);
 
-  struct covariance *cov = covariance_create (factor->n_vars, factor->vars,
+  struct covariance *cov = covariance_1pass_create (factor->n_vars, factor->vars,
                                              factor->wv, factor->exclude);
 
   for ( ; (c = casereader_read (r) ); case_unref (c))
@@ -1692,24 +1736,49 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
   idata->cov = covariance_calculate (cov);
 
+  if (idata->cov == NULL)
+    {
+      msg (MW, _("The dataset contains no complete observations. No analysis will be performed."));
+      covariance_destroy (cov);
+      goto finish;
+    }
+
   var_matrix = covariance_moments (cov, MOMENT_VARIANCE);
   mean_matrix = covariance_moments (cov, MOMENT_MEAN);
   idata->n = covariance_moments (cov, MOMENT_NONE);
+  
 
   if ( factor->method == METHOD_CORR)
     {
       idata->corr = correlation_from_covariance (idata->cov, var_matrix);
+      
       analysis_matrix = idata->corr;
     }
   else
     analysis_matrix = idata->cov;
 
+
+  if (factor->print & PRINT_DETERMINANT
+      || factor->print & PRINT_KMO)
+    {
+      int sign = 0;
+
+      const int size = idata->corr->size1;
+      gsl_permutation *p = gsl_permutation_calloc (size);
+      gsl_matrix *tmp = gsl_matrix_calloc (size, size);
+      gsl_matrix_memcpy (tmp, idata->corr);
+
+      gsl_linalg_LU_decomp (tmp, p, &sign);
+      idata->detR = gsl_linalg_LU_det (tmp, sign);
+      gsl_permutation_free (p);
+      gsl_matrix_free (tmp);
+    }
+
   if ( factor->print & PRINT_UNIVARIATE)
     {
+      const struct fmt_spec *wfmt = factor->wv ? var_get_print_format (factor->wv) : & F_8_0;
       const int nc = 4;
       int i;
-      const struct fmt_spec *wfmt = factor->wv ? var_get_print_format (factor->wv) : & F_8_0;
-
 
       const int heading_columns = 1;
       const int heading_rows = 1;
@@ -1717,6 +1786,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
       const int nr = heading_rows + factor->n_vars;
 
       struct tab_table *t = tab_create (nc, nr);
+      tab_set_format (t, RC_WEIGHT, wfmt);
       tab_title (t, _("Descriptive Statistics"));
 
       tab_headers (t, heading_columns, 0, heading_rows, 0);
@@ -1747,39 +1817,122 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
          const struct variable *v = factor->vars[i];
          tab_text (t, 0, i + heading_rows, TAB_LEFT | TAT_TITLE, var_to_string (v));
 
-         tab_double (t, 1, i + heading_rows, 0, gsl_matrix_get (mean_matrix, i, i), NULL);
-         tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (var_matrix, i, i)), NULL);
-         tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->n, i, i), wfmt);
+         tab_double (t, 1, i + heading_rows, 0, gsl_matrix_get (mean_matrix, i, i), NULL, RC_OTHER);
+         tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (var_matrix, i, i)), NULL, RC_OTHER);
+         tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->n, i, i), NULL, RC_WEIGHT);
        }
 
       tab_submit (t);
     }
 
+  if (factor->print & PRINT_KMO)
+    {
+      int i;
+      double sum_ssq_r = 0;
+      double sum_ssq_a = 0;
+
+      double df = factor->n_vars * ( factor->n_vars - 1) / 2;
+
+      double w = 0;
+
+
+      double xsq;
+
+      const int heading_columns = 2;
+      const int heading_rows = 0;
+
+      const int nr = heading_rows + 4;
+      const int nc = heading_columns + 1;
+
+      gsl_matrix *a, *x;
+
+      struct tab_table *t = tab_create (nc, nr);
+      tab_title (t, _("KMO and Bartlett's Test"));
+
+      x  = clone_matrix (idata->corr);
+      gsl_linalg_cholesky_decomp (x);
+      gsl_linalg_cholesky_invert (x);
+
+      a = anti_image (x);
+
+      for (i = 0; i < x->size1; ++i)
+       {
+         sum_ssq_r += ssq_od_n (x, i);
+         sum_ssq_a += ssq_od_n (a, i);
+       }
+
+      gsl_matrix_free (a);
+      gsl_matrix_free (x);
+
+      tab_headers (t, heading_columns, 0, heading_rows, 0);
+
+      /* Outline the box */
+      tab_box (t,
+              TAL_2, TAL_2,
+              -1, -1,
+              0, 0,
+              nc - 1, nr - 1);
+
+      tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
+
+      tab_text (t, 0, 0, TAT_TITLE | TAB_LEFT, _("Kaiser-Meyer-Olkin Measure of Sampling Adequacy"));
+
+      tab_double (t, 2, 0, 0, sum_ssq_r /  (sum_ssq_r + sum_ssq_a), NULL, RC_OTHER);
+
+      tab_text (t, 0, 1, TAT_TITLE | TAB_LEFT, _("Bartlett's Test of Sphericity"));
+
+      tab_text (t, 1, 1, TAT_TITLE, _("Approx. Chi-Square"));
+      tab_text (t, 1, 2, TAT_TITLE, _("df"));
+      tab_text (t, 1, 3, TAT_TITLE, _("Sig."));
+
+
+      /* The literature doesn't say what to do for the value of W when 
+        missing values are involved.  The best thing I can think of
+        is to take the mean average. */
+      w = 0;
+      for (i = 0; i < idata->n->size1; ++i)
+       w += gsl_matrix_get (idata->n, i, i);
+      w /= idata->n->size1;
+
+      xsq = w - 1 - (2 * factor->n_vars + 5) / 6.0;
+      xsq *= -log (idata->detR);
+
+      tab_double (t, 2, 1, 0, xsq, NULL, RC_OTHER);
+      tab_double (t, 2, 2, 0, df, NULL, RC_INTEGER);
+      tab_double (t, 2, 3, 0, gsl_cdf_chisq_Q (xsq, df), NULL, RC_PVALUE);
+      
+
+      tab_submit (t);
+    }
+
   show_correlation_matrix (factor, idata);
+  covariance_destroy (cov);
 
-#if 1
   {
+    gsl_matrix *am = matrix_dup (analysis_matrix);
     gsl_eigen_symmv_workspace *workspace = gsl_eigen_symmv_alloc (factor->n_vars);
     
-    gsl_eigen_symmv (matrix_dup (analysis_matrix), idata->eval, idata->evec, workspace);
+    gsl_eigen_symmv (am, idata->eval, idata->evec, workspace);
 
     gsl_eigen_symmv_free (workspace);
+    gsl_matrix_free (am);
   }
 
   gsl_eigen_symmv_sort (idata->eval, idata->evec, GSL_EIGEN_SORT_ABS_DESC);
-#endif
 
   idata->n_extractions = n_extracted_factors (factor, idata);
 
   if (idata->n_extractions == 0)
     {
-      msg (MW, _("The FACTOR criteria result in zero factors extracted. Therefore no analysis will be performed."));
+      msg (MW, _("The %s criteria result in zero factors extracted. Therefore no analysis will be performed."), "FACTOR");
       goto finish;
     }
 
   if (idata->n_extractions > factor->n_vars)
     {
-      msg (MW, _("The FACTOR criteria result in more factors than variables, which is not meaningful. No analysis will be performed."));
+      msg (MW, 
+          _("The %s criteria result in more factors than variables, which is not meaningful. No analysis will be performed."), 
+          "FACTOR");
       goto finish;
     }
     
@@ -1809,7 +1962,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
        gsl_vector_memcpy (initial_communalities, idata->msr);
 
-       for (i = 0; i < factor->iterations; ++i)
+       for (i = 0; i < factor->extraction_iterations; ++i)
          {
            double min, max;
            gsl_vector_memcpy (diff, idata->msr);
@@ -1876,6 +2029,8 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
 
 
+    gsl_matrix_free (factor_matrix);
+    gsl_vector_free (rotated_loadings);
     gsl_vector_free (initial_communalities);
     gsl_vector_free (extracted_communalities);
   }