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 static const struct variable *
  97 find_matrix_string_var (const struct dictionary *dict, const char *name)
  98 {
  99   const struct variable *var = dict_lookup_var (dict, name);
 100   if (var == NULL)
 101     {
 102       msg (ME, _("Matrix dataset lacks a variable called %s."), name);
 103       return NULL;
 104     }
 105   if (!var_is_alpha (var))
 106     {
 107       msg (ME, _("Matrix dataset variable %s should be of string type."), name);
 108       return NULL;
 109     }
 110   return var;
 111 }
 112
 113 struct matrix_reader *
 114 matrix_reader_create (const struct dictionary *dict,
 115                       struct casereader *in_reader)
 116 {
 117   const struct variable *varname = find_matrix_string_var (dict, "VARNAME_");
 118   const struct variable *rowtype = find_matrix_string_var (dict, "ROWTYPE_");
 119   if (!varname || !rowtype)
 120     return NULL;
 121
 122   for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
 123     {
 124       const struct variable *v = dict_get_var (dict, i);
 125       if (!var_is_numeric (v) && v != rowtype && v != varname)
 126         {
 127           msg (ME, _("Matrix dataset variable %s should be numeric."),
 128                var_get_name (v));
 129           return NULL;
 130         }
 131     }
 132
 133   size_t dvarcnt;
 134   const struct variable **dvars = NULL;
 135   dict_get_vars (dict, &dvars, &dvarcnt, DC_SCRATCH);
 136
 137   /* Continuous variables and split variables. */
 138   const struct variable **cvars = dvars + var_get_dict_index (varname) + 1;
 139   size_t n_cvars = dvarcnt - var_get_dict_index (varname) - 1;
 140   const struct variable **svars = dvars;
 141   size_t n_svars = var_get_dict_index (rowtype);
 142   if (!n_cvars)
 143     {
 144       msg (ME, _("Matrix dataset does not have any continuous variables."));
 145       free (dvars);
 146       return NULL;
 147     }
 148
 149   struct matrix_reader *mr = xmalloc (sizeof *mr);
 150   *mr = (struct matrix_reader) {
 151     .n_cvars = n_cvars,
 152     .cvars = xmemdup (cvars, n_cvars * sizeof *cvars),
 153     .rowtype = rowtype,
 154     .varname = varname,
 155     .dict = dict,
 156     .grouper = casegrouper_create_vars (in_reader, svars, n_svars)
 157   };
 158   free (dvars);
 159
 160   return mr;
 161 }
 162
 163 bool
 164 matrix_reader_destroy (struct matrix_reader *mr)
 165 {
 166   if (mr == NULL)
 167     return false;
 168   bool ret = casegrouper_destroy (mr->grouper);
 169   free (mr);
 170   return ret;
 171 }
 172
 173
 174 /*
 175    Allocates MATRIX if necessary,
 176    and populates row MROW, from the data in C corresponding to
 177    variables in VARS. N_VARS is the length of VARS.
 178 */
 179 static void
 180 matrix_fill_row (gsl_matrix **matrix,
 181       const struct ccase *c, int mrow,
 182       const struct variable **vars, size_t n_vars)
 183 {
 184   int col;
 185   if (*matrix == NULL)
 186     {
 187       *matrix = gsl_matrix_alloc (n_vars, n_vars);
 188       gsl_matrix_set_all (*matrix, SYSMIS);
 189     }
 190
 191   for (col = 0; col < n_vars; ++col)
 192     {
 193       const struct variable *cv = vars [col];
 194       double x = case_num (c, cv);
 195       assert (col  < (*matrix)->size2);
 196       assert (mrow < (*matrix)->size1);
 197       gsl_matrix_set (*matrix, mrow, col, x);
 198     }
 199 }
 200
 201 static int
 202 find_varname (const struct variable **vars, int n_vars,
 203               const char *varname)
 204 {
 205   for (int i = 0; i < n_vars; i++)
 206     if (!strcasecmp (var_get_name (vars[i]), varname))
 207       return i;
 208   return -1;
 209 }
 210
 211 struct substring
 212 matrix_reader_get_string (const struct ccase *c, const struct variable *var)
 213 {
 214   struct substring s = case_ss (c, var);
 215   ss_rtrim (&s, ss_cstr (CC_SPACES));
 216   return s;
 217 }
 218
 219 void
 220 matrix_reader_set_string (struct ccase *c, const struct variable *var,
 221                           struct substring src)
 222 {
 223   struct substring dst = case_ss (c, var);
 224   for (size_t i = 0; i < dst.length; i++)
 225     dst.string[i] = i < src.length ? src.string[i] : ' ';
 226 }
 227
 228 bool
 229 matrix_reader_next (struct matrix_material *mm, struct matrix_reader *mr,
 230                     struct casereader **groupp)
 231 {
 232   struct casereader *group;
 233   if (!casegrouper_get_next_group (mr->grouper, &group))
 234     {
 235       *mm = (struct matrix_material) MATRIX_MATERIAL_INIT;
 236       if (groupp)
 237         *groupp = NULL;
 238       return false;
 239     }
 240
 241   if (groupp)
 242     *groupp = casereader_clone (group);
 243
 244   const struct variable **vars = mr->cvars;
 245   size_t n_vars = mr->n_cvars;
 246
 247   *mm = (struct matrix_material) {
 248     .n = gsl_matrix_calloc (n_vars, n_vars),
 249     .mean_matrix = gsl_matrix_calloc (n_vars, n_vars),
 250     .var_matrix = gsl_matrix_calloc (n_vars, n_vars),
 251   };
 252
 253   struct matrix
 254     {
 255       const char *name;
 256       gsl_matrix **m;
 257       size_t good_rows;
 258       size_t bad_rows;
 259     };
 260   struct matrix matrices[] = {
 261     { .name = "CORR", .m = &mm->corr },
 262     { .name = "COV", .m = &mm->cov },
 263   };
 264   enum { N_MATRICES = 2 };
 265
 266   struct ccase *c;
 267   for (; (c = casereader_read (group)); case_unref (c))
 268     {
 269       struct substring rowtype = matrix_reader_get_string (c, mr->rowtype);
 270
 271       gsl_matrix *v
 272         = (ss_equals_case (rowtype, ss_cstr ("N")) ? mm->n
 273            : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? mm->mean_matrix
 274            : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? mm->var_matrix
 275            : NULL);
 276       if (v)
 277         {
 278           for (int x = 0; x < n_vars; ++x)
 279             {
 280               double n = case_num (c, vars[x]);
 281               if (v == mm->var_matrix)
 282                 n *= n;
 283               for (int y = 0; y < n_vars; ++y)
 284                 gsl_matrix_set (v, y, x, n);
 285             }
 286           continue;
 287         }
 288
 289       struct matrix *m = NULL;
 290       for (size_t i = 0; i < N_MATRICES; i++)
 291         if (ss_equals_case (rowtype, ss_cstr (matrices[i].name)))
 292           {
 293             m = &matrices[i];
 294             break;
 295           }
 296       if (m)
 297         {
 298           struct substring varname_raw = case_ss (c, mr->varname);
 299           struct substring varname = ss_cstr (
 300             recode_string (UTF8, dict_get_encoding (mr->dict),
 301                            varname_raw.string, varname_raw.length));
 302           ss_rtrim (&varname, ss_cstr (CC_SPACES));
 303           varname.string[varname.length] = '\0';
 304
 305           int y = find_varname (vars, n_vars, varname.string);
 306           if (y >= 0)
 307             {
 308               m->good_rows++;
 309               matrix_fill_row (m->m, c, y, vars, n_vars);
 310             }
 311           else
 312             m->bad_rows++;
 313           ss_dealloc (&varname);
 314         }
 315     }
 316   casereader_destroy (group);
 317
 318   for (size_t i = 0; i < N_MATRICES; i++)
 319     if (matrices[i].good_rows && matrices[i].good_rows != n_vars)
 320       msg (SW, _("%s matrix has %zu columns but %zu rows named variables "
 321                  "to be analyzed (and %zu rows named unknown variables)."),
 322            matrices[i].name, n_vars, matrices[i].good_rows,
 323            matrices[i].bad_rows);
 324
 325   return true;
 326 }
 327
 328 int
 329 cmd_debug_matrix_read (struct lexer *lexer UNUSED, struct dataset *ds)
 330 {
 331   struct matrix_reader *mr = matrix_reader_create (dataset_dict (ds),
 332                                                    proc_open (ds));
 333   if (!mr)
 334     return CMD_FAILURE;
 335
 336   struct pivot_table *pt = pivot_table_create ("Debug Matrix Reader");
 337
 338   enum mm_stat
 339     {
 340       MM_CORR,
 341       MM_COV,
 342       MM_N,
 343       MM_MEAN,
 344       MM_STDDEV,
 345     };
 346   const char *mm_stat_names[] = {
 347     [MM_CORR] = "Correlation",
 348     [MM_COV] = "Covariance",
 349     [MM_N] = "N",
 350     [MM_MEAN] = "Mean",
 351     [MM_STDDEV] = "Standard Deviation",
 352   };
 353   enum { N_STATS = sizeof mm_stat_names / sizeof *mm_stat_names };
 354   for (size_t i = 0; i < 2; i++)
 355     {
 356       struct pivot_dimension *d = pivot_dimension_create (
 357         pt,
 358         i ? PIVOT_AXIS_COLUMN : PIVOT_AXIS_ROW,
 359         i ? "Column" : "Row");
 360       if (!i)
 361         pivot_category_create_leaf_rc (d->root, pivot_value_new_text ("Value"),
 362                                        PIVOT_RC_CORRELATION);
 363       for (size_t j = 0; j < mr->n_cvars; j++)
 364         pivot_category_create_leaf_rc (
 365           d->root, pivot_value_new_variable (mr->cvars[j]),
 366           PIVOT_RC_CORRELATION);
 367     }
 368
 369   struct pivot_dimension *stat = pivot_dimension_create (pt, PIVOT_AXIS_ROW,
 370                                                          "Statistic");
 371   for (size_t i = 0; i < N_STATS; i++)
 372     pivot_category_create_leaf (stat->root,
 373                                 pivot_value_new_text (mm_stat_names[i]));
 374
 375   struct pivot_dimension *split = pivot_dimension_create (
 376     pt, PIVOT_AXIS_ROW, "Split");
 377
 378   int split_num = 0;
 379
 380   struct matrix_material mm = MATRIX_MATERIAL_INIT;
 381   while (matrix_reader_next (&mm, mr, NULL))
 382     {
 383       pivot_category_create_leaf (split->root,
 384                                   pivot_value_new_integer (split_num + 1));
 385
 386       const gsl_matrix *m[N_STATS] = {
 387         [MM_CORR] = mm.corr,
 388         [MM_COV] = mm.cov,
 389         [MM_N] = mm.n,
 390         [MM_MEAN] = mm.mean_matrix,
 391         [MM_STDDEV] = mm.var_matrix,
 392       };
 393
 394       for (size_t i = 0; i < N_STATS; i++)
 395         if (m[i])
 396           {
 397             if (i == MM_COV || i == MM_CORR)
 398               {
 399                 for (size_t y = 0; y < mr->n_cvars; y++)
 400                   for (size_t x = 0; x < mr->n_cvars; x++)
 401                     pivot_table_put4 (
 402                       pt, y + 1, x, i, split_num,
 403                       pivot_value_new_number (gsl_matrix_get (m[i], y, x)));
 404               }
 405             else
 406               for (size_t x = 0; x < mr->n_cvars; x++)
 407                 {
 408                   double n = gsl_matrix_get (m[i], 0, x);
 409                   if (i == MM_STDDEV)
 410                     n = sqrt (n);
 411                   pivot_table_put4 (pt, 0, x, i, split_num,
 412                                     pivot_value_new_number (n));
 413                 }
 414           }
 415
 416       split_num++;
 417       matrix_material_uninit (&mm);
 418     }
 419   pivot_table_submit (pt);
 420
 421   proc_commit (ds);
 422
 423   matrix_reader_destroy (mr);
 424   return CMD_SUCCESS;
 425 }