/* 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
struct covariance_matrix
{
struct design_matrix *cov;
+ struct design_matrix *ssize;
struct hsh_table *ca;
struct moments1 **m1;
struct moments **m;
const struct variable **v_variables;
+ const struct interaction_variable **interactions;
size_t n_variables;
+ size_t n_intr;
int n_pass;
int missing_handling;
enum mv_class missing_value;
result->ca = covariance_hsh_create (&result->n_variables);
result->m = NULL;
result->m1 = NULL;
+ result->n_intr = 0;
result->missing_handling = missing_handling;
result->missing_value = missing_value;
result->accumulate = (result->missing_handling == LISTWISE) ?
return result;
}
+void
+covariance_interaction_set (struct covariance_matrix *cov,
+ const struct interaction_variable **intr, size_t n_intr)
+{
+ cov->interactions = intr;
+ cov->n_intr = n_intr;
+}
+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.
*/
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 size_t
+get_n_rows_s (const struct variable *var)
+{
+ size_t result = 0;
+ if (var_is_numeric (var))
+ {
+ result++;
+ }
+ else
+ {
+ result += cat_get_n_categories (var) - 1;
+ }
+ return result;
+}
+static struct design_matrix *
+covariance_matrix_create_s (struct covariance_matrix *cov)
+{
+ struct variable **v_variables;
+ size_t n_variables;
+ size_t n_rows = 0;
+ size_t i;
+ size_t j;
+
+ n_variables = cov->n_variables + cov->n_intr;
+ v_variables = xnmalloc (n_variables, sizeof (*v_variables));
+ for (i = 0; i < cov->n_variables; i++)
+ {
+ v_variables[i] = cov->v_variables[i];
+ n_rows += get_n_rows_s (v_variables[i]);
+ }
+ for (j = 0; j < cov->n_intr; j++)
+ {
+ v_variables[i + j] = interaction_get_variable (cov->interactions[j]);
+ n_rows += get_n_rows_s (v_variables[i]);
+ }
+ return design_matrix_create (n_variables, v_variables, n_rows);
}
static void
assert (cov != NULL);
design_matrix_destroy (cov->cov);
+ design_matrix_destroy (cov->ssize);
hsh_destroy (cov->ca);
if (cov->n_pass == ONE_PASS)
{
col = design_matrix_var_to_column (cov, v2);
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);
+ tmp = design_matrix_get_element (cov, row, col);
+ design_matrix_set_element (cov, row, col, (val2->f - mean) * x + tmp);
+ design_matrix_set_element (cov, 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)
{
+ int width = var_get_width (v);
size_t col;
size_t i;
double y;
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;
}
- 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);
+ tmp = design_matrix_get_element (cov, row, col);
+ design_matrix_set_element (cov, row, col, x * y + tmp);
+ design_matrix_set_element (cov, col, row, x * y + tmp);
}
}
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;
}
row = design_matrix_var_to_column (cov, v1);
col = design_matrix_var_to_column (cov, v2);
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);
+ x += design_matrix_get_element (cov, col, row);
+ design_matrix_set_element (cov, row, col, x);
+ design_matrix_set_element (cov, col, row, x);
}
}
}
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;
}
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
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);
}
/*
if (a1 == NULL || a2 == NULL)
return 1;
-
+
return match_nodes (a1, a2->v1, a2->v2, a2->val1, a2->val2);
}
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))
{
}
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)
return cov->v_variables;
}
-
static void
-update_hash_entry (struct hsh_table *c,
- const struct variable *v1,
- const struct variable *v2,
- const union value *val1, const union value *val2,
- const struct interaction_value *i_val1,
- const struct interaction_value *i_val2)
+update_hash_entry_intr (struct hsh_table *c,
+ const struct variable *v1,
+ const struct variable *v2,
+ const union value *val1, const union value *val2,
+ const struct interaction_value *i_val1,
+ const struct interaction_value *i_val2)
{
struct covariance_accumulator *ca;
struct covariance_accumulator *new_entry;
new_entry->sum1 += ca->sum1;
new_entry->sum2 += ca->sum2;
/*
- If DOT_PRODUCT is null, CA was not already in the hash
- hable, so we don't free it because it was just inserted.
- If DOT_PRODUCT was not null, CA is already in the hash table.
- Unnecessary now, it must be freed here.
- */
+ If DOT_PRODUCT is null, CA was not already in the hash
+ hable, so we don't free it because it was just inserted.
+ If DOT_PRODUCT was not null, CA is already in the hash table.
+ Unnecessary now, it must be freed here.
+ */
+ free (ca);
+ }
+}
+
+static void
+update_hash_entry (struct hsh_table *c,
+ const struct variable *v1,
+ const struct variable *v2,
+ const union value *val1, const union value *val2)
+{
+ struct covariance_accumulator *ca;
+ struct covariance_accumulator *new_entry;
+
+ ca = get_new_covariance_accumulator (v1, v2, val1, val2);
+ ca->dot_product = update_product (ca->v1, ca->v2, ca->val1, ca->val2);
+ ca->sum1 = update_sum (ca->v1, ca->val1, 1.0);
+ ca->sum2 = update_sum (ca->v2, ca->val2, 1.0);
+ ca->ssize = 1.0;
+ new_entry = hsh_insert (c, ca);
+
+ if (new_entry != NULL)
+ {
+ new_entry->dot_product += ca->dot_product;
+ new_entry->ssize += 1.0;
+ new_entry->sum1 += ca->sum1;
+ new_entry->sum2 += ca->sum2;
+ /*
+ If DOT_PRODUCT is null, CA was not already in the hash
+ hable, so we don't free it because it was just inserted.
+ If DOT_PRODUCT was not null, CA is already in the hash table.
+ Unnecessary now, it must be freed here.
+ */
free (ca);
}
}
+static void
+inner_intr_loop (struct covariance_matrix *cov, const struct ccase *ccase, const struct variable *var1,
+ const union value *val1, const struct interaction_variable **i_var,
+ const struct interaction_value *i_val1, size_t j)
+{
+ struct variable *var2;
+ union value *val2;
+ struct interaction_value *i_val2;
+
+ var2 = interaction_get_variable (i_var[j]);
+ i_val2 = interaction_case_data (ccase, i_var[j]);
+ val2 = interaction_value_get (i_val2);
+
+ if (!var_is_value_missing (var2, val2, cov->missing_value))
+ {
+ update_hash_entry_intr (cov->ca, var1, var2, val1, val2, i_val1, i_val2);
+ }
+}
/*
Compute the covariance matrix in a single data-pass. Cases with
missing values are dropped pairwise, in other words, only if one of
const union value *val1;
const union value *val2;
const struct variable **v_variables;
+ const struct variable *var1;
+ const struct variable *var2;
struct interaction_value *i_val1 = NULL;
struct interaction_value *i_val2 = NULL;
for (i = 0; i < cov->n_variables; ++i)
{
- if (is_interaction (v_variables[i], i_var, n_intr))
- {
- i_val1 = interaction_case_data (ccase, v_variables[i], i_var, n_intr);
- val1 = interaction_value_get (i_val1);
- }
- else
+ var1 = v_variables[i];
+ val1 = case_data (ccase, var1);
+ if (!var_is_value_missing (var1, val1, cov->missing_value))
{
- val1 = case_data (ccase, v_variables[i]);
- }
- if (!var_is_value_missing (v_variables[i], val1, cov->missing_value))
- {
- cat_value_update (v_variables[i], val1);
- if (var_is_numeric (v_variables[i]))
+ cat_value_update (var1, val1);
+ if (var_is_numeric (var1))
cov->update_moments (cov, i, val1->f);
for (j = i; j < cov->n_variables; j++)
{
- if (is_interaction (v_variables[j], i_var, n_intr))
- {
- i_val2 = interaction_case_data (ccase, v_variables[j], i_var, n_intr);
- val2 = interaction_value_get (i_val2);
- }
- else
- {
- val2 = case_data (ccase, v_variables[j]);
- }
+ var2 = v_variables[j];
+ val2 = case_data (ccase, var2);
if (!var_is_value_missing
- (v_variables[j], val2, cov->missing_value))
+ (var2, val2, cov->missing_value))
{
- 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);
+ update_hash_entry (cov->ca, var1, var2, val1, val2);
}
}
+ for (j = 0; j < cov->n_intr; j++)
+ {
+ inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
+ }
+ }
+ }
+ for (i = 0; i < cov->n_intr; i++)
+ {
+ var1 = interaction_get_variable (i_var[i]);
+ i_val1 = interaction_case_data (ccase, i_var[i]);
+ val1 = interaction_value_get (i_val1);
+ cat_value_update (var1, val1);
+ if (!var_is_value_missing (var1, val1, cov->missing_value))
+ {
+ for (j = i; j < cov->n_intr; j++)
+ {
+ inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
+ }
}
}
}
for (i = 0; i < cov->n_variables; ++i)
{
- if (is_interaction (v_variables[i], i_var, n_intr))
- {
- i_val1 = interaction_case_data (ccase, v_variables[i], i_var, n_intr);
- val1 = interaction_value_get (i_val1);
- }
- else
- {
- val1 = case_data (ccase, v_variables[i]);
- }
+ val1 = case_data (ccase, v_variables[i]);
cat_value_update (v_variables[i], val1);
if (var_is_numeric (v_variables[i]))
cov->update_moments (cov, i, val1->f);
for (j = i; j < cov->n_variables; j++)
{
- if (is_interaction (v_variables[j], i_var, n_intr))
- {
- i_val2 = interaction_case_data (ccase, v_variables[j], i_var, n_intr);
- val2 = interaction_value_get (i_val2);
- }
- else
- {
- val2 = case_data (ccase, v_variables[j]);
- }
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);
+ val1, val2);
}
}
}
cov->accumulate (cov, ccase, (const struct interaction_variable **) aux, n_intr);
}
-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)
+/*
+ Return the value corresponding to subscript TARGET. If that value corresponds
+ to the origin, return NULL.
+ */
+static const union value *
+get_value_from_subscript (const struct design_matrix *dm, size_t target)
{
- size_t row;
- size_t col;
+ const union value *result = NULL;
+ const struct variable *var;
size_t i;
- const union value *tmp_val;
-
- assert (cov != NULL);
-
- row = design_matrix_var_to_column (cov, v1);
- if (var_is_alpha (v1))
+
+ var = design_matrix_col_to_var (dm, target);
+ if (var_is_numeric (var))
+ {
+ return NULL;
+ }
+ for (i = 0; i < cat_get_n_categories (var); i++)
{
- i = 0;
- tmp_val = cat_subscript_to_value (i, v1);
- while (compare_values_short (tmp_val, val1, v1))
+ result = cat_subscript_to_value (i, var);
+ if (dm_get_exact_subscript (dm, var, result) == target)
{
- i++;
- tmp_val = cat_subscript_to_value (i, v1);
+ return result;
}
- row += i;
- if (var_is_numeric (v2))
+ }
+ return NULL;
+}
+
+static bool
+is_covariance_contributor (const struct covariance_accumulator *ca, const struct design_matrix *dm,
+ size_t i, size_t j)
+{
+ 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))
+ {
+ if (var_get_dict_index (v2) == var_get_dict_index (ca->v2))
{
- col = design_matrix_var_to_column (cov, v2);
+ k = dm_get_exact_subscript (dm, v1, ca->val1);
+ if (k == i)
+ {
+ k = dm_get_exact_subscript (dm, v2, ca->val2);
+ if (k == j)
+ {
+ return true;
+ }
+ }
}
- else
+ }
+ 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))
{
- col = design_matrix_var_to_column (cov, v2);
- i = 0;
- tmp_val = cat_subscript_to_value (i, v1);
- while (compare_values_short (tmp_val, val1, v1))
+ k = dm_get_exact_subscript (dm, v1, ca->val2);
+ if (k == i)
{
- i++;
- tmp_val = cat_subscript_to_value (i, v1);
+ k = dm_get_exact_subscript (dm, v2, ca->val1);
+ if (k == j)
+ {
+ return true;
+ }
}
- col += i;
}
}
- else
+
+ 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;
+
+ assert ( cov != NULL);
+ var = design_matrix_col_to_var (cov->cov, i);
+ if (var != NULL)
{
- if (var_is_numeric (v2))
+ if (var_is_alpha (var))
{
- col = design_matrix_var_to_column (cov, v2);
+ val = get_value_from_subscript (cov->cov, i);
+ k = cat_value_find (var, val);
+ return cat_get_category_count (k, var);
}
else
{
- covariance_matrix_insert (cov, v2, v1, val2, val1, product);
+ k = 0;
+ while (cov->v_variables[k] != var && k < cov->n_variables)
+ {
+ k++;
+ }
+ if (k < cov->n_variables)
+ {
+ moments1_calculate (cov->m1[k], &n, &mean, NULL, NULL, NULL);
+ return mean * n;
+ }
}
}
-
- gsl_matrix_set (cov->m, row, col, product);
+
+ return 0.0;
}
-
-static struct design_matrix *
-covariance_accumulator_to_matrix (struct covariance_matrix *cov)
+static void
+update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_accumulator *ca)
{
+ 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))
+ {
+ var = design_matrix_col_to_var (dm, j);
+ if (var_get_dict_index (ca->v2) == var_get_dict_index (var))
+ {
+ tmp = design_matrix_get_element (dm, i, j);
+ tmp += ca->ssize;
+ design_matrix_set_element (dm, i, j, tmp);
+ }
+ }
+}
+static void
+covariance_accumulator_to_matrix (struct covariance_matrix *cov)
+{
+ size_t i;
+ size_t j;
+ double sum_i = 0.0;
+ double sum_j = 0.0;
+ double tmp = 0.0;
struct covariance_accumulator *entry;
- struct design_matrix *result = NULL;
struct hsh_iterator iter;
- result = covariance_matrix_create (cov->n_variables, cov->v_variables);
-
+ cov->cov = covariance_matrix_create_s (cov);
+ cov->ssize = covariance_matrix_create_s (cov);
entry = hsh_first (cov->ca, &iter);
-
while (entry != NULL)
{
- /*
- We compute the centered, un-normalized covariance matrix.
- */
- tmp = entry->dot_product - entry->sum1 * entry->sum2 / entry->ssize;
- covariance_matrix_insert (result, entry->v1, entry->v2, entry->val1,
- entry->val2, tmp);
entry = hsh_next (cov->ca, &iter);
}
- return result;
+
+ for (i = 0; i < design_matrix_get_n_cols (cov->cov); i++)
+ {
+ sum_i = get_sum (cov, i);
+ for (j = i; j < design_matrix_get_n_cols (cov->cov); j++)
+ {
+ sum_j = get_sum (cov, j);
+ entry = hsh_first (cov->ca, &iter);
+ while (entry != NULL)
+ {
+ update_ssize (cov->ssize, i, j, entry);
+ /*
+ We compute the centered, un-normalized covariance matrix.
+ */
+ if (is_covariance_contributor (entry, cov->cov, i, j))
+ {
+ 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);
+ }
+ }
}
{
if (cov->n_pass == ONE_PASS)
{
- cov->cov = covariance_accumulator_to_matrix (cov);
+ covariance_accumulator_to_matrix (cov);
}
}
}
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)
+{
+ return (design_matrix_get_element (c->cov, row, col));
+}
+