Change "union value" to dynamically allocate long strings.
[pspp-builds.git] / src / math / covariance-matrix.c
index bb17d99aa9c616d49aa1f7d3f3f62201d3508710..2a00a56d95d5c7a20a1a540bc714f539b7a9d58a 100644 (file)
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2008 Free Software Foundation, Inc.
+   Copyright (C) 2008, 2009 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -57,7 +57,6 @@ struct covariance_matrix
 {
   struct design_matrix *cov;
   struct design_matrix *ssize;
-  struct design_matrix *sums;
   struct hsh_table *ca;
   struct moments1 **m1;
   struct moments **m;
@@ -148,7 +147,25 @@ covariance_matrix_init (size_t n_variables,
 
   return result;
 }
-
+static size_t 
+get_n_rows (size_t n_variables, const struct variable *v_variables[])
+{
+  size_t i;
+  size_t result = 0;
+  for (i = 0; i < n_variables; i++)
+    {
+      if (var_is_numeric (v_variables[i]))
+       {
+         result++;
+       }
+      else if (var_is_alpha (v_variables[i]))
+       {
+         size_t n_categories = cat_get_n_categories (v_variables[i]);
+         result += n_categories - 1;
+       }
+    }
+  return result;
+}
 /*
   The covariances are stored in a DESIGN_MATRIX structure.
  */
@@ -156,8 +173,8 @@ struct design_matrix *
 covariance_matrix_create (size_t n_variables,
                          const struct variable *v_variables[])
 {
-  return design_matrix_create (n_variables, v_variables,
-                              (size_t) n_variables);
+  size_t n_rows = get_n_rows (n_variables, v_variables);
+  return design_matrix_create (n_variables, v_variables, n_rows);
 }
 
 static void
@@ -182,7 +199,6 @@ covariance_matrix_destroy (struct covariance_matrix *cov)
   assert (cov != NULL);
   design_matrix_destroy (cov->cov);
   design_matrix_destroy (cov->ssize);
-  design_matrix_destroy (cov->sums);
   hsh_destroy (cov->ca);
   if (cov->n_pass == ONE_PASS)
     {
@@ -227,6 +243,7 @@ static void
 column_iterate (struct design_matrix *cov, const struct variable *v,
                double ssize, double x, const union value *val1, size_t row)
 {
+  int width = var_get_width (v);
   size_t col;
   size_t i;
   double y;
@@ -239,7 +256,7 @@ column_iterate (struct design_matrix *cov, const struct variable *v,
       col += i;
       y = -1.0 * cat_get_category_count (i, v) / ssize;
       tmp_val = cat_subscript_to_value (i, v);
-      if (!compare_values_short (tmp_val, val1, v))
+      if (!value_equal (tmp_val, val1, width))
        {
          y += -1.0;
        }
@@ -274,7 +291,7 @@ covariance_pass_two (struct design_matrix *cov, double mean1, double mean2,
          row += i;
          x = -1.0 * cat_get_category_count (i, v1) / ssize;
          tmp_val = cat_subscript_to_value (i, v1);
-         if (!compare_values_short (tmp_val, val1, v1))
+         if (!value_equal (tmp_val, val1, var_get_width (v1)))
            {
              x += 1.0;
            }
@@ -354,10 +371,9 @@ covariance_accumulator_hash (const void *h, const void *aux)
     }
   if (var_is_alpha (v_max) && var_is_alpha (v_min))
     {
-      unsigned tmp = hsh_hash_bytes (val_max, var_get_width (v_max));
-      tmp ^= hsh_hash_bytes (val_min, var_get_width (v_min));
-      tmp += *n_vars * (*n_vars + 1 + idx_max) + idx_min;
-      return (size_t) tmp;
+      unsigned hash = value_hash (val_max, var_get_width (v_max), 0);
+      hash = value_hash (val_min, var_get_width (v_min), hash);
+      return hash_int (*n_vars * (*n_vars + 1 + idx_max) + idx_min, hash);
     }
   return -1u;
 }
@@ -384,6 +400,31 @@ covariance_accumulator_free (void *c_, const void *aux UNUSED)
   free (c);
 }
 
+static int 
+ordered_match_nodes (const struct covariance_accumulator *c, const struct variable *v1,
+                    const struct variable *v2, const union value *val1, const union value *val2)
+{
+  size_t result;
+  size_t m;
+
+  result = var_get_dict_index (v1) ^ var_get_dict_index (c->v1);
+  m = var_get_dict_index (v2) ^ var_get_dict_index (c->v2);
+  result = result|m;
+  if (var_is_alpha (v1))
+    {
+      result |= value_compare_3way (val1, c->val1, var_get_width (v1));
+      if (var_is_alpha (v2))
+       {
+         result |= value_compare_3way (val2, c->val2, var_get_width (v2));
+       }
+    }
+  else if (var_is_alpha (v2))
+    {
+      result |= value_compare_3way (val2, c->val2, var_get_width (v2));
+    }
+  return result;
+}
+  
 /*
   Hash comparison. Returns 0 for a match, or a non-zero int
   otherwise. The sign of a non-zero return value *should* indicate the
@@ -397,39 +438,12 @@ match_nodes (const struct covariance_accumulator *c,
             const struct variable *v1, const struct variable *v2,
             const union value *val1, const union value *val2)
 {
-  if (var_get_dict_index (v1) == var_get_dict_index (c->v1))
-    if (var_get_dict_index (v2) == var_get_dict_index (c->v2))
-      {
-       if (var_is_numeric (v1) && var_is_numeric (v2))
-         {
-           return 0;
-         }
-       if (var_is_numeric (v1) && var_is_alpha (v2))
-         {
-           if (!compare_values_short (val2, c->val2, v2))
-             {
-               return 0;
-             }
-         }
-       if (var_is_alpha (v1) && var_is_numeric (v2))
-         {
-           if (!compare_values_short (val1, c->val1, v1))
-             {
-               return 0;
-             }
-         }
-       if (var_is_alpha (v1) && var_is_alpha (v2))
-         {
-           if (!compare_values_short (val1, c->val1, v1))
-             {
-               if (!compare_values_short (val2, c->val2, v2))
-                 {
-                   return 0;
-                 }
-             }
-         }
-      }
-  return 1;
+  size_t n;
+  size_t m;
+
+  n = ordered_match_nodes (c, v1, v2, val1, val2);
+  m = ordered_match_nodes (c, v2, v1, val2, val1);
+  return (n & m);
 }
 
 /*
@@ -448,7 +462,7 @@ covariance_accumulator_compare (const void *a1_, const void *a2_,
 
   if (a1 == NULL || a2 == NULL)
     return 1;
-
+  
   return match_nodes (a1, a2->v1, a2->v2, a2->val1, a2->val2);
 }
 
@@ -460,7 +474,7 @@ hash_numeric_alpha (const struct variable *v1, const struct variable *v2,
   if (var_is_numeric (v1) && var_is_alpha (v2))
     {
       result = n_vars * ((n_vars + 1) + var_get_dict_index (v1))
-       + var_get_dict_index (v2) + hsh_hash_string (val->s);
+       + var_get_dict_index (v2) + value_hash (val, var_get_width (v2), 0);
     }
   else if (var_is_alpha (v1) && var_is_numeric (v2))
     {
@@ -488,13 +502,16 @@ update_product (const struct variable *v1, const struct variable *v2,
     }
   if (var_is_numeric (v1) && var_is_alpha (v2))
     {
-      return (val1->f);
+      return val1->f;
     }
   if (var_is_numeric (v2) && var_is_alpha (v1))
     {
-      update_product (v2, v1, val2, val1);
+      return val2->f;
+    }
+  else
+    {
+      return 0.0;
     }
-  return 0.0;
 }
 static double
 update_sum (const struct variable *var, const union value *val, double weight)
@@ -554,7 +571,7 @@ update_hash_entry (struct hsh_table *c,
   ca->sum2 = update_sum (ca->v2, ca->val2, iv_f2);
   ca->ssize = 1.0;
   new_entry = hsh_insert (c, ca);
-  
+
   if (new_entry != NULL)
     {
       new_entry->dot_product += ca->dot_product;
@@ -633,9 +650,6 @@ covariance_accumulate_pairwise (struct covariance_matrix *cov,
                {
                  update_hash_entry (cov->ca, v_variables[i], v_variables[j],
                                     val1, val2, i_val1, i_val2);
-                 if (j != i)
-                   update_hash_entry (cov->ca, v_variables[j],
-                                      v_variables[i], val2, val1, i_val2, i_val1);
                }
            }
        }
@@ -705,9 +719,6 @@ covariance_accumulate_listwise (struct covariance_matrix *cov,
            }
          update_hash_entry (cov->ca, v_variables[i], v_variables[j],
                             val1, val2, i_val1, i_val2);
-         if (j != i)
-           update_hash_entry (cov->ca, v_variables[j], v_variables[i],
-                              val2, val1, i_val2, i_val1);
        }
     }
 }
@@ -724,65 +735,34 @@ covariance_matrix_accumulate (struct covariance_matrix *cov,
 {
   cov->accumulate (cov, ccase, (const struct interaction_variable **) aux, n_intr);
 }
-/*
-  If VAR is categorical with d categories, its first category should
-  correspond to the origin in d-dimensional Euclidean space.
- */
-static bool
-is_origin (const struct variable *var, const union value *val)
-{
-  if (cat_value_find (var, val) == 0)
-    {
-      return true;
-    }
-  return false;
-}
 
 /*
-  Return the subscript of the column of the design matrix
-  corresponding to VAL. If VAR is categorical with d categories, its
-  first category should correspond to the origin in d-dimensional
-  Euclidean space, so there is no subscript for this value.
+  Return the value corresponding to subscript TARGET. If that value corresponds
+  to the origin, return NULL.
  */
-static size_t
-get_exact_subscript (const struct design_matrix *dm, const struct variable *var,
-                    const union value *val)
+static const union value *
+get_value_from_subscript (const struct design_matrix *dm, size_t target)
 {
-  size_t result;
-
-  if (is_origin (var, val))
+  const union value *result = NULL;
+  const struct variable *var;
+  size_t i;
+  
+  var = design_matrix_col_to_var (dm, target);
+  if (var_is_numeric (var))
     {
-      return -1u;
+      return NULL;
     }
-
-  result = design_matrix_var_to_column (dm, var);
-  if (var_is_alpha (var))
+  for (i = 0; i < cat_get_n_categories (var); i++)
     {
-      result += cat_value_find (var, val) - 1;
-    }
-  return result;
-}
-
-static void
-covariance_matrix_insert (struct design_matrix *cov,
-                         const struct variable *v1,
-                         const struct variable *v2, const union value *val1,
-                         const union value *val2, double product)
-{
-  size_t row;
-  size_t col;
-
-  assert (cov != NULL);
-
-  row = get_exact_subscript (cov, v1, val1);
-  col = get_exact_subscript (cov, v2, val2);
-  if (row != -1u && col != -1u)
-    {
-      design_matrix_set_element (cov, row, col, product);
+      result = cat_subscript_to_value (i, var);
+      if (dm_get_exact_subscript (dm, var, result) == target)
+       {
+         return result;
+       }
     }
+  return NULL;
 }
 
-
 static bool
 is_covariance_contributor (const struct covariance_accumulator *ca, const struct design_matrix *dm,
                           size_t i, size_t j)
@@ -790,18 +770,18 @@ is_covariance_contributor (const struct covariance_accumulator *ca, const struct
   size_t k;
   const struct variable *v1;
   const struct variable *v2;
-
+  
   assert (dm != NULL);
   v1 = design_matrix_col_to_var (dm, i);
+  v2 = design_matrix_col_to_var (dm, j);
   if (var_get_dict_index (v1) == var_get_dict_index(ca->v1))
     {
-      v2 = design_matrix_col_to_var (dm, j);
       if (var_get_dict_index (v2) == var_get_dict_index (ca->v2))
        {
-         k = get_exact_subscript (dm, v1, ca->val1);
+         k = dm_get_exact_subscript (dm, v1, ca->val1);
          if (k == i)
            {
-             k = get_exact_subscript (dm, v2, ca->val2);
+             k = dm_get_exact_subscript (dm, v2, ca->val2);
              if (k == j)
                {
                  return true;
@@ -809,16 +789,32 @@ is_covariance_contributor (const struct covariance_accumulator *ca, const struct
            }
        }
     }
+  else if (var_get_dict_index (v1) == var_get_dict_index (ca->v2))
+    {
+      if (var_get_dict_index (v2) == var_get_dict_index (ca->v1))
+       {
+         k = dm_get_exact_subscript (dm, v1, ca->val2);
+         if (k == i)
+           {
+             k = dm_get_exact_subscript (dm, v2, ca->val1);
+             if (k == j)
+               {
+                 return true;
+               }
+           }
+       }
+    }
+  
   return false;
 }
 static double
 get_sum (const struct covariance_matrix *cov, size_t i)
 {
   size_t k;
+  double mean;
+  double n;
   const struct variable *var;
   const union value *val = NULL;
-  struct covariance_accumulator ca;
-  struct covariance_accumulator *c;
 
   assert ( cov != NULL);
   var = design_matrix_col_to_var (cov->cov, i);
@@ -826,26 +822,28 @@ get_sum (const struct covariance_matrix *cov, size_t i)
     {
       if (var_is_alpha (var))
        {
-         k = design_matrix_var_to_column (cov->cov, var);
-         i -= k;
-         val = cat_subscript_to_value (i, var);
+         val = get_value_from_subscript (cov->cov, i);
+         k = cat_value_find (var, val);
+         return cat_get_category_count (k, var);
        }
-      ca.v1 = var;
-      ca.v2 = var;
-      ca.val1 = val;
-      ca.val2 = val;
-      c = (struct covariance_accumulator *) hsh_find (cov->ca, &ca);
-      if (c != NULL)
+      else
        {
-         return c->sum1;
+         k = 0;
+         while (var_get_dict_index (cov->v_variables[k]) != var_get_dict_index (var))
+           {
+             k++;
+           }
+         moments1_calculate (cov->m1[k], &n, &mean, NULL, NULL, NULL);
+         return mean * n;
        }
     }
+
   return 0.0;
 }
 static void
 update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_accumulator *ca)
 {
-  struct variable *var;
+  const struct variable *var;
   double tmp;
   var = design_matrix_col_to_var (dm, i);
   if (var_get_dict_index (ca->v1) == var_get_dict_index (var))
@@ -872,7 +870,12 @@ covariance_accumulator_to_matrix (struct covariance_matrix *cov)
 
   cov->cov = covariance_matrix_create (cov->n_variables, cov->v_variables);
   cov->ssize = covariance_matrix_create (cov->n_variables, cov->v_variables);
-  cov->sums = covariance_matrix_create (cov->n_variables, cov->v_variables);
+  entry = hsh_first (cov->ca, &iter);
+  while (entry != NULL)
+    {
+      entry = hsh_next (cov->ca, &iter);
+    }
+
   for (i = 0; i < design_matrix_get_n_cols (cov->cov); i++)
     {
       sum_i = get_sum (cov, i);
@@ -880,7 +883,6 @@ covariance_accumulator_to_matrix (struct covariance_matrix *cov)
        {
          sum_j = get_sum (cov, j);
          entry = hsh_first (cov->ca, &iter);
-         design_matrix_set_element (cov->sums, i, j, sum_i);     
          while (entry != NULL)
            {
              update_ssize (cov->ssize, i, j, entry);
@@ -889,14 +891,14 @@ covariance_accumulator_to_matrix (struct covariance_matrix *cov)
              */
              if (is_covariance_contributor (entry, cov->cov, i, j))
                {
-                 covariance_matrix_insert (cov->cov, entry->v1, entry->v2, entry->val1,
-                                           entry->val2, entry->dot_product);
+                 design_matrix_set_element (cov->cov, i, j, entry->dot_product);
                }
              entry = hsh_next (cov->ca, &iter);
            }
          tmp = design_matrix_get_element (cov->cov, i, j);
          tmp -= sum_i * sum_j / design_matrix_get_element (cov->ssize, i, j);
          design_matrix_set_element (cov->cov, i, j, tmp);
+         design_matrix_set_element (cov->cov, j, i, tmp);
        } 
     }
 }
@@ -923,6 +925,11 @@ covariance_to_design (const struct covariance_matrix *c)
     }
   return NULL;
 }
+size_t
+covariance_matrix_get_n_rows (const struct covariance_matrix *c)
+{
+  return design_matrix_get_n_rows (c->cov);
+}
 
 double 
 covariance_matrix_get_element (const struct covariance_matrix *c, size_t row, size_t col)