FACTOR: Add Anti-image matrix output
[pspp] / src / language / stats / factor.c
index 472c340fd7feab0feab6c2f89ecc22b0d148c80f..a1cd333389c7276bc5a6ed7c37b776f2daf583a1 100644 (file)
@@ -1,5 +1,6 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015, 2016 Free Software Foundation, Inc.
+   Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015,
+   2016, 2017 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
@@ -24,6 +25,7 @@
 #include <gsl/gsl_sort_vector.h>
 #include <gsl/gsl_cdf.h>
 
+#include "data/any-reader.h"
 #include "data/casegrouper.h"
 #include "data/casereader.h"
 #include "data/casewriter.h"
@@ -35,6 +37,8 @@
 #include "language/lexer/lexer.h"
 #include "language/lexer/value-parser.h"
 #include "language/lexer/variable-parser.h"
+#include "language/data-io/file-handle.h"
+#include "language/data-io/matrix-reader.h"
 #include "libpspp/cast.h"
 #include "libpspp/message.h"
 #include "libpspp/misc.h"
@@ -45,6 +49,7 @@
 #include "output/charts/scree.h"
 #include "output/tab.h"
 
+
 #include "gettext.h"
 #define _(msgid) gettext (msgid)
 #define N_(msgid) msgid
@@ -221,13 +226,13 @@ struct cmd_factor
   bool sort;
 };
 
+
 struct idata
 {
   /* Intermediate values used in calculation */
+  struct matrix_material mm;
 
-  const gsl_matrix *corr ;  /* The correlation matrix */
-  gsl_matrix *cov ;         /* The covariance matrix */
-  const gsl_matrix *n ;     /* Matrix of number of samples */
+  gsl_matrix *analysis_matrix; /* A pointer to either mm.corr or mm.cov */
 
   gsl_vector *eval ;  /* The eigenvalues */
   gsl_matrix *evec ;  /* The eigenvectors */
@@ -237,6 +242,10 @@ struct idata
   gsl_vector *msr ;  /* Multiple Squared Regressions */
 
   double detR;  /* The determinant of the correlation matrix */
+
+  gsl_matrix *ai_cov; /* The anti-image covariance matrix */
+  gsl_matrix *ai_cor; /* The anti-image correlation matrix */
+  struct covariance *cvm;
 };
 
 static struct idata *
@@ -259,39 +268,31 @@ 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));
+  gsl_matrix_free (id->ai_cov);
+  gsl_matrix_free (id->ai_cor);
 
   free (id);
 }
 
-
-static gsl_matrix *
-anti_image (const gsl_matrix *m)
+/* Return the sum of squares of all the elements in row J excluding column J */
+static double
+ssq_row_od_n (const gsl_matrix *m, int j)
 {
-  int i, j;
-  gsl_matrix *a;
+  int i;
+  double ss = 0;
   assert (m->size1 == m->size2);
 
-  a = gsl_matrix_alloc (m->size1, m->size2);
+  assert (j < m->size1);
 
   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);
-       }
+      if (i == j ) continue;
+      ss += pow2 (gsl_matrix_get (m, i, j));
     }
 
-  return a;
+  return ss;
 }
 
-
 /* Return the sum of all the elements excluding row N */
 static double
 ssq_od_n (const gsl_matrix *m, int n)
@@ -315,6 +316,58 @@ ssq_od_n (const gsl_matrix *m, int n)
 }
 
 
+static gsl_matrix *
+anti_image_corr (const gsl_matrix *m, const struct idata *idata)
+{
+  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 /= sqrt (gsl_matrix_get (m, i, i) *
+                      gsl_matrix_get (m, j, j));
+        }
+    }
+
+  for (i = 0; i < m->size1; ++i)
+    {
+      double r = ssq_row_od_n (idata->mm.corr, i);
+      double u = ssq_row_od_n (a, i);
+      gsl_matrix_set (a, i, i, r / (r + u));
+    }
+
+  return a;
+}
+
+static gsl_matrix *
+anti_image_cov (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;
+}
 
 #if 0
 static void
@@ -1002,11 +1055,14 @@ iterate_factor_matrix (const gsl_matrix *r, gsl_vector *communalities, gsl_matri
 
 static bool run_factor (struct dataset *ds, const struct cmd_factor *factor);
 
+static void do_factor_by_matrix (const struct cmd_factor *factor, struct idata *idata);
+
+
 
 int
 cmd_factor (struct lexer *lexer, struct dataset *ds)
 {
-  const struct dictionary *dict = dataset_dict (ds);
+  struct dictionary *dict = NULL;
   int n_iterations = 25;
   struct cmd_factor factor;
   factor.n_vars = 0;
@@ -1026,26 +1082,77 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
   factor.sort = false;
   factor.plot = 0;
   factor.rotation = ROT_VARIMAX;
+  factor.wv = NULL;
 
   factor.rconverge = 0.0001;
 
-  factor.wv = dict_get_weight (dict);
-
   lex_match (lexer, T_SLASH);
 
-  if (!lex_force_match_id (lexer, "VARIABLES"))
+  struct matrix_reader *mr = NULL;
+  struct casereader *matrix_reader = NULL;
+
+  if (lex_match_id (lexer, "VARIABLES"))
     {
-      goto error;
+      lex_match (lexer, T_EQUALS);
+      dict = dataset_dict (ds);
+      factor.wv = dict_get_weight (dict);
+
+      if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars,
+                                 PV_NO_DUPLICATE | PV_NUMERIC))
+       goto error;
     }
+  else if (lex_match_id (lexer, "MATRIX"))
+    {
+      lex_match (lexer, T_EQUALS);
+      if (! lex_force_match_id (lexer, "IN"))
+       goto error;
+      if (!lex_force_match (lexer, T_LPAREN))
+       {
+         goto error;
+       }
+      if (lex_match_id (lexer, "CORR"))
+       {
+       }
+      else if (lex_match_id (lexer, "COV"))
+       {
+       }
+      else
+       {
+         lex_error (lexer, _("Matrix input for %s must be either COV or CORR"), "FACTOR");
+         goto error;
+       }
+      if (! lex_force_match (lexer, T_EQUALS))
+       goto error;
+      if (lex_match (lexer, T_ASTERISK))
+       {
+         dict = dataset_dict (ds);
+         matrix_reader = casereader_clone (dataset_source (ds));
+       }
+      else
+       {
+         struct file_handle *fh = fh_parse (lexer, FH_REF_FILE, NULL);
+         if (fh == NULL)
+           goto error;
 
-  lex_match (lexer, T_EQUALS);
+         matrix_reader
+           = any_reader_open_and_decode (fh, NULL, &dict, NULL);
 
-  if (!parse_variables_const (lexer, dict, &factor.vars, &factor.n_vars,
-                             PV_NO_DUPLICATE | PV_NUMERIC))
-    goto error;
+         if (! (matrix_reader && dict))
+           {
+             goto error;
+           }
+       }
 
-  if (factor.n_vars < 2)
-    msg (MW, _("Factor analysis on a single variable is not useful."));
+      if (! lex_force_match (lexer, T_RPAREN))
+       goto error;
+
+      mr = create_matrix_reader_from_case_reader (dict, matrix_reader,
+                                                 &factor.vars, &factor.n_vars);
+    }
+  else
+    {
+      goto error;
+    }
 
   while (lex_token (lexer) != T_ENDCMD)
     {
@@ -1308,10 +1415,11 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
              else if (lex_match_id (lexer, "INV"))
                {
                }
+#endif
              else if (lex_match_id (lexer, "AIC"))
                {
+                 factor.print |= PRINT_AIC;
                }
-#endif
              else if (lex_match_id (lexer, "SIG"))
                {
                  factor.print |= PRINT_SIG;
@@ -1320,11 +1428,10 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
                {
                  factor.print |= PRINT_CORRELATION;
                }
-#if FACTOR_FULLY_IMPLEMENTED
              else if (lex_match_id (lexer, "COVARIANCE"))
                {
+                 factor.print |= PRINT_COVARIANCE;
                }
-#endif
              else if (lex_match_id (lexer, "ROTATION"))
                {
                  factor.print |= PRINT_ROTATION;
@@ -1408,13 +1515,43 @@ cmd_factor (struct lexer *lexer, struct dataset *ds)
   if ( factor.rotation == ROT_NONE )
     factor.print &= ~PRINT_ROTATION;
 
-  if ( ! run_factor (ds, &factor))
-    goto error;
+  if (factor.n_vars < 2)
+    msg (MW, _("Factor analysis on a single variable is not useful."));
+
+  if (factor.n_vars < 1)
+    {
+      msg (ME, _("Factor analysis without variables is not possible."));
+      goto error;
+    }
+
+  if (matrix_reader)
+    {
+      struct idata *id = idata_alloc (factor.n_vars);
+
+      while (next_matrix_from_reader (&id->mm, mr,
+                                     factor.vars, factor.n_vars))
+       {
+         do_factor_by_matrix (&factor, id);
+
+         gsl_matrix_free (id->mm.corr);
+         id->mm.corr = NULL;
+         gsl_matrix_free (id->mm.cov);
+         id->mm.cov = NULL;
+       }
+
+      idata_free (id);
+    }
+  else
+    if ( ! run_factor (ds, &factor))
+      goto error;
 
+
+  destroy_matrix_reader (mr);
   free (factor.vars);
   return CMD_SUCCESS;
 
  error:
+  destroy_matrix_reader (mr);
   free (factor.vars);
   return CMD_FAILURE;
 }
@@ -1474,14 +1611,14 @@ the_communality (const gsl_matrix *evec, const gsl_vector *eval, int n, int n_fa
 
 /* Return the communality of variable N, calculated to N_FACTORS */
 static double
-communality (struct idata *idata, int n, int n_factors)
+communality (const struct idata *idata, int n, int n_factors)
 {
   return the_communality (idata->evec, idata->eval, n, n_factors);
 }
 
 
 static void
-show_scree (const struct cmd_factor *f, struct idata *idata)
+show_scree (const struct cmd_factor *f, const struct idata *idata)
 {
   struct scree *s;
   const char *label ;
@@ -1566,7 +1703,7 @@ show_communalities (const struct cmd_factor * factor,
 
 
 static void
-show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const char *title, const gsl_matrix *fm)
+show_factor_matrix (const struct cmd_factor *factor, const struct idata *idata, const char *title, const gsl_matrix *fm)
 {
   int i;
 
@@ -1659,7 +1796,8 @@ show_factor_matrix (const struct cmd_factor *factor, struct idata *idata, const
 
 
 static void
-show_explained_variance (const struct cmd_factor * factor, struct idata *idata,
+show_explained_variance (const struct cmd_factor * factor,
+                        const struct idata *idata,
                         const gsl_vector *initial_eigenvalues,
                         const gsl_vector *extracted_eigenvalues,
                         const gsl_vector *rotated_loadings)
@@ -1887,13 +2025,88 @@ show_factor_correlation (const struct cmd_factor * factor, const gsl_matrix *fcm
   for (i = 0 ; i < fcm->size1; ++i)
     {
       for (j = 0 ; j < fcm->size2; ++j)
-       tab_double (t, heading_columns + i,  heading_rows +j, 0,
+       tab_double (t, heading_columns + j,  heading_rows + i, 0,
                    gsl_matrix_get (fcm, i, j), NULL, RC_OTHER);
     }
 
   tab_submit (t);
 }
 
+static void
+show_aic (const struct cmd_factor *factor, const struct idata *idata)
+{
+  struct tab_table *t ;
+  size_t i;
+
+  const int heading_rows = 1;
+  const int heading_columns = 2;
+
+  const int nc = heading_columns + factor->n_vars;
+  const int nr = heading_rows + 2 * factor->n_vars;
+
+  if ((factor->print & PRINT_AIC) == 0)
+    return;
+
+  t = tab_create (nc, nr);
+
+  tab_title (t, _("Anti-Image Matrices"));
+
+  tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
+
+  tab_headers (t, heading_columns, 0, heading_rows, 0);
+
+  tab_vline (t, TAL_2, 2, 0, nr - 1);
+
+  /* Outline the box */
+  tab_box (t,
+          TAL_2, TAL_2,
+          -1, -1,
+          0, 0,
+          nc - 1, nr - 1);
+
+  /* Vertical lines */
+  tab_box (t,
+          -1, -1,
+          -1, TAL_1,
+          heading_columns, 0,
+          nc - 1, nr - 1);
+
+
+  for (i = 0; i < factor->n_vars; ++i)
+    tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i]));
+
+  tab_text (t, 0, heading_rows, TAT_TITLE, _("Anti-image Covariance"));
+  tab_hline (t, TAL_1, 0, nc - 1, heading_rows + factor->n_vars);
+  tab_text (t, 0, heading_rows + factor->n_vars, TAT_TITLE, _("Anti-image Correlation"));
+
+  for (i = 0; i < factor->n_vars; ++i)
+    {
+      tab_text (t, 1, i + heading_rows, TAT_TITLE,
+               var_to_string (factor->vars[i]));
+
+      tab_text (t, 1, factor->n_vars + i + heading_rows, TAT_TITLE,
+               var_to_string (factor->vars[i]));
+    }
+
+  for (i = 0; i < factor->n_vars; ++i)
+    {
+      int j;
+      for (j = 0; j < factor->n_vars; ++j)
+       {
+         tab_double (t, heading_columns + i, heading_rows + j, 0,
+                     gsl_matrix_get (idata->ai_cov, i, j), NULL, RC_OTHER);
+       }
+
+
+      for (j = 0; j < factor->n_vars; ++j)
+       {
+         tab_double (t, heading_columns + i, factor->n_vars + heading_rows + j, 0,
+                     gsl_matrix_get (idata->ai_cor, i, j), NULL, RC_OTHER);
+       }
+    }
+
+  tab_submit (t);
+}
 
 static void
 show_correlation_matrix (const struct cmd_factor *factor, const struct idata *idata)
@@ -1983,7 +2196,7 @@ 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, RC_OTHER);
+               tab_double (t, heading_columns + j,  y + i, 0, gsl_matrix_get (idata->mm.corr, i, j), NULL, RC_OTHER);
            }
        }
 
@@ -1996,13 +2209,13 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id
            {
              for (j = 0; j < factor->n_vars; ++j)
                {
-                 double rho = gsl_matrix_get (idata->corr, i, j);
-                 double w = gsl_matrix_get (idata->n, i, j);
+                 double rho = gsl_matrix_get (idata->mm.corr, i, j);
+                 double w = gsl_matrix_get (idata->mm.n, i, j);
 
                  if (i == j)
                    continue;
 
-                 tab_double (t, heading_columns + i,  y + j, 0, significance_of_correlation (rho, w), NULL, RC_PVALUE);
+                 tab_double (t, heading_columns + j,  y + i, 0, significance_of_correlation (rho, w), NULL, RC_PVALUE);
                }
            }
        }
@@ -2018,59 +2231,174 @@ show_correlation_matrix (const struct cmd_factor *factor, const struct idata *id
   tab_submit (t);
 }
 
+static void
+show_covariance_matrix (const struct cmd_factor *factor, const struct idata *idata)
+{
+  struct tab_table *t ;
+  size_t i, j;
+  int y_pos_corr = -1;
+  int suffix_rows = 0;
+
+  const int heading_rows = 1;
+  const int heading_columns = 1;
+
+  int nc = heading_columns ;
+  int nr = heading_rows ;
+  int n_data_sets = 0;
+
+  if (factor->print & PRINT_COVARIANCE)
+    {
+      y_pos_corr = n_data_sets;
+      n_data_sets++;
+      nc = heading_columns + factor->n_vars;
+    }
+
+  nr += n_data_sets * factor->n_vars;
+
+  /* If the table would contain only headings, don't bother rendering it */
+  if (nr <= heading_rows && suffix_rows == 0)
+    return;
+
+  t = tab_create (nc, nr + suffix_rows);
+
+  tab_title (t, _("Covariance Matrix"));
+
+  tab_hline (t, TAL_1, 0, nc - 1, heading_rows);
+
+  if (nr > heading_rows)
+    {
+      tab_headers (t, heading_columns, 0, heading_rows, 0);
+
+      tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
+
+      /* Outline the box */
+      tab_box (t,
+              TAL_2, TAL_2,
+              -1, -1,
+              0, 0,
+              nc - 1, nr - 1);
+
+      /* Vertical lines */
+      tab_box (t,
+              -1, -1,
+              -1, TAL_1,
+              heading_columns, 0,
+              nc - 1, nr - 1);
+
+
+      for (i = 0; i < factor->n_vars; ++i)
+       tab_text (t, heading_columns + i, 0, TAT_TITLE, var_to_string (factor->vars[i]));
+
+
+      for (i = 0 ; i < n_data_sets; ++i)
+       {
+         int y = heading_rows + i * factor->n_vars;
+         size_t v;
+         for (v = 0; v < factor->n_vars; ++v)
+           tab_text (t, heading_columns -1, y + v, TAT_TITLE, var_to_string (factor->vars[v]));
+
+         tab_hline (t, TAL_1, 0, nc - 1, y);
+       }
+
+      if (factor->print & PRINT_COVARIANCE)
+       {
+         const double y = heading_rows + y_pos_corr;
+
+         for (i = 0; i < factor->n_vars; ++i)
+           {
+             for (j = 0; j < factor->n_vars; ++j)
+               tab_double (t, heading_columns + j,  y + i, 0, gsl_matrix_get (idata->mm.cov, i, j), NULL, RC_OTHER);
+           }
+       }
+    }
+
+  tab_submit (t);
+}
 
 
 static void
 do_factor (const struct cmd_factor *factor, struct casereader *r)
 {
   struct ccase *c;
-  const gsl_matrix *var_matrix;
-  const gsl_matrix *mean_matrix;
-
-  const gsl_matrix *analysis_matrix;
   struct idata *idata = idata_alloc (factor->n_vars);
 
-  struct covariance *cov = covariance_1pass_create (factor->n_vars, factor->vars,
-                                             factor->wv, factor->exclude);
+  idata->cvm = covariance_1pass_create (factor->n_vars, factor->vars,
+                                       factor->wv, factor->exclude, true);
 
   for ( ; (c = casereader_read (r) ); case_unref (c))
     {
-      covariance_accumulate (cov, c);
+      covariance_accumulate (idata->cvm, c);
     }
 
-  idata->cov = covariance_calculate (cov);
+  idata->mm.cov = covariance_calculate (idata->cvm);
 
-  if (idata->cov == NULL)
+  if (idata->mm.cov == NULL)
     {
       msg (MW, _("The dataset contains no complete observations. No analysis will be performed."));
-      covariance_destroy (cov);
+      covariance_destroy (idata->cvm);
       goto finish;
     }
 
-  var_matrix = covariance_moments (cov, MOMENT_VARIANCE);
-  mean_matrix = covariance_moments (cov, MOMENT_MEAN);
-  idata->n = covariance_moments (cov, MOMENT_NONE);
+  idata->mm.var_matrix = covariance_moments (idata->cvm, MOMENT_VARIANCE);
+  idata->mm.mean_matrix = covariance_moments (idata->cvm, MOMENT_MEAN);
+  idata->mm.n = covariance_moments (idata->cvm, MOMENT_NONE);
 
+  do_factor_by_matrix (factor, idata);
 
-  if ( factor->method == METHOD_CORR)
-    {
-      idata->corr = correlation_from_covariance (idata->cov, var_matrix);
+ finish:
+  gsl_matrix_free (idata->mm.corr);
+  gsl_matrix_free (idata->mm.cov);
+
+  idata_free (idata);
+  casereader_destroy (r);
+}
 
-      analysis_matrix = idata->corr;
+static void
+do_factor_by_matrix (const struct cmd_factor *factor, struct idata *idata)
+{
+  if (!idata->mm.cov && !idata->mm.corr)
+    {
+      msg (ME, _("The dataset has no complete covariance or correlation matrix."));
+      return;
     }
+
+  if (idata->mm.cov && !idata->mm.corr)
+    idata->mm.corr = correlation_from_covariance (idata->mm.cov, idata->mm.var_matrix);
+  if (idata->mm.corr && !idata->mm.cov)
+    idata->mm.cov = covariance_from_correlation (idata->mm.corr, idata->mm.var_matrix);
+  if (factor->method == METHOD_CORR)
+    idata->analysis_matrix = idata->mm.corr;
   else
-    analysis_matrix = idata->cov;
+    idata->analysis_matrix = idata->mm.cov;
+
+  gsl_matrix *r_inv;
+  r_inv  = clone_matrix (idata->mm.corr);
+  gsl_linalg_cholesky_decomp (r_inv);
+  gsl_linalg_cholesky_invert (r_inv);
 
+  idata->ai_cov = anti_image_cov (r_inv);
+  idata->ai_cor = anti_image_corr (r_inv, idata);
+
+  int i;
+  double sum_ssq_r = 0;
+  double sum_ssq_a = 0;
+  for (i = 0; i < r_inv->size1; ++i)
+    {
+      sum_ssq_r += ssq_od_n (r_inv, i);
+      sum_ssq_a += ssq_od_n (idata->ai_cov, i);
+    }
+
+  gsl_matrix_free (r_inv);
 
   if (factor->print & PRINT_DETERMINANT
       || factor->print & PRINT_KMO)
     {
       int sign = 0;
 
-      const int size = idata->corr->size1;
+      const int size = idata->mm.corr->size1;
       gsl_permutation *p = gsl_permutation_calloc (size);
       gsl_matrix *tmp = gsl_matrix_calloc (size, size);
-      gsl_matrix_memcpy (tmp, idata->corr);
+      gsl_matrix_memcpy (tmp, idata->mm.corr);
 
       gsl_linalg_LU_decomp (tmp, p, &sign);
       idata->detR = gsl_linalg_LU_det (tmp, sign);
@@ -2121,9 +2449,9 @@ 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, 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_double (t, 1, i + heading_rows, 0, gsl_matrix_get (idata->mm.mean_matrix, i, i), NULL, RC_OTHER);
+         tab_double (t, 2, i + heading_rows, 0, sqrt (gsl_matrix_get (idata->mm.var_matrix, i, i)), NULL, RC_OTHER);
+         tab_double (t, 3, i + heading_rows, 0, gsl_matrix_get (idata->mm.n, i, i), NULL, RC_WEIGHT);
        }
 
       tab_submit (t);
@@ -2132,10 +2460,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
   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 df = factor->n_vars * (factor->n_vars - 1) / 2;
 
       double w = 0;
 
@@ -2148,25 +2473,11 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
       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);
 
@@ -2194,9 +2505,9 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
         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;
+      for (i = 0; i < idata->mm.n->size1; ++i)
+       w += gsl_matrix_get (idata->mm.n, i, i);
+      w /= idata->mm.n->size1;
 
       xsq = w - 1 - (2 * factor->n_vars + 5) / 6.0;
       xsq *= -log (idata->detR);
@@ -2210,10 +2521,12 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
     }
 
   show_correlation_matrix (factor, idata);
-  covariance_destroy (cov);
+  show_covariance_matrix (factor, idata);
+  if (idata->cvm)
+    covariance_destroy (idata->cvm);
 
   {
-    gsl_matrix *am = matrix_dup (analysis_matrix);
+    gsl_matrix *am = matrix_dup (idata->analysis_matrix);
     gsl_eigen_symmv_workspace *workspace = gsl_eigen_symmv_alloc (factor->n_vars);
 
     gsl_eigen_symmv (am, idata->eval, idata->evec, workspace);
@@ -2229,7 +2542,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
   if (idata->n_extractions == 0)
     {
       msg (MW, _("The %s criteria result in zero factors extracted. Therefore no analysis will be performed."), "FACTOR");
-      goto finish;
+      return;
     }
 
   if (idata->n_extractions > factor->n_vars)
@@ -2237,7 +2550,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
       msg (MW,
           _("The %s criteria result in more factors than variables, which is not meaningful. No analysis will be performed."),
           "FACTOR");
-      goto finish;
+      return;
     }
 
   {
@@ -2256,11 +2569,11 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
     if ( factor->extraction == EXTRACTION_PAF)
       {
        gsl_vector *diff = gsl_vector_alloc (idata->msr->size);
-       struct smr_workspace *ws = ws_create (analysis_matrix);
+       struct smr_workspace *ws = ws_create (idata->analysis_matrix);
 
        for (i = 0 ; i < factor->n_vars ; ++i)
          {
-           double r2 = squared_multiple_correlation (analysis_matrix, i, ws);
+           double r2 = squared_multiple_correlation (idata->analysis_matrix, i, ws);
 
            gsl_vector_set (idata->msr, i, r2);
          }
@@ -2273,7 +2586,7 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
            double min, max;
            gsl_vector_memcpy (diff, idata->msr);
 
-           iterate_factor_matrix (analysis_matrix, idata->msr, factor_matrix, fmw);
+           iterate_factor_matrix (idata->analysis_matrix, idata->msr, factor_matrix, fmw);
 
            gsl_vector_sub (diff, idata->msr);
 
@@ -2296,16 +2609,16 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
 
        gsl_vector_memcpy (extracted_communalities, initial_communalities);
 
-       iterate_factor_matrix (analysis_matrix, extracted_communalities, factor_matrix, fmw);
+       iterate_factor_matrix (idata->analysis_matrix, extracted_communalities, factor_matrix, fmw);
 
 
        extracted_eigenvalues = idata->eval;
       }
 
 
+    show_aic (factor, idata);
     show_communalities (factor, initial_communalities, extracted_communalities);
 
-
     if ( factor->rotation != ROT_NONE)
       {
        rotated_factors = gsl_matrix_calloc (factor_matrix->size1, factor_matrix->size2);
@@ -2340,7 +2653,8 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
       {
        show_factor_matrix (factor, idata,
                            (factor->rotation == ROT_PROMAX) ? _("Structure Matrix") :
-                           (factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") : _("Rotated Factor Matrix")),
+                           (factor->extraction == EXTRACTION_PC ? _("Rotated Component Matrix") :
+                            _("Rotated Factor Matrix")),
                            rotated_factors);
 
        gsl_matrix_free (rotated_factors);
@@ -2357,12 +2671,6 @@ do_factor (const struct cmd_factor *factor, struct casereader *r)
     gsl_vector_free (initial_communalities);
     gsl_vector_free (extracted_communalities);
   }
-
- finish:
-
-  idata_free (idata);
-
-  casereader_destroy (r);
 }