Use cov->dim instead of cov->n_vars where appropriate

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

diff --git a/src/math/covariance.c b/src/math/covariance.c
index 1d908b3d1754ab9980f5717e968fb478598a14b5..d63159b96d5c3244b85ccd51436170b7e20ec5e9 100644 (file)
--- a/src/math/covariance.c
+++ b/src/math/covariance.c
@@ -29,6 +29,41 @@
 #define n_MOMENTS (MOMENT_VARIANCE + 1)
 
 
+/* Create a new matrix of NEW_SIZE x NEW_SIZE and copy the elements of
+   matrix IN into it.  IN must be a square matrix, and in normal usage
+   it will be smaller than NEW_SIZE.
+   IN is destroyed by this function.  The return value must be destroyed
+   when no longer required.
+*/
+static gsl_matrix *
+resize_matrix (gsl_matrix *in, size_t new_size)
+{
+  size_t i, j;
+
+  gsl_matrix *out = NULL;
+
+  assert (in->size1 == in->size2);
+
+  if (new_size <= in->size1)
+    return in;
+
+  out = gsl_matrix_calloc (new_size, new_size);
+
+  for (i = 0; i < in->size1; ++i)
+    {
+      for (j = 0; j < in->size2; ++j)
+       {
+         double x = gsl_matrix_get (in, i, j);
+
+         gsl_matrix_set (out, i, j, x);
+       }
+    }
+    
+  gsl_matrix_free (in);
+
+  return out;
+}
+
 struct covariance
 {
   /* The variables for which the covariance matrix is to be calculated. */
@@ -175,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) 
@@ -226,14 +261,14 @@ 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]);
 
       if ( var_is_value_missing (cov->vars[i], val1, cov->exclude))
        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]);
@@ -274,6 +309,13 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
       cov->n_cm = (cov->dim * (cov->dim - 1)  ) / 2;
       cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
 
+      /* Grow the moment matrices so that they're large enough to accommodate the
+        categorical elements */
+      for (i = 0; i < n_MOMENTS; ++i)
+       {
+         cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
+       }
+
       /* Divide the means by the number of samples */
       for (i = 0; i < cov->n_vars; ++i)
        {
@@ -281,18 +323,18 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
            {
              double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
              *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
-           }
+           }
        }
     }
 
-  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))
        continue;
 
-      for (j = 0 ; j < cov->n_vars; ++j)
+      for (j = 0 ; j < cov->dim; ++j)
        {
          int idx;
          double ss ;
@@ -345,14 +387,14 @@ 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))
        continue;
 
-      for (j = 0 ; j < cov->n_vars; ++j)
+      for (j = 0 ; j < cov->dim; ++j)
        {
          double pwr = 1.0;
          int idx;
@@ -389,12 +431,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);
@@ -403,7 +445,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);
@@ -417,9 +459,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);
@@ -448,9 +490,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);
@@ -463,9 +505,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)];
          
@@ -515,6 +557,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]);