+covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
+{
+ size_t i, j;
+ const double weight = cov->wv ? case_data (c, cov->wv)->f : 1.0;
+
+ assert (cov->passes == 2);
+ assert (cov->state >= 1);
+
+ if (! cov->pass_two_first_case_seen)
+ {
+ size_t m;
+ 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_df_total (cov->categoricals);
+
+ cov->n_cm = (cov->dim * (cov->dim - 1) ) / 2;
+ cov->cm = xcalloc (cov->n_cm, sizeof *cov->cm);
+
+ /* Grow the moment matrices so that they're large enough to accommodate the
+ categorical elements */
+ for (i = 0; i < n_MOMENTS; ++i)
+ {
+ cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
+ }
+
+ /* Populate the moments matrices with the categorical value elements */
+ for (i = cov->n_vars; i < cov->dim; ++i)
+ {
+ for (j = 0 ; j < cov->dim ; ++j) /* FIXME: This is WRONG !!! */
+ {
+ double w = categoricals_get_weight_by_subscript (cov->categoricals, i - cov->n_vars);
+
+ gsl_matrix_set (cov->moments[MOMENT_NONE], i, j, w);
+
+ w = categoricals_get_sum_by_subscript (cov->categoricals, i - cov->n_vars);
+
+ gsl_matrix_set (cov->moments[MOMENT_MEAN], i, j, w);
+ }
+ }
+
+ /* FIXME: This is WRONG!! It must be fixed to properly handle missing values. For
+ now it assumes there are none */
+ for (m = 0 ; m < n_MOMENTS; ++m)
+ {
+ for (i = 0 ; i < cov->dim ; ++i)
+ {
+ double x = gsl_matrix_get (cov->moments[m], i, cov->n_vars -1);
+ for (j = cov->n_vars; j < cov->dim; ++j)
+ {
+ gsl_matrix_set (cov->moments[m], i, j, x);
+ }
+ }
+ }
+
+ /* Divide the means by the number of samples */
+ for (i = 0; i < cov->dim; ++i)
+ {
+ for (j = 0; j < cov->dim; ++j)
+ {
+ double *x = gsl_matrix_ptr (cov->moments[MOMENT_MEAN], i, j);
+ *x /= gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
+ }
+ }
+ }
+
+ for (i = 0 ; i < cov->dim; ++i)
+ {
+ double v1 = get_val (cov, i, c);
+
+ if ( is_missing (cov, i, c))
+ continue;
+
+ for (j = 0 ; j < cov->dim; ++j)
+ {
+ int idx;
+ double ss ;
+ double v2 = get_val (cov, j, c);
+
+ const double s = pow2 (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j)) * weight;
+
+ if ( is_missing (cov, j, c))
+ continue;
+
+ {
+ double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
+ *x += s;
+ }
+
+ ss =
+ (v1 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
+ *
+ (v2 - gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
+ * weight
+ ;
+
+ idx = cm_idx (cov, i, j);
+ if (idx >= 0)
+ {
+ cov->cm [idx] += ss;
+ }
+ }
+ }
+
+ cov->pass_two_first_case_seen = true;
+}
+
+
+/* Call this function for every case in the data set.
+ After all cases have been passed, call covariance_calculate
+ */
+void