MATRIX DATA: Fully implement.
[pspp] / src / language / data-io / matrix-reader.c
index a7db15e029d8b9feac1d9111a8297edb87c6540b..cf7b524480334fa8a6e5568f5edb8470961b8c99 100644 (file)
 #include "matrix-reader.h"
 
 #include <stdbool.h>
-
-#include <libpspp/message.h>
-#include <libpspp/str.h>
-#include <data/casegrouper.h>
-#include <data/casereader.h>
-#include <data/dictionary.h>
-#include <data/variable.h>
-#include <data/data-out.h>
-#include <data/format.h>
+#include <math.h>
+
+#include "data/casegrouper.h"
+#include "data/casereader.h"
+#include "data/data-out.h"
+#include "data/dataset.h"
+#include "data/dictionary.h"
+#include "data/format.h"
+#include "data/variable.h"
+#include "language/command.h"
+#include "libpspp/i18n.h"
+#include "libpspp/message.h"
+#include "libpspp/str.h"
+#include "output/pivot-table.h"
 
 #include "gettext.h"
 #define _(msgid) gettext (msgid)
 #define N_(msgid) msgid
 
+struct lexer;
 
 /*
 This module interprets a "data matrix", typically generated by the command
@@ -77,16 +83,22 @@ s_0 ROWTYPE_   VARNAME_   v_0         v_1         v_2
 
 */
 
+void
+matrix_material_uninit (struct matrix_material *mm)
+{
+  gsl_matrix_free (mm->corr);
+  gsl_matrix_free (mm->cov);
+  gsl_matrix_free (mm->n);
+  gsl_matrix_free (mm->mean_matrix);
+  gsl_matrix_free (mm->var_matrix);
+}
+\f
 struct matrix_reader
 {
   const struct dictionary *dict;
   const struct variable *varname;
   const struct variable *rowtype;
   struct casegrouper *grouper;
-
-  gsl_matrix *n_vectors;
-  gsl_matrix *mean_vectors;
-  gsl_matrix *var_vectors;
 };
 
 struct matrix_reader *
@@ -177,7 +189,10 @@ matrix_fill_row (gsl_matrix **matrix,
 {
   int col;
   if (*matrix == NULL)
-    *matrix = gsl_matrix_alloc (n_vars, n_vars);
+    {
+      *matrix = gsl_matrix_alloc (n_vars, n_vars);
+      gsl_matrix_set_all (*matrix, SYSMIS);
+    }
 
   for (col = 0; col < n_vars; ++col)
     {
@@ -189,6 +204,16 @@ matrix_fill_row (gsl_matrix **matrix,
     }
 }
 
+static int
+find_varname (const struct variable **vars, int n_vars,
+              const char *varname)
+{
+  for (int i = 0; i < n_vars; i++)
+    if (!strcasecmp (var_get_name (vars[i]), varname))
+      return i;
+  return -1;
+}
+
 bool
 next_matrix_from_reader (struct matrix_material *mm,
                         struct matrix_reader *mr,
@@ -198,87 +223,191 @@ next_matrix_from_reader (struct matrix_material *mm,
 
   assert (vars);
 
-  gsl_matrix_free (mr->n_vectors);
-  gsl_matrix_free (mr->mean_vectors);
-  gsl_matrix_free (mr->var_vectors);
-
   if (!casegrouper_get_next_group (mr->grouper, &group))
-    return false;
-
-  mr->n_vectors    = gsl_matrix_alloc (n_vars, n_vars);
-  mr->mean_vectors = gsl_matrix_alloc (n_vars, n_vars);
-  mr->var_vectors  = gsl_matrix_alloc (n_vars, n_vars);
+    {
+      *mm = (struct matrix_material) MATRIX_MATERIAL_INIT;
+      return false;
+    }
 
-  mm->n = mr->n_vectors;
-  mm->mean_matrix = mr->mean_vectors;
-  mm->var_matrix = mr->var_vectors;
+  *mm = (struct matrix_material) {
+    .n = gsl_matrix_calloc (n_vars, n_vars),
+    .mean_matrix = gsl_matrix_calloc (n_vars, n_vars),
+    .var_matrix = gsl_matrix_calloc (n_vars, n_vars),
+  };
 
-  struct substring *var_names = XCALLOC (n_vars,  struct substring);
-  for (int i = 0; i < n_vars; ++i)
+  struct matrix
     {
-      ss_alloc_substring (var_names + i, ss_cstr (var_get_name (vars[i])));
-    }
+      const char *name;
+      gsl_matrix **m;
+      size_t good_rows;
+      size_t bad_rows;
+    };
+  struct matrix matrices[] = {
+    { .name = "CORR", .m = &mm->corr },
+    { .name = "COV", .m = &mm->cov },
+  };
+  enum { N_MATRICES = 2 };
 
   struct ccase *c;
   for (; (c = casereader_read (group)); case_unref (c))
     {
-      const union value *uv = case_data (c, mr->rowtype);
-      const char *row_type = CHAR_CAST (const char *, uv->s);
-      int col, row;
-      for (col = 0; col < n_vars; ++col)
-       {
-         const struct variable *cv = vars[col];
-         double x = case_data (c, cv)->f;
-         if (0 == strncasecmp (row_type, "N       ", 8))
-           for (row = 0; row < n_vars; ++row)
-             gsl_matrix_set (mr->n_vectors, row, col, x);
-         else if (0 == strncasecmp (row_type, "MEAN    ", 8))
-           for (row = 0; row < n_vars; ++row)
-             gsl_matrix_set (mr->mean_vectors, row, col, x);
-         else if (0 == strncasecmp (row_type, "STDDEV  ", 8))
-           for (row = 0; row < n_vars; ++row)
-             gsl_matrix_set (mr->var_vectors, row, col, x * x);
-       }
-
-      const char *enc = dict_get_encoding (mr->dict);
+      struct substring rowtype = case_ss (c, mr->rowtype);
+      ss_rtrim (&rowtype, ss_cstr (CC_SPACES));
+
+      gsl_matrix *v
+        = (ss_equals_case (rowtype, ss_cstr ("N")) ? mm->n
+           : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? mm->mean_matrix
+           : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? mm->var_matrix
+           : NULL);
+      if (v)
+        {
+          for (int x = 0; x < n_vars; ++x)
+            {
+              double n = case_num (c, vars[x]);
+              if (v == mm->var_matrix)
+                n *= n;
+              for (int y = 0; y < n_vars; ++y)
+                gsl_matrix_set (v, y, x, n);
+            }
+          continue;
+        }
+
+      struct matrix *m = NULL;
+      for (size_t i = 0; i < N_MATRICES; i++)
+        if (ss_equals_case (rowtype, ss_cstr (matrices[i].name)))
+          {
+            m = &matrices[i];
+            break;
+          }
+      if (m)
+        {
+          struct substring varname_raw = case_ss (c, mr->varname);
+          struct substring varname = ss_cstr (
+            recode_string (UTF8, dict_get_encoding (mr->dict),
+                           varname_raw.string, varname_raw.length));
+          ss_rtrim (&varname, ss_cstr (CC_SPACES));
+          varname.string[varname.length] = '\0';
+
+          int y = find_varname (vars, n_vars, varname.string);
+          if (y >= 0)
+            {
+              m->good_rows++;
+              matrix_fill_row (m->m, c, y, vars, n_vars);
+            }
+          else
+            m->bad_rows++;
+          ss_dealloc (&varname);
+        }
+    }
+  casereader_destroy (group);
 
-      const union value *uvv  = case_data (c, mr->varname);
-      int w = var_get_width (mr->varname);
+  for (size_t i = 0; i < N_MATRICES; i++)
+    if (matrices[i].good_rows && matrices[i].good_rows != n_vars)
+      msg (SW, _("%s matrix has %d columns but %zu rows named variables "
+                 "to be analyzed (and %zu rows named unknown variables)."),
+           matrices[i].name, n_vars, matrices[i].good_rows,
+           matrices[i].bad_rows);
 
-      struct fmt_spec fmt = { .type = FMT_A };
-      fmt.w = w;
-      char *vname = data_out (uvv, enc, &fmt, settings_get_fmt_settings ());
-      struct substring the_name = ss_cstr (vname);
+  return true;
+}
 
-      int mrow = -1;
-      for (int i = 0; i < n_vars; ++i)
-       {
-         if (ss_equals (var_names[i], the_name))
-           {
-             mrow = i;
-             break;
-           }
-       }
-      free (vname);
+int
+cmd_debug_matrix_read (struct lexer *lexer UNUSED, struct dataset *ds)
+{
+  const struct variable **vars;
+  size_t n_vars;
+  struct matrix_reader *mr = create_matrix_reader_from_case_reader (
+    dataset_dict (ds), proc_open (ds), &vars, &n_vars);
+  if (!mr)
+    return CMD_FAILURE;
 
-      if (mrow == -1)
-       continue;
+  struct pivot_table *pt = pivot_table_create ("Debug Matrix Reader");
 
-      if (0 == strncasecmp (row_type, "CORR    ", 8))
-       {
-         matrix_fill_row (&mm->corr, c, mrow, vars, n_vars);
-       }
-      else if (0 == strncasecmp (row_type, "COV     ", 8))
-       {
-         matrix_fill_row (&mm->cov, c, mrow, vars, n_vars);
-       }
+  enum mm_stat
+    {
+      MM_CORR,
+      MM_COV,
+      MM_N,
+      MM_MEAN,
+      MM_STDDEV,
+    };
+  const char *mm_stat_names[] = {
+    [MM_CORR] = "Correlation",
+    [MM_COV] = "Covariance",
+    [MM_N] = "N",
+    [MM_MEAN] = "Mean",
+    [MM_STDDEV] = "Standard Deviation",
+  };
+  enum { N_STATS = sizeof mm_stat_names / sizeof *mm_stat_names };
+  for (size_t i = 0; i < 2; i++)
+    {
+      struct pivot_dimension *d = pivot_dimension_create (
+        pt,
+        i ? PIVOT_AXIS_COLUMN : PIVOT_AXIS_ROW,
+        i ? "Column" : "Row");
+      if (!i)
+        pivot_category_create_leaf_rc (d->root, pivot_value_new_text ("Value"),
+                                       PIVOT_RC_CORRELATION);
+      for (size_t j = 0; j < n_vars; j++)
+        pivot_category_create_leaf_rc (
+          d->root, pivot_value_new_variable (vars[j]), PIVOT_RC_CORRELATION);
     }
 
-  casereader_destroy (group);
+  struct pivot_dimension *stat = pivot_dimension_create (pt, PIVOT_AXIS_ROW,
+                                                         "Statistic");
+  for (size_t i = 0; i < N_STATS; i++)
+    pivot_category_create_leaf (stat->root,
+                                pivot_value_new_text (mm_stat_names[i]));
 
-  for (int i = 0; i < n_vars; ++i)
-    ss_dealloc (var_names + i);
-  free (var_names);
+  struct pivot_dimension *split = pivot_dimension_create (
+    pt, PIVOT_AXIS_ROW, "Split");
 
-  return true;
+  int split_num = 0;
+
+  struct matrix_material mm = MATRIX_MATERIAL_INIT;
+  while (next_matrix_from_reader (&mm, mr, vars, n_vars))
+    {
+      pivot_category_create_leaf (split->root,
+                                  pivot_value_new_integer (split_num + 1));
+
+      const gsl_matrix *m[N_STATS] = {
+        [MM_CORR] = mm.corr,
+        [MM_COV] = mm.cov,
+        [MM_N] = mm.n,
+        [MM_MEAN] = mm.mean_matrix,
+        [MM_STDDEV] = mm.var_matrix,
+      };
+
+      for (size_t i = 0; i < N_STATS; i++)
+        if (m[i])
+          {
+            if (i == MM_COV || i == MM_CORR)
+              {
+                for (size_t y = 0; y < n_vars; y++)
+                  for (size_t x = 0; x < n_vars; x++)
+                    pivot_table_put4 (
+                      pt, y + 1, x, i, split_num,
+                      pivot_value_new_number (gsl_matrix_get (m[i], y, x)));
+              }
+            else
+              for (size_t x = 0; x < n_vars; x++)
+                {
+                  double n = gsl_matrix_get (m[i], 0, x);
+                  if (i == MM_STDDEV)
+                    n = sqrt (n);
+                  pivot_table_put4 (pt, 0, x, i, split_num,
+                                    pivot_value_new_number (n));
+                }
+          }
+
+      split_num++;
+      matrix_material_uninit (&mm);
+    }
+  pivot_table_submit (pt);
+
+  proc_commit (ds);
+
+  destroy_matrix_reader (mr);
+  free (vars);
+  return CMD_SUCCESS;
 }