X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance-matrix.c;h=95895eaeaca3c8337da4e55da49dc7f78f624419;hb=537fdeb3702c011e05d7826a8d556a7beeba2605;hp=972a3ff00726a12200bbdab00086180b0ec100df;hpb=428b694be3fe5ab2b4ba430f3a1d12d60a8e510c;p=pspp-builds.git

diff --git a/src/math/covariance-matrix.c b/src/math/covariance-matrix.c
index 972a3ff0..95895eae 100644
--- a/src/math/covariance-matrix.c
+++ b/src/math/covariance-matrix.c
@@ -56,6 +56,8 @@ struct covariance_accumulator
 struct covariance_matrix
 {
   struct design_matrix *cov;
+  struct design_matrix *ssize;
+  struct design_matrix *sums;
   struct hsh_table *ca;
   struct moments1 **m1;
   struct moments **m;
@@ -146,7 +148,25 @@ covariance_matrix_init (size_t n_variables,
 
   return result;
 }
-
+static size_t 
+get_n_rows (size_t n_variables, size_t *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 +174,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 +199,8 @@ covariance_matrix_destroy (struct covariance_matrix *cov)
 
   assert (cov != NULL);
   design_matrix_destroy (cov->cov);
+  design_matrix_destroy (cov->ssize);
+  design_matrix_destroy (cov->sums);
   hsh_destroy (cov->ca);
   if (cov->n_pass == ONE_PASS)
     {
@@ -215,9 +237,9 @@ 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,
@@ -239,9 +261,9 @@ column_iterate (struct design_matrix *cov, const struct variable *v,
 	{
 	  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);
     }
 }
 
@@ -302,9 +324,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);
     }
 }
 
@@ -767,8 +789,6 @@ covariance_matrix_insert (struct design_matrix *cov,
 {
   size_t row;
   size_t col;
-  size_t i;
-  const union value *tmp_val;
 
   assert (cov != NULL);
 
@@ -776,7 +796,7 @@ covariance_matrix_insert (struct design_matrix *cov,
   col = get_exact_subscript (cov, v2, val2);
   if (row != -1u && col != -1u)
     {
-      gsl_matrix_set (cov->m, row, col, product);
+      design_matrix_set_element (cov, row, col, product);
     }
 }
 
@@ -809,41 +829,94 @@ is_covariance_contributor (const struct covariance_accumulator *ca, const struct
     }
   return false;
 }
-  
-static struct design_matrix *
+static double
+get_sum (const struct covariance_matrix *cov, size_t i)
+{
+  size_t k;
+  const struct variable *var;
+  const union value *val = NULL;
+  struct covariance_accumulator ca;
+  struct covariance_accumulator *c;
+
+  assert ( cov != NULL);
+  var = design_matrix_col_to_var (cov->cov, i);
+  if (var != NULL)
+    {
+      if (var_is_alpha (var))
+	{
+	  k = design_matrix_var_to_column (cov->cov, var);
+	  i -= k;
+	  val = cat_subscript_to_value (i, var);
+	}
+      ca.v1 = var;
+      ca.v2 = var;
+      ca.val1 = val;
+      ca.val2 = val;
+      c = (struct covariance_accumulator *) hsh_find (cov->ca, &ca);
+      if (c != NULL)
+	{
+	  return c->sum1;
+	}
+    }
+  return 0.0;
+}
+static void
+update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_accumulator *ca)
+{
+  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 tmp;
+  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);
-
-  for (i = 0; i < design_matrix_get_n_cols (result); i++)
+  cov->cov = covariance_matrix_create (cov->n_variables, cov->v_variables);
+  cov->ssize = covariance_matrix_create (cov->n_variables, cov->v_variables);
+  cov->sums = covariance_matrix_create (cov->n_variables, cov->v_variables);
+  for (i = 0; i < design_matrix_get_n_cols (cov->cov); i++)
     {
-      for (j = i; j < design_matrix_get_n_cols (result); j++)
+      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);
-	  
+	  design_matrix_set_element (cov->sums, i, j, sum_i);	  
 	  while (entry != NULL)
 	    {
+	      update_ssize (cov->ssize, i, j, entry);
 	      /*
 		We compute the centered, un-normalized covariance matrix.
 	      */
-	      if (is_covariance_contributor (entry, result, i, j))
+	      if (is_covariance_contributor (entry, cov->cov, i, j))
 		{
-		  tmp = entry->dot_product - entry->sum1 * entry->sum2 / entry->ssize;		  
-		  covariance_matrix_insert (result, entry->v1, entry->v2, entry->val1,
-					    entry->val2, tmp);
+		  covariance_matrix_insert (cov->cov, entry->v1, entry->v2, entry->val1,
+					    entry->val2, 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);
 	} 
     }
-  return result;
 }
 
 
@@ -855,7 +928,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);
     }
 }
 
@@ -868,3 +941,10 @@ covariance_to_design (const struct covariance_matrix *c)
     }
   return NULL;
 }
+
+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));
+}
+