X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;ds=sidebyside;f=src%2Fmath%2Fcovariance.c;h=61f54f50752704e7e9717378574930f3a01b4752;hb=263613c4fe117bcf0702cc7d0d334068b2424d37;hp=d60d609e45681ffc7858abb83490b0b07b09fb46;hpb=a741cf4e8348e478ec5e8658cebeb80dc0376b1a;p=pspp

diff --git a/src/math/covariance.c b/src/math/covariance.c
index d60d609e45..61f54f5075 100644
--- a/src/math/covariance.c
+++ b/src/math/covariance.c
@@ -135,7 +135,7 @@ covariance_1pass_create (size_t n_vars, const struct variable **vars,
 			 const struct variable *weight, enum mv_class exclude)
 {
   size_t i;
-  struct covariance *cov = xmalloc (sizeof *cov);
+  struct covariance *cov = xzalloc (sizeof *cov);
 
   cov->passes = 1;
   cov->state = 0;
@@ -156,7 +156,9 @@ covariance_1pass_create (size_t n_vars, const struct variable **vars,
 
   cov->n_cm = (n_vars * (n_vars - 1)  ) / 2;
 
-  cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
+  if (cov->n_cm > 0)
+    cov->cm = xcalloc (sizeof *cov->cm, cov->n_cm);
+  cov->categoricals = NULL;
 
   return cov;
 }
@@ -265,6 +267,7 @@ get_val (const struct covariance *cov, int i, const struct ccase *c)
   return categoricals_get_binary_by_subscript (cov->categoricals, i - cov->n_vars, c);
 }
 
+#if 0
 void
 dump_matrix (const gsl_matrix *m)
 {
@@ -277,6 +280,7 @@ dump_matrix (const gsl_matrix *m)
       printf ("\n");
     }
 }
+#endif
 
 /* Call this function for every case in the data set */
 void
@@ -589,7 +593,6 @@ covariance_calculate_single_pass (struct covariance *cov)
 }
 
 
-
 /* 
    Return a pointer to gsl_matrix containing the pairwise covariances.
    The matrix remains owned by the COV object, and must not be freed.
@@ -613,6 +616,85 @@ covariance_calculate (struct covariance *cov)
     }
 }
 
+/*
+  Covariance computed without dividing by the sample size.
+ */
+static const 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);
+}
+
+static const gsl_matrix *
+covariance_calculate_single_pass_unnormalized (struct covariance *cov)
+{
+  size_t i, j;
+
+  for (i = 0 ; i < cov->dim; ++i)
+    {
+      for (j = 0 ; j < cov->dim; ++j)
+	{
+	  double *x = gsl_matrix_ptr (cov->moments[MOMENT_VARIANCE], i, j);
+	  *x -= pow2 (gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j))
+	    / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
+	}
+    }
+  for ( j = 0 ; j < cov->dim - 1; ++j)
+    {
+      for (i = j + 1 ; i < cov->dim; ++i)
+	{
+	  double *x = &cov->cm [cm_idx (cov, i, j)];
+	  
+	  *x -=
+	    gsl_matrix_get (cov->moments[MOMENT_MEAN], i, j) 
+	    *
+	    gsl_matrix_get (cov->moments[MOMENT_MEAN], j, i) 
+	  / gsl_matrix_get (cov->moments[MOMENT_NONE], i, j);
+	}
+    }
+
+  return cm_to_gsl (cov);
+}
+
+
+/* 
+   Return a pointer to gsl_matrix containing the pairwise covariances.
+   The matrix remains owned by the COV object, and must not be freed.
+   Call this function only after all data have been accumulated.
+*/
+const gsl_matrix *
+covariance_calculate_unnormalized (struct covariance *cov)
+{
+  assert ( cov->state > 0 );
+
+  switch (cov->passes)
+    {
+    case 1:
+      return covariance_calculate_single_pass_unnormalized (cov);  
+      break;
+    case 2:
+      return covariance_calculate_double_pass_unnormalized (cov);  
+      break;
+    default:
+      NOT_REACHED ();
+    }
+}
 
 
 
@@ -621,7 +703,7 @@ void
 covariance_destroy (struct covariance *cov)
 {
   size_t i;
-  free (cov->vars);
+
   categoricals_destroy (cov->categoricals);
 
   for (i = 0; i < n_MOMENTS; ++i)