X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance-matrix.c;h=2e83749c6ed27155235dac0f4bb37e9d98e212ea;hb=5d31af536467b5774fb04ebc998ef459e893f24c;hp=a33c33b243e46c07c8997a8825f27b00cfbc458b;hpb=d3181522c4d8e725d417cbae655e46f660fe168e;p=pspp-builds.git diff --git a/src/math/covariance-matrix.c b/src/math/covariance-matrix.c index a33c33b2..2e83749c 100644 --- a/src/math/covariance-matrix.c +++ b/src/math/covariance-matrix.c @@ -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,7 @@ struct covariance_matrix { struct design_matrix *cov; struct design_matrix *ssize; - struct design_matrix *means; + struct design_matrix *sums; struct hsh_table *ca; struct moments1 **m1; struct moments **m; @@ -148,7 +148,25 @@ covariance_matrix_init (size_t n_variables, return result; } - +static size_t +get_n_rows (size_t n_variables, size_t *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 +174,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 +200,7 @@ 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); + design_matrix_destroy (cov->sums); hsh_destroy (cov->ca); if (cov->n_pass == ONE_PASS) { @@ -219,9 +237,9 @@ covariance_update_categorical_numeric (struct design_matrix *cov, double mean, 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, @@ -243,9 +261,9 @@ column_iterate (struct design_matrix *cov, const struct variable *v, { 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); } } @@ -306,9 +324,9 @@ covariance_pass_two (struct design_matrix *cov, double mean1, double mean2, 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); } } @@ -354,10 +372,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 = hash_bytes (val_max, var_get_width (v_max), 0); + hash = hash_bytes (val_min, var_get_width (v_min), hash); + return hash_int (*n_vars * (*n_vars + 1 + idx_max) + idx_min, hash); } return -1u; } @@ -460,7 +477,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) + hash_string (val->s, 0); } else if (var_is_alpha (v1) && var_is_numeric (v2)) { @@ -778,7 +795,7 @@ covariance_matrix_insert (struct design_matrix *cov, col = get_exact_subscript (cov, v2, val2); if (row != -1u && col != -1u) { - gsl_matrix_set (cov->m, row, col, product); + design_matrix_set_element (cov, row, col, product); } } @@ -811,6 +828,37 @@ is_covariance_contributor (const struct covariance_accumulator *ca, const struct } 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_alpha (var)) + { + k = design_matrix_var_to_column (cov->cov, var); + i -= k; + val = cat_subscript_to_value (i, 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) + { + return c->sum1; + } + } + return 0.0; +} static void update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_accumulator *ca) { @@ -822,9 +870,9 @@ update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_ac 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 = design_matrix_get_element (dm, i, j); tmp += ca->ssize; - gsl_matrix_set (dm->m, i, j, tmp); + design_matrix_set_element (dm, i, j, tmp); } } } @@ -833,23 +881,26 @@ covariance_accumulator_to_matrix (struct covariance_matrix *cov) { size_t i; size_t j; - double tmp; + double sum_i = 0.0; + double sum_j = 0.0; + double tmp = 0.0; struct covariance_accumulator *entry; struct hsh_iterator iter; 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); + cov->sums = 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); - + design_matrix_set_element (cov->sums, i, j, sum_i); while (entry != NULL) { update_ssize (cov->ssize, i, j, entry); - /* We compute the centered, un-normalized covariance matrix. */ @@ -857,14 +908,12 @@ covariance_accumulator_to_matrix (struct covariance_matrix *cov) { covariance_matrix_insert (cov->cov, entry->v1, entry->v2, entry->val1, entry->val2, entry->dot_product); - covariance_matrix_insert (cov->cov, entry->v1, entry->v2, entry->val1, - entry->val2, entry->sum1 * entry->sum2); } 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); + 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); } } } @@ -891,3 +940,10 @@ covariance_to_design (const struct covariance_matrix *c) } return NULL; } + +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)); +} +