From: Jason H Stover <jhs@math.gcsu.edu>
Date: Thu, 6 Nov 2008 15:45:40 +0000 (-0500)
Subject: covariance-matrix.c: Added one-pass algorithm to compute covariance
X-Git-Tag: v0.7.1~50^2~18^2
X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?p=pspp-builds.git;a=commitdiff_plain;h=93cf6259300c4ddeef1df6e40b3958bb9e6802f5

covariance-matrix.c: Added one-pass algorithm to compute covariance
matrices.  New structure covariance_matrix contains the covariance
matrix, related data and member functions. Added different routines to
handle pairwise or listwise deletion of missing values.  Added many
new functions.

glm.q: Changed run_glm to use the new interface to the covariance
matrix. Removed unused function data_pass_one().
---

diff --git a/src/language/stats/glm.q b/src/language/stats/glm.q
index 024b5429..d6b4952c 100644
--- a/src/language/stats/glm.q
+++ b/src/language/stats/glm.q
@@ -158,69 +158,16 @@ coeff_init (pspp_linreg_cache * c, const struct design_matrix *cov)
   pspp_coeff_init (c->coeff, cov);
 }
 
-/* Encode categorical variables.
-   Returns number of valid cases. */
-static int
-data_pass_one (struct casereader *input,
-	       const struct variable **vars, size_t n_vars,
-	       struct moments_var **mom)
-{
-  int n_data;
-  struct ccase c;
-  size_t i;
-
-  for (i = 0; i < n_vars; i++)
-    {
-      mom[i] = xmalloc (sizeof (*mom[i]));
-      mom[i]->v = vars[i];
-      mom[i]->mean = xmalloc (sizeof (*mom[i]->mean));
-      mom[i]->variance = xmalloc (sizeof (*mom[i]->mean));
-      mom[i]->weight = xmalloc (sizeof (*mom[i]->weight));
-      mom[i]->m = moments1_create (MOMENT_VARIANCE);
-      if (var_is_alpha (vars[i]))
-	cat_stored_values_create (vars[i]);
-    }
-
-  n_data = 0;
-  for (; casereader_read (input, &c); case_destroy (&c))
-    {
-      /*
-         The second condition ensures the program will run even if
-         there is only one variable to act as both explanatory and
-         response.
-       */
-      for (i = 0; i < n_vars; i++)
-	{
-	  const union value *val = case_data (&c, vars[i]);
-	  if (var_is_alpha (vars[i]))
-	    cat_value_update (vars[i], val);
-	  else
-	    moments1_add (mom[i]->m, val->f, 1.0);
-	}
-      n_data++;
-    }
-  casereader_destroy (input);
-  for (i = 0; i < n_vars; i++)
-    {
-      if (var_is_numeric (mom[i]->v))
-	{
-	  moments1_calculate (mom[i]->m, mom[i]->weight, mom[i]->mean,
-			      mom[i]->variance, NULL, NULL);
-	}
-    }
-
-  return n_data;
-}
 
 static pspp_linreg_cache *
-fit_model (const struct design_matrix *cov, const struct moments1 **mom, 
+fit_model (const struct covariance_matrix *cov,
 	   const struct variable *dep_var, 
 	   const struct variable ** indep_vars, 
 	   size_t n_data, size_t n_indep)
 {
   pspp_linreg_cache *result = NULL;
   result = pspp_linreg_cache_alloc (dep_var, indep_vars, n_data, n_indep);
-  coeff_init (result, cov);
+  coeff_init (result, covariance_to_design (cov));
   pspp_linreg_with_cov (cov, result);  
   
   return result;
@@ -242,9 +189,7 @@ run_glm (struct casereader *input,
   size_t n_all_vars;
   size_t n_data;		/* Number of valid cases. */
   struct casereader *reader;
-  struct design_matrix *cov;
-  struct hsh_table *cov_hash;
-  struct moments1 **mom;
+  struct covariance_matrix *cov;
 
   if (!casereader_peek (input, 0, &c))
     {
@@ -277,9 +222,6 @@ run_glm (struct casereader *input,
       all_vars[i + n_dependent] = cmd->v_by[i];
     }
   n_indep = cmd->n_by;
-  mom = xnmalloc (n_all_vars, sizeof (*mom));
-  for (i = 0; i < n_all_vars; i++)
-    mom[i] = moments1_create (MOMENT_MEAN);
 
   reader = casereader_clone (input);
   reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
@@ -293,31 +235,25 @@ run_glm (struct casereader *input,
 	if (var_is_alpha (all_vars[i]))
 	  cat_stored_values_create (all_vars[i]);
       
-      cov_hash = covariance_hsh_create (n_all_vars);
+      cov = covariance_matrix_init (n_all_vars, all_vars, ONE_PASS, PAIRWISE, MV_ANY);
       reader = casereader_create_counter (reader, &row, -1);
       for (; casereader_read (reader, &c); case_destroy (&c))
 	{
 	  /* 
 	     Accumulate the covariance matrix.
 	  */
-	  covariance_accumulate (cov_hash, mom, &c, all_vars, n_all_vars);
+	  covariance_matrix_accumulate (cov, &c);
 	  n_data++;
 	}
-      cov = covariance_accumulator_to_matrix (cov_hash, mom, all_vars, n_all_vars, n_data);
+      covariance_matrix_compute (cov);
 
-      hsh_destroy (cov_hash);
       for (i = 0; i < n_dependent; i++)
 	{
-	  model = fit_model (cov, mom, v_dependent[i], indep_vars, n_data, n_indep);
+	  model = fit_model (cov, v_dependent[i], indep_vars, n_data, n_indep);
 	  pspp_linreg_cache_free (model);
 	}
 
       casereader_destroy (reader);
-      for (i = 0; i < n_all_vars; i++)
-	{
-	  moments1_destroy (mom[i]);
-	}
-      free (mom);
       covariance_matrix_destroy (cov);
     }
   else
diff --git a/src/math/covariance-matrix.c b/src/math/covariance-matrix.c
index 817f7ed3..fb67a64c 100644
--- a/src/math/covariance-matrix.c
+++ b/src/math/covariance-matrix.c
@@ -43,29 +43,151 @@ struct covariance_accumulator
 {
   const struct variable *v1;
   const struct variable *v2;
-  double product;
   const union value *val1;
   const union value *val2;
+  double dot_product;
+  double sum1;
+  double sum2;
+  double ssize;
 };
 
+
+
+struct covariance_matrix
+{
+  struct design_matrix *cov;
+  struct hsh_table *ca;
+  struct moments1 **m1;
+  struct moments **m;
+  const struct variable **v_variables;
+  size_t n_variables;
+  int n_pass;
+  int missing_handling;
+  enum mv_class missing_value;
+  void (*accumulate) (struct covariance_matrix *, const struct ccase *);
+  void (*update_moments) (struct covariance_matrix *, size_t, double);
+};
+
+static struct hsh_table *covariance_hsh_create (size_t);
 static hsh_hash_func covariance_accumulator_hash;
-static unsigned int hash_numeric_alpha (const struct variable *, const struct variable *, 
+static unsigned int hash_numeric_alpha (const struct variable *,
+					const struct variable *,
 					const union value *, size_t);
 static hsh_compare_func covariance_accumulator_compare;
 static hsh_free_func covariance_accumulator_free;
+static void update_moments1 (struct covariance_matrix *, size_t, double);
+static void update_moments2 (struct covariance_matrix *, size_t, double);
+static struct covariance_accumulator *get_new_covariance_accumulator (const
+								      struct
+								      variable
+								      *,
+								      const
+								      struct
+								      variable
+								      *,
+								      const
+								      union
+								      value *,
+								      const
+								      union
+								      value
+								      *);
+static void covariance_accumulate_listwise (struct covariance_matrix *,
+					    const struct ccase *);
+static void covariance_accumulate_pairwise (struct covariance_matrix *,
+					    const struct ccase *);
+
+struct covariance_matrix *
+covariance_matrix_init (size_t n_variables,
+			const struct variable *v_variables[], int n_pass,
+			int missing_handling, enum mv_class missing_value)
+{
+  size_t i;
+  struct covariance_matrix *result = NULL;
+
+  result = xmalloc (sizeof (*result));
+  result->cov = NULL;
+  result->ca = covariance_hsh_create (n_variables);
+  result->m = NULL;
+  result->m1 = NULL;
+  result->missing_handling = missing_handling;
+  result->missing_value = missing_value;
+  result->accumulate = (result->missing_handling == LISTWISE) ?
+    covariance_accumulate_listwise : covariance_accumulate_pairwise;
+  if (n_pass == ONE_PASS)
+    {
+      result->update_moments = update_moments1;
+      result->m1 = xnmalloc (n_variables, sizeof (*result->m1));
+      for (i = 0; i < n_variables; i++)
+	{
+	  result->m1[i] = moments1_create (MOMENT_MEAN);
+	}
+    }
+  else
+    {
+      result->update_moments = update_moments2;
+      result->m = xnmalloc (n_variables, sizeof (*result->m));
+      for (i = 0; i < n_variables; i++)
+	{
+	  result->m[i] = moments_create (MOMENT_MEAN);
+	}
+    }
+  result->v_variables = v_variables;
+  result->n_variables = n_variables;
+  result->n_pass = n_pass;
+
+  return result;
+}
 
 /*
   The covariances are stored in a DESIGN_MATRIX structure.
  */
 struct design_matrix *
-covariance_matrix_create (size_t n_variables, const struct variable *v_variables[])
+covariance_matrix_create (size_t n_variables,
+			  const struct variable *v_variables[])
 {
-  return design_matrix_create (n_variables, v_variables, (size_t) n_variables);
+  return design_matrix_create (n_variables, v_variables,
+			       (size_t) n_variables);
 }
 
-void covariance_matrix_destroy (struct design_matrix *x)
+static void
+update_moments1 (struct covariance_matrix *cov, size_t i, double x)
 {
-  design_matrix_destroy (x);
+  assert (cov->m1 != NULL);
+  moments1_add (cov->m1[i], x, 1.0);
+}
+
+static void
+update_moments2 (struct covariance_matrix *cov, size_t i, double x)
+{
+  assert (cov->m != NULL);
+  moments_pass_one (cov->m[i], x, 1.0);
+}
+
+void
+covariance_matrix_destroy (struct covariance_matrix *cov)
+{
+  size_t i;
+
+  assert (cov != NULL);
+  design_matrix_destroy (cov->cov);
+  hsh_destroy (cov->ca);
+  if (cov->n_pass == ONE_PASS)
+    {
+      for (i = 0; i < cov->n_variables; i++)
+	{
+	  moments1_destroy (cov->m1[i]);
+	}
+      free (cov->m1);
+    }
+  else
+    {
+      for (i = 0; i < cov->n_variables; i++)
+	{
+	  moments_destroy (cov->m[i]);
+	}
+      free (cov->m);
+    }
 }
 
 /*
@@ -74,12 +196,13 @@ void covariance_matrix_destroy (struct design_matrix *x)
  */
 static void
 covariance_update_categorical_numeric (struct design_matrix *cov, double mean,
-			  size_t row, 
-			  const struct variable *v2, double x, const union value *val2)
+				       size_t row,
+				       const struct variable *v2, double x,
+				       const union value *val2)
 {
   size_t col;
   double tmp;
-  
+
   assert (var_is_numeric (v2));
 
   col = design_matrix_var_to_column (cov, v2);
@@ -98,7 +221,7 @@ column_iterate (struct design_matrix *cov, const struct variable *v,
   double tmp;
   const union value *tmp_val;
 
-  col = design_matrix_var_to_column (cov, v);  
+  col = design_matrix_var_to_column (cov, v);
   for (i = 0; i < cat_get_n_categories (v) - 1; i++)
     {
       col += i;
@@ -113,14 +236,17 @@ column_iterate (struct design_matrix *cov, const struct variable *v,
       gsl_matrix_set (cov->m, col, row, x * y + tmp);
     }
 }
+
 /*
   Call this function in the second data pass. The central moments are
   MEAN1 and MEAN2. Any categorical variables should already have their
   values summarized in in its OBS_VALS element.
  */
-void covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
-			  double ssize, const struct variable *v1, 
-			  const struct variable *v2, const union value *val1, const union value *val2)
+void
+covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
+		     double ssize, const struct variable *v1,
+		     const struct variable *v2, const union value *val1,
+		     const union value *val2)
 {
   size_t row;
   size_t col;
@@ -142,7 +268,7 @@ void covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
 	    }
 	  if (var_is_numeric (v2))
 	    {
-	      covariance_update_categorical_numeric (cov, mean2, row, 
+	      covariance_update_categorical_numeric (cov, mean2, row,
 						     v2, x, val2);
 	    }
 	  else
@@ -155,17 +281,17 @@ void covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
   else if (var_is_alpha (v2))
     {
       /*
-	Reverse the orders of V1, V2, etc. and put ourselves back
-	in the previous IF scope.
+         Reverse the orders of V1, V2, etc. and put ourselves back
+         in the previous IF scope.
        */
       covariance_pass_two (cov, mean2, mean1, ssize, v2, v1, val2, val1);
     }
   else
     {
       /*
-	Both variables are numeric.
-      */
-      row = design_matrix_var_to_column (cov, v1);  
+         Both variables are numeric.
+       */
+      row = design_matrix_var_to_column (cov, v1);
       col = design_matrix_var_to_column (cov, v2);
       x = (val1->f - mean1) * (val2->f - mean2);
       x += gsl_matrix_get (cov->m, col, row);
@@ -187,11 +313,13 @@ covariance_accumulator_hash (const void *h, const void *aux)
   const union value *val_max;
 
   /*
-    Order everything by the variables' indices. This ensures we get the
-    same key regardless of the order in which the variables are stored
-    and passed around.
+     Order everything by the variables' indices. This ensures we get the
+     same key regardless of the order in which the variables are stored
+     and passed around.
    */
-  v_min = (var_get_dict_index (ca->v1) < var_get_dict_index (ca->v2)) ? ca->v1 : ca->v2;
+  v_min =
+    (var_get_dict_index (ca->v1) <
+     var_get_dict_index (ca->v2)) ? ca->v1 : ca->v2;
   v_max = (ca->v1 == v_min) ? ca->v2 : ca->v1;
 
   val_min = (v_min == ca->v1) ? ca->val1 : ca->val2;
@@ -215,12 +343,12 @@ covariance_accumulator_hash (const void *h, const void *aux)
   if (var_is_alpha (v_max) && var_is_alpha (v_min))
     {
       unsigned int tmp;
-      char *x = xnmalloc (1 + var_get_width (v_max) + var_get_width (v_min), sizeof (*x));
+      char *x =
+	xnmalloc (1 + var_get_width (v_max) + var_get_width (v_min),
+		  sizeof (*x));
       strncpy (x, val_max->s, var_get_width (v_max));
       strncat (x, val_min->s, var_get_width (v_min));
-      tmp = *n_vars * (*n_vars + 1 + idx_max)
-	+ idx_min
-	+ hsh_hash_string (x);
+      tmp = *n_vars * (*n_vars + 1 + idx_max) + idx_min + hsh_hash_string (x);
       free (x);
       return tmp;
     }
@@ -233,61 +361,67 @@ covariance_accumulator_hash (const void *h, const void *aux)
   in a single data pass. Call covariance_accumulate () for each case 
   in the data.
  */
-struct hsh_table *
+static struct hsh_table *
 covariance_hsh_create (size_t n_vars)
 {
-  return hsh_create (n_vars * (n_vars + 1) / 2, covariance_accumulator_compare, 
-		     covariance_accumulator_hash, covariance_accumulator_free, &n_vars);
+  return hsh_create (n_vars * n_vars, covariance_accumulator_compare,
+		     covariance_accumulator_hash, covariance_accumulator_free,
+		     &n_vars);
 }
 
-static void 
+static void
 covariance_accumulator_free (void *c_, const void *aux UNUSED)
 {
   struct covariance_accumulator *c = c_;
   assert (c != NULL);
   free (c);
 }
+
+/*
+  Hash comparison. Returns 0 for a match, or a non-zero int
+  otherwise. The sign of a non-zero return value *should* indicate the
+  position of C relative to the covariance_accumulator described by
+  the other arguments. But for now, it just returns 1 for any
+  non-match.  This should be changed when someone figures out how to
+  compute a sensible sign for the return value.
+ */
 static int
-match_nodes (const struct covariance_accumulator *c, const struct variable *v1,
-	     const struct variable *v2, const union value *val1,
-	     const union value *val2)
+match_nodes (const struct covariance_accumulator *c,
+	     const struct variable *v1, const struct variable *v2,
+	     const union value *val1, const union value *val2)
 {
-  if (var_get_dict_index (v1) == var_get_dict_index (c->v1) && 
-      var_get_dict_index (v2) == var_get_dict_index (c->v2))
-    {
-      if (var_is_numeric (v1) && var_is_numeric (v2))
-	{
-	  return 0;
-	}
-      if (var_is_numeric (v1) && var_is_alpha (v2))
-	{
-	  if (compare_values (val2, c->val2, v2))
-	    {
-	      return 0;
-	    }
-	}
-      if (var_is_alpha (v1) && var_is_numeric (v2))
-	{
-	  if (compare_values (val1, c->val1, v1))
-	    {
-	      return 0;
-	    }
-	}
-      if (var_is_alpha (v1) && var_is_alpha (v2))
-	{
-	  if (compare_values (val1, c->val1, v1))
-	    {
-	      if (compare_values (val2, c->val2, v2))
-		{
-		  return 0;
-		}
-	    }
-	}
-    }
-  else if (v2 == c->v1 && v1 == c->v2)
-    {
-      return -match_nodes (c, v2, v1, val2, val1);
-    }
+  if (var_get_dict_index (v1) == var_get_dict_index (c->v1))
+    if (var_get_dict_index (v2) == var_get_dict_index (c->v2))
+      {
+	if (var_is_numeric (v1) && var_is_numeric (v2))
+	  {
+	    return 0;
+	  }
+	if (var_is_numeric (v1) && var_is_alpha (v2))
+	  {
+	    if (compare_values (val2, c->val2, v2))
+	      {
+		return 0;
+	      }
+	  }
+	if (var_is_alpha (v1) && var_is_numeric (v2))
+	  {
+	    if (compare_values (val1, c->val1, v1))
+	      {
+		return 0;
+	      }
+	  }
+	if (var_is_alpha (v1) && var_is_alpha (v2))
+	  {
+	    if (compare_values (val1, c->val1, v1))
+	      {
+		if (compare_values (val2, c->val2, v2))
+		  {
+		    return 0;
+		  }
+	      }
+	  }
+      }
   return 1;
 }
 
@@ -296,10 +430,11 @@ match_nodes (const struct covariance_accumulator *c, const struct variable *v1,
   a struct hsh_table in src/libpspp/hash.c.
 */
 static int
-covariance_accumulator_compare (const void *a1_, const void *a2_, const void *aux UNUSED)
+covariance_accumulator_compare (const void *a1_, const void *a2_,
+				const void *aux UNUSED)
 {
-  const struct covariance_accumulator *a1 =  a1_;
-  const struct covariance_accumulator *a2 =  a2_;
+  const struct covariance_accumulator *a1 = a1_;
+  const struct covariance_accumulator *a2 = a2_;
 
   if (a1 == NULL && a2 == NULL)
     return 0;
@@ -311,7 +446,7 @@ covariance_accumulator_compare (const void *a1_, const void *a2_, const void *au
 }
 
 static unsigned int
-hash_numeric_alpha (const struct variable *v1, const struct variable *v2, 
+hash_numeric_alpha (const struct variable *v1, const struct variable *v2,
 		    const union value *val, size_t n_vars)
 {
   unsigned int result = -1u;
@@ -329,8 +464,8 @@ hash_numeric_alpha (const struct variable *v1, const struct variable *v2,
 
 
 static double
-update_product (const struct variable *v1, const struct variable *v2, const union value *val1,
-		const union value *val2)
+update_product (const struct variable *v1, const struct variable *v2,
+		const union value *val1, const union value *val2)
 {
   assert (v1 != NULL);
   assert (v2 != NULL);
@@ -354,60 +489,192 @@ update_product (const struct variable *v1, const struct variable *v2, const unio
     }
   return 0.0;
 }
+static double
+update_sum (const struct variable *var, const union value *val)
+{
+  assert (var != NULL);
+  assert (val != NULL);
+  if (var_is_alpha (var))
+    {
+      return 1.0;
+    }
+  return val->f;
+}
+static struct covariance_accumulator *
+get_new_covariance_accumulator (const struct variable *v1,
+				const struct variable *v2,
+				const union value *val1,
+				const union value *val2)
+{
+  if ((v1 != NULL) && (v2 != NULL) && (val1 != NULL) && (val2 != NULL))
+    {
+      struct covariance_accumulator *ca;
+      ca = xmalloc (sizeof (*ca));
+      ca->v1 = v1;
+      ca->v2 = v2;
+      ca->val1 = val1;
+      ca->val2 = val2;
+      return ca;
+    }
+  return NULL;
+}
+
+static const struct variable **
+get_covariance_variables (const struct covariance_matrix *cov)
+{
+  return cov->v_variables;
+}
+
+static void
+update_hash_entry (struct hsh_table *c,
+		   const struct variable *v1,
+		   const struct variable *v2,
+		   const union value *val1, const union value *val2)
+{
+  struct covariance_accumulator *ca;
+  struct covariance_accumulator *new_entry;
+
+
+  ca = get_new_covariance_accumulator (v1, v2, val1, val2);
+  ca->dot_product = update_product (ca->v1, ca->v2, ca->val1, ca->val2);
+  ca->sum1 = update_sum (ca->v1, ca->val1);
+  ca->sum2 = update_sum (ca->v2, ca->val2);
+  ca->ssize = 1.0;
+  new_entry = hsh_insert (c, ca);
+  if (new_entry != NULL)
+    {
+      new_entry->dot_product += ca->dot_product;
+      new_entry->ssize += 1.0;
+      new_entry->sum1 += ca->sum1;
+      new_entry->sum2 += ca->sum2;
+      /*
+         If DOT_PRODUCT is null, CA was not already in the hash
+         hable, so we don't free it because it was just inserted.
+         If DOT_PRODUCT was not null, CA is already in the hash table.
+         Unnecessary now, it must be freed here.
+       */
+      free (ca);
+    }
+}
+
 /*
-  Compute the covariance matrix in a single data-pass.
+  Compute the covariance matrix in a single data-pass. Cases with
+  missing values are dropped pairwise, in other words, only if one of
+  the two values necessary to accumulate the inner product is missing.
+
+  Do not call this function directly. Call it through the struct
+  covariance_matrix ACCUMULATE member function, for example,
+  cov->accumulate (cov, ccase).
  */
-void 
-covariance_accumulate (struct hsh_table *cov, struct moments1 **m,
-		       const struct ccase *ccase, const struct variable **vars,
-		       size_t n_vars)
+static void
+covariance_accumulate_pairwise (struct covariance_matrix *cov,
+				const struct ccase *ccase)
 {
   size_t i;
   size_t j;
-  const union value *val;
-  struct covariance_accumulator *ca;
-  struct covariance_accumulator *entry;
+  const union value *val1;
+  const union value *val2;
+  const struct variable **v_variables;
+
+  assert (cov != NULL);
+  assert (ccase != NULL);
 
-  assert (m != NULL);
+  v_variables = get_covariance_variables (cov);
+  assert (v_variables != NULL);
 
-  for (i = 0; i < n_vars; ++i)
+  for (i = 0; i < cov->n_variables; ++i)
     {
-      val = case_data (ccase, vars[i]);
-      if (var_is_alpha (vars[i]))
-	{
-	  cat_value_update (vars[i], val);
-	}
-      else
+      val1 = case_data (ccase, v_variables[i]);
+      if (!var_is_value_missing (v_variables[i], val1, cov->missing_value))
 	{
-	  moments1_add (m[i], val->f, 1.0);
-	}
-      for (j = i; j < n_vars; j++)
-	{
-	  ca = xmalloc (sizeof (*ca));
-	  ca->v1 = vars[i];
-	  ca->v2 = vars[j];
-	  ca->val1 = val;
-	  ca->val2 = case_data (ccase, ca->v2);
-	  ca->product = update_product (ca->v1, ca->v2, ca->val1, ca->val2);
-	  entry = hsh_insert (cov, ca);
-	  if (entry != NULL)
+	  cat_value_update (v_variables[i], val1);
+	  if (var_is_alpha (v_variables[i]))
+	    cov->update_moments (cov, i, val1->f);
+
+	  for (j = i; j < cov->n_variables; j++)
 	    {
-	      entry->product += ca->product;
-	      /*
-		If ENTRY is null, CA was not already in the hash
-		hable, so we don't free it because it was just inserted.
-		If ENTRY was not null, CA is already in the hash table.
-		Unnecessary now, it must be freed here.
-	      */
-	      free (ca);
+	      val2 = case_data (ccase, v_variables[j]);
+	      if (!var_is_value_missing
+		  (v_variables[j], val2, cov->missing_value))
+		{
+		  update_hash_entry (cov->ca, v_variables[i], v_variables[j],
+				     val1, val2);
+		  if (j != i)
+		    update_hash_entry (cov->ca, v_variables[j],
+				       v_variables[i], val2, val1);
+		}
 	    }
 	}
     }
 }
 
-static void 
-covariance_matrix_insert (struct design_matrix *cov, const struct variable *v1,
-			  const struct variable *v2, const union value *val1, 
+/*
+  Compute the covariance matrix in a single data-pass. Cases with
+  missing values are dropped listwise. In other words, if one of the
+  values for any variable in a case is missing, the entire case is
+  skipped. 
+
+  The caller must use a casefilter to remove the cases with missing
+  values before calling covariance_accumulate_listwise. This function
+  assumes that CCASE has already passed through this filter, and
+  contains no missing values.
+
+  Do not call this function directly. Call it through the struct
+  covariance_matrix ACCUMULATE member function, for example,
+  cov->accumulate (cov, ccase).
+ */
+static void
+covariance_accumulate_listwise (struct covariance_matrix *cov,
+				const struct ccase *ccase)
+{
+  size_t i;
+  size_t j;
+  const union value *val1;
+  const union value *val2;
+  const struct variable **v_variables;
+
+  assert (cov != NULL);
+  assert (ccase != NULL);
+
+  v_variables = get_covariance_variables (cov);
+  assert (v_variables != NULL);
+
+  for (i = 0; i < cov->n_variables; ++i)
+    {
+      val1 = case_data (ccase, v_variables[i]);
+      cat_value_update (v_variables[i], val1);
+      if (var_is_alpha (v_variables[i]))
+	cov->update_moments (cov, i, val1->f);
+
+      for (j = i; j < cov->n_variables; j++)
+	{
+	  val2 = case_data (ccase, v_variables[j]);
+	  update_hash_entry (cov->ca, v_variables[i], v_variables[j],
+			     val1, val2);
+	  if (j != i)
+	    update_hash_entry (cov->ca, v_variables[j], v_variables[i],
+			       val2, val1);
+	}
+    }
+}
+
+/*
+  Call this function during the data pass. Each case will be added to
+  a hash containing all values of the covariance matrix. After the
+  data have been passed, call covariance_matrix_compute to put the
+  values in the struct covariance_matrix.
+ */
+void
+covariance_matrix_accumulate (struct covariance_matrix *cov,
+			      const struct ccase *ccase)
+{
+  cov->accumulate (cov, ccase);
+}
+
+static void
+covariance_matrix_insert (struct design_matrix *cov,
+			  const struct variable *v1,
+			  const struct variable *v2, const union value *val1,
 			  const union value *val2, double product)
 {
   size_t row;
@@ -441,10 +708,10 @@ covariance_matrix_insert (struct design_matrix *cov, const struct variable *v1,
 	    {
 	      i++;
 	      tmp_val = cat_subscript_to_value (i, v1);
-	    } 
+	    }
 	  col += i;
 	}
-    }    
+    }
   else
     {
       if (var_is_numeric (v2))
@@ -457,80 +724,52 @@ covariance_matrix_insert (struct design_matrix *cov, const struct variable *v1,
 	}
     }
   gsl_matrix_set (cov->m, row, col, product);
-  gsl_matrix_set (cov->m, col, row, product);
-}
-
-static double
-get_center (const struct variable *v, const union value *val, 
-	    const struct variable **vars, const struct moments1 **m, size_t n_vars,
-	    size_t ssize)
-{
-  size_t i = 0;
-
-  while ((var_get_dict_index (vars[i]) != var_get_dict_index(v)) && (i < n_vars))
-    {
-      i++;
-    }  
-  if (var_is_numeric (v))
-    {
-      double mean;
-      moments1_calculate (m[i], NULL, &mean, NULL, NULL, NULL);
-      return mean;
-    }
-  else 
-    {
-      i = cat_value_find (v, val);
-      return (cat_get_category_count (i, v) / ssize);
-    }
-  return 0.0;
 }
 
-/*
-  Subtract the product of the means.
- */
-static double
-center_entry (const struct covariance_accumulator *ca, const struct variable **vars,
-	      const struct moments1 **m, size_t n_vars, size_t ssize)
-{
-  double m1;
-  double m2;
-  double result = 0.0;
-  
-  m1 = get_center (ca->v1, ca->val1, vars, m, n_vars, ssize);
-  m2 = get_center (ca->v2, ca->val2, vars, m, n_vars, ssize);
-  result = ca->product - ssize * m1 * m2;
-  return result;
-}
-
-/*
-  The first moments in M should be stored in the order corresponding
-  to the order of VARS. So, for example, VARS[0] has its moments in
-  M[0], VARS[1] has its moments in M[1], etc.
- */
-struct design_matrix *
-covariance_accumulator_to_matrix (struct hsh_table *cov, const struct moments1 **m,
-				  const struct variable **vars, size_t n_vars, size_t ssize)
+static struct design_matrix *
+covariance_accumulator_to_matrix (struct covariance_matrix *cov)
 {
   double tmp;
   struct covariance_accumulator *entry;
   struct design_matrix *result = NULL;
   struct hsh_iterator iter;
-  
-  result = covariance_matrix_create (n_vars, vars);
 
-  entry = hsh_first (cov, &iter);
-  
+  result = covariance_matrix_create (cov->n_variables, cov->v_variables);
+
+  entry = hsh_first (cov->ca, &iter);
+
   while (entry != NULL)
     {
       /*
-	We compute the centered, un-normalized covariance matrix.
+         We compute the centered, un-normalized covariance matrix.
        */
-      tmp = center_entry (entry, vars, m, n_vars, ssize);
+      tmp = entry->dot_product - entry->sum1 * entry->sum2 / entry->ssize;
       covariance_matrix_insert (result, entry->v1, entry->v2, entry->val1,
 				entry->val2, tmp);
-      entry = hsh_next (cov, &iter);
+      entry = hsh_next (cov->ca, &iter);
     }
-
   return result;
 }
 
+
+/*
+  Call this function after passing the data.
+ */
+void
+covariance_matrix_compute (struct covariance_matrix *cov)
+{
+  if (cov->n_pass == ONE_PASS)
+    {
+      cov->cov = covariance_accumulator_to_matrix (cov);
+    }
+}
+
+struct design_matrix *
+covariance_to_design (const struct covariance_matrix *c)
+{
+  if (c != NULL)
+    {
+      return c->cov;
+    }
+  return NULL;
+}
diff --git a/src/math/covariance-matrix.h b/src/math/covariance-matrix.h
index 573d3e7e..33a5d750 100644
--- a/src/math/covariance-matrix.h
+++ b/src/math/covariance-matrix.h
@@ -26,18 +26,33 @@
 struct moments1;
 struct ccase;
 struct hsh_table;
+struct covariance_matrix;
+enum
+{ ONE_PASS,
+  TWO_PASS
+};
 
-struct design_matrix *
-covariance_matrix_create (size_t, const struct variable *[]);
-
-void covariance_matrix_destroy (struct design_matrix *);
-
+/*
+  How to deal with missing values.
+ */
+enum
+{ LISTWISE,
+  PAIRWISE
+};
+struct design_matrix *covariance_matrix_create (size_t,
+						const struct variable *[]);
+
+void covariance_matrix_destroy (struct covariance_matrix *cov);
 void covariance_pass_two (struct design_matrix *, double,
-			  double, double, const struct variable *, 
-			  const struct variable *, const union value *, const union value *);
-void covariance_accumulate (struct hsh_table *, struct moments1 **,
-			    const struct ccase *, const struct variable **, size_t);
-struct hsh_table * covariance_hsh_create (size_t);
-struct design_matrix * covariance_accumulator_to_matrix (struct hsh_table *, const struct moments1 **,
-							 const struct variable **, size_t, size_t);
+			  double, double, const struct variable *,
+			  const struct variable *, const union value *,
+			  const union value *);
+void covariance_matrix_compute (struct covariance_matrix *);
+struct covariance_matrix *covariance_matrix_init (size_t,
+						  const struct variable *[],
+						  int, int, enum mv_class);
+void covariance_matrix_free (struct covariance_matrix *);
+void covariance_matrix_accumulate (struct covariance_matrix *,
+				   const struct ccase *);
+struct design_matrix *covariance_to_design (const struct covariance_matrix *);
 #endif