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