X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance-matrix.c;h=2ab5f406fade36931a94412bbf42196541b24874;hb=5282b76bcbd0501c1227b36e4d451f340388eca5;hp=c419ce02e8b2357a868526aff9147f172b1125d7;hpb=f05f98b6b0b7e81f77657a4313fccc259eb5c506;p=pspp-builds.git diff --git a/src/math/covariance-matrix.c b/src/math/covariance-matrix.c index c419ce02..2ab5f406 100644 --- a/src/math/covariance-matrix.c +++ b/src/math/covariance-matrix.c @@ -56,6 +56,8 @@ struct covariance_accumulator 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; @@ -179,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->sums); hsh_destroy (cov->ca); if (cov->n_pass == ONE_PASS) { @@ -528,7 +532,6 @@ get_covariance_variables (const struct covariance_matrix *cov) return cov->v_variables; } - static void update_hash_entry (struct hsh_table *c, const struct variable *v1, @@ -551,7 +554,7 @@ update_hash_entry (struct hsh_table *c, 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; @@ -559,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); } } @@ -721,6 +724,44 @@ covariance_matrix_accumulate (struct covariance_matrix *cov, { 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, @@ -730,77 +771,134 @@ 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; } @@ -812,7 +910,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); } }