1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2021 Free Software Foundation, Inc.
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.
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.
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/>. */
19 #include <gsl/gsl_blas.h>
20 #include <gsl/gsl_eigen.h>
21 #include <gsl/gsl_linalg.h>
22 #include <gsl/gsl_matrix.h>
23 #include <gsl/gsl_permutation.h>
24 #include <gsl/gsl_randist.h>
25 #include <gsl/gsl_vector.h>
30 #include "data/any-reader.h"
31 #include "data/any-writer.h"
32 #include "data/casereader.h"
33 #include "data/casewriter.h"
34 #include "data/data-in.h"
35 #include "data/data-out.h"
36 #include "data/dataset.h"
37 #include "data/dictionary.h"
38 #include "data/file-handle-def.h"
39 #include "language/command.h"
40 #include "language/data-io/data-reader.h"
41 #include "language/data-io/data-writer.h"
42 #include "language/data-io/file-handle.h"
43 #include "language/lexer/format-parser.h"
44 #include "language/lexer/lexer.h"
45 #include "language/lexer/variable-parser.h"
46 #include "libpspp/array.h"
47 #include "libpspp/assertion.h"
48 #include "libpspp/compiler.h"
49 #include "libpspp/hmap.h"
50 #include "libpspp/i18n.h"
51 #include "libpspp/misc.h"
52 #include "libpspp/str.h"
53 #include "libpspp/string-array.h"
54 #include "libpspp/stringi-set.h"
55 #include "libpspp/u8-line.h"
56 #include "math/random.h"
57 #include "output/driver.h"
58 #include "output/output-item.h"
59 #include "output/pivot-table.h"
61 #include "gl/c-ctype.h"
62 #include "gl/ftoastr.h"
63 #include "gl/minmax.h"
67 #define _(msgid) gettext (msgid)
68 #define N_(msgid) (msgid)
72 struct hmap_node hmap_node;
80 struct file_handle *outfile;
81 struct string_array variables;
82 struct string_array fnames;
83 struct string_array snames;
86 /* Collects and owns factors and splits. The individual msave_command
87 structs point to these but do not own them. */
88 struct matrix_expr **factors;
89 size_t n_factors, allocated_factors;
90 struct matrix_expr **splits;
91 size_t n_splits, allocated_splits;
93 /* Execution state. */
94 struct dictionary *dict;
95 struct casewriter *writer;
100 struct file_handle *file;
101 struct dfm_reader *reader;
107 struct file_handle *file;
108 struct dfm_writer *writer;
110 struct u8_line *held;
115 struct file_handle *file;
116 struct dataset *dataset;
118 /* Parameters from parsing the first SAVE command for 'file'. */
119 struct string_array variables;
120 struct matrix_expr *names;
121 struct stringi_set strings;
123 /* Results from the first (only) attempt to open this save_file. */
125 struct casewriter *writer;
126 struct dictionary *dict;
131 struct dataset *dataset;
132 struct session *session;
136 struct msave_common *common;
138 struct file_handle *prev_read_file;
139 struct read_file **read_files;
142 struct file_handle *prev_write_file;
143 struct write_file **write_files;
144 size_t n_write_files;
146 struct file_handle *prev_save_file;
147 struct save_file **save_files;
151 static struct matrix_var *
152 matrix_var_lookup (struct matrix_state *s, struct substring name)
154 struct matrix_var *var;
156 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
157 utf8_hash_case_substring (name, 0), &s->vars)
158 if (!utf8_sscasecmp (ss_cstr (var->name), name))
164 static struct matrix_var *
165 matrix_var_create (struct matrix_state *s, struct substring name)
167 struct matrix_var *var = xmalloc (sizeof *var);
168 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
169 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
174 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
176 gsl_matrix_free (var->value);
180 #define MATRIX_FUNCTIONS \
184 F(ARSIN, "ARSIN", m_m) \
185 F(ARTAN, "ARTAN", m_m) \
186 F(BLOCK, "BLOCK", m_any) \
187 F(CHOL, "CHOL", m_m) \
188 F(CMIN, "CMIN", m_m) \
189 F(CMAX, "CMAX", m_m) \
191 F(CSSQ, "CSSQ", m_m) \
192 F(CSUM, "CSUM", m_m) \
193 F(DESIGN, "DESIGN", m_m) \
195 F(DIAG, "DIAG", m_m) \
196 F(EVAL, "EVAL", m_m) \
198 F(GINV, "GINV", m_m) \
199 F(GRADE, "GRADE", m_m) \
200 F(GSCH, "GSCH", m_m) \
201 F(IDENT, "IDENT", IDENT) \
203 F(KRONEKER, "KRONEKER", m_mm) \
204 F(LG10, "LG10", m_m) \
206 F(MAGIC, "MAGIC", m_d) \
207 F(MAKE, "MAKE", m_ddd) \
208 F(MDIAG, "MDIAG", m_v) \
209 F(MMAX, "MMAX", d_m) \
210 F(MMIN, "MMIN", d_m) \
211 F(MOD, "MOD", m_md) \
212 F(MSSQ, "MSSQ", d_m) \
213 F(MSUM, "MSUM", d_m) \
214 F(NCOL, "NCOL", d_m) \
215 F(NROW, "NROW", d_m) \
216 F(RANK, "RANK", d_m) \
217 F(RESHAPE, "RESHAPE", m_mdd) \
218 F(RMAX, "RMAX", m_m) \
219 F(RMIN, "RMIN", m_m) \
221 F(RNKORDER, "RNKORDER", m_m) \
222 F(RSSQ, "RSSQ", m_m) \
223 F(RSUM, "RSUM", m_m) \
225 F(SOLVE, "SOLVE", m_mm) \
226 F(SQRT, "SQRT", m_m) \
227 F(SSCP, "SSCP", m_m) \
228 F(SVAL, "SVAL", m_m) \
229 F(SWEEP, "SWEEP", m_md) \
231 F(TRACE, "TRACE", d_m) \
232 F(TRANSPOS, "TRANSPOS", m_m) \
233 F(TRUNC, "TRUNC", m_m) \
234 F(UNIFORM, "UNIFORM", m_dd)
236 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
237 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
238 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
239 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
240 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
241 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
242 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
243 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
244 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
245 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
246 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
248 typedef gsl_matrix *matrix_proto_m_d (double);
249 typedef gsl_matrix *matrix_proto_m_dd (double, double);
250 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
251 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
252 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
253 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
254 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
255 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
256 typedef double matrix_proto_d_m (gsl_matrix *);
257 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
258 typedef gsl_matrix *matrix_proto_IDENT (double, double);
260 #define F(ENUM, STRING, PROTO) static matrix_proto_##PROTO matrix_eval_##ENUM;
264 /* Matrix expressions. */
271 #define F(ENUM, STRING, PROTO) MOP_F_##ENUM,
275 /* Elementwise and scalar arithmetic. */
276 MOP_NEGATE, /* unary - */
277 MOP_ADD_ELEMS, /* + */
278 MOP_SUB_ELEMS, /* - */
279 MOP_MUL_ELEMS, /* &* */
280 MOP_DIV_ELEMS, /* / and &/ */
281 MOP_EXP_ELEMS, /* &** */
283 MOP_SEQ_BY, /* a:b:c */
285 /* Matrix arithmetic. */
287 MOP_EXP_MAT, /* ** */
304 MOP_PASTE_HORZ, /* a, b, c, ... */
305 MOP_PASTE_VERT, /* a; b; c; ... */
309 MOP_VEC_INDEX, /* x(y) */
310 MOP_VEC_ALL, /* x(:) */
311 MOP_MAT_INDEX, /* x(y,z) */
312 MOP_ROW_INDEX, /* x(y,:) */
313 MOP_COL_INDEX, /* x(:,z) */
320 MOP_EOF, /* EOF('file') */
328 struct matrix_expr **subs;
333 struct matrix_var *variable;
334 struct read_file *eof;
339 matrix_expr_destroy (struct matrix_expr *e)
346 #define F(ENUM, STRING, PROTO) case MOP_F_##ENUM:
377 for (size_t i = 0; i < e->n_subs; i++)
378 matrix_expr_destroy (e->subs[i]);
390 static struct matrix_expr *
391 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
394 struct matrix_expr *e = xmalloc (sizeof *e);
395 *e = (struct matrix_expr) {
397 .subs = xmemdup (subs, n_subs * sizeof *subs),
403 static struct matrix_expr *
404 matrix_expr_create_0 (enum matrix_op op)
406 struct matrix_expr *sub;
407 return matrix_expr_create_subs (op, &sub, 0);
410 static struct matrix_expr *
411 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
413 return matrix_expr_create_subs (op, &sub, 1);
416 static struct matrix_expr *
417 matrix_expr_create_2 (enum matrix_op op,
418 struct matrix_expr *sub0, struct matrix_expr *sub1)
420 struct matrix_expr *subs[] = { sub0, sub1 };
421 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
424 static struct matrix_expr *
425 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
426 struct matrix_expr *sub1, struct matrix_expr *sub2)
428 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
429 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
432 static struct matrix_expr *
433 matrix_expr_create_number (double number)
435 struct matrix_expr *e = xmalloc (sizeof *e);
436 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
440 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
442 static struct matrix_expr *
443 matrix_parse_curly_comma (struct matrix_state *s)
445 struct matrix_expr *lhs = matrix_parse_or_xor (s);
449 while (lex_match (s->lexer, T_COMMA))
451 struct matrix_expr *rhs = matrix_parse_or_xor (s);
454 matrix_expr_destroy (lhs);
457 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
462 static struct matrix_expr *
463 matrix_parse_curly_semi (struct matrix_state *s)
465 if (lex_token (s->lexer) == T_RCURLY)
466 return matrix_expr_create_0 (MOP_EMPTY);
468 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
472 while (lex_match (s->lexer, T_SEMICOLON))
474 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
477 matrix_expr_destroy (lhs);
480 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
485 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
486 for (size_t Y = 0; Y < (M)->size1; Y++) \
487 for (size_t X = 0; X < (M)->size2; X++) \
488 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
491 is_vector (const gsl_matrix *m)
493 return m->size1 <= 1 || m->size2 <= 1;
497 to_vector (gsl_matrix *m)
499 return (m->size1 == 1
500 ? gsl_matrix_row (m, 0).vector
501 : gsl_matrix_column (m, 0).vector);
506 matrix_eval_ABS (gsl_matrix *m)
508 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
514 matrix_eval_ALL (gsl_matrix *m)
516 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
523 matrix_eval_ANY (gsl_matrix *m)
525 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
532 matrix_eval_ARSIN (gsl_matrix *m)
534 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
540 matrix_eval_ARTAN (gsl_matrix *m)
542 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
548 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
552 for (size_t i = 0; i < n; i++)
557 gsl_matrix *block = gsl_matrix_calloc (r, c);
559 for (size_t i = 0; i < n; i++)
561 for (size_t y = 0; y < m[i]->size1; y++)
562 for (size_t x = 0; x < m[i]->size2; x++)
563 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
571 matrix_eval_CHOL (gsl_matrix *m)
573 if (!gsl_linalg_cholesky_decomp1 (m))
575 for (size_t y = 0; y < m->size1; y++)
576 for (size_t x = y + 1; x < m->size2; x++)
577 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
579 for (size_t y = 0; y < m->size1; y++)
580 for (size_t x = 0; x < y; x++)
581 gsl_matrix_set (m, y, x, 0);
586 msg (SE, _("Input to CHOL function is not positive-definite."));
592 matrix_eval_col_extremum (gsl_matrix *m, bool min)
597 return gsl_matrix_alloc (1, 0);
599 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
600 for (size_t x = 0; x < m->size2; x++)
602 double ext = gsl_matrix_get (m, 0, x);
603 for (size_t y = 1; y < m->size1; y++)
605 double value = gsl_matrix_get (m, y, x);
606 if (min ? value < ext : value > ext)
609 gsl_matrix_set (cext, 0, x, ext);
615 matrix_eval_CMAX (gsl_matrix *m)
617 return matrix_eval_col_extremum (m, false);
621 matrix_eval_CMIN (gsl_matrix *m)
623 return matrix_eval_col_extremum (m, true);
627 matrix_eval_COS (gsl_matrix *m)
629 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
635 matrix_eval_col_sum (gsl_matrix *m, bool square)
640 return gsl_matrix_alloc (1, 0);
642 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
643 for (size_t x = 0; x < m->size2; x++)
646 for (size_t y = 0; y < m->size1; y++)
648 double d = gsl_matrix_get (m, y, x);
649 sum += square ? pow2 (d) : d;
651 gsl_matrix_set (result, 0, x, sum);
657 matrix_eval_CSSQ (gsl_matrix *m)
659 return matrix_eval_col_sum (m, true);
663 matrix_eval_CSUM (gsl_matrix *m)
665 return matrix_eval_col_sum (m, false);
669 compare_double_3way (const void *a_, const void *b_)
671 const double *a = a_;
672 const double *b = b_;
673 return *a < *b ? -1 : *a > *b;
677 matrix_eval_DESIGN (gsl_matrix *m)
679 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
680 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
681 gsl_matrix_transpose_memcpy (&m2, m);
683 for (size_t y = 0; y < m2.size1; y++)
684 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
686 size_t *n = xcalloc (m2.size1, sizeof *n);
688 for (size_t i = 0; i < m2.size1; i++)
690 double *row = tmp + m2.size2 * i;
691 for (size_t j = 0; j < m2.size2; )
694 for (k = j + 1; k < m2.size2; k++)
695 if (row[j] != row[k])
697 row[n[i]++] = row[j];
702 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
707 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
709 for (size_t i = 0; i < m->size2; i++)
714 const double *unique = tmp + m2.size2 * i;
715 for (size_t j = 0; j < n[i]; j++, x++)
717 double value = unique[j];
718 for (size_t y = 0; y < m->size1; y++)
719 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
730 matrix_eval_DET (gsl_matrix *m)
732 gsl_permutation *p = gsl_permutation_alloc (m->size1);
734 gsl_linalg_LU_decomp (m, p, &signum);
735 gsl_permutation_free (p);
736 return gsl_linalg_LU_det (m, signum);
740 matrix_eval_DIAG (gsl_matrix *m)
742 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
743 for (size_t i = 0; i < diag->size1; i++)
744 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
749 is_symmetric (const gsl_matrix *m)
751 if (m->size1 != m->size2)
754 for (size_t y = 0; y < m->size1; y++)
755 for (size_t x = 0; x < y; x++)
756 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
763 compare_double_desc (const void *a_, const void *b_)
765 const double *a = a_;
766 const double *b = b_;
767 return *a > *b ? -1 : *a < *b;
771 matrix_eval_EVAL (gsl_matrix *m)
773 if (!is_symmetric (m))
775 msg (SE, _("Argument of EVAL must be symmetric."));
779 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
780 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
781 gsl_vector v_eval = to_vector (eval);
782 gsl_eigen_symm (m, &v_eval, w);
783 gsl_eigen_symm_free (w);
785 assert (v_eval.stride == 1);
786 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
792 matrix_eval_EXP (gsl_matrix *m)
794 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
799 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
802 Charl Linssen <charl@itfromb.it>
806 matrix_eval_GINV (gsl_matrix *A)
811 gsl_matrix *tmp_mat = NULL;
814 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
815 tmp_mat = gsl_matrix_alloc (m, n);
816 gsl_matrix_transpose_memcpy (tmp_mat, A);
824 gsl_matrix *V = gsl_matrix_alloc (m, m);
825 gsl_vector *u = gsl_vector_alloc (m);
827 gsl_vector *tmp_vec = gsl_vector_alloc (m);
828 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
829 gsl_vector_free (tmp_vec);
832 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
833 gsl_matrix_set_zero (Sigma_pinv);
834 double cutoff = 1e-15 * gsl_vector_max (u);
836 for (size_t i = 0; i < m; ++i)
838 double x = gsl_vector_get (u, i);
839 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
842 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
843 gsl_matrix *U = gsl_matrix_calloc (n, n);
844 for (size_t i = 0; i < n; i++)
845 for (size_t j = 0; j < m; j++)
846 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
848 /* two dot products to obtain pseudoinverse */
849 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
850 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
855 A_pinv = gsl_matrix_alloc (n, m);
856 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
860 A_pinv = gsl_matrix_alloc (m, n);
861 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
864 gsl_matrix_free (tmp_mat);
865 gsl_matrix_free (tmp_mat2);
867 gsl_matrix_free (Sigma_pinv);
881 grade_compare_3way (const void *a_, const void *b_)
883 const struct grade *a = a_;
884 const struct grade *b = b_;
886 return (a->value < b->value ? -1
887 : a->value > b->value ? 1
895 matrix_eval_GRADE (gsl_matrix *m)
897 size_t n = m->size1 * m->size2;
898 struct grade *grades = xmalloc (n * sizeof *grades);
901 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
902 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
903 qsort (grades, n, sizeof *grades, grade_compare_3way);
905 for (size_t i = 0; i < n; i++)
906 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
914 dot (gsl_vector *a, gsl_vector *b)
917 for (size_t i = 0; i < a->size; i++)
918 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
923 norm2 (gsl_vector *v)
926 for (size_t i = 0; i < v->size; i++)
927 result += pow2 (gsl_vector_get (v, i));
934 return sqrt (norm2 (v));
938 matrix_eval_GSCH (gsl_matrix *v)
940 if (v->size2 < v->size1)
942 msg (SE, _("GSCH requires its argument to have at least as many columns "
943 "as rows, but it has dimensions %zu×%zu."),
947 if (!v->size1 || !v->size2)
950 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
952 for (size_t vx = 0; vx < v->size2; vx++)
954 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
955 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
957 gsl_vector_memcpy (&u_i, &v_i);
958 for (size_t j = 0; j < ux; j++)
960 gsl_vector u_j = gsl_matrix_column (u, j).vector;
961 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
962 for (size_t k = 0; k < u_i.size; k++)
963 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
964 - scale * gsl_vector_get (&u_j, k)));
967 double len = norm (&u_i);
970 gsl_vector_scale (&u_i, 1.0 / len);
971 if (++ux >= v->size1)
978 msg (SE, _("%zu×%zu argument to GSCH contains only "
979 "%zu linearly independent columns."),
980 v->size1, v->size2, ux);
990 matrix_eval_IDENT (double s1, double s2)
992 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
994 msg (SE, _("Arguments to IDENT must be non-negative integers."));
998 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
999 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1005 invert_matrix (gsl_matrix *x)
1007 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1009 gsl_linalg_LU_decomp (x, p, &signum);
1010 gsl_linalg_LU_invx (x, p);
1011 gsl_permutation_free (p);
1015 matrix_eval_INV (gsl_matrix *m)
1022 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1024 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1025 a->size2 * b->size2);
1027 for (size_t ar = 0; ar < a->size1; ar++)
1028 for (size_t br = 0; br < b->size1; br++, y++)
1031 for (size_t ac = 0; ac < a->size2; ac++)
1032 for (size_t bc = 0; bc < b->size2; bc++, x++)
1034 double av = gsl_matrix_get (a, ar, ac);
1035 double bv = gsl_matrix_get (b, br, bc);
1036 gsl_matrix_set (k, y, x, av * bv);
1043 matrix_eval_LG10 (gsl_matrix *m)
1045 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1051 matrix_eval_LN (gsl_matrix *m)
1053 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1059 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1061 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1064 for (size_t i = 1; i <= n * n; i++)
1066 gsl_matrix_set (m, y, x, i);
1068 size_t y1 = !y ? n - 1 : y - 1;
1069 size_t x1 = x + 1 >= n ? 0 : x + 1;
1070 if (gsl_matrix_get (m, y1, x1) == 0)
1076 y = y + 1 >= n ? 0 : y + 1;
1081 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1083 double a = gsl_matrix_get (m, y1, x1);
1084 double b = gsl_matrix_get (m, y2, x2);
1085 gsl_matrix_set (m, y1, x1, b);
1086 gsl_matrix_set (m, y2, x2, a);
1090 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1094 /* A. Umar, "On the Construction of Even Order Magic Squares",
1095 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1097 for (size_t i = 1; i <= n * n / 2; i++)
1099 gsl_matrix_set (m, y, x, i);
1109 for (size_t i = n * n; i > n * n / 2; i--)
1111 gsl_matrix_set (m, y, x, i);
1119 for (size_t y = 0; y < n; y++)
1120 for (size_t x = 0; x < n / 2; x++)
1122 unsigned int d = gsl_matrix_get (m, y, x);
1123 if (d % 2 != (y < n / 2))
1124 magic_exchange (m, y, x, y, n - x - 1);
1129 size_t x1 = n / 2 - 1;
1131 magic_exchange (m, y1, x1, y2, x1);
1132 magic_exchange (m, y1, x2, y2, x2);
1136 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1138 /* A. Umar, "On the Construction of Even Order Magic Squares",
1139 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1143 for (size_t i = 1; ; i++)
1145 gsl_matrix_set (m, y, x, i);
1146 if (++y == n / 2 - 1)
1158 for (size_t i = n * n; ; i--)
1160 gsl_matrix_set (m, y, x, i);
1161 if (++y == n / 2 - 1)
1170 for (size_t y = 0; y < n; y++)
1171 if (y != n / 2 - 1 && y != n / 2)
1172 for (size_t x = 0; x < n / 2; x++)
1174 unsigned int d = gsl_matrix_get (m, y, x);
1175 if (d % 2 != (y < n / 2))
1176 magic_exchange (m, y, x, y, n - x - 1);
1179 size_t a0 = (n * n - 2 * n) / 2 + 1;
1180 for (size_t i = 0; i < n / 2; i++)
1183 gsl_matrix_set (m, n / 2, i, a);
1184 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1186 for (size_t i = 0; i < n / 2; i++)
1188 size_t a = a0 + i + n / 2;
1189 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1190 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1192 for (size_t x = 1; x < n / 2; x += 2)
1193 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1194 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1195 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1196 size_t x1 = n / 2 - 2;
1197 size_t x2 = n / 2 + 1;
1198 size_t y1 = n / 2 - 2;
1199 size_t y2 = n / 2 + 1;
1200 magic_exchange (m, y1, x1, y2, x1);
1201 magic_exchange (m, y1, x2, y2, x2);
1205 matrix_eval_MAGIC (double n_)
1207 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1209 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1214 gsl_matrix *m = gsl_matrix_calloc (n, n);
1216 matrix_eval_MAGIC_odd (m, n);
1218 matrix_eval_MAGIC_singly_even (m, n);
1220 matrix_eval_MAGIC_doubly_even (m, n);
1225 matrix_eval_MAKE (double r, double c, double value)
1227 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1229 msg (SE, _("Size arguments to MAKE must be integers."));
1233 gsl_matrix *m = gsl_matrix_alloc (r, c);
1234 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1240 matrix_eval_MDIAG (gsl_vector *v)
1242 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1243 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1244 gsl_vector_memcpy (&diagonal, v);
1249 matrix_eval_MMAX (gsl_matrix *m)
1251 return gsl_matrix_max (m);
1255 matrix_eval_MMIN (gsl_matrix *m)
1257 return gsl_matrix_min (m);
1261 matrix_eval_MOD (gsl_matrix *m, double divisor)
1265 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1269 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1270 *d = fmod (*d, divisor);
1275 matrix_eval_MSSQ (gsl_matrix *m)
1278 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1284 matrix_eval_MSUM (gsl_matrix *m)
1287 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1293 matrix_eval_NCOL (gsl_matrix *m)
1299 matrix_eval_NROW (gsl_matrix *m)
1305 matrix_eval_RANK (gsl_matrix *m)
1307 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1308 gsl_linalg_QR_decomp (m, tau);
1309 gsl_vector_free (tau);
1311 return gsl_linalg_QRPT_rank (m, -1);
1315 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1317 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1319 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1324 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1326 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1327 "product of matrix dimensions."));
1331 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1334 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1336 gsl_matrix_set (dst, y1, x1, *d);
1347 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1352 return gsl_matrix_alloc (0, 1);
1354 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1355 for (size_t y = 0; y < m->size1; y++)
1357 double ext = gsl_matrix_get (m, y, 0);
1358 for (size_t x = 1; x < m->size2; x++)
1360 double value = gsl_matrix_get (m, y, x);
1361 if (min ? value < ext : value > ext)
1364 gsl_matrix_set (rext, y, 0, ext);
1370 matrix_eval_RMAX (gsl_matrix *m)
1372 return matrix_eval_row_extremum (m, false);
1376 matrix_eval_RMIN (gsl_matrix *m)
1378 return matrix_eval_row_extremum (m, true);
1382 matrix_eval_RND (gsl_matrix *m)
1384 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1396 rank_compare_3way (const void *a_, const void *b_)
1398 const struct rank *a = a_;
1399 const struct rank *b = b_;
1401 return a->value < b->value ? -1 : a->value > b->value;
1405 matrix_eval_RNKORDER (gsl_matrix *m)
1407 size_t n = m->size1 * m->size2;
1408 struct rank *ranks = xmalloc (n * sizeof *ranks);
1410 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1411 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1412 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1414 for (size_t i = 0; i < n; )
1417 for (j = i + 1; j < n; j++)
1418 if (ranks[i].value != ranks[j].value)
1421 double rank = (i + j + 1.0) / 2.0;
1422 for (size_t k = i; k < j; k++)
1423 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1434 matrix_eval_row_sum (gsl_matrix *m, bool square)
1439 return gsl_matrix_alloc (0, 1);
1441 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1442 for (size_t y = 0; y < m->size1; y++)
1445 for (size_t x = 0; x < m->size2; x++)
1447 double d = gsl_matrix_get (m, y, x);
1448 sum += square ? pow2 (d) : d;
1450 gsl_matrix_set (result, y, 0, sum);
1456 matrix_eval_RSSQ (gsl_matrix *m)
1458 return matrix_eval_row_sum (m, true);
1462 matrix_eval_RSUM (gsl_matrix *m)
1464 return matrix_eval_row_sum (m, false);
1468 matrix_eval_SIN (gsl_matrix *m)
1470 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1476 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1478 if (m1->size1 != m2->size1)
1480 msg (SE, _("SOLVE requires its arguments to have the same number of "
1481 "rows, but the first argument has dimensions %zu×%zu and "
1482 "the second %zu×%zu."),
1483 m1->size1, m1->size2,
1484 m2->size1, m2->size2);
1488 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1489 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1491 gsl_linalg_LU_decomp (m1, p, &signum);
1492 for (size_t i = 0; i < m2->size2; i++)
1494 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1495 gsl_vector xi = gsl_matrix_column (x, i).vector;
1496 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1498 gsl_permutation_free (p);
1503 matrix_eval_SQRT (gsl_matrix *m)
1505 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1509 msg (SE, _("Argument to SQRT must be nonnegative."));
1518 matrix_eval_SSCP (gsl_matrix *m)
1520 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1521 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1526 matrix_eval_SVAL (gsl_matrix *m)
1528 gsl_matrix *tmp_mat = NULL;
1529 if (m->size2 > m->size1)
1531 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1532 gsl_matrix_transpose_memcpy (tmp_mat, m);
1537 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1538 gsl_vector *S = gsl_vector_alloc (m->size2);
1539 gsl_vector *work = gsl_vector_alloc (m->size2);
1540 gsl_linalg_SV_decomp (m, V, S, work);
1542 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1543 for (size_t i = 0; i < m->size2; i++)
1544 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1546 gsl_matrix_free (V);
1547 gsl_vector_free (S);
1548 gsl_vector_free (work);
1549 gsl_matrix_free (tmp_mat);
1555 matrix_eval_SWEEP (gsl_matrix *m, double d)
1557 if (d < 1 || d > SIZE_MAX)
1559 msg (SE, _("Scalar argument to SWEEP must be integer."));
1563 if (k >= MIN (m->size1, m->size2))
1565 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1566 "equal to the smaller of the matrix argument's rows and "
1571 double m_kk = gsl_matrix_get (m, k, k);
1572 if (fabs (m_kk) > 1e-19)
1574 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1575 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1577 double m_ij = gsl_matrix_get (m, i, j);
1578 double m_ik = gsl_matrix_get (m, i, k);
1579 double m_kj = gsl_matrix_get (m, k, j);
1580 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1589 for (size_t i = 0; i < m->size1; i++)
1591 gsl_matrix_set (m, i, k, 0);
1592 gsl_matrix_set (m, k, i, 0);
1599 matrix_eval_TRACE (gsl_matrix *m)
1602 size_t n = MIN (m->size1, m->size2);
1603 for (size_t i = 0; i < n; i++)
1604 sum += gsl_matrix_get (m, i, i);
1609 matrix_eval_T (gsl_matrix *m)
1611 return matrix_eval_TRANSPOS (m);
1615 matrix_eval_TRANSPOS (gsl_matrix *m)
1617 if (m->size1 == m->size2)
1619 gsl_matrix_transpose (m);
1624 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1625 gsl_matrix_transpose_memcpy (t, m);
1631 matrix_eval_TRUNC (gsl_matrix *m)
1633 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1639 matrix_eval_UNIFORM (double r_, double c_)
1641 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1643 msg (SE, _("Arguments to UNIFORM must be integers."));
1648 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1650 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1654 gsl_matrix *m = gsl_matrix_alloc (r, c);
1655 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1656 *d = gsl_ran_flat (get_rng (), 0, 1);
1660 struct matrix_function
1664 size_t min_args, max_args;
1667 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1669 static const struct matrix_function *
1670 matrix_parse_function_name (const char *token)
1672 static const struct matrix_function functions[] =
1674 #define F(ENUM, STRING, PROTO) \
1675 { #ENUM, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1679 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1681 for (size_t i = 0; i < N_FUNCTIONS; i++)
1683 if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1684 return &functions[i];
1689 static struct read_file *
1690 read_file_create (struct matrix_state *s, struct file_handle *fh)
1692 for (size_t i = 0; i < s->n_read_files; i++)
1694 struct read_file *rf = s->read_files[i];
1702 struct read_file *rf = xmalloc (sizeof *rf);
1703 *rf = (struct read_file) { .file = fh };
1705 s->read_files = xrealloc (s->read_files,
1706 (s->n_read_files + 1) * sizeof *s->read_files);
1707 s->read_files[s->n_read_files++] = rf;
1711 static struct dfm_reader *
1712 read_file_open (struct read_file *rf)
1715 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1720 read_file_destroy (struct read_file *rf)
1724 fh_unref (rf->file);
1725 dfm_close_reader (rf->reader);
1726 free (rf->encoding);
1731 static struct write_file *
1732 write_file_create (struct matrix_state *s, struct file_handle *fh)
1734 for (size_t i = 0; i < s->n_write_files; i++)
1736 struct write_file *wf = s->write_files[i];
1744 struct write_file *wf = xmalloc (sizeof *wf);
1745 *wf = (struct write_file) { .file = fh };
1747 s->write_files = xrealloc (s->write_files,
1748 (s->n_write_files + 1) * sizeof *s->write_files);
1749 s->write_files[s->n_write_files++] = wf;
1753 static struct dfm_writer *
1754 write_file_open (struct write_file *wf)
1757 wf->writer = dfm_open_writer (wf->file, wf->encoding);
1762 write_file_destroy (struct write_file *wf)
1768 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
1769 wf->held->s.ss.length);
1770 u8_line_destroy (wf->held);
1774 fh_unref (wf->file);
1775 dfm_close_writer (wf->writer);
1776 free (wf->encoding);
1782 matrix_parse_function (struct matrix_state *s, const char *token,
1783 struct matrix_expr **exprp)
1786 if (lex_next_token (s->lexer, 1) != T_LPAREN)
1789 if (lex_match_id (s->lexer, "EOF"))
1792 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1796 if (!lex_force_match (s->lexer, T_RPAREN))
1802 struct read_file *rf = read_file_create (s, fh);
1804 struct matrix_expr *e = xmalloc (sizeof *e);
1805 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1810 const struct matrix_function *f = matrix_parse_function_name (token);
1814 lex_get_n (s->lexer, 2);
1816 struct matrix_expr *e = xmalloc (sizeof *e);
1817 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1819 size_t allocated_subs = 0;
1822 struct matrix_expr *sub = matrix_parse_expr (s);
1826 if (e->n_subs >= allocated_subs)
1827 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1828 e->subs[e->n_subs++] = sub;
1830 while (lex_match (s->lexer, T_COMMA));
1831 if (!lex_force_match (s->lexer, T_RPAREN))
1838 matrix_expr_destroy (e);
1842 static struct matrix_expr *
1843 matrix_parse_primary (struct matrix_state *s)
1845 if (lex_is_number (s->lexer))
1847 double number = lex_number (s->lexer);
1849 return matrix_expr_create_number (number);
1851 else if (lex_is_string (s->lexer))
1853 char string[sizeof (double)];
1854 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1858 memcpy (&number, string, sizeof number);
1859 return matrix_expr_create_number (number);
1861 else if (lex_match (s->lexer, T_LPAREN))
1863 struct matrix_expr *e = matrix_parse_or_xor (s);
1864 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1866 matrix_expr_destroy (e);
1871 else if (lex_match (s->lexer, T_LCURLY))
1873 struct matrix_expr *e = matrix_parse_curly_semi (s);
1874 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1876 matrix_expr_destroy (e);
1881 else if (lex_token (s->lexer) == T_ID)
1883 struct matrix_expr *retval;
1884 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1887 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1890 lex_error (s->lexer, _("Unknown variable %s."),
1891 lex_tokcstr (s->lexer));
1896 struct matrix_expr *e = xmalloc (sizeof *e);
1897 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1900 else if (lex_token (s->lexer) == T_ALL)
1902 struct matrix_expr *retval;
1903 if (matrix_parse_function (s, "ALL", &retval))
1907 lex_error (s->lexer, NULL);
1911 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1914 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1916 if (lex_match (s->lexer, T_COLON))
1923 *indexp = matrix_parse_expr (s);
1924 return *indexp != NULL;
1928 static struct matrix_expr *
1929 matrix_parse_postfix (struct matrix_state *s)
1931 struct matrix_expr *lhs = matrix_parse_primary (s);
1932 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1935 struct matrix_expr *i0;
1936 if (!matrix_parse_index_expr (s, &i0))
1938 matrix_expr_destroy (lhs);
1941 if (lex_match (s->lexer, T_RPAREN))
1943 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1944 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1945 else if (lex_match (s->lexer, T_COMMA))
1947 struct matrix_expr *i1;
1948 if (!matrix_parse_index_expr (s, &i1)
1949 || !lex_force_match (s->lexer, T_RPAREN))
1951 matrix_expr_destroy (lhs);
1952 matrix_expr_destroy (i0);
1953 matrix_expr_destroy (i1);
1956 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1957 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1958 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1963 lex_error_expecting (s->lexer, "`)'", "`,'");
1968 static struct matrix_expr *
1969 matrix_parse_unary (struct matrix_state *s)
1971 if (lex_match (s->lexer, T_DASH))
1973 struct matrix_expr *lhs = matrix_parse_unary (s);
1974 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1976 else if (lex_match (s->lexer, T_PLUS))
1977 return matrix_parse_unary (s);
1979 return matrix_parse_postfix (s);
1982 static struct matrix_expr *
1983 matrix_parse_seq (struct matrix_state *s)
1985 struct matrix_expr *start = matrix_parse_unary (s);
1986 if (!start || !lex_match (s->lexer, T_COLON))
1989 struct matrix_expr *end = matrix_parse_unary (s);
1992 matrix_expr_destroy (start);
1996 if (lex_match (s->lexer, T_COLON))
1998 struct matrix_expr *increment = matrix_parse_unary (s);
2001 matrix_expr_destroy (start);
2002 matrix_expr_destroy (end);
2005 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2008 return matrix_expr_create_2 (MOP_SEQ, start, end);
2011 static struct matrix_expr *
2012 matrix_parse_exp (struct matrix_state *s)
2014 struct matrix_expr *lhs = matrix_parse_seq (s);
2021 if (lex_match (s->lexer, T_EXP))
2023 else if (lex_match_phrase (s->lexer, "&**"))
2028 struct matrix_expr *rhs = matrix_parse_seq (s);
2031 matrix_expr_destroy (lhs);
2034 lhs = matrix_expr_create_2 (op, lhs, rhs);
2038 static struct matrix_expr *
2039 matrix_parse_mul_div (struct matrix_state *s)
2041 struct matrix_expr *lhs = matrix_parse_exp (s);
2048 if (lex_match (s->lexer, T_ASTERISK))
2050 else if (lex_match (s->lexer, T_SLASH))
2052 else if (lex_match_phrase (s->lexer, "&*"))
2054 else if (lex_match_phrase (s->lexer, "&/"))
2059 struct matrix_expr *rhs = matrix_parse_exp (s);
2062 matrix_expr_destroy (lhs);
2065 lhs = matrix_expr_create_2 (op, lhs, rhs);
2069 static struct matrix_expr *
2070 matrix_parse_add_sub (struct matrix_state *s)
2072 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2079 if (lex_match (s->lexer, T_PLUS))
2081 else if (lex_match (s->lexer, T_DASH))
2083 else if (lex_token (s->lexer) == T_NEG_NUM)
2088 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2091 matrix_expr_destroy (lhs);
2094 lhs = matrix_expr_create_2 (op, lhs, rhs);
2098 static struct matrix_expr *
2099 matrix_parse_relational (struct matrix_state *s)
2101 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2108 if (lex_match (s->lexer, T_GT))
2110 else if (lex_match (s->lexer, T_GE))
2112 else if (lex_match (s->lexer, T_LT))
2114 else if (lex_match (s->lexer, T_LE))
2116 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2118 else if (lex_match (s->lexer, T_NE))
2123 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2126 matrix_expr_destroy (lhs);
2129 lhs = matrix_expr_create_2 (op, lhs, rhs);
2133 static struct matrix_expr *
2134 matrix_parse_not (struct matrix_state *s)
2136 if (lex_match (s->lexer, T_NOT))
2138 struct matrix_expr *lhs = matrix_parse_not (s);
2139 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2142 return matrix_parse_relational (s);
2145 static struct matrix_expr *
2146 matrix_parse_and (struct matrix_state *s)
2148 struct matrix_expr *lhs = matrix_parse_not (s);
2152 while (lex_match (s->lexer, T_AND))
2154 struct matrix_expr *rhs = matrix_parse_not (s);
2157 matrix_expr_destroy (lhs);
2160 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2165 static struct matrix_expr *
2166 matrix_parse_or_xor (struct matrix_state *s)
2168 struct matrix_expr *lhs = matrix_parse_and (s);
2175 if (lex_match (s->lexer, T_OR))
2177 else if (lex_match_id (s->lexer, "XOR"))
2182 struct matrix_expr *rhs = matrix_parse_and (s);
2185 matrix_expr_destroy (lhs);
2188 lhs = matrix_expr_create_2 (op, lhs, rhs);
2192 static struct matrix_expr *
2193 matrix_parse_expr (struct matrix_state *s)
2195 return matrix_parse_or_xor (s);
2198 /* Expression evaluation. */
2201 matrix_op_eval (enum matrix_op op, double a, double b)
2205 case MOP_ADD_ELEMS: return a + b;
2206 case MOP_SUB_ELEMS: return a - b;
2207 case MOP_MUL_ELEMS: return a * b;
2208 case MOP_DIV_ELEMS: return a / b;
2209 case MOP_EXP_ELEMS: return pow (a, b);
2210 case MOP_GT: return a > b;
2211 case MOP_GE: return a >= b;
2212 case MOP_LT: return a < b;
2213 case MOP_LE: return a <= b;
2214 case MOP_EQ: return a == b;
2215 case MOP_NE: return a != b;
2216 case MOP_AND: return (a > 0) && (b > 0);
2217 case MOP_OR: return (a > 0) || (b > 0);
2218 case MOP_XOR: return (a > 0) != (b > 0);
2220 #define F(ENUM, STRING, PROTO) case MOP_F_##ENUM:
2229 case MOP_PASTE_HORZ:
2230 case MOP_PASTE_VERT:
2246 matrix_op_name (enum matrix_op op)
2250 case MOP_ADD_ELEMS: return "+";
2251 case MOP_SUB_ELEMS: return "-";
2252 case MOP_MUL_ELEMS: return "&*";
2253 case MOP_DIV_ELEMS: return "&/";
2254 case MOP_EXP_ELEMS: return "&**";
2255 case MOP_GT: return ">";
2256 case MOP_GE: return ">=";
2257 case MOP_LT: return "<";
2258 case MOP_LE: return "<=";
2259 case MOP_EQ: return "=";
2260 case MOP_NE: return "<>";
2261 case MOP_AND: return "AND";
2262 case MOP_OR: return "OR";
2263 case MOP_XOR: return "XOR";
2265 #define F(ENUM, STRING, PROTO) case MOP_F_##ENUM:
2274 case MOP_PASTE_HORZ:
2275 case MOP_PASTE_VERT:
2291 is_scalar (const gsl_matrix *m)
2293 return m->size1 == 1 && m->size2 == 1;
2297 to_scalar (const gsl_matrix *m)
2299 assert (is_scalar (m));
2300 return gsl_matrix_get (m, 0, 0);
2304 matrix_expr_evaluate_elementwise (enum matrix_op op,
2305 gsl_matrix *a, gsl_matrix *b)
2309 double be = to_scalar (b);
2310 for (size_t r = 0; r < a->size1; r++)
2311 for (size_t c = 0; c < a->size2; c++)
2313 double *ae = gsl_matrix_ptr (a, r, c);
2314 *ae = matrix_op_eval (op, *ae, be);
2318 else if (is_scalar (a))
2320 double ae = to_scalar (a);
2321 for (size_t r = 0; r < b->size1; r++)
2322 for (size_t c = 0; c < b->size2; c++)
2324 double *be = gsl_matrix_ptr (b, r, c);
2325 *be = matrix_op_eval (op, ae, *be);
2329 else if (a->size1 == b->size1 && a->size2 == b->size2)
2331 for (size_t r = 0; r < a->size1; r++)
2332 for (size_t c = 0; c < a->size2; c++)
2334 double *ae = gsl_matrix_ptr (a, r, c);
2335 double be = gsl_matrix_get (b, r, c);
2336 *ae = matrix_op_eval (op, *ae, be);
2342 msg (SE, _("Operands to %s must have the same dimensions or one "
2343 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
2344 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2350 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2352 if (is_scalar (a) || is_scalar (b))
2353 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2355 if (a->size2 != b->size1)
2357 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
2358 "not conformable for multiplication."),
2359 a->size1, a->size2, b->size1, b->size2);
2363 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2364 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2369 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2371 gsl_matrix *tmp = *a;
2377 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2380 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2381 swap_matrix (z, tmp);
2385 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2387 mul_matrix (x, *x, *x, tmp);
2391 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2394 if (x->size1 != x->size2)
2396 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2397 "the left-hand size, not one with dimensions %zu×%zu."),
2398 x->size1, x->size2);
2403 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2404 "right-hand side, not a matrix with dimensions %zu×%zu."),
2405 b->size1, b->size2);
2408 double bf = to_scalar (b);
2409 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2411 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2412 "or outside the valid range."), bf);
2417 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2419 gsl_matrix_set_identity (y);
2423 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2425 for (unsigned long int n = labs (bl); n > 1; n /= 2)
2428 mul_matrix (&y, x, y, &t);
2429 square_matrix (&x, &t);
2432 square_matrix (&x, &t);
2434 mul_matrix (&y, x, y, &t);
2438 /* Garbage collection.
2440 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2441 a permutation of them. We are returning one of them; that one must not be
2442 destroyed. We must not destroy 'x_' because the caller owns it. */
2444 gsl_matrix_free (y_);
2446 gsl_matrix_free (t_);
2452 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2455 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2457 msg (SE, _("All operands of : operator must be scalars."));
2461 long int start = to_scalar (start_);
2462 long int end = to_scalar (end_);
2463 long int by = by_ ? to_scalar (by_) : 1;
2467 msg (SE, _("The increment operand to : must be nonzero."));
2471 long int n = (end >= start && by > 0 ? (end - start + by) / by
2472 : end <= start && by < 0 ? (start - end - by) / -by
2474 gsl_matrix *m = gsl_matrix_alloc (1, n);
2475 for (long int i = 0; i < n; i++)
2476 gsl_matrix_set (m, 0, i, start + i * by);
2481 matrix_expr_evaluate_not (gsl_matrix *a)
2483 for (size_t r = 0; r < a->size1; r++)
2484 for (size_t c = 0; c < a->size2; c++)
2486 double *ae = gsl_matrix_ptr (a, r, c);
2493 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2495 if (a->size1 != b->size1)
2497 if (!a->size1 || !a->size2)
2499 else if (!b->size1 || !b->size2)
2502 msg (SE, _("All columns in a matrix must have the same number of rows, "
2503 "but this tries to paste matrices with %zu and %zu rows."),
2504 a->size1, b->size1);
2508 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2509 for (size_t y = 0; y < a->size1; y++)
2511 for (size_t x = 0; x < a->size2; x++)
2512 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2513 for (size_t x = 0; x < b->size2; x++)
2514 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2520 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2522 if (a->size2 != b->size2)
2524 if (!a->size1 || !a->size2)
2526 else if (!b->size1 || !b->size2)
2529 msg (SE, _("All rows in a matrix must have the same number of columns, "
2530 "but this tries to stack matrices with %zu and %zu columns."),
2531 a->size2, b->size2);
2535 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2536 for (size_t x = 0; x < a->size2; x++)
2538 for (size_t y = 0; y < a->size1; y++)
2539 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2540 for (size_t y = 0; y < b->size1; y++)
2541 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2547 matrix_to_vector (gsl_matrix *m)
2550 gsl_vector v = to_vector (m);
2551 assert (v.block == m->block || !v.block);
2555 gsl_matrix_free (m);
2556 return xmemdup (&v, sizeof v);
2570 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2573 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2574 enum index_type index_type, size_t other_size,
2575 struct index_vector *iv)
2584 msg (SE, _("Vector index must be scalar or vector, not a "
2586 m->size1, m->size2);
2590 msg (SE, _("Matrix row index must be scalar or vector, not a "
2592 m->size1, m->size2);
2596 msg (SE, _("Matrix column index must be scalar or vector, not a "
2598 m->size1, m->size2);
2604 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2605 *iv = (struct index_vector) {
2606 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2609 for (size_t i = 0; i < v.size; i++)
2611 double index = gsl_vector_get (&v, i);
2612 if (index < 1 || index >= size + 1)
2617 msg (SE, _("Index %g is out of range for vector "
2618 "with %zu elements."), index, size);
2622 msg (SE, _("%g is not a valid row index for "
2623 "a %zu×%zu matrix."),
2624 index, size, other_size);
2628 msg (SE, _("%g is not a valid column index for "
2629 "a %zu×%zu matrix."),
2630 index, other_size, size);
2637 iv->indexes[i] = index - 1;
2643 *iv = (struct index_vector) {
2644 .indexes = xnmalloc (size, sizeof *iv->indexes),
2647 for (size_t i = 0; i < size; i++)
2654 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2656 if (!is_vector (sm))
2658 msg (SE, _("Vector index operator may not be applied to "
2659 "a %zu×%zu matrix."),
2660 sm->size1, sm->size2);
2668 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2670 if (!matrix_expr_evaluate_vec_all (sm))
2673 gsl_vector sv = to_vector (sm);
2674 struct index_vector iv;
2675 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2678 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2679 sm->size1 == 1 ? iv.n : 1);
2680 gsl_vector dv = to_vector (dm);
2681 for (size_t dx = 0; dx < iv.n; dx++)
2683 size_t sx = iv.indexes[dx];
2684 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2692 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2695 struct index_vector iv0;
2696 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2699 struct index_vector iv1;
2700 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2707 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2708 for (size_t dy = 0; dy < iv0.n; dy++)
2710 size_t sy = iv0.indexes[dy];
2712 for (size_t dx = 0; dx < iv1.n; dx++)
2714 size_t sx = iv1.indexes[dx];
2715 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2723 #define F(ENUM, STRING, PROTO) \
2724 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
2725 const char *name, gsl_matrix *[], size_t, \
2726 matrix_proto_##PROTO *);
2731 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2733 if (!is_scalar (subs[index]))
2735 msg (SE, _("Function %s argument %zu must be a scalar, "
2736 "not a %zu×%zu matrix."),
2737 name, index + 1, subs[index]->size1, subs[index]->size2);
2744 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2746 if (!is_vector (subs[index]))
2748 msg (SE, _("Function %s argument %zu must be a vector, "
2749 "not a %zu×%zu matrix."),
2750 name, index + 1, subs[index]->size1, subs[index]->size2);
2757 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2759 for (size_t i = 0; i < n_subs; i++)
2761 if (!check_scalar_arg (name, subs, i))
2763 d[i] = to_scalar (subs[i]);
2769 matrix_expr_evaluate_m_d (const char *name,
2770 gsl_matrix *subs[], size_t n_subs,
2771 matrix_proto_m_d *f)
2773 assert (n_subs == 1);
2776 return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2780 matrix_expr_evaluate_m_dd (const char *name,
2781 gsl_matrix *subs[], size_t n_subs,
2782 matrix_proto_m_dd *f)
2784 assert (n_subs == 2);
2787 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2791 matrix_expr_evaluate_m_ddd (const char *name,
2792 gsl_matrix *subs[], size_t n_subs,
2793 matrix_proto_m_ddd *f)
2795 assert (n_subs == 3);
2798 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2802 matrix_expr_evaluate_m_m (const char *name UNUSED,
2803 gsl_matrix *subs[], size_t n_subs,
2804 matrix_proto_m_m *f)
2806 assert (n_subs == 1);
2811 matrix_expr_evaluate_m_md (const char *name UNUSED,
2812 gsl_matrix *subs[], size_t n_subs,
2813 matrix_proto_m_md *f)
2815 assert (n_subs == 2);
2816 if (!check_scalar_arg (name, subs, 1))
2818 return f (subs[0], to_scalar (subs[1]));
2822 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2823 gsl_matrix *subs[], size_t n_subs,
2824 matrix_proto_m_mdd *f)
2826 assert (n_subs == 3);
2827 if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2829 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2833 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2834 gsl_matrix *subs[], size_t n_subs,
2835 matrix_proto_m_mm *f)
2837 assert (n_subs == 2);
2838 return f (subs[0], subs[1]);
2842 matrix_expr_evaluate_m_v (const char *name,
2843 gsl_matrix *subs[], size_t n_subs,
2844 matrix_proto_m_v *f)
2846 assert (n_subs == 1);
2847 if (!check_vector_arg (name, subs, 0))
2849 gsl_vector v = to_vector (subs[0]);
2854 matrix_expr_evaluate_d_m (const char *name UNUSED,
2855 gsl_matrix *subs[], size_t n_subs,
2856 matrix_proto_d_m *f)
2858 assert (n_subs == 1);
2860 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2861 gsl_matrix_set (m, 0, 0, f (subs[0]));
2866 matrix_expr_evaluate_m_any (const char *name UNUSED,
2867 gsl_matrix *subs[], size_t n_subs,
2868 matrix_proto_m_any *f)
2870 return f (subs, n_subs);
2874 matrix_expr_evaluate_IDENT (const char *name,
2875 gsl_matrix *subs[], size_t n_subs,
2876 matrix_proto_IDENT *f)
2878 assert (n_subs <= 2);
2881 if (!to_scalar_args (name, subs, n_subs, d))
2884 return f (d[0], d[n_subs - 1]);
2888 matrix_expr_evaluate (const struct matrix_expr *e)
2890 if (e->op == MOP_NUMBER)
2892 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2893 gsl_matrix_set (m, 0, 0, e->number);
2896 else if (e->op == MOP_VARIABLE)
2898 const gsl_matrix *src = e->variable->value;
2901 msg (SE, _("Uninitialized variable %s used in expression."),
2906 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2907 gsl_matrix_memcpy (dst, src);
2910 else if (e->op == MOP_EOF)
2912 struct dfm_reader *reader = read_file_open (e->eof);
2913 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2914 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2918 enum { N_LOCAL = 3 };
2919 gsl_matrix *local_subs[N_LOCAL];
2920 gsl_matrix **subs = (e->n_subs < N_LOCAL
2922 : xmalloc (e->n_subs * sizeof *subs));
2924 for (size_t i = 0; i < e->n_subs; i++)
2926 subs[i] = matrix_expr_evaluate (e->subs[i]);
2929 for (size_t j = 0; j < i; j++)
2930 gsl_matrix_free (subs[j]);
2931 if (subs != local_subs)
2937 gsl_matrix *result = NULL;
2940 #define F(ENUM, STRING, PROTO) \
2941 case MOP_F_##ENUM: \
2942 result = matrix_expr_evaluate_##PROTO (STRING, subs, e->n_subs, \
2943 matrix_eval_##ENUM); \
2949 gsl_matrix_scale (subs[0], -1.0);
2967 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2971 result = matrix_expr_evaluate_not (subs[0]);
2975 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2979 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2983 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2987 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2990 case MOP_PASTE_HORZ:
2991 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2994 case MOP_PASTE_VERT:
2995 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2999 result = gsl_matrix_alloc (0, 0);
3003 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
3007 result = matrix_expr_evaluate_vec_all (subs[0]);
3011 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
3015 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
3019 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
3028 for (size_t i = 0; i < e->n_subs; i++)
3029 if (subs[i] != result)
3030 gsl_matrix_free (subs[i]);
3031 if (subs != local_subs)
3037 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
3040 gsl_matrix *m = matrix_expr_evaluate (e);
3046 msg (SE, _("Expression for %s must evaluate to scalar, "
3047 "not a %zu×%zu matrix."),
3048 context, m->size1, m->size2);
3049 gsl_matrix_free (m);
3054 gsl_matrix_free (m);
3059 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
3063 if (!matrix_expr_evaluate_scalar (e, context, &d))
3067 if (d < LONG_MIN || d > LONG_MAX)
3069 msg (SE, _("Expression for %s is outside the integer range."), context);
3076 struct matrix_lvalue
3078 struct matrix_var *var;
3079 struct matrix_expr *indexes[2];
3084 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
3088 for (size_t i = 0; i < lvalue->n_indexes; i++)
3089 matrix_expr_destroy (lvalue->indexes[i]);
3094 static struct matrix_lvalue *
3095 matrix_lvalue_parse (struct matrix_state *s)
3097 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3098 if (!lex_force_id (s->lexer))
3100 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3101 if (lex_next_token (s->lexer, 1) == T_LPAREN)
3105 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3110 lex_get_n (s->lexer, 2);
3112 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3114 lvalue->n_indexes++;
3116 if (lex_match (s->lexer, T_COMMA))
3118 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3120 lvalue->n_indexes++;
3122 if (!lex_force_match (s->lexer, T_RPAREN))
3128 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3134 matrix_lvalue_destroy (lvalue);
3139 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3140 enum index_type index_type, size_t other_size,
3141 struct index_vector *iv)
3146 m = matrix_expr_evaluate (e);
3153 bool ok = matrix_normalize_index_vector (m, size, index_type,
3155 gsl_matrix_free (m);
3160 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3161 struct index_vector *iv, gsl_matrix *sm)
3163 gsl_vector dv = to_vector (lvalue->var->value);
3165 /* Convert source matrix 'sm' to source vector 'sv'. */
3166 if (!is_vector (sm))
3168 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
3169 sm->size1, sm->size2);
3172 gsl_vector sv = to_vector (sm);
3173 if (iv->n != sv.size)
3175 msg (SE, _("Can't assign %zu-element vector "
3176 "to %zu-element subvector."), sv.size, iv->n);
3180 /* Assign elements. */
3181 for (size_t x = 0; x < iv->n; x++)
3182 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3187 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3188 struct index_vector *iv0,
3189 struct index_vector *iv1, gsl_matrix *sm)
3191 gsl_matrix *dm = lvalue->var->value;
3193 /* Convert source matrix 'sm' to source vector 'sv'. */
3194 if (iv0->n != sm->size1)
3196 msg (SE, _("Row index vector for assignment to %s has %zu elements "
3197 "but source matrix has %zu rows."),
3198 lvalue->var->name, iv0->n, sm->size1);
3201 if (iv1->n != sm->size2)
3203 msg (SE, _("Column index vector for assignment to %s has %zu elements "
3204 "but source matrix has %zu columns."),
3205 lvalue->var->name, iv1->n, sm->size2);
3209 /* Assign elements. */
3210 for (size_t y = 0; y < iv0->n; y++)
3212 size_t f0 = iv0->indexes[y];
3213 for (size_t x = 0; x < iv1->n; x++)
3215 size_t f1 = iv1->indexes[x];
3216 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3222 /* Takes ownership of M. */
3224 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3225 struct index_vector *iv0, struct index_vector *iv1,
3228 if (!lvalue->n_indexes)
3230 gsl_matrix_free (lvalue->var->value);
3231 lvalue->var->value = sm;
3236 bool ok = (lvalue->n_indexes == 1
3237 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3238 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3239 free (iv0->indexes);
3240 free (iv1->indexes);
3241 gsl_matrix_free (sm);
3247 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3248 struct index_vector *iv0,
3249 struct index_vector *iv1)
3251 *iv0 = INDEX_VECTOR_INIT;
3252 *iv1 = INDEX_VECTOR_INIT;
3253 if (!lvalue->n_indexes)
3256 /* Validate destination matrix exists and has the right shape. */
3257 gsl_matrix *dm = lvalue->var->value;
3260 msg (SE, _("Undefined variable %s."), lvalue->var->name);
3263 else if (dm->size1 == 0 || dm->size2 == 0)
3265 msg (SE, _("Cannot index %zu×%zu matrix."),
3266 dm->size1, dm->size2);
3269 else if (lvalue->n_indexes == 1)
3271 if (!is_vector (dm))
3273 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
3274 dm->size1, dm->size2, lvalue->var->name);
3277 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3278 MAX (dm->size1, dm->size2),
3283 assert (lvalue->n_indexes == 2);
3284 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3285 IV_ROW, dm->size2, iv0))
3288 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3289 IV_COLUMN, dm->size1, iv1))
3291 free (iv0->indexes);
3298 /* Takes ownership of M. */
3300 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3302 struct index_vector iv0, iv1;
3303 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3305 gsl_matrix_free (sm);
3309 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3313 /* Matrix command. */
3317 struct matrix_cmd **commands;
3321 static void matrix_cmds_uninit (struct matrix_cmds *);
3325 enum matrix_cmd_type
3348 struct compute_command
3350 struct matrix_lvalue *lvalue;
3351 struct matrix_expr *rvalue;
3355 struct print_command
3357 struct matrix_expr *expression;
3358 bool use_default_format;
3359 struct fmt_spec format;
3361 int space; /* -1 means NEWPAGE. */
3365 struct string_array literals; /* CLABELS/RLABELS. */
3366 struct matrix_expr *expr; /* CNAMES/RNAMES. */
3372 struct do_if_command
3376 struct matrix_expr *condition;
3377 struct matrix_cmds commands;
3387 struct matrix_var *var;
3388 struct matrix_expr *start;
3389 struct matrix_expr *end;
3390 struct matrix_expr *increment;
3392 /* Loop conditions. */
3393 struct matrix_expr *top_condition;
3394 struct matrix_expr *bottom_condition;
3397 struct matrix_cmds commands;
3401 struct display_command
3403 struct matrix_state *state;
3407 struct release_command
3409 struct matrix_var **vars;
3416 struct matrix_expr *expression;
3417 struct save_file *sf;
3423 struct read_file *rf;
3424 struct matrix_lvalue *dst;
3425 struct matrix_expr *size;
3427 enum fmt_type format;
3434 struct write_command
3436 struct write_file *wf;
3437 struct matrix_expr *expression;
3440 /* If this is nonnull, WRITE uses this format.
3442 If this is NULL, WRITE uses free-field format with as many
3443 digits of precision as needed. */
3444 struct fmt_spec *format;
3453 struct matrix_lvalue *dst;
3454 struct dataset *dataset;
3455 struct file_handle *file;
3457 struct var_syntax *vars;
3459 struct matrix_var *names;
3461 /* Treatment of missing values. */
3466 MGET_ERROR, /* Flag error on command. */
3467 MGET_ACCEPT, /* Accept without change, user-missing only. */
3468 MGET_OMIT, /* Drop this case. */
3469 MGET_RECODE /* Recode to 'substitute'. */
3472 double substitute; /* MGET_RECODE only. */
3478 struct msave_command
3480 struct msave_common *common;
3481 struct matrix_expr *expr;
3482 const char *rowtype;
3483 const struct matrix_expr *factors;
3484 const struct matrix_expr *splits;
3490 struct matrix_state *state;
3491 struct file_handle *file;
3493 struct stringi_set rowtypes;
3497 struct eigen_command
3499 struct matrix_expr *expr;
3500 struct matrix_var *evec;
3501 struct matrix_var *eval;
3505 struct setdiag_command
3507 struct matrix_var *dst;
3508 struct matrix_expr *expr;
3514 struct matrix_expr *expr;
3515 struct matrix_var *u;
3516 struct matrix_var *s;
3517 struct matrix_var *v;
3523 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3524 static bool matrix_cmd_execute (struct matrix_cmd *);
3525 static void matrix_cmd_destroy (struct matrix_cmd *);
3528 static struct matrix_cmd *
3529 matrix_parse_compute (struct matrix_state *s)
3531 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3532 *cmd = (struct matrix_cmd) {
3533 .type = MCMD_COMPUTE,
3534 .compute = { .lvalue = NULL },
3537 struct compute_command *compute = &cmd->compute;
3538 compute->lvalue = matrix_lvalue_parse (s);
3539 if (!compute->lvalue)
3542 if (!lex_force_match (s->lexer, T_EQUALS))
3545 compute->rvalue = matrix_parse_expr (s);
3546 if (!compute->rvalue)
3552 matrix_cmd_destroy (cmd);
3557 matrix_cmd_execute_compute (struct compute_command *compute)
3559 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3563 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3567 print_labels_destroy (struct print_labels *labels)
3571 string_array_destroy (&labels->literals);
3572 matrix_expr_destroy (labels->expr);
3578 parse_literal_print_labels (struct matrix_state *s,
3579 struct print_labels **labelsp)
3581 lex_match (s->lexer, T_EQUALS);
3582 print_labels_destroy (*labelsp);
3583 *labelsp = xzalloc (sizeof **labelsp);
3584 while (lex_token (s->lexer) != T_SLASH
3585 && lex_token (s->lexer) != T_ENDCMD
3586 && lex_token (s->lexer) != T_STOP)
3588 struct string label = DS_EMPTY_INITIALIZER;
3589 while (lex_token (s->lexer) != T_COMMA
3590 && lex_token (s->lexer) != T_SLASH
3591 && lex_token (s->lexer) != T_ENDCMD
3592 && lex_token (s->lexer) != T_STOP)
3594 if (!ds_is_empty (&label))
3595 ds_put_byte (&label, ' ');
3597 if (lex_token (s->lexer) == T_STRING)
3598 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3601 char *rep = lex_next_representation (s->lexer, 0, 0);
3602 ds_put_cstr (&label, rep);
3607 string_array_append_nocopy (&(*labelsp)->literals,
3608 ds_steal_cstr (&label));
3610 if (!lex_match (s->lexer, T_COMMA))
3616 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3618 lex_match (s->lexer, T_EQUALS);
3619 struct matrix_expr *e = matrix_parse_exp (s);
3623 print_labels_destroy (*labelsp);
3624 *labelsp = xzalloc (sizeof **labelsp);
3625 (*labelsp)->expr = e;
3629 static struct matrix_cmd *
3630 matrix_parse_print (struct matrix_state *s)
3632 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3633 *cmd = (struct matrix_cmd) {
3636 .use_default_format = true,
3640 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3643 for (size_t i = 0; ; i++)
3645 enum token_type t = lex_next_token (s->lexer, i);
3646 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3648 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3650 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3653 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3658 cmd->print.expression = matrix_parse_exp (s);
3659 if (!cmd->print.expression)
3663 while (lex_match (s->lexer, T_SLASH))
3665 if (lex_match_id (s->lexer, "FORMAT"))
3667 lex_match (s->lexer, T_EQUALS);
3668 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3670 cmd->print.use_default_format = false;
3672 else if (lex_match_id (s->lexer, "TITLE"))
3674 lex_match (s->lexer, T_EQUALS);
3675 if (!lex_force_string (s->lexer))
3677 free (cmd->print.title);
3678 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3681 else if (lex_match_id (s->lexer, "SPACE"))
3683 lex_match (s->lexer, T_EQUALS);
3684 if (lex_match_id (s->lexer, "NEWPAGE"))
3685 cmd->print.space = -1;
3686 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3688 cmd->print.space = lex_integer (s->lexer);
3694 else if (lex_match_id (s->lexer, "RLABELS"))
3695 parse_literal_print_labels (s, &cmd->print.rlabels);
3696 else if (lex_match_id (s->lexer, "CLABELS"))
3697 parse_literal_print_labels (s, &cmd->print.clabels);
3698 else if (lex_match_id (s->lexer, "RNAMES"))
3700 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3703 else if (lex_match_id (s->lexer, "CNAMES"))
3705 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3710 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3711 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3719 matrix_cmd_destroy (cmd);
3724 matrix_is_integer (const gsl_matrix *m)
3726 for (size_t y = 0; y < m->size1; y++)
3727 for (size_t x = 0; x < m->size2; x++)
3729 double d = gsl_matrix_get (m, y, x);
3737 matrix_max_magnitude (const gsl_matrix *m)
3740 for (size_t y = 0; y < m->size1; y++)
3741 for (size_t x = 0; x < m->size2; x++)
3743 double d = fabs (gsl_matrix_get (m, y, x));
3751 format_fits (struct fmt_spec format, double x)
3753 char *s = data_out (&(union value) { .f = x }, NULL,
3754 &format, settings_get_fmt_settings ());
3755 bool fits = *s != '*' && !strchr (s, 'E');
3760 static struct fmt_spec
3761 default_format (const gsl_matrix *m, int *log_scale)
3765 double max = matrix_max_magnitude (m);
3767 if (matrix_is_integer (m))
3768 for (int w = 1; w <= 12; w++)
3770 struct fmt_spec format = { .type = FMT_F, .w = w };
3771 if (format_fits (format, -max))
3775 if (max >= 1e9 || max <= 1e-4)
3778 snprintf (s, sizeof s, "%.1e", max);
3780 const char *e = strchr (s, 'e');
3782 *log_scale = atoi (e + 1);
3785 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3789 trimmed_string (double d)
3791 struct substring s = ss_buffer ((char *) &d, sizeof d);
3792 ss_rtrim (&s, ss_cstr (" "));
3793 return ss_xstrdup (s);
3796 static struct string_array *
3797 print_labels_get (const struct print_labels *labels, size_t n,
3798 const char *prefix, bool truncate)
3803 struct string_array *out = xzalloc (sizeof *out);
3804 if (labels->literals.n)
3805 string_array_clone (out, &labels->literals);
3806 else if (labels->expr)
3808 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3809 if (m && is_vector (m))
3811 gsl_vector v = to_vector (m);
3812 for (size_t i = 0; i < v.size; i++)
3813 string_array_append_nocopy (out, trimmed_string (
3814 gsl_vector_get (&v, i)));
3816 gsl_matrix_free (m);
3822 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3825 string_array_append (out, "");
3828 string_array_delete (out, out->n - 1);
3831 for (size_t i = 0; i < out->n; i++)
3833 char *s = out->strings[i];
3834 s[strnlen (s, 8)] = '\0';
3841 matrix_cmd_print_space (int space)
3844 output_item_submit (page_break_item_create ());
3845 for (int i = 0; i < space; i++)
3846 output_log ("%s", "");
3850 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3851 struct fmt_spec format, int log_scale)
3853 matrix_cmd_print_space (print->space);
3855 output_log ("%s", print->title);
3857 output_log (" 10 ** %d X", log_scale);
3859 struct string_array *clabels = print_labels_get (print->clabels,
3860 m->size2, "col", true);
3861 if (clabels && format.w < 8)
3863 struct string_array *rlabels = print_labels_get (print->rlabels,
3864 m->size1, "row", true);
3868 struct string line = DS_EMPTY_INITIALIZER;
3870 ds_put_byte_multiple (&line, ' ', 8);
3871 for (size_t x = 0; x < m->size2; x++)
3872 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3873 output_log_nocopy (ds_steal_cstr (&line));
3876 double scale = pow (10.0, log_scale);
3877 bool numeric = fmt_is_numeric (format.type);
3878 for (size_t y = 0; y < m->size1; y++)
3880 struct string line = DS_EMPTY_INITIALIZER;
3882 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3884 for (size_t x = 0; x < m->size2; x++)
3886 double f = gsl_matrix_get (m, y, x);
3888 ? data_out (&(union value) { .f = f / scale}, NULL,
3889 &format, settings_get_fmt_settings ())
3890 : trimmed_string (f));
3891 ds_put_format (&line, " %s", s);
3894 output_log_nocopy (ds_steal_cstr (&line));
3897 string_array_destroy (rlabels);
3899 string_array_destroy (clabels);
3904 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3905 const struct print_labels *print_labels, size_t n,
3906 const char *name, const char *prefix)
3908 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3910 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3911 for (size_t i = 0; i < n; i++)
3912 pivot_category_create_leaf (
3914 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3915 : pivot_value_new_integer (i + 1)));
3917 d->hide_all_labels = true;
3918 string_array_destroy (labels);
3923 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3924 struct fmt_spec format, int log_scale)
3926 struct pivot_table *table = pivot_table_create__ (
3927 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3930 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3932 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3933 N_("Columns"), "col");
3935 struct pivot_footnote *footnote = NULL;
3938 char *s = xasprintf ("× 10**%d", log_scale);
3939 footnote = pivot_table_create_footnote (
3940 table, pivot_value_new_user_text_nocopy (s));
3943 double scale = pow (10.0, log_scale);
3944 bool numeric = fmt_is_numeric (format.type);
3945 for (size_t y = 0; y < m->size1; y++)
3946 for (size_t x = 0; x < m->size2; x++)
3948 double f = gsl_matrix_get (m, y, x);
3949 struct pivot_value *v;
3952 v = pivot_value_new_number (f / scale);
3953 v->numeric.format = format;
3956 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3958 pivot_value_add_footnote (v, footnote);
3959 pivot_table_put2 (table, y, x, v);
3962 pivot_table_submit (table);
3966 matrix_cmd_execute_print (const struct print_command *print)
3968 if (print->expression)
3970 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3975 struct fmt_spec format = (print->use_default_format
3976 ? default_format (m, &log_scale)
3979 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3980 matrix_cmd_print_text (print, m, format, log_scale);
3982 matrix_cmd_print_tables (print, m, format, log_scale);
3984 gsl_matrix_free (m);
3988 matrix_cmd_print_space (print->space);
3990 output_log ("%s", print->title);
3997 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3998 const char *command_name,
3999 const char *stop1, const char *stop2)
4001 lex_end_of_command (s->lexer);
4002 lex_discard_rest_of_command (s->lexer);
4004 size_t allocated = 0;
4007 while (lex_token (s->lexer) == T_ENDCMD)
4010 if (lex_at_phrase (s->lexer, stop1)
4011 || (stop2 && lex_at_phrase (s->lexer, stop2)))
4014 if (lex_at_phrase (s->lexer, "END MATRIX"))
4016 msg (SE, _("Premature END MATRIX within %s."), command_name);
4020 struct matrix_cmd *cmd = matrix_parse_command (s);
4024 if (c->n >= allocated)
4025 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
4026 c->commands[c->n++] = cmd;
4031 matrix_parse_do_if_clause (struct matrix_state *s,
4032 struct do_if_command *ifc,
4033 bool parse_condition,
4034 size_t *allocated_clauses)
4036 if (ifc->n_clauses >= *allocated_clauses)
4037 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
4038 sizeof *ifc->clauses);
4039 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
4040 *c = (struct do_if_clause) { .condition = NULL };
4042 if (parse_condition)
4044 c->condition = matrix_parse_expr (s);
4049 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
4052 static struct matrix_cmd *
4053 matrix_parse_do_if (struct matrix_state *s)
4055 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4056 *cmd = (struct matrix_cmd) {
4058 .do_if = { .n_clauses = 0 }
4061 size_t allocated_clauses = 0;
4064 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
4067 while (lex_match_phrase (s->lexer, "ELSE IF"));
4069 if (lex_match_id (s->lexer, "ELSE")
4070 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
4073 if (!lex_match_phrase (s->lexer, "END IF"))
4078 matrix_cmd_destroy (cmd);
4083 matrix_cmd_execute_do_if (struct do_if_command *cmd)
4085 for (size_t i = 0; i < cmd->n_clauses; i++)
4087 struct do_if_clause *c = &cmd->clauses[i];
4091 if (!matrix_expr_evaluate_scalar (c->condition,
4092 i ? "ELSE IF" : "DO IF",
4097 for (size_t j = 0; j < c->commands.n; j++)
4098 if (!matrix_cmd_execute (c->commands.commands[j]))
4105 static struct matrix_cmd *
4106 matrix_parse_loop (struct matrix_state *s)
4108 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4109 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4111 struct loop_command *loop = &cmd->loop;
4112 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4114 struct substring name = lex_tokss (s->lexer);
4115 loop->var = matrix_var_lookup (s, name);
4117 loop->var = matrix_var_create (s, name);
4122 loop->start = matrix_parse_expr (s);
4123 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4126 loop->end = matrix_parse_expr (s);
4130 if (lex_match (s->lexer, T_BY))
4132 loop->increment = matrix_parse_expr (s);
4133 if (!loop->increment)
4138 if (lex_match_id (s->lexer, "IF"))
4140 loop->top_condition = matrix_parse_expr (s);
4141 if (!loop->top_condition)
4145 bool was_in_loop = s->in_loop;
4147 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4149 s->in_loop = was_in_loop;
4153 if (!lex_match_phrase (s->lexer, "END LOOP"))
4156 if (lex_match_id (s->lexer, "IF"))
4158 loop->bottom_condition = matrix_parse_expr (s);
4159 if (!loop->bottom_condition)
4166 matrix_cmd_destroy (cmd);
4171 matrix_cmd_execute_loop (struct loop_command *cmd)
4175 long int increment = 1;
4178 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4179 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4181 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4185 if (increment > 0 ? value > end
4186 : increment < 0 ? value < end
4191 int mxloops = settings_get_mxloops ();
4192 for (int i = 0; i < mxloops; i++)
4197 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4199 gsl_matrix_free (cmd->var->value);
4200 cmd->var->value = NULL;
4202 if (!cmd->var->value)
4203 cmd->var->value = gsl_matrix_alloc (1, 1);
4205 gsl_matrix_set (cmd->var->value, 0, 0, value);
4208 if (cmd->top_condition)
4211 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4216 for (size_t j = 0; j < cmd->commands.n; j++)
4217 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4220 if (cmd->bottom_condition)
4223 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4224 "END LOOP IF", &d) || d > 0)
4231 if (increment > 0 ? value > end : value < end)
4237 static struct matrix_cmd *
4238 matrix_parse_break (struct matrix_state *s)
4242 msg (SE, _("BREAK not inside LOOP."));
4246 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4247 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4251 static struct matrix_cmd *
4252 matrix_parse_display (struct matrix_state *s)
4254 while (lex_token (s->lexer) != T_ENDCMD)
4256 if (!lex_match_id (s->lexer, "DICTIONARY")
4257 && !lex_match_id (s->lexer, "STATUS"))
4259 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4264 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4265 *cmd = (struct matrix_cmd) { .type = MCMD_DISPLAY, .display = { s } };
4270 compare_matrix_var_pointers (const void *a_, const void *b_)
4272 const struct matrix_var *const *ap = a_;
4273 const struct matrix_var *const *bp = b_;
4274 const struct matrix_var *a = *ap;
4275 const struct matrix_var *b = *bp;
4276 return strcmp (a->name, b->name);
4280 matrix_cmd_execute_display (struct display_command *cmd)
4282 const struct matrix_state *s = cmd->state;
4284 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4285 struct pivot_dimension *attr_dimension
4286 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
4287 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
4288 N_("Rows"), N_("Columns"));
4289 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
4291 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4294 const struct matrix_var *var;
4295 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4297 vars[n_vars++] = var;
4298 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4300 struct pivot_dimension *rows = pivot_dimension_create (
4301 table, PIVOT_AXIS_ROW, N_("Variable"));
4302 for (size_t i = 0; i < n_vars; i++)
4304 const struct matrix_var *var = vars[i];
4305 pivot_category_create_leaf (
4306 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4308 size_t r = var->value->size1;
4309 size_t c = var->value->size2;
4310 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4311 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4312 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4315 pivot_table_submit (table);
4318 static struct matrix_cmd *
4319 matrix_parse_release (struct matrix_state *s)
4321 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4322 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4324 size_t allocated_vars = 0;
4325 while (lex_token (s->lexer) == T_ID)
4327 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4330 if (cmd->release.n_vars >= allocated_vars)
4331 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4332 sizeof *cmd->release.vars);
4333 cmd->release.vars[cmd->release.n_vars++] = var;
4336 lex_error (s->lexer, _("Variable name expected."));
4339 if (!lex_match (s->lexer, T_COMMA))
4347 matrix_cmd_execute_release (struct release_command *cmd)
4349 for (size_t i = 0; i < cmd->n_vars; i++)
4351 struct matrix_var *var = cmd->vars[i];
4352 gsl_matrix_free (var->value);
4357 static struct save_file *
4358 save_file_create (struct matrix_state *s, struct file_handle *fh,
4359 struct string_array *variables,
4360 struct matrix_expr *names,
4361 struct stringi_set *strings)
4363 for (size_t i = 0; i < s->n_save_files; i++)
4365 struct save_file *sf = s->save_files[i];
4366 if (fh_equal (sf->file, fh))
4370 string_array_destroy (variables);
4371 matrix_expr_destroy (names);
4372 stringi_set_destroy (strings);
4378 struct save_file *sf = xmalloc (sizeof *sf);
4379 *sf = (struct save_file) {
4381 .dataset = s->dataset,
4382 .variables = *variables,
4384 .strings = STRINGI_SET_INITIALIZER (sf->strings),
4387 stringi_set_swap (&sf->strings, strings);
4388 stringi_set_destroy (strings);
4390 s->save_files = xrealloc (s->save_files,
4391 (s->n_save_files + 1) * sizeof *s->save_files);
4392 s->save_files[s->n_save_files++] = sf;
4396 static struct casewriter *
4397 save_file_open (struct save_file *sf, gsl_matrix *m)
4399 if (sf->writer || sf->error)
4403 size_t n_variables = caseproto_get_n_widths (
4404 casewriter_get_proto (sf->writer));
4405 if (m->size2 != n_variables)
4407 msg (ME, _("The first SAVE to %s within this matrix program "
4408 "had %zu columns, so a %zu×%zu matrix cannot be "
4410 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
4411 n_variables, m->size1, m->size2);
4419 struct dictionary *dict = dict_create (get_default_encoding ());
4421 /* Fill 'names' with user-specified names if there were any, either from
4422 sf->variables or sf->names. */
4423 struct string_array names = { .n = 0 };
4424 if (sf->variables.n)
4425 string_array_clone (&names, &sf->variables);
4428 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
4429 if (nm && is_vector (nm))
4431 gsl_vector nv = to_vector (nm);
4432 for (size_t i = 0; i < nv.size; i++)
4434 char *name = trimmed_string (gsl_vector_get (&nv, i));
4435 if (dict_id_is_valid (dict, name, true))
4436 string_array_append_nocopy (&names, name);
4441 gsl_matrix_free (nm);
4444 struct stringi_set strings;
4445 stringi_set_clone (&strings, &sf->strings);
4447 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4452 name = names.strings[i];
4455 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4459 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4460 struct variable *var = dict_create_var (dict, name, width);
4463 msg (ME, _("Duplicate variable name %s in SAVE statement."),
4469 if (!stringi_set_is_empty (&strings))
4471 size_t count = stringi_set_count (&strings);
4472 const char *example = stringi_set_node_get_string (
4473 stringi_set_first (&strings));
4476 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
4477 "variable %s."), example);
4479 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
4480 "unknown variable: %s.",
4481 "The SAVE command STRINGS subcommand specifies %zu "
4482 "unknown variables, including %s.",
4487 stringi_set_destroy (&strings);
4488 string_array_destroy (&names);
4497 if (sf->file == fh_inline_file ())
4498 sf->writer = autopaging_writer_create (dict_get_proto (dict));
4500 sf->writer = any_writer_open (sf->file, dict);
4512 save_file_destroy (struct save_file *sf)
4516 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
4518 dataset_set_dict (sf->dataset, sf->dict);
4519 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
4523 casewriter_destroy (sf->writer);
4524 dict_unref (sf->dict);
4526 fh_unref (sf->file);
4527 string_array_destroy (&sf->variables);
4528 matrix_expr_destroy (sf->names);
4529 stringi_set_destroy (&sf->strings);
4534 static struct matrix_cmd *
4535 matrix_parse_save (struct matrix_state *s)
4537 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4538 *cmd = (struct matrix_cmd) {
4540 .save = { .expression = NULL },
4543 struct file_handle *fh = NULL;
4544 struct save_command *save = &cmd->save;
4546 struct string_array variables = STRING_ARRAY_INITIALIZER;
4547 struct matrix_expr *names = NULL;
4548 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
4550 save->expression = matrix_parse_exp (s);
4551 if (!save->expression)
4554 while (lex_match (s->lexer, T_SLASH))
4556 if (lex_match_id (s->lexer, "OUTFILE"))
4558 lex_match (s->lexer, T_EQUALS);
4561 fh = (lex_match (s->lexer, T_ASTERISK)
4563 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4567 else if (lex_match_id (s->lexer, "VARIABLES"))
4569 lex_match (s->lexer, T_EQUALS);
4573 struct dictionary *d = dict_create (get_default_encoding ());
4574 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4575 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4580 string_array_clear (&variables);
4581 variables = (struct string_array) {
4587 else if (lex_match_id (s->lexer, "NAMES"))
4589 lex_match (s->lexer, T_EQUALS);
4590 matrix_expr_destroy (names);
4591 names = matrix_parse_exp (s);
4595 else if (lex_match_id (s->lexer, "STRINGS"))
4597 lex_match (s->lexer, T_EQUALS);
4598 while (lex_token (s->lexer) == T_ID)
4600 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
4602 if (!lex_match (s->lexer, T_COMMA))
4608 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4616 if (s->prev_save_file)
4617 fh = fh_ref (s->prev_save_file);
4620 lex_sbc_missing ("OUTFILE");
4624 fh_unref (s->prev_save_file);
4625 s->prev_save_file = fh_ref (fh);
4627 if (variables.n && names)
4629 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4630 matrix_expr_destroy (names);
4634 save->sf = save_file_create (s, fh, &variables, names, &strings);
4638 string_array_destroy (&variables);
4639 matrix_expr_destroy (names);
4640 stringi_set_destroy (&strings);
4642 matrix_cmd_destroy (cmd);
4647 matrix_cmd_execute_save (const struct save_command *save)
4649 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4653 struct casewriter *writer = save_file_open (save->sf, m);
4656 gsl_matrix_free (m);
4660 const struct caseproto *proto = casewriter_get_proto (writer);
4661 for (size_t y = 0; y < m->size1; y++)
4663 struct ccase *c = case_create (proto);
4664 for (size_t x = 0; x < m->size2; x++)
4666 double d = gsl_matrix_get (m, y, x);
4667 int width = caseproto_get_width (proto, x);
4668 union value *value = case_data_rw_idx (c, x);
4672 memcpy (value->s, &d, width);
4674 casewriter_write (writer, c);
4676 gsl_matrix_free (m);
4679 static struct matrix_cmd *
4680 matrix_parse_read (struct matrix_state *s)
4682 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4683 *cmd = (struct matrix_cmd) {
4685 .read = { .format = FMT_F },
4688 struct file_handle *fh = NULL;
4689 char *encoding = NULL;
4690 struct read_command *read = &cmd->read;
4691 read->dst = matrix_lvalue_parse (s);
4696 int repetitions = 0;
4697 int record_width = 0;
4698 bool seen_format = false;
4699 while (lex_match (s->lexer, T_SLASH))
4701 if (lex_match_id (s->lexer, "FILE"))
4703 lex_match (s->lexer, T_EQUALS);
4706 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4710 else if (lex_match_id (s->lexer, "ENCODING"))
4712 lex_match (s->lexer, T_EQUALS);
4713 if (!lex_force_string (s->lexer))
4717 encoding = ss_xstrdup (lex_tokss (s->lexer));
4721 else if (lex_match_id (s->lexer, "FIELD"))
4723 lex_match (s->lexer, T_EQUALS);
4725 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4727 read->c1 = lex_integer (s->lexer);
4729 if (!lex_force_match (s->lexer, T_TO)
4730 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4732 read->c2 = lex_integer (s->lexer) + 1;
4735 record_width = read->c2 - read->c1;
4736 if (lex_match (s->lexer, T_BY))
4738 if (!lex_force_int_range (s->lexer, "BY", 1,
4739 read->c2 - read->c1))
4741 by = lex_integer (s->lexer);
4744 if (record_width % by)
4746 msg (SE, _("BY %d does not evenly divide record width %d."),
4754 else if (lex_match_id (s->lexer, "SIZE"))
4756 lex_match (s->lexer, T_EQUALS);
4757 matrix_expr_destroy (read->size);
4758 read->size = matrix_parse_exp (s);
4762 else if (lex_match_id (s->lexer, "MODE"))
4764 lex_match (s->lexer, T_EQUALS);
4765 if (lex_match_id (s->lexer, "RECTANGULAR"))
4766 read->symmetric = false;
4767 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4768 read->symmetric = true;
4771 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4775 else if (lex_match_id (s->lexer, "REREAD"))
4776 read->reread = true;
4777 else if (lex_match_id (s->lexer, "FORMAT"))
4781 lex_sbc_only_once ("FORMAT");
4786 lex_match (s->lexer, T_EQUALS);
4788 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4791 const char *p = lex_tokcstr (s->lexer);
4792 if (c_isdigit (p[0]))
4794 repetitions = atoi (p);
4795 p += strspn (p, "0123456789");
4796 if (!fmt_from_name (p, &read->format))
4798 lex_error (s->lexer, _("Unknown format %s."), p);
4803 else if (fmt_from_name (p, &read->format))
4807 struct fmt_spec format;
4808 if (!parse_format_specifier (s->lexer, &format))
4810 read->format = format.type;
4816 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4817 "REREAD", "FORMAT");
4824 lex_sbc_missing ("FIELD");
4828 if (!read->dst->n_indexes && !read->size)
4830 msg (SE, _("SIZE is required for reading data into a full matrix "
4831 "(as opposed to a submatrix)."));
4837 if (s->prev_read_file)
4838 fh = fh_ref (s->prev_read_file);
4841 lex_sbc_missing ("FILE");
4845 fh_unref (s->prev_read_file);
4846 s->prev_read_file = fh_ref (fh);
4848 read->rf = read_file_create (s, fh);
4852 free (read->rf->encoding);
4853 read->rf->encoding = encoding;
4857 /* Field width may be specified in multiple ways:
4860 2. The format on FORMAT.
4861 3. The repetition factor on FORMAT.
4863 (2) and (3) are mutually exclusive.
4865 If more than one of these is present, they must agree. If none of them is
4866 present, then free-field format is used.
4868 if (repetitions > record_width)
4870 msg (SE, _("%d repetitions cannot fit in record width %d."),
4871 repetitions, record_width);
4874 int w = (repetitions ? record_width / repetitions
4880 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4881 "which implies field width %d, "
4882 "but BY specifies field width %d."),
4883 repetitions, record_width, w, by);
4885 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4894 matrix_cmd_destroy (cmd);
4900 parse_error (const struct dfm_reader *reader, enum fmt_type format,
4901 struct substring data, size_t y, size_t x,
4902 int first_column, int last_column, char *error)
4904 int line_number = dfm_get_line_number (reader);
4905 struct msg_location *location = xmalloc (sizeof *location);
4906 *location = (struct msg_location) {
4907 .file_name = xstrdup (dfm_get_file_name (reader)),
4908 .first_line = line_number,
4909 .last_line = line_number + 1,
4910 .first_column = first_column,
4911 .last_column = last_column,
4913 struct msg *m = xmalloc (sizeof *m);
4915 .category = MSG_C_DATA,
4916 .severity = MSG_S_WARNING,
4917 .location = location,
4918 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
4919 "for matrix row %zu, column %zu: %s"),
4920 (int) data.length, data.string, fmt_name (format),
4921 y + 1, x + 1, error),
4929 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4930 gsl_matrix *m, struct substring p, size_t y, size_t x,
4931 const char *line_start)
4933 const char *input_encoding = dfm_reader_get_encoding (reader);
4936 if (fmt_is_numeric (read->format))
4939 error = data_in (p, input_encoding, read->format,
4940 settings_get_fmt_settings (), &v, 0, NULL);
4941 if (!error && v.f == SYSMIS)
4942 error = xstrdup (_("Matrix data may not contain missing value."));
4947 uint8_t s[sizeof (double)];
4948 union value v = { .s = s };
4949 error = data_in (p, input_encoding, read->format,
4950 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
4951 memcpy (&f, s, sizeof f);
4956 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
4957 int c2 = c1 + ss_utf8_count_columns (p) - 1;
4958 parse_error (reader, read->format, p, y, x, c1, c2, error);
4962 gsl_matrix_set (m, y, x, f);
4963 if (read->symmetric && x != y)
4964 gsl_matrix_set (m, x, y, f);
4969 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4970 struct substring *line, const char **startp)
4972 if (dfm_eof (reader))
4974 msg (SE, _("Unexpected end of file reading matrix data."));
4977 dfm_expand_tabs (reader);
4978 struct substring record = dfm_get_record (reader);
4979 /* XXX need to recode record into UTF-8 */
4980 *startp = record.string;
4981 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
4986 matrix_read (struct read_command *read, struct dfm_reader *reader,
4989 for (size_t y = 0; y < m->size1; y++)
4991 size_t nx = read->symmetric ? y + 1 : m->size2;
4993 struct substring line = ss_empty ();
4994 const char *line_start = line.string;
4995 for (size_t x = 0; x < nx; x++)
5002 ss_ltrim (&line, ss_cstr (" ,"));
5003 if (!ss_is_empty (line))
5005 if (!matrix_read_line (read, reader, &line, &line_start))
5007 dfm_forward_record (reader);
5010 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
5014 if (!matrix_read_line (read, reader, &line, &line_start))
5016 size_t fields_per_line = (read->c2 - read->c1) / read->w;
5017 int f = x % fields_per_line;
5018 if (f == fields_per_line - 1)
5019 dfm_forward_record (reader);
5021 p = ss_substr (line, read->w * f, read->w);
5024 matrix_read_set_field (read, reader, m, p, y, x, line_start);
5028 dfm_forward_record (reader);
5031 ss_ltrim (&line, ss_cstr (" ,"));
5032 if (!ss_is_empty (line))
5035 msg (SW, _("Trailing garbage on line \"%.*s\""),
5036 (int) line.length, line.string);
5043 matrix_cmd_execute_read (struct read_command *read)
5045 struct index_vector iv0, iv1;
5046 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
5049 size_t size[2] = { SIZE_MAX, SIZE_MAX };
5052 gsl_matrix *m = matrix_expr_evaluate (read->size);
5058 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5059 "not a %zu×%zu matrix."), m->size1, m->size2);
5060 gsl_matrix_free (m);
5066 gsl_vector v = to_vector (m);
5070 d[0] = gsl_vector_get (&v, 0);
5073 else if (v.size == 2)
5075 d[0] = gsl_vector_get (&v, 0);
5076 d[1] = gsl_vector_get (&v, 1);
5080 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5081 "not a %zu×%zu matrix."),
5082 m->size1, m->size2),
5083 gsl_matrix_free (m);
5088 gsl_matrix_free (m);
5090 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
5092 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
5093 "are outside valid range."),
5104 if (read->dst->n_indexes)
5106 size_t submatrix_size[2];
5107 if (read->dst->n_indexes == 2)
5109 submatrix_size[0] = iv0.n;
5110 submatrix_size[1] = iv1.n;
5112 else if (read->dst->var->value->size1 == 1)
5114 submatrix_size[0] = 1;
5115 submatrix_size[1] = iv0.n;
5119 submatrix_size[0] = iv0.n;
5120 submatrix_size[1] = 1;
5125 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
5127 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
5128 "differ from submatrix dimensions %zu×%zu."),
5130 submatrix_size[0], submatrix_size[1]);
5138 size[0] = submatrix_size[0];
5139 size[1] = submatrix_size[1];
5143 struct dfm_reader *reader = read_file_open (read->rf);
5145 dfm_reread_record (reader, 1);
5147 if (read->symmetric && size[0] != size[1])
5149 msg (SE, _("Cannot read non-square %zu×%zu matrix "
5150 "using READ with MODE=SYMMETRIC."),
5156 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
5157 matrix_read (read, reader, tmp);
5158 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
5161 static struct matrix_cmd *
5162 matrix_parse_write (struct matrix_state *s)
5164 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5165 *cmd = (struct matrix_cmd) {
5169 struct file_handle *fh = NULL;
5170 char *encoding = NULL;
5171 struct write_command *write = &cmd->write;
5172 write->expression = matrix_parse_exp (s);
5173 if (!write->expression)
5177 int repetitions = 0;
5178 int record_width = 0;
5179 enum fmt_type format = FMT_F;
5180 bool has_format = false;
5181 while (lex_match (s->lexer, T_SLASH))
5183 if (lex_match_id (s->lexer, "OUTFILE"))
5185 lex_match (s->lexer, T_EQUALS);
5188 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5192 else if (lex_match_id (s->lexer, "ENCODING"))
5194 lex_match (s->lexer, T_EQUALS);
5195 if (!lex_force_string (s->lexer))
5199 encoding = ss_xstrdup (lex_tokss (s->lexer));
5203 else if (lex_match_id (s->lexer, "FIELD"))
5205 lex_match (s->lexer, T_EQUALS);
5207 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5209 write->c1 = lex_integer (s->lexer);
5211 if (!lex_force_match (s->lexer, T_TO)
5212 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5214 write->c2 = lex_integer (s->lexer) + 1;
5217 record_width = write->c2 - write->c1;
5218 if (lex_match (s->lexer, T_BY))
5220 if (!lex_force_int_range (s->lexer, "BY", 1,
5221 write->c2 - write->c1))
5223 by = lex_integer (s->lexer);
5226 if (record_width % by)
5228 msg (SE, _("BY %d does not evenly divide record width %d."),
5236 else if (lex_match_id (s->lexer, "MODE"))
5238 lex_match (s->lexer, T_EQUALS);
5239 if (lex_match_id (s->lexer, "RECTANGULAR"))
5240 write->triangular = false;
5241 else if (lex_match_id (s->lexer, "TRIANGULAR"))
5242 write->triangular = true;
5245 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5249 else if (lex_match_id (s->lexer, "HOLD"))
5251 else if (lex_match_id (s->lexer, "FORMAT"))
5253 if (has_format || write->format)
5255 lex_sbc_only_once ("FORMAT");
5259 lex_match (s->lexer, T_EQUALS);
5261 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5264 const char *p = lex_tokcstr (s->lexer);
5265 if (c_isdigit (p[0]))
5267 repetitions = atoi (p);
5268 p += strspn (p, "0123456789");
5269 if (!fmt_from_name (p, &format))
5271 lex_error (s->lexer, _("Unknown format %s."), p);
5277 else if (fmt_from_name (p, &format))
5284 struct fmt_spec spec;
5285 if (!parse_format_specifier (s->lexer, &spec))
5287 write->format = xmemdup (&spec, sizeof spec);
5292 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5300 lex_sbc_missing ("FIELD");
5306 if (s->prev_write_file)
5307 fh = fh_ref (s->prev_write_file);
5310 lex_sbc_missing ("OUTFILE");
5314 fh_unref (s->prev_write_file);
5315 s->prev_write_file = fh_ref (fh);
5317 write->wf = write_file_create (s, fh);
5321 free (write->wf->encoding);
5322 write->wf->encoding = encoding;
5326 /* Field width may be specified in multiple ways:
5329 2. The format on FORMAT.
5330 3. The repetition factor on FORMAT.
5332 (2) and (3) are mutually exclusive.
5334 If more than one of these is present, they must agree. If none of them is
5335 present, then free-field format is used.
5337 if (repetitions > record_width)
5339 msg (SE, _("%d repetitions cannot fit in record width %d."),
5340 repetitions, record_width);
5343 int w = (repetitions ? record_width / repetitions
5344 : write->format ? write->format->w
5349 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5350 "which implies field width %d, "
5351 "but BY specifies field width %d."),
5352 repetitions, record_width, w, by);
5354 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5358 if (w && !write->format)
5360 write->format = xmalloc (sizeof *write->format);
5361 *write->format = (struct fmt_spec) { .type = format, .w = w };
5363 if (!fmt_check_output (write->format))
5367 if (write->format && fmt_var_width (write->format) > sizeof (double))
5369 char s[FMT_STRING_LEN_MAX + 1];
5370 fmt_to_string (write->format, s);
5371 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
5372 s, sizeof (double));
5380 matrix_cmd_destroy (cmd);
5385 matrix_cmd_execute_write (struct write_command *write)
5387 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5391 if (write->triangular && m->size1 != m->size2)
5393 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5394 "the matrix to be written has dimensions %zu×%zu."),
5395 m->size1, m->size2);
5396 gsl_matrix_free (m);
5400 struct dfm_writer *writer = write_file_open (write->wf);
5401 if (!writer || !m->size1)
5403 gsl_matrix_free (m);
5407 const struct fmt_settings *settings = settings_get_fmt_settings ();
5408 struct u8_line *line = write->wf->held;
5409 for (size_t y = 0; y < m->size1; y++)
5413 line = xmalloc (sizeof *line);
5414 u8_line_init (line);
5416 size_t nx = write->triangular ? y + 1 : m->size2;
5418 for (size_t x = 0; x < nx; x++)
5421 double f = gsl_matrix_get (m, y, x);
5425 if (fmt_is_numeric (write->format->type))
5428 v.s = (uint8_t *) &f;
5429 s = data_out (&v, NULL, write->format, settings);
5433 s = xmalloc (DBL_BUFSIZE_BOUND);
5434 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
5435 >= DBL_BUFSIZE_BOUND)
5438 size_t len = strlen (s);
5439 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5440 if (width + x0 > write->c2)
5442 dfm_put_record_utf8 (writer, line->s.ss.string,
5444 u8_line_clear (line);
5447 u8_line_put (line, x0, x0 + width, s, len);
5450 x0 += write->format ? write->format->w : width + 1;
5453 if (y + 1 >= m->size1 && write->hold)
5455 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5456 u8_line_clear (line);
5460 u8_line_destroy (line);
5464 write->wf->held = line;
5466 gsl_matrix_free (m);
5469 static struct matrix_cmd *
5470 matrix_parse_get (struct matrix_state *s)
5472 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5473 *cmd = (struct matrix_cmd) {
5476 .dataset = s->dataset,
5477 .user = { .treatment = MGET_ERROR },
5478 .system = { .treatment = MGET_ERROR },
5482 struct get_command *get = &cmd->get;
5483 get->dst = matrix_lvalue_parse (s);
5487 while (lex_match (s->lexer, T_SLASH))
5489 if (lex_match_id (s->lexer, "FILE"))
5491 lex_match (s->lexer, T_EQUALS);
5493 fh_unref (get->file);
5494 if (lex_match (s->lexer, T_ASTERISK))
5498 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5503 else if (lex_match_id (s->lexer, "ENCODING"))
5505 lex_match (s->lexer, T_EQUALS);
5506 if (!lex_force_string (s->lexer))
5509 free (get->encoding);
5510 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5514 else if (lex_match_id (s->lexer, "VARIABLES"))
5516 lex_match (s->lexer, T_EQUALS);
5520 lex_sbc_only_once ("VARIABLES");
5524 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
5527 else if (lex_match_id (s->lexer, "NAMES"))
5529 lex_match (s->lexer, T_EQUALS);
5530 if (!lex_force_id (s->lexer))
5533 struct substring name = lex_tokss (s->lexer);
5534 get->names = matrix_var_lookup (s, name);
5536 get->names = matrix_var_create (s, name);
5539 else if (lex_match_id (s->lexer, "MISSING"))
5541 lex_match (s->lexer, T_EQUALS);
5542 if (lex_match_id (s->lexer, "ACCEPT"))
5543 get->user.treatment = MGET_ACCEPT;
5544 else if (lex_match_id (s->lexer, "OMIT"))
5545 get->user.treatment = MGET_OMIT;
5546 else if (lex_is_number (s->lexer))
5548 get->user.treatment = MGET_RECODE;
5549 get->user.substitute = lex_number (s->lexer);
5554 lex_error (s->lexer, NULL);
5558 else if (lex_match_id (s->lexer, "SYSMIS"))
5560 lex_match (s->lexer, T_EQUALS);
5561 if (lex_match_id (s->lexer, "OMIT"))
5562 get->system.treatment = MGET_OMIT;
5563 else if (lex_is_number (s->lexer))
5565 get->system.treatment = MGET_RECODE;
5566 get->system.substitute = lex_number (s->lexer);
5571 lex_error (s->lexer, NULL);
5577 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5578 "MISSING", "SYSMIS");
5583 if (get->user.treatment != MGET_ACCEPT)
5584 get->system.treatment = MGET_ERROR;
5589 matrix_cmd_destroy (cmd);
5594 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
5595 const struct dictionary *dict)
5597 struct variable **vars;
5602 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
5603 &vars, &n_vars, PV_NUMERIC))
5608 n_vars = dict_get_var_cnt (dict);
5609 vars = xnmalloc (n_vars, sizeof *vars);
5610 for (size_t i = 0; i < n_vars; i++)
5612 struct variable *var = dict_get_var (dict, i);
5613 if (!var_is_numeric (var))
5615 msg (SE, _("GET: Variable %s is not numeric."),
5616 var_get_name (var));
5626 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
5627 for (size_t i = 0; i < n_vars; i++)
5629 char s[sizeof (double)];
5631 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
5632 memcpy (&f, s, sizeof f);
5633 gsl_matrix_set (names, i, 0, f);
5636 gsl_matrix_free (get->names->value);
5637 get->names->value = names;
5641 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5642 long long int casenum = 1;
5644 for (struct ccase *c = casereader_read (reader); c;
5645 c = casereader_read (reader), casenum++)
5647 if (n_rows >= m->size1)
5649 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5650 for (size_t y = 0; y < n_rows; y++)
5651 for (size_t x = 0; x < n_vars; x++)
5652 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5653 gsl_matrix_free (m);
5658 for (size_t x = 0; x < n_vars; x++)
5660 const struct variable *var = vars[x];
5661 double d = case_num (c, var);
5664 if (get->system.treatment == MGET_RECODE)
5665 d = get->system.substitute;
5666 else if (get->system.treatment == MGET_OMIT)
5670 msg (SE, _("GET: Variable %s in case %lld "
5671 "is system-missing."),
5672 var_get_name (var), casenum);
5676 else if (var_is_num_missing (var, d, MV_USER))
5678 if (get->user.treatment == MGET_RECODE)
5679 d = get->user.substitute;
5680 else if (get->user.treatment == MGET_OMIT)
5682 else if (get->user.treatment != MGET_ACCEPT)
5684 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5686 var_get_name (var), casenum, d);
5690 gsl_matrix_set (m, n_rows, x, d);
5701 matrix_lvalue_evaluate_and_assign (get->dst, m);
5704 gsl_matrix_free (m);
5709 matrix_open_casereader (const char *command_name,
5710 struct file_handle *file, const char *encoding,
5711 struct dataset *dataset,
5712 struct casereader **readerp, struct dictionary **dictp)
5716 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
5717 return *readerp != NULL;
5721 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
5723 msg (ME, _("The %s command cannot read an empty active file."),
5727 *readerp = proc_open (dataset);
5728 *dictp = dict_ref (dataset_dict (dataset));
5734 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
5735 struct casereader *reader, struct dictionary *dict)
5738 casereader_destroy (reader);
5740 proc_commit (dataset);
5744 matrix_cmd_execute_get (struct get_command *get)
5746 struct casereader *r;
5747 struct dictionary *d;
5748 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
5751 matrix_cmd_execute_get__ (get, r, d);
5752 matrix_close_casereader (get->file, get->dataset, r, d);
5757 match_rowtype (struct lexer *lexer)
5759 static const char *rowtypes[] = {
5760 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5762 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5764 for (size_t i = 0; i < n_rowtypes; i++)
5765 if (lex_match_id (lexer, rowtypes[i]))
5768 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5773 parse_var_names (struct lexer *lexer, struct string_array *sa)
5775 lex_match (lexer, T_EQUALS);
5777 string_array_clear (sa);
5779 struct dictionary *dict = dict_create (get_default_encoding ());
5782 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5783 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5788 for (size_t i = 0; i < n_names; i++)
5789 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
5790 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
5792 msg (SE, _("Variable name %s is reserved."), names[i]);
5793 for (size_t j = 0; j < n_names; j++)
5799 string_array_clear (sa);
5800 sa->strings = names;
5801 sa->n = sa->allocated = n_names;
5807 msave_common_destroy (struct msave_common *common)
5811 fh_unref (common->outfile);
5812 string_array_destroy (&common->variables);
5813 string_array_destroy (&common->fnames);
5814 string_array_destroy (&common->snames);
5816 for (size_t i = 0; i < common->n_factors; i++)
5817 matrix_expr_destroy (common->factors[i]);
5818 free (common->factors);
5820 for (size_t i = 0; i < common->n_splits; i++)
5821 matrix_expr_destroy (common->splits[i]);
5822 free (common->splits);
5824 dict_unref (common->dict);
5825 casewriter_destroy (common->writer);
5832 compare_variables (const char *keyword,
5833 const struct string_array *new,
5834 const struct string_array *old)
5840 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5841 "on the first MSAVE within MATRIX."), keyword);
5844 else if (!string_array_equal_case (old, new))
5846 msg (SE, _("%s must specify the same variables each time within "
5847 "a given MATRIX."), keyword);
5854 static struct matrix_cmd *
5855 matrix_parse_msave (struct matrix_state *s)
5857 struct msave_common *common = xmalloc (sizeof *common);
5858 *common = (struct msave_common) { .outfile = NULL };
5860 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5861 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5863 struct matrix_expr *splits = NULL;
5864 struct matrix_expr *factors = NULL;
5866 struct msave_command *msave = &cmd->msave;
5867 msave->expr = matrix_parse_exp (s);
5871 while (lex_match (s->lexer, T_SLASH))
5873 if (lex_match_id (s->lexer, "TYPE"))
5875 lex_match (s->lexer, T_EQUALS);
5877 msave->rowtype = match_rowtype (s->lexer);
5878 if (!msave->rowtype)
5881 else if (lex_match_id (s->lexer, "OUTFILE"))
5883 lex_match (s->lexer, T_EQUALS);
5885 fh_unref (common->outfile);
5886 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5887 if (!common->outfile)
5890 else if (lex_match_id (s->lexer, "VARIABLES"))
5892 if (!parse_var_names (s->lexer, &common->variables))
5895 else if (lex_match_id (s->lexer, "FNAMES"))
5897 if (!parse_var_names (s->lexer, &common->fnames))
5900 else if (lex_match_id (s->lexer, "SNAMES"))
5902 if (!parse_var_names (s->lexer, &common->snames))
5905 else if (lex_match_id (s->lexer, "SPLIT"))
5907 lex_match (s->lexer, T_EQUALS);
5909 matrix_expr_destroy (splits);
5910 splits = matrix_parse_exp (s);
5914 else if (lex_match_id (s->lexer, "FACTOR"))
5916 lex_match (s->lexer, T_EQUALS);
5918 matrix_expr_destroy (factors);
5919 factors = matrix_parse_exp (s);
5925 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5926 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5930 if (!msave->rowtype)
5932 lex_sbc_missing ("TYPE");
5938 if (common->fnames.n && !factors)
5940 msg (SE, _("FNAMES requires FACTOR."));
5943 if (common->snames.n && !splits)
5945 msg (SE, _("SNAMES requires SPLIT."));
5948 if (!common->outfile)
5950 lex_sbc_missing ("OUTFILE");
5957 if (common->outfile)
5959 if (!fh_equal (common->outfile, s->common->outfile))
5961 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5962 "within a single MATRIX command."));
5966 if (!compare_variables ("VARIABLES",
5967 &common->variables, &s->common->variables)
5968 || !compare_variables ("FNAMES", &common->fnames, &s->common->fnames)
5969 || !compare_variables ("SNAMES", &common->snames, &s->common->snames))
5971 msave_common_destroy (common);
5973 msave->common = s->common;
5975 struct msave_common *c = s->common;
5978 if (c->n_factors >= c->allocated_factors)
5979 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
5980 sizeof *c->factors);
5981 c->factors[c->n_factors++] = factors;
5983 if (c->n_factors > 0)
5984 msave->factors = c->factors[c->n_factors - 1];
5988 if (c->n_splits >= c->allocated_splits)
5989 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
5991 c->splits[c->n_splits++] = splits;
5993 if (c->n_splits > 0)
5994 msave->splits = c->splits[c->n_splits - 1];
5999 matrix_expr_destroy (splits);
6000 matrix_expr_destroy (factors);
6001 msave_common_destroy (common);
6002 matrix_cmd_destroy (cmd);
6007 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
6009 gsl_matrix *m = matrix_expr_evaluate (e);
6015 msg (SE, _("%s expression must evaluate to vector, "
6016 "not a %zu×%zu matrix."),
6017 name, m->size1, m->size2);
6018 gsl_matrix_free (m);
6022 return matrix_to_vector (m);
6026 msave_add_vars (struct dictionary *d, const struct string_array *vars)
6028 for (size_t i = 0; i < vars->n; i++)
6029 if (!dict_create_var (d, vars->strings[i], 0))
6030 return vars->strings[i];
6034 static struct dictionary *
6035 msave_create_dict (const struct msave_common *common)
6037 struct dictionary *dict = dict_create (get_default_encoding ());
6039 const char *dup_split = msave_add_vars (dict, &common->snames);
6042 /* Should not be possible because the parser ensures that the names are
6047 dict_create_var_assert (dict, "ROWTYPE_", 8);
6049 const char *dup_factor = msave_add_vars (dict, &common->fnames);
6052 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
6056 dict_create_var_assert (dict, "VARNAME_", 8);
6058 const char *dup_var = msave_add_vars (dict, &common->variables);
6061 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
6073 matrix_cmd_execute_msave (struct msave_command *msave)
6075 struct msave_common *common = msave->common;
6076 gsl_matrix *m = NULL;
6077 gsl_vector *factors = NULL;
6078 gsl_vector *splits = NULL;
6080 m = matrix_expr_evaluate (msave->expr);
6084 if (!common->variables.n)
6085 for (size_t i = 0; i < m->size2; i++)
6086 string_array_append_nocopy (&common->variables,
6087 xasprintf ("COL%zu", i + 1));
6088 else if (m->size2 != common->variables.n)
6091 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
6092 m->size2, common->variables.n);
6098 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
6102 if (!common->fnames.n)
6103 for (size_t i = 0; i < factors->size; i++)
6104 string_array_append_nocopy (&common->fnames,
6105 xasprintf ("FAC%zu", i + 1));
6106 else if (factors->size != common->fnames.n)
6108 msg (SE, _("There are %zu factor variables, "
6109 "but %zu factor values were supplied."),
6110 common->fnames.n, factors->size);
6117 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
6121 if (!common->snames.n)
6122 for (size_t i = 0; i < splits->size; i++)
6123 string_array_append_nocopy (&common->snames,
6124 xasprintf ("SPL%zu", i + 1));
6125 else if (splits->size != common->snames.n)
6127 msg (SE, _("There are %zu split variables, "
6128 "but %zu split values were supplied."),
6129 common->snames.n, splits->size);
6134 if (!common->writer)
6136 struct dictionary *dict = msave_create_dict (common);
6140 common->writer = any_writer_open (common->outfile, dict);
6141 if (!common->writer)
6147 common->dict = dict;
6150 bool matrix = (!strcmp (msave->rowtype, "COV")
6151 || !strcmp (msave->rowtype, "CORR"));
6152 for (size_t y = 0; y < m->size1; y++)
6154 struct ccase *c = case_create (dict_get_proto (common->dict));
6157 /* Split variables */
6159 for (size_t i = 0; i < splits->size; i++)
6160 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
6163 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6164 msave->rowtype, ' ');
6168 for (size_t i = 0; i < factors->size; i++)
6169 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
6172 const char *varname_ = (matrix && y < common->variables.n
6173 ? common->variables.strings[y]
6175 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6178 /* Continuous variables. */
6179 for (size_t x = 0; x < m->size2; x++)
6180 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
6181 casewriter_write (common->writer, c);
6185 gsl_matrix_free (m);
6186 gsl_vector_free (factors);
6187 gsl_vector_free (splits);
6190 static struct matrix_cmd *
6191 matrix_parse_mget (struct matrix_state *s)
6193 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6194 *cmd = (struct matrix_cmd) {
6198 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
6202 struct mget_command *mget = &cmd->mget;
6204 lex_match (s->lexer, T_SLASH);
6205 while (lex_token (s->lexer) != T_ENDCMD)
6207 if (lex_match_id (s->lexer, "FILE"))
6209 lex_match (s->lexer, T_EQUALS);
6211 fh_unref (mget->file);
6212 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6216 else if (lex_match_id (s->lexer, "ENCODING"))
6218 lex_match (s->lexer, T_EQUALS);
6219 if (!lex_force_string (s->lexer))
6222 free (mget->encoding);
6223 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
6227 else if (lex_match_id (s->lexer, "TYPE"))
6229 lex_match (s->lexer, T_EQUALS);
6230 while (lex_token (s->lexer) != T_SLASH
6231 && lex_token (s->lexer) != T_ENDCMD)
6233 const char *rowtype = match_rowtype (s->lexer);
6237 stringi_set_insert (&mget->rowtypes, rowtype);
6242 lex_error_expecting (s->lexer, "FILE", "TYPE");
6245 lex_match (s->lexer, T_SLASH);
6250 matrix_cmd_destroy (cmd);
6254 static const struct variable *
6255 get_a8_var (const struct dictionary *d, const char *name)
6257 const struct variable *v = dict_lookup_var (d, name);
6260 msg (SE, _("Matrix data file lacks %s variable."), name);
6263 if (var_get_width (v) != 8)
6265 msg (SE, _("%s variable in matrix data file must be "
6266 "8-byte string, but it has width %d."),
6267 name, var_get_width (v));
6274 var_changed (const struct ccase *ca, const struct ccase *cb,
6275 const struct variable *var)
6278 ? !value_equal (case_data (ca, var), case_data (cb, var),
6279 var_get_width (var))
6284 vars_changed (const struct ccase *ca, const struct ccase *cb,
6285 const struct dictionary *d,
6286 size_t first_var, size_t n_vars)
6288 for (size_t i = 0; i < n_vars; i++)
6290 const struct variable *v = dict_get_var (d, first_var + i);
6291 if (var_changed (ca, cb, v))
6298 vars_all_missing (const struct ccase *c, const struct dictionary *d,
6299 size_t first_var, size_t n_vars)
6301 for (size_t i = 0; i < n_vars; i++)
6302 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
6308 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6309 const struct dictionary *d,
6310 const struct variable *rowtype_var,
6311 const struct stringi_set *accepted_rowtypes,
6312 struct matrix_state *s,
6313 size_t ss, size_t sn, size_t si,
6314 size_t fs, size_t fn, size_t fi,
6315 size_t cs, size_t cn,
6316 struct pivot_table *pt,
6317 struct pivot_dimension *var_dimension)
6322 /* Is this a matrix for pooled data, either where there are no factor
6323 variables or the factor variables are missing? */
6324 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
6326 struct substring rowtype = case_ss (rows[0], rowtype_var);
6327 ss_rtrim (&rowtype, ss_cstr (" "));
6328 if (!stringi_set_is_empty (accepted_rowtypes)
6329 && !stringi_set_contains_len (accepted_rowtypes,
6330 rowtype.string, rowtype.length))
6333 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
6334 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
6335 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
6336 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
6337 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
6338 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
6342 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
6343 (int) rowtype.length, rowtype.string);
6347 struct string name = DS_EMPTY_INITIALIZER;
6348 ds_put_cstr (&name, prefix);
6350 ds_put_format (&name, "F%zu", fi);
6352 ds_put_format (&name, "S%zu", si);
6354 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6356 mv = matrix_var_create (s, ds_ss (&name));
6359 msg (SW, _("Matrix data file contains variable with existing name %s."),
6361 goto exit_free_name;
6364 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6365 size_t n_missing = 0;
6366 for (size_t y = 0; y < n_rows; y++)
6368 for (size_t x = 0; x < cn; x++)
6370 struct variable *var = dict_get_var (d, cs + x);
6371 double value = case_num (rows[y], var);
6372 if (var_is_num_missing (var, value, MV_ANY))
6377 gsl_matrix_set (m, y, x, value);
6381 int var_index = pivot_category_create_leaf (
6382 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
6383 double values[] = { n_rows, cn };
6384 for (size_t j = 0; j < sn; j++)
6386 struct variable *var = dict_get_var (d, ss + j);
6387 const union value *value = case_data (rows[0], var);
6388 pivot_table_put2 (pt, j, var_index,
6389 pivot_value_new_var_value (var, value));
6391 for (size_t j = 0; j < fn; j++)
6393 struct variable *var = dict_get_var (d, fs + j);
6394 const union value sysmis = { .f = SYSMIS };
6395 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
6396 pivot_table_put2 (pt, j + sn, var_index,
6397 pivot_value_new_var_value (var, value));
6399 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
6400 pivot_table_put2 (pt, j + sn + fn, var_index,
6401 pivot_value_new_integer (values[j]));
6404 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6405 "value, which was treated as zero.",
6406 "Matrix data file variable %s contains %zu missing "
6407 "values, which were treated as zero.", n_missing),
6408 ds_cstr (&name), n_missing);
6415 for (size_t y = 0; y < n_rows; y++)
6416 case_unref (rows[y]);
6420 matrix_cmd_execute_mget__ (struct mget_command *mget,
6421 struct casereader *r, const struct dictionary *d)
6423 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6424 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6425 if (!rowtype_ || !varname_)
6428 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6430 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6433 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6435 msg (SE, _("Matrix data file contains no continuous variables."));
6439 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6441 const struct variable *v = dict_get_var (d, i);
6442 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6445 _("Matrix data file contains unexpected string variable %s."),
6451 /* SPLIT variables. */
6453 size_t sn = var_get_dict_index (rowtype_);
6454 struct ccase *sc = NULL;
6457 /* FACTOR variables. */
6458 size_t fs = var_get_dict_index (rowtype_) + 1;
6459 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6460 struct ccase *fc = NULL;
6463 /* Continuous variables. */
6464 size_t cs = var_get_dict_index (varname_) + 1;
6465 size_t cn = dict_get_var_cnt (d) - cs;
6466 struct ccase *cc = NULL;
6469 struct pivot_table *pt = pivot_table_create (
6470 N_("Matrix Variables Created by MGET"));
6471 struct pivot_dimension *attr_dimension = pivot_dimension_create (
6472 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
6473 struct pivot_dimension *var_dimension = pivot_dimension_create (
6474 pt, PIVOT_AXIS_ROW, N_("Variable"));
6477 struct pivot_category *splits = pivot_category_create_group (
6478 attr_dimension->root, N_("Split Values"));
6479 for (size_t i = 0; i < sn; i++)
6480 pivot_category_create_leaf (splits, pivot_value_new_variable (
6481 dict_get_var (d, ss + i)));
6485 struct pivot_category *factors = pivot_category_create_group (
6486 attr_dimension->root, N_("Factors"));
6487 for (size_t i = 0; i < fn; i++)
6488 pivot_category_create_leaf (factors, pivot_value_new_variable (
6489 dict_get_var (d, fs + i)));
6491 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
6492 N_("Rows"), N_("Columns"));
6495 struct ccase **rows = NULL;
6496 size_t allocated_rows = 0;
6500 while ((c = casereader_read (r)) != NULL)
6502 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
6512 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
6513 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
6514 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
6517 if (change != NOTHING_CHANGED)
6519 matrix_mget_commit_var (rows, n_rows, d,
6520 rowtype_, &mget->rowtypes,
6531 if (n_rows >= allocated_rows)
6532 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6535 if (change == SPLITS_CHANGED)
6541 /* Reset the factor number, if there are factors. */
6545 if (row_has_factors)
6551 else if (change == FACTORS_CHANGED)
6553 if (row_has_factors)
6559 matrix_mget_commit_var (rows, n_rows, d,
6560 rowtype_, &mget->rowtypes,
6572 if (var_dimension->n_leaves)
6573 pivot_table_submit (pt);
6575 pivot_table_unref (pt);
6579 matrix_cmd_execute_mget (struct mget_command *mget)
6581 struct casereader *r;
6582 struct dictionary *d;
6583 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
6584 mget->state->dataset, &r, &d))
6586 matrix_cmd_execute_mget__ (mget, r, d);
6587 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
6592 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6594 if (!lex_force_id (s->lexer))
6597 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6599 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6604 static struct matrix_cmd *
6605 matrix_parse_eigen (struct matrix_state *s)
6607 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6608 *cmd = (struct matrix_cmd) {
6610 .eigen = { .expr = NULL }
6613 struct eigen_command *eigen = &cmd->eigen;
6614 if (!lex_force_match (s->lexer, T_LPAREN))
6616 eigen->expr = matrix_parse_expr (s);
6618 || !lex_force_match (s->lexer, T_COMMA)
6619 || !matrix_parse_dst_var (s, &eigen->evec)
6620 || !lex_force_match (s->lexer, T_COMMA)
6621 || !matrix_parse_dst_var (s, &eigen->eval)
6622 || !lex_force_match (s->lexer, T_RPAREN))
6628 matrix_cmd_destroy (cmd);
6633 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6635 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6638 if (!is_symmetric (A))
6640 msg (SE, _("Argument of EIGEN must be symmetric."));
6641 gsl_matrix_free (A);
6645 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6646 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6647 gsl_vector v_eval = to_vector (eval);
6648 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6649 gsl_eigen_symmv (A, &v_eval, evec, w);
6650 gsl_eigen_symmv_free (w);
6652 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6654 gsl_matrix_free (A);
6656 gsl_matrix_free (eigen->eval->value);
6657 eigen->eval->value = eval;
6659 gsl_matrix_free (eigen->evec->value);
6660 eigen->evec->value = evec;
6663 static struct matrix_cmd *
6664 matrix_parse_setdiag (struct matrix_state *s)
6666 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6667 *cmd = (struct matrix_cmd) {
6668 .type = MCMD_SETDIAG,
6669 .setdiag = { .dst = NULL }
6672 struct setdiag_command *setdiag = &cmd->setdiag;
6673 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6676 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6679 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6684 if (!lex_force_match (s->lexer, T_COMMA))
6687 setdiag->expr = matrix_parse_expr (s);
6691 if (!lex_force_match (s->lexer, T_RPAREN))
6697 matrix_cmd_destroy (cmd);
6702 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6704 gsl_matrix *dst = setdiag->dst->value;
6707 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6708 setdiag->dst->name);
6712 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6716 size_t n = MIN (dst->size1, dst->size2);
6717 if (is_scalar (src))
6719 double d = to_scalar (src);
6720 for (size_t i = 0; i < n; i++)
6721 gsl_matrix_set (dst, i, i, d);
6723 else if (is_vector (src))
6725 gsl_vector v = to_vector (src);
6726 for (size_t i = 0; i < n && i < v.size; i++)
6727 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6730 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6731 "not a %zu×%zu matrix."),
6732 src->size1, src->size2);
6733 gsl_matrix_free (src);
6736 static struct matrix_cmd *
6737 matrix_parse_svd (struct matrix_state *s)
6739 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6740 *cmd = (struct matrix_cmd) {
6742 .svd = { .expr = NULL }
6745 struct svd_command *svd = &cmd->svd;
6746 if (!lex_force_match (s->lexer, T_LPAREN))
6748 svd->expr = matrix_parse_expr (s);
6750 || !lex_force_match (s->lexer, T_COMMA)
6751 || !matrix_parse_dst_var (s, &svd->u)
6752 || !lex_force_match (s->lexer, T_COMMA)
6753 || !matrix_parse_dst_var (s, &svd->s)
6754 || !lex_force_match (s->lexer, T_COMMA)
6755 || !matrix_parse_dst_var (s, &svd->v)
6756 || !lex_force_match (s->lexer, T_RPAREN))
6762 matrix_cmd_destroy (cmd);
6767 matrix_cmd_execute_svd (struct svd_command *svd)
6769 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6773 if (m->size1 >= m->size2)
6776 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6777 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6778 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6779 gsl_vector *work = gsl_vector_alloc (A->size2);
6780 gsl_linalg_SV_decomp (A, V, &Sv, work);
6781 gsl_vector_free (work);
6783 matrix_var_set (svd->u, A);
6784 matrix_var_set (svd->s, S);
6785 matrix_var_set (svd->v, V);
6789 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6790 gsl_matrix_transpose_memcpy (At, m);
6791 gsl_matrix_free (m);
6793 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6794 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6795 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6796 gsl_vector *work = gsl_vector_alloc (At->size2);
6797 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6798 gsl_vector_free (work);
6800 matrix_var_set (svd->v, At);
6801 matrix_var_set (svd->s, St);
6802 matrix_var_set (svd->u, Vt);
6807 matrix_cmd_execute (struct matrix_cmd *cmd)
6812 matrix_cmd_execute_compute (&cmd->compute);
6816 matrix_cmd_execute_print (&cmd->print);
6820 return matrix_cmd_execute_do_if (&cmd->do_if);
6823 matrix_cmd_execute_loop (&cmd->loop);
6830 matrix_cmd_execute_display (&cmd->display);
6834 matrix_cmd_execute_release (&cmd->release);
6838 matrix_cmd_execute_save (&cmd->save);
6842 matrix_cmd_execute_read (&cmd->read);
6846 matrix_cmd_execute_write (&cmd->write);
6850 matrix_cmd_execute_get (&cmd->get);
6854 matrix_cmd_execute_msave (&cmd->msave);
6858 matrix_cmd_execute_mget (&cmd->mget);
6862 matrix_cmd_execute_eigen (&cmd->eigen);
6866 matrix_cmd_execute_setdiag (&cmd->setdiag);
6870 matrix_cmd_execute_svd (&cmd->svd);
6878 matrix_cmds_uninit (struct matrix_cmds *cmds)
6880 for (size_t i = 0; i < cmds->n; i++)
6881 matrix_cmd_destroy (cmds->commands[i]);
6882 free (cmds->commands);
6886 matrix_cmd_destroy (struct matrix_cmd *cmd)
6894 matrix_lvalue_destroy (cmd->compute.lvalue);
6895 matrix_expr_destroy (cmd->compute.rvalue);
6899 matrix_expr_destroy (cmd->print.expression);
6900 free (cmd->print.title);
6901 print_labels_destroy (cmd->print.rlabels);
6902 print_labels_destroy (cmd->print.clabels);
6906 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6908 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6909 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6911 free (cmd->do_if.clauses);
6915 matrix_expr_destroy (cmd->loop.start);
6916 matrix_expr_destroy (cmd->loop.end);
6917 matrix_expr_destroy (cmd->loop.increment);
6918 matrix_expr_destroy (cmd->loop.top_condition);
6919 matrix_expr_destroy (cmd->loop.bottom_condition);
6920 matrix_cmds_uninit (&cmd->loop.commands);
6930 free (cmd->release.vars);
6934 matrix_expr_destroy (cmd->save.expression);
6938 matrix_lvalue_destroy (cmd->read.dst);
6939 matrix_expr_destroy (cmd->read.size);
6943 matrix_expr_destroy (cmd->write.expression);
6944 free (cmd->write.format);
6948 matrix_lvalue_destroy (cmd->get.dst);
6949 fh_unref (cmd->get.file);
6950 free (cmd->get.encoding);
6951 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
6955 matrix_expr_destroy (cmd->msave.expr);
6959 fh_unref (cmd->mget.file);
6960 stringi_set_destroy (&cmd->mget.rowtypes);
6964 matrix_expr_destroy (cmd->eigen.expr);
6968 matrix_expr_destroy (cmd->setdiag.expr);
6972 matrix_expr_destroy (cmd->svd.expr);
6978 struct matrix_command_name
6981 struct matrix_cmd *(*parse) (struct matrix_state *);
6984 static const struct matrix_command_name *
6985 matrix_parse_command_name (struct lexer *lexer)
6987 static const struct matrix_command_name commands[] = {
6988 { "COMPUTE", matrix_parse_compute },
6989 { "CALL EIGEN", matrix_parse_eigen },
6990 { "CALL SETDIAG", matrix_parse_setdiag },
6991 { "CALL SVD", matrix_parse_svd },
6992 { "PRINT", matrix_parse_print },
6993 { "DO IF", matrix_parse_do_if },
6994 { "LOOP", matrix_parse_loop },
6995 { "BREAK", matrix_parse_break },
6996 { "READ", matrix_parse_read },
6997 { "WRITE", matrix_parse_write },
6998 { "GET", matrix_parse_get },
6999 { "SAVE", matrix_parse_save },
7000 { "MGET", matrix_parse_mget },
7001 { "MSAVE", matrix_parse_msave },
7002 { "DISPLAY", matrix_parse_display },
7003 { "RELEASE", matrix_parse_release },
7005 static size_t n = sizeof commands / sizeof *commands;
7007 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
7008 if (lex_match_phrase (lexer, c->name))
7013 static struct matrix_cmd *
7014 matrix_parse_command (struct matrix_state *s)
7016 size_t nesting_level = SIZE_MAX;
7018 struct matrix_cmd *c = NULL;
7019 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
7021 lex_error (s->lexer, _("Unknown matrix command."));
7022 else if (!cmd->parse)
7023 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
7027 nesting_level = output_open_group (
7028 group_item_create_nocopy (utf8_to_title (cmd->name),
7029 utf8_to_title (cmd->name)));
7034 lex_end_of_command (s->lexer);
7035 lex_discard_rest_of_command (s->lexer);
7036 if (nesting_level != SIZE_MAX)
7037 output_close_groups (nesting_level);
7043 cmd_matrix (struct lexer *lexer, struct dataset *ds)
7045 if (!lex_force_match (lexer, T_ENDCMD))
7048 struct matrix_state state = {
7050 .session = dataset_session (ds),
7052 .vars = HMAP_INITIALIZER (state.vars),
7057 while (lex_match (lexer, T_ENDCMD))
7059 if (lex_token (lexer) == T_STOP)
7061 msg (SE, _("Unexpected end of input expecting matrix command."));
7065 if (lex_match_phrase (lexer, "END MATRIX"))
7068 struct matrix_cmd *c = matrix_parse_command (&state);
7071 matrix_cmd_execute (c);
7072 matrix_cmd_destroy (c);
7076 struct matrix_var *var, *next;
7077 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
7080 gsl_matrix_free (var->value);
7081 hmap_delete (&state.vars, &var->hmap_node);
7084 hmap_destroy (&state.vars);
7085 msave_common_destroy (state.common);
7086 fh_unref (state.prev_read_file);
7087 for (size_t i = 0; i < state.n_read_files; i++)
7088 read_file_destroy (state.read_files[i]);
7089 free (state.read_files);
7090 fh_unref (state.prev_write_file);
7091 for (size_t i = 0; i < state.n_write_files; i++)
7092 write_file_destroy (state.write_files[i]);
7093 free (state.write_files);
7094 fh_unref (state.prev_save_file);
7095 for (size_t i = 0; i < state.n_save_files; i++)
7096 save_file_destroy (state.save_files[i]);
7097 free (state.save_files);