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;
88 /* Execution state. */
89 struct dictionary *dict;
90 struct casewriter *writer;
95 struct file_handle *file;
96 struct dfm_reader *reader;
102 struct file_handle *file;
103 struct dfm_writer *writer;
105 struct u8_line *held;
110 struct file_handle *file;
112 /* Parameters from parsing the first SAVE command for 'file'. */
113 struct string_array variables;
114 struct matrix_expr *names;
115 struct stringi_set strings;
117 /* Results from the first attempt to open this save_file. */
119 struct casewriter *writer;
124 struct dataset *dataset;
125 struct session *session;
129 struct msave_common *common;
131 struct file_handle *prev_read_file;
132 struct read_file **read_files;
135 struct file_handle *prev_write_file;
136 struct write_file **write_files;
137 size_t n_write_files;
139 struct file_handle *prev_save_file;
140 struct save_file **save_files;
144 static struct matrix_var *
145 matrix_var_lookup (struct matrix_state *s, struct substring name)
147 struct matrix_var *var;
149 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
150 utf8_hash_case_substring (name, 0), &s->vars)
151 if (!utf8_sscasecmp (ss_cstr (var->name), name))
157 static struct matrix_var *
158 matrix_var_create (struct matrix_state *s, struct substring name)
160 struct matrix_var *var = xmalloc (sizeof *var);
161 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
162 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
167 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
169 gsl_matrix_free (var->value);
173 #define MATRIX_FUNCTIONS \
229 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
230 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
231 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
232 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
233 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
234 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
235 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
236 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
237 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
238 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
239 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
241 typedef gsl_matrix *matrix_proto_m_d (double);
242 typedef gsl_matrix *matrix_proto_m_dd (double, double);
243 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
244 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
245 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
246 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
247 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
248 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
249 typedef double matrix_proto_d_m (gsl_matrix *);
250 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
251 typedef gsl_matrix *matrix_proto_IDENT (double, double);
253 #define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
257 /* Matrix expressions. */
264 #define F(NAME, PROTOTYPE) MOP_F_##NAME,
268 /* Elementwise and scalar arithmetic. */
269 MOP_NEGATE, /* unary - */
270 MOP_ADD_ELEMS, /* + */
271 MOP_SUB_ELEMS, /* - */
272 MOP_MUL_ELEMS, /* &* */
273 MOP_DIV_ELEMS, /* / and &/ */
274 MOP_EXP_ELEMS, /* &** */
276 MOP_SEQ_BY, /* a:b:c */
278 /* Matrix arithmetic. */
280 MOP_EXP_MAT, /* ** */
297 MOP_PASTE_HORZ, /* a, b, c, ... */
298 MOP_PASTE_VERT, /* a; b; c; ... */
302 MOP_VEC_INDEX, /* x(y) */
303 MOP_VEC_ALL, /* x(:) */
304 MOP_MAT_INDEX, /* x(y,z) */
305 MOP_ROW_INDEX, /* x(y,:) */
306 MOP_COL_INDEX, /* x(:,z) */
313 MOP_EOF, /* EOF('file') */
321 struct matrix_expr **subs;
326 struct matrix_var *variable;
327 struct read_file *eof;
332 matrix_expr_destroy (struct matrix_expr *e)
339 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
370 for (size_t i = 0; i < e->n_subs; i++)
371 matrix_expr_destroy (e->subs[i]);
383 static struct matrix_expr *
384 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
387 struct matrix_expr *e = xmalloc (sizeof *e);
388 *e = (struct matrix_expr) {
390 .subs = xmemdup (subs, n_subs * sizeof *subs),
396 static struct matrix_expr *
397 matrix_expr_create_0 (enum matrix_op op)
399 struct matrix_expr *sub;
400 return matrix_expr_create_subs (op, &sub, 0);
403 static struct matrix_expr *
404 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
406 return matrix_expr_create_subs (op, &sub, 1);
409 static struct matrix_expr *
410 matrix_expr_create_2 (enum matrix_op op,
411 struct matrix_expr *sub0, struct matrix_expr *sub1)
413 struct matrix_expr *subs[] = { sub0, sub1 };
414 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
417 static struct matrix_expr *
418 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
419 struct matrix_expr *sub1, struct matrix_expr *sub2)
421 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
422 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
425 static struct matrix_expr *
426 matrix_expr_create_number (double number)
428 struct matrix_expr *e = xmalloc (sizeof *e);
429 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
433 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
435 static struct matrix_expr *
436 matrix_parse_curly_comma (struct matrix_state *s)
438 struct matrix_expr *lhs = matrix_parse_or_xor (s);
442 while (lex_match (s->lexer, T_COMMA))
444 struct matrix_expr *rhs = matrix_parse_or_xor (s);
447 matrix_expr_destroy (lhs);
450 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
455 static struct matrix_expr *
456 matrix_parse_curly_semi (struct matrix_state *s)
458 if (lex_token (s->lexer) == T_RCURLY)
459 return matrix_expr_create_0 (MOP_EMPTY);
461 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
465 while (lex_match (s->lexer, T_SEMICOLON))
467 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
470 matrix_expr_destroy (lhs);
473 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
478 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
479 for (size_t Y = 0; Y < (M)->size1; Y++) \
480 for (size_t X = 0; X < (M)->size2; X++) \
481 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
484 is_vector (const gsl_matrix *m)
486 return m->size1 <= 1 || m->size2 <= 1;
490 to_vector (gsl_matrix *m)
492 return (m->size1 == 1
493 ? gsl_matrix_row (m, 0).vector
494 : gsl_matrix_column (m, 0).vector);
499 matrix_eval_ABS (gsl_matrix *m)
501 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
507 matrix_eval_ALL (gsl_matrix *m)
509 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
516 matrix_eval_ANY (gsl_matrix *m)
518 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
525 matrix_eval_ARSIN (gsl_matrix *m)
527 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
533 matrix_eval_ARTAN (gsl_matrix *m)
535 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
541 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
545 for (size_t i = 0; i < n; i++)
550 gsl_matrix *block = gsl_matrix_calloc (r, c);
552 for (size_t i = 0; i < n; i++)
554 for (size_t y = 0; y < m[i]->size1; y++)
555 for (size_t x = 0; x < m[i]->size2; x++)
556 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
564 matrix_eval_CHOL (gsl_matrix *m)
566 if (!gsl_linalg_cholesky_decomp1 (m))
568 for (size_t y = 0; y < m->size1; y++)
569 for (size_t x = y + 1; x < m->size2; x++)
570 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
572 for (size_t y = 0; y < m->size1; y++)
573 for (size_t x = 0; x < y; x++)
574 gsl_matrix_set (m, y, x, 0);
579 msg (SE, _("Input to CHOL function is not positive-definite."));
585 matrix_eval_col_extremum (gsl_matrix *m, bool min)
590 return gsl_matrix_alloc (1, 0);
592 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
593 for (size_t x = 0; x < m->size2; x++)
595 double ext = gsl_matrix_get (m, 0, x);
596 for (size_t y = 1; y < m->size1; y++)
598 double value = gsl_matrix_get (m, y, x);
599 if (min ? value < ext : value > ext)
602 gsl_matrix_set (cext, 0, x, ext);
608 matrix_eval_CMAX (gsl_matrix *m)
610 return matrix_eval_col_extremum (m, false);
614 matrix_eval_CMIN (gsl_matrix *m)
616 return matrix_eval_col_extremum (m, true);
620 matrix_eval_COS (gsl_matrix *m)
622 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
628 matrix_eval_col_sum (gsl_matrix *m, bool square)
633 return gsl_matrix_alloc (1, 0);
635 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
636 for (size_t x = 0; x < m->size2; x++)
639 for (size_t y = 0; y < m->size1; y++)
641 double d = gsl_matrix_get (m, y, x);
642 sum += square ? pow2 (d) : d;
644 gsl_matrix_set (result, 0, x, sum);
650 matrix_eval_CSSQ (gsl_matrix *m)
652 return matrix_eval_col_sum (m, true);
656 matrix_eval_CSUM (gsl_matrix *m)
658 return matrix_eval_col_sum (m, false);
662 compare_double_3way (const void *a_, const void *b_)
664 const double *a = a_;
665 const double *b = b_;
666 return *a < *b ? -1 : *a > *b;
670 matrix_eval_DESIGN (gsl_matrix *m)
672 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
673 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
674 gsl_matrix_transpose_memcpy (&m2, m);
676 for (size_t y = 0; y < m2.size1; y++)
677 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
679 size_t *n = xcalloc (m2.size1, sizeof *n);
681 for (size_t i = 0; i < m2.size1; i++)
683 double *row = tmp + m2.size2 * i;
684 for (size_t j = 0; j < m2.size2; )
687 for (k = j + 1; k < m2.size2; k++)
688 if (row[j] != row[k])
690 row[n[i]++] = row[j];
695 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
700 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
702 for (size_t i = 0; i < m->size2; i++)
707 const double *unique = tmp + m2.size2 * i;
708 for (size_t j = 0; j < n[i]; j++, x++)
710 double value = unique[j];
711 for (size_t y = 0; y < m->size1; y++)
712 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
723 matrix_eval_DET (gsl_matrix *m)
725 gsl_permutation *p = gsl_permutation_alloc (m->size1);
727 gsl_linalg_LU_decomp (m, p, &signum);
728 gsl_permutation_free (p);
729 return gsl_linalg_LU_det (m, signum);
733 matrix_eval_DIAG (gsl_matrix *m)
735 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
736 for (size_t i = 0; i < diag->size1; i++)
737 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
742 is_symmetric (const gsl_matrix *m)
744 if (m->size1 != m->size2)
747 for (size_t y = 0; y < m->size1; y++)
748 for (size_t x = 0; x < y; x++)
749 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
756 compare_double_desc (const void *a_, const void *b_)
758 const double *a = a_;
759 const double *b = b_;
760 return *a > *b ? -1 : *a < *b;
764 matrix_eval_EVAL (gsl_matrix *m)
766 if (!is_symmetric (m))
768 msg (SE, _("Argument of EVAL must be symmetric."));
772 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
773 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
774 gsl_vector v_eval = to_vector (eval);
775 gsl_eigen_symm (m, &v_eval, w);
776 gsl_eigen_symm_free (w);
778 assert (v_eval.stride == 1);
779 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
785 matrix_eval_EXP (gsl_matrix *m)
787 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
792 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
795 Charl Linssen <charl@itfromb.it>
799 matrix_eval_GINV (gsl_matrix *A)
804 gsl_matrix *tmp_mat = NULL;
807 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
808 tmp_mat = gsl_matrix_alloc (m, n);
809 gsl_matrix_transpose_memcpy (tmp_mat, A);
817 gsl_matrix *V = gsl_matrix_alloc (m, m);
818 gsl_vector *u = gsl_vector_alloc (m);
820 gsl_vector *tmp_vec = gsl_vector_alloc (m);
821 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
822 gsl_vector_free (tmp_vec);
825 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
826 gsl_matrix_set_zero (Sigma_pinv);
827 double cutoff = 1e-15 * gsl_vector_max (u);
829 for (size_t i = 0; i < m; ++i)
831 double x = gsl_vector_get (u, i);
832 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
835 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
836 gsl_matrix *U = gsl_matrix_calloc (n, n);
837 for (size_t i = 0; i < n; i++)
838 for (size_t j = 0; j < m; j++)
839 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
841 /* two dot products to obtain pseudoinverse */
842 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
843 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
848 A_pinv = gsl_matrix_alloc (n, m);
849 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
853 A_pinv = gsl_matrix_alloc (m, n);
854 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
857 gsl_matrix_free (tmp_mat);
858 gsl_matrix_free (tmp_mat2);
860 gsl_matrix_free (Sigma_pinv);
874 grade_compare_3way (const void *a_, const void *b_)
876 const struct grade *a = a_;
877 const struct grade *b = b_;
879 return (a->value < b->value ? -1
880 : a->value > b->value ? 1
888 matrix_eval_GRADE (gsl_matrix *m)
890 size_t n = m->size1 * m->size2;
891 struct grade *grades = xmalloc (n * sizeof *grades);
894 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
895 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
896 qsort (grades, n, sizeof *grades, grade_compare_3way);
898 for (size_t i = 0; i < n; i++)
899 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
907 dot (gsl_vector *a, gsl_vector *b)
910 for (size_t i = 0; i < a->size; i++)
911 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
916 norm2 (gsl_vector *v)
919 for (size_t i = 0; i < v->size; i++)
920 result += pow2 (gsl_vector_get (v, i));
927 return sqrt (norm2 (v));
931 matrix_eval_GSCH (gsl_matrix *v)
933 if (v->size2 < v->size1)
935 msg (SE, _("GSCH requires its argument to have at least as many columns "
936 "as rows, but it has dimensions %zu×%zu."),
940 if (!v->size1 || !v->size2)
943 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
945 for (size_t vx = 0; vx < v->size2; vx++)
947 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
948 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
950 gsl_vector_memcpy (&u_i, &v_i);
951 for (size_t j = 0; j < ux; j++)
953 gsl_vector u_j = gsl_matrix_column (u, j).vector;
954 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
955 for (size_t k = 0; k < u_i.size; k++)
956 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
957 - scale * gsl_vector_get (&u_j, k)));
960 double len = norm (&u_i);
963 gsl_vector_scale (&u_i, 1.0 / len);
964 if (++ux >= v->size1)
971 msg (SE, _("%zu×%zu argument to GSCH contains only "
972 "%zu linearly independent columns."),
973 v->size1, v->size2, ux);
983 matrix_eval_IDENT (double s1, double s2)
985 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
987 msg (SE, _("Arguments to IDENT must be non-negative integers."));
991 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
992 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
998 invert_matrix (gsl_matrix *x)
1000 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1002 gsl_linalg_LU_decomp (x, p, &signum);
1003 gsl_linalg_LU_invx (x, p);
1004 gsl_permutation_free (p);
1008 matrix_eval_INV (gsl_matrix *m)
1015 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1017 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1018 a->size2 * b->size2);
1020 for (size_t ar = 0; ar < a->size1; ar++)
1021 for (size_t br = 0; br < b->size1; br++, y++)
1024 for (size_t ac = 0; ac < a->size2; ac++)
1025 for (size_t bc = 0; bc < b->size2; bc++, x++)
1027 double av = gsl_matrix_get (a, ar, ac);
1028 double bv = gsl_matrix_get (b, br, bc);
1029 gsl_matrix_set (k, y, x, av * bv);
1036 matrix_eval_LG10 (gsl_matrix *m)
1038 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1044 matrix_eval_LN (gsl_matrix *m)
1046 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1052 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1054 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1057 for (size_t i = 1; i <= n * n; i++)
1059 gsl_matrix_set (m, y, x, i);
1061 size_t y1 = !y ? n - 1 : y - 1;
1062 size_t x1 = x + 1 >= n ? 0 : x + 1;
1063 if (gsl_matrix_get (m, y1, x1) == 0)
1069 y = y + 1 >= n ? 0 : y + 1;
1074 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1076 double a = gsl_matrix_get (m, y1, x1);
1077 double b = gsl_matrix_get (m, y2, x2);
1078 gsl_matrix_set (m, y1, x1, b);
1079 gsl_matrix_set (m, y2, x2, a);
1083 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1087 /* A. Umar, "On the Construction of Even Order Magic Squares",
1088 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1090 for (size_t i = 1; i <= n * n / 2; i++)
1092 gsl_matrix_set (m, y, x, i);
1102 for (size_t i = n * n; i > n * n / 2; i--)
1104 gsl_matrix_set (m, y, x, i);
1112 for (size_t y = 0; y < n; y++)
1113 for (size_t x = 0; x < n / 2; x++)
1115 unsigned int d = gsl_matrix_get (m, y, x);
1116 if (d % 2 != (y < n / 2))
1117 magic_exchange (m, y, x, y, n - x - 1);
1122 size_t x1 = n / 2 - 1;
1124 magic_exchange (m, y1, x1, y2, x1);
1125 magic_exchange (m, y1, x2, y2, x2);
1129 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1131 /* A. Umar, "On the Construction of Even Order Magic Squares",
1132 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1136 for (size_t i = 1; ; i++)
1138 gsl_matrix_set (m, y, x, i);
1139 if (++y == n / 2 - 1)
1151 for (size_t i = n * n; ; i--)
1153 gsl_matrix_set (m, y, x, i);
1154 if (++y == n / 2 - 1)
1163 for (size_t y = 0; y < n; y++)
1164 if (y != n / 2 - 1 && y != n / 2)
1165 for (size_t x = 0; x < n / 2; x++)
1167 unsigned int d = gsl_matrix_get (m, y, x);
1168 if (d % 2 != (y < n / 2))
1169 magic_exchange (m, y, x, y, n - x - 1);
1172 size_t a0 = (n * n - 2 * n) / 2 + 1;
1173 for (size_t i = 0; i < n / 2; i++)
1176 gsl_matrix_set (m, n / 2, i, a);
1177 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1179 for (size_t i = 0; i < n / 2; i++)
1181 size_t a = a0 + i + n / 2;
1182 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1183 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1185 for (size_t x = 1; x < n / 2; x += 2)
1186 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1187 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1188 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1189 size_t x1 = n / 2 - 2;
1190 size_t x2 = n / 2 + 1;
1191 size_t y1 = n / 2 - 2;
1192 size_t y2 = n / 2 + 1;
1193 magic_exchange (m, y1, x1, y2, x1);
1194 magic_exchange (m, y1, x2, y2, x2);
1198 matrix_eval_MAGIC (double n_)
1200 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1202 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1207 gsl_matrix *m = gsl_matrix_calloc (n, n);
1209 matrix_eval_MAGIC_odd (m, n);
1211 matrix_eval_MAGIC_singly_even (m, n);
1213 matrix_eval_MAGIC_doubly_even (m, n);
1218 matrix_eval_MAKE (double r, double c, double value)
1220 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1222 msg (SE, _("Size arguments to MAKE must be integers."));
1226 gsl_matrix *m = gsl_matrix_alloc (r, c);
1227 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1233 matrix_eval_MDIAG (gsl_vector *v)
1235 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1236 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1237 gsl_vector_memcpy (&diagonal, v);
1242 matrix_eval_MMAX (gsl_matrix *m)
1244 return gsl_matrix_max (m);
1248 matrix_eval_MMIN (gsl_matrix *m)
1250 return gsl_matrix_min (m);
1254 matrix_eval_MOD (gsl_matrix *m, double divisor)
1258 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1262 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1263 *d = fmod (*d, divisor);
1268 matrix_eval_MSSQ (gsl_matrix *m)
1271 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1277 matrix_eval_MSUM (gsl_matrix *m)
1280 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1286 matrix_eval_NCOL (gsl_matrix *m)
1292 matrix_eval_NROW (gsl_matrix *m)
1298 matrix_eval_RANK (gsl_matrix *m)
1300 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1301 gsl_linalg_QR_decomp (m, tau);
1302 gsl_vector_free (tau);
1304 return gsl_linalg_QRPT_rank (m, -1);
1308 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1310 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1312 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1317 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1319 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1320 "product of matrix dimensions."));
1324 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1327 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1329 gsl_matrix_set (dst, y1, x1, *d);
1340 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1345 return gsl_matrix_alloc (0, 1);
1347 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1348 for (size_t y = 0; y < m->size1; y++)
1350 double ext = gsl_matrix_get (m, y, 0);
1351 for (size_t x = 1; x < m->size2; x++)
1353 double value = gsl_matrix_get (m, y, x);
1354 if (min ? value < ext : value > ext)
1357 gsl_matrix_set (rext, y, 0, ext);
1363 matrix_eval_RMAX (gsl_matrix *m)
1365 return matrix_eval_row_extremum (m, false);
1369 matrix_eval_RMIN (gsl_matrix *m)
1371 return matrix_eval_row_extremum (m, true);
1375 matrix_eval_RND (gsl_matrix *m)
1377 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1389 rank_compare_3way (const void *a_, const void *b_)
1391 const struct rank *a = a_;
1392 const struct rank *b = b_;
1394 return a->value < b->value ? -1 : a->value > b->value;
1398 matrix_eval_RNKORDER (gsl_matrix *m)
1400 size_t n = m->size1 * m->size2;
1401 struct rank *ranks = xmalloc (n * sizeof *ranks);
1403 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1404 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1405 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1407 for (size_t i = 0; i < n; )
1410 for (j = i + 1; j < n; j++)
1411 if (ranks[i].value != ranks[j].value)
1414 double rank = (i + j + 1.0) / 2.0;
1415 for (size_t k = i; k < j; k++)
1416 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1427 matrix_eval_row_sum (gsl_matrix *m, bool square)
1432 return gsl_matrix_alloc (0, 1);
1434 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1435 for (size_t y = 0; y < m->size1; y++)
1438 for (size_t x = 0; x < m->size2; x++)
1440 double d = gsl_matrix_get (m, y, x);
1441 sum += square ? pow2 (d) : d;
1443 gsl_matrix_set (result, y, 0, sum);
1449 matrix_eval_RSSQ (gsl_matrix *m)
1451 return matrix_eval_row_sum (m, true);
1455 matrix_eval_RSUM (gsl_matrix *m)
1457 return matrix_eval_row_sum (m, false);
1461 matrix_eval_SIN (gsl_matrix *m)
1463 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1469 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1471 if (m1->size1 != m2->size1)
1473 msg (SE, _("SOLVE requires its arguments to have the same number of "
1474 "rows, but the first argument has dimensions %zu×%zu and "
1475 "the second %zu×%zu."),
1476 m1->size1, m1->size2,
1477 m2->size1, m2->size2);
1481 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1482 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1484 gsl_linalg_LU_decomp (m1, p, &signum);
1485 for (size_t i = 0; i < m2->size2; i++)
1487 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1488 gsl_vector xi = gsl_matrix_column (x, i).vector;
1489 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1491 gsl_permutation_free (p);
1496 matrix_eval_SQRT (gsl_matrix *m)
1498 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1502 msg (SE, _("Argument to SQRT must be nonnegative."));
1511 matrix_eval_SSCP (gsl_matrix *m)
1513 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1514 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1519 matrix_eval_SVAL (gsl_matrix *m)
1521 gsl_matrix *tmp_mat = NULL;
1522 if (m->size2 > m->size1)
1524 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1525 gsl_matrix_transpose_memcpy (tmp_mat, m);
1530 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1531 gsl_vector *S = gsl_vector_alloc (m->size2);
1532 gsl_vector *work = gsl_vector_alloc (m->size2);
1533 gsl_linalg_SV_decomp (m, V, S, work);
1535 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1536 for (size_t i = 0; i < m->size2; i++)
1537 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1539 gsl_matrix_free (V);
1540 gsl_vector_free (S);
1541 gsl_vector_free (work);
1542 gsl_matrix_free (tmp_mat);
1548 matrix_eval_SWEEP (gsl_matrix *m, double d)
1550 if (d < 1 || d > SIZE_MAX)
1552 msg (SE, _("Scalar argument to SWEEP must be integer."));
1556 if (k >= MIN (m->size1, m->size2))
1558 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1559 "equal to the smaller of the matrix argument's rows and "
1564 double m_kk = gsl_matrix_get (m, k, k);
1565 if (fabs (m_kk) > 1e-19)
1567 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1568 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1570 double m_ij = gsl_matrix_get (m, i, j);
1571 double m_ik = gsl_matrix_get (m, i, k);
1572 double m_kj = gsl_matrix_get (m, k, j);
1573 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1582 for (size_t i = 0; i < m->size1; i++)
1584 gsl_matrix_set (m, i, k, 0);
1585 gsl_matrix_set (m, k, i, 0);
1592 matrix_eval_TRACE (gsl_matrix *m)
1595 size_t n = MIN (m->size1, m->size2);
1596 for (size_t i = 0; i < n; i++)
1597 sum += gsl_matrix_get (m, i, i);
1602 matrix_eval_T (gsl_matrix *m)
1604 return matrix_eval_TRANSPOS (m);
1608 matrix_eval_TRANSPOS (gsl_matrix *m)
1610 if (m->size1 == m->size2)
1612 gsl_matrix_transpose (m);
1617 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1618 gsl_matrix_transpose_memcpy (t, m);
1624 matrix_eval_TRUNC (gsl_matrix *m)
1626 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1632 matrix_eval_UNIFORM (double r_, double c_)
1634 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1636 msg (SE, _("Arguments to UNIFORM must be integers."));
1641 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1643 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1647 gsl_matrix *m = gsl_matrix_alloc (r, c);
1648 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1649 *d = gsl_ran_flat (get_rng (), 0, 1);
1653 struct matrix_function
1657 size_t min_args, max_args;
1660 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1662 static const struct matrix_function *
1663 matrix_parse_function_name (const char *token)
1665 static const struct matrix_function functions[] =
1667 #define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1671 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1673 for (size_t i = 0; i < N_FUNCTIONS; i++)
1675 if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1676 return &functions[i];
1681 static struct read_file *
1682 read_file_create (struct matrix_state *s, struct file_handle *fh)
1684 for (size_t i = 0; i < s->n_read_files; i++)
1686 struct read_file *rf = s->read_files[i];
1694 struct read_file *rf = xmalloc (sizeof *rf);
1695 *rf = (struct read_file) { .file = fh };
1697 s->read_files = xrealloc (s->read_files,
1698 (s->n_read_files + 1) * sizeof *s->read_files);
1699 s->read_files[s->n_read_files++] = rf;
1703 static struct dfm_reader *
1704 read_file_open (struct read_file *rf)
1707 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1712 read_file_destroy (struct read_file *rf)
1716 fh_unref (rf->file);
1717 dfm_close_reader (rf->reader);
1718 free (rf->encoding);
1723 static struct write_file *
1724 write_file_create (struct matrix_state *s, struct file_handle *fh)
1726 for (size_t i = 0; i < s->n_write_files; i++)
1728 struct write_file *wf = s->write_files[i];
1736 struct write_file *wf = xmalloc (sizeof *wf);
1737 *wf = (struct write_file) { .file = fh };
1739 s->write_files = xrealloc (s->write_files,
1740 (s->n_write_files + 1) * sizeof *s->write_files);
1741 s->write_files[s->n_write_files++] = wf;
1745 static struct dfm_writer *
1746 write_file_open (struct write_file *wf)
1749 wf->writer = dfm_open_writer (wf->file, wf->encoding);
1754 write_file_destroy (struct write_file *wf)
1760 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
1761 wf->held->s.ss.length);
1762 u8_line_destroy (wf->held);
1766 fh_unref (wf->file);
1767 dfm_close_writer (wf->writer);
1768 free (wf->encoding);
1774 matrix_parse_function (struct matrix_state *s, const char *token,
1775 struct matrix_expr **exprp)
1778 if (lex_next_token (s->lexer, 1) != T_LPAREN)
1781 if (lex_match_id (s->lexer, "EOF"))
1784 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1788 if (!lex_force_match (s->lexer, T_RPAREN))
1794 struct read_file *rf = read_file_create (s, fh);
1796 struct matrix_expr *e = xmalloc (sizeof *e);
1797 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1802 const struct matrix_function *f = matrix_parse_function_name (token);
1806 lex_get_n (s->lexer, 2);
1808 struct matrix_expr *e = xmalloc (sizeof *e);
1809 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1811 size_t allocated_subs = 0;
1814 struct matrix_expr *sub = matrix_parse_expr (s);
1818 if (e->n_subs >= allocated_subs)
1819 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1820 e->subs[e->n_subs++] = sub;
1822 while (lex_match (s->lexer, T_COMMA));
1823 if (!lex_force_match (s->lexer, T_RPAREN))
1830 matrix_expr_destroy (e);
1834 static struct matrix_expr *
1835 matrix_parse_primary (struct matrix_state *s)
1837 if (lex_is_number (s->lexer))
1839 double number = lex_number (s->lexer);
1841 return matrix_expr_create_number (number);
1843 else if (lex_is_string (s->lexer))
1845 char string[sizeof (double)];
1846 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1850 memcpy (&number, string, sizeof number);
1851 return matrix_expr_create_number (number);
1853 else if (lex_match (s->lexer, T_LPAREN))
1855 struct matrix_expr *e = matrix_parse_or_xor (s);
1856 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1858 matrix_expr_destroy (e);
1863 else if (lex_match (s->lexer, T_LCURLY))
1865 struct matrix_expr *e = matrix_parse_curly_semi (s);
1866 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1868 matrix_expr_destroy (e);
1873 else if (lex_token (s->lexer) == T_ID)
1875 struct matrix_expr *retval;
1876 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1879 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1882 lex_error (s->lexer, _("Unknown variable %s."),
1883 lex_tokcstr (s->lexer));
1888 struct matrix_expr *e = xmalloc (sizeof *e);
1889 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1892 else if (lex_token (s->lexer) == T_ALL)
1894 struct matrix_expr *retval;
1895 if (matrix_parse_function (s, "ALL", &retval))
1899 lex_error (s->lexer, NULL);
1903 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1906 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1908 if (lex_match (s->lexer, T_COLON))
1915 *indexp = matrix_parse_expr (s);
1916 return *indexp != NULL;
1920 static struct matrix_expr *
1921 matrix_parse_postfix (struct matrix_state *s)
1923 struct matrix_expr *lhs = matrix_parse_primary (s);
1924 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1927 struct matrix_expr *i0;
1928 if (!matrix_parse_index_expr (s, &i0))
1930 matrix_expr_destroy (lhs);
1933 if (lex_match (s->lexer, T_RPAREN))
1935 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1936 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1937 else if (lex_match (s->lexer, T_COMMA))
1939 struct matrix_expr *i1;
1940 if (!matrix_parse_index_expr (s, &i1)
1941 || !lex_force_match (s->lexer, T_RPAREN))
1943 matrix_expr_destroy (lhs);
1944 matrix_expr_destroy (i0);
1945 matrix_expr_destroy (i1);
1948 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1949 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1950 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1955 lex_error_expecting (s->lexer, "`)'", "`,'");
1960 static struct matrix_expr *
1961 matrix_parse_unary (struct matrix_state *s)
1963 if (lex_match (s->lexer, T_DASH))
1965 struct matrix_expr *lhs = matrix_parse_unary (s);
1966 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1968 else if (lex_match (s->lexer, T_PLUS))
1969 return matrix_parse_unary (s);
1971 return matrix_parse_postfix (s);
1974 static struct matrix_expr *
1975 matrix_parse_seq (struct matrix_state *s)
1977 struct matrix_expr *start = matrix_parse_unary (s);
1978 if (!start || !lex_match (s->lexer, T_COLON))
1981 struct matrix_expr *end = matrix_parse_unary (s);
1984 matrix_expr_destroy (start);
1988 if (lex_match (s->lexer, T_COLON))
1990 struct matrix_expr *increment = matrix_parse_unary (s);
1993 matrix_expr_destroy (start);
1994 matrix_expr_destroy (end);
1997 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2000 return matrix_expr_create_2 (MOP_SEQ, start, end);
2003 static struct matrix_expr *
2004 matrix_parse_exp (struct matrix_state *s)
2006 struct matrix_expr *lhs = matrix_parse_seq (s);
2013 if (lex_match (s->lexer, T_EXP))
2015 else if (lex_match_phrase (s->lexer, "&**"))
2020 struct matrix_expr *rhs = matrix_parse_seq (s);
2023 matrix_expr_destroy (lhs);
2026 lhs = matrix_expr_create_2 (op, lhs, rhs);
2030 static struct matrix_expr *
2031 matrix_parse_mul_div (struct matrix_state *s)
2033 struct matrix_expr *lhs = matrix_parse_exp (s);
2040 if (lex_match (s->lexer, T_ASTERISK))
2042 else if (lex_match (s->lexer, T_SLASH))
2044 else if (lex_match_phrase (s->lexer, "&*"))
2046 else if (lex_match_phrase (s->lexer, "&/"))
2051 struct matrix_expr *rhs = matrix_parse_exp (s);
2054 matrix_expr_destroy (lhs);
2057 lhs = matrix_expr_create_2 (op, lhs, rhs);
2061 static struct matrix_expr *
2062 matrix_parse_add_sub (struct matrix_state *s)
2064 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2071 if (lex_match (s->lexer, T_PLUS))
2073 else if (lex_match (s->lexer, T_DASH))
2075 else if (lex_token (s->lexer) == T_NEG_NUM)
2080 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2083 matrix_expr_destroy (lhs);
2086 lhs = matrix_expr_create_2 (op, lhs, rhs);
2090 static struct matrix_expr *
2091 matrix_parse_relational (struct matrix_state *s)
2093 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2100 if (lex_match (s->lexer, T_GT))
2102 else if (lex_match (s->lexer, T_GE))
2104 else if (lex_match (s->lexer, T_LT))
2106 else if (lex_match (s->lexer, T_LE))
2108 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2110 else if (lex_match (s->lexer, T_NE))
2115 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2118 matrix_expr_destroy (lhs);
2121 lhs = matrix_expr_create_2 (op, lhs, rhs);
2125 static struct matrix_expr *
2126 matrix_parse_not (struct matrix_state *s)
2128 if (lex_match (s->lexer, T_NOT))
2130 struct matrix_expr *lhs = matrix_parse_not (s);
2131 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2134 return matrix_parse_relational (s);
2137 static struct matrix_expr *
2138 matrix_parse_and (struct matrix_state *s)
2140 struct matrix_expr *lhs = matrix_parse_not (s);
2144 while (lex_match (s->lexer, T_AND))
2146 struct matrix_expr *rhs = matrix_parse_not (s);
2149 matrix_expr_destroy (lhs);
2152 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2157 static struct matrix_expr *
2158 matrix_parse_or_xor (struct matrix_state *s)
2160 struct matrix_expr *lhs = matrix_parse_and (s);
2167 if (lex_match (s->lexer, T_OR))
2169 else if (lex_match_id (s->lexer, "XOR"))
2174 struct matrix_expr *rhs = matrix_parse_and (s);
2177 matrix_expr_destroy (lhs);
2180 lhs = matrix_expr_create_2 (op, lhs, rhs);
2184 static struct matrix_expr *
2185 matrix_parse_expr (struct matrix_state *s)
2187 return matrix_parse_or_xor (s);
2190 /* Expression evaluation. */
2193 matrix_op_eval (enum matrix_op op, double a, double b)
2197 case MOP_ADD_ELEMS: return a + b;
2198 case MOP_SUB_ELEMS: return a - b;
2199 case MOP_MUL_ELEMS: return a * b;
2200 case MOP_DIV_ELEMS: return a / b;
2201 case MOP_EXP_ELEMS: return pow (a, b);
2202 case MOP_GT: return a > b;
2203 case MOP_GE: return a >= b;
2204 case MOP_LT: return a < b;
2205 case MOP_LE: return a <= b;
2206 case MOP_EQ: return a == b;
2207 case MOP_NE: return a != b;
2208 case MOP_AND: return (a > 0) && (b > 0);
2209 case MOP_OR: return (a > 0) || (b > 0);
2210 case MOP_XOR: return (a > 0) != (b > 0);
2212 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2221 case MOP_PASTE_HORZ:
2222 case MOP_PASTE_VERT:
2238 matrix_op_name (enum matrix_op op)
2242 case MOP_ADD_ELEMS: return "+";
2243 case MOP_SUB_ELEMS: return "-";
2244 case MOP_MUL_ELEMS: return "&*";
2245 case MOP_DIV_ELEMS: return "&/";
2246 case MOP_EXP_ELEMS: return "&**";
2247 case MOP_GT: return ">";
2248 case MOP_GE: return ">=";
2249 case MOP_LT: return "<";
2250 case MOP_LE: return "<=";
2251 case MOP_EQ: return "=";
2252 case MOP_NE: return "<>";
2253 case MOP_AND: return "AND";
2254 case MOP_OR: return "OR";
2255 case MOP_XOR: return "XOR";
2257 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2266 case MOP_PASTE_HORZ:
2267 case MOP_PASTE_VERT:
2283 is_scalar (const gsl_matrix *m)
2285 return m->size1 == 1 && m->size2 == 1;
2289 to_scalar (const gsl_matrix *m)
2291 assert (is_scalar (m));
2292 return gsl_matrix_get (m, 0, 0);
2296 matrix_expr_evaluate_elementwise (enum matrix_op op,
2297 gsl_matrix *a, gsl_matrix *b)
2301 double be = to_scalar (b);
2302 for (size_t r = 0; r < a->size1; r++)
2303 for (size_t c = 0; c < a->size2; c++)
2305 double *ae = gsl_matrix_ptr (a, r, c);
2306 *ae = matrix_op_eval (op, *ae, be);
2310 else if (is_scalar (a))
2312 double ae = to_scalar (a);
2313 for (size_t r = 0; r < b->size1; r++)
2314 for (size_t c = 0; c < b->size2; c++)
2316 double *be = gsl_matrix_ptr (b, r, c);
2317 *be = matrix_op_eval (op, ae, *be);
2321 else if (a->size1 == b->size1 && a->size2 == b->size2)
2323 for (size_t r = 0; r < a->size1; r++)
2324 for (size_t c = 0; c < a->size2; c++)
2326 double *ae = gsl_matrix_ptr (a, r, c);
2327 double be = gsl_matrix_get (b, r, c);
2328 *ae = matrix_op_eval (op, *ae, be);
2334 msg (SE, _("Operands to %s must have the same dimensions or one "
2335 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
2336 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2342 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2344 if (is_scalar (a) || is_scalar (b))
2345 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2347 if (a->size2 != b->size1)
2349 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
2350 "not conformable for multiplication."),
2351 a->size1, a->size2, b->size1, b->size2);
2355 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2356 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2361 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2363 gsl_matrix *tmp = *a;
2369 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2372 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2373 swap_matrix (z, tmp);
2377 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2379 mul_matrix (x, *x, *x, tmp);
2383 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2386 if (x->size1 != x->size2)
2388 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2389 "the left-hand size, not one with dimensions %zu×%zu."),
2390 x->size1, x->size2);
2395 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2396 "right-hand side, not a matrix with dimensions %zu×%zu."),
2397 b->size1, b->size2);
2400 double bf = to_scalar (b);
2401 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2403 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2404 "or outside the valid range."), bf);
2409 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2411 gsl_matrix_set_identity (y);
2415 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2417 for (unsigned long int n = labs (bl); n > 1; n /= 2)
2420 mul_matrix (&y, x, y, &t);
2421 square_matrix (&x, &t);
2424 square_matrix (&x, &t);
2426 mul_matrix (&y, x, y, &t);
2430 /* Garbage collection.
2432 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2433 a permutation of them. We are returning one of them; that one must not be
2434 destroyed. We must not destroy 'x_' because the caller owns it. */
2436 gsl_matrix_free (y_);
2438 gsl_matrix_free (t_);
2444 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2447 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2449 msg (SE, _("All operands of : operator must be scalars."));
2453 long int start = to_scalar (start_);
2454 long int end = to_scalar (end_);
2455 long int by = by_ ? to_scalar (by_) : 1;
2459 msg (SE, _("The increment operand to : must be nonzero."));
2463 long int n = (end >= start && by > 0 ? (end - start + by) / by
2464 : end <= start && by < 0 ? (start - end - by) / -by
2466 gsl_matrix *m = gsl_matrix_alloc (1, n);
2467 for (long int i = 0; i < n; i++)
2468 gsl_matrix_set (m, 0, i, start + i * by);
2473 matrix_expr_evaluate_not (gsl_matrix *a)
2475 for (size_t r = 0; r < a->size1; r++)
2476 for (size_t c = 0; c < a->size2; c++)
2478 double *ae = gsl_matrix_ptr (a, r, c);
2485 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2487 if (a->size1 != b->size1)
2489 if (!a->size1 || !a->size2)
2491 else if (!b->size1 || !b->size2)
2494 msg (SE, _("All columns in a matrix must have the same number of rows, "
2495 "but this tries to paste matrices with %zu and %zu rows."),
2496 a->size1, b->size1);
2500 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2501 for (size_t y = 0; y < a->size1; y++)
2503 for (size_t x = 0; x < a->size2; x++)
2504 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2505 for (size_t x = 0; x < b->size2; x++)
2506 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2512 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2514 if (a->size2 != b->size2)
2516 if (!a->size1 || !a->size2)
2518 else if (!b->size1 || !b->size2)
2521 msg (SE, _("All rows in a matrix must have the same number of columns, "
2522 "but this tries to stack matrices with %zu and %zu columns."),
2523 a->size2, b->size2);
2527 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2528 for (size_t x = 0; x < a->size2; x++)
2530 for (size_t y = 0; y < a->size1; y++)
2531 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2532 for (size_t y = 0; y < b->size1; y++)
2533 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2539 matrix_to_vector (gsl_matrix *m)
2542 gsl_vector v = to_vector (m);
2543 assert (v.block == m->block || !v.block);
2547 return xmemdup (&v, sizeof v);
2561 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2564 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2565 enum index_type index_type, size_t other_size,
2566 struct index_vector *iv)
2575 msg (SE, _("Vector index must be scalar or vector, not a "
2577 m->size1, m->size2);
2581 msg (SE, _("Matrix row index must be scalar or vector, not a "
2583 m->size1, m->size2);
2587 msg (SE, _("Matrix column index must be scalar or vector, not a "
2589 m->size1, m->size2);
2595 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2596 *iv = (struct index_vector) {
2597 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2600 for (size_t i = 0; i < v.size; i++)
2602 double index = gsl_vector_get (&v, i);
2603 if (index < 1 || index >= size + 1)
2608 msg (SE, _("Index %g is out of range for vector "
2609 "with %zu elements."), index, size);
2613 msg (SE, _("%g is not a valid row index for "
2614 "a %zu×%zu matrix."),
2615 index, size, other_size);
2619 msg (SE, _("%g is not a valid column index for "
2620 "a %zu×%zu matrix."),
2621 index, other_size, size);
2628 iv->indexes[i] = index - 1;
2634 *iv = (struct index_vector) {
2635 .indexes = xnmalloc (size, sizeof *iv->indexes),
2638 for (size_t i = 0; i < size; i++)
2645 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2647 if (!is_vector (sm))
2649 msg (SE, _("Vector index operator may not be applied to "
2650 "a %zu×%zu matrix."),
2651 sm->size1, sm->size2);
2659 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2661 if (!matrix_expr_evaluate_vec_all (sm))
2664 gsl_vector sv = to_vector (sm);
2665 struct index_vector iv;
2666 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2669 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2670 sm->size1 == 1 ? iv.n : 1);
2671 gsl_vector dv = to_vector (dm);
2672 for (size_t dx = 0; dx < iv.n; dx++)
2674 size_t sx = iv.indexes[dx];
2675 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2683 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2686 struct index_vector iv0;
2687 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2690 struct index_vector iv1;
2691 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2698 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2699 for (size_t dy = 0; dy < iv0.n; dy++)
2701 size_t sy = iv0.indexes[dy];
2703 for (size_t dx = 0; dx < iv1.n; dx++)
2705 size_t sx = iv1.indexes[dx];
2706 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2714 #define F(NAME, PROTOTYPE) \
2715 static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2716 const char *name, gsl_matrix *[], size_t, \
2717 matrix_proto_##PROTOTYPE *);
2722 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2724 if (!is_scalar (subs[index]))
2726 msg (SE, _("Function %s argument %zu must be a scalar, "
2727 "not a %zu×%zu matrix."),
2728 name, index + 1, subs[index]->size1, subs[index]->size2);
2735 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2737 if (!is_vector (subs[index]))
2739 msg (SE, _("Function %s argument %zu must be a vector, "
2740 "not a %zu×%zu matrix."),
2741 name, index + 1, subs[index]->size1, subs[index]->size2);
2748 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2750 for (size_t i = 0; i < n_subs; i++)
2752 if (!check_scalar_arg (name, subs, i))
2754 d[i] = to_scalar (subs[i]);
2760 matrix_expr_evaluate_m_d (const char *name,
2761 gsl_matrix *subs[], size_t n_subs,
2762 matrix_proto_m_d *f)
2764 assert (n_subs == 1);
2767 return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2771 matrix_expr_evaluate_m_dd (const char *name,
2772 gsl_matrix *subs[], size_t n_subs,
2773 matrix_proto_m_dd *f)
2775 assert (n_subs == 2);
2778 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2782 matrix_expr_evaluate_m_ddd (const char *name,
2783 gsl_matrix *subs[], size_t n_subs,
2784 matrix_proto_m_ddd *f)
2786 assert (n_subs == 3);
2789 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2793 matrix_expr_evaluate_m_m (const char *name UNUSED,
2794 gsl_matrix *subs[], size_t n_subs,
2795 matrix_proto_m_m *f)
2797 assert (n_subs == 1);
2802 matrix_expr_evaluate_m_md (const char *name UNUSED,
2803 gsl_matrix *subs[], size_t n_subs,
2804 matrix_proto_m_md *f)
2806 assert (n_subs == 2);
2807 if (!check_scalar_arg (name, subs, 1))
2809 return f (subs[0], to_scalar (subs[1]));
2813 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2814 gsl_matrix *subs[], size_t n_subs,
2815 matrix_proto_m_mdd *f)
2817 assert (n_subs == 3);
2818 if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2820 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2824 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2825 gsl_matrix *subs[], size_t n_subs,
2826 matrix_proto_m_mm *f)
2828 assert (n_subs == 2);
2829 return f (subs[0], subs[1]);
2833 matrix_expr_evaluate_m_v (const char *name,
2834 gsl_matrix *subs[], size_t n_subs,
2835 matrix_proto_m_v *f)
2837 assert (n_subs == 1);
2838 if (!check_vector_arg (name, subs, 0))
2840 gsl_vector v = to_vector (subs[0]);
2845 matrix_expr_evaluate_d_m (const char *name UNUSED,
2846 gsl_matrix *subs[], size_t n_subs,
2847 matrix_proto_d_m *f)
2849 assert (n_subs == 1);
2851 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2852 gsl_matrix_set (m, 0, 0, f (subs[0]));
2857 matrix_expr_evaluate_m_any (const char *name UNUSED,
2858 gsl_matrix *subs[], size_t n_subs,
2859 matrix_proto_m_any *f)
2861 return f (subs, n_subs);
2865 matrix_expr_evaluate_IDENT (const char *name,
2866 gsl_matrix *subs[], size_t n_subs,
2867 matrix_proto_IDENT *f)
2869 assert (n_subs <= 2);
2872 if (!to_scalar_args (name, subs, n_subs, d))
2875 return f (d[0], d[n_subs - 1]);
2879 matrix_expr_evaluate (const struct matrix_expr *e)
2881 if (e->op == MOP_NUMBER)
2883 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2884 gsl_matrix_set (m, 0, 0, e->number);
2887 else if (e->op == MOP_VARIABLE)
2889 const gsl_matrix *src = e->variable->value;
2892 msg (SE, _("Uninitialized variable %s used in expression."),
2897 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2898 gsl_matrix_memcpy (dst, src);
2901 else if (e->op == MOP_EOF)
2903 struct dfm_reader *reader = read_file_open (e->eof);
2904 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2905 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2909 enum { N_LOCAL = 3 };
2910 gsl_matrix *local_subs[N_LOCAL];
2911 gsl_matrix **subs = (e->n_subs < N_LOCAL
2913 : xmalloc (e->n_subs * sizeof *subs));
2915 for (size_t i = 0; i < e->n_subs; i++)
2917 subs[i] = matrix_expr_evaluate (e->subs[i]);
2920 for (size_t j = 0; j < i; j++)
2921 gsl_matrix_free (subs[j]);
2922 if (subs != local_subs)
2928 gsl_matrix *result = NULL;
2931 #define F(NAME, PROTOTYPE) \
2932 case MOP_F_##NAME: \
2933 result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \
2935 matrix_eval_##NAME); \
2941 gsl_matrix_scale (subs[0], -1.0);
2959 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2963 result = matrix_expr_evaluate_not (subs[0]);
2967 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2971 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2975 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2979 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2982 case MOP_PASTE_HORZ:
2983 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2986 case MOP_PASTE_VERT:
2987 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2991 result = gsl_matrix_alloc (0, 0);
2995 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
2999 result = matrix_expr_evaluate_vec_all (subs[0]);
3003 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
3007 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
3011 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
3020 for (size_t i = 0; i < e->n_subs; i++)
3021 if (subs[i] != result)
3022 gsl_matrix_free (subs[i]);
3023 if (subs != local_subs)
3029 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
3032 gsl_matrix *m = matrix_expr_evaluate (e);
3038 msg (SE, _("Expression for %s must evaluate to scalar, "
3039 "not a %zu×%zu matrix."),
3040 context, m->size1, m->size2);
3041 gsl_matrix_free (m);
3046 gsl_matrix_free (m);
3051 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
3055 if (!matrix_expr_evaluate_scalar (e, context, &d))
3059 if (d < LONG_MIN || d > LONG_MAX)
3061 msg (SE, _("Expression for %s is outside the integer range."), context);
3068 struct matrix_lvalue
3070 struct matrix_var *var;
3071 struct matrix_expr *indexes[2];
3076 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
3080 for (size_t i = 0; i < lvalue->n_indexes; i++)
3081 matrix_expr_destroy (lvalue->indexes[i]);
3086 static struct matrix_lvalue *
3087 matrix_lvalue_parse (struct matrix_state *s)
3089 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3090 if (!lex_force_id (s->lexer))
3092 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3093 if (lex_next_token (s->lexer, 1) == T_LPAREN)
3097 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3102 lex_get_n (s->lexer, 2);
3104 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3106 lvalue->n_indexes++;
3108 if (lex_match (s->lexer, T_COMMA))
3110 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3112 lvalue->n_indexes++;
3114 if (!lex_force_match (s->lexer, T_RPAREN))
3120 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3126 matrix_lvalue_destroy (lvalue);
3131 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3132 enum index_type index_type, size_t other_size,
3133 struct index_vector *iv)
3138 m = matrix_expr_evaluate (e);
3145 bool ok = matrix_normalize_index_vector (m, size, index_type,
3147 gsl_matrix_free (m);
3152 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3153 struct index_vector *iv, gsl_matrix *sm)
3155 gsl_vector dv = to_vector (lvalue->var->value);
3157 /* Convert source matrix 'sm' to source vector 'sv'. */
3158 if (!is_vector (sm))
3160 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
3161 sm->size1, sm->size2);
3164 gsl_vector sv = to_vector (sm);
3165 if (iv->n != sv.size)
3167 msg (SE, _("Can't assign %zu-element vector "
3168 "to %zu-element subvector."), sv.size, iv->n);
3172 /* Assign elements. */
3173 for (size_t x = 0; x < iv->n; x++)
3174 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3179 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3180 struct index_vector *iv0,
3181 struct index_vector *iv1, gsl_matrix *sm)
3183 gsl_matrix *dm = lvalue->var->value;
3185 /* Convert source matrix 'sm' to source vector 'sv'. */
3186 if (iv0->n != sm->size1)
3188 msg (SE, _("Row index vector for assignment to %s has %zu elements "
3189 "but source matrix has %zu rows."),
3190 lvalue->var->name, iv0->n, sm->size1);
3193 if (iv1->n != sm->size2)
3195 msg (SE, _("Column index vector for assignment to %s has %zu elements "
3196 "but source matrix has %zu columns."),
3197 lvalue->var->name, iv1->n, sm->size2);
3201 /* Assign elements. */
3202 for (size_t y = 0; y < iv0->n; y++)
3204 size_t f0 = iv0->indexes[y];
3205 for (size_t x = 0; x < iv1->n; x++)
3207 size_t f1 = iv1->indexes[x];
3208 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3214 /* Takes ownership of M. */
3216 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3217 struct index_vector *iv0, struct index_vector *iv1,
3220 if (!lvalue->n_indexes)
3222 gsl_matrix_free (lvalue->var->value);
3223 lvalue->var->value = sm;
3228 bool ok = (lvalue->n_indexes == 1
3229 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3230 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3231 free (iv0->indexes);
3232 free (iv1->indexes);
3233 gsl_matrix_free (sm);
3239 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3240 struct index_vector *iv0,
3241 struct index_vector *iv1)
3243 *iv0 = INDEX_VECTOR_INIT;
3244 *iv1 = INDEX_VECTOR_INIT;
3245 if (!lvalue->n_indexes)
3248 /* Validate destination matrix exists and has the right shape. */
3249 gsl_matrix *dm = lvalue->var->value;
3252 msg (SE, _("Undefined variable %s."), lvalue->var->name);
3255 else if (dm->size1 == 0 || dm->size2 == 0)
3257 msg (SE, _("Cannot index %zu×%zu matrix."),
3258 dm->size1, dm->size2);
3261 else if (lvalue->n_indexes == 1)
3263 if (!is_vector (dm))
3265 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
3266 dm->size1, dm->size2, lvalue->var->name);
3269 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3270 MAX (dm->size1, dm->size2),
3275 assert (lvalue->n_indexes == 2);
3276 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3277 IV_ROW, dm->size2, iv0))
3280 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3281 IV_COLUMN, dm->size1, iv1))
3283 free (iv0->indexes);
3290 /* Takes ownership of M. */
3292 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3294 struct index_vector iv0, iv1;
3295 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3297 gsl_matrix_free (sm);
3301 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3305 /* Matrix command. */
3309 struct matrix_cmd **commands;
3313 static void matrix_cmds_uninit (struct matrix_cmds *);
3317 enum matrix_cmd_type
3340 struct compute_command
3342 struct matrix_lvalue *lvalue;
3343 struct matrix_expr *rvalue;
3347 struct print_command
3349 struct matrix_expr *expression;
3350 bool use_default_format;
3351 struct fmt_spec format;
3353 int space; /* -1 means NEWPAGE. */
3357 struct string_array literals; /* CLABELS/RLABELS. */
3358 struct matrix_expr *expr; /* CNAMES/RNAMES. */
3364 struct do_if_command
3368 struct matrix_expr *condition;
3369 struct matrix_cmds commands;
3379 struct matrix_var *var;
3380 struct matrix_expr *start;
3381 struct matrix_expr *end;
3382 struct matrix_expr *increment;
3384 /* Loop conditions. */
3385 struct matrix_expr *top_condition;
3386 struct matrix_expr *bottom_condition;
3389 struct matrix_cmds commands;
3393 struct display_command
3395 struct matrix_state *state;
3401 struct release_command
3403 struct matrix_var **vars;
3410 struct matrix_expr *expression;
3411 struct save_file *sf;
3417 struct read_file *rf;
3418 struct matrix_lvalue *dst;
3419 struct matrix_expr *size;
3421 enum fmt_type format;
3428 struct write_command
3430 struct write_file *wf;
3431 struct matrix_expr *expression;
3434 /* If this is nonnull, WRITE uses this format.
3436 If this is NULL, WRITE uses free-field format with as many
3437 digits of precision as needed. */
3438 struct fmt_spec *format;
3447 struct matrix_lvalue *dst;
3448 struct dataset *dataset;
3449 struct file_handle *file;
3451 struct var_syntax *vars;
3453 struct matrix_var *names;
3455 /* Treatment of missing values. */
3460 MGET_ERROR, /* Flag error on command. */
3461 MGET_ACCEPT, /* Accept without change, user-missing only. */
3462 MGET_OMIT, /* Drop this case. */
3463 MGET_RECODE /* Recode to 'substitute'. */
3466 double substitute; /* MGET_RECODE only. */
3472 struct msave_command
3474 struct msave_common *common;
3476 struct matrix_expr *expr;
3477 const char *rowtype;
3478 struct matrix_expr *factors;
3479 struct matrix_expr *splits;
3485 struct matrix_state *state;
3486 struct file_handle *file;
3487 struct stringi_set rowtypes;
3491 struct eigen_command
3493 struct matrix_expr *expr;
3494 struct matrix_var *evec;
3495 struct matrix_var *eval;
3499 struct setdiag_command
3501 struct matrix_var *dst;
3502 struct matrix_expr *expr;
3508 struct matrix_expr *expr;
3509 struct matrix_var *u;
3510 struct matrix_var *s;
3511 struct matrix_var *v;
3517 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3518 static bool matrix_cmd_execute (struct matrix_cmd *);
3519 static void matrix_cmd_destroy (struct matrix_cmd *);
3522 static struct matrix_cmd *
3523 matrix_parse_compute (struct matrix_state *s)
3525 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3526 *cmd = (struct matrix_cmd) {
3527 .type = MCMD_COMPUTE,
3528 .compute = { .lvalue = NULL },
3531 struct compute_command *compute = &cmd->compute;
3532 compute->lvalue = matrix_lvalue_parse (s);
3533 if (!compute->lvalue)
3536 if (!lex_force_match (s->lexer, T_EQUALS))
3539 compute->rvalue = matrix_parse_expr (s);
3540 if (!compute->rvalue)
3546 matrix_cmd_destroy (cmd);
3551 matrix_cmd_execute_compute (struct compute_command *compute)
3553 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3557 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3561 print_labels_destroy (struct print_labels *labels)
3565 string_array_destroy (&labels->literals);
3566 matrix_expr_destroy (labels->expr);
3572 parse_literal_print_labels (struct matrix_state *s,
3573 struct print_labels **labelsp)
3575 lex_match (s->lexer, T_EQUALS);
3576 print_labels_destroy (*labelsp);
3577 *labelsp = xzalloc (sizeof **labelsp);
3578 while (lex_token (s->lexer) != T_SLASH
3579 && lex_token (s->lexer) != T_ENDCMD
3580 && lex_token (s->lexer) != T_STOP)
3582 struct string label = DS_EMPTY_INITIALIZER;
3583 while (lex_token (s->lexer) != T_COMMA
3584 && lex_token (s->lexer) != T_SLASH
3585 && lex_token (s->lexer) != T_ENDCMD
3586 && lex_token (s->lexer) != T_STOP)
3588 if (!ds_is_empty (&label))
3589 ds_put_byte (&label, ' ');
3591 if (lex_token (s->lexer) == T_STRING)
3592 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3595 char *rep = lex_next_representation (s->lexer, 0, 0);
3596 ds_put_cstr (&label, rep);
3601 string_array_append_nocopy (&(*labelsp)->literals,
3602 ds_steal_cstr (&label));
3604 if (!lex_match (s->lexer, T_COMMA))
3610 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3612 lex_match (s->lexer, T_EQUALS);
3613 struct matrix_expr *e = matrix_parse_exp (s);
3617 print_labels_destroy (*labelsp);
3618 *labelsp = xzalloc (sizeof **labelsp);
3619 (*labelsp)->expr = e;
3623 static struct matrix_cmd *
3624 matrix_parse_print (struct matrix_state *s)
3626 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3627 *cmd = (struct matrix_cmd) {
3630 .use_default_format = true,
3634 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3637 for (size_t i = 0; ; i++)
3639 enum token_type t = lex_next_token (s->lexer, i);
3640 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3642 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3644 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3647 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3652 cmd->print.expression = matrix_parse_exp (s);
3653 if (!cmd->print.expression)
3657 while (lex_match (s->lexer, T_SLASH))
3659 if (lex_match_id (s->lexer, "FORMAT"))
3661 lex_match (s->lexer, T_EQUALS);
3662 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3664 cmd->print.use_default_format = false;
3666 else if (lex_match_id (s->lexer, "TITLE"))
3668 lex_match (s->lexer, T_EQUALS);
3669 if (!lex_force_string (s->lexer))
3671 free (cmd->print.title);
3672 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3675 else if (lex_match_id (s->lexer, "SPACE"))
3677 lex_match (s->lexer, T_EQUALS);
3678 if (lex_match_id (s->lexer, "NEWPAGE"))
3679 cmd->print.space = -1;
3680 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3682 cmd->print.space = lex_integer (s->lexer);
3688 else if (lex_match_id (s->lexer, "RLABELS"))
3689 parse_literal_print_labels (s, &cmd->print.rlabels);
3690 else if (lex_match_id (s->lexer, "CLABELS"))
3691 parse_literal_print_labels (s, &cmd->print.clabels);
3692 else if (lex_match_id (s->lexer, "RNAMES"))
3694 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3697 else if (lex_match_id (s->lexer, "CNAMES"))
3699 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3704 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3705 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3713 matrix_cmd_destroy (cmd);
3718 matrix_is_integer (const gsl_matrix *m)
3720 for (size_t y = 0; y < m->size1; y++)
3721 for (size_t x = 0; x < m->size2; x++)
3723 double d = gsl_matrix_get (m, y, x);
3731 matrix_max_magnitude (const gsl_matrix *m)
3734 for (size_t y = 0; y < m->size1; y++)
3735 for (size_t x = 0; x < m->size2; x++)
3737 double d = fabs (gsl_matrix_get (m, y, x));
3745 format_fits (struct fmt_spec format, double x)
3747 char *s = data_out (&(union value) { .f = x }, NULL,
3748 &format, settings_get_fmt_settings ());
3749 bool fits = *s != '*' && !strchr (s, 'E');
3754 static struct fmt_spec
3755 default_format (const gsl_matrix *m, int *log_scale)
3759 double max = matrix_max_magnitude (m);
3761 if (matrix_is_integer (m))
3762 for (int w = 1; w <= 12; w++)
3764 struct fmt_spec format = { .type = FMT_F, .w = w };
3765 if (format_fits (format, -max))
3769 if (max >= 1e9 || max <= 1e-4)
3772 snprintf (s, sizeof s, "%.1e", max);
3774 const char *e = strchr (s, 'e');
3776 *log_scale = atoi (e + 1);
3779 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3783 trimmed_string (double d)
3785 struct substring s = ss_buffer ((char *) &d, sizeof d);
3786 ss_rtrim (&s, ss_cstr (" "));
3787 return ss_xstrdup (s);
3790 static struct string_array *
3791 print_labels_get (const struct print_labels *labels, size_t n,
3792 const char *prefix, bool truncate)
3797 struct string_array *out = xzalloc (sizeof *out);
3798 if (labels->literals.n)
3799 string_array_clone (out, &labels->literals);
3800 else if (labels->expr)
3802 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3803 if (m && is_vector (m))
3805 gsl_vector v = to_vector (m);
3806 for (size_t i = 0; i < v.size; i++)
3807 string_array_append_nocopy (out, trimmed_string (
3808 gsl_vector_get (&v, i)));
3810 gsl_matrix_free (m);
3816 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3819 string_array_append (out, "");
3822 string_array_delete (out, out->n - 1);
3825 for (size_t i = 0; i < out->n; i++)
3827 char *s = out->strings[i];
3828 s[strnlen (s, 8)] = '\0';
3835 matrix_cmd_print_space (int space)
3838 output_item_submit (page_break_item_create ());
3839 for (int i = 0; i < space; i++)
3840 output_log ("%s", "");
3844 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3845 struct fmt_spec format, int log_scale)
3847 matrix_cmd_print_space (print->space);
3849 output_log ("%s", print->title);
3851 output_log (" 10 ** %d X", log_scale);
3853 struct string_array *clabels = print_labels_get (print->clabels,
3854 m->size2, "col", true);
3855 if (clabels && format.w < 8)
3857 struct string_array *rlabels = print_labels_get (print->rlabels,
3858 m->size1, "row", true);
3862 struct string line = DS_EMPTY_INITIALIZER;
3864 ds_put_byte_multiple (&line, ' ', 8);
3865 for (size_t x = 0; x < m->size2; x++)
3866 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3867 output_log_nocopy (ds_steal_cstr (&line));
3870 double scale = pow (10.0, log_scale);
3871 bool numeric = fmt_is_numeric (format.type);
3872 for (size_t y = 0; y < m->size1; y++)
3874 struct string line = DS_EMPTY_INITIALIZER;
3876 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3878 for (size_t x = 0; x < m->size2; x++)
3880 double f = gsl_matrix_get (m, y, x);
3882 ? data_out (&(union value) { .f = f / scale}, NULL,
3883 &format, settings_get_fmt_settings ())
3884 : trimmed_string (f));
3885 ds_put_format (&line, " %s", s);
3888 output_log_nocopy (ds_steal_cstr (&line));
3891 string_array_destroy (rlabels);
3893 string_array_destroy (clabels);
3898 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3899 const struct print_labels *print_labels, size_t n,
3900 const char *name, const char *prefix)
3902 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3904 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3905 for (size_t i = 0; i < n; i++)
3906 pivot_category_create_leaf (
3908 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3909 : pivot_value_new_integer (i + 1)));
3911 d->hide_all_labels = true;
3912 string_array_destroy (labels);
3917 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3918 struct fmt_spec format, int log_scale)
3920 struct pivot_table *table = pivot_table_create__ (
3921 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3924 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3926 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3927 N_("Columns"), "col");
3929 struct pivot_footnote *footnote = NULL;
3932 char *s = xasprintf ("× 10**%d", log_scale);
3933 footnote = pivot_table_create_footnote (
3934 table, pivot_value_new_user_text_nocopy (s));
3937 double scale = pow (10.0, log_scale);
3938 bool numeric = fmt_is_numeric (format.type);
3939 for (size_t y = 0; y < m->size1; y++)
3940 for (size_t x = 0; x < m->size2; x++)
3942 double f = gsl_matrix_get (m, y, x);
3943 struct pivot_value *v;
3946 v = pivot_value_new_number (f / scale);
3947 v->numeric.format = format;
3950 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3952 pivot_value_add_footnote (v, footnote);
3953 pivot_table_put2 (table, y, x, v);
3956 pivot_table_submit (table);
3960 matrix_cmd_execute_print (const struct print_command *print)
3962 if (print->expression)
3964 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3969 struct fmt_spec format = (print->use_default_format
3970 ? default_format (m, &log_scale)
3973 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3974 matrix_cmd_print_text (print, m, format, log_scale);
3976 matrix_cmd_print_tables (print, m, format, log_scale);
3978 gsl_matrix_free (m);
3982 matrix_cmd_print_space (print->space);
3984 output_log ("%s", print->title);
3991 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3992 const char *command_name,
3993 const char *stop1, const char *stop2)
3995 lex_end_of_command (s->lexer);
3996 lex_discard_rest_of_command (s->lexer);
3998 size_t allocated = 0;
4001 while (lex_token (s->lexer) == T_ENDCMD)
4004 if (lex_at_phrase (s->lexer, stop1)
4005 || (stop2 && lex_at_phrase (s->lexer, stop2)))
4008 if (lex_at_phrase (s->lexer, "END MATRIX"))
4010 msg (SE, _("Premature END MATRIX within %s."), command_name);
4014 struct matrix_cmd *cmd = matrix_parse_command (s);
4018 if (c->n >= allocated)
4019 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
4020 c->commands[c->n++] = cmd;
4025 matrix_parse_do_if_clause (struct matrix_state *s,
4026 struct do_if_command *ifc,
4027 bool parse_condition,
4028 size_t *allocated_clauses)
4030 if (ifc->n_clauses >= *allocated_clauses)
4031 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
4032 sizeof *ifc->clauses);
4033 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
4034 *c = (struct do_if_clause) { .condition = NULL };
4036 if (parse_condition)
4038 c->condition = matrix_parse_expr (s);
4043 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
4046 static struct matrix_cmd *
4047 matrix_parse_do_if (struct matrix_state *s)
4049 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4050 *cmd = (struct matrix_cmd) {
4052 .do_if = { .n_clauses = 0 }
4055 size_t allocated_clauses = 0;
4058 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
4061 while (lex_match_phrase (s->lexer, "ELSE IF"));
4063 if (lex_match_id (s->lexer, "ELSE")
4064 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
4067 if (!lex_match_phrase (s->lexer, "END IF"))
4072 matrix_cmd_destroy (cmd);
4077 matrix_cmd_execute_do_if (struct do_if_command *cmd)
4079 for (size_t i = 0; i < cmd->n_clauses; i++)
4081 struct do_if_clause *c = &cmd->clauses[i];
4085 if (!matrix_expr_evaluate_scalar (c->condition,
4086 i ? "ELSE IF" : "DO IF",
4091 for (size_t j = 0; j < c->commands.n; j++)
4092 if (!matrix_cmd_execute (c->commands.commands[j]))
4099 static struct matrix_cmd *
4100 matrix_parse_loop (struct matrix_state *s)
4102 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4103 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4105 struct loop_command *loop = &cmd->loop;
4106 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4108 struct substring name = lex_tokss (s->lexer);
4109 loop->var = matrix_var_lookup (s, name);
4111 loop->var = matrix_var_create (s, name);
4116 loop->start = matrix_parse_expr (s);
4117 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4120 loop->end = matrix_parse_expr (s);
4124 if (lex_match (s->lexer, T_BY))
4126 loop->increment = matrix_parse_expr (s);
4127 if (!loop->increment)
4132 if (lex_match_id (s->lexer, "IF"))
4134 loop->top_condition = matrix_parse_expr (s);
4135 if (!loop->top_condition)
4139 bool was_in_loop = s->in_loop;
4141 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4143 s->in_loop = was_in_loop;
4147 if (!lex_match_phrase (s->lexer, "END LOOP"))
4150 if (lex_match_id (s->lexer, "IF"))
4152 loop->bottom_condition = matrix_parse_expr (s);
4153 if (!loop->bottom_condition)
4160 matrix_cmd_destroy (cmd);
4165 matrix_cmd_execute_loop (struct loop_command *cmd)
4169 long int increment = 1;
4172 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4173 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4175 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4179 if (increment > 0 ? value > end
4180 : increment < 0 ? value < end
4185 int mxloops = settings_get_mxloops ();
4186 for (int i = 0; i < mxloops; i++)
4191 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4193 gsl_matrix_free (cmd->var->value);
4194 cmd->var->value = NULL;
4196 if (!cmd->var->value)
4197 cmd->var->value = gsl_matrix_alloc (1, 1);
4199 gsl_matrix_set (cmd->var->value, 0, 0, value);
4202 if (cmd->top_condition)
4205 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4210 for (size_t j = 0; j < cmd->commands.n; j++)
4211 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4214 if (cmd->bottom_condition)
4217 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4218 "END LOOP IF", &d) || d > 0)
4225 if (increment > 0 ? value > end : value < end)
4231 static struct matrix_cmd *
4232 matrix_parse_break (struct matrix_state *s)
4236 msg (SE, _("BREAK not inside LOOP."));
4240 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4241 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4245 static struct matrix_cmd *
4246 matrix_parse_display (struct matrix_state *s)
4248 bool dictionary = false;
4249 bool status = false;
4250 while (lex_token (s->lexer) == T_ID)
4252 if (lex_match_id (s->lexer, "DICTIONARY"))
4254 else if (lex_match_id (s->lexer, "STATUS"))
4258 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4262 if (!dictionary && !status)
4263 dictionary = status = true;
4265 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4266 *cmd = (struct matrix_cmd) {
4267 .type = MCMD_DISPLAY,
4270 .dictionary = dictionary,
4278 compare_matrix_var_pointers (const void *a_, const void *b_)
4280 const struct matrix_var *const *ap = a_;
4281 const struct matrix_var *const *bp = b_;
4282 const struct matrix_var *a = *ap;
4283 const struct matrix_var *b = *bp;
4284 return strcmp (a->name, b->name);
4288 matrix_cmd_execute_display (struct display_command *cmd)
4290 const struct matrix_state *s = cmd->state;
4292 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4293 pivot_dimension_create (
4294 table, PIVOT_AXIS_COLUMN, N_("Property"),
4295 N_("Rows"), N_("Columns"), N_("Size (kB)"));
4297 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4300 const struct matrix_var *var;
4301 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4303 vars[n_vars++] = var;
4304 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4306 struct pivot_dimension *rows = pivot_dimension_create (
4307 table, PIVOT_AXIS_ROW, N_("Variable"));
4308 for (size_t i = 0; i < n_vars; i++)
4310 const struct matrix_var *var = vars[i];
4311 pivot_category_create_leaf (
4312 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4314 size_t r = var->value->size1;
4315 size_t c = var->value->size2;
4316 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4317 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4318 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4320 pivot_table_submit (table);
4323 static struct matrix_cmd *
4324 matrix_parse_release (struct matrix_state *s)
4326 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4327 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4329 size_t allocated_vars = 0;
4330 while (lex_token (s->lexer) == T_ID)
4332 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4335 if (cmd->release.n_vars >= allocated_vars)
4336 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4337 sizeof *cmd->release.vars);
4338 cmd->release.vars[cmd->release.n_vars++] = var;
4341 lex_error (s->lexer, _("Variable name expected."));
4344 if (!lex_match (s->lexer, T_COMMA))
4352 matrix_cmd_execute_release (struct release_command *cmd)
4354 for (size_t i = 0; i < cmd->n_vars; i++)
4356 struct matrix_var *var = cmd->vars[i];
4357 gsl_matrix_free (var->value);
4362 static struct save_file *
4363 save_file_create (struct matrix_state *s, struct file_handle *fh,
4364 struct string_array *variables,
4365 struct matrix_expr *names,
4366 struct stringi_set *strings)
4368 for (size_t i = 0; i < s->n_save_files; i++)
4370 struct save_file *sf = s->save_files[i];
4375 string_array_destroy (variables);
4376 matrix_expr_destroy (names);
4377 stringi_set_destroy (strings);
4383 struct save_file *sf = xmalloc (sizeof *sf);
4384 *sf = (struct save_file) {
4386 .variables = *variables,
4388 .strings = STRINGI_SET_INITIALIZER (sf->strings),
4391 stringi_set_swap (&sf->strings, strings);
4392 stringi_set_destroy (strings);
4394 s->save_files = xrealloc (s->save_files,
4395 (s->n_save_files + 1) * sizeof *s->save_files);
4396 s->save_files[s->n_save_files++] = sf;
4400 static struct casewriter *
4401 save_file_open (struct save_file *sf, gsl_matrix *m)
4403 if (sf->writer || sf->error)
4407 size_t n_variables = caseproto_get_n_widths (
4408 casewriter_get_proto (sf->writer));
4409 if (m->size2 != n_variables)
4411 msg (ME, _("The first SAVE to %s within this matrix program "
4412 "had %zu columns, so a %zu×%zu matrix cannot be "
4414 fh_get_name (sf->file), n_variables, m->size1, m->size2);
4422 struct dictionary *dict = dict_create (get_default_encoding ());
4424 /* Fill 'names' with user-specified names if there were any, either from
4425 sf->variables or sf->names. */
4426 struct string_array names = { .n = 0 };
4427 if (sf->variables.n)
4428 string_array_clone (&names, &sf->variables);
4431 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
4432 if (nm && is_vector (nm))
4434 gsl_vector nv = to_vector (nm);
4435 for (size_t i = 0; i < nv.size; i++)
4437 char *name = trimmed_string (gsl_vector_get (&nv, i));
4438 if (dict_id_is_valid (dict, name, true))
4439 string_array_append_nocopy (&names, name);
4444 gsl_matrix_free (nm);
4447 struct stringi_set strings;
4448 stringi_set_clone (&strings, &sf->strings);
4450 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4455 name = names.strings[i];
4458 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4462 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4463 struct variable *var = dict_create_var (dict, name, width);
4466 msg (ME, _("Duplicate variable name %s in SAVE statement."),
4472 if (!stringi_set_is_empty (&strings))
4474 const char *example = stringi_set_node_get_string (
4475 stringi_set_first (&strings));
4476 msg (ME, ngettext ("STRINGS specified a variable %s, but no variable "
4477 "with that name was found on SAVE.",
4478 "STRINGS specified %2$zu variables, including %1$s, "
4479 "whose names were not found on SAVE.",
4480 stringi_set_count (&strings)),
4481 example, stringi_set_count (&strings));
4484 stringi_set_destroy (&strings);
4485 string_array_destroy (&names);
4494 sf->writer = any_writer_open (sf->file, dict);
4505 save_file_destroy (struct save_file *sf)
4509 fh_unref (sf->file);
4510 string_array_destroy (&sf->variables);
4511 matrix_expr_destroy (sf->names);
4512 stringi_set_destroy (&sf->strings);
4513 casewriter_destroy (sf->writer);
4518 static struct matrix_cmd *
4519 matrix_parse_save (struct matrix_state *s)
4521 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4522 *cmd = (struct matrix_cmd) {
4524 .save = { .expression = NULL },
4527 struct file_handle *fh = NULL;
4528 struct save_command *save = &cmd->save;
4530 struct string_array variables = STRING_ARRAY_INITIALIZER;
4531 struct matrix_expr *names = NULL;
4532 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
4534 save->expression = matrix_parse_exp (s);
4535 if (!save->expression)
4538 while (lex_match (s->lexer, T_SLASH))
4540 if (lex_match_id (s->lexer, "OUTFILE"))
4542 lex_match (s->lexer, T_EQUALS);
4545 fh = (lex_match (s->lexer, T_ASTERISK)
4547 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4551 else if (lex_match_id (s->lexer, "VARIABLES"))
4553 lex_match (s->lexer, T_EQUALS);
4557 struct dictionary *d = dict_create (get_default_encoding ());
4558 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4559 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4564 string_array_clear (&variables);
4565 variables = (struct string_array) {
4571 else if (lex_match_id (s->lexer, "NAMES"))
4573 lex_match (s->lexer, T_EQUALS);
4574 matrix_expr_destroy (names);
4575 names = matrix_parse_exp (s);
4579 else if (lex_match_id (s->lexer, "STRINGS"))
4581 lex_match (s->lexer, T_EQUALS);
4582 while (lex_token (s->lexer) == T_ID)
4584 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
4586 if (!lex_match (s->lexer, T_COMMA))
4592 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4600 if (s->prev_save_file)
4601 fh = fh_ref (s->prev_save_file);
4604 lex_sbc_missing ("OUTFILE");
4608 fh_unref (s->prev_save_file);
4609 s->prev_save_file = fh_ref (fh);
4611 if (variables.n && names)
4613 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4614 matrix_expr_destroy (names);
4618 save->sf = save_file_create (s, fh, &variables, names, &strings);
4624 string_array_destroy (&variables);
4625 matrix_expr_destroy (names);
4626 stringi_set_destroy (&strings);
4628 matrix_cmd_destroy (cmd);
4633 matrix_cmd_execute_save (const struct save_command *save)
4635 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4639 struct casewriter *writer = save_file_open (save->sf, m);
4642 gsl_matrix_free (m);
4646 const struct caseproto *proto = casewriter_get_proto (writer);
4647 for (size_t y = 0; y < m->size1; y++)
4649 struct ccase *c = case_create (proto);
4650 for (size_t x = 0; x < m->size2; x++)
4652 double d = gsl_matrix_get (m, y, x);
4653 int width = caseproto_get_width (proto, x);
4654 union value *value = case_data_rw_idx (c, x);
4658 memcpy (value->s, &d, width);
4660 casewriter_write (writer, c);
4662 gsl_matrix_free (m);
4665 static struct matrix_cmd *
4666 matrix_parse_read (struct matrix_state *s)
4668 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4669 *cmd = (struct matrix_cmd) {
4671 .read = { .format = FMT_F },
4674 struct file_handle *fh = NULL;
4675 char *encoding = NULL;
4676 struct read_command *read = &cmd->read;
4677 read->dst = matrix_lvalue_parse (s);
4682 int repetitions = 0;
4683 int record_width = 0;
4684 bool seen_format = false;
4685 while (lex_match (s->lexer, T_SLASH))
4687 if (lex_match_id (s->lexer, "FILE"))
4689 lex_match (s->lexer, T_EQUALS);
4692 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4696 else if (lex_match_id (s->lexer, "ENCODING"))
4698 lex_match (s->lexer, T_EQUALS);
4699 if (!lex_force_string (s->lexer))
4703 encoding = ss_xstrdup (lex_tokss (s->lexer));
4707 else if (lex_match_id (s->lexer, "FIELD"))
4709 lex_match (s->lexer, T_EQUALS);
4711 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4713 read->c1 = lex_integer (s->lexer);
4715 if (!lex_force_match (s->lexer, T_TO)
4716 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4718 read->c2 = lex_integer (s->lexer) + 1;
4721 record_width = read->c2 - read->c1;
4722 if (lex_match (s->lexer, T_BY))
4724 if (!lex_force_int_range (s->lexer, "BY", 1,
4725 read->c2 - read->c1))
4727 by = lex_integer (s->lexer);
4730 if (record_width % by)
4732 msg (SE, _("BY %d does not evenly divide record width %d."),
4740 else if (lex_match_id (s->lexer, "SIZE"))
4742 lex_match (s->lexer, T_EQUALS);
4743 matrix_expr_destroy (read->size);
4744 read->size = matrix_parse_exp (s);
4748 else if (lex_match_id (s->lexer, "MODE"))
4750 lex_match (s->lexer, T_EQUALS);
4751 if (lex_match_id (s->lexer, "RECTANGULAR"))
4752 read->symmetric = false;
4753 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4754 read->symmetric = true;
4757 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4761 else if (lex_match_id (s->lexer, "REREAD"))
4762 read->reread = true;
4763 else if (lex_match_id (s->lexer, "FORMAT"))
4767 lex_sbc_only_once ("FORMAT");
4772 lex_match (s->lexer, T_EQUALS);
4774 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4777 const char *p = lex_tokcstr (s->lexer);
4778 if (c_isdigit (p[0]))
4780 repetitions = atoi (p);
4781 p += strspn (p, "0123456789");
4782 if (!fmt_from_name (p, &read->format))
4784 lex_error (s->lexer, _("Unknown format %s."), p);
4789 else if (fmt_from_name (p, &read->format))
4793 struct fmt_spec format;
4794 if (!parse_format_specifier (s->lexer, &format))
4796 read->format = format.type;
4802 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4803 "REREAD", "FORMAT");
4810 lex_sbc_missing ("FIELD");
4814 if (!read->dst->n_indexes && !read->size)
4816 msg (SE, _("SIZE is required for reading data into a full matrix "
4817 "(as opposed to a submatrix)."));
4823 if (s->prev_read_file)
4824 fh = fh_ref (s->prev_read_file);
4827 lex_sbc_missing ("FILE");
4831 fh_unref (s->prev_read_file);
4832 s->prev_read_file = fh_ref (fh);
4834 read->rf = read_file_create (s, fh);
4838 free (read->rf->encoding);
4839 read->rf->encoding = encoding;
4843 /* Field width may be specified in multiple ways:
4846 2. The format on FORMAT.
4847 3. The repetition factor on FORMAT.
4849 (2) and (3) are mutually exclusive.
4851 If more than one of these is present, they must agree. If none of them is
4852 present, then free-field format is used.
4854 if (repetitions > record_width)
4856 msg (SE, _("%d repetitions cannot fit in record width %d."),
4857 repetitions, record_width);
4860 int w = (repetitions ? record_width / repetitions
4866 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4867 "which implies field width %d, "
4868 "but BY specifies field width %d."),
4869 repetitions, record_width, w, by);
4871 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4880 matrix_cmd_destroy (cmd);
4886 parse_error (const struct dfm_reader *reader, enum fmt_type format,
4887 struct substring data, size_t y, size_t x,
4888 int first_column, int last_column, char *error)
4890 int line_number = dfm_get_line_number (reader);
4891 struct msg_location *location = xmalloc (sizeof *location);
4892 *location = (struct msg_location) {
4893 .file_name = xstrdup (dfm_get_file_name (reader)),
4894 .first_line = line_number,
4895 .last_line = line_number + 1,
4896 .first_column = first_column,
4897 .last_column = last_column,
4899 struct msg *m = xmalloc (sizeof *m);
4901 .category = MSG_C_DATA,
4902 .severity = MSG_S_WARNING,
4903 .location = location,
4904 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
4905 "for matrix row %zu, column %zu: %s"),
4906 (int) data.length, data.string, fmt_name (format),
4907 y + 1, x + 1, error),
4915 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4916 gsl_matrix *m, struct substring p, size_t y, size_t x,
4917 const char *line_start)
4919 const char *input_encoding = dfm_reader_get_encoding (reader);
4922 if (fmt_is_numeric (read->format))
4925 error = data_in (p, input_encoding, read->format,
4926 settings_get_fmt_settings (), &v, 0, NULL);
4927 if (!error && v.f == SYSMIS)
4928 error = xstrdup (_("Matrix data may not contain missing value."));
4933 uint8_t s[sizeof (double)];
4934 union value v = { .s = s };
4935 error = data_in (p, input_encoding, read->format,
4936 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
4937 memcpy (&f, s, sizeof f);
4942 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
4943 int c2 = c1 + ss_utf8_count_columns (p) - 1;
4944 parse_error (reader, read->format, p, y, x, c1, c2, error);
4948 gsl_matrix_set (m, y, x, f);
4949 if (read->symmetric && x != y)
4950 gsl_matrix_set (m, x, y, f);
4955 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4956 struct substring *line, const char **startp)
4958 if (dfm_eof (reader))
4960 msg (SE, _("Unexpected end of file reading matrix data."));
4963 dfm_expand_tabs (reader);
4964 struct substring record = dfm_get_record (reader);
4965 /* XXX need to recode record into UTF-8 */
4966 *startp = record.string;
4967 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
4972 matrix_read (struct read_command *read, struct dfm_reader *reader,
4975 for (size_t y = 0; y < m->size1; y++)
4977 size_t nx = read->symmetric ? y + 1 : m->size2;
4979 struct substring line = ss_empty ();
4980 const char *line_start = line.string;
4981 for (size_t x = 0; x < nx; x++)
4988 ss_ltrim (&line, ss_cstr (" ,"));
4989 if (!ss_is_empty (line))
4991 if (!matrix_read_line (read, reader, &line, &line_start))
4993 dfm_forward_record (reader);
4996 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
5000 if (!matrix_read_line (read, reader, &line, &line_start))
5002 size_t fields_per_line = (read->c2 - read->c1) / read->w;
5003 int f = x % fields_per_line;
5004 if (f == fields_per_line - 1)
5005 dfm_forward_record (reader);
5007 p = ss_substr (line, read->w * f, read->w);
5010 matrix_read_set_field (read, reader, m, p, y, x, line_start);
5014 dfm_forward_record (reader);
5017 ss_ltrim (&line, ss_cstr (" ,"));
5018 if (!ss_is_empty (line))
5021 msg (SW, _("Trailing garbage on line \"%.*s\""),
5022 (int) line.length, line.string);
5029 matrix_cmd_execute_read (struct read_command *read)
5031 struct index_vector iv0, iv1;
5032 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
5035 size_t size[2] = { SIZE_MAX, SIZE_MAX };
5038 gsl_matrix *m = matrix_expr_evaluate (read->size);
5044 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5045 "not a %zu×%zu matrix."), m->size1, m->size2);
5046 gsl_matrix_free (m);
5052 gsl_vector v = to_vector (m);
5056 d[0] = gsl_vector_get (&v, 0);
5059 else if (v.size == 2)
5061 d[0] = gsl_vector_get (&v, 0);
5062 d[1] = gsl_vector_get (&v, 1);
5066 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5067 "not a %zu×%zu matrix."),
5068 m->size1, m->size2),
5069 gsl_matrix_free (m);
5074 gsl_matrix_free (m);
5076 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
5078 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
5079 "are outside valid range."),
5090 if (read->dst->n_indexes)
5092 size_t submatrix_size[2];
5093 if (read->dst->n_indexes == 2)
5095 submatrix_size[0] = iv0.n;
5096 submatrix_size[1] = iv1.n;
5098 else if (read->dst->var->value->size1 == 1)
5100 submatrix_size[0] = 1;
5101 submatrix_size[1] = iv0.n;
5105 submatrix_size[0] = iv0.n;
5106 submatrix_size[1] = 1;
5111 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
5113 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
5114 "differ from submatrix dimensions %zu×%zu."),
5116 submatrix_size[0], submatrix_size[1]);
5124 size[0] = submatrix_size[0];
5125 size[1] = submatrix_size[1];
5129 struct dfm_reader *reader = read_file_open (read->rf);
5131 dfm_reread_record (reader, 1);
5133 if (read->symmetric && size[0] != size[1])
5135 msg (SE, _("Cannot read non-square %zu×%zu matrix "
5136 "using READ with MODE=SYMMETRIC."),
5142 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
5143 matrix_read (read, reader, tmp);
5144 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
5147 static struct matrix_cmd *
5148 matrix_parse_write (struct matrix_state *s)
5150 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5151 *cmd = (struct matrix_cmd) {
5155 struct file_handle *fh = NULL;
5156 char *encoding = NULL;
5157 struct write_command *write = &cmd->write;
5158 write->expression = matrix_parse_exp (s);
5159 if (!write->expression)
5163 int repetitions = 0;
5164 int record_width = 0;
5165 enum fmt_type format = FMT_F;
5166 bool has_format = false;
5167 while (lex_match (s->lexer, T_SLASH))
5169 if (lex_match_id (s->lexer, "OUTFILE"))
5171 lex_match (s->lexer, T_EQUALS);
5174 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5178 else if (lex_match_id (s->lexer, "ENCODING"))
5180 lex_match (s->lexer, T_EQUALS);
5181 if (!lex_force_string (s->lexer))
5185 encoding = ss_xstrdup (lex_tokss (s->lexer));
5189 else if (lex_match_id (s->lexer, "FIELD"))
5191 lex_match (s->lexer, T_EQUALS);
5193 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5195 write->c1 = lex_integer (s->lexer);
5197 if (!lex_force_match (s->lexer, T_TO)
5198 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5200 write->c2 = lex_integer (s->lexer) + 1;
5203 record_width = write->c2 - write->c1;
5204 if (lex_match (s->lexer, T_BY))
5206 if (!lex_force_int_range (s->lexer, "BY", 1,
5207 write->c2 - write->c1))
5209 by = lex_integer (s->lexer);
5212 if (record_width % by)
5214 msg (SE, _("BY %d does not evenly divide record width %d."),
5222 else if (lex_match_id (s->lexer, "MODE"))
5224 lex_match (s->lexer, T_EQUALS);
5225 if (lex_match_id (s->lexer, "RECTANGULAR"))
5226 write->triangular = false;
5227 else if (lex_match_id (s->lexer, "TRIANGULAR"))
5228 write->triangular = true;
5231 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5235 else if (lex_match_id (s->lexer, "HOLD"))
5237 else if (lex_match_id (s->lexer, "FORMAT"))
5239 if (has_format || write->format)
5241 lex_sbc_only_once ("FORMAT");
5245 lex_match (s->lexer, T_EQUALS);
5247 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5250 const char *p = lex_tokcstr (s->lexer);
5251 if (c_isdigit (p[0]))
5253 repetitions = atoi (p);
5254 p += strspn (p, "0123456789");
5255 if (!fmt_from_name (p, &format))
5257 lex_error (s->lexer, _("Unknown format %s."), p);
5263 else if (fmt_from_name (p, &format))
5270 struct fmt_spec spec;
5271 if (!parse_format_specifier (s->lexer, &spec))
5273 write->format = xmemdup (&spec, sizeof spec);
5278 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5286 lex_sbc_missing ("FIELD");
5292 if (s->prev_write_file)
5293 fh = fh_ref (s->prev_write_file);
5296 lex_sbc_missing ("OUTFILE");
5300 fh_unref (s->prev_write_file);
5301 s->prev_write_file = fh_ref (fh);
5303 write->wf = write_file_create (s, fh);
5307 free (write->wf->encoding);
5308 write->wf->encoding = encoding;
5312 /* Field width may be specified in multiple ways:
5315 2. The format on FORMAT.
5316 3. The repetition factor on FORMAT.
5318 (2) and (3) are mutually exclusive.
5320 If more than one of these is present, they must agree. If none of them is
5321 present, then free-field format is used.
5323 if (repetitions > record_width)
5325 msg (SE, _("%d repetitions cannot fit in record width %d."),
5326 repetitions, record_width);
5329 int w = (repetitions ? record_width / repetitions
5330 : write->format ? write->format->w
5335 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5336 "which implies field width %d, "
5337 "but BY specifies field width %d."),
5338 repetitions, record_width, w, by);
5340 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5344 if (w && !write->format)
5346 write->format = xmalloc (sizeof *write->format);
5347 *write->format = (struct fmt_spec) { .type = format, .w = w };
5349 if (!fmt_check_output (write->format))
5353 if (write->format && fmt_var_width (write->format) > sizeof (double))
5355 char s[FMT_STRING_LEN_MAX + 1];
5356 fmt_to_string (write->format, s);
5357 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
5358 s, sizeof (double));
5366 matrix_cmd_destroy (cmd);
5371 matrix_cmd_execute_write (struct write_command *write)
5373 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5377 if (write->triangular && m->size1 != m->size2)
5379 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5380 "the matrix to be written has dimensions %zu×%zu."),
5381 m->size1, m->size2);
5382 gsl_matrix_free (m);
5386 struct dfm_writer *writer = write_file_open (write->wf);
5387 if (!writer || !m->size1)
5389 gsl_matrix_free (m);
5393 const struct fmt_settings *settings = settings_get_fmt_settings ();
5394 struct u8_line *line = write->wf->held;
5395 for (size_t y = 0; y < m->size1; y++)
5399 line = xmalloc (sizeof *line);
5400 u8_line_init (line);
5402 size_t nx = write->triangular ? y + 1 : m->size2;
5404 for (size_t x = 0; x < nx; x++)
5407 double f = gsl_matrix_get (m, y, x);
5411 if (fmt_is_numeric (write->format->type))
5414 v.s = (uint8_t *) &f;
5415 s = data_out (&v, NULL, write->format, settings);
5419 s = xmalloc (DBL_BUFSIZE_BOUND);
5420 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
5421 >= DBL_BUFSIZE_BOUND)
5424 size_t len = strlen (s);
5425 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5426 if (width + x0 > write->c2)
5428 dfm_put_record_utf8 (writer, line->s.ss.string,
5430 u8_line_clear (line);
5433 u8_line_put (line, x0, x0 + width, s, len);
5436 x0 += write->format ? write->format->w : width + 1;
5439 if (y + 1 >= m->size1 && write->hold)
5441 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5442 u8_line_clear (line);
5446 u8_line_destroy (line);
5450 write->wf->held = line;
5452 gsl_matrix_free (m);
5455 static struct matrix_cmd *
5456 matrix_parse_get (struct matrix_state *s)
5458 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5459 *cmd = (struct matrix_cmd) {
5462 .dataset = s->dataset,
5463 .user = { .treatment = MGET_ERROR },
5464 .system = { .treatment = MGET_ERROR },
5468 struct get_command *get = &cmd->get;
5469 get->dst = matrix_lvalue_parse (s);
5473 while (lex_match (s->lexer, T_SLASH))
5475 if (lex_match_id (s->lexer, "FILE"))
5477 lex_match (s->lexer, T_EQUALS);
5479 fh_unref (get->file);
5480 if (lex_match (s->lexer, T_ASTERISK))
5484 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5489 else if (lex_match_id (s->lexer, "ENCODING"))
5491 lex_match (s->lexer, T_EQUALS);
5492 if (!lex_force_string (s->lexer))
5495 free (get->encoding);
5496 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5500 else if (lex_match_id (s->lexer, "VARIABLES"))
5502 lex_match (s->lexer, T_EQUALS);
5506 lex_sbc_only_once ("VARIABLES");
5510 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
5513 else if (lex_match_id (s->lexer, "NAMES"))
5515 lex_match (s->lexer, T_EQUALS);
5516 if (!lex_force_id (s->lexer))
5519 struct substring name = lex_tokss (s->lexer);
5520 get->names = matrix_var_lookup (s, name);
5522 get->names = matrix_var_create (s, name);
5525 else if (lex_match_id (s->lexer, "MISSING"))
5527 lex_match (s->lexer, T_EQUALS);
5528 if (lex_match_id (s->lexer, "ACCEPT"))
5529 get->user.treatment = MGET_ACCEPT;
5530 else if (lex_match_id (s->lexer, "OMIT"))
5531 get->user.treatment = MGET_OMIT;
5532 else if (lex_is_number (s->lexer))
5534 get->user.treatment = MGET_RECODE;
5535 get->user.substitute = lex_number (s->lexer);
5540 lex_error (s->lexer, NULL);
5544 else if (lex_match_id (s->lexer, "SYSMIS"))
5546 lex_match (s->lexer, T_EQUALS);
5547 if (lex_match_id (s->lexer, "OMIT"))
5548 get->system.treatment = MGET_OMIT;
5549 else if (lex_is_number (s->lexer))
5551 get->system.treatment = MGET_RECODE;
5552 get->system.substitute = lex_number (s->lexer);
5557 lex_error (s->lexer, NULL);
5563 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5564 "MISSING", "SYSMIS");
5569 if (get->user.treatment != MGET_ACCEPT)
5570 get->system.treatment = MGET_ERROR;
5575 matrix_cmd_destroy (cmd);
5580 matrix_cmd_execute_get__ (struct get_command *get,
5581 const struct dictionary *dict,
5582 struct casereader *reader)
5584 struct variable **vars;
5589 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
5590 &vars, &n_vars, PV_NUMERIC))
5595 n_vars = dict_get_var_cnt (dict);
5596 vars = xnmalloc (n_vars, sizeof *vars);
5597 for (size_t i = 0; i < n_vars; i++)
5599 struct variable *var = dict_get_var (dict, i);
5600 if (!var_is_numeric (var))
5602 msg (SE, _("GET: Variable %s is not numeric."),
5603 var_get_name (var));
5613 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
5614 for (size_t i = 0; i < n_vars; i++)
5616 char s[sizeof (double)];
5618 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
5619 memcpy (&f, s, sizeof f);
5620 gsl_matrix_set (names, i, 0, f);
5623 gsl_matrix_free (get->names->value);
5624 get->names->value = names;
5628 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5629 long long int casenum = 1;
5631 for (struct ccase *c = casereader_read (reader); c;
5632 c = casereader_read (reader), casenum++)
5634 if (n_rows >= m->size1)
5636 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5637 for (size_t y = 0; y < n_rows; y++)
5638 for (size_t x = 0; x < n_vars; x++)
5639 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5640 gsl_matrix_free (m);
5645 for (size_t x = 0; x < n_vars; x++)
5647 const struct variable *var = vars[x];
5648 double d = case_num (c, var);
5651 if (get->system.treatment == MGET_RECODE)
5652 d = get->system.substitute;
5653 else if (get->system.treatment == MGET_OMIT)
5657 msg (SE, _("GET: Variable %s in case %lld "
5658 "is system-missing."),
5659 var_get_name (var), casenum);
5663 else if (var_is_num_missing (var, d, MV_USER))
5665 if (get->user.treatment == MGET_RECODE)
5666 d = get->user.substitute;
5667 else if (get->user.treatment == MGET_OMIT)
5669 else if (get->user.treatment != MGET_ACCEPT)
5671 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5673 var_get_name (var), casenum, d);
5677 gsl_matrix_set (m, n_rows, x, d);
5688 matrix_lvalue_evaluate_and_assign (get->dst, m);
5691 gsl_matrix_free (m);
5696 matrix_cmd_execute_get (struct get_command *get)
5698 struct dictionary *dict;
5699 struct casereader *reader;
5702 reader = any_reader_open_and_decode (get->file, get->encoding,
5709 if (dict_get_var_cnt (dataset_dict (get->dataset)) == 0)
5711 msg (ME, _("GET cannot read empty active file."));
5714 reader = proc_open (get->dataset);
5715 dict = dict_ref (dataset_dict (get->dataset));
5718 matrix_cmd_execute_get__ (get, dict, reader);
5721 casereader_destroy (reader);
5723 proc_commit (get->dataset);
5727 match_rowtype (struct lexer *lexer)
5729 static const char *rowtypes[] = {
5730 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5732 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5734 for (size_t i = 0; i < n_rowtypes; i++)
5735 if (lex_match_id (lexer, rowtypes[i]))
5738 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5743 parse_var_names (struct lexer *lexer, struct string_array *sa)
5745 lex_match (lexer, T_EQUALS);
5747 string_array_clear (sa);
5749 struct dictionary *dict = dict_create (get_default_encoding ());
5750 dict_create_var_assert (dict, "ROWTYPE_", 8);
5751 dict_create_var_assert (dict, "VARNAME_", 8);
5754 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5755 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5760 string_array_clear (sa);
5761 sa->strings = names;
5762 sa->n = sa->allocated = n_names;
5768 msave_common_uninit (struct msave_common *common)
5772 fh_unref (common->outfile);
5773 string_array_destroy (&common->variables);
5774 string_array_destroy (&common->fnames);
5775 string_array_destroy (&common->snames);
5780 compare_variables (const char *keyword,
5781 const struct string_array *new,
5782 const struct string_array *old)
5788 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5789 "on the first MSAVE within MATRIX."), keyword);
5792 else if (!string_array_equal_case (old, new))
5794 msg (SE, _("%s must specify the same variables each time within "
5795 "a given MATRIX."), keyword);
5802 static struct matrix_cmd *
5803 matrix_parse_msave (struct matrix_state *s)
5805 struct msave_common common = { .outfile = NULL };
5806 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5807 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5809 struct msave_command *msave = &cmd->msave;
5810 if (lex_token (s->lexer) == T_ID)
5811 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5812 msave->expr = matrix_parse_exp (s);
5816 while (lex_match (s->lexer, T_SLASH))
5818 if (lex_match_id (s->lexer, "TYPE"))
5820 lex_match (s->lexer, T_EQUALS);
5822 msave->rowtype = match_rowtype (s->lexer);
5823 if (!msave->rowtype)
5826 else if (lex_match_id (s->lexer, "OUTFILE"))
5828 lex_match (s->lexer, T_EQUALS);
5830 fh_unref (common.outfile);
5831 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5832 if (!common.outfile)
5835 else if (lex_match_id (s->lexer, "VARIABLES"))
5837 if (!parse_var_names (s->lexer, &common.variables))
5840 else if (lex_match_id (s->lexer, "FNAMES"))
5842 if (!parse_var_names (s->lexer, &common.fnames))
5845 else if (lex_match_id (s->lexer, "SNAMES"))
5847 if (!parse_var_names (s->lexer, &common.snames))
5850 else if (lex_match_id (s->lexer, "SPLIT"))
5852 lex_match (s->lexer, T_EQUALS);
5854 matrix_expr_destroy (msave->splits);
5855 msave->splits = matrix_parse_exp (s);
5859 else if (lex_match_id (s->lexer, "FACTOR"))
5861 lex_match (s->lexer, T_EQUALS);
5863 matrix_expr_destroy (msave->factors);
5864 msave->factors = matrix_parse_exp (s);
5865 if (!msave->factors)
5870 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5871 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5875 if (!msave->rowtype)
5877 lex_sbc_missing ("TYPE");
5880 common.has_splits = msave->splits || common.snames.n;
5881 common.has_factors = msave->factors || common.fnames.n;
5883 struct msave_common *c = s->common ? s->common : &common;
5884 if (c->fnames.n && !msave->factors)
5886 msg (SE, _("FNAMES requires FACTOR."));
5889 if (c->snames.n && !msave->splits)
5891 msg (SE, _("SNAMES requires SPLIT."));
5894 if (c->has_factors && !common.has_factors)
5896 msg (SE, _("%s is required because it was present on the first "
5897 "MSAVE in this MATRIX command."), "FACTOR");
5900 if (c->has_splits && !common.has_splits)
5902 msg (SE, _("%s is required because it was present on the first "
5903 "MSAVE in this MATRIX command."), "SPLIT");
5909 if (!common.outfile)
5911 lex_sbc_missing ("OUTFILE");
5914 s->common = xmemdup (&common, sizeof common);
5920 bool same = common.outfile == s->common->outfile;
5921 fh_unref (common.outfile);
5924 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5925 "within a single MATRIX command."));
5929 if (!compare_variables ("VARIABLES",
5930 &common.variables, &s->common->variables)
5931 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5932 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5934 msave_common_uninit (&common);
5936 msave->common = s->common;
5937 if (!msave->varname_)
5938 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5942 msave_common_uninit (&common);
5943 matrix_cmd_destroy (cmd);
5948 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5950 gsl_matrix *m = matrix_expr_evaluate (e);
5956 msg (SE, _("%s expression must evaluate to vector, "
5957 "not a %zu×%zu matrix."),
5958 name, m->size1, m->size2);
5959 gsl_matrix_free (m);
5963 return matrix_to_vector (m);
5967 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5969 for (size_t i = 0; i < vars->n; i++)
5970 if (!dict_create_var (d, vars->strings[i], 0))
5971 return vars->strings[i];
5975 static struct dictionary *
5976 msave_create_dict (const struct msave_common *common)
5978 struct dictionary *dict = dict_create (get_default_encoding ());
5980 const char *dup_split = msave_add_vars (dict, &common->snames);
5983 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5987 dict_create_var_assert (dict, "ROWTYPE_", 8);
5989 const char *dup_factor = msave_add_vars (dict, &common->fnames);
5992 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5996 dict_create_var_assert (dict, "VARNAME_", 8);
5998 const char *dup_var = msave_add_vars (dict, &common->variables);
6001 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
6013 matrix_cmd_execute_msave (struct msave_command *msave)
6015 struct msave_common *common = msave->common;
6016 gsl_matrix *m = NULL;
6017 gsl_vector *factors = NULL;
6018 gsl_vector *splits = NULL;
6020 m = matrix_expr_evaluate (msave->expr);
6024 if (!common->variables.n)
6025 for (size_t i = 0; i < m->size2; i++)
6026 string_array_append_nocopy (&common->variables,
6027 xasprintf ("COL%zu", i + 1));
6029 if (m->size2 != common->variables.n)
6031 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
6032 m->size2, common->variables.n);
6038 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
6042 if (!common->fnames.n)
6043 for (size_t i = 0; i < factors->size; i++)
6044 string_array_append_nocopy (&common->fnames,
6045 xasprintf ("FAC%zu", i + 1));
6047 if (factors->size != common->fnames.n)
6049 msg (SE, _("There are %zu factor variables, "
6050 "but %zu split values were supplied."),
6051 common->fnames.n, factors->size);
6058 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
6062 if (!common->fnames.n)
6063 for (size_t i = 0; i < splits->size; i++)
6064 string_array_append_nocopy (&common->fnames,
6065 xasprintf ("SPL%zu", i + 1));
6067 if (splits->size != common->fnames.n)
6069 msg (SE, _("There are %zu split variables, "
6070 "but %zu split values were supplied."),
6071 common->fnames.n, splits->size);
6076 if (!common->writer)
6078 struct dictionary *dict = msave_create_dict (common);
6082 common->writer = any_writer_open (common->outfile, dict);
6083 if (!common->writer)
6089 common->dict = dict;
6092 for (size_t y = 0; y < m->size1; y++)
6094 struct ccase *c = case_create (dict_get_proto (common->dict));
6097 /* Split variables */
6099 for (size_t i = 0; i < splits->size; i++)
6100 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
6103 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6104 msave->rowtype, ' ');
6108 for (size_t i = 0; i < factors->size; i++)
6109 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
6112 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6113 msave->varname_, ' ');
6115 /* Continuous variables. */
6116 for (size_t x = 0; x < m->size2; x++)
6117 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
6118 casewriter_write (common->writer, c);
6122 gsl_matrix_free (m);
6123 gsl_vector_free (factors);
6124 gsl_vector_free (splits);
6127 static struct matrix_cmd *
6128 matrix_parse_mget (struct matrix_state *s)
6130 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6131 *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
6133 struct mget_command *mget = &cmd->mget;
6137 if (lex_match_id (s->lexer, "FILE"))
6139 lex_match (s->lexer, T_EQUALS);
6141 fh_unref (mget->file);
6142 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6146 else if (lex_match_id (s->lexer, "TYPE"))
6148 lex_match (s->lexer, T_EQUALS);
6149 while (lex_token (s->lexer) != T_SLASH
6150 && lex_token (s->lexer) != T_ENDCMD)
6152 const char *rowtype = match_rowtype (s->lexer);
6156 stringi_set_insert (&mget->rowtypes, rowtype);
6161 lex_error_expecting (s->lexer, "FILE", "TYPE");
6164 if (lex_token (s->lexer) == T_ENDCMD)
6167 if (!lex_force_match (s->lexer, T_SLASH))
6173 matrix_cmd_destroy (cmd);
6177 static const struct variable *
6178 get_a8_var (const struct dictionary *d, const char *name)
6180 const struct variable *v = dict_lookup_var (d, name);
6183 msg (SE, _("Matrix data file lacks %s variable."), name);
6186 if (var_get_width (v) != 8)
6188 msg (SE, _("%s variable in matrix data file must be "
6189 "8-byte string, but it has width %d."),
6190 name, var_get_width (v));
6197 is_rowtype (const union value *v, const char *rowtype)
6199 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
6200 ss_rtrim (&vs, ss_cstr (" "));
6201 return ss_equals_case (vs, ss_cstr (rowtype));
6205 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6206 const struct dictionary *d,
6207 const struct variable *rowtype_var,
6208 struct matrix_state *s, size_t si, size_t fi,
6209 size_t cs, size_t cn)
6214 const union value *rowtype_ = case_data (rows[0], rowtype_var);
6215 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
6216 : is_rowtype (rowtype_, "CORR") ? "CR"
6217 : is_rowtype (rowtype_, "MEAN") ? "MN"
6218 : is_rowtype (rowtype_, "STDDEV") ? "SD"
6219 : is_rowtype (rowtype_, "N") ? "NC"
6222 struct string name = DS_EMPTY_INITIALIZER;
6223 ds_put_cstr (&name, name_prefix);
6225 ds_put_format (&name, "F%zu", fi);
6227 ds_put_format (&name, "S%zu", si);
6229 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6231 mv = matrix_var_create (s, ds_ss (&name));
6234 msg (SW, _("Matrix data file contains variable with existing name %s."),
6239 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6240 size_t n_missing = 0;
6241 for (size_t y = 0; y < n_rows; y++)
6243 for (size_t x = 0; x < cn; x++)
6245 struct variable *var = dict_get_var (d, cs + x);
6246 double value = case_num (rows[y], var);
6247 if (var_is_num_missing (var, value, MV_ANY))
6252 gsl_matrix_set (m, y, x, value);
6257 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6258 "value, which was treated as zero.",
6259 "Matrix data file variable %s contains %zu missing "
6260 "values, which were treated as zero.", n_missing),
6261 ds_cstr (&name), n_missing);
6266 for (size_t y = 0; y < n_rows; y++)
6267 case_unref (rows[y]);
6271 var_changed (const struct ccase *ca, const struct ccase *cb,
6272 const struct variable *var)
6274 return !value_equal (case_data (ca, var), case_data (cb, var),
6275 var_get_width (var));
6279 vars_changed (const struct ccase *ca, const struct ccase *cb,
6280 const struct dictionary *d,
6281 size_t first_var, size_t n_vars)
6283 for (size_t i = 0; i < n_vars; i++)
6285 const struct variable *v = dict_get_var (d, first_var + i);
6286 if (var_changed (ca, cb, v))
6293 matrix_cmd_execute_mget (struct mget_command *mget)
6295 struct dictionary *d;
6296 struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6301 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6302 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6303 if (!rowtype_ || !varname_)
6306 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6308 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6311 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6313 msg (SE, _("Matrix data file contains no continuous variables."));
6317 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6319 const struct variable *v = dict_get_var (d, i);
6320 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6323 _("Matrix data file contains unexpected string variable %s."),
6329 /* SPLIT variables. */
6331 size_t sn = var_get_dict_index (rowtype_);
6332 struct ccase *sc = NULL;
6335 /* FACTOR variables. */
6336 size_t fs = var_get_dict_index (rowtype_) + 1;
6337 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6338 struct ccase *fc = NULL;
6341 /* Continuous variables. */
6342 size_t cs = var_get_dict_index (varname_) + 1;
6343 size_t cn = dict_get_var_cnt (d) - cs;
6344 struct ccase *cc = NULL;
6347 struct ccase **rows = NULL;
6348 size_t allocated_rows = 0;
6352 while ((c = casereader_read (r)) != NULL)
6354 bool sd = vars_changed (sc, c, d, ss, sn);
6355 bool fd = sd || vars_changed (fc, c, d, fs, fn);
6356 bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6371 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6372 mget->state, si, fi, cs, cn);
6378 if (n_rows >= allocated_rows)
6379 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6382 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6383 mget->state, si, fi, cs, cn);
6387 casereader_destroy (r);
6391 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6393 if (!lex_force_id (s->lexer))
6396 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6398 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6403 static struct matrix_cmd *
6404 matrix_parse_eigen (struct matrix_state *s)
6406 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6407 *cmd = (struct matrix_cmd) {
6409 .eigen = { .expr = NULL }
6412 struct eigen_command *eigen = &cmd->eigen;
6413 if (!lex_force_match (s->lexer, T_LPAREN))
6415 eigen->expr = matrix_parse_expr (s);
6417 || !lex_force_match (s->lexer, T_COMMA)
6418 || !matrix_parse_dst_var (s, &eigen->evec)
6419 || !lex_force_match (s->lexer, T_COMMA)
6420 || !matrix_parse_dst_var (s, &eigen->eval)
6421 || !lex_force_match (s->lexer, T_RPAREN))
6427 matrix_cmd_destroy (cmd);
6432 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6434 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6437 if (!is_symmetric (A))
6439 msg (SE, _("Argument of EIGEN must be symmetric."));
6440 gsl_matrix_free (A);
6444 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6445 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6446 gsl_vector v_eval = to_vector (eval);
6447 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6448 gsl_eigen_symmv (A, &v_eval, evec, w);
6449 gsl_eigen_symmv_free (w);
6451 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6453 gsl_matrix_free (A);
6455 gsl_matrix_free (eigen->eval->value);
6456 eigen->eval->value = eval;
6458 gsl_matrix_free (eigen->evec->value);
6459 eigen->evec->value = evec;
6462 static struct matrix_cmd *
6463 matrix_parse_setdiag (struct matrix_state *s)
6465 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6466 *cmd = (struct matrix_cmd) {
6467 .type = MCMD_SETDIAG,
6468 .setdiag = { .dst = NULL }
6471 struct setdiag_command *setdiag = &cmd->setdiag;
6472 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6475 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6478 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6483 if (!lex_force_match (s->lexer, T_COMMA))
6486 setdiag->expr = matrix_parse_expr (s);
6490 if (!lex_force_match (s->lexer, T_RPAREN))
6496 matrix_cmd_destroy (cmd);
6501 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6503 gsl_matrix *dst = setdiag->dst->value;
6506 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6507 setdiag->dst->name);
6511 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6515 size_t n = MIN (dst->size1, dst->size2);
6516 if (is_scalar (src))
6518 double d = to_scalar (src);
6519 for (size_t i = 0; i < n; i++)
6520 gsl_matrix_set (dst, i, i, d);
6522 else if (is_vector (src))
6524 gsl_vector v = to_vector (src);
6525 for (size_t i = 0; i < n && i < v.size; i++)
6526 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6529 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6530 "not a %zu×%zu matrix."),
6531 src->size1, src->size2);
6532 gsl_matrix_free (src);
6535 static struct matrix_cmd *
6536 matrix_parse_svd (struct matrix_state *s)
6538 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6539 *cmd = (struct matrix_cmd) {
6541 .svd = { .expr = NULL }
6544 struct svd_command *svd = &cmd->svd;
6545 if (!lex_force_match (s->lexer, T_LPAREN))
6547 svd->expr = matrix_parse_expr (s);
6549 || !lex_force_match (s->lexer, T_COMMA)
6550 || !matrix_parse_dst_var (s, &svd->u)
6551 || !lex_force_match (s->lexer, T_COMMA)
6552 || !matrix_parse_dst_var (s, &svd->s)
6553 || !lex_force_match (s->lexer, T_COMMA)
6554 || !matrix_parse_dst_var (s, &svd->v)
6555 || !lex_force_match (s->lexer, T_RPAREN))
6561 matrix_cmd_destroy (cmd);
6566 matrix_cmd_execute_svd (struct svd_command *svd)
6568 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6572 if (m->size1 >= m->size2)
6575 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6576 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6577 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6578 gsl_vector *work = gsl_vector_alloc (A->size2);
6579 gsl_linalg_SV_decomp (A, V, &Sv, work);
6580 gsl_vector_free (work);
6582 matrix_var_set (svd->u, A);
6583 matrix_var_set (svd->s, S);
6584 matrix_var_set (svd->v, V);
6588 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6589 gsl_matrix_transpose_memcpy (At, m);
6590 gsl_matrix_free (m);
6592 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6593 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6594 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6595 gsl_vector *work = gsl_vector_alloc (At->size2);
6596 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6597 gsl_vector_free (work);
6599 matrix_var_set (svd->v, At);
6600 matrix_var_set (svd->s, St);
6601 matrix_var_set (svd->u, Vt);
6606 matrix_cmd_execute (struct matrix_cmd *cmd)
6611 matrix_cmd_execute_compute (&cmd->compute);
6615 matrix_cmd_execute_print (&cmd->print);
6619 return matrix_cmd_execute_do_if (&cmd->do_if);
6622 matrix_cmd_execute_loop (&cmd->loop);
6629 matrix_cmd_execute_display (&cmd->display);
6633 matrix_cmd_execute_release (&cmd->release);
6637 matrix_cmd_execute_save (&cmd->save);
6641 matrix_cmd_execute_read (&cmd->read);
6645 matrix_cmd_execute_write (&cmd->write);
6649 matrix_cmd_execute_get (&cmd->get);
6653 matrix_cmd_execute_msave (&cmd->msave);
6657 matrix_cmd_execute_mget (&cmd->mget);
6661 matrix_cmd_execute_eigen (&cmd->eigen);
6665 matrix_cmd_execute_setdiag (&cmd->setdiag);
6669 matrix_cmd_execute_svd (&cmd->svd);
6677 matrix_cmds_uninit (struct matrix_cmds *cmds)
6679 for (size_t i = 0; i < cmds->n; i++)
6680 matrix_cmd_destroy (cmds->commands[i]);
6681 free (cmds->commands);
6685 matrix_cmd_destroy (struct matrix_cmd *cmd)
6693 matrix_lvalue_destroy (cmd->compute.lvalue);
6694 matrix_expr_destroy (cmd->compute.rvalue);
6698 matrix_expr_destroy (cmd->print.expression);
6699 free (cmd->print.title);
6700 print_labels_destroy (cmd->print.rlabels);
6701 print_labels_destroy (cmd->print.clabels);
6705 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6707 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6708 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6710 free (cmd->do_if.clauses);
6714 matrix_expr_destroy (cmd->loop.start);
6715 matrix_expr_destroy (cmd->loop.end);
6716 matrix_expr_destroy (cmd->loop.increment);
6717 matrix_expr_destroy (cmd->loop.top_condition);
6718 matrix_expr_destroy (cmd->loop.bottom_condition);
6719 matrix_cmds_uninit (&cmd->loop.commands);
6729 free (cmd->release.vars);
6733 matrix_expr_destroy (cmd->save.expression);
6737 matrix_lvalue_destroy (cmd->read.dst);
6738 matrix_expr_destroy (cmd->read.size);
6742 matrix_expr_destroy (cmd->write.expression);
6743 free (cmd->write.format);
6747 matrix_lvalue_destroy (cmd->get.dst);
6748 fh_unref (cmd->get.file);
6749 free (cmd->get.encoding);
6750 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
6754 free (cmd->msave.varname_);
6755 matrix_expr_destroy (cmd->msave.expr);
6756 matrix_expr_destroy (cmd->msave.factors);
6757 matrix_expr_destroy (cmd->msave.splits);
6761 fh_unref (cmd->mget.file);
6762 stringi_set_destroy (&cmd->mget.rowtypes);
6766 matrix_expr_destroy (cmd->eigen.expr);
6770 matrix_expr_destroy (cmd->setdiag.expr);
6774 matrix_expr_destroy (cmd->svd.expr);
6780 struct matrix_command_name
6783 struct matrix_cmd *(*parse) (struct matrix_state *);
6786 static const struct matrix_command_name *
6787 matrix_parse_command_name (struct lexer *lexer)
6789 static const struct matrix_command_name commands[] = {
6790 { "COMPUTE", matrix_parse_compute },
6791 { "CALL EIGEN", matrix_parse_eigen },
6792 { "CALL SETDIAG", matrix_parse_setdiag },
6793 { "CALL SVD", matrix_parse_svd },
6794 { "PRINT", matrix_parse_print },
6795 { "DO IF", matrix_parse_do_if },
6796 { "LOOP", matrix_parse_loop },
6797 { "BREAK", matrix_parse_break },
6798 { "READ", matrix_parse_read },
6799 { "WRITE", matrix_parse_write },
6800 { "GET", matrix_parse_get },
6801 { "SAVE", matrix_parse_save },
6802 { "MGET", matrix_parse_mget },
6803 { "MSAVE", matrix_parse_msave },
6804 { "DISPLAY", matrix_parse_display },
6805 { "RELEASE", matrix_parse_release },
6807 static size_t n = sizeof commands / sizeof *commands;
6809 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6810 if (lex_match_phrase (lexer, c->name))
6815 static struct matrix_cmd *
6816 matrix_parse_command (struct matrix_state *s)
6818 size_t nesting_level = SIZE_MAX;
6820 struct matrix_cmd *c = NULL;
6821 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6823 lex_error (s->lexer, _("Unknown matrix command."));
6824 else if (!cmd->parse)
6825 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6829 nesting_level = output_open_group (
6830 group_item_create_nocopy (utf8_to_title (cmd->name),
6831 utf8_to_title (cmd->name)));
6836 lex_end_of_command (s->lexer);
6837 lex_discard_rest_of_command (s->lexer);
6838 if (nesting_level != SIZE_MAX)
6839 output_close_groups (nesting_level);
6845 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6847 if (!lex_force_match (lexer, T_ENDCMD))
6850 struct matrix_state state = {
6852 .session = dataset_session (ds),
6854 .vars = HMAP_INITIALIZER (state.vars),
6859 while (lex_match (lexer, T_ENDCMD))
6861 if (lex_token (lexer) == T_STOP)
6863 msg (SE, _("Unexpected end of input expecting matrix command."));
6867 if (lex_match_phrase (lexer, "END MATRIX"))
6870 struct matrix_cmd *c = matrix_parse_command (&state);
6873 matrix_cmd_execute (c);
6874 matrix_cmd_destroy (c);
6878 struct matrix_var *var, *next;
6879 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6882 gsl_matrix_free (var->value);
6883 hmap_delete (&state.vars, &var->hmap_node);
6886 hmap_destroy (&state.vars);
6889 dict_unref (state.common->dict);
6890 casewriter_destroy (state.common->writer);
6891 free (state.common);
6893 fh_unref (state.prev_read_file);
6894 for (size_t i = 0; i < state.n_read_files; i++)
6895 read_file_destroy (state.read_files[i]);
6896 free (state.read_files);
6897 fh_unref (state.prev_write_file);
6898 for (size_t i = 0; i < state.n_write_files; i++)
6899 write_file_destroy (state.write_files[i]);
6900 free (state.write_files);
6901 fh_unref (state.prev_save_file);
6902 for (size_t i = 0; i < state.n_save_files; i++)
6903 save_file_destroy (state.save_files[i]);
6904 free (state.save_files);