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;
111 struct dataset *dataset;
113 /* Parameters from parsing the first SAVE command for 'file'. */
114 struct string_array variables;
115 struct matrix_expr *names;
116 struct stringi_set strings;
118 /* Results from the first (only) attempt to open this save_file. */
120 struct casewriter *writer;
121 struct dictionary *dict;
126 struct dataset *dataset;
127 struct session *session;
131 struct msave_common *common;
133 struct file_handle *prev_read_file;
134 struct read_file **read_files;
137 struct file_handle *prev_write_file;
138 struct write_file **write_files;
139 size_t n_write_files;
141 struct file_handle *prev_save_file;
142 struct save_file **save_files;
146 static struct matrix_var *
147 matrix_var_lookup (struct matrix_state *s, struct substring name)
149 struct matrix_var *var;
151 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
152 utf8_hash_case_substring (name, 0), &s->vars)
153 if (!utf8_sscasecmp (ss_cstr (var->name), name))
159 static struct matrix_var *
160 matrix_var_create (struct matrix_state *s, struct substring name)
162 struct matrix_var *var = xmalloc (sizeof *var);
163 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
164 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
169 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
171 gsl_matrix_free (var->value);
175 #define MATRIX_FUNCTIONS \
231 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
232 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
233 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
234 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
235 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
236 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
237 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
238 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
239 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
240 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
241 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
243 typedef gsl_matrix *matrix_proto_m_d (double);
244 typedef gsl_matrix *matrix_proto_m_dd (double, double);
245 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
246 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
247 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
248 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
249 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
250 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
251 typedef double matrix_proto_d_m (gsl_matrix *);
252 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
253 typedef gsl_matrix *matrix_proto_IDENT (double, double);
255 #define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
259 /* Matrix expressions. */
266 #define F(NAME, PROTOTYPE) MOP_F_##NAME,
270 /* Elementwise and scalar arithmetic. */
271 MOP_NEGATE, /* unary - */
272 MOP_ADD_ELEMS, /* + */
273 MOP_SUB_ELEMS, /* - */
274 MOP_MUL_ELEMS, /* &* */
275 MOP_DIV_ELEMS, /* / and &/ */
276 MOP_EXP_ELEMS, /* &** */
278 MOP_SEQ_BY, /* a:b:c */
280 /* Matrix arithmetic. */
282 MOP_EXP_MAT, /* ** */
299 MOP_PASTE_HORZ, /* a, b, c, ... */
300 MOP_PASTE_VERT, /* a; b; c; ... */
304 MOP_VEC_INDEX, /* x(y) */
305 MOP_VEC_ALL, /* x(:) */
306 MOP_MAT_INDEX, /* x(y,z) */
307 MOP_ROW_INDEX, /* x(y,:) */
308 MOP_COL_INDEX, /* x(:,z) */
315 MOP_EOF, /* EOF('file') */
323 struct matrix_expr **subs;
328 struct matrix_var *variable;
329 struct read_file *eof;
334 matrix_expr_destroy (struct matrix_expr *e)
341 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
372 for (size_t i = 0; i < e->n_subs; i++)
373 matrix_expr_destroy (e->subs[i]);
385 static struct matrix_expr *
386 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
389 struct matrix_expr *e = xmalloc (sizeof *e);
390 *e = (struct matrix_expr) {
392 .subs = xmemdup (subs, n_subs * sizeof *subs),
398 static struct matrix_expr *
399 matrix_expr_create_0 (enum matrix_op op)
401 struct matrix_expr *sub;
402 return matrix_expr_create_subs (op, &sub, 0);
405 static struct matrix_expr *
406 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
408 return matrix_expr_create_subs (op, &sub, 1);
411 static struct matrix_expr *
412 matrix_expr_create_2 (enum matrix_op op,
413 struct matrix_expr *sub0, struct matrix_expr *sub1)
415 struct matrix_expr *subs[] = { sub0, sub1 };
416 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
419 static struct matrix_expr *
420 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
421 struct matrix_expr *sub1, struct matrix_expr *sub2)
423 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
424 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
427 static struct matrix_expr *
428 matrix_expr_create_number (double number)
430 struct matrix_expr *e = xmalloc (sizeof *e);
431 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
435 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
437 static struct matrix_expr *
438 matrix_parse_curly_comma (struct matrix_state *s)
440 struct matrix_expr *lhs = matrix_parse_or_xor (s);
444 while (lex_match (s->lexer, T_COMMA))
446 struct matrix_expr *rhs = matrix_parse_or_xor (s);
449 matrix_expr_destroy (lhs);
452 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
457 static struct matrix_expr *
458 matrix_parse_curly_semi (struct matrix_state *s)
460 if (lex_token (s->lexer) == T_RCURLY)
461 return matrix_expr_create_0 (MOP_EMPTY);
463 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
467 while (lex_match (s->lexer, T_SEMICOLON))
469 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
472 matrix_expr_destroy (lhs);
475 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
480 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
481 for (size_t Y = 0; Y < (M)->size1; Y++) \
482 for (size_t X = 0; X < (M)->size2; X++) \
483 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
486 is_vector (const gsl_matrix *m)
488 return m->size1 <= 1 || m->size2 <= 1;
492 to_vector (gsl_matrix *m)
494 return (m->size1 == 1
495 ? gsl_matrix_row (m, 0).vector
496 : gsl_matrix_column (m, 0).vector);
501 matrix_eval_ABS (gsl_matrix *m)
503 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
509 matrix_eval_ALL (gsl_matrix *m)
511 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
518 matrix_eval_ANY (gsl_matrix *m)
520 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
527 matrix_eval_ARSIN (gsl_matrix *m)
529 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
535 matrix_eval_ARTAN (gsl_matrix *m)
537 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
543 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
547 for (size_t i = 0; i < n; i++)
552 gsl_matrix *block = gsl_matrix_calloc (r, c);
554 for (size_t i = 0; i < n; i++)
556 for (size_t y = 0; y < m[i]->size1; y++)
557 for (size_t x = 0; x < m[i]->size2; x++)
558 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
566 matrix_eval_CHOL (gsl_matrix *m)
568 if (!gsl_linalg_cholesky_decomp1 (m))
570 for (size_t y = 0; y < m->size1; y++)
571 for (size_t x = y + 1; x < m->size2; x++)
572 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
574 for (size_t y = 0; y < m->size1; y++)
575 for (size_t x = 0; x < y; x++)
576 gsl_matrix_set (m, y, x, 0);
581 msg (SE, _("Input to CHOL function is not positive-definite."));
587 matrix_eval_col_extremum (gsl_matrix *m, bool min)
592 return gsl_matrix_alloc (1, 0);
594 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
595 for (size_t x = 0; x < m->size2; x++)
597 double ext = gsl_matrix_get (m, 0, x);
598 for (size_t y = 1; y < m->size1; y++)
600 double value = gsl_matrix_get (m, y, x);
601 if (min ? value < ext : value > ext)
604 gsl_matrix_set (cext, 0, x, ext);
610 matrix_eval_CMAX (gsl_matrix *m)
612 return matrix_eval_col_extremum (m, false);
616 matrix_eval_CMIN (gsl_matrix *m)
618 return matrix_eval_col_extremum (m, true);
622 matrix_eval_COS (gsl_matrix *m)
624 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
630 matrix_eval_col_sum (gsl_matrix *m, bool square)
635 return gsl_matrix_alloc (1, 0);
637 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
638 for (size_t x = 0; x < m->size2; x++)
641 for (size_t y = 0; y < m->size1; y++)
643 double d = gsl_matrix_get (m, y, x);
644 sum += square ? pow2 (d) : d;
646 gsl_matrix_set (result, 0, x, sum);
652 matrix_eval_CSSQ (gsl_matrix *m)
654 return matrix_eval_col_sum (m, true);
658 matrix_eval_CSUM (gsl_matrix *m)
660 return matrix_eval_col_sum (m, false);
664 compare_double_3way (const void *a_, const void *b_)
666 const double *a = a_;
667 const double *b = b_;
668 return *a < *b ? -1 : *a > *b;
672 matrix_eval_DESIGN (gsl_matrix *m)
674 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
675 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
676 gsl_matrix_transpose_memcpy (&m2, m);
678 for (size_t y = 0; y < m2.size1; y++)
679 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
681 size_t *n = xcalloc (m2.size1, sizeof *n);
683 for (size_t i = 0; i < m2.size1; i++)
685 double *row = tmp + m2.size2 * i;
686 for (size_t j = 0; j < m2.size2; )
689 for (k = j + 1; k < m2.size2; k++)
690 if (row[j] != row[k])
692 row[n[i]++] = row[j];
697 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
702 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
704 for (size_t i = 0; i < m->size2; i++)
709 const double *unique = tmp + m2.size2 * i;
710 for (size_t j = 0; j < n[i]; j++, x++)
712 double value = unique[j];
713 for (size_t y = 0; y < m->size1; y++)
714 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
725 matrix_eval_DET (gsl_matrix *m)
727 gsl_permutation *p = gsl_permutation_alloc (m->size1);
729 gsl_linalg_LU_decomp (m, p, &signum);
730 gsl_permutation_free (p);
731 return gsl_linalg_LU_det (m, signum);
735 matrix_eval_DIAG (gsl_matrix *m)
737 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
738 for (size_t i = 0; i < diag->size1; i++)
739 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
744 is_symmetric (const gsl_matrix *m)
746 if (m->size1 != m->size2)
749 for (size_t y = 0; y < m->size1; y++)
750 for (size_t x = 0; x < y; x++)
751 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
758 compare_double_desc (const void *a_, const void *b_)
760 const double *a = a_;
761 const double *b = b_;
762 return *a > *b ? -1 : *a < *b;
766 matrix_eval_EVAL (gsl_matrix *m)
768 if (!is_symmetric (m))
770 msg (SE, _("Argument of EVAL must be symmetric."));
774 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
775 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
776 gsl_vector v_eval = to_vector (eval);
777 gsl_eigen_symm (m, &v_eval, w);
778 gsl_eigen_symm_free (w);
780 assert (v_eval.stride == 1);
781 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
787 matrix_eval_EXP (gsl_matrix *m)
789 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
794 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
797 Charl Linssen <charl@itfromb.it>
801 matrix_eval_GINV (gsl_matrix *A)
806 gsl_matrix *tmp_mat = NULL;
809 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
810 tmp_mat = gsl_matrix_alloc (m, n);
811 gsl_matrix_transpose_memcpy (tmp_mat, A);
819 gsl_matrix *V = gsl_matrix_alloc (m, m);
820 gsl_vector *u = gsl_vector_alloc (m);
822 gsl_vector *tmp_vec = gsl_vector_alloc (m);
823 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
824 gsl_vector_free (tmp_vec);
827 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
828 gsl_matrix_set_zero (Sigma_pinv);
829 double cutoff = 1e-15 * gsl_vector_max (u);
831 for (size_t i = 0; i < m; ++i)
833 double x = gsl_vector_get (u, i);
834 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
837 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
838 gsl_matrix *U = gsl_matrix_calloc (n, n);
839 for (size_t i = 0; i < n; i++)
840 for (size_t j = 0; j < m; j++)
841 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
843 /* two dot products to obtain pseudoinverse */
844 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
845 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
850 A_pinv = gsl_matrix_alloc (n, m);
851 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
855 A_pinv = gsl_matrix_alloc (m, n);
856 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
859 gsl_matrix_free (tmp_mat);
860 gsl_matrix_free (tmp_mat2);
862 gsl_matrix_free (Sigma_pinv);
876 grade_compare_3way (const void *a_, const void *b_)
878 const struct grade *a = a_;
879 const struct grade *b = b_;
881 return (a->value < b->value ? -1
882 : a->value > b->value ? 1
890 matrix_eval_GRADE (gsl_matrix *m)
892 size_t n = m->size1 * m->size2;
893 struct grade *grades = xmalloc (n * sizeof *grades);
896 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
897 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
898 qsort (grades, n, sizeof *grades, grade_compare_3way);
900 for (size_t i = 0; i < n; i++)
901 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
909 dot (gsl_vector *a, gsl_vector *b)
912 for (size_t i = 0; i < a->size; i++)
913 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
918 norm2 (gsl_vector *v)
921 for (size_t i = 0; i < v->size; i++)
922 result += pow2 (gsl_vector_get (v, i));
929 return sqrt (norm2 (v));
933 matrix_eval_GSCH (gsl_matrix *v)
935 if (v->size2 < v->size1)
937 msg (SE, _("GSCH requires its argument to have at least as many columns "
938 "as rows, but it has dimensions %zu×%zu."),
942 if (!v->size1 || !v->size2)
945 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
947 for (size_t vx = 0; vx < v->size2; vx++)
949 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
950 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
952 gsl_vector_memcpy (&u_i, &v_i);
953 for (size_t j = 0; j < ux; j++)
955 gsl_vector u_j = gsl_matrix_column (u, j).vector;
956 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
957 for (size_t k = 0; k < u_i.size; k++)
958 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
959 - scale * gsl_vector_get (&u_j, k)));
962 double len = norm (&u_i);
965 gsl_vector_scale (&u_i, 1.0 / len);
966 if (++ux >= v->size1)
973 msg (SE, _("%zu×%zu argument to GSCH contains only "
974 "%zu linearly independent columns."),
975 v->size1, v->size2, ux);
985 matrix_eval_IDENT (double s1, double s2)
987 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
989 msg (SE, _("Arguments to IDENT must be non-negative integers."));
993 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
994 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1000 invert_matrix (gsl_matrix *x)
1002 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1004 gsl_linalg_LU_decomp (x, p, &signum);
1005 gsl_linalg_LU_invx (x, p);
1006 gsl_permutation_free (p);
1010 matrix_eval_INV (gsl_matrix *m)
1017 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1019 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1020 a->size2 * b->size2);
1022 for (size_t ar = 0; ar < a->size1; ar++)
1023 for (size_t br = 0; br < b->size1; br++, y++)
1026 for (size_t ac = 0; ac < a->size2; ac++)
1027 for (size_t bc = 0; bc < b->size2; bc++, x++)
1029 double av = gsl_matrix_get (a, ar, ac);
1030 double bv = gsl_matrix_get (b, br, bc);
1031 gsl_matrix_set (k, y, x, av * bv);
1038 matrix_eval_LG10 (gsl_matrix *m)
1040 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1046 matrix_eval_LN (gsl_matrix *m)
1048 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1054 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1056 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1059 for (size_t i = 1; i <= n * n; i++)
1061 gsl_matrix_set (m, y, x, i);
1063 size_t y1 = !y ? n - 1 : y - 1;
1064 size_t x1 = x + 1 >= n ? 0 : x + 1;
1065 if (gsl_matrix_get (m, y1, x1) == 0)
1071 y = y + 1 >= n ? 0 : y + 1;
1076 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1078 double a = gsl_matrix_get (m, y1, x1);
1079 double b = gsl_matrix_get (m, y2, x2);
1080 gsl_matrix_set (m, y1, x1, b);
1081 gsl_matrix_set (m, y2, x2, a);
1085 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1089 /* A. Umar, "On the Construction of Even Order Magic Squares",
1090 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1092 for (size_t i = 1; i <= n * n / 2; i++)
1094 gsl_matrix_set (m, y, x, i);
1104 for (size_t i = n * n; i > n * n / 2; i--)
1106 gsl_matrix_set (m, y, x, i);
1114 for (size_t y = 0; y < n; y++)
1115 for (size_t x = 0; x < n / 2; x++)
1117 unsigned int d = gsl_matrix_get (m, y, x);
1118 if (d % 2 != (y < n / 2))
1119 magic_exchange (m, y, x, y, n - x - 1);
1124 size_t x1 = n / 2 - 1;
1126 magic_exchange (m, y1, x1, y2, x1);
1127 magic_exchange (m, y1, x2, y2, x2);
1131 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1133 /* A. Umar, "On the Construction of Even Order Magic Squares",
1134 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1138 for (size_t i = 1; ; i++)
1140 gsl_matrix_set (m, y, x, i);
1141 if (++y == n / 2 - 1)
1153 for (size_t i = n * n; ; i--)
1155 gsl_matrix_set (m, y, x, i);
1156 if (++y == n / 2 - 1)
1165 for (size_t y = 0; y < n; y++)
1166 if (y != n / 2 - 1 && y != n / 2)
1167 for (size_t x = 0; x < n / 2; x++)
1169 unsigned int d = gsl_matrix_get (m, y, x);
1170 if (d % 2 != (y < n / 2))
1171 magic_exchange (m, y, x, y, n - x - 1);
1174 size_t a0 = (n * n - 2 * n) / 2 + 1;
1175 for (size_t i = 0; i < n / 2; i++)
1178 gsl_matrix_set (m, n / 2, i, a);
1179 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1181 for (size_t i = 0; i < n / 2; i++)
1183 size_t a = a0 + i + n / 2;
1184 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1185 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1187 for (size_t x = 1; x < n / 2; x += 2)
1188 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1189 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1190 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1191 size_t x1 = n / 2 - 2;
1192 size_t x2 = n / 2 + 1;
1193 size_t y1 = n / 2 - 2;
1194 size_t y2 = n / 2 + 1;
1195 magic_exchange (m, y1, x1, y2, x1);
1196 magic_exchange (m, y1, x2, y2, x2);
1200 matrix_eval_MAGIC (double n_)
1202 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1204 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1209 gsl_matrix *m = gsl_matrix_calloc (n, n);
1211 matrix_eval_MAGIC_odd (m, n);
1213 matrix_eval_MAGIC_singly_even (m, n);
1215 matrix_eval_MAGIC_doubly_even (m, n);
1220 matrix_eval_MAKE (double r, double c, double value)
1222 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1224 msg (SE, _("Size arguments to MAKE must be integers."));
1228 gsl_matrix *m = gsl_matrix_alloc (r, c);
1229 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1235 matrix_eval_MDIAG (gsl_vector *v)
1237 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1238 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1239 gsl_vector_memcpy (&diagonal, v);
1244 matrix_eval_MMAX (gsl_matrix *m)
1246 return gsl_matrix_max (m);
1250 matrix_eval_MMIN (gsl_matrix *m)
1252 return gsl_matrix_min (m);
1256 matrix_eval_MOD (gsl_matrix *m, double divisor)
1260 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1264 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1265 *d = fmod (*d, divisor);
1270 matrix_eval_MSSQ (gsl_matrix *m)
1273 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1279 matrix_eval_MSUM (gsl_matrix *m)
1282 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1288 matrix_eval_NCOL (gsl_matrix *m)
1294 matrix_eval_NROW (gsl_matrix *m)
1300 matrix_eval_RANK (gsl_matrix *m)
1302 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1303 gsl_linalg_QR_decomp (m, tau);
1304 gsl_vector_free (tau);
1306 return gsl_linalg_QRPT_rank (m, -1);
1310 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1312 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1314 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1319 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1321 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1322 "product of matrix dimensions."));
1326 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1329 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1331 gsl_matrix_set (dst, y1, x1, *d);
1342 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1347 return gsl_matrix_alloc (0, 1);
1349 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1350 for (size_t y = 0; y < m->size1; y++)
1352 double ext = gsl_matrix_get (m, y, 0);
1353 for (size_t x = 1; x < m->size2; x++)
1355 double value = gsl_matrix_get (m, y, x);
1356 if (min ? value < ext : value > ext)
1359 gsl_matrix_set (rext, y, 0, ext);
1365 matrix_eval_RMAX (gsl_matrix *m)
1367 return matrix_eval_row_extremum (m, false);
1371 matrix_eval_RMIN (gsl_matrix *m)
1373 return matrix_eval_row_extremum (m, true);
1377 matrix_eval_RND (gsl_matrix *m)
1379 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1391 rank_compare_3way (const void *a_, const void *b_)
1393 const struct rank *a = a_;
1394 const struct rank *b = b_;
1396 return a->value < b->value ? -1 : a->value > b->value;
1400 matrix_eval_RNKORDER (gsl_matrix *m)
1402 size_t n = m->size1 * m->size2;
1403 struct rank *ranks = xmalloc (n * sizeof *ranks);
1405 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1406 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1407 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1409 for (size_t i = 0; i < n; )
1412 for (j = i + 1; j < n; j++)
1413 if (ranks[i].value != ranks[j].value)
1416 double rank = (i + j + 1.0) / 2.0;
1417 for (size_t k = i; k < j; k++)
1418 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1429 matrix_eval_row_sum (gsl_matrix *m, bool square)
1434 return gsl_matrix_alloc (0, 1);
1436 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1437 for (size_t y = 0; y < m->size1; y++)
1440 for (size_t x = 0; x < m->size2; x++)
1442 double d = gsl_matrix_get (m, y, x);
1443 sum += square ? pow2 (d) : d;
1445 gsl_matrix_set (result, y, 0, sum);
1451 matrix_eval_RSSQ (gsl_matrix *m)
1453 return matrix_eval_row_sum (m, true);
1457 matrix_eval_RSUM (gsl_matrix *m)
1459 return matrix_eval_row_sum (m, false);
1463 matrix_eval_SIN (gsl_matrix *m)
1465 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1471 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1473 if (m1->size1 != m2->size1)
1475 msg (SE, _("SOLVE requires its arguments to have the same number of "
1476 "rows, but the first argument has dimensions %zu×%zu and "
1477 "the second %zu×%zu."),
1478 m1->size1, m1->size2,
1479 m2->size1, m2->size2);
1483 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1484 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1486 gsl_linalg_LU_decomp (m1, p, &signum);
1487 for (size_t i = 0; i < m2->size2; i++)
1489 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1490 gsl_vector xi = gsl_matrix_column (x, i).vector;
1491 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1493 gsl_permutation_free (p);
1498 matrix_eval_SQRT (gsl_matrix *m)
1500 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1504 msg (SE, _("Argument to SQRT must be nonnegative."));
1513 matrix_eval_SSCP (gsl_matrix *m)
1515 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1516 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1521 matrix_eval_SVAL (gsl_matrix *m)
1523 gsl_matrix *tmp_mat = NULL;
1524 if (m->size2 > m->size1)
1526 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1527 gsl_matrix_transpose_memcpy (tmp_mat, m);
1532 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1533 gsl_vector *S = gsl_vector_alloc (m->size2);
1534 gsl_vector *work = gsl_vector_alloc (m->size2);
1535 gsl_linalg_SV_decomp (m, V, S, work);
1537 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1538 for (size_t i = 0; i < m->size2; i++)
1539 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1541 gsl_matrix_free (V);
1542 gsl_vector_free (S);
1543 gsl_vector_free (work);
1544 gsl_matrix_free (tmp_mat);
1550 matrix_eval_SWEEP (gsl_matrix *m, double d)
1552 if (d < 1 || d > SIZE_MAX)
1554 msg (SE, _("Scalar argument to SWEEP must be integer."));
1558 if (k >= MIN (m->size1, m->size2))
1560 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1561 "equal to the smaller of the matrix argument's rows and "
1566 double m_kk = gsl_matrix_get (m, k, k);
1567 if (fabs (m_kk) > 1e-19)
1569 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1570 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1572 double m_ij = gsl_matrix_get (m, i, j);
1573 double m_ik = gsl_matrix_get (m, i, k);
1574 double m_kj = gsl_matrix_get (m, k, j);
1575 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1584 for (size_t i = 0; i < m->size1; i++)
1586 gsl_matrix_set (m, i, k, 0);
1587 gsl_matrix_set (m, k, i, 0);
1594 matrix_eval_TRACE (gsl_matrix *m)
1597 size_t n = MIN (m->size1, m->size2);
1598 for (size_t i = 0; i < n; i++)
1599 sum += gsl_matrix_get (m, i, i);
1604 matrix_eval_T (gsl_matrix *m)
1606 return matrix_eval_TRANSPOS (m);
1610 matrix_eval_TRANSPOS (gsl_matrix *m)
1612 if (m->size1 == m->size2)
1614 gsl_matrix_transpose (m);
1619 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1620 gsl_matrix_transpose_memcpy (t, m);
1626 matrix_eval_TRUNC (gsl_matrix *m)
1628 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1634 matrix_eval_UNIFORM (double r_, double c_)
1636 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1638 msg (SE, _("Arguments to UNIFORM must be integers."));
1643 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1645 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1649 gsl_matrix *m = gsl_matrix_alloc (r, c);
1650 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1651 *d = gsl_ran_flat (get_rng (), 0, 1);
1655 struct matrix_function
1659 size_t min_args, max_args;
1662 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1664 static const struct matrix_function *
1665 matrix_parse_function_name (const char *token)
1667 static const struct matrix_function functions[] =
1669 #define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1673 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1675 for (size_t i = 0; i < N_FUNCTIONS; i++)
1677 if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1678 return &functions[i];
1683 static struct read_file *
1684 read_file_create (struct matrix_state *s, struct file_handle *fh)
1686 for (size_t i = 0; i < s->n_read_files; i++)
1688 struct read_file *rf = s->read_files[i];
1696 struct read_file *rf = xmalloc (sizeof *rf);
1697 *rf = (struct read_file) { .file = fh };
1699 s->read_files = xrealloc (s->read_files,
1700 (s->n_read_files + 1) * sizeof *s->read_files);
1701 s->read_files[s->n_read_files++] = rf;
1705 static struct dfm_reader *
1706 read_file_open (struct read_file *rf)
1709 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1714 read_file_destroy (struct read_file *rf)
1718 fh_unref (rf->file);
1719 dfm_close_reader (rf->reader);
1720 free (rf->encoding);
1725 static struct write_file *
1726 write_file_create (struct matrix_state *s, struct file_handle *fh)
1728 for (size_t i = 0; i < s->n_write_files; i++)
1730 struct write_file *wf = s->write_files[i];
1738 struct write_file *wf = xmalloc (sizeof *wf);
1739 *wf = (struct write_file) { .file = fh };
1741 s->write_files = xrealloc (s->write_files,
1742 (s->n_write_files + 1) * sizeof *s->write_files);
1743 s->write_files[s->n_write_files++] = wf;
1747 static struct dfm_writer *
1748 write_file_open (struct write_file *wf)
1751 wf->writer = dfm_open_writer (wf->file, wf->encoding);
1756 write_file_destroy (struct write_file *wf)
1762 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
1763 wf->held->s.ss.length);
1764 u8_line_destroy (wf->held);
1768 fh_unref (wf->file);
1769 dfm_close_writer (wf->writer);
1770 free (wf->encoding);
1776 matrix_parse_function (struct matrix_state *s, const char *token,
1777 struct matrix_expr **exprp)
1780 if (lex_next_token (s->lexer, 1) != T_LPAREN)
1783 if (lex_match_id (s->lexer, "EOF"))
1786 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1790 if (!lex_force_match (s->lexer, T_RPAREN))
1796 struct read_file *rf = read_file_create (s, fh);
1798 struct matrix_expr *e = xmalloc (sizeof *e);
1799 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1804 const struct matrix_function *f = matrix_parse_function_name (token);
1808 lex_get_n (s->lexer, 2);
1810 struct matrix_expr *e = xmalloc (sizeof *e);
1811 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1813 size_t allocated_subs = 0;
1816 struct matrix_expr *sub = matrix_parse_expr (s);
1820 if (e->n_subs >= allocated_subs)
1821 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1822 e->subs[e->n_subs++] = sub;
1824 while (lex_match (s->lexer, T_COMMA));
1825 if (!lex_force_match (s->lexer, T_RPAREN))
1832 matrix_expr_destroy (e);
1836 static struct matrix_expr *
1837 matrix_parse_primary (struct matrix_state *s)
1839 if (lex_is_number (s->lexer))
1841 double number = lex_number (s->lexer);
1843 return matrix_expr_create_number (number);
1845 else if (lex_is_string (s->lexer))
1847 char string[sizeof (double)];
1848 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1852 memcpy (&number, string, sizeof number);
1853 return matrix_expr_create_number (number);
1855 else if (lex_match (s->lexer, T_LPAREN))
1857 struct matrix_expr *e = matrix_parse_or_xor (s);
1858 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1860 matrix_expr_destroy (e);
1865 else if (lex_match (s->lexer, T_LCURLY))
1867 struct matrix_expr *e = matrix_parse_curly_semi (s);
1868 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1870 matrix_expr_destroy (e);
1875 else if (lex_token (s->lexer) == T_ID)
1877 struct matrix_expr *retval;
1878 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1881 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1884 lex_error (s->lexer, _("Unknown variable %s."),
1885 lex_tokcstr (s->lexer));
1890 struct matrix_expr *e = xmalloc (sizeof *e);
1891 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1894 else if (lex_token (s->lexer) == T_ALL)
1896 struct matrix_expr *retval;
1897 if (matrix_parse_function (s, "ALL", &retval))
1901 lex_error (s->lexer, NULL);
1905 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1908 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1910 if (lex_match (s->lexer, T_COLON))
1917 *indexp = matrix_parse_expr (s);
1918 return *indexp != NULL;
1922 static struct matrix_expr *
1923 matrix_parse_postfix (struct matrix_state *s)
1925 struct matrix_expr *lhs = matrix_parse_primary (s);
1926 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1929 struct matrix_expr *i0;
1930 if (!matrix_parse_index_expr (s, &i0))
1932 matrix_expr_destroy (lhs);
1935 if (lex_match (s->lexer, T_RPAREN))
1937 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1938 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1939 else if (lex_match (s->lexer, T_COMMA))
1941 struct matrix_expr *i1;
1942 if (!matrix_parse_index_expr (s, &i1)
1943 || !lex_force_match (s->lexer, T_RPAREN))
1945 matrix_expr_destroy (lhs);
1946 matrix_expr_destroy (i0);
1947 matrix_expr_destroy (i1);
1950 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1951 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1952 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1957 lex_error_expecting (s->lexer, "`)'", "`,'");
1962 static struct matrix_expr *
1963 matrix_parse_unary (struct matrix_state *s)
1965 if (lex_match (s->lexer, T_DASH))
1967 struct matrix_expr *lhs = matrix_parse_unary (s);
1968 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1970 else if (lex_match (s->lexer, T_PLUS))
1971 return matrix_parse_unary (s);
1973 return matrix_parse_postfix (s);
1976 static struct matrix_expr *
1977 matrix_parse_seq (struct matrix_state *s)
1979 struct matrix_expr *start = matrix_parse_unary (s);
1980 if (!start || !lex_match (s->lexer, T_COLON))
1983 struct matrix_expr *end = matrix_parse_unary (s);
1986 matrix_expr_destroy (start);
1990 if (lex_match (s->lexer, T_COLON))
1992 struct matrix_expr *increment = matrix_parse_unary (s);
1995 matrix_expr_destroy (start);
1996 matrix_expr_destroy (end);
1999 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2002 return matrix_expr_create_2 (MOP_SEQ, start, end);
2005 static struct matrix_expr *
2006 matrix_parse_exp (struct matrix_state *s)
2008 struct matrix_expr *lhs = matrix_parse_seq (s);
2015 if (lex_match (s->lexer, T_EXP))
2017 else if (lex_match_phrase (s->lexer, "&**"))
2022 struct matrix_expr *rhs = matrix_parse_seq (s);
2025 matrix_expr_destroy (lhs);
2028 lhs = matrix_expr_create_2 (op, lhs, rhs);
2032 static struct matrix_expr *
2033 matrix_parse_mul_div (struct matrix_state *s)
2035 struct matrix_expr *lhs = matrix_parse_exp (s);
2042 if (lex_match (s->lexer, T_ASTERISK))
2044 else if (lex_match (s->lexer, T_SLASH))
2046 else if (lex_match_phrase (s->lexer, "&*"))
2048 else if (lex_match_phrase (s->lexer, "&/"))
2053 struct matrix_expr *rhs = matrix_parse_exp (s);
2056 matrix_expr_destroy (lhs);
2059 lhs = matrix_expr_create_2 (op, lhs, rhs);
2063 static struct matrix_expr *
2064 matrix_parse_add_sub (struct matrix_state *s)
2066 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2073 if (lex_match (s->lexer, T_PLUS))
2075 else if (lex_match (s->lexer, T_DASH))
2077 else if (lex_token (s->lexer) == T_NEG_NUM)
2082 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2085 matrix_expr_destroy (lhs);
2088 lhs = matrix_expr_create_2 (op, lhs, rhs);
2092 static struct matrix_expr *
2093 matrix_parse_relational (struct matrix_state *s)
2095 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2102 if (lex_match (s->lexer, T_GT))
2104 else if (lex_match (s->lexer, T_GE))
2106 else if (lex_match (s->lexer, T_LT))
2108 else if (lex_match (s->lexer, T_LE))
2110 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2112 else if (lex_match (s->lexer, T_NE))
2117 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2120 matrix_expr_destroy (lhs);
2123 lhs = matrix_expr_create_2 (op, lhs, rhs);
2127 static struct matrix_expr *
2128 matrix_parse_not (struct matrix_state *s)
2130 if (lex_match (s->lexer, T_NOT))
2132 struct matrix_expr *lhs = matrix_parse_not (s);
2133 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2136 return matrix_parse_relational (s);
2139 static struct matrix_expr *
2140 matrix_parse_and (struct matrix_state *s)
2142 struct matrix_expr *lhs = matrix_parse_not (s);
2146 while (lex_match (s->lexer, T_AND))
2148 struct matrix_expr *rhs = matrix_parse_not (s);
2151 matrix_expr_destroy (lhs);
2154 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2159 static struct matrix_expr *
2160 matrix_parse_or_xor (struct matrix_state *s)
2162 struct matrix_expr *lhs = matrix_parse_and (s);
2169 if (lex_match (s->lexer, T_OR))
2171 else if (lex_match_id (s->lexer, "XOR"))
2176 struct matrix_expr *rhs = matrix_parse_and (s);
2179 matrix_expr_destroy (lhs);
2182 lhs = matrix_expr_create_2 (op, lhs, rhs);
2186 static struct matrix_expr *
2187 matrix_parse_expr (struct matrix_state *s)
2189 return matrix_parse_or_xor (s);
2192 /* Expression evaluation. */
2195 matrix_op_eval (enum matrix_op op, double a, double b)
2199 case MOP_ADD_ELEMS: return a + b;
2200 case MOP_SUB_ELEMS: return a - b;
2201 case MOP_MUL_ELEMS: return a * b;
2202 case MOP_DIV_ELEMS: return a / b;
2203 case MOP_EXP_ELEMS: return pow (a, b);
2204 case MOP_GT: return a > b;
2205 case MOP_GE: return a >= b;
2206 case MOP_LT: return a < b;
2207 case MOP_LE: return a <= b;
2208 case MOP_EQ: return a == b;
2209 case MOP_NE: return a != b;
2210 case MOP_AND: return (a > 0) && (b > 0);
2211 case MOP_OR: return (a > 0) || (b > 0);
2212 case MOP_XOR: return (a > 0) != (b > 0);
2214 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2223 case MOP_PASTE_HORZ:
2224 case MOP_PASTE_VERT:
2240 matrix_op_name (enum matrix_op op)
2244 case MOP_ADD_ELEMS: return "+";
2245 case MOP_SUB_ELEMS: return "-";
2246 case MOP_MUL_ELEMS: return "&*";
2247 case MOP_DIV_ELEMS: return "&/";
2248 case MOP_EXP_ELEMS: return "&**";
2249 case MOP_GT: return ">";
2250 case MOP_GE: return ">=";
2251 case MOP_LT: return "<";
2252 case MOP_LE: return "<=";
2253 case MOP_EQ: return "=";
2254 case MOP_NE: return "<>";
2255 case MOP_AND: return "AND";
2256 case MOP_OR: return "OR";
2257 case MOP_XOR: return "XOR";
2259 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2268 case MOP_PASTE_HORZ:
2269 case MOP_PASTE_VERT:
2285 is_scalar (const gsl_matrix *m)
2287 return m->size1 == 1 && m->size2 == 1;
2291 to_scalar (const gsl_matrix *m)
2293 assert (is_scalar (m));
2294 return gsl_matrix_get (m, 0, 0);
2298 matrix_expr_evaluate_elementwise (enum matrix_op op,
2299 gsl_matrix *a, gsl_matrix *b)
2303 double be = to_scalar (b);
2304 for (size_t r = 0; r < a->size1; r++)
2305 for (size_t c = 0; c < a->size2; c++)
2307 double *ae = gsl_matrix_ptr (a, r, c);
2308 *ae = matrix_op_eval (op, *ae, be);
2312 else if (is_scalar (a))
2314 double ae = to_scalar (a);
2315 for (size_t r = 0; r < b->size1; r++)
2316 for (size_t c = 0; c < b->size2; c++)
2318 double *be = gsl_matrix_ptr (b, r, c);
2319 *be = matrix_op_eval (op, ae, *be);
2323 else if (a->size1 == b->size1 && a->size2 == b->size2)
2325 for (size_t r = 0; r < a->size1; r++)
2326 for (size_t c = 0; c < a->size2; c++)
2328 double *ae = gsl_matrix_ptr (a, r, c);
2329 double be = gsl_matrix_get (b, r, c);
2330 *ae = matrix_op_eval (op, *ae, be);
2336 msg (SE, _("Operands to %s must have the same dimensions or one "
2337 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
2338 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2344 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2346 if (is_scalar (a) || is_scalar (b))
2347 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2349 if (a->size2 != b->size1)
2351 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
2352 "not conformable for multiplication."),
2353 a->size1, a->size2, b->size1, b->size2);
2357 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2358 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2363 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2365 gsl_matrix *tmp = *a;
2371 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2374 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2375 swap_matrix (z, tmp);
2379 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2381 mul_matrix (x, *x, *x, tmp);
2385 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2388 if (x->size1 != x->size2)
2390 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2391 "the left-hand size, not one with dimensions %zu×%zu."),
2392 x->size1, x->size2);
2397 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2398 "right-hand side, not a matrix with dimensions %zu×%zu."),
2399 b->size1, b->size2);
2402 double bf = to_scalar (b);
2403 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2405 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2406 "or outside the valid range."), bf);
2411 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2413 gsl_matrix_set_identity (y);
2417 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2419 for (unsigned long int n = labs (bl); n > 1; n /= 2)
2422 mul_matrix (&y, x, y, &t);
2423 square_matrix (&x, &t);
2426 square_matrix (&x, &t);
2428 mul_matrix (&y, x, y, &t);
2432 /* Garbage collection.
2434 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2435 a permutation of them. We are returning one of them; that one must not be
2436 destroyed. We must not destroy 'x_' because the caller owns it. */
2438 gsl_matrix_free (y_);
2440 gsl_matrix_free (t_);
2446 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2449 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2451 msg (SE, _("All operands of : operator must be scalars."));
2455 long int start = to_scalar (start_);
2456 long int end = to_scalar (end_);
2457 long int by = by_ ? to_scalar (by_) : 1;
2461 msg (SE, _("The increment operand to : must be nonzero."));
2465 long int n = (end >= start && by > 0 ? (end - start + by) / by
2466 : end <= start && by < 0 ? (start - end - by) / -by
2468 gsl_matrix *m = gsl_matrix_alloc (1, n);
2469 for (long int i = 0; i < n; i++)
2470 gsl_matrix_set (m, 0, i, start + i * by);
2475 matrix_expr_evaluate_not (gsl_matrix *a)
2477 for (size_t r = 0; r < a->size1; r++)
2478 for (size_t c = 0; c < a->size2; c++)
2480 double *ae = gsl_matrix_ptr (a, r, c);
2487 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2489 if (a->size1 != b->size1)
2491 if (!a->size1 || !a->size2)
2493 else if (!b->size1 || !b->size2)
2496 msg (SE, _("All columns in a matrix must have the same number of rows, "
2497 "but this tries to paste matrices with %zu and %zu rows."),
2498 a->size1, b->size1);
2502 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2503 for (size_t y = 0; y < a->size1; y++)
2505 for (size_t x = 0; x < a->size2; x++)
2506 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2507 for (size_t x = 0; x < b->size2; x++)
2508 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2514 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2516 if (a->size2 != b->size2)
2518 if (!a->size1 || !a->size2)
2520 else if (!b->size1 || !b->size2)
2523 msg (SE, _("All rows in a matrix must have the same number of columns, "
2524 "but this tries to stack matrices with %zu and %zu columns."),
2525 a->size2, b->size2);
2529 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2530 for (size_t x = 0; x < a->size2; x++)
2532 for (size_t y = 0; y < a->size1; y++)
2533 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2534 for (size_t y = 0; y < b->size1; y++)
2535 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2541 matrix_to_vector (gsl_matrix *m)
2544 gsl_vector v = to_vector (m);
2545 assert (v.block == m->block || !v.block);
2549 return xmemdup (&v, sizeof v);
2563 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2566 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2567 enum index_type index_type, size_t other_size,
2568 struct index_vector *iv)
2577 msg (SE, _("Vector index must be scalar or vector, not a "
2579 m->size1, m->size2);
2583 msg (SE, _("Matrix row index must be scalar or vector, not a "
2585 m->size1, m->size2);
2589 msg (SE, _("Matrix column index must be scalar or vector, not a "
2591 m->size1, m->size2);
2597 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2598 *iv = (struct index_vector) {
2599 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2602 for (size_t i = 0; i < v.size; i++)
2604 double index = gsl_vector_get (&v, i);
2605 if (index < 1 || index >= size + 1)
2610 msg (SE, _("Index %g is out of range for vector "
2611 "with %zu elements."), index, size);
2615 msg (SE, _("%g is not a valid row index for "
2616 "a %zu×%zu matrix."),
2617 index, size, other_size);
2621 msg (SE, _("%g is not a valid column index for "
2622 "a %zu×%zu matrix."),
2623 index, other_size, size);
2630 iv->indexes[i] = index - 1;
2636 *iv = (struct index_vector) {
2637 .indexes = xnmalloc (size, sizeof *iv->indexes),
2640 for (size_t i = 0; i < size; i++)
2647 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2649 if (!is_vector (sm))
2651 msg (SE, _("Vector index operator may not be applied to "
2652 "a %zu×%zu matrix."),
2653 sm->size1, sm->size2);
2661 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2663 if (!matrix_expr_evaluate_vec_all (sm))
2666 gsl_vector sv = to_vector (sm);
2667 struct index_vector iv;
2668 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2671 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2672 sm->size1 == 1 ? iv.n : 1);
2673 gsl_vector dv = to_vector (dm);
2674 for (size_t dx = 0; dx < iv.n; dx++)
2676 size_t sx = iv.indexes[dx];
2677 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2685 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2688 struct index_vector iv0;
2689 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2692 struct index_vector iv1;
2693 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2700 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2701 for (size_t dy = 0; dy < iv0.n; dy++)
2703 size_t sy = iv0.indexes[dy];
2705 for (size_t dx = 0; dx < iv1.n; dx++)
2707 size_t sx = iv1.indexes[dx];
2708 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2716 #define F(NAME, PROTOTYPE) \
2717 static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2718 const char *name, gsl_matrix *[], size_t, \
2719 matrix_proto_##PROTOTYPE *);
2724 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2726 if (!is_scalar (subs[index]))
2728 msg (SE, _("Function %s argument %zu must be a scalar, "
2729 "not a %zu×%zu matrix."),
2730 name, index + 1, subs[index]->size1, subs[index]->size2);
2737 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2739 if (!is_vector (subs[index]))
2741 msg (SE, _("Function %s argument %zu must be a vector, "
2742 "not a %zu×%zu matrix."),
2743 name, index + 1, subs[index]->size1, subs[index]->size2);
2750 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2752 for (size_t i = 0; i < n_subs; i++)
2754 if (!check_scalar_arg (name, subs, i))
2756 d[i] = to_scalar (subs[i]);
2762 matrix_expr_evaluate_m_d (const char *name,
2763 gsl_matrix *subs[], size_t n_subs,
2764 matrix_proto_m_d *f)
2766 assert (n_subs == 1);
2769 return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2773 matrix_expr_evaluate_m_dd (const char *name,
2774 gsl_matrix *subs[], size_t n_subs,
2775 matrix_proto_m_dd *f)
2777 assert (n_subs == 2);
2780 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2784 matrix_expr_evaluate_m_ddd (const char *name,
2785 gsl_matrix *subs[], size_t n_subs,
2786 matrix_proto_m_ddd *f)
2788 assert (n_subs == 3);
2791 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2795 matrix_expr_evaluate_m_m (const char *name UNUSED,
2796 gsl_matrix *subs[], size_t n_subs,
2797 matrix_proto_m_m *f)
2799 assert (n_subs == 1);
2804 matrix_expr_evaluate_m_md (const char *name UNUSED,
2805 gsl_matrix *subs[], size_t n_subs,
2806 matrix_proto_m_md *f)
2808 assert (n_subs == 2);
2809 if (!check_scalar_arg (name, subs, 1))
2811 return f (subs[0], to_scalar (subs[1]));
2815 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2816 gsl_matrix *subs[], size_t n_subs,
2817 matrix_proto_m_mdd *f)
2819 assert (n_subs == 3);
2820 if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2822 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2826 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2827 gsl_matrix *subs[], size_t n_subs,
2828 matrix_proto_m_mm *f)
2830 assert (n_subs == 2);
2831 return f (subs[0], subs[1]);
2835 matrix_expr_evaluate_m_v (const char *name,
2836 gsl_matrix *subs[], size_t n_subs,
2837 matrix_proto_m_v *f)
2839 assert (n_subs == 1);
2840 if (!check_vector_arg (name, subs, 0))
2842 gsl_vector v = to_vector (subs[0]);
2847 matrix_expr_evaluate_d_m (const char *name UNUSED,
2848 gsl_matrix *subs[], size_t n_subs,
2849 matrix_proto_d_m *f)
2851 assert (n_subs == 1);
2853 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2854 gsl_matrix_set (m, 0, 0, f (subs[0]));
2859 matrix_expr_evaluate_m_any (const char *name UNUSED,
2860 gsl_matrix *subs[], size_t n_subs,
2861 matrix_proto_m_any *f)
2863 return f (subs, n_subs);
2867 matrix_expr_evaluate_IDENT (const char *name,
2868 gsl_matrix *subs[], size_t n_subs,
2869 matrix_proto_IDENT *f)
2871 assert (n_subs <= 2);
2874 if (!to_scalar_args (name, subs, n_subs, d))
2877 return f (d[0], d[n_subs - 1]);
2881 matrix_expr_evaluate (const struct matrix_expr *e)
2883 if (e->op == MOP_NUMBER)
2885 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2886 gsl_matrix_set (m, 0, 0, e->number);
2889 else if (e->op == MOP_VARIABLE)
2891 const gsl_matrix *src = e->variable->value;
2894 msg (SE, _("Uninitialized variable %s used in expression."),
2899 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2900 gsl_matrix_memcpy (dst, src);
2903 else if (e->op == MOP_EOF)
2905 struct dfm_reader *reader = read_file_open (e->eof);
2906 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2907 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2911 enum { N_LOCAL = 3 };
2912 gsl_matrix *local_subs[N_LOCAL];
2913 gsl_matrix **subs = (e->n_subs < N_LOCAL
2915 : xmalloc (e->n_subs * sizeof *subs));
2917 for (size_t i = 0; i < e->n_subs; i++)
2919 subs[i] = matrix_expr_evaluate (e->subs[i]);
2922 for (size_t j = 0; j < i; j++)
2923 gsl_matrix_free (subs[j]);
2924 if (subs != local_subs)
2930 gsl_matrix *result = NULL;
2933 #define F(NAME, PROTOTYPE) \
2934 case MOP_F_##NAME: \
2935 result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \
2937 matrix_eval_##NAME); \
2943 gsl_matrix_scale (subs[0], -1.0);
2961 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2965 result = matrix_expr_evaluate_not (subs[0]);
2969 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2973 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2977 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2981 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2984 case MOP_PASTE_HORZ:
2985 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2988 case MOP_PASTE_VERT:
2989 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2993 result = gsl_matrix_alloc (0, 0);
2997 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
3001 result = matrix_expr_evaluate_vec_all (subs[0]);
3005 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
3009 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
3013 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
3022 for (size_t i = 0; i < e->n_subs; i++)
3023 if (subs[i] != result)
3024 gsl_matrix_free (subs[i]);
3025 if (subs != local_subs)
3031 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
3034 gsl_matrix *m = matrix_expr_evaluate (e);
3040 msg (SE, _("Expression for %s must evaluate to scalar, "
3041 "not a %zu×%zu matrix."),
3042 context, m->size1, m->size2);
3043 gsl_matrix_free (m);
3048 gsl_matrix_free (m);
3053 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
3057 if (!matrix_expr_evaluate_scalar (e, context, &d))
3061 if (d < LONG_MIN || d > LONG_MAX)
3063 msg (SE, _("Expression for %s is outside the integer range."), context);
3070 struct matrix_lvalue
3072 struct matrix_var *var;
3073 struct matrix_expr *indexes[2];
3078 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
3082 for (size_t i = 0; i < lvalue->n_indexes; i++)
3083 matrix_expr_destroy (lvalue->indexes[i]);
3088 static struct matrix_lvalue *
3089 matrix_lvalue_parse (struct matrix_state *s)
3091 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3092 if (!lex_force_id (s->lexer))
3094 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3095 if (lex_next_token (s->lexer, 1) == T_LPAREN)
3099 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3104 lex_get_n (s->lexer, 2);
3106 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3108 lvalue->n_indexes++;
3110 if (lex_match (s->lexer, T_COMMA))
3112 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3114 lvalue->n_indexes++;
3116 if (!lex_force_match (s->lexer, T_RPAREN))
3122 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3128 matrix_lvalue_destroy (lvalue);
3133 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3134 enum index_type index_type, size_t other_size,
3135 struct index_vector *iv)
3140 m = matrix_expr_evaluate (e);
3147 bool ok = matrix_normalize_index_vector (m, size, index_type,
3149 gsl_matrix_free (m);
3154 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3155 struct index_vector *iv, gsl_matrix *sm)
3157 gsl_vector dv = to_vector (lvalue->var->value);
3159 /* Convert source matrix 'sm' to source vector 'sv'. */
3160 if (!is_vector (sm))
3162 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
3163 sm->size1, sm->size2);
3166 gsl_vector sv = to_vector (sm);
3167 if (iv->n != sv.size)
3169 msg (SE, _("Can't assign %zu-element vector "
3170 "to %zu-element subvector."), sv.size, iv->n);
3174 /* Assign elements. */
3175 for (size_t x = 0; x < iv->n; x++)
3176 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3181 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3182 struct index_vector *iv0,
3183 struct index_vector *iv1, gsl_matrix *sm)
3185 gsl_matrix *dm = lvalue->var->value;
3187 /* Convert source matrix 'sm' to source vector 'sv'. */
3188 if (iv0->n != sm->size1)
3190 msg (SE, _("Row index vector for assignment to %s has %zu elements "
3191 "but source matrix has %zu rows."),
3192 lvalue->var->name, iv0->n, sm->size1);
3195 if (iv1->n != sm->size2)
3197 msg (SE, _("Column index vector for assignment to %s has %zu elements "
3198 "but source matrix has %zu columns."),
3199 lvalue->var->name, iv1->n, sm->size2);
3203 /* Assign elements. */
3204 for (size_t y = 0; y < iv0->n; y++)
3206 size_t f0 = iv0->indexes[y];
3207 for (size_t x = 0; x < iv1->n; x++)
3209 size_t f1 = iv1->indexes[x];
3210 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3216 /* Takes ownership of M. */
3218 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3219 struct index_vector *iv0, struct index_vector *iv1,
3222 if (!lvalue->n_indexes)
3224 gsl_matrix_free (lvalue->var->value);
3225 lvalue->var->value = sm;
3230 bool ok = (lvalue->n_indexes == 1
3231 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3232 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3233 free (iv0->indexes);
3234 free (iv1->indexes);
3235 gsl_matrix_free (sm);
3241 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3242 struct index_vector *iv0,
3243 struct index_vector *iv1)
3245 *iv0 = INDEX_VECTOR_INIT;
3246 *iv1 = INDEX_VECTOR_INIT;
3247 if (!lvalue->n_indexes)
3250 /* Validate destination matrix exists and has the right shape. */
3251 gsl_matrix *dm = lvalue->var->value;
3254 msg (SE, _("Undefined variable %s."), lvalue->var->name);
3257 else if (dm->size1 == 0 || dm->size2 == 0)
3259 msg (SE, _("Cannot index %zu×%zu matrix."),
3260 dm->size1, dm->size2);
3263 else if (lvalue->n_indexes == 1)
3265 if (!is_vector (dm))
3267 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
3268 dm->size1, dm->size2, lvalue->var->name);
3271 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3272 MAX (dm->size1, dm->size2),
3277 assert (lvalue->n_indexes == 2);
3278 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3279 IV_ROW, dm->size2, iv0))
3282 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3283 IV_COLUMN, dm->size1, iv1))
3285 free (iv0->indexes);
3292 /* Takes ownership of M. */
3294 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3296 struct index_vector iv0, iv1;
3297 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3299 gsl_matrix_free (sm);
3303 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3307 /* Matrix command. */
3311 struct matrix_cmd **commands;
3315 static void matrix_cmds_uninit (struct matrix_cmds *);
3319 enum matrix_cmd_type
3342 struct compute_command
3344 struct matrix_lvalue *lvalue;
3345 struct matrix_expr *rvalue;
3349 struct print_command
3351 struct matrix_expr *expression;
3352 bool use_default_format;
3353 struct fmt_spec format;
3355 int space; /* -1 means NEWPAGE. */
3359 struct string_array literals; /* CLABELS/RLABELS. */
3360 struct matrix_expr *expr; /* CNAMES/RNAMES. */
3366 struct do_if_command
3370 struct matrix_expr *condition;
3371 struct matrix_cmds commands;
3381 struct matrix_var *var;
3382 struct matrix_expr *start;
3383 struct matrix_expr *end;
3384 struct matrix_expr *increment;
3386 /* Loop conditions. */
3387 struct matrix_expr *top_condition;
3388 struct matrix_expr *bottom_condition;
3391 struct matrix_cmds commands;
3395 struct display_command
3397 struct matrix_state *state;
3403 struct release_command
3405 struct matrix_var **vars;
3412 struct matrix_expr *expression;
3413 struct save_file *sf;
3419 struct read_file *rf;
3420 struct matrix_lvalue *dst;
3421 struct matrix_expr *size;
3423 enum fmt_type format;
3430 struct write_command
3432 struct write_file *wf;
3433 struct matrix_expr *expression;
3436 /* If this is nonnull, WRITE uses this format.
3438 If this is NULL, WRITE uses free-field format with as many
3439 digits of precision as needed. */
3440 struct fmt_spec *format;
3449 struct matrix_lvalue *dst;
3450 struct dataset *dataset;
3451 struct file_handle *file;
3453 struct var_syntax *vars;
3455 struct matrix_var *names;
3457 /* Treatment of missing values. */
3462 MGET_ERROR, /* Flag error on command. */
3463 MGET_ACCEPT, /* Accept without change, user-missing only. */
3464 MGET_OMIT, /* Drop this case. */
3465 MGET_RECODE /* Recode to 'substitute'. */
3468 double substitute; /* MGET_RECODE only. */
3474 struct msave_command
3476 struct msave_common *common;
3478 struct matrix_expr *expr;
3479 const char *rowtype;
3480 struct matrix_expr *factors;
3481 struct matrix_expr *splits;
3487 struct matrix_state *state;
3488 struct file_handle *file;
3490 struct stringi_set rowtypes;
3494 struct eigen_command
3496 struct matrix_expr *expr;
3497 struct matrix_var *evec;
3498 struct matrix_var *eval;
3502 struct setdiag_command
3504 struct matrix_var *dst;
3505 struct matrix_expr *expr;
3511 struct matrix_expr *expr;
3512 struct matrix_var *u;
3513 struct matrix_var *s;
3514 struct matrix_var *v;
3520 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3521 static bool matrix_cmd_execute (struct matrix_cmd *);
3522 static void matrix_cmd_destroy (struct matrix_cmd *);
3525 static struct matrix_cmd *
3526 matrix_parse_compute (struct matrix_state *s)
3528 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3529 *cmd = (struct matrix_cmd) {
3530 .type = MCMD_COMPUTE,
3531 .compute = { .lvalue = NULL },
3534 struct compute_command *compute = &cmd->compute;
3535 compute->lvalue = matrix_lvalue_parse (s);
3536 if (!compute->lvalue)
3539 if (!lex_force_match (s->lexer, T_EQUALS))
3542 compute->rvalue = matrix_parse_expr (s);
3543 if (!compute->rvalue)
3549 matrix_cmd_destroy (cmd);
3554 matrix_cmd_execute_compute (struct compute_command *compute)
3556 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3560 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3564 print_labels_destroy (struct print_labels *labels)
3568 string_array_destroy (&labels->literals);
3569 matrix_expr_destroy (labels->expr);
3575 parse_literal_print_labels (struct matrix_state *s,
3576 struct print_labels **labelsp)
3578 lex_match (s->lexer, T_EQUALS);
3579 print_labels_destroy (*labelsp);
3580 *labelsp = xzalloc (sizeof **labelsp);
3581 while (lex_token (s->lexer) != T_SLASH
3582 && lex_token (s->lexer) != T_ENDCMD
3583 && lex_token (s->lexer) != T_STOP)
3585 struct string label = DS_EMPTY_INITIALIZER;
3586 while (lex_token (s->lexer) != T_COMMA
3587 && lex_token (s->lexer) != T_SLASH
3588 && lex_token (s->lexer) != T_ENDCMD
3589 && lex_token (s->lexer) != T_STOP)
3591 if (!ds_is_empty (&label))
3592 ds_put_byte (&label, ' ');
3594 if (lex_token (s->lexer) == T_STRING)
3595 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3598 char *rep = lex_next_representation (s->lexer, 0, 0);
3599 ds_put_cstr (&label, rep);
3604 string_array_append_nocopy (&(*labelsp)->literals,
3605 ds_steal_cstr (&label));
3607 if (!lex_match (s->lexer, T_COMMA))
3613 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3615 lex_match (s->lexer, T_EQUALS);
3616 struct matrix_expr *e = matrix_parse_exp (s);
3620 print_labels_destroy (*labelsp);
3621 *labelsp = xzalloc (sizeof **labelsp);
3622 (*labelsp)->expr = e;
3626 static struct matrix_cmd *
3627 matrix_parse_print (struct matrix_state *s)
3629 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3630 *cmd = (struct matrix_cmd) {
3633 .use_default_format = true,
3637 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3640 for (size_t i = 0; ; i++)
3642 enum token_type t = lex_next_token (s->lexer, i);
3643 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3645 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3647 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3650 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3655 cmd->print.expression = matrix_parse_exp (s);
3656 if (!cmd->print.expression)
3660 while (lex_match (s->lexer, T_SLASH))
3662 if (lex_match_id (s->lexer, "FORMAT"))
3664 lex_match (s->lexer, T_EQUALS);
3665 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3667 cmd->print.use_default_format = false;
3669 else if (lex_match_id (s->lexer, "TITLE"))
3671 lex_match (s->lexer, T_EQUALS);
3672 if (!lex_force_string (s->lexer))
3674 free (cmd->print.title);
3675 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3678 else if (lex_match_id (s->lexer, "SPACE"))
3680 lex_match (s->lexer, T_EQUALS);
3681 if (lex_match_id (s->lexer, "NEWPAGE"))
3682 cmd->print.space = -1;
3683 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3685 cmd->print.space = lex_integer (s->lexer);
3691 else if (lex_match_id (s->lexer, "RLABELS"))
3692 parse_literal_print_labels (s, &cmd->print.rlabels);
3693 else if (lex_match_id (s->lexer, "CLABELS"))
3694 parse_literal_print_labels (s, &cmd->print.clabels);
3695 else if (lex_match_id (s->lexer, "RNAMES"))
3697 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3700 else if (lex_match_id (s->lexer, "CNAMES"))
3702 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3707 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3708 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3716 matrix_cmd_destroy (cmd);
3721 matrix_is_integer (const gsl_matrix *m)
3723 for (size_t y = 0; y < m->size1; y++)
3724 for (size_t x = 0; x < m->size2; x++)
3726 double d = gsl_matrix_get (m, y, x);
3734 matrix_max_magnitude (const gsl_matrix *m)
3737 for (size_t y = 0; y < m->size1; y++)
3738 for (size_t x = 0; x < m->size2; x++)
3740 double d = fabs (gsl_matrix_get (m, y, x));
3748 format_fits (struct fmt_spec format, double x)
3750 char *s = data_out (&(union value) { .f = x }, NULL,
3751 &format, settings_get_fmt_settings ());
3752 bool fits = *s != '*' && !strchr (s, 'E');
3757 static struct fmt_spec
3758 default_format (const gsl_matrix *m, int *log_scale)
3762 double max = matrix_max_magnitude (m);
3764 if (matrix_is_integer (m))
3765 for (int w = 1; w <= 12; w++)
3767 struct fmt_spec format = { .type = FMT_F, .w = w };
3768 if (format_fits (format, -max))
3772 if (max >= 1e9 || max <= 1e-4)
3775 snprintf (s, sizeof s, "%.1e", max);
3777 const char *e = strchr (s, 'e');
3779 *log_scale = atoi (e + 1);
3782 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3786 trimmed_string (double d)
3788 struct substring s = ss_buffer ((char *) &d, sizeof d);
3789 ss_rtrim (&s, ss_cstr (" "));
3790 return ss_xstrdup (s);
3793 static struct string_array *
3794 print_labels_get (const struct print_labels *labels, size_t n,
3795 const char *prefix, bool truncate)
3800 struct string_array *out = xzalloc (sizeof *out);
3801 if (labels->literals.n)
3802 string_array_clone (out, &labels->literals);
3803 else if (labels->expr)
3805 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3806 if (m && is_vector (m))
3808 gsl_vector v = to_vector (m);
3809 for (size_t i = 0; i < v.size; i++)
3810 string_array_append_nocopy (out, trimmed_string (
3811 gsl_vector_get (&v, i)));
3813 gsl_matrix_free (m);
3819 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3822 string_array_append (out, "");
3825 string_array_delete (out, out->n - 1);
3828 for (size_t i = 0; i < out->n; i++)
3830 char *s = out->strings[i];
3831 s[strnlen (s, 8)] = '\0';
3838 matrix_cmd_print_space (int space)
3841 output_item_submit (page_break_item_create ());
3842 for (int i = 0; i < space; i++)
3843 output_log ("%s", "");
3847 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3848 struct fmt_spec format, int log_scale)
3850 matrix_cmd_print_space (print->space);
3852 output_log ("%s", print->title);
3854 output_log (" 10 ** %d X", log_scale);
3856 struct string_array *clabels = print_labels_get (print->clabels,
3857 m->size2, "col", true);
3858 if (clabels && format.w < 8)
3860 struct string_array *rlabels = print_labels_get (print->rlabels,
3861 m->size1, "row", true);
3865 struct string line = DS_EMPTY_INITIALIZER;
3867 ds_put_byte_multiple (&line, ' ', 8);
3868 for (size_t x = 0; x < m->size2; x++)
3869 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3870 output_log_nocopy (ds_steal_cstr (&line));
3873 double scale = pow (10.0, log_scale);
3874 bool numeric = fmt_is_numeric (format.type);
3875 for (size_t y = 0; y < m->size1; y++)
3877 struct string line = DS_EMPTY_INITIALIZER;
3879 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3881 for (size_t x = 0; x < m->size2; x++)
3883 double f = gsl_matrix_get (m, y, x);
3885 ? data_out (&(union value) { .f = f / scale}, NULL,
3886 &format, settings_get_fmt_settings ())
3887 : trimmed_string (f));
3888 ds_put_format (&line, " %s", s);
3891 output_log_nocopy (ds_steal_cstr (&line));
3894 string_array_destroy (rlabels);
3896 string_array_destroy (clabels);
3901 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3902 const struct print_labels *print_labels, size_t n,
3903 const char *name, const char *prefix)
3905 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3907 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3908 for (size_t i = 0; i < n; i++)
3909 pivot_category_create_leaf (
3911 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3912 : pivot_value_new_integer (i + 1)));
3914 d->hide_all_labels = true;
3915 string_array_destroy (labels);
3920 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3921 struct fmt_spec format, int log_scale)
3923 struct pivot_table *table = pivot_table_create__ (
3924 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3927 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3929 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3930 N_("Columns"), "col");
3932 struct pivot_footnote *footnote = NULL;
3935 char *s = xasprintf ("× 10**%d", log_scale);
3936 footnote = pivot_table_create_footnote (
3937 table, pivot_value_new_user_text_nocopy (s));
3940 double scale = pow (10.0, log_scale);
3941 bool numeric = fmt_is_numeric (format.type);
3942 for (size_t y = 0; y < m->size1; y++)
3943 for (size_t x = 0; x < m->size2; x++)
3945 double f = gsl_matrix_get (m, y, x);
3946 struct pivot_value *v;
3949 v = pivot_value_new_number (f / scale);
3950 v->numeric.format = format;
3953 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3955 pivot_value_add_footnote (v, footnote);
3956 pivot_table_put2 (table, y, x, v);
3959 pivot_table_submit (table);
3963 matrix_cmd_execute_print (const struct print_command *print)
3965 if (print->expression)
3967 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3972 struct fmt_spec format = (print->use_default_format
3973 ? default_format (m, &log_scale)
3976 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3977 matrix_cmd_print_text (print, m, format, log_scale);
3979 matrix_cmd_print_tables (print, m, format, log_scale);
3981 gsl_matrix_free (m);
3985 matrix_cmd_print_space (print->space);
3987 output_log ("%s", print->title);
3994 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3995 const char *command_name,
3996 const char *stop1, const char *stop2)
3998 lex_end_of_command (s->lexer);
3999 lex_discard_rest_of_command (s->lexer);
4001 size_t allocated = 0;
4004 while (lex_token (s->lexer) == T_ENDCMD)
4007 if (lex_at_phrase (s->lexer, stop1)
4008 || (stop2 && lex_at_phrase (s->lexer, stop2)))
4011 if (lex_at_phrase (s->lexer, "END MATRIX"))
4013 msg (SE, _("Premature END MATRIX within %s."), command_name);
4017 struct matrix_cmd *cmd = matrix_parse_command (s);
4021 if (c->n >= allocated)
4022 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
4023 c->commands[c->n++] = cmd;
4028 matrix_parse_do_if_clause (struct matrix_state *s,
4029 struct do_if_command *ifc,
4030 bool parse_condition,
4031 size_t *allocated_clauses)
4033 if (ifc->n_clauses >= *allocated_clauses)
4034 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
4035 sizeof *ifc->clauses);
4036 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
4037 *c = (struct do_if_clause) { .condition = NULL };
4039 if (parse_condition)
4041 c->condition = matrix_parse_expr (s);
4046 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
4049 static struct matrix_cmd *
4050 matrix_parse_do_if (struct matrix_state *s)
4052 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4053 *cmd = (struct matrix_cmd) {
4055 .do_if = { .n_clauses = 0 }
4058 size_t allocated_clauses = 0;
4061 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
4064 while (lex_match_phrase (s->lexer, "ELSE IF"));
4066 if (lex_match_id (s->lexer, "ELSE")
4067 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
4070 if (!lex_match_phrase (s->lexer, "END IF"))
4075 matrix_cmd_destroy (cmd);
4080 matrix_cmd_execute_do_if (struct do_if_command *cmd)
4082 for (size_t i = 0; i < cmd->n_clauses; i++)
4084 struct do_if_clause *c = &cmd->clauses[i];
4088 if (!matrix_expr_evaluate_scalar (c->condition,
4089 i ? "ELSE IF" : "DO IF",
4094 for (size_t j = 0; j < c->commands.n; j++)
4095 if (!matrix_cmd_execute (c->commands.commands[j]))
4102 static struct matrix_cmd *
4103 matrix_parse_loop (struct matrix_state *s)
4105 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4106 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4108 struct loop_command *loop = &cmd->loop;
4109 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4111 struct substring name = lex_tokss (s->lexer);
4112 loop->var = matrix_var_lookup (s, name);
4114 loop->var = matrix_var_create (s, name);
4119 loop->start = matrix_parse_expr (s);
4120 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4123 loop->end = matrix_parse_expr (s);
4127 if (lex_match (s->lexer, T_BY))
4129 loop->increment = matrix_parse_expr (s);
4130 if (!loop->increment)
4135 if (lex_match_id (s->lexer, "IF"))
4137 loop->top_condition = matrix_parse_expr (s);
4138 if (!loop->top_condition)
4142 bool was_in_loop = s->in_loop;
4144 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4146 s->in_loop = was_in_loop;
4150 if (!lex_match_phrase (s->lexer, "END LOOP"))
4153 if (lex_match_id (s->lexer, "IF"))
4155 loop->bottom_condition = matrix_parse_expr (s);
4156 if (!loop->bottom_condition)
4163 matrix_cmd_destroy (cmd);
4168 matrix_cmd_execute_loop (struct loop_command *cmd)
4172 long int increment = 1;
4175 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4176 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4178 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4182 if (increment > 0 ? value > end
4183 : increment < 0 ? value < end
4188 int mxloops = settings_get_mxloops ();
4189 for (int i = 0; i < mxloops; i++)
4194 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4196 gsl_matrix_free (cmd->var->value);
4197 cmd->var->value = NULL;
4199 if (!cmd->var->value)
4200 cmd->var->value = gsl_matrix_alloc (1, 1);
4202 gsl_matrix_set (cmd->var->value, 0, 0, value);
4205 if (cmd->top_condition)
4208 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4213 for (size_t j = 0; j < cmd->commands.n; j++)
4214 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4217 if (cmd->bottom_condition)
4220 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4221 "END LOOP IF", &d) || d > 0)
4228 if (increment > 0 ? value > end : value < end)
4234 static struct matrix_cmd *
4235 matrix_parse_break (struct matrix_state *s)
4239 msg (SE, _("BREAK not inside LOOP."));
4243 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4244 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4248 static struct matrix_cmd *
4249 matrix_parse_display (struct matrix_state *s)
4251 bool dictionary = false;
4252 bool status = false;
4253 while (lex_token (s->lexer) == T_ID)
4255 if (lex_match_id (s->lexer, "DICTIONARY"))
4257 else if (lex_match_id (s->lexer, "STATUS"))
4261 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4265 if (!dictionary && !status)
4266 dictionary = status = true;
4268 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4269 *cmd = (struct matrix_cmd) {
4270 .type = MCMD_DISPLAY,
4273 .dictionary = dictionary,
4281 compare_matrix_var_pointers (const void *a_, const void *b_)
4283 const struct matrix_var *const *ap = a_;
4284 const struct matrix_var *const *bp = b_;
4285 const struct matrix_var *a = *ap;
4286 const struct matrix_var *b = *bp;
4287 return strcmp (a->name, b->name);
4291 matrix_cmd_execute_display (struct display_command *cmd)
4293 const struct matrix_state *s = cmd->state;
4295 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4296 pivot_dimension_create (
4297 table, PIVOT_AXIS_COLUMN, N_("Property"),
4298 N_("Rows"), N_("Columns"), N_("Size (kB)"));
4300 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4303 const struct matrix_var *var;
4304 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4306 vars[n_vars++] = var;
4307 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4309 struct pivot_dimension *rows = pivot_dimension_create (
4310 table, PIVOT_AXIS_ROW, N_("Variable"));
4311 for (size_t i = 0; i < n_vars; i++)
4313 const struct matrix_var *var = vars[i];
4314 pivot_category_create_leaf (
4315 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4317 size_t r = var->value->size1;
4318 size_t c = var->value->size2;
4319 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4320 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4321 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4324 pivot_table_submit (table);
4327 static struct matrix_cmd *
4328 matrix_parse_release (struct matrix_state *s)
4330 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4331 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4333 size_t allocated_vars = 0;
4334 while (lex_token (s->lexer) == T_ID)
4336 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4339 if (cmd->release.n_vars >= allocated_vars)
4340 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4341 sizeof *cmd->release.vars);
4342 cmd->release.vars[cmd->release.n_vars++] = var;
4345 lex_error (s->lexer, _("Variable name expected."));
4348 if (!lex_match (s->lexer, T_COMMA))
4356 matrix_cmd_execute_release (struct release_command *cmd)
4358 for (size_t i = 0; i < cmd->n_vars; i++)
4360 struct matrix_var *var = cmd->vars[i];
4361 gsl_matrix_free (var->value);
4366 static struct save_file *
4367 save_file_create (struct matrix_state *s, struct file_handle *fh,
4368 struct string_array *variables,
4369 struct matrix_expr *names,
4370 struct stringi_set *strings)
4372 for (size_t i = 0; i < s->n_save_files; i++)
4374 struct save_file *sf = s->save_files[i];
4375 if (fh_equal (sf->file, fh))
4379 string_array_destroy (variables);
4380 matrix_expr_destroy (names);
4381 stringi_set_destroy (strings);
4387 struct save_file *sf = xmalloc (sizeof *sf);
4388 *sf = (struct save_file) {
4390 .dataset = s->dataset,
4391 .variables = *variables,
4393 .strings = STRINGI_SET_INITIALIZER (sf->strings),
4396 stringi_set_swap (&sf->strings, strings);
4397 stringi_set_destroy (strings);
4399 s->save_files = xrealloc (s->save_files,
4400 (s->n_save_files + 1) * sizeof *s->save_files);
4401 s->save_files[s->n_save_files++] = sf;
4405 static struct casewriter *
4406 save_file_open (struct save_file *sf, gsl_matrix *m)
4408 if (sf->writer || sf->error)
4412 size_t n_variables = caseproto_get_n_widths (
4413 casewriter_get_proto (sf->writer));
4414 if (m->size2 != n_variables)
4416 msg (ME, _("The first SAVE to %s within this matrix program "
4417 "had %zu columns, so a %zu×%zu matrix cannot be "
4419 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
4420 n_variables, m->size1, m->size2);
4428 struct dictionary *dict = dict_create (get_default_encoding ());
4430 /* Fill 'names' with user-specified names if there were any, either from
4431 sf->variables or sf->names. */
4432 struct string_array names = { .n = 0 };
4433 if (sf->variables.n)
4434 string_array_clone (&names, &sf->variables);
4437 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
4438 if (nm && is_vector (nm))
4440 gsl_vector nv = to_vector (nm);
4441 for (size_t i = 0; i < nv.size; i++)
4443 char *name = trimmed_string (gsl_vector_get (&nv, i));
4444 if (dict_id_is_valid (dict, name, true))
4445 string_array_append_nocopy (&names, name);
4450 gsl_matrix_free (nm);
4453 struct stringi_set strings;
4454 stringi_set_clone (&strings, &sf->strings);
4456 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4461 name = names.strings[i];
4464 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4468 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4469 struct variable *var = dict_create_var (dict, name, width);
4472 msg (ME, _("Duplicate variable name %s in SAVE statement."),
4478 if (!stringi_set_is_empty (&strings))
4480 size_t count = stringi_set_count (&strings);
4481 const char *example = stringi_set_node_get_string (
4482 stringi_set_first (&strings));
4485 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
4486 "variable %s."), example);
4488 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
4489 "unknown variable: %s.",
4490 "The SAVE command STRINGS subcommand specifies %zu "
4491 "unknown variables, including %s.",
4496 stringi_set_destroy (&strings);
4497 string_array_destroy (&names);
4506 if (sf->file == fh_inline_file ())
4507 sf->writer = autopaging_writer_create (dict_get_proto (dict));
4509 sf->writer = any_writer_open (sf->file, dict);
4521 save_file_destroy (struct save_file *sf)
4525 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
4527 dataset_set_dict (sf->dataset, sf->dict);
4528 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
4532 casewriter_destroy (sf->writer);
4533 dict_unref (sf->dict);
4535 fh_unref (sf->file);
4536 string_array_destroy (&sf->variables);
4537 matrix_expr_destroy (sf->names);
4538 stringi_set_destroy (&sf->strings);
4543 static struct matrix_cmd *
4544 matrix_parse_save (struct matrix_state *s)
4546 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4547 *cmd = (struct matrix_cmd) {
4549 .save = { .expression = NULL },
4552 struct file_handle *fh = NULL;
4553 struct save_command *save = &cmd->save;
4555 struct string_array variables = STRING_ARRAY_INITIALIZER;
4556 struct matrix_expr *names = NULL;
4557 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
4559 save->expression = matrix_parse_exp (s);
4560 if (!save->expression)
4563 while (lex_match (s->lexer, T_SLASH))
4565 if (lex_match_id (s->lexer, "OUTFILE"))
4567 lex_match (s->lexer, T_EQUALS);
4570 fh = (lex_match (s->lexer, T_ASTERISK)
4572 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4576 else if (lex_match_id (s->lexer, "VARIABLES"))
4578 lex_match (s->lexer, T_EQUALS);
4582 struct dictionary *d = dict_create (get_default_encoding ());
4583 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4584 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4589 string_array_clear (&variables);
4590 variables = (struct string_array) {
4596 else if (lex_match_id (s->lexer, "NAMES"))
4598 lex_match (s->lexer, T_EQUALS);
4599 matrix_expr_destroy (names);
4600 names = matrix_parse_exp (s);
4604 else if (lex_match_id (s->lexer, "STRINGS"))
4606 lex_match (s->lexer, T_EQUALS);
4607 while (lex_token (s->lexer) == T_ID)
4609 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
4611 if (!lex_match (s->lexer, T_COMMA))
4617 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4625 if (s->prev_save_file)
4626 fh = fh_ref (s->prev_save_file);
4629 lex_sbc_missing ("OUTFILE");
4633 fh_unref (s->prev_save_file);
4634 s->prev_save_file = fh_ref (fh);
4636 if (variables.n && names)
4638 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4639 matrix_expr_destroy (names);
4643 save->sf = save_file_create (s, fh, &variables, names, &strings);
4647 string_array_destroy (&variables);
4648 matrix_expr_destroy (names);
4649 stringi_set_destroy (&strings);
4651 matrix_cmd_destroy (cmd);
4656 matrix_cmd_execute_save (const struct save_command *save)
4658 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4662 struct casewriter *writer = save_file_open (save->sf, m);
4665 gsl_matrix_free (m);
4669 const struct caseproto *proto = casewriter_get_proto (writer);
4670 for (size_t y = 0; y < m->size1; y++)
4672 struct ccase *c = case_create (proto);
4673 for (size_t x = 0; x < m->size2; x++)
4675 double d = gsl_matrix_get (m, y, x);
4676 int width = caseproto_get_width (proto, x);
4677 union value *value = case_data_rw_idx (c, x);
4681 memcpy (value->s, &d, width);
4683 casewriter_write (writer, c);
4685 gsl_matrix_free (m);
4688 static struct matrix_cmd *
4689 matrix_parse_read (struct matrix_state *s)
4691 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4692 *cmd = (struct matrix_cmd) {
4694 .read = { .format = FMT_F },
4697 struct file_handle *fh = NULL;
4698 char *encoding = NULL;
4699 struct read_command *read = &cmd->read;
4700 read->dst = matrix_lvalue_parse (s);
4705 int repetitions = 0;
4706 int record_width = 0;
4707 bool seen_format = false;
4708 while (lex_match (s->lexer, T_SLASH))
4710 if (lex_match_id (s->lexer, "FILE"))
4712 lex_match (s->lexer, T_EQUALS);
4715 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4719 else if (lex_match_id (s->lexer, "ENCODING"))
4721 lex_match (s->lexer, T_EQUALS);
4722 if (!lex_force_string (s->lexer))
4726 encoding = ss_xstrdup (lex_tokss (s->lexer));
4730 else if (lex_match_id (s->lexer, "FIELD"))
4732 lex_match (s->lexer, T_EQUALS);
4734 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4736 read->c1 = lex_integer (s->lexer);
4738 if (!lex_force_match (s->lexer, T_TO)
4739 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4741 read->c2 = lex_integer (s->lexer) + 1;
4744 record_width = read->c2 - read->c1;
4745 if (lex_match (s->lexer, T_BY))
4747 if (!lex_force_int_range (s->lexer, "BY", 1,
4748 read->c2 - read->c1))
4750 by = lex_integer (s->lexer);
4753 if (record_width % by)
4755 msg (SE, _("BY %d does not evenly divide record width %d."),
4763 else if (lex_match_id (s->lexer, "SIZE"))
4765 lex_match (s->lexer, T_EQUALS);
4766 matrix_expr_destroy (read->size);
4767 read->size = matrix_parse_exp (s);
4771 else if (lex_match_id (s->lexer, "MODE"))
4773 lex_match (s->lexer, T_EQUALS);
4774 if (lex_match_id (s->lexer, "RECTANGULAR"))
4775 read->symmetric = false;
4776 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4777 read->symmetric = true;
4780 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4784 else if (lex_match_id (s->lexer, "REREAD"))
4785 read->reread = true;
4786 else if (lex_match_id (s->lexer, "FORMAT"))
4790 lex_sbc_only_once ("FORMAT");
4795 lex_match (s->lexer, T_EQUALS);
4797 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4800 const char *p = lex_tokcstr (s->lexer);
4801 if (c_isdigit (p[0]))
4803 repetitions = atoi (p);
4804 p += strspn (p, "0123456789");
4805 if (!fmt_from_name (p, &read->format))
4807 lex_error (s->lexer, _("Unknown format %s."), p);
4812 else if (fmt_from_name (p, &read->format))
4816 struct fmt_spec format;
4817 if (!parse_format_specifier (s->lexer, &format))
4819 read->format = format.type;
4825 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4826 "REREAD", "FORMAT");
4833 lex_sbc_missing ("FIELD");
4837 if (!read->dst->n_indexes && !read->size)
4839 msg (SE, _("SIZE is required for reading data into a full matrix "
4840 "(as opposed to a submatrix)."));
4846 if (s->prev_read_file)
4847 fh = fh_ref (s->prev_read_file);
4850 lex_sbc_missing ("FILE");
4854 fh_unref (s->prev_read_file);
4855 s->prev_read_file = fh_ref (fh);
4857 read->rf = read_file_create (s, fh);
4861 free (read->rf->encoding);
4862 read->rf->encoding = encoding;
4866 /* Field width may be specified in multiple ways:
4869 2. The format on FORMAT.
4870 3. The repetition factor on FORMAT.
4872 (2) and (3) are mutually exclusive.
4874 If more than one of these is present, they must agree. If none of them is
4875 present, then free-field format is used.
4877 if (repetitions > record_width)
4879 msg (SE, _("%d repetitions cannot fit in record width %d."),
4880 repetitions, record_width);
4883 int w = (repetitions ? record_width / repetitions
4889 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4890 "which implies field width %d, "
4891 "but BY specifies field width %d."),
4892 repetitions, record_width, w, by);
4894 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4903 matrix_cmd_destroy (cmd);
4909 parse_error (const struct dfm_reader *reader, enum fmt_type format,
4910 struct substring data, size_t y, size_t x,
4911 int first_column, int last_column, char *error)
4913 int line_number = dfm_get_line_number (reader);
4914 struct msg_location *location = xmalloc (sizeof *location);
4915 *location = (struct msg_location) {
4916 .file_name = xstrdup (dfm_get_file_name (reader)),
4917 .first_line = line_number,
4918 .last_line = line_number + 1,
4919 .first_column = first_column,
4920 .last_column = last_column,
4922 struct msg *m = xmalloc (sizeof *m);
4924 .category = MSG_C_DATA,
4925 .severity = MSG_S_WARNING,
4926 .location = location,
4927 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
4928 "for matrix row %zu, column %zu: %s"),
4929 (int) data.length, data.string, fmt_name (format),
4930 y + 1, x + 1, error),
4938 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4939 gsl_matrix *m, struct substring p, size_t y, size_t x,
4940 const char *line_start)
4942 const char *input_encoding = dfm_reader_get_encoding (reader);
4945 if (fmt_is_numeric (read->format))
4948 error = data_in (p, input_encoding, read->format,
4949 settings_get_fmt_settings (), &v, 0, NULL);
4950 if (!error && v.f == SYSMIS)
4951 error = xstrdup (_("Matrix data may not contain missing value."));
4956 uint8_t s[sizeof (double)];
4957 union value v = { .s = s };
4958 error = data_in (p, input_encoding, read->format,
4959 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
4960 memcpy (&f, s, sizeof f);
4965 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
4966 int c2 = c1 + ss_utf8_count_columns (p) - 1;
4967 parse_error (reader, read->format, p, y, x, c1, c2, error);
4971 gsl_matrix_set (m, y, x, f);
4972 if (read->symmetric && x != y)
4973 gsl_matrix_set (m, x, y, f);
4978 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4979 struct substring *line, const char **startp)
4981 if (dfm_eof (reader))
4983 msg (SE, _("Unexpected end of file reading matrix data."));
4986 dfm_expand_tabs (reader);
4987 struct substring record = dfm_get_record (reader);
4988 /* XXX need to recode record into UTF-8 */
4989 *startp = record.string;
4990 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
4995 matrix_read (struct read_command *read, struct dfm_reader *reader,
4998 for (size_t y = 0; y < m->size1; y++)
5000 size_t nx = read->symmetric ? y + 1 : m->size2;
5002 struct substring line = ss_empty ();
5003 const char *line_start = line.string;
5004 for (size_t x = 0; x < nx; x++)
5011 ss_ltrim (&line, ss_cstr (" ,"));
5012 if (!ss_is_empty (line))
5014 if (!matrix_read_line (read, reader, &line, &line_start))
5016 dfm_forward_record (reader);
5019 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
5023 if (!matrix_read_line (read, reader, &line, &line_start))
5025 size_t fields_per_line = (read->c2 - read->c1) / read->w;
5026 int f = x % fields_per_line;
5027 if (f == fields_per_line - 1)
5028 dfm_forward_record (reader);
5030 p = ss_substr (line, read->w * f, read->w);
5033 matrix_read_set_field (read, reader, m, p, y, x, line_start);
5037 dfm_forward_record (reader);
5040 ss_ltrim (&line, ss_cstr (" ,"));
5041 if (!ss_is_empty (line))
5044 msg (SW, _("Trailing garbage on line \"%.*s\""),
5045 (int) line.length, line.string);
5052 matrix_cmd_execute_read (struct read_command *read)
5054 struct index_vector iv0, iv1;
5055 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
5058 size_t size[2] = { SIZE_MAX, SIZE_MAX };
5061 gsl_matrix *m = matrix_expr_evaluate (read->size);
5067 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5068 "not a %zu×%zu matrix."), m->size1, m->size2);
5069 gsl_matrix_free (m);
5075 gsl_vector v = to_vector (m);
5079 d[0] = gsl_vector_get (&v, 0);
5082 else if (v.size == 2)
5084 d[0] = gsl_vector_get (&v, 0);
5085 d[1] = gsl_vector_get (&v, 1);
5089 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5090 "not a %zu×%zu matrix."),
5091 m->size1, m->size2),
5092 gsl_matrix_free (m);
5097 gsl_matrix_free (m);
5099 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
5101 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
5102 "are outside valid range."),
5113 if (read->dst->n_indexes)
5115 size_t submatrix_size[2];
5116 if (read->dst->n_indexes == 2)
5118 submatrix_size[0] = iv0.n;
5119 submatrix_size[1] = iv1.n;
5121 else if (read->dst->var->value->size1 == 1)
5123 submatrix_size[0] = 1;
5124 submatrix_size[1] = iv0.n;
5128 submatrix_size[0] = iv0.n;
5129 submatrix_size[1] = 1;
5134 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
5136 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
5137 "differ from submatrix dimensions %zu×%zu."),
5139 submatrix_size[0], submatrix_size[1]);
5147 size[0] = submatrix_size[0];
5148 size[1] = submatrix_size[1];
5152 struct dfm_reader *reader = read_file_open (read->rf);
5154 dfm_reread_record (reader, 1);
5156 if (read->symmetric && size[0] != size[1])
5158 msg (SE, _("Cannot read non-square %zu×%zu matrix "
5159 "using READ with MODE=SYMMETRIC."),
5165 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
5166 matrix_read (read, reader, tmp);
5167 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
5170 static struct matrix_cmd *
5171 matrix_parse_write (struct matrix_state *s)
5173 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5174 *cmd = (struct matrix_cmd) {
5178 struct file_handle *fh = NULL;
5179 char *encoding = NULL;
5180 struct write_command *write = &cmd->write;
5181 write->expression = matrix_parse_exp (s);
5182 if (!write->expression)
5186 int repetitions = 0;
5187 int record_width = 0;
5188 enum fmt_type format = FMT_F;
5189 bool has_format = false;
5190 while (lex_match (s->lexer, T_SLASH))
5192 if (lex_match_id (s->lexer, "OUTFILE"))
5194 lex_match (s->lexer, T_EQUALS);
5197 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5201 else if (lex_match_id (s->lexer, "ENCODING"))
5203 lex_match (s->lexer, T_EQUALS);
5204 if (!lex_force_string (s->lexer))
5208 encoding = ss_xstrdup (lex_tokss (s->lexer));
5212 else if (lex_match_id (s->lexer, "FIELD"))
5214 lex_match (s->lexer, T_EQUALS);
5216 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5218 write->c1 = lex_integer (s->lexer);
5220 if (!lex_force_match (s->lexer, T_TO)
5221 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5223 write->c2 = lex_integer (s->lexer) + 1;
5226 record_width = write->c2 - write->c1;
5227 if (lex_match (s->lexer, T_BY))
5229 if (!lex_force_int_range (s->lexer, "BY", 1,
5230 write->c2 - write->c1))
5232 by = lex_integer (s->lexer);
5235 if (record_width % by)
5237 msg (SE, _("BY %d does not evenly divide record width %d."),
5245 else if (lex_match_id (s->lexer, "MODE"))
5247 lex_match (s->lexer, T_EQUALS);
5248 if (lex_match_id (s->lexer, "RECTANGULAR"))
5249 write->triangular = false;
5250 else if (lex_match_id (s->lexer, "TRIANGULAR"))
5251 write->triangular = true;
5254 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5258 else if (lex_match_id (s->lexer, "HOLD"))
5260 else if (lex_match_id (s->lexer, "FORMAT"))
5262 if (has_format || write->format)
5264 lex_sbc_only_once ("FORMAT");
5268 lex_match (s->lexer, T_EQUALS);
5270 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5273 const char *p = lex_tokcstr (s->lexer);
5274 if (c_isdigit (p[0]))
5276 repetitions = atoi (p);
5277 p += strspn (p, "0123456789");
5278 if (!fmt_from_name (p, &format))
5280 lex_error (s->lexer, _("Unknown format %s."), p);
5286 else if (fmt_from_name (p, &format))
5293 struct fmt_spec spec;
5294 if (!parse_format_specifier (s->lexer, &spec))
5296 write->format = xmemdup (&spec, sizeof spec);
5301 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5309 lex_sbc_missing ("FIELD");
5315 if (s->prev_write_file)
5316 fh = fh_ref (s->prev_write_file);
5319 lex_sbc_missing ("OUTFILE");
5323 fh_unref (s->prev_write_file);
5324 s->prev_write_file = fh_ref (fh);
5326 write->wf = write_file_create (s, fh);
5330 free (write->wf->encoding);
5331 write->wf->encoding = encoding;
5335 /* Field width may be specified in multiple ways:
5338 2. The format on FORMAT.
5339 3. The repetition factor on FORMAT.
5341 (2) and (3) are mutually exclusive.
5343 If more than one of these is present, they must agree. If none of them is
5344 present, then free-field format is used.
5346 if (repetitions > record_width)
5348 msg (SE, _("%d repetitions cannot fit in record width %d."),
5349 repetitions, record_width);
5352 int w = (repetitions ? record_width / repetitions
5353 : write->format ? write->format->w
5358 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5359 "which implies field width %d, "
5360 "but BY specifies field width %d."),
5361 repetitions, record_width, w, by);
5363 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5367 if (w && !write->format)
5369 write->format = xmalloc (sizeof *write->format);
5370 *write->format = (struct fmt_spec) { .type = format, .w = w };
5372 if (!fmt_check_output (write->format))
5376 if (write->format && fmt_var_width (write->format) > sizeof (double))
5378 char s[FMT_STRING_LEN_MAX + 1];
5379 fmt_to_string (write->format, s);
5380 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
5381 s, sizeof (double));
5389 matrix_cmd_destroy (cmd);
5394 matrix_cmd_execute_write (struct write_command *write)
5396 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5400 if (write->triangular && m->size1 != m->size2)
5402 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5403 "the matrix to be written has dimensions %zu×%zu."),
5404 m->size1, m->size2);
5405 gsl_matrix_free (m);
5409 struct dfm_writer *writer = write_file_open (write->wf);
5410 if (!writer || !m->size1)
5412 gsl_matrix_free (m);
5416 const struct fmt_settings *settings = settings_get_fmt_settings ();
5417 struct u8_line *line = write->wf->held;
5418 for (size_t y = 0; y < m->size1; y++)
5422 line = xmalloc (sizeof *line);
5423 u8_line_init (line);
5425 size_t nx = write->triangular ? y + 1 : m->size2;
5427 for (size_t x = 0; x < nx; x++)
5430 double f = gsl_matrix_get (m, y, x);
5434 if (fmt_is_numeric (write->format->type))
5437 v.s = (uint8_t *) &f;
5438 s = data_out (&v, NULL, write->format, settings);
5442 s = xmalloc (DBL_BUFSIZE_BOUND);
5443 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
5444 >= DBL_BUFSIZE_BOUND)
5447 size_t len = strlen (s);
5448 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5449 if (width + x0 > write->c2)
5451 dfm_put_record_utf8 (writer, line->s.ss.string,
5453 u8_line_clear (line);
5456 u8_line_put (line, x0, x0 + width, s, len);
5459 x0 += write->format ? write->format->w : width + 1;
5462 if (y + 1 >= m->size1 && write->hold)
5464 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5465 u8_line_clear (line);
5469 u8_line_destroy (line);
5473 write->wf->held = line;
5475 gsl_matrix_free (m);
5478 static struct matrix_cmd *
5479 matrix_parse_get (struct matrix_state *s)
5481 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5482 *cmd = (struct matrix_cmd) {
5485 .dataset = s->dataset,
5486 .user = { .treatment = MGET_ERROR },
5487 .system = { .treatment = MGET_ERROR },
5491 struct get_command *get = &cmd->get;
5492 get->dst = matrix_lvalue_parse (s);
5496 while (lex_match (s->lexer, T_SLASH))
5498 if (lex_match_id (s->lexer, "FILE"))
5500 lex_match (s->lexer, T_EQUALS);
5502 fh_unref (get->file);
5503 if (lex_match (s->lexer, T_ASTERISK))
5507 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5512 else if (lex_match_id (s->lexer, "ENCODING"))
5514 lex_match (s->lexer, T_EQUALS);
5515 if (!lex_force_string (s->lexer))
5518 free (get->encoding);
5519 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5523 else if (lex_match_id (s->lexer, "VARIABLES"))
5525 lex_match (s->lexer, T_EQUALS);
5529 lex_sbc_only_once ("VARIABLES");
5533 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
5536 else if (lex_match_id (s->lexer, "NAMES"))
5538 lex_match (s->lexer, T_EQUALS);
5539 if (!lex_force_id (s->lexer))
5542 struct substring name = lex_tokss (s->lexer);
5543 get->names = matrix_var_lookup (s, name);
5545 get->names = matrix_var_create (s, name);
5548 else if (lex_match_id (s->lexer, "MISSING"))
5550 lex_match (s->lexer, T_EQUALS);
5551 if (lex_match_id (s->lexer, "ACCEPT"))
5552 get->user.treatment = MGET_ACCEPT;
5553 else if (lex_match_id (s->lexer, "OMIT"))
5554 get->user.treatment = MGET_OMIT;
5555 else if (lex_is_number (s->lexer))
5557 get->user.treatment = MGET_RECODE;
5558 get->user.substitute = lex_number (s->lexer);
5563 lex_error (s->lexer, NULL);
5567 else if (lex_match_id (s->lexer, "SYSMIS"))
5569 lex_match (s->lexer, T_EQUALS);
5570 if (lex_match_id (s->lexer, "OMIT"))
5571 get->system.treatment = MGET_OMIT;
5572 else if (lex_is_number (s->lexer))
5574 get->system.treatment = MGET_RECODE;
5575 get->system.substitute = lex_number (s->lexer);
5580 lex_error (s->lexer, NULL);
5586 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5587 "MISSING", "SYSMIS");
5592 if (get->user.treatment != MGET_ACCEPT)
5593 get->system.treatment = MGET_ERROR;
5598 matrix_cmd_destroy (cmd);
5603 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
5604 const struct dictionary *dict)
5606 struct variable **vars;
5611 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
5612 &vars, &n_vars, PV_NUMERIC))
5617 n_vars = dict_get_var_cnt (dict);
5618 vars = xnmalloc (n_vars, sizeof *vars);
5619 for (size_t i = 0; i < n_vars; i++)
5621 struct variable *var = dict_get_var (dict, i);
5622 if (!var_is_numeric (var))
5624 msg (SE, _("GET: Variable %s is not numeric."),
5625 var_get_name (var));
5635 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
5636 for (size_t i = 0; i < n_vars; i++)
5638 char s[sizeof (double)];
5640 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
5641 memcpy (&f, s, sizeof f);
5642 gsl_matrix_set (names, i, 0, f);
5645 gsl_matrix_free (get->names->value);
5646 get->names->value = names;
5650 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5651 long long int casenum = 1;
5653 for (struct ccase *c = casereader_read (reader); c;
5654 c = casereader_read (reader), casenum++)
5656 if (n_rows >= m->size1)
5658 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5659 for (size_t y = 0; y < n_rows; y++)
5660 for (size_t x = 0; x < n_vars; x++)
5661 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5662 gsl_matrix_free (m);
5667 for (size_t x = 0; x < n_vars; x++)
5669 const struct variable *var = vars[x];
5670 double d = case_num (c, var);
5673 if (get->system.treatment == MGET_RECODE)
5674 d = get->system.substitute;
5675 else if (get->system.treatment == MGET_OMIT)
5679 msg (SE, _("GET: Variable %s in case %lld "
5680 "is system-missing."),
5681 var_get_name (var), casenum);
5685 else if (var_is_num_missing (var, d, MV_USER))
5687 if (get->user.treatment == MGET_RECODE)
5688 d = get->user.substitute;
5689 else if (get->user.treatment == MGET_OMIT)
5691 else if (get->user.treatment != MGET_ACCEPT)
5693 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5695 var_get_name (var), casenum, d);
5699 gsl_matrix_set (m, n_rows, x, d);
5710 matrix_lvalue_evaluate_and_assign (get->dst, m);
5713 gsl_matrix_free (m);
5718 matrix_open_casereader (const char *command_name,
5719 struct file_handle *file, const char *encoding,
5720 struct dataset *dataset,
5721 struct casereader **readerp, struct dictionary **dictp)
5725 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
5726 return *readerp != NULL;
5730 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
5732 msg (ME, _("The %s command cannot read empty active file."),
5736 *readerp = proc_open (dataset);
5737 *dictp = dict_ref (dataset_dict (dataset));
5743 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
5744 struct casereader *reader, struct dictionary *dict)
5747 casereader_destroy (reader);
5749 proc_commit (dataset);
5753 matrix_cmd_execute_get (struct get_command *get)
5755 struct casereader *r;
5756 struct dictionary *d;
5757 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
5760 matrix_cmd_execute_get__ (get, r, d);
5761 matrix_close_casereader (get->file, get->dataset, r, d);
5766 match_rowtype (struct lexer *lexer)
5768 static const char *rowtypes[] = {
5769 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5771 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5773 for (size_t i = 0; i < n_rowtypes; i++)
5774 if (lex_match_id (lexer, rowtypes[i]))
5777 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5782 parse_var_names (struct lexer *lexer, struct string_array *sa)
5784 lex_match (lexer, T_EQUALS);
5786 string_array_clear (sa);
5788 struct dictionary *dict = dict_create (get_default_encoding ());
5789 dict_create_var_assert (dict, "ROWTYPE_", 8);
5790 dict_create_var_assert (dict, "VARNAME_", 8);
5793 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5794 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5799 string_array_clear (sa);
5800 sa->strings = names;
5801 sa->n = sa->allocated = n_names;
5807 msave_common_uninit (struct msave_common *common)
5811 fh_unref (common->outfile);
5812 string_array_destroy (&common->variables);
5813 string_array_destroy (&common->fnames);
5814 string_array_destroy (&common->snames);
5819 compare_variables (const char *keyword,
5820 const struct string_array *new,
5821 const struct string_array *old)
5827 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5828 "on the first MSAVE within MATRIX."), keyword);
5831 else if (!string_array_equal_case (old, new))
5833 msg (SE, _("%s must specify the same variables each time within "
5834 "a given MATRIX."), keyword);
5841 static struct matrix_cmd *
5842 matrix_parse_msave (struct matrix_state *s)
5844 struct msave_common common = { .outfile = NULL };
5845 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5846 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5848 struct msave_command *msave = &cmd->msave;
5849 if (lex_token (s->lexer) == T_ID)
5850 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5851 msave->expr = matrix_parse_exp (s);
5855 while (lex_match (s->lexer, T_SLASH))
5857 if (lex_match_id (s->lexer, "TYPE"))
5859 lex_match (s->lexer, T_EQUALS);
5861 msave->rowtype = match_rowtype (s->lexer);
5862 if (!msave->rowtype)
5865 else if (lex_match_id (s->lexer, "OUTFILE"))
5867 lex_match (s->lexer, T_EQUALS);
5869 fh_unref (common.outfile);
5870 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5871 if (!common.outfile)
5874 else if (lex_match_id (s->lexer, "VARIABLES"))
5876 if (!parse_var_names (s->lexer, &common.variables))
5879 else if (lex_match_id (s->lexer, "FNAMES"))
5881 if (!parse_var_names (s->lexer, &common.fnames))
5884 else if (lex_match_id (s->lexer, "SNAMES"))
5886 if (!parse_var_names (s->lexer, &common.snames))
5889 else if (lex_match_id (s->lexer, "SPLIT"))
5891 lex_match (s->lexer, T_EQUALS);
5893 matrix_expr_destroy (msave->splits);
5894 msave->splits = matrix_parse_exp (s);
5898 else if (lex_match_id (s->lexer, "FACTOR"))
5900 lex_match (s->lexer, T_EQUALS);
5902 matrix_expr_destroy (msave->factors);
5903 msave->factors = matrix_parse_exp (s);
5904 if (!msave->factors)
5909 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5910 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5914 if (!msave->rowtype)
5916 lex_sbc_missing ("TYPE");
5919 common.has_splits = msave->splits || common.snames.n;
5920 common.has_factors = msave->factors || common.fnames.n;
5922 struct msave_common *c = s->common ? s->common : &common;
5923 if (c->fnames.n && !msave->factors)
5925 msg (SE, _("FNAMES requires FACTOR."));
5928 if (c->snames.n && !msave->splits)
5930 msg (SE, _("SNAMES requires SPLIT."));
5933 if (c->has_factors && !common.has_factors)
5935 msg (SE, _("%s is required because it was present on the first "
5936 "MSAVE in this MATRIX command."), "FACTOR");
5939 if (c->has_splits && !common.has_splits)
5941 msg (SE, _("%s is required because it was present on the first "
5942 "MSAVE in this MATRIX command."), "SPLIT");
5948 if (!common.outfile)
5950 lex_sbc_missing ("OUTFILE");
5953 s->common = xmemdup (&common, sizeof common);
5959 bool same = common.outfile == s->common->outfile;
5960 fh_unref (common.outfile);
5963 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5964 "within a single MATRIX command."));
5968 if (!compare_variables ("VARIABLES",
5969 &common.variables, &s->common->variables)
5970 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5971 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5973 msave_common_uninit (&common);
5975 msave->common = s->common;
5976 if (!msave->varname_)
5977 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5981 msave_common_uninit (&common);
5982 matrix_cmd_destroy (cmd);
5987 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5989 gsl_matrix *m = matrix_expr_evaluate (e);
5995 msg (SE, _("%s expression must evaluate to vector, "
5996 "not a %zu×%zu matrix."),
5997 name, m->size1, m->size2);
5998 gsl_matrix_free (m);
6002 return matrix_to_vector (m);
6006 msave_add_vars (struct dictionary *d, const struct string_array *vars)
6008 for (size_t i = 0; i < vars->n; i++)
6009 if (!dict_create_var (d, vars->strings[i], 0))
6010 return vars->strings[i];
6014 static struct dictionary *
6015 msave_create_dict (const struct msave_common *common)
6017 struct dictionary *dict = dict_create (get_default_encoding ());
6019 const char *dup_split = msave_add_vars (dict, &common->snames);
6022 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
6026 dict_create_var_assert (dict, "ROWTYPE_", 8);
6028 const char *dup_factor = msave_add_vars (dict, &common->fnames);
6031 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
6035 dict_create_var_assert (dict, "VARNAME_", 8);
6037 const char *dup_var = msave_add_vars (dict, &common->variables);
6040 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
6052 matrix_cmd_execute_msave (struct msave_command *msave)
6054 struct msave_common *common = msave->common;
6055 gsl_matrix *m = NULL;
6056 gsl_vector *factors = NULL;
6057 gsl_vector *splits = NULL;
6059 m = matrix_expr_evaluate (msave->expr);
6063 if (!common->variables.n)
6064 for (size_t i = 0; i < m->size2; i++)
6065 string_array_append_nocopy (&common->variables,
6066 xasprintf ("COL%zu", i + 1));
6068 if (m->size2 != common->variables.n)
6070 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
6071 m->size2, common->variables.n);
6077 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
6081 if (!common->fnames.n)
6082 for (size_t i = 0; i < factors->size; i++)
6083 string_array_append_nocopy (&common->fnames,
6084 xasprintf ("FAC%zu", i + 1));
6086 if (factors->size != common->fnames.n)
6088 msg (SE, _("There are %zu factor variables, "
6089 "but %zu split values were supplied."),
6090 common->fnames.n, factors->size);
6097 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
6101 if (!common->fnames.n)
6102 for (size_t i = 0; i < splits->size; i++)
6103 string_array_append_nocopy (&common->fnames,
6104 xasprintf ("SPL%zu", i + 1));
6106 if (splits->size != common->fnames.n)
6108 msg (SE, _("There are %zu split variables, "
6109 "but %zu split values were supplied."),
6110 common->fnames.n, splits->size);
6115 if (!common->writer)
6117 struct dictionary *dict = msave_create_dict (common);
6121 common->writer = any_writer_open (common->outfile, dict);
6122 if (!common->writer)
6128 common->dict = dict;
6131 for (size_t y = 0; y < m->size1; y++)
6133 struct ccase *c = case_create (dict_get_proto (common->dict));
6136 /* Split variables */
6138 for (size_t i = 0; i < splits->size; i++)
6139 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
6142 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6143 msave->rowtype, ' ');
6147 for (size_t i = 0; i < factors->size; i++)
6148 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
6151 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6152 msave->varname_, ' ');
6154 /* Continuous variables. */
6155 for (size_t x = 0; x < m->size2; x++)
6156 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
6157 casewriter_write (common->writer, c);
6161 gsl_matrix_free (m);
6162 gsl_vector_free (factors);
6163 gsl_vector_free (splits);
6166 static struct matrix_cmd *
6167 matrix_parse_mget (struct matrix_state *s)
6169 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6170 *cmd = (struct matrix_cmd) {
6174 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
6178 struct mget_command *mget = &cmd->mget;
6180 while (lex_token (s->lexer) != T_ENDCMD)
6182 if (lex_match_id (s->lexer, "FILE"))
6184 lex_match (s->lexer, T_EQUALS);
6186 fh_unref (mget->file);
6187 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6191 else if (lex_match_id (s->lexer, "ENCODING"))
6193 lex_match (s->lexer, T_EQUALS);
6194 if (!lex_force_string (s->lexer))
6197 free (mget->encoding);
6198 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
6202 else if (lex_match_id (s->lexer, "TYPE"))
6204 lex_match (s->lexer, T_EQUALS);
6205 while (lex_token (s->lexer) != T_SLASH
6206 && lex_token (s->lexer) != T_ENDCMD)
6208 const char *rowtype = match_rowtype (s->lexer);
6212 stringi_set_insert (&mget->rowtypes, rowtype);
6217 lex_error_expecting (s->lexer, "FILE", "TYPE");
6220 lex_match (s->lexer, T_SLASH);
6225 matrix_cmd_destroy (cmd);
6229 static const struct variable *
6230 get_a8_var (const struct dictionary *d, const char *name)
6232 const struct variable *v = dict_lookup_var (d, name);
6235 msg (SE, _("Matrix data file lacks %s variable."), name);
6238 if (var_get_width (v) != 8)
6240 msg (SE, _("%s variable in matrix data file must be "
6241 "8-byte string, but it has width %d."),
6242 name, var_get_width (v));
6249 is_rowtype (const union value *v, const char *rowtype)
6251 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
6252 ss_rtrim (&vs, ss_cstr (" "));
6253 return ss_equals_case (vs, ss_cstr (rowtype));
6257 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6258 const struct dictionary *d,
6259 const struct variable *rowtype_var,
6260 struct matrix_state *s, size_t si, size_t fi,
6261 size_t cs, size_t cn,
6262 struct pivot_table *pt,
6263 struct pivot_dimension *var_dimension)
6268 const union value *rowtype_ = case_data (rows[0], rowtype_var);
6269 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
6270 : is_rowtype (rowtype_, "CORR") ? "CR"
6271 : is_rowtype (rowtype_, "MEAN") ? "MN"
6272 : is_rowtype (rowtype_, "STDDEV") ? "SD"
6273 : is_rowtype (rowtype_, "N") ? "NC"
6276 struct string name = DS_EMPTY_INITIALIZER;
6277 ds_put_cstr (&name, name_prefix);
6279 ds_put_format (&name, "F%zu", fi);
6281 ds_put_format (&name, "S%zu", si);
6283 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6285 mv = matrix_var_create (s, ds_ss (&name));
6288 msg (SW, _("Matrix data file contains variable with existing name %s."),
6293 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6294 size_t n_missing = 0;
6295 for (size_t y = 0; y < n_rows; y++)
6297 for (size_t x = 0; x < cn; x++)
6299 struct variable *var = dict_get_var (d, cs + x);
6300 double value = case_num (rows[y], var);
6301 if (var_is_num_missing (var, value, MV_ANY))
6306 gsl_matrix_set (m, y, x, value);
6310 int var_index = pivot_category_create_leaf (
6311 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
6312 double values[] = { n_rows, cn };
6313 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
6314 pivot_table_put2 (pt, j, var_index, pivot_value_new_integer (values[j]));
6317 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6318 "value, which was treated as zero.",
6319 "Matrix data file variable %s contains %zu missing "
6320 "values, which were treated as zero.", n_missing),
6321 ds_cstr (&name), n_missing);
6326 for (size_t y = 0; y < n_rows; y++)
6327 case_unref (rows[y]);
6331 var_changed (const struct ccase *ca, const struct ccase *cb,
6332 const struct variable *var)
6335 ? !value_equal (case_data (ca, var), case_data (cb, var),
6336 var_get_width (var))
6341 vars_changed (const struct ccase *ca, const struct ccase *cb,
6342 const struct dictionary *d,
6343 size_t first_var, size_t n_vars)
6345 for (size_t i = 0; i < n_vars; i++)
6347 const struct variable *v = dict_get_var (d, first_var + i);
6348 if (var_changed (ca, cb, v))
6355 matrix_cmd_execute_mget__ (struct mget_command *mget,
6356 struct casereader *r, const struct dictionary *d,
6357 struct pivot_table *pt,
6358 struct pivot_dimension *var_dimension)
6360 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6361 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6362 if (!rowtype_ || !varname_)
6365 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6367 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6370 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6372 msg (SE, _("Matrix data file contains no continuous variables."));
6376 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6378 const struct variable *v = dict_get_var (d, i);
6379 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6382 _("Matrix data file contains unexpected string variable %s."),
6388 /* SPLIT variables. */
6390 size_t sn = var_get_dict_index (rowtype_);
6391 struct ccase *sc = NULL;
6394 /* FACTOR variables. */
6395 size_t fs = var_get_dict_index (rowtype_) + 1;
6396 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6397 struct ccase *fc = NULL;
6400 /* Continuous variables. */
6401 size_t cs = var_get_dict_index (varname_) + 1;
6402 size_t cn = dict_get_var_cnt (d) - cs;
6403 struct ccase *cc = NULL;
6406 struct ccase **rows = NULL;
6407 size_t allocated_rows = 0;
6411 while ((c = casereader_read (r)) != NULL)
6421 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
6422 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
6423 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
6426 if (change != NOTHING_CHANGED)
6428 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6429 mget->state, si, fi, cs, cn,
6436 if (n_rows >= allocated_rows)
6437 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6440 if (change == SPLITS_CHANGED)
6446 /* Reset the factor number, if there are factors. */
6454 else if (change == FACTORS_CHANGED)
6461 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6462 mget->state, si, fi, cs, cn,
6472 matrix_cmd_execute_mget (struct mget_command *mget)
6474 struct casereader *r;
6475 struct dictionary *d;
6476 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
6477 mget->state->dataset, &r, &d))
6479 struct pivot_table *pt = pivot_table_create (
6480 N_("Matrix Variables Created by MGET"));
6481 pivot_dimension_create (pt, PIVOT_AXIS_COLUMN, N_("Dimension"),
6482 N_("Rows"), N_("Columns"));
6483 struct pivot_dimension *var_dimension = pivot_dimension_create (
6484 pt, PIVOT_AXIS_ROW, N_("Variable"));
6485 matrix_cmd_execute_mget__ (mget, r, d, pt, var_dimension);
6486 if (var_dimension->n_leaves)
6487 pivot_table_submit (pt);
6489 pivot_table_unref (pt);
6491 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
6496 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6498 if (!lex_force_id (s->lexer))
6501 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6503 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6508 static struct matrix_cmd *
6509 matrix_parse_eigen (struct matrix_state *s)
6511 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6512 *cmd = (struct matrix_cmd) {
6514 .eigen = { .expr = NULL }
6517 struct eigen_command *eigen = &cmd->eigen;
6518 if (!lex_force_match (s->lexer, T_LPAREN))
6520 eigen->expr = matrix_parse_expr (s);
6522 || !lex_force_match (s->lexer, T_COMMA)
6523 || !matrix_parse_dst_var (s, &eigen->evec)
6524 || !lex_force_match (s->lexer, T_COMMA)
6525 || !matrix_parse_dst_var (s, &eigen->eval)
6526 || !lex_force_match (s->lexer, T_RPAREN))
6532 matrix_cmd_destroy (cmd);
6537 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6539 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6542 if (!is_symmetric (A))
6544 msg (SE, _("Argument of EIGEN must be symmetric."));
6545 gsl_matrix_free (A);
6549 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6550 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6551 gsl_vector v_eval = to_vector (eval);
6552 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6553 gsl_eigen_symmv (A, &v_eval, evec, w);
6554 gsl_eigen_symmv_free (w);
6556 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6558 gsl_matrix_free (A);
6560 gsl_matrix_free (eigen->eval->value);
6561 eigen->eval->value = eval;
6563 gsl_matrix_free (eigen->evec->value);
6564 eigen->evec->value = evec;
6567 static struct matrix_cmd *
6568 matrix_parse_setdiag (struct matrix_state *s)
6570 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6571 *cmd = (struct matrix_cmd) {
6572 .type = MCMD_SETDIAG,
6573 .setdiag = { .dst = NULL }
6576 struct setdiag_command *setdiag = &cmd->setdiag;
6577 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6580 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6583 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6588 if (!lex_force_match (s->lexer, T_COMMA))
6591 setdiag->expr = matrix_parse_expr (s);
6595 if (!lex_force_match (s->lexer, T_RPAREN))
6601 matrix_cmd_destroy (cmd);
6606 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6608 gsl_matrix *dst = setdiag->dst->value;
6611 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6612 setdiag->dst->name);
6616 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6620 size_t n = MIN (dst->size1, dst->size2);
6621 if (is_scalar (src))
6623 double d = to_scalar (src);
6624 for (size_t i = 0; i < n; i++)
6625 gsl_matrix_set (dst, i, i, d);
6627 else if (is_vector (src))
6629 gsl_vector v = to_vector (src);
6630 for (size_t i = 0; i < n && i < v.size; i++)
6631 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6634 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6635 "not a %zu×%zu matrix."),
6636 src->size1, src->size2);
6637 gsl_matrix_free (src);
6640 static struct matrix_cmd *
6641 matrix_parse_svd (struct matrix_state *s)
6643 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6644 *cmd = (struct matrix_cmd) {
6646 .svd = { .expr = NULL }
6649 struct svd_command *svd = &cmd->svd;
6650 if (!lex_force_match (s->lexer, T_LPAREN))
6652 svd->expr = matrix_parse_expr (s);
6654 || !lex_force_match (s->lexer, T_COMMA)
6655 || !matrix_parse_dst_var (s, &svd->u)
6656 || !lex_force_match (s->lexer, T_COMMA)
6657 || !matrix_parse_dst_var (s, &svd->s)
6658 || !lex_force_match (s->lexer, T_COMMA)
6659 || !matrix_parse_dst_var (s, &svd->v)
6660 || !lex_force_match (s->lexer, T_RPAREN))
6666 matrix_cmd_destroy (cmd);
6671 matrix_cmd_execute_svd (struct svd_command *svd)
6673 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6677 if (m->size1 >= m->size2)
6680 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6681 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6682 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6683 gsl_vector *work = gsl_vector_alloc (A->size2);
6684 gsl_linalg_SV_decomp (A, V, &Sv, work);
6685 gsl_vector_free (work);
6687 matrix_var_set (svd->u, A);
6688 matrix_var_set (svd->s, S);
6689 matrix_var_set (svd->v, V);
6693 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6694 gsl_matrix_transpose_memcpy (At, m);
6695 gsl_matrix_free (m);
6697 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6698 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6699 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6700 gsl_vector *work = gsl_vector_alloc (At->size2);
6701 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6702 gsl_vector_free (work);
6704 matrix_var_set (svd->v, At);
6705 matrix_var_set (svd->s, St);
6706 matrix_var_set (svd->u, Vt);
6711 matrix_cmd_execute (struct matrix_cmd *cmd)
6716 matrix_cmd_execute_compute (&cmd->compute);
6720 matrix_cmd_execute_print (&cmd->print);
6724 return matrix_cmd_execute_do_if (&cmd->do_if);
6727 matrix_cmd_execute_loop (&cmd->loop);
6734 matrix_cmd_execute_display (&cmd->display);
6738 matrix_cmd_execute_release (&cmd->release);
6742 matrix_cmd_execute_save (&cmd->save);
6746 matrix_cmd_execute_read (&cmd->read);
6750 matrix_cmd_execute_write (&cmd->write);
6754 matrix_cmd_execute_get (&cmd->get);
6758 matrix_cmd_execute_msave (&cmd->msave);
6762 matrix_cmd_execute_mget (&cmd->mget);
6766 matrix_cmd_execute_eigen (&cmd->eigen);
6770 matrix_cmd_execute_setdiag (&cmd->setdiag);
6774 matrix_cmd_execute_svd (&cmd->svd);
6782 matrix_cmds_uninit (struct matrix_cmds *cmds)
6784 for (size_t i = 0; i < cmds->n; i++)
6785 matrix_cmd_destroy (cmds->commands[i]);
6786 free (cmds->commands);
6790 matrix_cmd_destroy (struct matrix_cmd *cmd)
6798 matrix_lvalue_destroy (cmd->compute.lvalue);
6799 matrix_expr_destroy (cmd->compute.rvalue);
6803 matrix_expr_destroy (cmd->print.expression);
6804 free (cmd->print.title);
6805 print_labels_destroy (cmd->print.rlabels);
6806 print_labels_destroy (cmd->print.clabels);
6810 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6812 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6813 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6815 free (cmd->do_if.clauses);
6819 matrix_expr_destroy (cmd->loop.start);
6820 matrix_expr_destroy (cmd->loop.end);
6821 matrix_expr_destroy (cmd->loop.increment);
6822 matrix_expr_destroy (cmd->loop.top_condition);
6823 matrix_expr_destroy (cmd->loop.bottom_condition);
6824 matrix_cmds_uninit (&cmd->loop.commands);
6834 free (cmd->release.vars);
6838 matrix_expr_destroy (cmd->save.expression);
6842 matrix_lvalue_destroy (cmd->read.dst);
6843 matrix_expr_destroy (cmd->read.size);
6847 matrix_expr_destroy (cmd->write.expression);
6848 free (cmd->write.format);
6852 matrix_lvalue_destroy (cmd->get.dst);
6853 fh_unref (cmd->get.file);
6854 free (cmd->get.encoding);
6855 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
6859 free (cmd->msave.varname_);
6860 matrix_expr_destroy (cmd->msave.expr);
6861 matrix_expr_destroy (cmd->msave.factors);
6862 matrix_expr_destroy (cmd->msave.splits);
6866 fh_unref (cmd->mget.file);
6867 stringi_set_destroy (&cmd->mget.rowtypes);
6871 matrix_expr_destroy (cmd->eigen.expr);
6875 matrix_expr_destroy (cmd->setdiag.expr);
6879 matrix_expr_destroy (cmd->svd.expr);
6885 struct matrix_command_name
6888 struct matrix_cmd *(*parse) (struct matrix_state *);
6891 static const struct matrix_command_name *
6892 matrix_parse_command_name (struct lexer *lexer)
6894 static const struct matrix_command_name commands[] = {
6895 { "COMPUTE", matrix_parse_compute },
6896 { "CALL EIGEN", matrix_parse_eigen },
6897 { "CALL SETDIAG", matrix_parse_setdiag },
6898 { "CALL SVD", matrix_parse_svd },
6899 { "PRINT", matrix_parse_print },
6900 { "DO IF", matrix_parse_do_if },
6901 { "LOOP", matrix_parse_loop },
6902 { "BREAK", matrix_parse_break },
6903 { "READ", matrix_parse_read },
6904 { "WRITE", matrix_parse_write },
6905 { "GET", matrix_parse_get },
6906 { "SAVE", matrix_parse_save },
6907 { "MGET", matrix_parse_mget },
6908 { "MSAVE", matrix_parse_msave },
6909 { "DISPLAY", matrix_parse_display },
6910 { "RELEASE", matrix_parse_release },
6912 static size_t n = sizeof commands / sizeof *commands;
6914 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6915 if (lex_match_phrase (lexer, c->name))
6920 static struct matrix_cmd *
6921 matrix_parse_command (struct matrix_state *s)
6923 size_t nesting_level = SIZE_MAX;
6925 struct matrix_cmd *c = NULL;
6926 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6928 lex_error (s->lexer, _("Unknown matrix command."));
6929 else if (!cmd->parse)
6930 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6934 nesting_level = output_open_group (
6935 group_item_create_nocopy (utf8_to_title (cmd->name),
6936 utf8_to_title (cmd->name)));
6941 lex_end_of_command (s->lexer);
6942 lex_discard_rest_of_command (s->lexer);
6943 if (nesting_level != SIZE_MAX)
6944 output_close_groups (nesting_level);
6950 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6952 if (!lex_force_match (lexer, T_ENDCMD))
6955 struct matrix_state state = {
6957 .session = dataset_session (ds),
6959 .vars = HMAP_INITIALIZER (state.vars),
6964 while (lex_match (lexer, T_ENDCMD))
6966 if (lex_token (lexer) == T_STOP)
6968 msg (SE, _("Unexpected end of input expecting matrix command."));
6972 if (lex_match_phrase (lexer, "END MATRIX"))
6975 struct matrix_cmd *c = matrix_parse_command (&state);
6978 matrix_cmd_execute (c);
6979 matrix_cmd_destroy (c);
6983 struct matrix_var *var, *next;
6984 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6987 gsl_matrix_free (var->value);
6988 hmap_delete (&state.vars, &var->hmap_node);
6991 hmap_destroy (&state.vars);
6994 dict_unref (state.common->dict);
6995 casewriter_destroy (state.common->writer);
6996 free (state.common);
6998 fh_unref (state.prev_read_file);
6999 for (size_t i = 0; i < state.n_read_files; i++)
7000 read_file_destroy (state.read_files[i]);
7001 free (state.read_files);
7002 fh_unref (state.prev_write_file);
7003 for (size_t i = 0; i < state.n_write_files; i++)
7004 write_file_destroy (state.write_files[i]);
7005 free (state.write_files);
7006 fh_unref (state.prev_save_file);
7007 for (size_t i = 0; i < state.n_save_files; i++)
7008 save_file_destroy (state.save_files[i]);
7009 free (state.save_files);