+ x = (val1->f - mean1) * (val2->f - mean2);
+ x += design_matrix_get_element (cov, col, row);
+ design_matrix_set_element (cov, row, col, x);
+ design_matrix_set_element (cov, col, row, x);
+ }
+}
+
+static unsigned int
+covariance_accumulator_hash (const void *h, const void *aux)
+{
+ struct covariance_accumulator *ca = (struct covariance_accumulator *) h;
+ size_t *n_vars = (size_t *) aux;
+ size_t idx_max;
+ size_t idx_min;
+ const struct variable *v_min;
+ const struct variable *v_max;
+ const union value *val_min;
+ const union value *val_max;
+
+ /*
+ Order everything by the variables' indices. This ensures we get the
+ same key regardless of the order in which the variables are stored
+ and passed around.
+ */
+ v_min =
+ (var_get_dict_index (ca->v1) <
+ var_get_dict_index (ca->v2)) ? ca->v1 : ca->v2;
+ v_max = (ca->v1 == v_min) ? ca->v2 : ca->v1;
+
+ val_min = (v_min == ca->v1) ? ca->val1 : ca->val2;
+ val_max = (ca->val1 == val_min) ? ca->val2 : ca->val1;
+
+ idx_min = var_get_dict_index (v_min);
+ idx_max = var_get_dict_index (v_max);
+
+ if (var_is_numeric (v_max) && var_is_numeric (v_min))
+ {
+ return (*n_vars * idx_max + idx_min);
+ }
+ if (var_is_numeric (v_max) && var_is_alpha (v_min))
+ {
+ return hash_numeric_alpha (v_max, v_min, val_min, *n_vars);
+ }
+ if (var_is_alpha (v_max) && var_is_numeric (v_min))
+ {
+ return (hash_numeric_alpha (v_min, v_max, val_max, *n_vars));
+ }
+ if (var_is_alpha (v_max) && var_is_alpha (v_min))
+ {
+ unsigned hash = value_hash (val_max, var_get_width (v_max), 0);
+ hash = value_hash (val_min, var_get_width (v_min), hash);
+ return hash_int (*n_vars * (*n_vars + 1 + idx_max) + idx_min, hash);
+ }
+ return -1u;
+}
+
+/*
+ Make a hash table consisting of struct covariance_accumulators.
+ This allows the accumulation of the elements of a covariance matrix
+ in a single data pass. Call covariance_accumulate () for each case
+ in the data.
+ */
+static struct hsh_table *
+covariance_hsh_create (size_t *n_vars)
+{
+ return hsh_create (*n_vars * *n_vars, covariance_accumulator_compare,
+ covariance_accumulator_hash, covariance_accumulator_free,
+ n_vars);
+}
+
+static void
+covariance_accumulator_free (void *c_, const void *aux UNUSED)
+{
+ struct covariance_accumulator *c = c_;
+ assert (c != NULL);
+ free (c);
+}
+
+static int
+ordered_match_nodes (const struct covariance_accumulator *c, const struct variable *v1,
+ const struct variable *v2, const union value *val1, const union value *val2)
+{
+ size_t result;
+ size_t m;
+
+ result = var_get_dict_index (v1) ^ var_get_dict_index (c->v1);
+ m = var_get_dict_index (v2) ^ var_get_dict_index (c->v2);
+ result = result|m;
+ if (var_is_alpha (v1))
+ {
+ result |= value_compare_3way (val1, c->val1, var_get_width (v1));
+ if (var_is_alpha (v2))
+ {
+ result |= value_compare_3way (val2, c->val2, var_get_width (v2));
+ }
+ }
+ else if (var_is_alpha (v2))
+ {
+ result |= value_compare_3way (val2, c->val2, var_get_width (v2));
+ }
+ return result;
+}
+
+/*
+ Hash comparison. Returns 0 for a match, or a non-zero int
+ otherwise. The sign of a non-zero return value *should* indicate the
+ position of C relative to the covariance_accumulator described by
+ the other arguments. But for now, it just returns 1 for any
+ non-match. This should be changed when someone figures out how to
+ compute a sensible sign for the return value.
+ */
+static int
+match_nodes (const struct covariance_accumulator *c,
+ const struct variable *v1, const struct variable *v2,
+ const union value *val1, const union value *val2)
+{
+ size_t n;
+ size_t m;
+
+ n = ordered_match_nodes (c, v1, v2, val1, val2);
+ m = ordered_match_nodes (c, v2, v1, val2, val1);
+ return (n & m);
+}
+
+/*
+ This function is meant to be used as a comparison function for
+ a struct hsh_table in src/libpspp/hash.c.
+*/
+static int
+covariance_accumulator_compare (const void *a1_, const void *a2_,
+ const void *aux UNUSED)
+{
+ const struct covariance_accumulator *a1 = a1_;
+ const struct covariance_accumulator *a2 = a2_;
+
+ if (a1 == NULL && a2 == NULL)
+ return 0;
+
+ if (a1 == NULL || a2 == NULL)
+ return 1;
+
+ return match_nodes (a1, a2->v1, a2->v2, a2->val1, a2->val2);
+}
+
+static unsigned int
+hash_numeric_alpha (const struct variable *v1, const struct variable *v2,
+ const union value *val, size_t n_vars)
+{
+ unsigned int result = -1u;
+ 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) + value_hash (val, var_get_width (v2), 0);
+ }
+ else if (var_is_alpha (v1) && var_is_numeric (v2))
+ {
+ result = hash_numeric_alpha (v2, v1, val, n_vars);
+ }
+ return result;
+}
+
+
+static double
+update_product (const struct variable *v1, const struct variable *v2,
+ const union value *val1, const union value *val2)
+{
+ assert (v1 != NULL);
+ assert (v2 != NULL);
+ assert (val1 != NULL);
+ assert (val2 != NULL);
+ if (var_is_alpha (v1) && var_is_alpha (v2))
+ {
+ return 1.0;
+ }
+ if (var_is_numeric (v1) && var_is_numeric (v2))
+ {
+ return (val1->f * val2->f);
+ }
+ if (var_is_numeric (v1) && var_is_alpha (v2))
+ {
+ return val1->f;
+ }
+ if (var_is_numeric (v2) && var_is_alpha (v1))
+ {
+ return val2->f;
+ }
+ else
+ {
+ return 0.0;
+ }
+}
+static double
+update_sum (const struct variable *var, const union value *val, double weight)
+{
+ assert (var != NULL);
+ assert (val != NULL);
+ if (var_is_alpha (var))
+ {
+ return weight;
+ }
+ return val->f;
+}
+static struct covariance_accumulator *
+get_new_covariance_accumulator (const struct variable *v1,
+ const struct variable *v2,
+ const union value *val1,
+ const union value *val2)
+{
+ if ((v1 != NULL) && (v2 != NULL) && (val1 != NULL) && (val2 != NULL))
+ {
+ struct covariance_accumulator *ca;
+ ca = xmalloc (sizeof (*ca));
+ ca->v1 = v1;
+ ca->v2 = v2;
+ ca->val1 = val1;
+ ca->val2 = val2;
+ return ca;
+ }
+ 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);
+ }
+ }