X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Fcovariance.c;h=3ec90efbba05188c68bddd54ce53bc2a8e7a5d0f;hb=afa29b53daef68497af3b44291c0d93c9f85a1d5;hp=dc3166923d7ddff670d903274caef1dfe548beb4;hpb=f550aee00a62fe1d8baf62d83cd7efef6cc2ee92;p=pspp

diff --git a/src/math/covariance.c b/src/math/covariance.c
index dc3166923d..3ec90efbba 100644
--- a/src/math/covariance.c
+++ b/src/math/covariance.c
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2009 Free Software Foundation, Inc.
+   Copyright (C) 2009, 2010, 2011 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -16,15 +16,19 @@
 
 #include <config.h>
 
-#include <libpspp/assertion.h>
-#include "covariance.h"
-#include <gl/xalloc.h>
-#include "moments.h"
+#include "math/covariance.h"
+
 #include <gsl/gsl_matrix.h>
-#include <data/case.h>
-#include <data/variable.h>
-#include <libpspp/misc.h>
-#include "categoricals.h"
+
+#include "data/case.h"
+#include "data/variable.h"
+#include "libpspp/assertion.h"
+#include "libpspp/misc.h"
+#include "math/categoricals.h"
+#include "math/interaction.h"
+#include "math/moments.h"
+
+#include "gl/xalloc.h"
 
 #define n_MOMENTS (MOMENT_VARIANCE + 1)
 
@@ -68,7 +72,7 @@ struct covariance
 {
   /* The variables for which the covariance matrix is to be calculated. */
   size_t n_vars;
-  const struct variable **vars;
+  const struct variable *const *vars;
 
   /* Categorical variables. */
   struct categoricals *categoricals;
@@ -109,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;
 };
 
 
@@ -131,11 +137,11 @@ covariance_moments (const struct covariance *cov, int m)
 /* Create a covariance struct.
  */
 struct covariance *
-covariance_1pass_create (size_t n_vars, const struct variable **vars,
+covariance_1pass_create (size_t n_vars, const struct variable *const *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 +162,8 @@ 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);
+
+  cov->cm = xcalloc (cov->n_cm, sizeof *cov->cm);
   cov->categoricals = NULL;
 
   return cov;
@@ -169,8 +176,8 @@ covariance_1pass_create (size_t n_vars, const struct variable **vars,
   until then.
  */
 struct covariance *
-covariance_2pass_create (size_t n_vars, const struct variable **vars,
-			 size_t n_catvars, const struct variable **catvars, 
+covariance_2pass_create (size_t n_vars, const struct variable *const *vars,
+			 struct categoricals *cats,
 			 const struct variable *wv, enum mv_class exclude)
 {
   size_t i;
@@ -196,7 +203,8 @@ covariance_2pass_create (size_t n_vars, const struct variable **vars,
   cov->n_cm = -1;
   cov->cm = NULL;
 
-  cov->categoricals = categoricals_create (catvars, n_catvars, wv, exclude);
+  cov->categoricals = cats;
+  cov->unnormalised = NULL;
 
   return cov;
 }
@@ -243,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);
 
@@ -263,9 +271,10 @@ 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
 void
 dump_matrix (const gsl_matrix *m)
 {
@@ -278,6 +287,7 @@ dump_matrix (const gsl_matrix *m)
       printf ("\n");
     }
 }
+#endif
 
 /* Call this function for every case in the data set */
 void
@@ -293,7 +303,8 @@ covariance_accumulate_pass1 (struct covariance *cov, const struct ccase *c)
       cov->state = 1;
     }
 
-  categoricals_update (cov->categoricals, c);
+  if (cov->categoricals)
+    categoricals_update (cov->categoricals, c);
 
   for (i = 0 ; i < cov->dim; ++i)
     {
@@ -339,11 +350,16 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
       assert (cov->state == 1);
       cov->state = 2;
 
-      cov->dim = cov->n_vars +
-	categoricals_total (cov->categoricals) - categoricals_get_n_variables (cov->categoricals);
+      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 (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 */
@@ -352,8 +368,6 @@ covariance_accumulate_pass2 (struct covariance *cov, const struct ccase *c)
 	  cov->moments[i] = resize_matrix (cov->moments[i], cov->dim);
 	}
 
-      categoricals_done (cov->categoricals);
-
       /* Populate the moments matrices with the categorical value elements */
       for (i = cov->n_vars; i < cov->dim; ++i)
 	{
@@ -521,7 +535,7 @@ cm_to_gsl (struct covariance *cov)
 }
 
 
-static const gsl_matrix *
+static gsl_matrix *
 covariance_calculate_double_pass (struct covariance *cov)
 {
   size_t i, j;
@@ -545,7 +559,7 @@ covariance_calculate_double_pass (struct covariance *cov)
   return  cm_to_gsl (cov);
 }
 
-static const gsl_matrix *
+static gsl_matrix *
 covariance_calculate_single_pass (struct covariance *cov)
 {
   size_t i, j;
@@ -590,15 +604,16 @@ 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.
-   Call this function only after all data have been accumulated.
-*/
-const gsl_matrix *
+/* 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.
+
+   Call this function only after all data have been accumulated.  */
+gsl_matrix *
 covariance_calculate (struct covariance *cov)
 {
-  assert ( cov->state > 0 );
+  if ( cov->state <= 0 )
+    return NULL;
 
   switch (cov->passes)
     {
@@ -616,33 +631,16 @@ covariance_calculate (struct covariance *cov)
 /*
   Covariance computed without dividing by the sample size.
  */
-static const gsl_matrix *
+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);
 }
 
-static const gsl_matrix *
+static gsl_matrix *
 covariance_calculate_single_pass_unnormalized (struct covariance *cov)
 {
   size_t i, j;
-  size_t m;
 
   for (i = 0 ; i < cov->dim; ++i)
     {
@@ -671,29 +669,42 @@ covariance_calculate_single_pass_unnormalized (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.
-   Call this function only after all data have been accumulated.
-*/
+/* Return a pointer to gsl_matrix containing the pairwise covariances.  The
+   returned matrix is owned by the structure, 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 );
+  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
+   The return value is owned by the COV
+*/
+const struct categoricals *
+covariance_get_categoricals (const struct covariance *cov)
+{
+  return cov->categoricals;
+}
 
 
 /* Destroy the COV object */
@@ -707,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);
+    }
+}