New function covariance_calculate_unnormalized

[pspp-builds.git] / src / math / covariance.c

diff --git a/src/math/covariance.c b/src/math/covariance.c
index dbcf4f65efb8d8b5c586c3e90afb61b73327da33..dc3166923d7ddff670d903274caef1dfe548beb4 100644 (file)
--- a/src/math/covariance.c
+++ b/src/math/covariance.c
@@ -157,6 +157,7 @@ covariance_1pass_create (size_t n_vars, const struct variable **vars,
   cov->n_cm = (n_vars * (n_vars - 1)  ) / 2;
 
   cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
+  cov->categoricals = NULL;
 
   return cov;
 }
@@ -192,7 +193,7 @@ 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, exclude);
@@ -265,7 +266,7 @@ get_val (const struct covariance *cov, int i, const struct ccase *c)
   return categoricals_get_binary_by_subscript (cov->categoricals, i - cov->n_vars, c);
 }
 
-static void
+void
 dump_matrix (const gsl_matrix *m)
 {
   size_t i, j;
@@ -334,10 +335,13 @@ 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;
 
-      cov->dim = cov->n_vars + categoricals_total (cov->categoricals);
+      cov->dim = cov->n_vars +
+       categoricals_total (cov->categoricals) - categoricals_get_n_variables (cov->categoricals);
+
       cov->n_cm = (cov->dim * (cov->dim - 1)  ) / 2;
       cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
 
@@ -350,10 +354,39 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 
       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);
@@ -396,7 +429,6 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
            {
              cov->cm [idx] += ss;
            }
-
        }
     }
 
@@ -558,7 +590,6 @@ covariance_calculate_single_pass (struct covariance *cov)
 }
 
 
-
 /* 
    Return a pointer to gsl_matrix containing the pairwise covariances.
    The matrix remains owned by the COV object, and must not be freed.
@@ -582,6 +613,86 @@ covariance_calculate (struct covariance *cov)
     }
 }
 
+/*
+  Covariance computed without dividing by the sample size.
+ */
+static const gsl_matrix *
+covariance_calculate_double_pass_unnormalized (struct covariance *cov)
+{
+  size_t i, j;
+  for (i = 0 ; i < cov->dim; ++i)
+    {
+      for (j = 0 ; j < cov->dim; ++j)
+       {
+         int idx;
+         double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
+
+         idx = cm_idx (cov, i, j);
+         if ( idx >= 0)
+           {
+             x = &cov->cm [idx];
+           }
+       }
+    }
+
+  return  cm_to_gsl (cov);
+}
+
+static const gsl_matrix *
+covariance_calculate_single_pass_unnormalized (struct covariance *cov)
+{
+  size_t i, j;
+  size_t m;
+
+  for (i = 0 ; i < cov->dim; ++i)
+    {
+      for (j = 0 ; j < cov->dim; ++j)
+       {
+         double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
+         *x -= pow2 (gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
+           / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
+       }
+    }
+  for ( j = 0 ; j < cov->dim - 1; ++j)
+    {
+      for (i = j + 1 ; i < cov->dim; ++i)
+       {
+         double *x = &cov->cm [cm_idx (cov, i, j)];
+         
+         *x -=
+           gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j) 
+           *
+           gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i) 
+         / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
+       }
+    }
+
+  return cm_to_gsl (cov);
+}
+
+
+/* 
+   Return a pointer to gsl_matrix containing the pairwise covariances.
+   The matrix remains owned by the COV object, and must not be freed.
+   Call this function only after all data have been accumulated.
+*/
+const gsl_matrix *
+covariance_calculate_unnormalized (struct covariance *cov)
+{
+  assert ( cov->state > 0 );
+
+  switch (cov->passes)
+    {
+    case 1:
+      return covariance_calculate_single_pass_unnormalized (cov);  
+      break;
+    case 2:
+      return covariance_calculate_double_pass_unnormalized (cov);  
+      break;
+    default:
+      NOT_REACHED ();
+    }
+}
 
 
 
@@ -590,7 +701,7 @@ void
 covariance_destroy (struct covariance *cov)
 {
   size_t i;
-  free (cov->vars);
+
   categoricals_destroy (cov->categoricals);
 
   for (i = 0; i < n_MOMENTS; ++i)