X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance-matrix.c;h=2a00a56d95d5c7a20a1a540bc714f539b7a9d58a;hb=5c3291dc396b795696e94f47780308fd7ace6fc4;hp=c419ce02e8b2357a868526aff9147f172b1125d7;hpb=f05f98b6b0b7e81f77657a4313fccc259eb5c506;p=pspp-builds.git

diff --git a/src/math/covariance-matrix.c b/src/math/covariance-matrix.c
index c419ce02..2a00a56d 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
@@ -56,6 +56,7 @@ struct covariance_accumulator
 struct covariance_matrix
 {
   struct design_matrix *cov;
+  struct design_matrix *ssize;
   struct hsh_table *ca;
   struct moments1 **m1;
   struct moments **m;
@@ -146,7 +147,25 @@ covariance_matrix_init (size_t n_variables,
 
   return result;
 }
-
+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.
  */
@@ -154,8 +173,8 @@ 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 void
@@ -179,6 +198,7 @@ covariance_matrix_destroy (struct covariance_matrix *cov)
 
   assert (cov != NULL);
   design_matrix_destroy (cov->cov);
+  design_matrix_destroy (cov->ssize);
   hsh_destroy (cov->ca);
   if (cov->n_pass == ONE_PASS)
     {
@@ -215,14 +235,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;
@@ -235,13 +256,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);
     }
 }
 
@@ -270,7 +291,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;
 	    }
@@ -302,9 +323,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);
     }
 }
 
@@ -350,10 +371,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;
 }
@@ -380,6 +400,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
@@ -393,39 +438,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);
 }
 
 /*
@@ -444,7 +462,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);
 }
 
@@ -456,7 +474,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))
     {
@@ -484,13 +502,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)
@@ -528,7 +549,6 @@ 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,
@@ -559,11 +579,11 @@ update_hash_entry (struct hsh_table *c,
       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.
-       */
+	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);
     }
 }
@@ -630,9 +650,6 @@ covariance_accumulate_pairwise (struct covariance_matrix *cov,
 		{
 		  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);
 		}
 	    }
 	}
@@ -702,9 +719,6 @@ covariance_accumulate_listwise (struct covariance_matrix *cov,
 	    }
 	  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);
 	}
     }
 }
@@ -722,85 +736,171 @@ covariance_matrix_accumulate (struct covariance_matrix *cov,
   cov->accumulate (cov, ccase, (const struct interaction_variable **) aux, n_intr);
 }
 
-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)
+/*
+  Return the value corresponding to subscript TARGET. If that value corresponds
+  to the origin, return NULL.
+ */
+static const union value *
+get_value_from_subscript (const struct design_matrix *dm, size_t target)
 {
-  size_t row;
-  size_t col;
+  const union value *result = NULL;
+  const struct variable *var;
   size_t i;
-  const union value *tmp_val;
-
-  assert (cov != NULL);
-
-  row = design_matrix_var_to_column (cov, v1);
-  if (var_is_alpha (v1))
+  
+  var = design_matrix_col_to_var (dm, target);
+  if (var_is_numeric (var))
     {
-      i = 0;
-      tmp_val = cat_subscript_to_value (i, v1);
-      while (compare_values_short (tmp_val, val1, v1))
+      return NULL;
+    }
+  for (i = 0; i < cat_get_n_categories (var); i++)
+    {
+      result = cat_subscript_to_value (i, var);
+      if (dm_get_exact_subscript (dm, var, result) == target)
 	{
-	  i++;
-	  tmp_val = cat_subscript_to_value (i, v1);
+	  return result;
 	}
-      row += i;
-      if (var_is_numeric (v2))
+    }
+  return NULL;
+}
+
+static bool
+is_covariance_contributor (const struct covariance_accumulator *ca, const struct design_matrix *dm,
+			   size_t i, size_t j)
+{
+  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))
+    {
+      if (var_get_dict_index (v2) == var_get_dict_index (ca->v2))
 	{
-	  col = design_matrix_var_to_column (cov, v2);
+	  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
+    }
+  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))
 	{
-	  col = design_matrix_var_to_column (cov, v2);
-	  i = 0;
-	  tmp_val = cat_subscript_to_value (i, v1);
-	  while (compare_values_short (tmp_val, val1, v1))
+	  k = dm_get_exact_subscript (dm, v1, ca->val2);
+	  if (k == i)
 	    {
-	      i++;
-	      tmp_val = cat_subscript_to_value (i, v1);
+	      k = dm_get_exact_subscript (dm, v2, ca->val1);
+	      if (k == j)
+		{
+		  return true;
+		}
 	    }
-	  col += i;
 	}
     }
-  else
+  
+  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_numeric (v2))
+      if (var_is_alpha (var))
 	{
-	  col = design_matrix_var_to_column (cov, v2);
+	  val = get_value_from_subscript (cov->cov, i);
+	  k = cat_value_find (var, val);
+	  return cat_get_category_count (k, var);
 	}
       else
 	{
-	  covariance_matrix_insert (cov, v2, v1, val2, val1, product);
+	  k = 0;
+	  while (var_get_dict_index (cov->v_variables[k]) != var_get_dict_index (var))
+	    {
+	      k++;
+	    }
+	  moments1_calculate (cov->m1[k], &n, &mean, NULL, NULL, NULL);
+	  return mean * n;
 	}
     }
 
-  gsl_matrix_set (cov->m, row, col, product);
+  return 0.0;
 }
-
-static struct design_matrix *
-covariance_accumulator_to_matrix (struct covariance_matrix *cov)
+static void
+update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_accumulator *ca)
 {
+  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))
+    {
+      var = design_matrix_col_to_var (dm, j);
+      if (var_get_dict_index (ca->v2) == var_get_dict_index (var))
+	{
+	  tmp = design_matrix_get_element (dm, i, j);
+	  tmp += ca->ssize;
+	  design_matrix_set_element (dm, i, j, tmp);
+	}
+    }
+}
+static void
+covariance_accumulator_to_matrix (struct covariance_matrix *cov)
+{
+  size_t i;
+  size_t j;
+  double sum_i = 0.0;
+  double sum_j = 0.0;
+  double tmp = 0.0;
   struct covariance_accumulator *entry;
-  struct design_matrix *result = NULL;
   struct hsh_iterator iter;
 
-  result = covariance_matrix_create (cov->n_variables, cov->v_variables);
-
+  cov->cov = covariance_matrix_create (cov->n_variables, cov->v_variables);
+  cov->ssize = 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.
-       */
-      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->ca, &iter);
     }
-  return result;
+
+  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))
+		{
+		  design_matrix_set_element (cov->cov, i, j, entry->dot_product);
+		}
+	      entry = hsh_next (cov->ca, &iter);
+	    }
+	  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);
+	} 
+    }
 }
 
 
@@ -812,7 +912,7 @@ covariance_matrix_compute (struct covariance_matrix *cov)
 {
   if (cov->n_pass == ONE_PASS)
     {
-      cov->cov = covariance_accumulator_to_matrix (cov);
+      covariance_accumulator_to_matrix (cov);
     }
 }
 
@@ -825,3 +925,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));
+}
+