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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2017 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 "data/case.h"
  20 #include "data/casereader.h"
  21 #include "data/casewriter.h"
  22 #include "data/dataset.h"
  23 #include "data/dictionary.h"
  24 #include "data/format.h"
  25 #include "data/transformations.h"
  26 #include "data/variable.h"
  27 #include "language/command.h"
  28 #include "language/data-io/data-parser.h"
  29 #include "language/data-io/data-reader.h"
  30 #include "language/data-io/file-handle.h"
  31 #include "language/data-io/inpt-pgm.h"
  32 #include "language/data-io/placement-parser.h"
  33 #include "language/lexer/lexer.h"
  34 #include "language/lexer/variable-parser.h"
  35 #include "libpspp/i18n.h"
  36 #include "libpspp/message.h"
  37 #include "libpspp/misc.h"
  38
  39 #include "gl/xsize.h"
  40 #include "gl/xalloc.h"
  41
  42 #include "gettext.h"
  43 #define _(msgid) gettext (msgid)
  44 \f
  45 /* DATA LIST transformation data. */
  46 struct data_list_trns
  47   {
  48     struct data_parser *parser; /* Parser. */
  49     struct dfm_reader *reader;  /* Data file reader. */
  50     struct variable *end;       /* Variable specified on END subcommand. */
  51   };
  52
  53 static trns_free_func data_list_trns_free;
  54 static trns_proc_func data_list_trns_proc;
  55
  56 enum diagonal
  57   {
  58     DIAGONAL,
  59     NO_DIAGONAL
  60   };
  61
  62 enum triangle
  63   {
  64     LOWER,
  65     UPPER,
  66     FULL
  67   };
  68
  69 static const int ROWTYPE_WIDTH = 8;
  70
  71 struct matrix_format
  72 {
  73   enum triangle triangle;
  74   enum diagonal diagonal;
  75   const struct variable *rowtype;
  76   const struct variable *varname;
  77   int n_continuous_vars;
  78   struct variable **split_vars;
  79   size_t n_split_vars;
  80 };
  81
  82 /*
  83 valid rowtype_ values:
  84   CORR,
  85   COV,
  86   MAT,
  87
  88
  89   MSE,
  90   DFE,
  91   MEAN,
  92   STDDEV (or SD),
  93   N_VECTOR (or N),
  94   N_SCALAR,
  95   N_MATRIX,
  96   COUNT,
  97   PROX.
  98 */
  99
 100 /* Sets the value of OUTCASE which corresponds to VNAME
 101    to the value STR.  VNAME must be of type string.
 102  */
 103 static void
 104 set_varname_column (struct ccase *outcase, const struct variable *vname,
 105      const char *str)
 106 {
 107   int len = var_get_width (vname);
 108   uint8_t *s = value_str_rw (case_data_rw (outcase, vname), len);
 109
 110   strncpy ((char *) s, str, len);
 111 }
 112
 113 static void
 114 blank_varname_column (struct ccase *outcase, const struct variable *vname)
 115 {
 116   int len = var_get_width (vname);
 117   uint8_t *s = value_str_rw (case_data_rw (outcase, vname), len);
 118
 119   memset (s, ' ', len);
 120 }
 121
 122 static struct casereader *
 123 preprocess (struct casereader *casereader0, const struct dictionary *dict, void *aux)
 124 {
 125   struct matrix_format *mformat = aux;
 126   const struct caseproto *proto = casereader_get_proto (casereader0);
 127   struct casewriter *writer;
 128   writer = autopaging_writer_create (proto);
 129   struct ccase *prev_case = NULL;
 130   double **matrices = NULL;
 131   size_t n_splits = 0;
 132
 133   const size_t sizeof_matrix =
 134     sizeof (double) * mformat->n_continuous_vars * mformat->n_continuous_vars;
 135
 136
 137   /* Make an initial pass to populate our temporary matrix */
 138   struct casereader *pass0 = casereader_clone (casereader0);
 139   struct ccase *c;
 140   unsigned int prev_split_hash = 1;
 141   int row = (mformat->triangle == LOWER && mformat->diagonal == NO_DIAGONAL) ? 1 : 0;
 142   for (; (c = casereader_read (pass0)) != NULL; case_unref (c))
 143     {
 144       int s;
 145       unsigned int split_hash = 0;
 146       for (s = 0; s < mformat->n_split_vars; ++s)
 147         {
 148           const struct variable *svar = mformat->split_vars[s];
 149           const union value *sv = case_data (c, svar);
 150           split_hash = value_hash (sv, var_get_width (svar), split_hash);
 151         }
 152
 153       if (matrices == NULL || prev_split_hash != split_hash)
 154         {
 155           row = (mformat->triangle == LOWER && mformat->diagonal == NO_DIAGONAL) ?
 156             1 : 0;
 157
 158           n_splits++;
 159           matrices = xrealloc (matrices, sizeof (double*)  * n_splits);
 160           matrices[n_splits - 1] = xmalloc (sizeof_matrix);
 161         }
 162
 163       prev_split_hash = split_hash;
 164
 165       int c_offset = (mformat->triangle == UPPER) ? row : 0;
 166       if (mformat->triangle == UPPER && mformat->diagonal == NO_DIAGONAL)
 167         c_offset++;
 168       const union value *v = case_data (c, mformat->rowtype);
 169       const char *val = (const char *) value_str (v, ROWTYPE_WIDTH);
 170       if (0 == strncasecmp (val, "corr    ", ROWTYPE_WIDTH) ||
 171           0 == strncasecmp (val, "cov     ", ROWTYPE_WIDTH))
 172         {
 173           if (row >= mformat->n_continuous_vars)
 174             {
 175               msg (SE,
 176                    _("There are %d variable declared but the data has at least %d matrix rows."),
 177                    mformat->n_continuous_vars, row + 1);
 178               goto error;
 179             }
 180           int col;
 181           for (col = c_offset; col < mformat->n_continuous_vars; ++col)
 182             {
 183               const struct variable *var =
 184                 dict_get_var (dict,
 185                               1 + col - c_offset +
 186                               var_get_dict_index (mformat->varname));
 187
 188               double e = case_data (c, var)->f;
 189               if (e == SYSMIS)
 190                 continue;
 191
 192               /* Fill in the lower triangle */
 193               (matrices[n_splits-1])[col + mformat->n_continuous_vars * row] = e;
 194
 195               if (mformat->triangle != FULL)
 196                 /* Fill in the upper triangle */
 197                 (matrices[n_splits-1]) [row + mformat->n_continuous_vars * col] = e;
 198             }
 199           row++;
 200         }
 201     }
 202   casereader_destroy (pass0);
 203
 204   /* Now make a second pass to fill in the other triangle from our
 205      temporary matrix */
 206   const int idx = var_get_dict_index (mformat->varname);
 207   row = 0;
 208
 209   prev_split_hash = 1;
 210   n_splits = 0;
 211   for (; (c = casereader_read (casereader0)) != NULL; prev_case = c)
 212     {
 213       int s;
 214       unsigned int split_hash = 0;
 215       for (s = 0; s < mformat->n_split_vars; ++s)
 216         {
 217           const struct variable *svar = mformat->split_vars[s];
 218           const union value *sv = case_data (c, svar);
 219           split_hash = value_hash (sv, var_get_width (svar), split_hash);
 220         }
 221       if (prev_split_hash != split_hash)
 222         {
 223           n_splits++;
 224           row = 0;
 225         }
 226
 227       prev_split_hash = split_hash;
 228
 229       case_unref (prev_case);
 230       struct ccase *outcase = case_create (proto);
 231       case_copy (outcase, 0, c, 0, caseproto_get_n_widths (proto));
 232       const union value *v = case_data (c, mformat->rowtype);
 233       const char *val = (const char *) value_str (v, ROWTYPE_WIDTH);
 234       if (0 == strncasecmp (val, "corr    ", ROWTYPE_WIDTH) ||
 235           0 == strncasecmp (val, "cov     ", ROWTYPE_WIDTH))
 236         {
 237           int col;
 238           const struct variable *var = dict_get_var (dict, idx + 1 + row);
 239           set_varname_column (outcase, mformat->varname, var_get_name (var));
 240           value_copy (case_data_rw (outcase, mformat->rowtype), v, ROWTYPE_WIDTH);
 241
 242           for (col = 0; col < mformat->n_continuous_vars; ++col)
 243             {
 244               union value *dest_val =
 245                 case_data_rw_idx (outcase,
 246                                   1 + col + var_get_dict_index (mformat->varname));
 247               dest_val->f = (matrices[n_splits - 1])[col + mformat->n_continuous_vars * row];
 248               if (col == row && mformat->diagonal == NO_DIAGONAL)
 249                 dest_val->f = 1.0;
 250             }
 251           row++;
 252         }
 253       else
 254         {
 255           blank_varname_column (outcase, mformat->varname);
 256         }
 257
 258       /* Special case for SD and N_VECTOR: Rewrite as STDDEV and N respectively */
 259       if (0 == strncasecmp (val, "sd      ", ROWTYPE_WIDTH))
 260         {
 261           value_copy_buf_rpad (case_data_rw (outcase, mformat->rowtype), ROWTYPE_WIDTH,
 262                                (uint8_t *) "STDDEV", 6, ' ');
 263         }
 264       else if (0 == strncasecmp (val, "n_vector", ROWTYPE_WIDTH))
 265         {
 266           value_copy_buf_rpad (case_data_rw (outcase, mformat->rowtype), ROWTYPE_WIDTH,
 267                                (uint8_t *) "N", 1, ' ');
 268         }
 269
 270       casewriter_write (writer, outcase);
 271     }
 272
 273   /* If NODIAGONAL is specified, then a final case must be written */
 274   if (mformat->diagonal == NO_DIAGONAL)
 275     {
 276       int col;
 277       struct ccase *outcase = case_create (proto);
 278
 279       if (prev_case)
 280         case_copy (outcase, 0, prev_case, 0, caseproto_get_n_widths (proto));
 281
 282       const struct variable *var = dict_get_var (dict, idx + 1 + row);
 283       set_varname_column (outcase, mformat->varname, var_get_name (var));
 284
 285       for (col = 0; col < mformat->n_continuous_vars; ++col)
 286         {
 287           union value *dest_val =
 288             case_data_rw_idx (outcase, 1 + col +
 289                               var_get_dict_index (mformat->varname));
 290           dest_val->f = (matrices[n_splits - 1]) [col + mformat->n_continuous_vars * row];
 291           if (col == row && mformat->diagonal == NO_DIAGONAL)
 292             dest_val->f = 1.0;
 293         }
 294
 295       casewriter_write (writer, outcase);
 296     }
 297
 298
 299   if (prev_case)
 300     case_unref (prev_case);
 301
 302   int i;
 303   for (i = 0 ; i < n_splits; ++i)
 304     free (matrices[i]);
 305   free (matrices);
 306   struct casereader *reader1 = casewriter_make_reader (writer);
 307   casereader_destroy (casereader0);
 308   return reader1;
 309
 310
 311 error:
 312   if (prev_case)
 313     case_unref (prev_case);
 314
 315   for (i = 0 ; i < n_splits; ++i)
 316     free (matrices[i]);
 317   free (matrices);
 318   casereader_destroy (casereader0);
 319   return NULL;
 320 }
 321
 322 int
 323 cmd_matrix (struct lexer *lexer, struct dataset *ds)
 324 {
 325   struct dictionary *dict;
 326   struct data_parser *parser;
 327   struct dfm_reader *reader;
 328   struct file_handle *fh = NULL;
 329   char *encoding = NULL;
 330   struct matrix_format mformat;
 331   int i;
 332   size_t n_names;
 333   char **names = NULL;
 334
 335   mformat.triangle = LOWER;
 336   mformat.diagonal = DIAGONAL;
 337   mformat.n_split_vars = 0;
 338   mformat.split_vars = NULL;
 339
 340   dict = (in_input_program ()
 341           ? dataset_dict (ds)
 342           : dict_create (get_default_encoding ()));
 343   parser = data_parser_create (dict);
 344   reader = NULL;
 345
 346   data_parser_set_type (parser, DP_DELIMITED);
 347   data_parser_set_warn_missing_fields (parser, false);
 348   data_parser_set_span (parser, false);
 349
 350   mformat.rowtype = dict_create_var (dict, "ROWTYPE_", ROWTYPE_WIDTH);
 351
 352   mformat.n_continuous_vars = 0;
 353   mformat.n_split_vars = 0;
 354
 355   if (! lex_force_match_id (lexer, "VARIABLES"))
 356     goto error;
 357
 358   lex_match (lexer, T_EQUALS);
 359
 360   if (! parse_mixed_vars (lexer, dict, &names, &n_names, PV_NO_DUPLICATE))
 361     {
 362       int i;
 363       for (i = 0; i < n_names; ++i)
 364         free (names[i]);
 365       free (names);
 366       goto error;
 367     }
 368
 369   int longest_name = 0;
 370   for (i = 0; i < n_names; ++i)
 371     {
 372       maximize_int (&longest_name, strlen (names[i]));
 373     }
 374
 375   mformat.varname = dict_create_var (dict, "VARNAME_",
 376                                      8 * DIV_RND_UP (longest_name, 8));
 377
 378   for (i = 0; i < n_names; ++i)
 379     {
 380       if (0 == strcasecmp (names[i], "ROWTYPE_"))
 381         {
 382           const struct fmt_spec fmt = fmt_for_input (FMT_A, 8, 0);
 383           data_parser_add_delimited_field (parser,
 384                                            &fmt,
 385                                            var_get_case_index (mformat.rowtype),
 386                                            "ROWTYPE_");
 387         }
 388       else
 389         {
 390           const struct fmt_spec fmt = fmt_for_input (FMT_F, 10, 4);
 391           struct variable *v = dict_create_var (dict, names[i], 0);
 392           var_set_both_formats (v, &fmt);
 393           data_parser_add_delimited_field (parser,
 394                                            &fmt,
 395                                            var_get_case_index (mformat.varname) +
 396                                            ++mformat.n_continuous_vars,
 397                                            names[i]);
 398         }
 399     }
 400   for (i = 0; i < n_names; ++i)
 401     free (names[i]);
 402   free (names);
 403
 404   while (lex_token (lexer) != T_ENDCMD)
 405     {
 406       if (! lex_force_match (lexer, T_SLASH))
 407         goto error;
 408
 409       if (lex_match_id (lexer, "FORMAT"))
 410         {
 411           lex_match (lexer, T_EQUALS);
 412
 413           while (lex_token (lexer) != T_SLASH && (lex_token (lexer) != T_ENDCMD))
 414             {
 415               if (lex_match_id (lexer, "LIST"))
 416                 {
 417                   data_parser_set_span (parser, false);
 418                 }
 419               else if (lex_match_id (lexer, "FREE"))
 420                 {
 421                   data_parser_set_span (parser, true);
 422                 }
 423               else if (lex_match_id (lexer, "UPPER"))
 424                 {
 425                   mformat.triangle = UPPER;
 426                 }
 427               else if (lex_match_id (lexer, "LOWER"))
 428                 {
 429                   mformat.triangle = LOWER;
 430                 }
 431               else if (lex_match_id (lexer, "FULL"))
 432                 {
 433                   mformat.triangle = FULL;
 434                 }
 435               else if (lex_match_id (lexer, "DIAGONAL"))
 436                 {
 437                   mformat.diagonal = DIAGONAL;
 438                 }
 439               else if (lex_match_id (lexer, "NODIAGONAL"))
 440                 {
 441                   mformat.diagonal = NO_DIAGONAL;
 442                 }
 443               else
 444                 {
 445                   lex_error (lexer, NULL);
 446                   goto error;
 447                 }
 448             }
 449         }
 450       else if (lex_match_id (lexer, "FILE"))
 451         {
 452           lex_match (lexer, T_EQUALS);
 453           fh_unref (fh);
 454           fh = fh_parse (lexer, FH_REF_FILE | FH_REF_INLINE, NULL);
 455           if (fh == NULL)
 456             goto error;
 457         }
 458       else if (lex_match_id (lexer, "SPLIT"))
 459         {
 460           lex_match (lexer, T_EQUALS);
 461           if (! parse_variables (lexer, dict, &mformat.split_vars, &mformat.n_split_vars, 0))
 462             {
 463               free (mformat.split_vars);
 464               goto error;
 465             }
 466           int i;
 467           for (i = 0; i < mformat.n_split_vars; ++i)
 468             {
 469               const struct fmt_spec fmt = fmt_for_input (FMT_F, 4, 0);
 470               var_set_both_formats (mformat.split_vars[i], &fmt);
 471             }
 472           dict_reorder_vars (dict, mformat.split_vars, mformat.n_split_vars);
 473           mformat.n_continuous_vars -= mformat.n_split_vars;
 474         }
 475       else
 476         {
 477           lex_error (lexer, NULL);
 478           goto error;
 479         }
 480     }
 481
 482   if (mformat.diagonal == NO_DIAGONAL && mformat.triangle == FULL)
 483     {
 484       msg (SE, _("FORMAT = FULL and FORMAT = NODIAGONAL are mutually exclusive."));
 485       goto error;
 486     }
 487
 488   if (fh == NULL)
 489     fh = fh_inline_file ();
 490   fh_set_default_handle (fh);
 491
 492   if (!data_parser_any_fields (parser))
 493     {
 494       msg (SE, _("At least one variable must be specified."));
 495       goto error;
 496     }
 497
 498   if (lex_end_of_command (lexer) != CMD_SUCCESS)
 499     goto error;
 500
 501   reader = dfm_open_reader (fh, lexer, encoding);
 502   if (reader == NULL)
 503     goto error;
 504
 505   if (in_input_program ())
 506     {
 507       struct data_list_trns *trns = xmalloc (sizeof *trns);
 508       trns->parser = parser;
 509       trns->reader = reader;
 510       trns->end = NULL;
 511       add_transformation (ds, data_list_trns_proc, data_list_trns_free, trns);
 512     }
 513   else
 514     {
 515       data_parser_make_active_file (parser, ds, reader, dict, preprocess, &mformat);
 516     }
 517
 518   fh_unref (fh);
 519   free (encoding);
 520   free (mformat.split_vars);
 521
 522   return CMD_DATA_LIST;
 523
 524  error:
 525   data_parser_destroy (parser);
 526   if (!in_input_program ())
 527     dict_destroy (dict);
 528   fh_unref (fh);
 529   free (encoding);
 530   free (mformat.split_vars);
 531   return CMD_CASCADING_FAILURE;
 532 }
 533
 534 \f
 535 /* Input procedure. */
 536
 537 /* Destroys DATA LIST transformation TRNS.
 538    Returns true if successful, false if an I/O error occurred. */
 539 static bool
 540 data_list_trns_free (void *trns_)
 541 {
 542   struct data_list_trns *trns = trns_;
 543   data_parser_destroy (trns->parser);
 544   dfm_close_reader (trns->reader);
 545   free (trns);
 546   return true;
 547 }
 548
 549 /* Handle DATA LIST transformation TRNS, parsing data into *C. */
 550 static int
 551 data_list_trns_proc (void *trns_, struct ccase **c, casenumber case_num UNUSED)
 552 {
 553   struct data_list_trns *trns = trns_;
 554   int retval;
 555
 556   *c = case_unshare (*c);
 557   if (data_parser_parse (trns->parser, trns->reader, *c))
 558     retval = TRNS_CONTINUE;
 559   else if (dfm_reader_error (trns->reader) || dfm_eof (trns->reader) > 1)
 560     {
 561       /* An I/O error, or encountering end of file for a second
 562          time, should be escalated into a more serious error. */
 563       retval = TRNS_ERROR;
 564     }
 565   else
 566     retval = TRNS_END_FILE;
 567
 568   /* If there was an END subcommand handle it. */
 569   if (trns->end != NULL)
 570     {
 571       double *end = &case_data_rw (*c, trns->end)->f;
 572       if (retval == TRNS_END_FILE)
 573         {
 574           *end = 1.0;
 575           retval = TRNS_CONTINUE;
 576         }
 577       else
 578         *end = 0.0;
 579     }
 580
 581   return retval;
 582 }
 583