X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance-matrix.c;h=82f4b4479ad18b2903339ded476c340075217c2a;hb=1288f8123ea04686150b79b9e6a5d5d78c420353;hp=b6411fb3d78262e2450f3cbbdd1bccd1844f9184;hpb=a1feacc810212d66d78dc95269396e36b01a7072;p=pspp-builds.git diff --git a/src/math/covariance-matrix.c b/src/math/covariance-matrix.c index b6411fb3..82f4b447 100644 --- a/src/math/covariance-matrix.c +++ b/src/math/covariance-matrix.c @@ -26,7 +26,6 @@ #include #include #include -#include #include #include #include @@ -57,6 +56,8 @@ struct covariance_accumulator struct covariance_matrix { struct design_matrix *cov; + struct design_matrix *ssize; + struct design_matrix *means; struct hsh_table *ca; struct moments1 **m1; struct moments **m; @@ -65,11 +66,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 +98,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 +116,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 +143,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; @@ -172,6 +181,8 @@ covariance_matrix_destroy (struct covariance_matrix *cov) assert (cov != NULL); design_matrix_destroy (cov->cov); + design_matrix_destroy (cov->ssize); + design_matrix_destroy (cov->means); hsh_destroy (cov->ca); if (cov->n_pass == ONE_PASS) { @@ -228,7 +239,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 +274,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 +354,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 +369,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 @@ -400,23 +406,23 @@ match_nodes (const struct covariance_accumulator *c, } if (var_is_numeric (v1) && var_is_alpha (v2)) { - if (compare_values_short (val2, c->val2, 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)) + 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 (val1, c->val1, v1)) { - if (compare_values_short (val2, c->val2, v2)) + if (!compare_values_short (val2, c->val2, v2)) { return 0; } @@ -491,13 +497,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 +536,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 +562,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 +582,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 +602,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 +619,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 +659,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 +679,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 +716,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) +{ + 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) { - cov->accumulate (cov, ccase); + 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 @@ -680,76 +771,136 @@ 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); + cov->cov = covariance_matrix_create (cov->n_variables, cov->v_variables); + cov->ssize = covariance_matrix_create (cov->n_variables, cov->v_variables); + cov->means = covariance_matrix_create (cov->n_variables, cov->v_variables); + 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)) + { - entry = hsh_first (cov->ca, &iter); + 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 -= gsl_matrix_get (cov->means->m, i, j) / gsl_matrix_get (cov->ssize->m, i, j); + gsl_matrix_set (cov->cov->m, i, j, tmp); - 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; } @@ -761,7 +912,7 @@ covariance_matrix_compute (struct covariance_matrix *cov) { if (cov->n_pass == ONE_PASS) { - cov->cov = covariance_accumulator_to_matrix (cov); + covariance_accumulator_to_matrix (cov); } }