X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance.c;h=c247148ee9f88d64a29edf36e087d79ab21e95cd;hb=5055ec22100c84153e76662bae23063408a6a935;hp=e9f29a299f9b94b31ec196d8eb3a7424bc672e2e;hpb=8539828ab3d09acbc9e0bf7242db21d472f24e83;p=pspp-builds.git

diff --git a/src/math/covariance.c b/src/math/covariance.c
index e9f29a29..c247148e 100644
--- a/src/math/covariance.c
+++ b/src/math/covariance.c
@@ -192,10 +192,10 @@ covariance_2pass_create (size_t n_vars, const struct variable **vars,
 
   cov->exclude = exclude;
 
-  cov->n_cm = - 1;
+  cov->n_cm = -1;
   cov->cm = NULL;
 
-  cov->categoricals = categoricals_create (catvars, n_catvars, wv);
+  cov->categoricals = categoricals_create (catvars, n_catvars, wv, exclude);
 
   return cov;
 }
@@ -210,16 +210,16 @@ static int
 cm_idx (const struct covariance *cov, int i, int j)
 {
   int as;
-  const int n2j = cov->n_vars - 2 - j;
-  const int nj = cov->n_vars - 2 ;
+  const int n2j = cov->dim - 2 - j;
+  const int nj = cov->dim - 2 ;
   
   assert (i >= 0);
-  assert (j < cov->n_vars);
+  assert (j < cov->dim);
 
   if ( i == 0)
     return -1;
 
-  if (j >= cov->n_vars - 1)
+  if (j >= cov->dim - 1)
     return -1;
 
   if ( i <= j) 
@@ -232,7 +232,40 @@ cm_idx (const struct covariance *cov, int i, int j)
   return i - 1 + as;
 }
 
-static void
+
+/*
+  Returns true iff the variable corresponding to the Ith element of the covariance matrix 
+   has a missing value for case C
+*/
+static bool
+is_missing (const struct covariance *cov, int i, const struct ccase *c)
+{
+  const struct variable *var = i < cov->n_vars ?
+    cov->vars[i] : 
+    categoricals_get_variable_by_subscript (cov->categoricals, i - cov->n_vars);
+
+  const union value *val = case_data (c, var);
+
+  return var_is_value_missing (var, val, cov->exclude);
+}
+
+
+static double
+get_val (const struct covariance *cov, int i, const struct ccase *c)
+{
+  if ( i < cov->n_vars)
+    {
+      const struct variable *var = cov->vars[i];
+
+      const union value *val = case_data (c, var);
+
+      return val->f;
+    }
+
+  return categoricals_get_binary_by_subscript (cov->categoricals, i - cov->n_vars, c);
+}
+
+void
 dump_matrix (const gsl_matrix *m)
 {
   size_t i, j;
@@ -261,19 +294,18 @@ covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
 
   categoricals_update (cov->categoricals, c);
 
-  for (i = 0 ; i < cov->n_vars; ++i)
+  for (i = 0 ; i < cov->dim; ++i)
     {
-      const union value *val1 = case_data (c, cov->vars[i]);
+      double v1 = get_val (cov, i, c);
 
-      if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
+      if ( is_missing (cov, i, c))
 	continue;
 
-      for (j = 0 ; j < cov->n_vars; ++j)
+      for (j = 0 ; j < cov->dim; ++j)
 	{
 	  double pwr = 1.0;
-	  const union value *val2 = case_data (c, cov->vars[j]);
 
-	  if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
+	  if ( is_missing (cov, j, c))
 	    continue;
 
 	  for (m = 0 ; m <= MOMENT_MEAN; ++m)
@@ -281,7 +313,7 @@ covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
 	      double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 
 	      *x += pwr * weight;
-	      pwr *= val1->f;
+	      pwr *= v1;
 	    }
 	}
     }
@@ -302,6 +334,7 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 
   if (! cov->pass_two_first_case_seen)
     {
+      size_t m;
       assert (cov->state == 1);
       cov->state = 2;
 
@@ -316,10 +349,41 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 	  cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
 	}
 
+      categoricals_done (cov->categoricals);
+
+      /* Populate the moments matrices with the categorical value elements */
+      for (i = cov->n_vars; i < cov->dim; ++i)
+	{
+	  for (j = 0 ; j < cov->dim ; ++j) /* FIXME: This is WRONG !!! */
+	    {
+	      double w = categoricals_get_weight_by_subscript (cov->categoricals, i - cov->n_vars);
+
+	      gsl_matrix_set (cov->moments[MOMENT_NONE], i, j, w);
+
+	      w = categoricals_get_sum_by_subscript (cov->categoricals, i - cov->n_vars);
+
+	      gsl_matrix_set (cov->moments[MOMENT_MEAN], i, j, w);
+	    }
+	}
+
+      /* FIXME: This is WRONG!!  It must be fixed to properly handle missing values.  For
+       now it assumes there are none */
+      for (m = 0 ; m < n_MOMENTS; ++m)
+	{
+	  for (i = 0 ; i < cov->dim ; ++i)
+	    {
+	      double x = gsl_matrix_get (cov->moments[m], i, cov->n_vars -1);
+	      for (j = cov->n_vars; j < cov->dim; ++j)
+		{
+		  gsl_matrix_set (cov->moments[m], i, j, x);
+		}
+	    }
+	}
+
       /* Divide the means by the number of samples */
-      for (i = 0; i < cov->n_vars; ++i)
+      for (i = 0; i < cov->dim; ++i)
 	{
-	  for (j = 0; j < cov->n_vars; ++j)
+	  for (j = 0; j < cov->dim; ++j)
 	    {
 	      double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
 	      *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
@@ -327,22 +391,22 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 	}
     }
 
-  for (i = 0 ; i < cov->n_vars; ++i)
+  for (i = 0 ; i < cov->dim; ++i)
     {
-      const union value *val1 = case_data (c, cov->vars[i]);
+      double v1 = get_val (cov, i, c);
 
-      if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
+      if ( is_missing (cov, i, c))
 	continue;
 
-      for (j = 0 ; j < cov->n_vars; ++j)
+      for (j = 0 ; j < cov->dim; ++j)
 	{
 	  int idx;
 	  double ss ;
-	  const union value *val2 = case_data (c, cov->vars[j]);
+	  double v2 = get_val (cov, j, c);
 
-	  const double s = pow2 (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
+	  const double s = pow2 (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
 
-	  if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
+	  if ( is_missing (cov, j, c))
 	    continue;
 
 	  {
@@ -351,9 +415,9 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 	  }
 
 	  ss = 
-	    (val1->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
+	    (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 	    * 
-	    (val2->f - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
+	    (v2 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
 	    * weight
 	    ;
 
@@ -362,7 +426,6 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 	    {
 	      cov->cm [idx] += ss;
 	    }
-
 	}
     }
 
@@ -387,20 +450,20 @@ covariance_accumulate (struct covariance *cov, const struct ccase *c)
       cov->state = 1;
     }
 
-  for (i = 0 ; i < cov->n_vars; ++i)
+  for (i = 0 ; i < cov->dim; ++i)
     {
       const union value *val1 = case_data (c, cov->vars[i]);
 
-      if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
+      if ( is_missing (cov, i, c))
 	continue;
 
-      for (j = 0 ; j < cov->n_vars; ++j)
+      for (j = 0 ; j < cov->dim; ++j)
 	{
 	  double pwr = 1.0;
 	  int idx;
 	  const union value *val2 = case_data (c, cov->vars[j]);
 
-	  if ( var_is_value_missing (cov->vars[j], val2, cov->exclude))
+	  if ( is_missing (cov, j, c))
 	    continue;
 
 	  idx = cm_idx (cov, i, j);
@@ -431,12 +494,12 @@ static gsl_matrix *
 cm_to_gsl (struct covariance *cov)
 {
   int i, j;
-  gsl_matrix *m = gsl_matrix_calloc (cov->n_vars, cov->n_vars);
+  gsl_matrix *m = gsl_matrix_calloc (cov->dim, cov->dim);
 
   /* Copy the non-diagonal elements from cov->cm */
-  for ( j = 0 ; j < cov->n_vars - 1; ++j)
+  for ( j = 0 ; j < cov->dim - 1; ++j)
     {
-      for (i = j+1 ; i < cov->n_vars; ++i)
+      for (i = j+1 ; i < cov->dim; ++i)
 	{
 	  double x = cov->cm [cm_idx (cov, i, j)];
 	  gsl_matrix_set (m, i, j, x);
@@ -445,7 +508,7 @@ cm_to_gsl (struct covariance *cov)
     }
 
   /* Copy the diagonal elements from cov->moments[2] */
-  for (j = 0 ; j < cov->n_vars ; ++j)
+  for (j = 0 ; j < cov->dim ; ++j)
     {
       double sigma = gsl_matrix_get (cov->moments[2], j, j);
       gsl_matrix_set (m, j, j, sigma);
@@ -459,9 +522,9 @@ static const gsl_matrix *
 covariance_calculate_double_pass (struct covariance *cov)
 {
   size_t i, j;
-  for (i = 0 ; i < cov->n_vars; ++i)
+  for (i = 0 ; i < cov->dim; ++i)
     {
-      for (j = 0 ; j < cov->n_vars; ++j)
+      for (j = 0 ; j < cov->dim; ++j)
 	{
 	  int idx;
 	  double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
@@ -490,9 +553,9 @@ covariance_calculate_single_pass (struct covariance *cov)
       /* Divide the moments by the number of samples */
       if ( m > 0)
 	{
-	  for (i = 0 ; i < cov->n_vars; ++i)
+	  for (i = 0 ; i < cov->dim; ++i)
 	    {
-	      for (j = 0 ; j < cov->n_vars; ++j)
+	      for (j = 0 ; j < cov->dim; ++j)
 		{
 		  double *x = gsl_matrix_ptr (cov->moments[m], i, j);
 		  *x /= gsl_matrix_get (cov->moments[0], i, j);
@@ -505,9 +568,9 @@ covariance_calculate_single_pass (struct covariance *cov)
     }
 
   /* Centre the moments */
-  for ( j = 0 ; j < cov->n_vars - 1; ++j)
+  for ( j = 0 ; j < cov->dim - 1; ++j)
     {
-      for (i = j + 1 ; i < cov->n_vars; ++i)
+      for (i = j + 1 ; i < cov->dim; ++i)
 	{
 	  double *x = &cov->cm [cm_idx (cov, i, j)];
 	  
@@ -557,6 +620,7 @@ covariance_destroy (struct covariance *cov)
 {
   size_t i;
   free (cov->vars);
+  categoricals_destroy (cov->categoricals);
 
   for (i = 0; i < n_MOMENTS; ++i)
     gsl_matrix_free (cov->moments[i]);