Set categoricals to NULL in one-pass covariance; avoid freeing NULL in categoricals_d...

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

diff --git a/src/math/covariance.c b/src/math/covariance.c
index d63159b96d5c3244b85ccd51436170b7e20ec5e9..1b5a238558ac773fa30b7016bc83bf521b406e1c 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,10 +193,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;
 }
@@ -232,7 +233,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;
@@ -263,17 +297,16 @@ covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
 
   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->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 +314,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,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);
 
@@ -316,10 +352,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);
@@ -329,20 +396,20 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 
   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->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 +418,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 +429,6 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
            {
              cov->cm [idx] += ss;
            }
-
        }
     }
 
@@ -391,7 +457,7 @@ covariance_accumulate (struct covariance *cov, const struct ccase *c)
     {
       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->dim; ++j)
@@ -400,7 +466,7 @@ covariance_accumulate (struct covariance *cov, const struct ccase *c)
          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);