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;
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)
{
}
if (var_is_alpha (v_max) && var_is_alpha (v_min))
{
- unsigned int tmp;
- char *x =
- xnmalloc (1 + var_get_width (v_max) + var_get_width (v_min),
- sizeof (*x));
- strncpy (x, val_max->s, var_get_width (v_max));
- strncat (x, val_min->s, var_get_width (v_min));
- tmp = *n_vars * (*n_vars + 1 + idx_max) + idx_min + hsh_hash_string (x);
- free (x);
- return tmp;
+ 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;
}
return -1u;
}
return cov->v_variables;
}
-
static void
update_hash_entry (struct hsh_table *c,
const struct variable *v1,
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;
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);
}
}
{
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.
+ */
+static size_t
+get_exact_subscript (const struct design_matrix *dm, const struct variable *var,
+ const union value *val)
+{
+ size_t result;
+
+ if (is_origin (var, val))
+ {
+ return -1u;
+ }
+
+ result = design_matrix_var_to_column (dm, var);
+ if (var_is_alpha (var))
+ {
+ result += cat_value_find (var, val) - 1;
+ }
+ return result;
+}
static void
covariance_matrix_insert (struct design_matrix *cov,
{
size_t row;
size_t col;
- size_t i;
- const union value *tmp_val;
assert (cov != NULL);
- row = design_matrix_var_to_column (cov, v1);
- if (var_is_alpha (v1))
+ row = get_exact_subscript (cov, v1, val1);
+ col = get_exact_subscript (cov, v2, val2);
+ if (row != -1u && col != -1u)
{
- i = 0;
- tmp_val = cat_subscript_to_value (i, v1);
- while (compare_values_short (tmp_val, val1, v1))
- {
- i++;
- tmp_val = cat_subscript_to_value (i, v1);
- }
- row += i;
- if (var_is_numeric (v2))
- {
- col = design_matrix_var_to_column (cov, v2);
- }
- else
+ gsl_matrix_set (cov->m, row, col, product);
+ }
+}
+
+
+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);
+ 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))
{
- 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 = get_exact_subscript (dm, v1, ca->val1);
+ if (k == i)
{
- i++;
- tmp_val = cat_subscript_to_value (i, v1);
+ k = get_exact_subscript (dm, v2, ca->val2);
+ 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;
+ 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);
+ if (var != NULL)
{
- if (var_is_numeric (v2))
+ if (var_is_alpha (var))
{
- col = design_matrix_var_to_column (cov, v2);
+ k = design_matrix_var_to_column (cov->cov, var);
+ i -= k;
+ val = cat_subscript_to_value (i, var);
}
- else
+ ca.v1 = var;
+ ca.v2 = var;
+ ca.val1 = val;
+ ca.val2 = val;
+ c = (struct covariance_accumulator *) hsh_find (cov->ca, &ca);
+ if (c != NULL)
{
- covariance_matrix_insert (cov, v2, v1, val2, val1, product);
+ return c->sum1;
}
}
-
- 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)
{
+ 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 = gsl_matrix_get (dm->m, i, j);
+ tmp += ca->ssize;
+ gsl_matrix_set (dm->m, 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);
-
- entry = hsh_first (cov->ca, &iter);
-
- while (entry != NULL)
+ 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);
+ for (i = 0; i < design_matrix_get_n_cols (cov->cov); i++)
{
- /*
- 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);
+ 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);
+ gsl_matrix_set (cov->sums->m, i, j, sum_i);
+ 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))
+ {
+ covariance_matrix_insert (cov->cov, entry->v1, entry->v2, entry->val1,
+ entry->val2, entry->dot_product);
+ }
+ entry = hsh_next (cov->ca, &iter);
+ }
+ tmp = gsl_matrix_get (cov->cov->m, i, j);
+ tmp -= sum_i * sum_j / gsl_matrix_get (cov->ssize->m, i, j);
+ gsl_matrix_set (cov->cov->m, i, j, tmp);
+ }
}
- return result;
}
{
if (cov->n_pass == ONE_PASS)
{
- cov->cov = covariance_accumulator_to_matrix (cov);
+ covariance_accumulator_to_matrix (cov);
}
}
projects / pspp-builds.git / blobdiff