pintos-os.org Git - pspp/blob - src/language/data-io/matrix-reader.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2017, 2019 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include "matrix-reader.h"
  20
  21 #include <stdbool.h>
  22 #include <math.h>
  23
  24 #include "data/casegrouper.h"
  25 #include "data/casereader.h"
  26 #include "data/data-out.h"
  27 #include "data/dataset.h"
  28 #include "data/dictionary.h"
  29 #include "data/format.h"
  30 #include "data/variable.h"
  31 #include "language/command.h"
  32 #include "libpspp/i18n.h"
  33 #include "libpspp/message.h"
  34 #include "libpspp/str.h"
  35 #include "output/pivot-table.h"
  36
  37 #include "gettext.h"
  38 #define _(msgid) gettext (msgid)
  39 #define N_(msgid) msgid
  40
  41 struct lexer;
  42
  43 /*
  44 This module interprets a "data matrix", typically generated by the command
  45 MATRIX DATA.  The dictionary of such a matrix takes the form:
  46
  47  s_0, s_1, ... s_m, ROWTYPE_, VARNAME_, v_0, v_1, .... v_n
  48
  49 where s_0, s_1 ... s_m are the variables defining the splits, and
  50 v_0, v_1 ... v_n are the continuous variables.
  51
  52 m >= 0; n >= 0
  53
  54 The ROWTYPE_ variable is of type A8.
  55 The VARNAME_ variable is a string type whose width is not predetermined.
  56 The variables s_x are of type F4.0 (although this reader accepts any type),
  57 and v_x are of any numeric type.
  58
  59 The values of the ROWTYPE_ variable are in the set {MEAN, STDDEV, N, CORR, COV}
  60 and determine the purpose of that case.
  61 The values of the VARNAME_ variable must correspond to the names of the varibles
  62 in {v_0, v_1 ... v_n} and indicate the rows of the correlation or covariance
  63 matrices.
  64
  65
  66
  67 A typical example is as follows:
  68
  69 s_0 ROWTYPE_   VARNAME_   v_0         v_1         v_2
  70
  71 0   MEAN                5.0000       4.0000       3.0000
  72 0   STDDEV              1.0000       2.0000       3.0000
  73 0   N                   9.0000       9.0000       9.0000
  74 0   CORR       V1       1.0000        .6000        .7000
  75 0   CORR       V2        .6000       1.0000        .8000
  76 0   CORR       V3        .7000        .8000       1.0000
  77 1   MEAN                9.0000       8.0000       7.0000
  78 1   STDDEV              5.0000       6.0000       7.0000
  79 1   N                   9.0000       9.0000       9.0000
  80 1   CORR       V1       1.0000        .4000        .3000
  81 1   CORR       V2        .4000       1.0000        .2000
  82 1   CORR       V3        .3000        .2000       1.0000
  83
  84 */
  85
  86 void
  87 matrix_material_uninit (struct matrix_material *mm)
  88 {
  89   gsl_matrix_free (mm->corr);
  90   gsl_matrix_free (mm->cov);
  91   gsl_matrix_free (mm->n);
  92   gsl_matrix_free (mm->mean_matrix);
  93   gsl_matrix_free (mm->var_matrix);
  94 }
  95 \f
  96 struct matrix_reader
  97 {
  98   const struct dictionary *dict;
  99   const struct variable *varname;
 100   const struct variable *rowtype;
 101   struct casegrouper *grouper;
 102 };
 103
 104 struct matrix_reader *
 105 create_matrix_reader_from_case_reader (const struct dictionary *dict, struct casereader *in_reader,
 106                                        const struct variable ***vars, size_t *n_vars)
 107 {
 108   struct matrix_reader *mr = xzalloc (sizeof *mr);
 109
 110   mr->varname = dict_lookup_var (dict, "varname_");
 111   mr->dict = dict;
 112   if (mr->varname == NULL)
 113     {
 114       msg (ME, _("Matrix dataset lacks a variable called %s."), "VARNAME_");
 115       free (mr);
 116       return NULL;
 117     }
 118
 119   if (!var_is_alpha (mr->varname))
 120     {
 121       msg (ME, _("Matrix dataset variable %s should be of string type."),
 122            "VARNAME_");
 123       free (mr);
 124       return NULL;
 125     }
 126
 127   mr->rowtype = dict_lookup_var (dict, "rowtype_");
 128   if (mr->rowtype == NULL)
 129     {
 130       msg (ME, _("Matrix dataset lacks a variable called %s."), "ROWTYPE_");
 131       free (mr);
 132       return NULL;
 133     }
 134
 135   if (!var_is_alpha (mr->rowtype))
 136     {
 137       msg (ME, _("Matrix dataset variable %s should be of string type."),
 138            "ROWTYPE_");
 139       free (mr);
 140       return NULL;
 141     }
 142
 143   size_t dvarcnt;
 144   const struct variable **dvars = NULL;
 145   dict_get_vars (dict, &dvars, &dvarcnt, DC_SCRATCH);
 146
 147   if (n_vars)
 148     *n_vars = dvarcnt - var_get_dict_index (mr->varname) - 1;
 149
 150   if (vars)
 151     {
 152       int i;
 153       *vars = xcalloc (*n_vars, sizeof (struct variable **));
 154
 155       for (i = 0; i < *n_vars; ++i)
 156         {
 157           (*vars)[i] = dvars[i + var_get_dict_index (mr->varname) + 1];
 158         }
 159     }
 160
 161   /* All the variables before ROWTYPE_ (if any) are split variables */
 162   mr->grouper = casegrouper_create_vars (in_reader, dvars, var_get_dict_index (mr->rowtype));
 163
 164   free (dvars);
 165
 166   return mr;
 167 }
 168
 169 bool
 170 destroy_matrix_reader (struct matrix_reader *mr)
 171 {
 172   if (mr == NULL)
 173     return false;
 174   bool ret = casegrouper_destroy (mr->grouper);
 175   free (mr);
 176   return ret;
 177 }
 178
 179
 180 /*
 181    Allocates MATRIX if necessary,
 182    and populates row MROW, from the data in C corresponding to
 183    variables in VARS. N_VARS is the length of VARS.
 184 */
 185 static void
 186 matrix_fill_row (gsl_matrix **matrix,
 187       const struct ccase *c, int mrow,
 188       const struct variable **vars, size_t n_vars)
 189 {
 190   int col;
 191   if (*matrix == NULL)
 192     {
 193       *matrix = gsl_matrix_alloc (n_vars, n_vars);
 194       gsl_matrix_set_all (*matrix, SYSMIS);
 195     }
 196
 197   for (col = 0; col < n_vars; ++col)
 198     {
 199       const struct variable *cv = vars [col];
 200       double x = case_num (c, cv);
 201       assert (col  < (*matrix)->size2);
 202       assert (mrow < (*matrix)->size1);
 203       gsl_matrix_set (*matrix, mrow, col, x);
 204     }
 205 }
 206
 207 static int
 208 find_varname (const struct variable **vars, int n_vars,
 209               const char *varname)
 210 {
 211   for (int i = 0; i < n_vars; i++)
 212     if (!strcasecmp (var_get_name (vars[i]), varname))
 213       return i;
 214   return -1;
 215 }
 216
 217 bool
 218 next_matrix_from_reader (struct matrix_material *mm,
 219                          struct matrix_reader *mr,
 220                          const struct variable **vars, int n_vars)
 221 {
 222   struct casereader *group;
 223
 224   assert (vars);
 225
 226   if (!casegrouper_get_next_group (mr->grouper, &group))
 227     {
 228       *mm = (struct matrix_material) MATRIX_MATERIAL_INIT;
 229       return false;
 230     }
 231
 232   *mm = (struct matrix_material) {
 233     .n = gsl_matrix_calloc (n_vars, n_vars),
 234     .mean_matrix = gsl_matrix_calloc (n_vars, n_vars),
 235     .var_matrix = gsl_matrix_calloc (n_vars, n_vars),
 236   };
 237
 238   struct matrix
 239     {
 240       const char *name;
 241       gsl_matrix **m;
 242       size_t good_rows;
 243       size_t bad_rows;
 244     };
 245   struct matrix matrices[] = {
 246     { .name = "CORR", .m = &mm->corr },
 247     { .name = "COV", .m = &mm->cov },
 248   };
 249   enum { N_MATRICES = 2 };
 250
 251   struct ccase *c;
 252   for (; (c = casereader_read (group)); case_unref (c))
 253     {
 254       struct substring rowtype = case_ss (c, mr->rowtype);
 255       ss_rtrim (&rowtype, ss_cstr (CC_SPACES));
 256
 257       gsl_matrix *v
 258         = (ss_equals_case (rowtype, ss_cstr ("N")) ? mm->n
 259            : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? mm->mean_matrix
 260            : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? mm->var_matrix
 261            : NULL);
 262       if (v)
 263         {
 264           for (int x = 0; x < n_vars; ++x)
 265             {
 266               double n = case_num (c, vars[x]);
 267               if (v == mm->var_matrix)
 268                 n *= n;
 269               for (int y = 0; y < n_vars; ++y)
 270                 gsl_matrix_set (v, y, x, n);
 271             }
 272           continue;
 273         }
 274
 275       struct matrix *m = NULL;
 276       for (size_t i = 0; i < N_MATRICES; i++)
 277         if (ss_equals_case (rowtype, ss_cstr (matrices[i].name)))
 278           {
 279             m = &matrices[i];
 280             break;
 281           }
 282       if (m)
 283         {
 284           struct substring varname_raw = case_ss (c, mr->varname);
 285           struct substring varname = ss_cstr (
 286             recode_string (UTF8, dict_get_encoding (mr->dict),
 287                            varname_raw.string, varname_raw.length));
 288           ss_rtrim (&varname, ss_cstr (CC_SPACES));
 289           varname.string[varname.length] = '\0';
 290
 291           int y = find_varname (vars, n_vars, varname.string);
 292           if (y >= 0)
 293             {
 294               m->good_rows++;
 295               matrix_fill_row (m->m, c, y, vars, n_vars);
 296             }
 297           else
 298             m->bad_rows++;
 299           ss_dealloc (&varname);
 300         }
 301     }
 302   casereader_destroy (group);
 303
 304   for (size_t i = 0; i < N_MATRICES; i++)
 305     if (matrices[i].good_rows && matrices[i].good_rows != n_vars)
 306       msg (SW, _("%s matrix has %d columns but %zu rows named variables "
 307                  "to be analyzed (and %zu rows named unknown variables)."),
 308            matrices[i].name, n_vars, matrices[i].good_rows,
 309            matrices[i].bad_rows);
 310
 311   return true;
 312 }
 313
 314 int
 315 cmd_debug_matrix_read (struct lexer *lexer UNUSED, struct dataset *ds)
 316 {
 317   const struct variable **vars;
 318   size_t n_vars;
 319   struct matrix_reader *mr = create_matrix_reader_from_case_reader (
 320     dataset_dict (ds), proc_open (ds), &vars, &n_vars);
 321   if (!mr)
 322     return CMD_FAILURE;
 323
 324   struct pivot_table *pt = pivot_table_create ("Debug Matrix Reader");
 325
 326   enum mm_stat
 327     {
 328       MM_CORR,
 329       MM_COV,
 330       MM_N,
 331       MM_MEAN,
 332       MM_STDDEV,
 333     };
 334   const char *mm_stat_names[] = {
 335     [MM_CORR] = "Correlation",
 336     [MM_COV] = "Covariance",
 337     [MM_N] = "N",
 338     [MM_MEAN] = "Mean",
 339     [MM_STDDEV] = "Standard Deviation",
 340   };
 341   enum { N_STATS = sizeof mm_stat_names / sizeof *mm_stat_names };
 342   for (size_t i = 0; i < 2; i++)
 343     {
 344       struct pivot_dimension *d = pivot_dimension_create (
 345         pt,
 346         i ? PIVOT_AXIS_COLUMN : PIVOT_AXIS_ROW,
 347         i ? "Column" : "Row");
 348       if (!i)
 349         pivot_category_create_leaf_rc (d->root, pivot_value_new_text ("Value"),
 350                                        PIVOT_RC_CORRELATION);
 351       for (size_t j = 0; j < n_vars; j++)
 352         pivot_category_create_leaf_rc (
 353           d->root, pivot_value_new_variable (vars[j]), PIVOT_RC_CORRELATION);
 354     }
 355
 356   struct pivot_dimension *stat = pivot_dimension_create (pt, PIVOT_AXIS_ROW,
 357                                                          "Statistic");
 358   for (size_t i = 0; i < N_STATS; i++)
 359     pivot_category_create_leaf (stat->root,
 360                                 pivot_value_new_text (mm_stat_names[i]));
 361
 362   struct pivot_dimension *split = pivot_dimension_create (
 363     pt, PIVOT_AXIS_ROW, "Split");
 364
 365   int split_num = 0;
 366
 367   struct matrix_material mm = MATRIX_MATERIAL_INIT;
 368   while (next_matrix_from_reader (&mm, mr, vars, n_vars))
 369     {
 370       pivot_category_create_leaf (split->root,
 371                                   pivot_value_new_integer (split_num + 1));
 372
 373       const gsl_matrix *m[N_STATS] = {
 374         [MM_CORR] = mm.corr,
 375         [MM_COV] = mm.cov,
 376         [MM_N] = mm.n,
 377         [MM_MEAN] = mm.mean_matrix,
 378         [MM_STDDEV] = mm.var_matrix,
 379       };
 380
 381       for (size_t i = 0; i < N_STATS; i++)
 382         if (m[i])
 383           {
 384             if (i == MM_COV || i == MM_CORR)
 385               {
 386                 for (size_t y = 0; y < n_vars; y++)
 387                   for (size_t x = 0; x < n_vars; x++)
 388                     pivot_table_put4 (
 389                       pt, y + 1, x, i, split_num,
 390                       pivot_value_new_number (gsl_matrix_get (m[i], y, x)));
 391               }
 392             else
 393               for (size_t x = 0; x < n_vars; x++)
 394                 {
 395                   double n = gsl_matrix_get (m[i], 0, x);
 396                   if (i == MM_STDDEV)
 397                     n = sqrt (n);
 398                   pivot_table_put4 (pt, 0, x, i, split_num,
 399                                     pivot_value_new_number (n));
 400                 }
 401           }
 402
 403       split_num++;
 404       matrix_material_uninit (&mm);
 405     }
 406   pivot_table_submit (pt);
 407
 408   proc_commit (ds);
 409
 410   destroy_matrix_reader (mr);
 411   free (vars);
 412   return CMD_SUCCESS;
 413 }