refactoring
[pspp] / 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       {
851         struct variable *var = dict_create_var_assert (dict, names[i], 0);
852         var_set_measure (var, MEASURE_SCALE);
853       }
854
855   for (size_t i = 0; i < n_names; ++i)
856     free (names[i]);
857   free (names);
858
859   if (dict_lookup_var (dict, "VARNAME_"))
860     {
861       msg (SE, _("VARIABLES may not include VARNAME_."));
862       dict_unref (dict);
863       return NULL;
864     }
865   return dict;
866 }
867
868 static bool
869 parse_matrix_data_subvars (struct lexer *lexer, struct dictionary *dict,
870                            bool *taken_vars,
871                            struct variable ***vars, size_t **indexes,
872                            size_t *n_vars)
873 {
874   if (!parse_variables (lexer, dict, vars, n_vars, 0))
875     return false;
876
877   *indexes = xnmalloc (*n_vars, sizeof **indexes);
878   for (size_t i = 0; i < *n_vars; i++)
879     {
880       struct variable *v = (*vars)[i];
881       if (!strcasecmp (var_get_name (v), "ROWTYPE_"))
882         {
883           msg (SE, _("ROWTYPE_ is not allowed on SPLIT or FACTORS."));
884           goto error;
885         }
886       (*indexes)[i] = var_get_dict_index (v);
887
888       bool *tv = &taken_vars[var_get_dict_index (v)];
889       if (*tv)
890         {
891           msg (SE, _("%s may not appear on both SPLIT and FACTORS."),
892                var_get_name (v));
893           goto error;
894         }
895       *tv = true;
896
897       var_set_measure (v, MEASURE_NOMINAL);
898       var_set_both_formats (v, &(struct fmt_spec) { .type = FMT_F, .w = 4 });
899     }
900   return true;
901
902 error:
903   free (*vars);
904   *vars = NULL;
905   *n_vars = 0;
906   free (*indexes);
907   *indexes = NULL;
908   return false;
909 }
910
911 int
912 cmd_matrix_data (struct lexer *lexer, struct dataset *ds)
913 {
914   struct dictionary *dict = parse_matrix_data_variables (lexer);
915   if (!dict)
916     return CMD_FAILURE;
917
918   size_t n_input_vars = dict_get_n_vars (dict);
919   struct variable **input_vars = xnmalloc (n_input_vars, sizeof *input_vars);
920   for (size_t i = 0; i < n_input_vars; i++)
921     input_vars[i] = dict_get_var (dict, i);
922
923   int varname_width = 8;
924   for (size_t i = 0; i < n_input_vars; i++)
925     {
926       int w = strlen (var_get_name (input_vars[i]));
927       varname_width = MAX (w, varname_width);
928     }
929
930   struct variable *rowtype = dict_lookup_var (dict, "ROWTYPE_");
931   bool input_rowtype = rowtype != NULL;
932   if (!rowtype)
933     rowtype = dict_create_var_assert (dict, "ROWTYPE_", 8);
934
935   struct matrix_format mf = {
936     .input_rowtype = input_rowtype,
937     .input_vars = input_vars,
938     .n_input_vars = n_input_vars,
939
940     .rowtype = rowtype,
941     .varname = dict_create_var_assert (dict, "VARNAME_", varname_width),
942
943     .triangle = LOWER,
944     .diagonal = DIAGONAL,
945     .n = -1,
946     .cells = -1,
947   };
948
949   bool *taken_vars = XCALLOC (n_input_vars, bool);
950   if (input_rowtype)
951     taken_vars[var_get_dict_index (rowtype)] = true;
952
953   struct file_handle *fh = NULL;
954   while (lex_token (lexer) != T_ENDCMD)
955     {
956       if (!lex_force_match (lexer, T_SLASH))
957         goto error;
958
959       if (lex_match_id (lexer, "N"))
960         {
961           lex_match (lexer, T_EQUALS);
962
963           if (!lex_force_int_range (lexer, "N", 0, INT_MAX))
964             goto error;
965
966           mf.n = lex_integer (lexer);
967           lex_get (lexer);
968         }
969       else if (lex_match_id (lexer, "FORMAT"))
970         {
971           lex_match (lexer, T_EQUALS);
972
973           while (lex_token (lexer) != T_SLASH && lex_token (lexer) != T_ENDCMD)
974             {
975               if (lex_match_id (lexer, "LIST"))
976                 mf.span = false;
977               else if (lex_match_id (lexer, "FREE"))
978                 mf.span = true;
979               else if (lex_match_id (lexer, "UPPER"))
980                 mf.triangle = UPPER;
981               else if (lex_match_id (lexer, "LOWER"))
982                 mf.triangle = LOWER;
983               else if (lex_match_id (lexer, "FULL"))
984                 mf.triangle = FULL;
985               else if (lex_match_id (lexer, "DIAGONAL"))
986                 mf.diagonal = DIAGONAL;
987               else if (lex_match_id (lexer, "NODIAGONAL"))
988                 mf.diagonal = NO_DIAGONAL;
989               else
990                 {
991                   lex_error (lexer, NULL);
992                   goto error;
993                 }
994             }
995         }
996       else if (lex_match_id (lexer, "FILE"))
997         {
998           lex_match (lexer, T_EQUALS);
999           fh_unref (fh);
1000           fh = fh_parse (lexer, FH_REF_FILE | FH_REF_INLINE, NULL);
1001           if (!fh)
1002             goto error;
1003         }
1004       else if (!mf.n_svars && lex_match_id (lexer, "SPLIT"))
1005         {
1006           lex_match (lexer, T_EQUALS);
1007           if (!mf.input_rowtype
1008               && lex_token (lexer) == T_ID
1009               && !dict_lookup_var (dict, lex_tokcstr (lexer)))
1010             {
1011               mf.svars = xmalloc (sizeof *mf.svars);
1012               mf.svars[0] = dict_create_var_assert (dict, lex_tokcstr (lexer),
1013                                                     0);
1014               var_set_measure (mf.svars[0], MEASURE_NOMINAL);
1015               var_set_both_formats (
1016                 mf.svars[0], &(struct fmt_spec) { .type = FMT_F, .w = 4 });
1017               mf.n_svars = 1;
1018               lex_get (lexer);
1019             }
1020           else if (!parse_matrix_data_subvars (lexer, dict, taken_vars,
1021                                                &mf.svars, &mf.svar_indexes,
1022                                                &mf.n_svars))
1023             goto error;
1024         }
1025       else if (!mf.n_fvars && lex_match_id (lexer, "FACTORS"))
1026         {
1027           lex_match (lexer, T_EQUALS);
1028           if (!parse_matrix_data_subvars (lexer, dict, taken_vars,
1029                                           &mf.fvars, &mf.fvar_indexes,
1030                                           &mf.n_fvars))
1031             goto error;
1032         }
1033       else if (lex_match_id (lexer, "CELLS"))
1034         {
1035           if (mf.input_rowtype)
1036             msg (SW, _("CELLS is ignored when VARIABLES includes ROWTYPE_"));
1037
1038           lex_match (lexer, T_EQUALS);
1039
1040           if (!lex_force_int_range (lexer, "CELLS", 0, INT_MAX))
1041             goto error;
1042
1043           mf.cells = lex_integer (lexer);
1044           lex_get (lexer);
1045         }
1046       else if (lex_match_id (lexer, "CONTENTS"))
1047         {
1048           lex_match (lexer, T_EQUALS);
1049
1050           size_t allocated_contents = mf.n_contents;
1051           bool in_parens = false;
1052           for (;;)
1053             {
1054               bool open = !in_parens && lex_match (lexer, T_LPAREN);
1055               enum rowtype rt;
1056               if (!rowtype_parse (lexer, &rt))
1057                 {
1058                   if (open || in_parens || (lex_token (lexer) != T_ENDCMD
1059                                             && lex_token (lexer) != T_SLASH))
1060                     {
1061                       lex_error (lexer, _("Row type keyword expected."));
1062                       goto error;
1063                     }
1064                   break;
1065                 }
1066
1067               if (open)
1068                 in_parens = true;
1069
1070               if (in_parens)
1071                 mf.factor_rowtype_mask |= 1u << rt;
1072               else
1073                 mf.pooled_rowtype_mask |= 1u << rt;
1074
1075               bool close = in_parens && lex_match (lexer, T_RPAREN);
1076               if (close)
1077                 in_parens = false;
1078
1079               if (mf.n_contents >= allocated_contents)
1080                 mf.contents = x2nrealloc (mf.contents, &allocated_contents,
1081                                           sizeof *mf.contents);
1082               mf.contents[mf.n_contents++] = (struct content) {
1083                 .open = open, .rowtype = rt, .close = close
1084               };
1085             }
1086         }
1087       else
1088         {
1089           lex_error (lexer, NULL);
1090           goto error;
1091         }
1092     }
1093   if (mf.diagonal == NO_DIAGONAL && mf.triangle == FULL)
1094     {
1095       msg (SE, _("FORMAT=FULL and FORMAT=NODIAGONAL are mutually exclusive."));
1096       goto error;
1097     }
1098   if (!mf.input_rowtype)
1099     {
1100       if (mf.cells < 0)
1101         {
1102           if (mf.n_fvars)
1103             {
1104               msg (SE, _("CELLS is required when factor variables are specified "
1105                          "and VARIABLES does not include ROWTYPE_."));
1106               goto error;
1107             }
1108           mf.cells = 1;
1109         }
1110
1111       if (!mf.n_contents)
1112         {
1113           msg (SW, _("CONTENTS was not specified and VARIABLES does not "
1114                      "include ROWTYPE_.  Assuming CONTENTS=CORR."));
1115
1116           mf.n_contents = 1;
1117           mf.contents = xmalloc (sizeof *mf.contents);
1118           *mf.contents = (struct content) { .rowtype = C_CORR };
1119         }
1120     }
1121   mf.cvars = xmalloc (mf.n_input_vars * sizeof *mf.cvars);
1122   for (size_t i = 0; i < mf.n_input_vars; i++)
1123     if (!taken_vars[i])
1124       {
1125         struct variable *v = input_vars[i];
1126         mf.cvars[mf.n_cvars++] = v;
1127         var_set_both_formats (v, &(struct fmt_spec) { .type = FMT_F, .w = 10,
1128                                                       .d = 4 });
1129       }
1130   if (!mf.n_cvars)
1131     {
1132       msg (SE, _("At least one continuous variable is required."));
1133       goto error;
1134     }
1135   if (mf.input_rowtype)
1136     {
1137       for (size_t i = 0; i < mf.n_cvars; i++)
1138         if (mf.cvars[i] != input_vars[n_input_vars - mf.n_cvars + i])
1139           {
1140             msg (SE, _("VARIABLES includes ROWTYPE_ but the continuous "
1141                        "variables are not the last ones on VARIABLES."));
1142             goto error;
1143           }
1144     }
1145   unsigned int rowtype_mask = mf.pooled_rowtype_mask | mf.factor_rowtype_mask;
1146   if (rowtype_mask & (1u << C_N) && mf.n >= 0)
1147     {
1148       msg (SE, _("Cannot specify N on CONTENTS along with the N subcommand."));
1149       goto error;
1150     }
1151
1152   struct variable **order = xnmalloc (dict_get_n_vars (dict), sizeof *order);
1153   size_t n_order = 0;
1154   for (size_t i = 0; i < mf.n_svars; i++)
1155     order[n_order++] = mf.svars[i];
1156   order[n_order++] = mf.rowtype;
1157   for (size_t i = 0; i < mf.n_fvars; i++)
1158     order[n_order++] = mf.fvars[i];
1159   order[n_order++] = mf.varname;
1160   for (size_t i = 0; i < mf.n_cvars; i++)
1161     order[n_order++] = mf.cvars[i];
1162   assert (n_order == dict_get_n_vars (dict));
1163   dict_reorder_vars (dict, order, n_order);
1164   free (order);
1165
1166   dict_set_split_vars (dict, mf.svars, mf.n_svars, SPLIT_LAYERED);
1167
1168   schedule_matrices (&mf);
1169
1170   if (fh == NULL)
1171     fh = fh_inline_file ();
1172
1173   if (lex_end_of_command (lexer) != CMD_SUCCESS)
1174     goto error;
1175
1176   struct dfm_reader *reader = dfm_open_reader (fh, lexer, NULL);
1177   if (reader == NULL)
1178     goto error;
1179
1180   struct casewriter *writer = autopaging_writer_create (dict_get_proto (dict));
1181   if (mf.input_rowtype)
1182     parse_data_with_rowtype (&mf, reader, writer);
1183   else
1184     parse_data_without_rowtype (&mf, reader, writer);
1185   dfm_close_reader (reader);
1186
1187   dataset_set_dict (ds, dict);
1188   dataset_set_source (ds, casewriter_make_reader (writer));
1189
1190   matrix_format_uninit (&mf);
1191   free (taken_vars);
1192   fh_unref (fh);
1193
1194   return CMD_SUCCESS;
1195
1196  error:
1197   matrix_format_uninit (&mf);
1198   free (taken_vars);
1199   dict_unref (dict);
1200   fh_unref (fh);
1201   return CMD_FAILURE;
1202 }