X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance.c;h=3ec90efbba05188c68bddd54ce53bc2a8e7a5d0f;hb=ca4012bcf0f8790ceb8539b55bbc296d0802d5d7;hp=446ee6621a1e11881c37e5dba04d8d645c1610a8;hpb=81579d9e9f994fb2908f50af41c3eb033d216e58;p=pspp diff --git a/src/math/covariance.c b/src/math/covariance.c index 446ee6621a..3ec90efbba 100644 --- a/src/math/covariance.c +++ b/src/math/covariance.c @@ -25,6 +25,7 @@ #include "libpspp/assertion.h" #include "libpspp/misc.h" #include "math/categoricals.h" +#include "math/interaction.h" #include "math/moments.h" #include "gl/xalloc.h" @@ -112,6 +113,8 @@ struct covariance /* Flags indicating that the first case has been seen */ bool pass_one_first_case_seen; bool pass_two_first_case_seen; + + gsl_matrix *unnormalised; }; @@ -159,8 +162,8 @@ covariance_1pass_create (size_t n_vars, const struct variable *const *vars, cov->n_cm = (n_vars * (n_vars - 1) ) / 2; - if (cov->n_cm > 0) - cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm); + + cov->cm = xcalloc (cov->n_cm, sizeof *cov->cm); cov->categoricals = NULL; return cov; @@ -201,6 +204,7 @@ covariance_2pass_create (size_t n_vars, const struct variable *const *vars, cov->cm = NULL; cov->categoricals = cats; + cov->unnormalised = NULL; return cov; } @@ -247,7 +251,7 @@ is_missing (const struct covariance *cov, int i, const struct ccase *c) { const struct variable *var = i < cov->n_vars ? cov->vars[i] : - categoricals_get_variable_by_subscript (cov->categoricals, i - cov->n_vars); + categoricals_get_interaction_by_subscript (cov->categoricals, i - cov->n_vars)->vars[0]; const union value *val = case_data (c, var); @@ -267,7 +271,7 @@ get_val (const struct covariance *cov, int i, const struct ccase *c) return val->f; } - return categoricals_get_binary_by_subscript (cov->categoricals, i - cov->n_vars, c); + return categoricals_get_effects_code_for_case (cov->categoricals, i - cov->n_vars, c); } #if 0 @@ -346,14 +350,16 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c) assert (cov->state == 1); cov->state = 2; + if (cov->categoricals) + categoricals_done (cov->categoricals); + cov->dim = cov->n_vars; if (cov->categoricals) - cov->dim += categoricals_total (cov->categoricals) - - categoricals_get_n_variables (cov->categoricals); + cov->dim += categoricals_df_total (cov->categoricals); cov->n_cm = (cov->dim * (cov->dim - 1) ) / 2; - cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm); + cov->cm = xcalloc (cov->n_cm, sizeof *cov->cm); /* Grow the moment matrices so that they're large enough to accommodate the categorical elements */ @@ -362,9 +368,6 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c) cov->moments[i] = resize_matrix (cov->moments[i], cov->dim); } - if (cov->categoricals) - categoricals_done (cov->categoricals); - /* Populate the moments matrices with the categorical value elements */ for (i = cov->n_vars; i < cov->dim; ++i) { @@ -631,22 +634,6 @@ covariance_calculate (struct covariance *cov) static gsl_matrix * covariance_calculate_double_pass_unnormalized (struct covariance *cov) { - size_t i, j; - for (i = 0 ; i < cov->dim; ++i) - { - for (j = 0 ; j < cov->dim; ++j) - { - int idx; - double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j); - - idx = cm_idx (cov, i, j); - if ( idx >= 0) - { - x = &cov->cm [idx]; - } - } - } - return cm_to_gsl (cov); } @@ -683,27 +670,31 @@ covariance_calculate_single_pass_unnormalized (struct covariance *cov) /* Return a pointer to gsl_matrix containing the pairwise covariances. The - caller owns the returned matrix and must free it when it is no longer - needed. + returned matrix is owned by the structure, and must not be freed. Call this function only after all data have been accumulated. */ -gsl_matrix * +const gsl_matrix * covariance_calculate_unnormalized (struct covariance *cov) { if ( cov->state <= 0 ) return NULL; + if (cov->unnormalised != NULL) + return cov->unnormalised; + switch (cov->passes) { case 1: - return covariance_calculate_single_pass_unnormalized (cov); + cov->unnormalised = covariance_calculate_single_pass_unnormalized (cov); break; case 2: - return covariance_calculate_double_pass_unnormalized (cov); + cov->unnormalised = covariance_calculate_double_pass_unnormalized (cov); break; default: NOT_REACHED (); } + + return cov->unnormalised; } /* Function to access the categoricals used by COV @@ -727,7 +718,101 @@ covariance_destroy (struct covariance *cov) for (i = 0; i < n_MOMENTS; ++i) gsl_matrix_free (cov->moments[i]); + gsl_matrix_free (cov->unnormalised); free (cov->moments); free (cov->cm); free (cov); } + +size_t +covariance_dim (const struct covariance * cov) +{ + return (cov->dim); +} + + + +/* + Routines to assist debugging. + The following are not thoroughly tested and in certain respects + unreliable. They should only be + used for aids to development. Not as user accessible code. +*/ + +#include "libpspp/str.h" +#include "output/tab.h" +#include "data/format.h" + + +/* Create a table which can be populated with the encodings for + the covariance matrix COV */ +struct tab_table * +covariance_dump_enc_header (const struct covariance *cov, int length) +{ + struct tab_table *t = tab_create (cov->dim, length); + int n; + int i; + + tab_title (t, "Covariance Encoding"); + + tab_box (t, + TAL_2, TAL_2, 0, 0, + 0, 0, tab_nc (t) - 1, tab_nr (t) - 1); + + tab_hline (t, TAL_2, 0, tab_nc (t) - 1, 1); + + + for (i = 0 ; i < cov->n_vars; ++i) + { + tab_text (t, i, 0, TAT_TITLE, var_get_name (cov->vars[i])); + tab_vline (t, TAL_1, i + 1, 0, tab_nr (t) - 1); + } + + n = 0; + while (i < cov->dim) + { + struct string str; + int idx = i - cov->n_vars; + const struct interaction *iact = + categoricals_get_interaction_by_subscript (cov->categoricals, idx); + int df; + + ds_init_empty (&str); + interaction_to_string (iact, &str); + + df = categoricals_df (cov->categoricals, n); + + tab_joint_text (t, + i, 0, + i + df - 1, 0, + TAT_TITLE, ds_cstr (&str)); + + if (i + df < tab_nr (t) - 1) + tab_vline (t, TAL_1, i + df, 0, tab_nr (t) - 1); + + i += df; + n++; + ds_destroy (&str); + } + + return t; +} + + +/* + Append table T, which should have been returned by covariance_dump_enc_header + with an entry corresponding to case C for the covariance matrix COV + */ +void +covariance_dump_enc (const struct covariance *cov, const struct ccase *c, + struct tab_table *t) +{ + static int row = 0; + int i; + ++row; + for (i = 0 ; i < cov->dim; ++i) + { + double v = get_val (cov, i, c); + tab_double (t, i, row, 0, v, i < cov->n_vars ? NULL : &F_8_0, RC_OTHER); + } +}