+ return NULL;
+}
+
+static const struct variable **
+get_covariance_variables (const struct covariance_matrix *cov)
+{
+ return cov->v_variables;
+}
+
+static void
+update_hash_entry_intr (struct hsh_table *c,
+ const struct variable *v1,
+ const struct variable *v2,
+ 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->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;
+ new_entry->ssize += 1.0;
+ 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.
+ */
+ free (ca);
+ }
+}
+
+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)
+{
+ struct covariance_accumulator *ca;
+ struct covariance_accumulator *new_entry;
+
+ 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, 1.0);
+ ca->sum2 = update_sum (ca->v2, ca->val2, 1.0);
+ ca->ssize = 1.0;
+ new_entry = hsh_insert (c, ca);
+
+ if (new_entry != NULL)
+ {
+ new_entry->dot_product += ca->dot_product;
+ new_entry->ssize += 1.0;
+ 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.
+ */
+ free (ca);
+ }
+}
+
+static void
+inner_intr_loop (struct covariance_matrix *cov, const struct ccase *ccase, const struct variable *var1,
+ const union value *val1, const struct interaction_variable **i_var,
+ const struct interaction_value *i_val1, size_t j)
+{
+ struct variable *var2;
+ union value *val2;
+ struct interaction_value *i_val2;
+
+ var2 = interaction_get_variable (i_var[j]);
+ i_val2 = interaction_case_data (ccase, i_var[j]);
+ val2 = interaction_value_get (i_val2);
+
+ if (!var_is_value_missing (var2, val2, cov->missing_value))
+ {
+ update_hash_entry_intr (cov->ca, var1, var2, val1, val2, i_val1, i_val2);
+ }
+}
+/*
+ Compute the covariance matrix in a single data-pass. Cases with
+ missing values are dropped pairwise, in other words, only if one of
+ the two values necessary to accumulate the inner product is missing.
+
+ Do not call this function directly. Call it through the struct
+ covariance_matrix ACCUMULATE member function, for example,
+ cov->accumulate (cov, ccase).
+ */
+static void
+covariance_accumulate_pairwise (struct covariance_matrix *cov,
+ 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;
+ const struct variable *var1;
+ const struct variable *var2;
+ struct interaction_value *i_val1 = NULL;
+ struct interaction_value *i_val2 = NULL;
+
+ assert (cov != NULL);
+ assert (ccase != NULL);
+
+ v_variables = get_covariance_variables (cov);
+ assert (v_variables != NULL);
+
+ for (i = 0; i < cov->n_variables; ++i)
+ {
+ var1 = v_variables[i];
+ val1 = case_data (ccase, var1);
+ if (!var_is_value_missing (var1, val1, cov->missing_value))
+ {
+ cat_value_update (var1, val1);
+ if (var_is_numeric (var1))
+ cov->update_moments (cov, i, val1->f);
+
+ for (j = i; j < cov->n_variables; j++)
+ {
+ var2 = v_variables[j];
+ val2 = case_data (ccase, var2);
+ if (!var_is_value_missing
+ (var2, val2, cov->missing_value))
+ {
+ update_hash_entry (cov->ca, var1, var2, val1, val2);
+ }
+ }
+ for (j = 0; j < cov->n_intr; j++)
+ {
+ inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
+ }
+ }
+ }
+ for (i = 0; i < cov->n_intr; i++)
+ {
+ var1 = interaction_get_variable (i_var[i]);
+ i_val1 = interaction_case_data (ccase, i_var[i]);
+ val1 = interaction_value_get (i_val1);
+ cat_value_update (var1, val1);
+ if (!var_is_value_missing (var1, val1, cov->missing_value))
+ {
+ for (j = i; j < cov->n_intr; j++)
+ {
+ inner_intr_loop (cov, ccase, var1, val1, i_var, i_val1, j);
+ }
+ }
+ }
+}
+
+/*
+ Compute the covariance matrix in a single data-pass. Cases with
+ missing values are dropped listwise. In other words, if one of the
+ values for any variable in a case is missing, the entire case is
+ skipped.
+
+ The caller must use a casefilter to remove the cases with missing
+ values before calling covariance_accumulate_listwise. This function
+ assumes that CCASE has already passed through this filter, and
+ contains no missing values.
+
+ Do not call this function directly. Call it through the struct
+ covariance_matrix ACCUMULATE member function, for example,
+ cov->accumulate (cov, ccase).
+ */
+static void
+covariance_accumulate_listwise (struct covariance_matrix *cov,
+ 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);
+
+ v_variables = get_covariance_variables (cov);
+ assert (v_variables != NULL);
+
+ for (i = 0; i < cov->n_variables; ++i)
+ {
+ 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++)
+ {
+ update_hash_entry (cov->ca, v_variables[i], v_variables[j],
+ val1, val2);
+ }
+ }
+}
+
+/*
+ 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.
+ */
+void
+covariance_matrix_accumulate (struct covariance_matrix *cov,
+ const struct ccase *ccase, void **aux, size_t n_intr)
+{
+ cov->accumulate (cov, ccase, (const struct interaction_variable **) aux, n_intr);
+}
+
+/*
+ Return the value corresponding to subscript TARGET. If that value corresponds
+ to the origin, return NULL.
+ */
+static const union value *
+get_value_from_subscript (const struct design_matrix *dm, size_t target)
+{
+ const union value *result = NULL;
+ const struct variable *var;
+ size_t i;
+
+ var = design_matrix_col_to_var (dm, target);
+ if (var_is_numeric (var))
+ {
+ return NULL;
+ }
+ for (i = 0; i < cat_get_n_categories (var); i++)
+ {
+ result = cat_subscript_to_value (i, var);
+ if (dm_get_exact_subscript (dm, var, result) == target)
+ {
+ return result;
+ }
+ }
+ return NULL;
+}
+
+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);
+ v2 = design_matrix_col_to_var (dm, j);
+ if (var_get_dict_index (v1) == var_get_dict_index(ca->v1))
+ {
+ if (var_get_dict_index (v2) == var_get_dict_index (ca->v2))
+ {
+ k = dm_get_exact_subscript (dm, v1, ca->val1);
+ if (k == i)
+ {
+ k = dm_get_exact_subscript (dm, v2, ca->val2);
+ if (k == j)
+ {
+ return true;
+ }
+ }
+ }
+ }
+ else if (var_get_dict_index (v1) == var_get_dict_index (ca->v2))
+ {
+ if (var_get_dict_index (v2) == var_get_dict_index (ca->v1))
+ {
+ k = dm_get_exact_subscript (dm, v1, ca->val2);
+ if (k == i)
+ {
+ k = dm_get_exact_subscript (dm, v2, ca->val1);
+ if (k == j)
+ {
+ return true;
+ }
+ }
+ }
+ }
+
+ return false;
+}
+static double
+get_sum (const struct covariance_matrix *cov, size_t i)
+{
+ size_t k;
+ double mean;
+ double n;
+ const struct variable *var;
+ const union value *val = NULL;
+
+ assert ( cov != NULL);
+ var = design_matrix_col_to_var (cov->cov, i);
+ if (var != NULL)
+ {
+ if (var_is_alpha (var))
+ {
+ val = get_value_from_subscript (cov->cov, i);
+ k = cat_value_find (var, val);
+ return cat_get_category_count (k, var);
+ }
+ else
+ {
+ k = 0;
+ while (cov->v_variables[k] != var && k < cov->n_variables)
+ {
+ k++;
+ }
+ if (k < cov->n_variables)
+ {
+ moments1_calculate (cov->m1[k], &n, &mean, NULL, NULL, NULL);
+ return mean * n;
+ }
+ }
+ }
+
+ return 0.0;
+}
+static void
+update_ssize (struct design_matrix *dm, size_t i, size_t j, struct covariance_accumulator *ca)
+{
+ const 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 = design_matrix_get_element (dm, i, j);
+ tmp += ca->ssize;
+ design_matrix_set_element (dm, 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 hsh_iterator iter;
+
+ cov->cov = covariance_matrix_create_s (cov);
+ cov->ssize = covariance_matrix_create_s (cov);
+ entry = hsh_first (cov->ca, &iter);
+ while (entry != NULL)
+ {
+ entry = hsh_next (cov->ca, &iter);
+ }
+
+ 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))
+ {
+ design_matrix_set_element (cov->cov, i, j, entry->dot_product);
+ }
+ entry = hsh_next (cov->ca, &iter);
+ }
+ 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);
+ design_matrix_set_element (cov->cov, j, i, tmp);
+ }
+ }
+}
+
+
+/*
+ Call this function after passing the data.
+ */
+void
+covariance_matrix_compute (struct covariance_matrix *cov)
+{
+ if (cov->n_pass == ONE_PASS)
+ {
+ covariance_accumulator_to_matrix (cov);
+ }
+}
+
+struct design_matrix *
+covariance_to_design (const struct covariance_matrix *c)
+{
+ if (c != NULL)
+ {
+ return c->cov;
+ }
+ return NULL;
+}
+size_t
+covariance_matrix_get_n_rows (const struct covariance_matrix *c)
+{
+ return design_matrix_get_n_rows (c->cov);
+}
+
+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));