X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance-matrix.c;h=89660ba9c1e3f5fbc1f4168af2ef8af549a43c7d;hb=8829e1f148c279db0f19b9e3bd746ace07d2d7f1;hp=a33c33b243e46c07c8997a8825f27b00cfbc458b;hpb=d3181522c4d8e725d417cbae655e46f660fe168e;p=pspp-builds.git

diff --git a/src/math/covariance-matrix.c b/src/math/covariance-matrix.c
index a33c33b2..89660ba9 100644
--- a/src/math/covariance-matrix.c
+++ b/src/math/covariance-matrix.c
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2008 Free Software Foundation, Inc.
+   Copyright (C) 2008, 2009 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
@@ -57,12 +57,13 @@ struct covariance_matrix
 {
   struct design_matrix *cov;
   struct design_matrix *ssize;
-  struct design_matrix *means;
   struct hsh_table *ca;
   struct moments1 **m1;
   struct moments **m;
   const struct variable **v_variables;
+  const struct interaction_variable **interactions;
   size_t n_variables;
+  size_t n_intr;
   int n_pass;
   int missing_handling;
   enum mv_class missing_value;
@@ -120,6 +121,7 @@ covariance_matrix_init (size_t n_variables,
   result->ca = covariance_hsh_create (&result->n_variables);
   result->m = NULL;
   result->m1 = NULL;
+  result->n_intr = 0;
   result->missing_handling = missing_handling;
   result->missing_value = missing_value;
   result->accumulate = (result->missing_handling == LISTWISE) ?
@@ -148,7 +150,33 @@ covariance_matrix_init (size_t n_variables,
 
   return result;
 }
+void
+covariance_interaction_set (struct covariance_matrix *cov, 
+			    const struct interaction_variable **intr, size_t n_intr)
+{
+  cov->interactions = intr;
+  cov->n_intr = n_intr;
+}
 
+static size_t 
+get_n_rows (size_t n_variables, const struct variable *v_variables[])
+{
+  size_t i;
+  size_t result = 0;
+  for (i = 0; i < n_variables; i++)
+    {
+      if (var_is_numeric (v_variables[i]))
+	{
+	  result++;
+	}
+      else if (var_is_alpha (v_variables[i]))
+	{
+	  size_t n_categories = cat_get_n_categories (v_variables[i]);
+	  result += n_categories - 1;
+	}
+    }
+  return result;
+}
 /*
   The covariances are stored in a DESIGN_MATRIX structure.
  */
@@ -156,8 +184,46 @@ struct design_matrix *
 covariance_matrix_create (size_t n_variables,
 			  const struct variable *v_variables[])
 {
-  return design_matrix_create (n_variables, v_variables,
-			       (size_t) n_variables);
+  size_t n_rows = get_n_rows (n_variables, v_variables);
+  return design_matrix_create (n_variables, v_variables, n_rows);
+}
+
+static size_t 
+get_n_rows_s (const struct variable *var)
+{
+  size_t result = 0;
+  if (var_is_numeric (var))
+    {
+      result++;
+    }
+  else
+    {
+      result += cat_get_n_categories (var) - 1;
+    }
+  return result;
+}
+static struct design_matrix *
+covariance_matrix_create_s (struct covariance_matrix *cov)
+{
+  struct variable **v_variables;
+  size_t n_variables;
+  size_t n_rows = 0;
+  size_t i;
+  size_t j;
+
+  n_variables = cov->n_variables + cov->n_intr;
+  v_variables = xnmalloc (n_variables, sizeof (*v_variables));
+  for (i = 0; i < cov->n_variables; i++)
+    {
+      v_variables[i] = cov->v_variables[i];
+      n_rows += get_n_rows_s (v_variables[i]);
+    }
+  for (j = 0; j < cov->n_intr; j++)
+    {
+      v_variables[i + j] = interaction_get_variable (cov->interactions[j]);
+      n_rows += get_n_rows_s (v_variables[i]);
+    }
+  return design_matrix_create (n_variables, v_variables, n_rows);
 }
 
 static void
@@ -182,7 +248,6 @@ covariance_matrix_destroy (struct covariance_matrix *cov)
   assert (cov != NULL);
   design_matrix_destroy (cov->cov);
   design_matrix_destroy (cov->ssize);
-  design_matrix_destroy (cov->means);
   hsh_destroy (cov->ca);
   if (cov->n_pass == ONE_PASS)
     {
@@ -219,14 +284,15 @@ covariance_update_categorical_numeric (struct design_matrix *cov, double mean,
 
   col = design_matrix_var_to_column (cov, v2);
   assert (val2 != NULL);
-  tmp = gsl_matrix_get (cov->m, row, col);
-  gsl_matrix_set (cov->m, row, col, (val2->f - mean) * x + tmp);
-  gsl_matrix_set (cov->m, col, row, (val2->f - mean) * x + tmp);
+  tmp = design_matrix_get_element (cov, row, col);
+  design_matrix_set_element (cov, row, col, (val2->f - mean) * x + tmp);
+  design_matrix_set_element (cov, col, row, (val2->f - mean) * x + tmp);
 }
 static void
 column_iterate (struct design_matrix *cov, const struct variable *v,
 		double ssize, double x, const union value *val1, size_t row)
 {
+  int width = var_get_width (v);
   size_t col;
   size_t i;
   double y;
@@ -239,13 +305,13 @@ column_iterate (struct design_matrix *cov, const struct variable *v,
       col += i;
       y = -1.0 * cat_get_category_count (i, v) / ssize;
       tmp_val = cat_subscript_to_value (i, v);
-      if (!compare_values_short (tmp_val, val1, v))
+      if (!value_equal (tmp_val, val1, width))
 	{
 	  y += -1.0;
 	}
-      tmp = gsl_matrix_get (cov->m, row, col);
-      gsl_matrix_set (cov->m, row, col, x * y + tmp);
-      gsl_matrix_set (cov->m, col, row, x * y + tmp);
+      tmp = design_matrix_get_element (cov, row, col);
+      design_matrix_set_element (cov, row, col, x * y + tmp);
+      design_matrix_set_element (cov, col, row, x * y + tmp);
     }
 }
 
@@ -274,7 +340,7 @@ covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
 	  row += i;
 	  x = -1.0 * cat_get_category_count (i, v1) / ssize;
 	  tmp_val = cat_subscript_to_value (i, v1);
-	  if (!compare_values_short (tmp_val, val1, v1))
+	  if (!value_equal (tmp_val, val1, var_get_width (v1)))
 	    {
 	      x += 1.0;
 	    }
@@ -306,9 +372,9 @@ covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
       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);
-      gsl_matrix_set (cov->m, row, col, x);
-      gsl_matrix_set (cov->m, col, row, x);
+      x += design_matrix_get_element (cov, col, row);
+      design_matrix_set_element (cov, row, col, x);
+      design_matrix_set_element (cov, col, row, x);
     }
 }
 
@@ -354,10 +420,9 @@ covariance_accumulator_hash (const void *h, const void *aux)
     }
   if (var_is_alpha (v_max) && var_is_alpha (v_min))
     {
-      unsigned tmp = hsh_hash_bytes (val_max, var_get_width (v_max));
-      tmp ^= hsh_hash_bytes (val_min, var_get_width (v_min));
-      tmp += *n_vars * (*n_vars + 1 + idx_max) + idx_min;
-      return (size_t) tmp;
+      unsigned hash = value_hash (val_max, var_get_width (v_max), 0);
+      hash = value_hash (val_min, var_get_width (v_min), hash);
+      return hash_int (*n_vars * (*n_vars + 1 + idx_max) + idx_min, hash);
     }
   return -1u;
 }
@@ -384,6 +449,31 @@ covariance_accumulator_free (void *c_, const void *aux UNUSED)
   free (c);
 }
 
+static int 
+ordered_match_nodes (const struct covariance_accumulator *c, const struct variable *v1,
+		     const struct variable *v2, const union value *val1, const union value *val2)
+{
+  size_t result;
+  size_t m;
+
+  result = var_get_dict_index (v1) ^ var_get_dict_index (c->v1);
+  m = var_get_dict_index (v2) ^ var_get_dict_index (c->v2);
+  result = result|m;
+  if (var_is_alpha (v1))
+    {
+      result |= value_compare_3way (val1, c->val1, var_get_width (v1));
+      if (var_is_alpha (v2))
+	{
+	  result |= value_compare_3way (val2, c->val2, var_get_width (v2));
+	}
+    }
+  else if (var_is_alpha (v2))
+    {
+      result |= value_compare_3way (val2, c->val2, var_get_width (v2));
+    }
+  return result;
+}
+  
 /*
   Hash comparison. Returns 0 for a match, or a non-zero int
   otherwise. The sign of a non-zero return value *should* indicate the
@@ -397,39 +487,12 @@ 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))
-    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_short (val2, c->val2, v2))
-	      {
-		return 0;
-	      }
-	  }
-	if (var_is_alpha (v1) && var_is_numeric (v2))
-	  {
-	    if (!compare_values_short (val1, c->val1, v1))
-	      {
-		return 0;
-	      }
-	  }
-	if (var_is_alpha (v1) && var_is_alpha (v2))
-	  {
-	    if (!compare_values_short (val1, c->val1, v1))
-	      {
-		if (!compare_values_short (val2, c->val2, v2))
-		  {
-		    return 0;
-		  }
-	      }
-	  }
-      }
-  return 1;
+  size_t n;
+  size_t m;
+
+  n = ordered_match_nodes (c, v1, v2, val1, val2);
+  m = ordered_match_nodes (c, v2, v1, val2, val1);
+  return (n & m);
 }
 
 /*
@@ -448,7 +511,7 @@ covariance_accumulator_compare (const void *a1_, const void *a2_,
 
   if (a1 == NULL || a2 == NULL)
     return 1;
-
+  
   return match_nodes (a1, a2->v1, a2->v2, a2->val1, a2->val2);
 }
 
@@ -460,7 +523,7 @@ hash_numeric_alpha (const struct variable *v1, const struct variable *v2,
   if (var_is_numeric (v1) && var_is_alpha (v2))
     {
       result = n_vars * ((n_vars + 1) + var_get_dict_index (v1))
-	+ var_get_dict_index (v2) + hsh_hash_string (val->s);
+	+ var_get_dict_index (v2) + value_hash (val, var_get_width (v2), 0);
     }
   else if (var_is_alpha (v1) && var_is_numeric (v2))
     {
@@ -488,13 +551,16 @@ update_product (const struct variable *v1, const struct variable *v2,
     }
   if (var_is_numeric (v1) && var_is_alpha (v2))
     {
-      return (val1->f);
+      return val1->f;
     }
   if (var_is_numeric (v2) && var_is_alpha (v1))
     {
-      update_product (v2, v1, val2, val1);
+      return val2->f;
+    }
+  else
+    {
+      return 0.0;
     }
-  return 0.0;
 }
 static double
 update_sum (const struct variable *var, const union value *val, double weight)
@@ -533,12 +599,12 @@ get_covariance_variables (const struct covariance_matrix *cov)
 }
 
 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, 
-		   const struct interaction_value *i_val1,
-		   const struct interaction_value *i_val2)
+update_hash_entry_intr (struct hsh_table *c,
+			const struct variable *v1,
+			const struct variable *v2,
+			const union value *val1, const union value *val2, 
+			const struct interaction_value *i_val1,
+			const struct interaction_value *i_val2)
 {
   struct covariance_accumulator *ca;
   struct covariance_accumulator *new_entry;
@@ -554,7 +620,39 @@ update_hash_entry (struct hsh_table *c,
   ca->sum2 = update_sum (ca->v2, ca->val2, iv_f2);
   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);
+    }
+}
+
+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, 1.0);
+  ca->sum2 = update_sum (ca->v2, ca->val2, 1.0);
+  ca->ssize = 1.0;
+  new_entry = hsh_insert (c, ca);
+
   if (new_entry != NULL)
     {
       new_entry->dot_product += ca->dot_product;
@@ -571,6 +669,24 @@ update_hash_entry (struct hsh_table *c,
     }
 }
 
+static void
+inner_intr_loop (struct covariance_matrix *cov, const struct ccase  *ccase, const struct variable *var1,
+		 const union value *val1, const struct interaction_variable **i_var, 
+		 const struct interaction_value *i_val1, size_t j)
+{
+  struct variable *var2;
+  union value *val2;
+  struct interaction_value *i_val2;
+
+  var2 = interaction_get_variable (i_var[j]);
+  i_val2 = interaction_case_data (ccase, i_var[j]);
+  val2 = interaction_value_get (i_val2);
+  
+  if (!var_is_value_missing (var2, val2, cov->missing_value))
+    {
+      update_hash_entry_intr (cov->ca, var1, var2, val1, val2, i_val1, i_val2);
+    }
+}	 
 /*
   Compute the covariance matrix in a single data-pass. Cases with
   missing values are dropped pairwise, in other words, only if one of
@@ -591,6 +707,8 @@ covariance_accumulate_pairwise (struct covariance_matrix *cov,
   const union value *val1;
   const union value *val2;
   const struct variable **v_variables;
+  const struct variable *var1;
+  const struct variable *var2;
   struct interaction_value *i_val1 = NULL;
   struct interaction_value *i_val2 = NULL;
 
@@ -602,42 +720,42 @@ covariance_accumulate_pairwise (struct covariance_matrix *cov,
 
   for (i = 0; i < cov->n_variables; ++i)
     {
-      if (is_interaction (v_variables[i], i_var, n_intr))
-	{
-	  i_val1 = interaction_case_data (ccase, v_variables[i], i_var, n_intr);
-	  val1 = interaction_value_get (i_val1);
-	}
-      else
-	{
-	  val1 = case_data (ccase, v_variables[i]);
-	}
-      if (!var_is_value_missing (v_variables[i], val1, cov->missing_value))
+      var1 = v_variables[i];
+      val1 = case_data (ccase, var1);
+      if (!var_is_value_missing (var1, val1, cov->missing_value))
 	{
-	  cat_value_update (v_variables[i], val1);
-	  if (var_is_numeric (v_variables[i]))
+	  cat_value_update (var1, val1);
+	  if (var_is_numeric (var1))
 	    cov->update_moments (cov, i, val1->f);
 
 	  for (j = i; j < cov->n_variables; j++)
 	    {
-	      if (is_interaction (v_variables[j], i_var, n_intr))
-		{
-		  i_val2 = interaction_case_data (ccase, v_variables[j], i_var, n_intr);
-		  val2 = interaction_value_get (i_val2);
-		}
-	      else
-		{
-		  val2 = case_data (ccase, v_variables[j]);
-		}
+	      var2 = v_variables[j];
+	      val2 = case_data (ccase, var2);
 	      if (!var_is_value_missing
-		  (v_variables[j], val2, cov->missing_value))
+		  (var2, val2, cov->missing_value))
 		{
-		  update_hash_entry (cov->ca, v_variables[i], v_variables[j],
-				     val1, val2, i_val1, i_val2);
-		  if (j != i)
-		    update_hash_entry (cov->ca, v_variables[j],
-				       v_variables[i], val2, val1, i_val2, i_val1);
+		  update_hash_entry (cov->ca, var1, var2, val1, val2);
 		}
 	    }
+	  for (j = 0; j < cov->n_intr; j++)
+	    {
+	      inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
+	    }
+	}
+    }
+  for (i = 0; i < cov->n_intr; i++)
+    {
+      var1 = interaction_get_variable (i_var[i]);
+      i_val1 = interaction_case_data (ccase, i_var[i]);
+      val1 = interaction_value_get (i_val1);
+      cat_value_update (var1, val1);
+      if (!var_is_value_missing (var1, val1, cov->missing_value))
+	{
+	  for (j = i; j < cov->n_intr; j++)
+	    {
+	      inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
+	    }
 	}
     }
 }
@@ -679,35 +797,15 @@ covariance_accumulate_listwise (struct covariance_matrix *cov,
 
   for (i = 0; i < cov->n_variables; ++i)
     {
-      if (is_interaction (v_variables[i], i_var, n_intr))
-	{
-	  i_val1 = interaction_case_data (ccase, v_variables[i], i_var, n_intr);
-	  val1 = interaction_value_get (i_val1);
-	}
-      else
-	{
-	  val1 = case_data (ccase, v_variables[i]);
-	}
+      val1 = case_data (ccase, v_variables[i]);
       cat_value_update (v_variables[i], val1);
       if (var_is_numeric (v_variables[i]))
 	cov->update_moments (cov, i, val1->f);
 
       for (j = i; j < cov->n_variables; j++)
 	{
-	  if (is_interaction (v_variables[j], i_var, n_intr))
-	    {
-	      i_val2 = interaction_case_data (ccase, v_variables[j], i_var, n_intr);
-	      val2 = interaction_value_get (i_val2);
-	    }
-	  else
-	    {
-	      val2 = case_data (ccase, v_variables[j]);
-	    }
 	  update_hash_entry (cov->ca, v_variables[i], v_variables[j],
-			     val1, val2, i_val1, i_val2);
-	  if (j != i)
-	    update_hash_entry (cov->ca, v_variables[j], v_variables[i],
-			       val2, val1, i_val2, i_val1);
+			     val1, val2);
 	}
     }
 }
@@ -724,65 +822,34 @@ covariance_matrix_accumulate (struct covariance_matrix *cov,
 {
   cov->accumulate (cov, ccase, (const struct interaction_variable **) aux, n_intr);
 }
-/*
-  If VAR is categorical with d categories, its first category should
-  correspond to the origin in d-dimensional Euclidean space.
- */
-static bool
-is_origin (const struct variable *var, const union value *val)
-{
-  if (cat_value_find (var, val) == 0)
-    {
-      return true;
-    }
-  return false;
-}
 
 /*
-  Return the subscript of the column of the design matrix
-  corresponding to VAL. If VAR is categorical with d categories, its
-  first category should correspond to the origin in d-dimensional
-  Euclidean space, so there is no subscript for this value.
+  Return the value corresponding to subscript TARGET. If that value corresponds
+  to the origin, return NULL.
  */
-static size_t
-get_exact_subscript (const struct design_matrix *dm, const struct variable *var,
-		     const union value *val)
+static const union value *
+get_value_from_subscript (const struct design_matrix *dm, size_t target)
 {
-  size_t result;
-
-  if (is_origin (var, val))
-    {
-      return -1u;
-    }
-
-  result = design_matrix_var_to_column (dm, var);
-  if (var_is_alpha (var))
+  const union value *result = NULL;
+  const struct variable *var;
+  size_t i;
+  
+  var = design_matrix_col_to_var (dm, target);
+  if (var_is_numeric (var))
     {
-      result += cat_value_find (var, val) - 1;
+      return NULL;
     }
-  return result;
-}
-
-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;
-  size_t col;
-
-  assert (cov != NULL);
-
-  row = get_exact_subscript (cov, v1, val1);
-  col = get_exact_subscript (cov, v2, val2);
-  if (row != -1u && col != -1u)
+  for (i = 0; i < cat_get_n_categories (var); i++)
     {
-      gsl_matrix_set (cov->m, row, col, product);
+      result = cat_subscript_to_value (i, var);
+      if (dm_get_exact_subscript (dm, var, result) == target)
+	{
+	  return result;
+	}
     }
+  return NULL;
 }
 
-
 static bool
 is_covariance_contributor (const struct covariance_accumulator *ca, const struct design_matrix *dm,
 			   size_t i, size_t j)
@@ -790,18 +857,33 @@ is_covariance_contributor (const struct covariance_accumulator *ca, const struct
   size_t k;
   const struct variable *v1;
   const struct variable *v2;
-
+  
   assert (dm != NULL);
   v1 = design_matrix_col_to_var (dm, i);
+  v2 = design_matrix_col_to_var (dm, j);
   if (var_get_dict_index (v1) == var_get_dict_index(ca->v1))
     {
-      v2 = design_matrix_col_to_var (dm, j);
       if (var_get_dict_index (v2) == var_get_dict_index (ca->v2))
 	{
-	  k = get_exact_subscript (dm, v1, ca->val1);
+	  k = dm_get_exact_subscript (dm, v1, ca->val1);
+	  if (k == i)
+	    {
+	      k = dm_get_exact_subscript (dm, v2, ca->val2);
+	      if (k == j)
+		{
+		  return true;
+		}
+	    }
+	}
+    }
+  else if (var_get_dict_index (v1) == var_get_dict_index (ca->v2))
+    {
+      if (var_get_dict_index (v2) == var_get_dict_index (ca->v1))
+	{
+	  k = dm_get_exact_subscript (dm, v1, ca->val2);
 	  if (k == i)
 	    {
-	      k = get_exact_subscript (dm, v2, ca->val2);
+	      k = dm_get_exact_subscript (dm, v2, ca->val1);
 	      if (k == j)
 		{
 		  return true;
@@ -809,12 +891,49 @@ is_covariance_contributor (const struct covariance_accumulator *ca, const struct
 	    }
 	}
     }
+  
   return false;
 }
+static double
+get_sum (const struct covariance_matrix *cov, size_t i)
+{
+  size_t k;
+  double mean;
+  double n;
+  const struct variable *var;
+  const union value *val = NULL;
+
+  assert ( cov != NULL);
+  var = design_matrix_col_to_var (cov->cov, i);
+  if (var != NULL)
+    {
+      if (var_is_alpha (var))
+	{
+	  val = get_value_from_subscript (cov->cov, i);
+	  k = cat_value_find (var, val);
+	  return cat_get_category_count (k, var);
+	}
+      else
+	{
+	  k = 0;
+	  while (cov->v_variables[k] != var && k  < cov->n_variables)
+	    {
+	      k++;
+	    }
+	  if (k < cov->n_variables)
+	    {
+	      moments1_calculate (cov->m1[k], &n, &mean, NULL, NULL, NULL);
+	      return mean * n;
+	    }
+	}
+    }
+      
+  return 0.0;
+}
 static void
 update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_accumulator *ca)
 {
-  struct variable *var;
+  const struct variable *var;
   double tmp;
   var = design_matrix_col_to_var (dm, i);
   if (var_get_dict_index (ca->v1) == var_get_dict_index (var))
@@ -822,9 +941,9 @@ update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_ac
       var = design_matrix_col_to_var (dm, j);
       if (var_get_dict_index (ca->v2) == var_get_dict_index (var))
 	{
-	  tmp = gsl_matrix_get (dm->m, i, j);
+	  tmp = design_matrix_get_element (dm, i, j);
 	  tmp += ca->ssize;
-	  gsl_matrix_set (dm->m, i, j, tmp);
+	  design_matrix_set_element (dm, i, j, tmp);
 	}
     }
 }
@@ -833,38 +952,43 @@ covariance_accumulator_to_matrix (struct covariance_matrix *cov)
 {
   size_t i;
   size_t j;
-  double tmp;
+  double sum_i = 0.0;
+  double sum_j = 0.0;
+  double tmp = 0.0;
   struct covariance_accumulator *entry;
   struct hsh_iterator iter;
 
-  cov->cov = covariance_matrix_create (cov->n_variables, cov->v_variables);
-  cov->ssize = covariance_matrix_create (cov->n_variables, cov->v_variables);
-  cov->means = covariance_matrix_create (cov->n_variables, cov->v_variables);
+  cov->cov = covariance_matrix_create_s (cov);
+  cov->ssize = covariance_matrix_create_s (cov);
+  entry = hsh_first (cov->ca, &iter);
+  while (entry != NULL)
+    {
+      entry = hsh_next (cov->ca, &iter);
+    }
+
   for (i = 0; i < design_matrix_get_n_cols (cov->cov); i++)
     {
+      sum_i = get_sum (cov, i);
       for (j = i; j < design_matrix_get_n_cols (cov->cov); j++)
 	{
+	  sum_j = get_sum (cov, j);
 	  entry = hsh_first (cov->ca, &iter);
-	  
 	  while (entry != NULL)
 	    {
 	      update_ssize (cov->ssize, i, j, entry);
-
 	      /*
 		We compute the centered, un-normalized covariance matrix.
 	      */
 	      if (is_covariance_contributor (entry, cov->cov, i, j))
 		{
-		  covariance_matrix_insert (cov->cov, entry->v1, entry->v2, entry->val1,
-					    entry->val2, entry->dot_product);
-		  covariance_matrix_insert (cov->cov, entry->v1, entry->v2, entry->val1,
-					    entry->val2, entry->sum1 * entry->sum2);
+		  design_matrix_set_element (cov->cov, i, j, entry->dot_product);
 		}
 	      entry = hsh_next (cov->ca, &iter);
 	    }
-	  tmp = gsl_matrix_get (cov->cov->m, i, j);
-	  tmp -= gsl_matrix_get (cov->means->m, i, j) / gsl_matrix_get (cov->ssize->m, i, j);
-	  gsl_matrix_set (cov->cov->m, i, j, tmp);
+	  tmp = design_matrix_get_element (cov->cov, i, j);
+	  tmp -= sum_i * sum_j / design_matrix_get_element (cov->ssize, i, j);
+	  design_matrix_set_element (cov->cov, i, j, tmp);
+	  design_matrix_set_element (cov->cov, j, i, tmp);
 	} 
     }
 }
@@ -891,3 +1015,15 @@ covariance_to_design (const struct covariance_matrix *c)
     }
   return NULL;
 }
+size_t
+covariance_matrix_get_n_rows (const struct covariance_matrix *c)
+{
+  return design_matrix_get_n_rows (c->cov);
+}
+
+double 
+covariance_matrix_get_element (const struct covariance_matrix *c, size_t row, size_t col)
+{
+  return (design_matrix_get_element (c->cov, row, col));
+}
+