Whitespace changes only
[pspp] / src / math / covariance-matrix.c
index 85998dfe49d80e1e20081f2bb503ce15bfc230e9..f929a370cf5470456ba0c33018a7f292d2e4e3a4 100644 (file)
@@ -39,80 +39,96 @@ void covariance_matrix_destroy (struct design_matrix *x)
 }
 
 /*
-  Update the covariance matrix with the new entries, assuming that V1 
-  is categorical and V2 is numeric.
+  Update the covariance matrix with the new entries, assuming that ROW
+  corresponds to a categorical variable and V2 is numeric.
  */
 static void
 covariance_update_categorical_numeric (struct design_matrix *cov, double mean,
-                         double weight, double ssize, const struct variable *v1
-                         const struct variable *v2, const union value *val1, const union value *val2)
+                         size_t row
+                         const struct variable *v2, double x, const union value *val2)
 {
-  double x;
-  size_t i;
   size_t col;
-  size_t row;
+  double tmp;
   
-  assert (var_is_alpha (v1));
   assert (var_is_numeric (v2));
 
-  row = design_matrix_var_to_column (cov, v1);  
   col = design_matrix_var_to_column (cov, v2);
-  for (i = 0; i < cat_get_n_categories (v1); i++)
+  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);
+}
+static void
+column_iterate (struct design_matrix *cov, const struct variable *v,
+               double ssize, double x, const union value *val1, size_t row)
+{
+  size_t col;
+  size_t i;
+  double y;
+  double tmp;
+  const union value *tmp_val;
+
+  col = design_matrix_var_to_column (cov, v);  
+  for (i = 0; i < cat_get_n_categories (v) - 1; i++)
     {
-      row += i;
-      x = -1.0 * cat_get_n_categories (v1) / ssize;
-      if (i == cat_value_find (v1, val1))
+      col += i;
+      y = -1.0 * cat_get_category_count (i, v) / ssize;
+      tmp_val = cat_subscript_to_value (i, v);
+      if (compare_values (tmp_val, val1, v))
        {
-         x += 1.0;
+         y += -1.0;
        }
-      assert (val2 != NULL);
-      gsl_matrix_set (cov->m, row, col, (val2->f - mean) * x * weight);
+      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);
     }
 }
-
 /*
-  Call this function in the first data pass. The central moments are
+  Call this function in the second data pass. The central moments are
   MEAN1 and MEAN2. Any categorical variables should already have their
   values summarized in in its OBS_VALS element.
  */
-void covariance_pass_one (struct design_matrix *cov, double mean1, double mean2,
-                         double weight, double ssize, const struct variable *v1, 
+void covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
+                         double ssize, const struct variable *v1, 
                          const struct variable *v2, const union value *val1, const union value *val2)
 {
   size_t row;
   size_t col;
   size_t i;
   double x;
-  double y;
+  const union value *tmp_val;
 
   if (var_is_alpha (v1))
     {
-      if (var_is_numeric (v2))
-       {
-         covariance_update_categorical_numeric (cov, mean2, weight, ssize, v1, 
-                                                v2, val1, val2);
-       }
-      else
+      row = design_matrix_var_to_column (cov, v1);
+      for (i = 0; i < cat_get_n_categories (v1) - 1; i++)
        {
-         row = design_matrix_var_to_column (cov, v1);  
-         col = design_matrix_var_to_column (cov, v2);
-         for (i = 0; i < cat_get_n_categories (v2); i++)
+         row += i;
+         x = -1.0 * cat_get_category_count (i, v1) / ssize;
+         tmp_val = cat_subscript_to_value (i, v1);
+         if (compare_values (tmp_val, val1, v1))
            {
-             col += i;
-             y = -1.0 * cat_get_n_categories (v2) / ssize;
-             if (i == cat_value_find (v2, val2))
-               {
-                 y += 1.0;
-               }
-             gsl_matrix_set (cov->m, row, col, x * y * weight);
-             gsl_matrix_set (cov->m, col, row, x * y * weight);
+             x += 1.0;
+           }
+         if (var_is_numeric (v2))
+           {
+             covariance_update_categorical_numeric (cov, mean2, row, 
+                                                    v2, x, val2);
+           }
+         else
+           {
+             column_iterate (cov, v1, ssize, x, val1, row);
+             column_iterate (cov, v2, ssize, x, val2, row);
            }
        }
     }
   else if (var_is_alpha (v2))
     {
-      covariance_update_categorical_numeric (cov, mean1, weight, ssize, v2, 
-                                            v1, val2, val1);
+      /*
+       Reverse the orders of V1, V2, etc. and put ourselves back
+       in the previous IF scope.
+       */
+      covariance_pass_two (cov, mean2, mean1, ssize, v2, v1, val2, val1);
     }
   else
     {
@@ -121,7 +137,8 @@ void covariance_pass_one (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) * weight;
+      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);
     }