X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance-matrix.c;h=972a3ff00726a12200bbdab00086180b0ec100df;hb=428b694be3fe5ab2b4ba430f3a1d12d60a8e510c;hp=4f14a7c60517f7b97379d0e963a8eff0c74fce1d;hpb=1870f69a33f783ca49bf1827ba13a0463dbc189a;p=pspp-builds.git diff --git a/src/math/covariance-matrix.c b/src/math/covariance-matrix.c index 4f14a7c6..972a3ff0 100644 --- a/src/math/covariance-matrix.c +++ b/src/math/covariance-matrix.c @@ -26,7 +26,6 @@ #include #include #include -#include #include #include #include @@ -65,11 +64,14 @@ struct covariance_matrix int n_pass; int missing_handling; enum mv_class missing_value; - void (*accumulate) (struct covariance_matrix *, const struct ccase *); + void (*accumulate) (struct covariance_matrix *, const struct ccase *, + const struct interaction_variable **, size_t); void (*update_moments) (struct covariance_matrix *, size_t, double); }; -static struct hsh_table *covariance_hsh_create (size_t); + + +static struct hsh_table *covariance_hsh_create (size_t *); static hsh_hash_func covariance_accumulator_hash; static unsigned int hash_numeric_alpha (const struct variable *, const struct variable *, @@ -94,9 +96,13 @@ static struct covariance_accumulator *get_new_covariance_accumulator (const value *); static void covariance_accumulate_listwise (struct covariance_matrix *, - const struct ccase *); + const struct ccase *, + const struct interaction_variable **, + size_t); static void covariance_accumulate_pairwise (struct covariance_matrix *, - const struct ccase *); + const struct ccase *, + const struct interaction_variable **, + size_t); struct covariance_matrix * covariance_matrix_init (size_t n_variables, @@ -108,7 +114,8 @@ covariance_matrix_init (size_t n_variables, result = xmalloc (sizeof (*result)); result->cov = NULL; - result->ca = covariance_hsh_create (n_variables); + result->n_variables = n_variables; + result->ca = covariance_hsh_create (&result->n_variables); result->m = NULL; result->m1 = NULL; result->missing_handling = missing_handling; @@ -134,7 +141,7 @@ covariance_matrix_init (size_t n_variables, } } result->v_variables = v_variables; - result->n_variables = n_variables; + result->n_pass = n_pass; return result; @@ -228,7 +235,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 (!compare_values_short (tmp_val, val1, v)) { y += -1.0; } @@ -263,7 +270,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 (!compare_values_short (tmp_val, val1, v1)) { x += 1.0; } @@ -343,15 +350,10 @@ covariance_accumulator_hash (const void *h, const void *aux) } 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; } @@ -363,11 +365,11 @@ covariance_accumulator_hash (const void *h, const void *aux) in the data. */ static struct hsh_table * -covariance_hsh_create (size_t n_vars) +covariance_hsh_create (size_t *n_vars) { - return hsh_create (n_vars * n_vars, covariance_accumulator_compare, + return hsh_create (*n_vars * *n_vars, covariance_accumulator_compare, covariance_accumulator_hash, covariance_accumulator_free, - &n_vars); + n_vars); } static void @@ -414,7 +416,7 @@ match_nodes (const struct covariance_accumulator *c, } if (var_is_alpha (v1) && var_is_alpha (v2)) { - if (compare_values_short (val1, c->val1, v1)) + if (!compare_values_short (val1, c->val1, v1)) { if (!compare_values_short (val2, c->val2, v2)) { @@ -491,13 +493,13 @@ update_product (const struct variable *v1, const struct variable *v2, return 0.0; } static double -update_sum (const struct variable *var, const union value *val) +update_sum (const struct variable *var, const union value *val, double weight) { assert (var != NULL); assert (val != NULL); if (var_is_alpha (var)) { - return 1.0; + return weight; } return val->f; } @@ -530,18 +532,25 @@ 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 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; + double iv_f1; + double iv_f2; - + iv_f1 = interaction_value_get_nonzero_entry (i_val1); + iv_f2 = interaction_value_get_nonzero_entry (i_val2); 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); - ca->sum2 = update_sum (ca->v2, ca->val2); + ca->dot_product *= iv_f1 * iv_f2; + ca->sum1 = update_sum (ca->v1, ca->val1, iv_f1); + 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; @@ -549,11 +558,11 @@ update_hash_entry (struct hsh_table *c, 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); } } @@ -569,13 +578,17 @@ update_hash_entry (struct hsh_table *c, */ static void covariance_accumulate_pairwise (struct covariance_matrix *cov, - const struct ccase *ccase) + const struct ccase *ccase, + const struct interaction_variable **i_var, + size_t n_intr) { size_t i; size_t j; const union value *val1; const union value *val2; const struct variable **v_variables; + struct interaction_value *i_val1 = NULL; + struct interaction_value *i_val2 = NULL; assert (cov != NULL); assert (ccase != NULL); @@ -585,7 +598,15 @@ covariance_accumulate_pairwise (struct covariance_matrix *cov, for (i = 0; i < cov->n_variables; ++i) { - val1 = case_data (ccase, v_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]); + } if (!var_is_value_missing (v_variables[i], val1, cov->missing_value)) { cat_value_update (v_variables[i], val1); @@ -594,15 +615,23 @@ covariance_accumulate_pairwise (struct covariance_matrix *cov, for (j = i; j < cov->n_variables; j++) { - val2 = case_data (ccase, v_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]); + } if (!var_is_value_missing (v_variables[j], val2, cov->missing_value)) { update_hash_entry (cov->ca, v_variables[i], v_variables[j], - val1, val2); + val1, val2, i_val1, i_val2); if (j != i) update_hash_entry (cov->ca, v_variables[j], - v_variables[i], val2, val1); + v_variables[i], val2, val1, i_val2, i_val1); } } } @@ -626,13 +655,17 @@ covariance_accumulate_pairwise (struct covariance_matrix *cov, */ static void covariance_accumulate_listwise (struct covariance_matrix *cov, - const struct ccase *ccase) + const struct ccase *ccase, + const struct interaction_variable **i_var, + size_t n_intr) { size_t i; size_t j; const union value *val1; const union value *val2; const struct variable **v_variables; + struct interaction_value *i_val1 = NULL; + struct interaction_value *i_val2 = NULL; assert (cov != NULL); assert (ccase != NULL); @@ -642,19 +675,35 @@ covariance_accumulate_listwise (struct covariance_matrix *cov, for (i = 0; i < cov->n_variables; ++i) { - val1 = case_data (ccase, v_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]); + } 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++) { - val2 = case_data (ccase, v_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); + val1, val2, i_val1, i_val2); if (j != i) update_hash_entry (cov->ca, v_variables[j], v_variables[i], - val2, val1); + val2, val1, i_val2, i_val1); } } } @@ -663,13 +712,51 @@ covariance_accumulate_listwise (struct covariance_matrix *cov, Call this function during the data pass. Each case will be added to a hash containing all values of the covariance matrix. After the data have been passed, call covariance_matrix_compute to put the - values in the struct covariance_matrix. + values in the struct covariance_matrix. */ void covariance_matrix_accumulate (struct covariance_matrix *cov, - const struct ccase *ccase) + const struct ccase *ccase, void **aux, size_t n_intr) +{ + 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) { - cov->accumulate (cov, ccase); + 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 @@ -685,51 +772,49 @@ covariance_matrix_insert (struct design_matrix *cov, 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 - { - 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)) - { - i++; - tmp_val = cat_subscript_to_value (i, v1); - } - col += i; - } + gsl_matrix_set (cov->m, row, col, product); } - else +} + + +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)) { - if (var_is_numeric (v2)) - { - col = design_matrix_var_to_column (cov, v2); - } - else + v2 = design_matrix_col_to_var (dm, j); + if (var_get_dict_index (v2) == var_get_dict_index (ca->v2)) { - covariance_matrix_insert (cov, v2, v1, val2, val1, product); + k = get_exact_subscript (dm, v1, ca->val1); + if (k == i) + { + k = get_exact_subscript (dm, v2, ca->val2); + if (k == j) + { + return true; + } + } } } - gsl_matrix_set (cov->m, row, col, product); + return false; } - + static struct design_matrix * covariance_accumulator_to_matrix (struct covariance_matrix *cov) { + size_t i; + size_t j; double tmp; struct covariance_accumulator *entry; struct design_matrix *result = NULL; @@ -737,17 +822,26 @@ covariance_accumulator_to_matrix (struct covariance_matrix *cov) result = covariance_matrix_create (cov->n_variables, cov->v_variables); - entry = hsh_first (cov->ca, &iter); - - while (entry != NULL) + for (i = 0; i < design_matrix_get_n_cols (result); 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); + for (j = i; j < design_matrix_get_n_cols (result); j++) + { + entry = hsh_first (cov->ca, &iter); + + while (entry != NULL) + { + /* + We compute the centered, un-normalized covariance matrix. + */ + if (is_covariance_contributor (entry, result, i, j)) + { + 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; }