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 <gsl/gsl_matrix.h>
  20 #include <gsl/gsl_vector.h>
  21
  22 #include "data/case.h"
  23 #include "data/casereader.h"
  24 #include "data/casewriter.h"
  25 #include "data/data-in.h"
  26 #include "data/dataset.h"
  27 #include "data/dictionary.h"
  28 #include "data/format.h"
  29 #include "data/short-names.h"
  30 #include "data/transformations.h"
  31 #include "data/variable.h"
  32 #include "language/command.h"
  33 #include "language/data-io/data-parser.h"
  34 #include "language/data-io/data-reader.h"
  35 #include "language/data-io/file-handle.h"
  36 #include "language/data-io/inpt-pgm.h"
  37 #include "language/data-io/placement-parser.h"
  38 #include "language/lexer/lexer.h"
  39 #include "language/lexer/variable-parser.h"
  40 #include "libpspp/assertion.h"
  41 #include "libpspp/i18n.h"
  42 #include "libpspp/message.h"
  43 #include "libpspp/str.h"
  44
  45 #include "gl/c-ctype.h"
  46 #include "gl/minmax.h"
  47 #include "gl/xsize.h"
  48 #include "gl/xalloc.h"
  49
  50 #include "gettext.h"
  51 #define _(msgid) gettext (msgid)
  52 \f
  53 #define ROWTYPES                                \
  54     /* Matrix row types. */                     \
  55     RT(CORR,     2)                             \
  56     RT(COV,      2)                             \
  57     RT(MAT,      2)                             \
  58     RT(N_MATRIX, 2)                             \
  59     RT(PROX,     2)                             \
  60                                                 \
  61     /* Vector row types. */                     \
  62     RT(COUNT,    1)                             \
  63     RT(DFE,      1)                             \
  64     RT(MEAN,     1)                             \
  65     RT(MSE,      1)                             \
  66     RT(STDDEV,   1)                             \
  67     RT(N, 1)                                    \
  68                                                 \
  69     /* Scalar row types. */                     \
  70     RT(N_SCALAR, 0)
  71
  72 enum rowtype
  73   {
  74 #define RT(NAME, DIMS) C_##NAME,
  75     ROWTYPES
  76 #undef RT
  77   };
  78
  79 enum
  80   {
  81 #define RT(NAME, DIMS) +1
  82     N_ROWTYPES = ROWTYPES
  83 #undef RT
  84   };
  85 verify (N_ROWTYPES < 32);
  86
  87 /* Returns the number of dimensions in the indexes for row type RT.  A matrix
  88    has 2 dimensions, a vector has 1, a scalar has 0. */
  89 static int
  90 rowtype_dimensions (enum rowtype rt)
  91 {
  92   static const int rowtype_dims[N_ROWTYPES] = {
  93 #define RT(NAME, DIMS) [C_##NAME] = DIMS,
  94     ROWTYPES
  95 #undef RT
  96   };
  97   return rowtype_dims[rt];
  98 }
  99
 100 static struct substring
 101 rowtype_name (enum rowtype rt)
 102 {
 103   static const struct substring rowtype_names[N_ROWTYPES] = {
 104 #define RT(NAME, DIMS) [C_##NAME] = SS_LITERAL_INITIALIZER (#NAME),
 105     ROWTYPES
 106 #undef RT
 107   };
 108
 109   return rowtype_names[rt];
 110 }
 111
 112 static bool
 113 rowtype_from_string (struct substring token, enum rowtype *rt)
 114 {
 115   ss_trim (&token, ss_cstr (CC_SPACES));
 116   for (size_t i = 0; i < N_ROWTYPES; i++)
 117     if (lex_id_match (rowtype_name (i), token))
 118       {
 119         *rt = i;
 120         return true;
 121       }
 122
 123   if (lex_id_match (ss_cstr ("N_VECTOR"), token))
 124     {
 125       *rt = C_N;
 126       return true;
 127     }
 128   else if (lex_id_match (ss_cstr ("SD"), token))
 129     {
 130       *rt = C_STDDEV;
 131       return true;
 132     }
 133
 134   return false;
 135 }
 136
 137 static bool
 138 rowtype_parse (struct lexer *lexer, enum rowtype *rt)
 139 {
 140   bool parsed = (lex_token (lexer) == T_ID
 141                  && rowtype_from_string (lex_tokss (lexer), rt));
 142   if (parsed)
 143     lex_get (lexer);
 144   return parsed;
 145 }
 146 \f
 147 struct matrix_format
 148   {
 149     bool span;
 150     enum triangle
 151       {
 152         LOWER,
 153         UPPER,
 154         FULL
 155       }
 156     triangle;
 157     enum diagonal
 158       {
 159         DIAGONAL,
 160         NO_DIAGONAL
 161       }
 162     diagonal;
 163
 164     bool input_rowtype;
 165     struct variable **input_vars;
 166     size_t n_input_vars;
 167
 168     /* How to read matrices with each possible number of dimensions (0=scalar,
 169        1=vector, 2=matrix). */
 170     struct matrix_sched
 171       {
 172         /* Number of rows and columns in the matrix: (1,1) for a scalar, (1,n) for
 173            a vector, (n,n) for a matrix. */
 174         int nr, nc;
 175
 176         /* Rows of data to read and the number of columns in each.  Because we
 177            often read just a triangle and sometimes omit the diagonal, 'n_rp' can
 178            be less than 'nr' and 'rp[i]->y' isn't always 'y'. */
 179         struct row_sched
 180           {
 181             /* The y-value of the row inside the matrix. */
 182             int y;
 183
 184             /* first and last (exclusive) columns to read in this row. */
 185             int x0, x1;
 186           }
 187           *rp;
 188         size_t n_rp;
 189       }
 190     ms[3];
 191
 192     struct variable *rowtype;
 193     struct variable *varname;
 194     struct variable **cvars;
 195     int n_cvars;
 196     struct variable **svars;
 197     size_t *svar_indexes;
 198     size_t n_svars;
 199     struct variable **fvars;
 200     size_t *fvar_indexes;
 201     size_t n_fvars;
 202     int cells;
 203     int n;
 204
 205     unsigned int pooled_rowtype_mask;
 206     unsigned int factor_rowtype_mask;
 207
 208     struct content
 209       {
 210         bool open;
 211         enum rowtype rowtype;
 212         bool close;
 213       }
 214     *contents;
 215     size_t n_contents;
 216   };
 217
 218 static void
 219 matrix_format_uninit (struct matrix_format *mf)
 220 {
 221   free (mf->input_vars);
 222   for (int i = 0; i < 3; i++)
 223     free (mf->ms[i].rp);
 224   free (mf->cvars);
 225   free (mf->svars);
 226   free (mf->svar_indexes);
 227   free (mf->fvars);
 228   free (mf->fvar_indexes);
 229   free (mf->contents);
 230 }
 231
 232 static void
 233 set_string (struct ccase *outcase, const struct variable *var,
 234             struct substring src)
 235 {
 236   struct substring dst = case_ss (outcase, var);
 237   for (size_t i = 0; i < dst.length; i++)
 238     dst.string[i] = i < src.length ? src.string[i] : ' ';
 239 }
 240
 241 static void
 242 parse_msg (struct dfm_reader *reader, const struct substring *token,
 243            char *text, enum msg_severity severity)
 244 {
 245   int first_column = 0;
 246   if (token)
 247     {
 248       struct substring line = dfm_get_record (reader);
 249       if (token->string >= line.string && token->string < ss_end (line))
 250         first_column = ss_pointer_to_position (line, token->string) + 1;
 251     }
 252
 253   int line_number = dfm_get_line_number (reader);
 254   struct msg_location *location = xmalloc (sizeof *location);
 255   *location = (struct msg_location) {
 256     .file_name = xstrdup (dfm_get_file_name (reader)),
 257     .first_line = line_number,
 258     .last_line = line_number + 1,
 259     .first_column = first_column,
 260     .last_column = first_column ? first_column + token->length : 0,
 261   };
 262   struct msg *m = xmalloc (sizeof *m);
 263   *m = (struct msg) {
 264     .category = MSG_C_DATA,
 265     .severity = severity,
 266     .location = location,
 267     .text = text,
 268   };
 269   msg_emit (m);
 270 }
 271
 272 static void PRINTF_FORMAT (3, 4)
 273 parse_warning (struct dfm_reader *reader, const struct substring *token,
 274                const char *format, ...)
 275 {
 276   va_list args;
 277   va_start (args, format);
 278   parse_msg (reader, token, xvasprintf (format, args), MSG_S_WARNING);
 279   va_end (args);
 280 }
 281
 282 static void PRINTF_FORMAT (3, 4)
 283 parse_error (struct dfm_reader *reader, const struct substring *token,
 284              const char *format, ...)
 285 {
 286   va_list args;
 287   va_start (args, format);
 288   parse_msg (reader, token, xvasprintf (format, args), MSG_S_ERROR);
 289   va_end (args);
 290 }
 291
 292 /* Advance to beginning of next token. */
 293 static bool
 294 more_tokens (struct substring *p, struct dfm_reader *r)
 295 {
 296   for (;;)
 297     {
 298       ss_ltrim (p, ss_cstr (CC_SPACES ","));
 299       if (p->length)
 300         return true;
 301
 302       dfm_forward_record (r);
 303       if (dfm_eof (r))
 304         return false;
 305       *p = dfm_get_record (r);
 306     }
 307 }
 308
 309 static bool
 310 next_token (struct substring *p, struct dfm_reader *r, struct substring *token)
 311 {
 312   if (!more_tokens (p, r))
 313     return false;
 314
 315   /* Collect token. */
 316   int c = ss_first (*p);
 317   if (c == '\'' || c == '"')
 318     {
 319       ss_advance (p, 1);
 320       ss_get_until (p, c, token);
 321     }
 322   else
 323     {
 324       size_t n = 1;
 325       for (;;)
 326         {
 327           c = ss_at (*p, n);
 328           if (c == EOF
 329               || ss_find_byte (ss_cstr (CC_SPACES ","), c) != SIZE_MAX
 330               || ((c == '+' || c == '-')
 331                   && ss_find_byte (ss_cstr ("dDeE"),
 332                                    ss_at (*p, n - 1)) == SIZE_MAX))
 333             break;
 334           n++;
 335         }
 336       ss_get_bytes (p, n, token);
 337     }
 338   return true;
 339 }
 340
 341 static bool
 342 next_number (struct substring *p, struct dfm_reader *r, double *d)
 343 {
 344   struct substring token;
 345   if (!next_token (p, r, &token))
 346     return false;
 347
 348   union value v;
 349   char *error = data_in (token, dfm_reader_get_encoding (r), FMT_F,
 350                          settings_get_fmt_settings (), &v, 0, NULL);
 351   if (error)
 352     {
 353       parse_error (r, &token, "%s", error);
 354       free (error);
 355     }
 356   *d = v.f;
 357   return true;
 358 }
 359
 360 static bool
 361 next_rowtype (struct substring *p, struct dfm_reader *r, enum rowtype *rt)
 362 {
 363   struct substring token;
 364   if (!next_token (p, r, &token))
 365     return false;
 366
 367   if (rowtype_from_string (token, rt))
 368     return true;
 369
 370   parse_error (r, &token, _("Unknown row type \"%.*s\"."),
 371                (int) token.length, token.string);
 372   return false;
 373 }
 374
 375 struct read_matrix_params
 376   {
 377     /* Adjustments to first and last row to read. */
 378     int dy0, dy1;
 379
 380     /* Left and right columns to read in first row, inclusive.
 381        For x1, INT_MAX is the rightmost column. */
 382     int x0, x1;
 383
 384     /* Adjustment to x0 and x1 for each subsequent row we read.  Each of these
 385        is 0 to keep it the same or -1 or +1 to adjust it by that much. */
 386     int dx0, dx1;
 387   };
 388
 389 static const struct read_matrix_params *
 390 get_read_matrix_params (const struct matrix_format *mf)
 391 {
 392   if (mf->triangle == FULL)
 393     {
 394       /* 1 2 3 4
 395          2 1 5 6
 396          3 5 1 7
 397          4 6 7 1 */
 398       static const struct read_matrix_params rmp = { 0, 0, 0, INT_MAX, 0, 0 };
 399       return &rmp;
 400     }
 401   else if (mf->triangle == LOWER)
 402     {
 403       if (mf->diagonal == DIAGONAL)
 404         {
 405           /* 1 . . .
 406              2 1 . .
 407              3 5 1 .
 408              4 6 7 1 */
 409           static const struct read_matrix_params rmp = { 0, 0, 0, 0, 0, 1 };
 410           return &rmp;
 411         }
 412       else
 413         {
 414           /* . . . .
 415              2 . . .
 416              3 5 . .
 417              4 6 7 . */
 418           static const struct read_matrix_params rmp = { 1, 0, 0, 0, 0, 1 };
 419           return &rmp;
 420         }
 421     }
 422   else if (mf->triangle == UPPER)
 423     {
 424       if (mf->diagonal == DIAGONAL)
 425         {
 426           /* 1 2 3 4
 427              . 1 5 6
 428              . . 1 7
 429              . . . 1 */
 430           static const struct read_matrix_params rmp = { 0, 0, 0, INT_MAX, 1, 0 };
 431           return &rmp;
 432         }
 433       else
 434         {
 435           /* . 2 3 4
 436              . . 5 6
 437              . . . 7
 438              . . . . */
 439           static const struct read_matrix_params rmp = { 0, -1, 1, INT_MAX, 1, 0 };
 440           return &rmp;
 441         }
 442     }
 443   else
 444     NOT_REACHED ();
 445 }
 446
 447 static void
 448 schedule_matrices (struct matrix_format *mf)
 449 {
 450   struct matrix_sched *ms0 = &mf->ms[0];
 451   ms0->nr = 1;
 452   ms0->nc = 1;
 453   ms0->rp = xmalloc (sizeof *ms0->rp);
 454   ms0->rp[0] = (struct row_sched) { .y = 0, .x0 = 0, .x1 = 1 };
 455   ms0->n_rp = 1;
 456
 457   struct matrix_sched *ms1 = &mf->ms[1];
 458   ms1->nr = 1;
 459   ms1->nc = mf->n_cvars;
 460   ms1->rp = xmalloc (sizeof *ms1->rp);
 461   ms1->rp[0] = (struct row_sched) { .y = 0, .x0 = 0, .x1 = mf->n_cvars };
 462   ms1->n_rp = 1;
 463
 464   struct matrix_sched *ms2 = &mf->ms[2];
 465   ms2->nr = mf->n_cvars;
 466   ms2->nc = mf->n_cvars;
 467   ms2->rp = xmalloc (mf->n_cvars * sizeof *ms2->rp);
 468   ms2->n_rp = 0;
 469
 470   const struct read_matrix_params *rmp = get_read_matrix_params (mf);
 471   int x0 = rmp->x0;
 472   int x1 = rmp->x1 < mf->n_cvars ? rmp->x1 : mf->n_cvars - 1;
 473   int y0 = rmp->dy0;
 474   int y1 = (int) mf->n_cvars + rmp->dy1;
 475   for (int y = y0; y < y1; y++)
 476     {
 477       assert (x0 >= 0 && x0 < mf->n_cvars);
 478       assert (x1 >= 0 && x1 < mf->n_cvars);
 479       assert (x1 >= x0);
 480
 481       ms2->rp[ms2->n_rp++] = (struct row_sched) {
 482         .y = y, .x0 = x0, .x1 = x1 + 1
 483       };
 484
 485       x0 += rmp->dx0;
 486       x1 += rmp->dx1;
 487     }
 488 }
 489
 490 static bool
 491 read_id_columns (const struct matrix_format *mf,
 492                  struct substring *p, struct dfm_reader *r,
 493                  double *d, enum rowtype *rt)
 494 {
 495   for (size_t i = 0; mf->input_vars[i] != mf->cvars[0]; i++)
 496     if (!(mf->input_vars[i] == mf->rowtype
 497           ? next_rowtype (p, r, rt)
 498           : next_number (p, r, &d[i])))
 499       return false;
 500   return true;
 501 }
 502
 503 static bool
 504 equal_id_columns (const struct matrix_format *mf,
 505                   const double *a, const double *b)
 506 {
 507   for (size_t i = 0; mf->input_vars[i] != mf->cvars[0]; i++)
 508     if (mf->input_vars[i] != mf->rowtype && a[i] != b[i])
 509       return false;
 510   return true;
 511 }
 512
 513 static bool
 514 equal_split_columns (const struct matrix_format *mf,
 515                      const double *a, const double *b)
 516 {
 517   for (size_t i = 0; i < mf->n_svars; i++)
 518     {
 519       size_t idx = mf->svar_indexes[i];
 520       if (a[idx] != b[idx])
 521         return false;
 522     }
 523   return true;
 524 }
 525
 526 static bool
 527 is_pooled (const struct matrix_format *mf, const double *d)
 528 {
 529   for (size_t i = 0; i < mf->n_fvars; i++)
 530     if (d[mf->fvar_indexes[i]] != SYSMIS)
 531       return false;
 532   return true;
 533 }
 534
 535 static void
 536 matrix_sched_init (const struct matrix_format *mf, enum rowtype rt,
 537                    gsl_matrix *m)
 538 {
 539   int n_dims = rowtype_dimensions (rt);
 540   const struct matrix_sched *ms = &mf->ms[n_dims];
 541   double diagonal = n_dims < 2 || rt != C_CORR ? SYSMIS : 1.0;
 542   for (size_t y = 0; y < ms->nr; y++)
 543     for (size_t x = 0; x < ms->nc; x++)
 544       gsl_matrix_set (m, y, x, y == x ? diagonal : SYSMIS);
 545 }
 546
 547 static void
 548 matrix_sched_output (const struct matrix_format *mf, enum rowtype rt,
 549                      gsl_matrix *m, const double *d, int split_num,
 550                      struct casewriter *w)
 551 {
 552   int n_dims = rowtype_dimensions (rt);
 553   const struct matrix_sched *ms = &mf->ms[n_dims];
 554
 555   if (rt == C_N_SCALAR)
 556     {
 557       for (size_t x = 1; x < mf->n_cvars; x++)
 558         gsl_matrix_set (m, 0, x, gsl_matrix_get (m, 0, 0));
 559       rt = C_N;
 560     }
 561
 562   for (int y = 0; y < ms->nr; y++)
 563     {
 564       struct ccase *c = case_create (casewriter_get_proto (w));
 565       for (size_t i = 0; mf->input_vars[i] != mf->cvars[0]; i++)
 566         if (mf->input_vars[i] != mf->rowtype)
 567           *case_num_rw (c, mf->input_vars[i]) = d[i];
 568       if (mf->n_svars && !mf->svar_indexes)
 569         *case_num_rw (c, mf->svars[0]) = split_num;
 570       set_string (c, mf->rowtype, rowtype_name (rt));
 571       const char *varname = n_dims == 2 ? var_get_name (mf->cvars[y]) : "";
 572       set_string (c, mf->varname, ss_cstr (varname));
 573       for (int x = 0; x < mf->n_cvars; x++)
 574         *case_num_rw (c, mf->cvars[x]) = gsl_matrix_get (m, y, x);
 575       casewriter_write (w, c);
 576     }
 577 }
 578
 579 static void
 580 matrix_sched_output_n (const struct matrix_format *mf, double n,
 581                        gsl_matrix *m, const double *d, int split_num,
 582                        struct casewriter *w)
 583 {
 584   gsl_matrix_set (m, 0, 0, n);
 585   matrix_sched_output (mf, C_N_SCALAR, m, d, split_num, w);
 586 }
 587
 588 static void
 589 check_eol (const struct matrix_format *mf, struct substring *p,
 590            struct dfm_reader *r)
 591 {
 592   if (!mf->span)
 593     {
 594       ss_ltrim (p, ss_cstr (CC_SPACES ","));
 595       if (p->length)
 596         {
 597           parse_error (r, p, _("Extraneous data expecting end of line."));
 598           p->length = 0;
 599         }
 600     }
 601 }
 602
 603 static void
 604 parse_data_with_rowtype (const struct matrix_format *mf,
 605                          struct dfm_reader *r, struct casewriter *w)
 606 {
 607   if (dfm_eof (r))
 608     return;
 609   struct substring p = dfm_get_record (r);
 610
 611   double *prev = NULL;
 612   gsl_matrix *m = gsl_matrix_alloc (mf->n_cvars, mf->n_cvars);
 613
 614   double *d = xnmalloc (mf->n_input_vars, sizeof *d);
 615   enum rowtype rt;
 616
 617   double *d_next = xnmalloc (mf->n_input_vars, sizeof *d_next);
 618
 619   if (!read_id_columns (mf, &p, r, d, &rt))
 620     goto exit;
 621   for (;;)
 622     {
 623       /* If this has rowtype N but there was an N subcommand, then the
 624          subcommand takes precedence, so we will suppress outputting this
 625          record.  We still need to parse it, though, so we can't skip other
 626          work. */
 627       bool suppress_output = mf->n >= 0 && (rt == C_N || rt == C_N_SCALAR);
 628       if (suppress_output)
 629         parse_error (r, NULL, _("N record is not allowed with N subcommand.  "
 630                                 "Ignoring N record."));
 631
 632       /* If there's an N subcommand, and this is a new split, then output an N
 633          record. */
 634       if (mf->n >= 0 && (!prev || !equal_split_columns (mf, prev, d)))
 635         {
 636           matrix_sched_output_n (mf, mf->n, m, d, 0, w);
 637
 638           if (!prev)
 639             prev = xnmalloc (mf->n_input_vars, sizeof *prev);
 640           memcpy (prev, d, mf->n_input_vars * sizeof *prev);
 641         }
 642
 643       /* Usually users don't provide the CONTENTS subcommand with ROWTYPE_, but
 644          if they did then warn if ROWTYPE_ is an unexpected type. */
 645       if (mf->factor_rowtype_mask || mf->pooled_rowtype_mask)
 646         {
 647           const char *name = rowtype_name (rt).string;
 648           if (is_pooled (mf, d))
 649             {
 650               if (!((1u << rt) & mf->pooled_rowtype_mask))
 651                 parse_warning (r, NULL, _("Data contains pooled row type %s not "
 652                                           "included in CONTENTS."), name);
 653             }
 654           else
 655             {
 656               if (!((1u << rt) & mf->factor_rowtype_mask))
 657                 parse_warning (r, NULL, _("Data contains with-factors row type "
 658                                           "%s not included in CONTENTS."), name);
 659             }
 660         }
 661
 662       /* Initialize the matrix to be filled-in. */
 663       int n_dims = rowtype_dimensions (rt);
 664       const struct matrix_sched *ms = &mf->ms[n_dims];
 665       matrix_sched_init (mf, rt, m);
 666
 667       enum rowtype rt_next;
 668       bool eof;
 669
 670       size_t n_rows;
 671       for (n_rows = 1; ; n_rows++)
 672         {
 673           if (n_rows <= ms->n_rp)
 674             {
 675               const struct row_sched *rs = &ms->rp[n_rows - 1];
 676               size_t y = rs->y;
 677               for (size_t x = rs->x0; x < rs->x1; x++)
 678                 {
 679                   double e;
 680                   if (!next_number (&p, r, &e))
 681                     goto exit;
 682                   gsl_matrix_set (m, y, x, e);
 683                   if (n_dims == 2 && mf->triangle != FULL)
 684                     gsl_matrix_set (m, x, y, e);
 685                 }
 686               check_eol (mf, &p, r);
 687             }
 688           else
 689             {
 690               /* Suppress bad input data.  We'll issue an error later. */
 691               p.length = 0;
 692             }
 693
 694           eof = (!more_tokens (&p, r)
 695                  || !read_id_columns (mf, &p, r, d_next, &rt_next));
 696           if (eof)
 697             break;
 698
 699           if (!equal_id_columns (mf, d, d_next) || rt_next != rt)
 700             break;
 701         }
 702       if (!suppress_output)
 703         matrix_sched_output (mf, rt, m, d, 0, w);
 704
 705       if (n_rows != ms->n_rp)
 706         parse_error (r, NULL,
 707                      _("Matrix %s had %zu rows but %zu rows were expected."),
 708                      rowtype_name (rt).string, n_rows, ms->n_rp);
 709       if (eof)
 710         break;
 711
 712       double *d_tmp = d;
 713       d = d_next;
 714       d_next = d_tmp;
 715
 716       rt = rt_next;
 717     }
 718
 719 exit:
 720   free (prev);
 721   gsl_matrix_free (m);
 722   free (d);
 723   free (d_next);
 724 }
 725
 726 static void
 727 parse_matrix_without_rowtype (const struct matrix_format *mf,
 728                               struct substring *p, struct dfm_reader *r,
 729                               gsl_matrix *m, enum rowtype rowtype, bool pooled,
 730                               int split_num, struct casewriter *w)
 731 {
 732   int n_dims = rowtype_dimensions (rowtype);
 733   const struct matrix_sched *ms = &mf->ms[n_dims];
 734
 735   double *d = xnmalloc (mf->n_input_vars, sizeof *d);
 736   matrix_sched_init (mf, rowtype, m);
 737   for (size_t i = 0; i < ms->n_rp; i++)
 738     {
 739       int y = ms->rp[i].y;
 740       int k = 0;
 741       int h = 0;
 742       for (size_t j = 0; j < mf->n_input_vars; j++)
 743         {
 744           const struct variable *iv = mf->input_vars[j];
 745           if (k < mf->n_cvars && iv == mf->cvars[k])
 746             {
 747               if (k < ms->rp[i].x1 - ms->rp[i].x0)
 748                 {
 749                   double e;
 750                   if (!next_number (p, r, &e))
 751                     goto exit;
 752
 753                   int x = k + ms->rp[i].x0;
 754                   gsl_matrix_set (m, y, x, e);
 755                   if (n_dims == 2 && mf->triangle != FULL)
 756                     gsl_matrix_set (m, x, y, e);
 757                 }
 758               k++;
 759               continue;
 760             }
 761           if (h < mf->n_fvars && iv == mf->fvars[h])
 762             {
 763               h++;
 764               if (pooled)
 765                 {
 766                   d[j] = SYSMIS;
 767                   continue;
 768                 }
 769             }
 770
 771           double e;
 772           if (!next_number (p, r, &e))
 773             goto exit;
 774           d[j] = e;
 775         }
 776       check_eol (mf, p, r);
 777     }
 778
 779   matrix_sched_output (mf, rowtype, m, d, split_num, w);
 780 exit:
 781   free (d);
 782 }
 783
 784 static void
 785 parse_data_without_rowtype (const struct matrix_format *mf,
 786                             struct dfm_reader *r, struct casewriter *w)
 787 {
 788   if (dfm_eof (r))
 789     return;
 790   struct substring p = dfm_get_record (r);
 791
 792   gsl_matrix *m = gsl_matrix_alloc (mf->n_cvars, mf->n_cvars);
 793
 794   int split_num = 1;
 795   do
 796     {
 797       for (size_t i = 0; i < mf->n_contents; )
 798         {
 799           size_t j = i;
 800           if (mf->contents[i].open)
 801             while (!mf->contents[j].close)
 802               j++;
 803
 804           if (mf->contents[i].open)
 805             {
 806               for (size_t k = 0; k < mf->cells; k++)
 807                 for (size_t h = i; h <= j; h++)
 808                   parse_matrix_without_rowtype (mf, &p, r, m,
 809                                                 mf->contents[h].rowtype, false,
 810                                                 split_num, w);
 811             }
 812           else
 813             parse_matrix_without_rowtype (mf, &p, r, m, mf->contents[i].rowtype,
 814                                           true, split_num, w);
 815           i = j + 1;
 816         }
 817
 818       split_num++;
 819     }
 820   while (more_tokens (&p, r));
 821
 822   gsl_matrix_free (m);
 823 }
 824
 825 /* Parses VARIABLES=varnames for MATRIX DATA and returns a dictionary with the
 826    named variables in it. */
 827 static struct dictionary *
 828 parse_matrix_data_variables (struct lexer *lexer)
 829 {
 830   if (!lex_force_match_id (lexer, "VARIABLES"))
 831     return NULL;
 832   lex_match (lexer, T_EQUALS);
 833
 834   struct dictionary *dict = dict_create (get_default_encoding ());
 835
 836   size_t n_names = 0;
 837   char **names = NULL;
 838   if (!parse_DATA_LIST_vars (lexer, dict, &names, &n_names, PV_NO_DUPLICATE))
 839     {
 840       dict_unref (dict);
 841       return NULL;
 842     }
 843
 844   for (size_t i = 0; i < n_names; i++)
 845     if (!strcasecmp (names[i], "ROWTYPE_"))
 846       dict_create_var_assert (dict, "ROWTYPE_", 8);
 847     else
 848       dict_create_var_assert (dict, names[i], 0);
 849
 850   for (size_t i = 0; i < n_names; ++i)
 851     free (names[i]);
 852   free (names);
 853
 854   if (dict_lookup_var (dict, "VARNAME_"))
 855     {
 856       msg (SE, _("VARIABLES may not include VARNAME_."));
 857       dict_unref (dict);
 858       return NULL;
 859     }
 860   return dict;
 861 }
 862
 863 static bool
 864 parse_matrix_data_subvars (struct lexer *lexer, struct dictionary *dict,
 865                            bool *taken_vars,
 866                            struct variable ***vars, size_t **indexes,
 867                            size_t *n_vars)
 868 {
 869   if (!parse_variables (lexer, dict, vars, n_vars, 0))
 870     return false;
 871
 872   *indexes = xnmalloc (*n_vars, sizeof **indexes);
 873   for (size_t i = 0; i < *n_vars; i++)
 874     {
 875       struct variable *v = (*vars)[i];
 876       if (!strcasecmp (var_get_name (v), "ROWTYPE_"))
 877         {
 878           msg (SE, _("ROWTYPE_ is not allowed on SPLIT or FACTORS."));
 879           goto error;
 880         }
 881       (*indexes)[i] = var_get_dict_index (v);
 882
 883       bool *tv = &taken_vars[var_get_dict_index (v)];
 884       if (*tv)
 885         {
 886           msg (SE, _("%s may not appear on both SPLIT and FACTORS."),
 887                var_get_name (v));
 888           goto error;
 889         }
 890       *tv = true;
 891
 892       var_set_both_formats (v, &(struct fmt_spec) { .type = FMT_F, .w = 4 });
 893     }
 894   return true;
 895
 896 error:
 897   free (*vars);
 898   *vars = NULL;
 899   *n_vars = 0;
 900   free (*indexes);
 901   *indexes = NULL;
 902   return false;
 903 }
 904
 905 int
 906 cmd_matrix_data (struct lexer *lexer, struct dataset *ds)
 907 {
 908   struct dictionary *dict = parse_matrix_data_variables (lexer);
 909   if (!dict)
 910     return CMD_FAILURE;
 911
 912   size_t n_input_vars = dict_get_var_cnt (dict);
 913   struct variable **input_vars = xnmalloc (n_input_vars, sizeof *input_vars);
 914   for (size_t i = 0; i < n_input_vars; i++)
 915     input_vars[i] = dict_get_var (dict, i);
 916
 917   int varname_width = 8;
 918   for (size_t i = 0; i < n_input_vars; i++)
 919     {
 920       int w = strlen (var_get_name (input_vars[i]));
 921       varname_width = MAX (w, varname_width);
 922     }
 923
 924   struct variable *rowtype = dict_lookup_var (dict, "ROWTYPE_");
 925   bool input_rowtype = rowtype != NULL;
 926   if (!rowtype)
 927     rowtype = dict_create_var_assert (dict, "ROWTYPE_", 8);
 928
 929   struct matrix_format mf = {
 930     .input_rowtype = input_rowtype,
 931     .input_vars = input_vars,
 932     .n_input_vars = n_input_vars,
 933
 934     .rowtype = rowtype,
 935     .varname = dict_create_var_assert (dict, "VARNAME_", varname_width),
 936
 937     .triangle = LOWER,
 938     .diagonal = DIAGONAL,
 939     .n = -1,
 940     .cells = -1,
 941   };
 942
 943   bool *taken_vars = XCALLOC (n_input_vars, bool);
 944   if (input_rowtype)
 945     taken_vars[var_get_dict_index (rowtype)] = true;
 946
 947   struct file_handle *fh = NULL;
 948   while (lex_token (lexer) != T_ENDCMD)
 949     {
 950       if (!lex_force_match (lexer, T_SLASH))
 951         goto error;
 952
 953       if (lex_match_id (lexer, "N"))
 954         {
 955           lex_match (lexer, T_EQUALS);
 956
 957           if (!lex_force_int_range (lexer, "N", 0, INT_MAX))
 958             goto error;
 959
 960           mf.n = lex_integer (lexer);
 961           lex_get (lexer);
 962         }
 963       else if (lex_match_id (lexer, "FORMAT"))
 964         {
 965           lex_match (lexer, T_EQUALS);
 966
 967           while (lex_token (lexer) != T_SLASH && lex_token (lexer) != T_ENDCMD)
 968             {
 969               if (lex_match_id (lexer, "LIST"))
 970                 mf.span = false;
 971               else if (lex_match_id (lexer, "FREE"))
 972                 mf.span = true;
 973               else if (lex_match_id (lexer, "UPPER"))
 974                 mf.triangle = UPPER;
 975               else if (lex_match_id (lexer, "LOWER"))
 976                 mf.triangle = LOWER;
 977               else if (lex_match_id (lexer, "FULL"))
 978                 mf.triangle = FULL;
 979               else if (lex_match_id (lexer, "DIAGONAL"))
 980                 mf.diagonal = DIAGONAL;
 981               else if (lex_match_id (lexer, "NODIAGONAL"))
 982                 mf.diagonal = NO_DIAGONAL;
 983               else
 984                 {
 985                   lex_error (lexer, NULL);
 986                   goto error;
 987                 }
 988             }
 989         }
 990       else if (lex_match_id (lexer, "FILE"))
 991         {
 992           lex_match (lexer, T_EQUALS);
 993           fh_unref (fh);
 994           fh = fh_parse (lexer, FH_REF_FILE | FH_REF_INLINE, NULL);
 995           if (!fh)
 996             goto error;
 997         }
 998       else if (!mf.n_svars && lex_match_id (lexer, "SPLIT"))
 999         {
1000           lex_match (lexer, T_EQUALS);
1001           if (!mf.input_rowtype
1002               && lex_token (lexer) == T_ID
1003               && !dict_lookup_var (dict, lex_tokcstr (lexer)))
1004             {
1005               mf.svars = xmalloc (sizeof *mf.svars);
1006               mf.svars[0] = dict_create_var_assert (dict, lex_tokcstr (lexer),
1007                                                     0);
1008               var_set_both_formats (
1009                 mf.svars[0], &(struct fmt_spec) { .type = FMT_F, .w = 4 });
1010               mf.n_svars = 1;
1011               lex_get (lexer);
1012             }
1013           else if (!parse_matrix_data_subvars (lexer, dict, taken_vars,
1014                                                &mf.svars, &mf.svar_indexes,
1015                                                &mf.n_svars))
1016             goto error;
1017         }
1018       else if (!mf.n_fvars && lex_match_id (lexer, "FACTORS"))
1019         {
1020           lex_match (lexer, T_EQUALS);
1021           if (!parse_matrix_data_subvars (lexer, dict, taken_vars,
1022                                           &mf.fvars, &mf.fvar_indexes,
1023                                           &mf.n_fvars))
1024             goto error;
1025         }
1026       else if (lex_match_id (lexer, "CELLS"))
1027         {
1028           if (mf.input_rowtype)
1029             msg (SW, _("CELLS is ignored when VARIABLES includes ROWTYPE_"));
1030
1031           lex_match (lexer, T_EQUALS);
1032
1033           if (!lex_force_int_range (lexer, "CELLS", 0, INT_MAX))
1034             goto error;
1035
1036           mf.cells = lex_integer (lexer);
1037           lex_get (lexer);
1038         }
1039       else if (lex_match_id (lexer, "CONTENTS"))
1040         {
1041           lex_match (lexer, T_EQUALS);
1042
1043           size_t allocated_contents = mf.n_contents;
1044           bool in_parens = false;
1045           for (;;)
1046             {
1047               bool open = !in_parens && lex_match (lexer, T_LPAREN);
1048               enum rowtype rt;
1049               if (!rowtype_parse (lexer, &rt))
1050                 {
1051                   if (open || in_parens || (lex_token (lexer) != T_ENDCMD
1052                                             && lex_token (lexer) != T_SLASH))
1053                     {
1054                       lex_error (lexer, _("Row type keyword expected."));
1055                       goto error;
1056                     }
1057                   break;
1058                 }
1059
1060               if (open)
1061                 in_parens = true;
1062
1063               if (in_parens)
1064                 mf.factor_rowtype_mask |= 1u << rt;
1065               else
1066                 mf.pooled_rowtype_mask |= 1u << rt;
1067
1068               bool close = in_parens && lex_match (lexer, T_RPAREN);
1069               if (close)
1070                 in_parens = false;
1071
1072               if (mf.n_contents >= allocated_contents)
1073                 mf.contents = x2nrealloc (mf.contents, &allocated_contents,
1074                                           sizeof *mf.contents);
1075               mf.contents[mf.n_contents++] = (struct content) {
1076                 .open = open, .rowtype = rt, .close = close
1077               };
1078             }
1079         }
1080       else
1081         {
1082           lex_error (lexer, NULL);
1083           goto error;
1084         }
1085     }
1086   if (mf.diagonal == NO_DIAGONAL && mf.triangle == FULL)
1087     {
1088       msg (SE, _("FORMAT=FULL and FORMAT=NODIAGONAL are mutually exclusive."));
1089       goto error;
1090     }
1091   if (!mf.input_rowtype)
1092     {
1093       if (mf.cells < 0)
1094         {
1095           if (mf.n_fvars)
1096             {
1097               msg (SE, _("CELLS is required when factor variables are specified "
1098                          "and VARIABLES does not include ROWTYPE_."));
1099               goto error;
1100             }
1101           mf.cells = 1;
1102         }
1103
1104       if (!mf.n_contents)
1105         {
1106           msg (SW, _("CONTENTS was not specified and VARIABLES does not "
1107                      "include ROWTYPE_.  Assuming CONTENTS=CORR."));
1108
1109           mf.n_contents = 1;
1110           mf.contents = xmalloc (sizeof *mf.contents);
1111           *mf.contents = (struct content) { .rowtype = C_CORR };
1112         }
1113     }
1114   mf.cvars = xmalloc (mf.n_input_vars * sizeof *mf.cvars);
1115   for (size_t i = 0; i < mf.n_input_vars; i++)
1116     if (!taken_vars[i])
1117       {
1118         struct variable *v = input_vars[i];
1119         mf.cvars[mf.n_cvars++] = v;
1120         var_set_both_formats (v, &(struct fmt_spec) { .type = FMT_F, .w = 10,
1121                                                       .d = 4 });
1122       }
1123   if (!mf.n_cvars)
1124     {
1125       msg (SE, _("At least one continuous variable is required."));
1126       goto error;
1127     }
1128   if (mf.input_rowtype)
1129     {
1130       for (size_t i = 0; i < mf.n_cvars; i++)
1131         if (mf.cvars[i] != input_vars[n_input_vars - mf.n_cvars + i])
1132           {
1133             msg (SE, _("VARIABLES includes ROWTYPE_ but the continuous "
1134                        "variables are not the last ones on VARIABLES."));
1135             goto error;
1136           }
1137     }
1138   unsigned int rowtype_mask = mf.pooled_rowtype_mask | mf.factor_rowtype_mask;
1139   if (rowtype_mask & (1u << C_N) && mf.n >= 0)
1140     {
1141       msg (SE, _("Cannot specify N on CONTENTS along with the N subcommand."));
1142       goto error;
1143     }
1144
1145   struct variable **order = xnmalloc (dict_get_var_cnt (dict), sizeof *order);
1146   size_t n_order = 0;
1147   for (size_t i = 0; i < mf.n_svars; i++)
1148     order[n_order++] = mf.svars[i];
1149   order[n_order++] = mf.rowtype;
1150   for (size_t i = 0; i < mf.n_fvars; i++)
1151     order[n_order++] = mf.fvars[i];
1152   order[n_order++] = mf.varname;
1153   for (size_t i = 0; i < mf.n_cvars; i++)
1154     order[n_order++] = mf.cvars[i];
1155   assert (n_order == dict_get_var_cnt (dict));
1156   dict_reorder_vars (dict, order, n_order);
1157   free (order);
1158
1159   dict_set_split_vars (dict, mf.svars, mf.n_svars);
1160
1161   schedule_matrices (&mf);
1162
1163   if (fh == NULL)
1164     fh = fh_inline_file ();
1165
1166   if (lex_end_of_command (lexer) != CMD_SUCCESS)
1167     goto error;
1168
1169   struct dfm_reader *reader = dfm_open_reader (fh, lexer, NULL);
1170   if (reader == NULL)
1171     goto error;
1172
1173   struct casewriter *writer = autopaging_writer_create (dict_get_proto (dict));
1174   if (mf.input_rowtype)
1175     parse_data_with_rowtype (&mf, reader, writer);
1176   else
1177     parse_data_without_rowtype (&mf, reader, writer);
1178   dfm_close_reader (reader);
1179
1180   dataset_set_dict (ds, dict);
1181   dataset_set_source (ds, casewriter_make_reader (writer));
1182
1183   matrix_format_uninit (&mf);
1184   free (taken_vars);
1185   fh_unref (fh);
1186
1187   return CMD_SUCCESS;
1188
1189  error:
1190   matrix_format_uninit (&mf);
1191   free (taken_vars);
1192   dict_unref (dict);
1193   fh_unref (fh);
1194   return CMD_FAILURE;
1195 }