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/minmax.h"
66 #define _(msgid) gettext (msgid)
67 #define N_(msgid) (msgid)
71 struct hmap_node hmap_node;
79 struct file_handle *outfile;
80 struct string_array variables;
81 struct string_array fnames;
82 struct string_array snames;
87 /* Execution state. */
88 struct dictionary *dict;
89 struct casewriter *writer;
94 struct file_handle *file;
95 struct dfm_reader *reader;
101 struct dataset *dataset;
102 struct session *session;
106 struct file_handle *prev_save_outfile;
107 struct file_handle *prev_write_outfile;
108 struct msave_common *common;
110 struct file_handle *prev_read_file;
111 struct read_file **read_files;
115 static struct matrix_var *
116 matrix_var_lookup (struct matrix_state *s, struct substring name)
118 struct matrix_var *var;
120 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
121 utf8_hash_case_substring (name, 0), &s->vars)
122 if (!utf8_sscasecmp (ss_cstr (var->name), name))
128 static struct matrix_var *
129 matrix_var_create (struct matrix_state *s, struct substring name)
131 struct matrix_var *var = xmalloc (sizeof *var);
132 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
133 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
138 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
140 gsl_matrix_free (var->value);
144 #define MATRIX_FUNCTIONS \
200 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
201 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
202 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
203 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
204 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
205 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
206 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
207 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
208 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
209 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
210 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
212 typedef gsl_matrix *matrix_proto_m_d (double);
213 typedef gsl_matrix *matrix_proto_m_dd (double, double);
214 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
215 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
216 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
217 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
218 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
219 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
220 typedef double matrix_proto_d_m (gsl_matrix *);
221 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
222 typedef gsl_matrix *matrix_proto_IDENT (double, double);
224 #define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
228 /* Matrix expressions. */
235 #define F(NAME, PROTOTYPE) MOP_F_##NAME,
239 /* Elementwise and scalar arithmetic. */
240 MOP_NEGATE, /* unary - */
241 MOP_ADD_ELEMS, /* + */
242 MOP_SUB_ELEMS, /* - */
243 MOP_MUL_ELEMS, /* &* */
244 MOP_DIV_ELEMS, /* / and &/ */
245 MOP_EXP_ELEMS, /* &** */
247 MOP_SEQ_BY, /* a:b:c */
249 /* Matrix arithmetic. */
251 MOP_EXP_MAT, /* ** */
268 MOP_PASTE_HORZ, /* a, b, c, ... */
269 MOP_PASTE_VERT, /* a; b; c; ... */
273 MOP_VEC_INDEX, /* x(y) */
274 MOP_VEC_ALL, /* x(:) */
275 MOP_MAT_INDEX, /* x(y,z) */
276 MOP_ROW_INDEX, /* x(y,:) */
277 MOP_COL_INDEX, /* x(:,z) */
284 MOP_EOF, /* EOF('file') */
292 struct matrix_expr **subs;
297 struct matrix_var *variable;
298 struct read_file *eof;
303 matrix_expr_destroy (struct matrix_expr *e)
310 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
341 for (size_t i = 0; i < e->n_subs; i++)
342 matrix_expr_destroy (e->subs[i]);
354 static struct matrix_expr *
355 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
358 struct matrix_expr *e = xmalloc (sizeof *e);
359 *e = (struct matrix_expr) {
361 .subs = xmemdup (subs, n_subs * sizeof *subs),
367 static struct matrix_expr *
368 matrix_expr_create_0 (enum matrix_op op)
370 struct matrix_expr *sub;
371 return matrix_expr_create_subs (op, &sub, 0);
374 static struct matrix_expr *
375 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
377 return matrix_expr_create_subs (op, &sub, 1);
380 static struct matrix_expr *
381 matrix_expr_create_2 (enum matrix_op op,
382 struct matrix_expr *sub0, struct matrix_expr *sub1)
384 struct matrix_expr *subs[] = { sub0, sub1 };
385 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
388 static struct matrix_expr *
389 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
390 struct matrix_expr *sub1, struct matrix_expr *sub2)
392 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
393 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
396 static struct matrix_expr *
397 matrix_expr_create_number (double number)
399 struct matrix_expr *e = xmalloc (sizeof *e);
400 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
404 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
406 static struct matrix_expr *
407 matrix_parse_curly_comma (struct matrix_state *s)
409 struct matrix_expr *lhs = matrix_parse_or_xor (s);
413 while (lex_match (s->lexer, T_COMMA))
415 struct matrix_expr *rhs = matrix_parse_or_xor (s);
418 matrix_expr_destroy (lhs);
421 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
426 static struct matrix_expr *
427 matrix_parse_curly_semi (struct matrix_state *s)
429 if (lex_token (s->lexer) == T_RCURLY)
430 return matrix_expr_create_0 (MOP_EMPTY);
432 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
436 while (lex_match (s->lexer, T_SEMICOLON))
438 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
441 matrix_expr_destroy (lhs);
444 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
449 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
450 for (size_t Y = 0; Y < (M)->size1; Y++) \
451 for (size_t X = 0; X < (M)->size2; X++) \
452 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
455 is_vector (const gsl_matrix *m)
457 return m->size1 <= 1 || m->size2 <= 1;
461 to_vector (gsl_matrix *m)
463 return (m->size1 == 1
464 ? gsl_matrix_row (m, 0).vector
465 : gsl_matrix_column (m, 0).vector);
470 matrix_eval_ABS (gsl_matrix *m)
472 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
478 matrix_eval_ALL (gsl_matrix *m)
480 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
487 matrix_eval_ANY (gsl_matrix *m)
489 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
496 matrix_eval_ARSIN (gsl_matrix *m)
498 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
504 matrix_eval_ARTAN (gsl_matrix *m)
506 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
512 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
516 for (size_t i = 0; i < n; i++)
521 gsl_matrix *block = gsl_matrix_calloc (r, c);
523 for (size_t i = 0; i < n; i++)
525 for (size_t y = 0; y < m[i]->size1; y++)
526 for (size_t x = 0; x < m[i]->size2; x++)
527 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
535 matrix_eval_CHOL (gsl_matrix *m)
537 if (!gsl_linalg_cholesky_decomp1 (m))
539 for (size_t y = 0; y < m->size1; y++)
540 for (size_t x = y + 1; x < m->size2; x++)
541 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
543 for (size_t y = 0; y < m->size1; y++)
544 for (size_t x = 0; x < y; x++)
545 gsl_matrix_set (m, y, x, 0);
550 msg (SE, _("Input to CHOL function is not positive-definite."));
556 matrix_eval_col_extremum (gsl_matrix *m, bool min)
561 return gsl_matrix_alloc (1, 0);
563 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
564 for (size_t x = 0; x < m->size2; x++)
566 double ext = gsl_matrix_get (m, 0, x);
567 for (size_t y = 1; y < m->size1; y++)
569 double value = gsl_matrix_get (m, y, x);
570 if (min ? value < ext : value > ext)
573 gsl_matrix_set (cext, 0, x, ext);
579 matrix_eval_CMAX (gsl_matrix *m)
581 return matrix_eval_col_extremum (m, false);
585 matrix_eval_CMIN (gsl_matrix *m)
587 return matrix_eval_col_extremum (m, true);
591 matrix_eval_COS (gsl_matrix *m)
593 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
599 matrix_eval_col_sum (gsl_matrix *m, bool square)
604 return gsl_matrix_alloc (1, 0);
606 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
607 for (size_t x = 0; x < m->size2; x++)
610 for (size_t y = 0; y < m->size1; y++)
612 double d = gsl_matrix_get (m, y, x);
613 sum += square ? pow2 (d) : d;
615 gsl_matrix_set (result, 0, x, sum);
621 matrix_eval_CSSQ (gsl_matrix *m)
623 return matrix_eval_col_sum (m, true);
627 matrix_eval_CSUM (gsl_matrix *m)
629 return matrix_eval_col_sum (m, false);
633 compare_double_3way (const void *a_, const void *b_)
635 const double *a = a_;
636 const double *b = b_;
637 return *a < *b ? -1 : *a > *b;
641 matrix_eval_DESIGN (gsl_matrix *m)
643 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
644 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
645 gsl_matrix_transpose_memcpy (&m2, m);
647 for (size_t y = 0; y < m2.size1; y++)
648 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
650 size_t *n = xcalloc (m2.size1, sizeof *n);
652 for (size_t i = 0; i < m2.size1; i++)
654 double *row = tmp + m2.size2 * i;
655 for (size_t j = 0; j < m2.size2; )
658 for (k = j + 1; k < m2.size2; k++)
659 if (row[j] != row[k])
661 row[n[i]++] = row[j];
666 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
671 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
673 for (size_t i = 0; i < m->size2; i++)
678 const double *unique = tmp + m2.size2 * i;
679 for (size_t j = 0; j < n[i]; j++, x++)
681 double value = unique[j];
682 for (size_t y = 0; y < m->size1; y++)
683 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
694 matrix_eval_DET (gsl_matrix *m)
696 gsl_permutation *p = gsl_permutation_alloc (m->size1);
698 gsl_linalg_LU_decomp (m, p, &signum);
699 gsl_permutation_free (p);
700 return gsl_linalg_LU_det (m, signum);
704 matrix_eval_DIAG (gsl_matrix *m)
706 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
707 for (size_t i = 0; i < diag->size1; i++)
708 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
713 is_symmetric (const gsl_matrix *m)
715 if (m->size1 != m->size2)
718 for (size_t y = 0; y < m->size1; y++)
719 for (size_t x = 0; x < y; x++)
720 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
727 compare_double_desc (const void *a_, const void *b_)
729 const double *a = a_;
730 const double *b = b_;
731 return *a > *b ? -1 : *a < *b;
735 matrix_eval_EVAL (gsl_matrix *m)
737 if (!is_symmetric (m))
739 msg (SE, _("Argument of EVAL must be symmetric."));
743 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
744 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
745 gsl_vector v_eval = to_vector (eval);
746 gsl_eigen_symm (m, &v_eval, w);
747 gsl_eigen_symm_free (w);
749 assert (v_eval.stride == 1);
750 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
756 matrix_eval_EXP (gsl_matrix *m)
758 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
763 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
766 Charl Linssen <charl@itfromb.it>
770 matrix_eval_GINV (gsl_matrix *A)
775 gsl_matrix *tmp_mat = NULL;
778 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
779 tmp_mat = gsl_matrix_alloc (m, n);
780 gsl_matrix_transpose_memcpy (tmp_mat, A);
788 gsl_matrix *V = gsl_matrix_alloc (m, m);
789 gsl_vector *u = gsl_vector_alloc (m);
791 gsl_vector *tmp_vec = gsl_vector_alloc (m);
792 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
793 gsl_vector_free (tmp_vec);
796 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
797 gsl_matrix_set_zero (Sigma_pinv);
798 double cutoff = 1e-15 * gsl_vector_max (u);
800 for (size_t i = 0; i < m; ++i)
802 double x = gsl_vector_get (u, i);
803 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
806 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
807 gsl_matrix *U = gsl_matrix_calloc (n, n);
808 for (size_t i = 0; i < n; i++)
809 for (size_t j = 0; j < m; j++)
810 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
812 /* two dot products to obtain pseudoinverse */
813 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
814 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
819 A_pinv = gsl_matrix_alloc (n, m);
820 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
824 A_pinv = gsl_matrix_alloc (m, n);
825 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
828 gsl_matrix_free (tmp_mat);
829 gsl_matrix_free (tmp_mat2);
831 gsl_matrix_free (Sigma_pinv);
845 grade_compare_3way (const void *a_, const void *b_)
847 const struct grade *a = a_;
848 const struct grade *b = b_;
850 return (a->value < b->value ? -1
851 : a->value > b->value ? 1
859 matrix_eval_GRADE (gsl_matrix *m)
861 size_t n = m->size1 * m->size2;
862 struct grade *grades = xmalloc (n * sizeof *grades);
865 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
866 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
867 qsort (grades, n, sizeof *grades, grade_compare_3way);
869 for (size_t i = 0; i < n; i++)
870 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
878 dot (gsl_vector *a, gsl_vector *b)
881 for (size_t i = 0; i < a->size; i++)
882 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
887 norm2 (gsl_vector *v)
890 for (size_t i = 0; i < v->size; i++)
891 result += pow2 (gsl_vector_get (v, i));
898 return sqrt (norm2 (v));
902 matrix_eval_GSCH (gsl_matrix *v)
904 if (v->size2 < v->size1)
906 msg (SE, _("GSCH requires its argument to have at least as many columns "
907 "as rows, but it has dimensions %zu×%zu."),
911 if (!v->size1 || !v->size2)
914 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
916 for (size_t vx = 0; vx < v->size2; vx++)
918 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
919 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
921 gsl_vector_memcpy (&u_i, &v_i);
922 for (size_t j = 0; j < ux; j++)
924 gsl_vector u_j = gsl_matrix_column (u, j).vector;
925 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
926 for (size_t k = 0; k < u_i.size; k++)
927 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
928 - scale * gsl_vector_get (&u_j, k)));
931 double len = norm (&u_i);
934 gsl_vector_scale (&u_i, 1.0 / len);
935 if (++ux >= v->size1)
942 msg (SE, _("%zu×%zu argument to GSCH contains only "
943 "%zu linearly independent columns."),
944 v->size1, v->size2, ux);
954 matrix_eval_IDENT (double s1, double s2)
956 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
958 msg (SE, _("Arguments to IDENT must be non-negative integers."));
962 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
963 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
969 invert_matrix (gsl_matrix *x)
971 gsl_permutation *p = gsl_permutation_alloc (x->size1);
973 gsl_linalg_LU_decomp (x, p, &signum);
974 gsl_linalg_LU_invx (x, p);
975 gsl_permutation_free (p);
979 matrix_eval_INV (gsl_matrix *m)
986 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
988 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
989 a->size2 * b->size2);
991 for (size_t ar = 0; ar < a->size1; ar++)
992 for (size_t br = 0; br < b->size1; br++, y++)
995 for (size_t ac = 0; ac < a->size2; ac++)
996 for (size_t bc = 0; bc < b->size2; bc++, x++)
998 double av = gsl_matrix_get (a, ar, ac);
999 double bv = gsl_matrix_get (b, br, bc);
1000 gsl_matrix_set (k, y, x, av * bv);
1007 matrix_eval_LG10 (gsl_matrix *m)
1009 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1015 matrix_eval_LN (gsl_matrix *m)
1017 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1023 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1025 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1028 for (size_t i = 1; i <= n * n; i++)
1030 gsl_matrix_set (m, y, x, i);
1032 size_t y1 = !y ? n - 1 : y - 1;
1033 size_t x1 = x + 1 >= n ? 0 : x + 1;
1034 if (gsl_matrix_get (m, y1, x1) == 0)
1040 y = y + 1 >= n ? 0 : y + 1;
1045 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1047 double a = gsl_matrix_get (m, y1, x1);
1048 double b = gsl_matrix_get (m, y2, x2);
1049 gsl_matrix_set (m, y1, x1, b);
1050 gsl_matrix_set (m, y2, x2, a);
1054 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1058 /* A. Umar, "On the Construction of Even Order Magic Squares",
1059 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1061 for (size_t i = 1; i <= n * n / 2; i++)
1063 gsl_matrix_set (m, y, x, i);
1073 for (size_t i = n * n; i > n * n / 2; i--)
1075 gsl_matrix_set (m, y, x, i);
1083 for (size_t y = 0; y < n; y++)
1084 for (size_t x = 0; x < n / 2; x++)
1086 unsigned int d = gsl_matrix_get (m, y, x);
1087 if (d % 2 != (y < n / 2))
1088 magic_exchange (m, y, x, y, n - x - 1);
1093 size_t x1 = n / 2 - 1;
1095 magic_exchange (m, y1, x1, y2, x1);
1096 magic_exchange (m, y1, x2, y2, x2);
1100 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1102 /* A. Umar, "On the Construction of Even Order Magic Squares",
1103 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1107 for (size_t i = 1; ; i++)
1109 gsl_matrix_set (m, y, x, i);
1110 if (++y == n / 2 - 1)
1122 for (size_t i = n * n; ; i--)
1124 gsl_matrix_set (m, y, x, i);
1125 if (++y == n / 2 - 1)
1134 for (size_t y = 0; y < n; y++)
1135 if (y != n / 2 - 1 && y != n / 2)
1136 for (size_t x = 0; x < n / 2; x++)
1138 unsigned int d = gsl_matrix_get (m, y, x);
1139 if (d % 2 != (y < n / 2))
1140 magic_exchange (m, y, x, y, n - x - 1);
1143 size_t a0 = (n * n - 2 * n) / 2 + 1;
1144 for (size_t i = 0; i < n / 2; i++)
1147 gsl_matrix_set (m, n / 2, i, a);
1148 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1150 for (size_t i = 0; i < n / 2; i++)
1152 size_t a = a0 + i + n / 2;
1153 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1154 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1156 for (size_t x = 1; x < n / 2; x += 2)
1157 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1158 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1159 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1160 size_t x1 = n / 2 - 2;
1161 size_t x2 = n / 2 + 1;
1162 size_t y1 = n / 2 - 2;
1163 size_t y2 = n / 2 + 1;
1164 magic_exchange (m, y1, x1, y2, x1);
1165 magic_exchange (m, y1, x2, y2, x2);
1169 matrix_eval_MAGIC (double n_)
1171 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1173 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1178 gsl_matrix *m = gsl_matrix_calloc (n, n);
1180 matrix_eval_MAGIC_odd (m, n);
1182 matrix_eval_MAGIC_singly_even (m, n);
1184 matrix_eval_MAGIC_doubly_even (m, n);
1189 matrix_eval_MAKE (double r, double c, double value)
1191 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1193 msg (SE, _("Size arguments to MAKE must be integers."));
1197 gsl_matrix *m = gsl_matrix_alloc (r, c);
1198 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1204 matrix_eval_MDIAG (gsl_vector *v)
1206 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1207 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1208 gsl_vector_memcpy (&diagonal, v);
1213 matrix_eval_MMAX (gsl_matrix *m)
1215 return gsl_matrix_max (m);
1219 matrix_eval_MMIN (gsl_matrix *m)
1221 return gsl_matrix_min (m);
1225 matrix_eval_MOD (gsl_matrix *m, double divisor)
1229 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1233 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1234 *d = fmod (*d, divisor);
1239 matrix_eval_MSSQ (gsl_matrix *m)
1242 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1248 matrix_eval_MSUM (gsl_matrix *m)
1251 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1257 matrix_eval_NCOL (gsl_matrix *m)
1263 matrix_eval_NROW (gsl_matrix *m)
1269 matrix_eval_RANK (gsl_matrix *m)
1271 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1272 gsl_linalg_QR_decomp (m, tau);
1273 gsl_vector_free (tau);
1275 return gsl_linalg_QRPT_rank (m, -1);
1279 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1281 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1283 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1288 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1290 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1291 "product of matrix dimensions."));
1295 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1298 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1300 gsl_matrix_set (dst, y1, x1, *d);
1311 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1316 return gsl_matrix_alloc (0, 1);
1318 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1319 for (size_t y = 0; y < m->size1; y++)
1321 double ext = gsl_matrix_get (m, y, 0);
1322 for (size_t x = 1; x < m->size2; x++)
1324 double value = gsl_matrix_get (m, y, x);
1325 if (min ? value < ext : value > ext)
1328 gsl_matrix_set (rext, y, 0, ext);
1334 matrix_eval_RMAX (gsl_matrix *m)
1336 return matrix_eval_row_extremum (m, false);
1340 matrix_eval_RMIN (gsl_matrix *m)
1342 return matrix_eval_row_extremum (m, true);
1346 matrix_eval_RND (gsl_matrix *m)
1348 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1360 rank_compare_3way (const void *a_, const void *b_)
1362 const struct rank *a = a_;
1363 const struct rank *b = b_;
1365 return a->value < b->value ? -1 : a->value > b->value;
1369 matrix_eval_RNKORDER (gsl_matrix *m)
1371 size_t n = m->size1 * m->size2;
1372 struct rank *ranks = xmalloc (n * sizeof *ranks);
1374 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1375 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1376 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1378 for (size_t i = 0; i < n; )
1381 for (j = i + 1; j < n; j++)
1382 if (ranks[i].value != ranks[j].value)
1385 double rank = (i + j + 1.0) / 2.0;
1386 for (size_t k = i; k < j; k++)
1387 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1398 matrix_eval_row_sum (gsl_matrix *m, bool square)
1403 return gsl_matrix_alloc (0, 1);
1405 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1406 for (size_t y = 0; y < m->size1; y++)
1409 for (size_t x = 0; x < m->size2; x++)
1411 double d = gsl_matrix_get (m, y, x);
1412 sum += square ? pow2 (d) : d;
1414 gsl_matrix_set (result, y, 0, sum);
1420 matrix_eval_RSSQ (gsl_matrix *m)
1422 return matrix_eval_row_sum (m, true);
1426 matrix_eval_RSUM (gsl_matrix *m)
1428 return matrix_eval_row_sum (m, false);
1432 matrix_eval_SIN (gsl_matrix *m)
1434 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1440 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1442 if (m1->size1 != m2->size1)
1444 msg (SE, _("SOLVE requires its arguments to have the same number of "
1445 "rows, but the first argument has dimensions %zu×%zu and "
1446 "the second %zu×%zu."),
1447 m1->size1, m1->size2,
1448 m2->size1, m2->size2);
1452 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1453 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1455 gsl_linalg_LU_decomp (m1, p, &signum);
1456 for (size_t i = 0; i < m2->size2; i++)
1458 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1459 gsl_vector xi = gsl_matrix_column (x, i).vector;
1460 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1462 gsl_permutation_free (p);
1467 matrix_eval_SQRT (gsl_matrix *m)
1469 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1473 msg (SE, _("Argument to SQRT must be nonnegative."));
1482 matrix_eval_SSCP (gsl_matrix *m)
1484 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1485 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1490 matrix_eval_SVAL (gsl_matrix *m)
1492 gsl_matrix *tmp_mat = NULL;
1493 if (m->size2 > m->size1)
1495 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1496 gsl_matrix_transpose_memcpy (tmp_mat, m);
1501 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1502 gsl_vector *S = gsl_vector_alloc (m->size2);
1503 gsl_vector *work = gsl_vector_alloc (m->size2);
1504 gsl_linalg_SV_decomp (m, V, S, work);
1506 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1507 for (size_t i = 0; i < m->size2; i++)
1508 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1510 gsl_matrix_free (V);
1511 gsl_vector_free (S);
1512 gsl_vector_free (work);
1513 gsl_matrix_free (tmp_mat);
1519 matrix_eval_SWEEP (gsl_matrix *m, double d)
1521 if (d < 1 || d > SIZE_MAX)
1523 msg (SE, _("Scalar argument to SWEEP must be integer."));
1527 if (k >= MIN (m->size1, m->size2))
1529 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1530 "equal to the smaller of the matrix argument's rows and "
1535 double m_kk = gsl_matrix_get (m, k, k);
1536 if (fabs (m_kk) > 1e-19)
1538 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1539 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1541 double m_ij = gsl_matrix_get (m, i, j);
1542 double m_ik = gsl_matrix_get (m, i, k);
1543 double m_kj = gsl_matrix_get (m, k, j);
1544 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1553 for (size_t i = 0; i < m->size1; i++)
1555 gsl_matrix_set (m, i, k, 0);
1556 gsl_matrix_set (m, k, i, 0);
1563 matrix_eval_TRACE (gsl_matrix *m)
1566 size_t n = MIN (m->size1, m->size2);
1567 for (size_t i = 0; i < n; i++)
1568 sum += gsl_matrix_get (m, i, i);
1573 matrix_eval_T (gsl_matrix *m)
1575 return matrix_eval_TRANSPOS (m);
1579 matrix_eval_TRANSPOS (gsl_matrix *m)
1581 if (m->size1 == m->size2)
1583 gsl_matrix_transpose (m);
1588 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1589 gsl_matrix_transpose_memcpy (t, m);
1595 matrix_eval_TRUNC (gsl_matrix *m)
1597 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1603 matrix_eval_UNIFORM (double r_, double c_)
1605 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1607 msg (SE, _("Arguments to UNIFORM must be integers."));
1612 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1614 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1618 gsl_matrix *m = gsl_matrix_alloc (r, c);
1619 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1620 *d = gsl_ran_flat (get_rng (), 0, 1);
1624 struct matrix_function
1628 size_t min_args, max_args;
1631 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1633 static const struct matrix_function *
1634 matrix_parse_function_name (const char *token)
1636 static const struct matrix_function functions[] =
1638 #define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1642 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1644 for (size_t i = 0; i < N_FUNCTIONS; i++)
1646 if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1647 return &functions[i];
1652 static struct read_file *
1653 read_file_create (struct matrix_state *s, struct file_handle *fh)
1655 for (size_t i = 0; i < s->n_read_files; i++)
1657 struct read_file *rf = s->read_files[i];
1665 struct read_file *rf = xmalloc (sizeof *rf);
1666 *rf = (struct read_file) { .file = fh };
1668 s->read_files = xrealloc (s->read_files,
1669 (s->n_read_files + 1) * sizeof *s->read_files);
1670 s->read_files[s->n_read_files++] = rf;
1674 static struct dfm_reader *
1675 read_file_open (struct read_file *rf)
1678 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1683 read_file_destroy (struct read_file *rf)
1687 fh_unref (rf->file);
1688 dfm_close_reader (rf->reader);
1689 free (rf->encoding);
1695 matrix_parse_function (struct matrix_state *s, const char *token,
1696 struct matrix_expr **exprp)
1699 if (lex_next_token (s->lexer, 1) != T_LPAREN)
1702 if (lex_match_id (s->lexer, "EOF"))
1705 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1709 if (!lex_force_match (s->lexer, T_RPAREN))
1715 struct read_file *rf = read_file_create (s, fh);
1717 struct matrix_expr *e = xmalloc (sizeof *e);
1718 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1723 const struct matrix_function *f = matrix_parse_function_name (token);
1727 lex_get_n (s->lexer, 2);
1729 struct matrix_expr *e = xmalloc (sizeof *e);
1730 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1732 size_t allocated_subs = 0;
1735 struct matrix_expr *sub = matrix_parse_expr (s);
1739 if (e->n_subs >= allocated_subs)
1740 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1741 e->subs[e->n_subs++] = sub;
1743 while (lex_match (s->lexer, T_COMMA));
1744 if (!lex_force_match (s->lexer, T_RPAREN))
1751 matrix_expr_destroy (e);
1755 static struct matrix_expr *
1756 matrix_parse_primary (struct matrix_state *s)
1758 if (lex_is_number (s->lexer))
1760 double number = lex_number (s->lexer);
1762 return matrix_expr_create_number (number);
1764 else if (lex_is_string (s->lexer))
1766 char string[sizeof (double)];
1767 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1771 memcpy (&number, string, sizeof number);
1772 return matrix_expr_create_number (number);
1774 else if (lex_match (s->lexer, T_LPAREN))
1776 struct matrix_expr *e = matrix_parse_or_xor (s);
1777 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1779 matrix_expr_destroy (e);
1784 else if (lex_match (s->lexer, T_LCURLY))
1786 struct matrix_expr *e = matrix_parse_curly_semi (s);
1787 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1789 matrix_expr_destroy (e);
1794 else if (lex_token (s->lexer) == T_ID)
1796 struct matrix_expr *retval;
1797 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1800 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1803 lex_error (s->lexer, _("Unknown variable %s."),
1804 lex_tokcstr (s->lexer));
1809 struct matrix_expr *e = xmalloc (sizeof *e);
1810 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1813 else if (lex_token (s->lexer) == T_ALL)
1815 struct matrix_expr *retval;
1816 if (matrix_parse_function (s, "ALL", &retval))
1820 lex_error (s->lexer, NULL);
1824 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1827 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1829 if (lex_match (s->lexer, T_COLON))
1836 *indexp = matrix_parse_expr (s);
1837 return *indexp != NULL;
1841 static struct matrix_expr *
1842 matrix_parse_postfix (struct matrix_state *s)
1844 struct matrix_expr *lhs = matrix_parse_primary (s);
1845 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1848 struct matrix_expr *i0;
1849 if (!matrix_parse_index_expr (s, &i0))
1851 matrix_expr_destroy (lhs);
1854 if (lex_match (s->lexer, T_RPAREN))
1856 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1857 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1858 else if (lex_match (s->lexer, T_COMMA))
1860 struct matrix_expr *i1;
1861 if (!matrix_parse_index_expr (s, &i1)
1862 || !lex_force_match (s->lexer, T_RPAREN))
1864 matrix_expr_destroy (lhs);
1865 matrix_expr_destroy (i0);
1866 matrix_expr_destroy (i1);
1869 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1870 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1871 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1876 lex_error_expecting (s->lexer, "`)'", "`,'");
1881 static struct matrix_expr *
1882 matrix_parse_unary (struct matrix_state *s)
1884 if (lex_match (s->lexer, T_DASH))
1886 struct matrix_expr *lhs = matrix_parse_unary (s);
1887 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1889 else if (lex_match (s->lexer, T_PLUS))
1890 return matrix_parse_unary (s);
1892 return matrix_parse_postfix (s);
1895 static struct matrix_expr *
1896 matrix_parse_seq (struct matrix_state *s)
1898 struct matrix_expr *start = matrix_parse_unary (s);
1899 if (!start || !lex_match (s->lexer, T_COLON))
1902 struct matrix_expr *end = matrix_parse_unary (s);
1905 matrix_expr_destroy (start);
1909 if (lex_match (s->lexer, T_COLON))
1911 struct matrix_expr *increment = matrix_parse_unary (s);
1914 matrix_expr_destroy (start);
1915 matrix_expr_destroy (end);
1918 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
1921 return matrix_expr_create_2 (MOP_SEQ, start, end);
1924 static struct matrix_expr *
1925 matrix_parse_exp (struct matrix_state *s)
1927 struct matrix_expr *lhs = matrix_parse_seq (s);
1934 if (lex_match (s->lexer, T_EXP))
1936 else if (lex_match_phrase (s->lexer, "&**"))
1941 struct matrix_expr *rhs = matrix_parse_seq (s);
1944 matrix_expr_destroy (lhs);
1947 lhs = matrix_expr_create_2 (op, lhs, rhs);
1951 static struct matrix_expr *
1952 matrix_parse_mul_div (struct matrix_state *s)
1954 struct matrix_expr *lhs = matrix_parse_exp (s);
1961 if (lex_match (s->lexer, T_ASTERISK))
1963 else if (lex_match (s->lexer, T_SLASH))
1965 else if (lex_match_phrase (s->lexer, "&*"))
1967 else if (lex_match_phrase (s->lexer, "&/"))
1972 struct matrix_expr *rhs = matrix_parse_exp (s);
1975 matrix_expr_destroy (lhs);
1978 lhs = matrix_expr_create_2 (op, lhs, rhs);
1982 static struct matrix_expr *
1983 matrix_parse_add_sub (struct matrix_state *s)
1985 struct matrix_expr *lhs = matrix_parse_mul_div (s);
1992 if (lex_match (s->lexer, T_PLUS))
1994 else if (lex_match (s->lexer, T_DASH))
1996 else if (lex_token (s->lexer) == T_NEG_NUM)
2001 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2004 matrix_expr_destroy (lhs);
2007 lhs = matrix_expr_create_2 (op, lhs, rhs);
2011 static struct matrix_expr *
2012 matrix_parse_relational (struct matrix_state *s)
2014 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2021 if (lex_match (s->lexer, T_GT))
2023 else if (lex_match (s->lexer, T_GE))
2025 else if (lex_match (s->lexer, T_LT))
2027 else if (lex_match (s->lexer, T_LE))
2029 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2031 else if (lex_match (s->lexer, T_NE))
2036 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2039 matrix_expr_destroy (lhs);
2042 lhs = matrix_expr_create_2 (op, lhs, rhs);
2046 static struct matrix_expr *
2047 matrix_parse_not (struct matrix_state *s)
2049 if (lex_match (s->lexer, T_NOT))
2051 struct matrix_expr *lhs = matrix_parse_not (s);
2052 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2055 return matrix_parse_relational (s);
2058 static struct matrix_expr *
2059 matrix_parse_and (struct matrix_state *s)
2061 struct matrix_expr *lhs = matrix_parse_not (s);
2065 while (lex_match (s->lexer, T_AND))
2067 struct matrix_expr *rhs = matrix_parse_not (s);
2070 matrix_expr_destroy (lhs);
2073 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2078 static struct matrix_expr *
2079 matrix_parse_or_xor (struct matrix_state *s)
2081 struct matrix_expr *lhs = matrix_parse_and (s);
2088 if (lex_match (s->lexer, T_OR))
2090 else if (lex_match_id (s->lexer, "XOR"))
2095 struct matrix_expr *rhs = matrix_parse_and (s);
2098 matrix_expr_destroy (lhs);
2101 lhs = matrix_expr_create_2 (op, lhs, rhs);
2105 static struct matrix_expr *
2106 matrix_parse_expr (struct matrix_state *s)
2108 return matrix_parse_or_xor (s);
2111 /* Expression evaluation. */
2114 matrix_op_eval (enum matrix_op op, double a, double b)
2118 case MOP_ADD_ELEMS: return a + b;
2119 case MOP_SUB_ELEMS: return a - b;
2120 case MOP_MUL_ELEMS: return a * b;
2121 case MOP_DIV_ELEMS: return a / b;
2122 case MOP_EXP_ELEMS: return pow (a, b);
2123 case MOP_GT: return a > b;
2124 case MOP_GE: return a >= b;
2125 case MOP_LT: return a < b;
2126 case MOP_LE: return a <= b;
2127 case MOP_EQ: return a == b;
2128 case MOP_NE: return a != b;
2129 case MOP_AND: return (a > 0) && (b > 0);
2130 case MOP_OR: return (a > 0) || (b > 0);
2131 case MOP_XOR: return (a > 0) != (b > 0);
2133 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2142 case MOP_PASTE_HORZ:
2143 case MOP_PASTE_VERT:
2159 matrix_op_name (enum matrix_op op)
2163 case MOP_ADD_ELEMS: return "+";
2164 case MOP_SUB_ELEMS: return "-";
2165 case MOP_MUL_ELEMS: return "&*";
2166 case MOP_DIV_ELEMS: return "&/";
2167 case MOP_EXP_ELEMS: return "&**";
2168 case MOP_GT: return ">";
2169 case MOP_GE: return ">=";
2170 case MOP_LT: return "<";
2171 case MOP_LE: return "<=";
2172 case MOP_EQ: return "=";
2173 case MOP_NE: return "<>";
2174 case MOP_AND: return "AND";
2175 case MOP_OR: return "OR";
2176 case MOP_XOR: return "XOR";
2178 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2187 case MOP_PASTE_HORZ:
2188 case MOP_PASTE_VERT:
2204 is_scalar (const gsl_matrix *m)
2206 return m->size1 == 1 && m->size2 == 1;
2210 to_scalar (const gsl_matrix *m)
2212 assert (is_scalar (m));
2213 return gsl_matrix_get (m, 0, 0);
2217 matrix_expr_evaluate_elementwise (enum matrix_op op,
2218 gsl_matrix *a, gsl_matrix *b)
2222 double be = to_scalar (b);
2223 for (size_t r = 0; r < a->size1; r++)
2224 for (size_t c = 0; c < a->size2; c++)
2226 double *ae = gsl_matrix_ptr (a, r, c);
2227 *ae = matrix_op_eval (op, *ae, be);
2231 else if (is_scalar (a))
2233 double ae = to_scalar (a);
2234 for (size_t r = 0; r < b->size1; r++)
2235 for (size_t c = 0; c < b->size2; c++)
2237 double *be = gsl_matrix_ptr (b, r, c);
2238 *be = matrix_op_eval (op, ae, *be);
2242 else if (a->size1 == b->size1 && a->size2 == b->size2)
2244 for (size_t r = 0; r < a->size1; r++)
2245 for (size_t c = 0; c < a->size2; c++)
2247 double *ae = gsl_matrix_ptr (a, r, c);
2248 double be = gsl_matrix_get (b, r, c);
2249 *ae = matrix_op_eval (op, *ae, be);
2255 msg (SE, _("Operands to %s must have the same dimensions or one "
2256 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
2257 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2263 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2265 if (is_scalar (a) || is_scalar (b))
2266 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2268 if (a->size2 != b->size1)
2270 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
2271 "not conformable for multiplication."),
2272 a->size1, a->size2, b->size1, b->size2);
2276 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2277 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2282 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2284 gsl_matrix *tmp = *a;
2290 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2293 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2294 swap_matrix (z, tmp);
2298 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2300 mul_matrix (x, *x, *x, tmp);
2304 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2307 if (x->size1 != x->size2)
2309 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2310 "the left-hand size, not one with dimensions %zu×%zu."),
2311 x->size1, x->size2);
2316 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2317 "right-hand side, not a matrix with dimensions %zu×%zu."),
2318 b->size1, b->size2);
2321 double bf = to_scalar (b);
2322 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2324 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2325 "or outside the valid range."), bf);
2330 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2332 gsl_matrix_set_identity (y);
2336 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2338 for (unsigned long int n = labs (bl); n > 1; n /= 2)
2341 mul_matrix (&y, x, y, &t);
2342 square_matrix (&x, &t);
2345 square_matrix (&x, &t);
2347 mul_matrix (&y, x, y, &t);
2351 /* Garbage collection.
2353 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2354 a permutation of them. We are returning one of them; that one must not be
2355 destroyed. We must not destroy 'x_' because the caller owns it. */
2357 gsl_matrix_free (y_);
2359 gsl_matrix_free (t_);
2365 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2368 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2370 msg (SE, _("All operands of : operator must be scalars."));
2374 long int start = to_scalar (start_);
2375 long int end = to_scalar (end_);
2376 long int by = by_ ? to_scalar (by_) : 1;
2380 msg (SE, _("The increment operand to : must be nonzero."));
2384 long int n = (end >= start && by > 0 ? (end - start + by) / by
2385 : end < start && by < 0 ? (start - end - by) / -by
2387 gsl_matrix *m = gsl_matrix_alloc (1, n);
2388 for (long int i = 0; i < n; i++)
2389 gsl_matrix_set (m, 0, i, start + i * by);
2394 matrix_expr_evaluate_not (gsl_matrix *a)
2396 for (size_t r = 0; r < a->size1; r++)
2397 for (size_t c = 0; c < a->size2; c++)
2399 double *ae = gsl_matrix_ptr (a, r, c);
2406 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2408 if (a->size1 != b->size1)
2410 if (!a->size1 || !a->size2)
2412 else if (!b->size1 || !b->size2)
2415 msg (SE, _("All columns in a matrix must have the same number of rows, "
2416 "but this tries to paste matrices with %zu and %zu rows."),
2417 a->size1, b->size1);
2421 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2422 for (size_t y = 0; y < a->size1; y++)
2424 for (size_t x = 0; x < a->size2; x++)
2425 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2426 for (size_t x = 0; x < b->size2; x++)
2427 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2433 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2435 if (a->size2 != b->size2)
2437 if (!a->size1 || !a->size2)
2439 else if (!b->size1 || !b->size2)
2442 msg (SE, _("All rows in a matrix must have the same number of columns, "
2443 "but this tries to stack matrices with %zu and %zu columns."),
2444 a->size2, b->size2);
2448 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2449 for (size_t x = 0; x < a->size2; x++)
2451 for (size_t y = 0; y < a->size1; y++)
2452 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2453 for (size_t y = 0; y < b->size1; y++)
2454 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2460 matrix_to_vector (gsl_matrix *m)
2463 gsl_vector v = to_vector (m);
2464 assert (v.block == m->block || !v.block);
2468 return xmemdup (&v, sizeof v);
2482 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2485 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2486 enum index_type index_type, size_t other_size,
2487 struct index_vector *iv)
2496 msg (SE, _("Vector index must be scalar or vector, not a "
2498 m->size1, m->size2);
2502 msg (SE, _("Matrix row index must be scalar or vector, not a "
2504 m->size1, m->size2);
2508 msg (SE, _("Matrix column index must be scalar or vector, not a "
2510 m->size1, m->size2);
2516 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2517 *iv = (struct index_vector) {
2518 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2521 for (size_t i = 0; i < v.size; i++)
2523 double index = gsl_vector_get (&v, i);
2524 if (index < 1 || index >= size + 1)
2529 msg (SE, _("Index %g is out of range for vector "
2530 "with %zu elements."), index, size);
2534 msg (SE, _("%g is not a valid row index for "
2535 "a %zu×%zu matrix."),
2536 index, size, other_size);
2540 msg (SE, _("%g is not a valid column index for "
2541 "a %zu×%zu matrix."),
2542 index, other_size, size);
2549 iv->indexes[i] = index - 1;
2555 *iv = (struct index_vector) {
2556 .indexes = xnmalloc (size, sizeof *iv->indexes),
2559 for (size_t i = 0; i < size; i++)
2566 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2568 if (!is_vector (sm))
2570 msg (SE, _("Vector index operator may not be applied to "
2571 "a %zu×%zu matrix."),
2572 sm->size1, sm->size2);
2580 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2582 if (!matrix_expr_evaluate_vec_all (sm))
2585 gsl_vector sv = to_vector (sm);
2586 struct index_vector iv;
2587 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2590 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2591 sm->size1 == 1 ? iv.n : 1);
2592 gsl_vector dv = to_vector (dm);
2593 for (size_t dx = 0; dx < iv.n; dx++)
2595 size_t sx = iv.indexes[dx];
2596 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2604 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2607 struct index_vector iv0;
2608 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2611 struct index_vector iv1;
2612 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2619 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2620 for (size_t dy = 0; dy < iv0.n; dy++)
2622 size_t sy = iv0.indexes[dy];
2624 for (size_t dx = 0; dx < iv1.n; dx++)
2626 size_t sx = iv1.indexes[dx];
2627 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2635 #define F(NAME, PROTOTYPE) \
2636 static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2637 const char *name, gsl_matrix *[], size_t, \
2638 matrix_proto_##PROTOTYPE *);
2643 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2645 if (!is_scalar (subs[index]))
2647 msg (SE, _("Function %s argument %zu must be a scalar, "
2648 "not a %zu×%zu matrix."),
2649 name, index + 1, subs[index]->size1, subs[index]->size2);
2656 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2658 if (!is_vector (subs[index]))
2660 msg (SE, _("Function %s argument %zu must be a vector, "
2661 "not a %zu×%zu matrix."),
2662 name, index + 1, subs[index]->size1, subs[index]->size2);
2669 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2671 for (size_t i = 0; i < n_subs; i++)
2673 if (!check_scalar_arg (name, subs, i))
2675 d[i] = to_scalar (subs[i]);
2681 matrix_expr_evaluate_m_d (const char *name,
2682 gsl_matrix *subs[], size_t n_subs,
2683 matrix_proto_m_d *f)
2685 assert (n_subs == 1);
2688 return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2692 matrix_expr_evaluate_m_dd (const char *name,
2693 gsl_matrix *subs[], size_t n_subs,
2694 matrix_proto_m_dd *f)
2696 assert (n_subs == 2);
2699 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2703 matrix_expr_evaluate_m_ddd (const char *name,
2704 gsl_matrix *subs[], size_t n_subs,
2705 matrix_proto_m_ddd *f)
2707 assert (n_subs == 3);
2710 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2714 matrix_expr_evaluate_m_m (const char *name UNUSED,
2715 gsl_matrix *subs[], size_t n_subs,
2716 matrix_proto_m_m *f)
2718 assert (n_subs == 1);
2723 matrix_expr_evaluate_m_md (const char *name UNUSED,
2724 gsl_matrix *subs[], size_t n_subs,
2725 matrix_proto_m_md *f)
2727 assert (n_subs == 2);
2728 if (!check_scalar_arg (name, subs, 1))
2730 return f (subs[0], to_scalar (subs[1]));
2734 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2735 gsl_matrix *subs[], size_t n_subs,
2736 matrix_proto_m_mdd *f)
2738 assert (n_subs == 3);
2739 if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2741 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2745 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2746 gsl_matrix *subs[], size_t n_subs,
2747 matrix_proto_m_mm *f)
2749 assert (n_subs == 2);
2750 return f (subs[0], subs[1]);
2754 matrix_expr_evaluate_m_v (const char *name,
2755 gsl_matrix *subs[], size_t n_subs,
2756 matrix_proto_m_v *f)
2758 assert (n_subs == 1);
2759 if (!check_vector_arg (name, subs, 0))
2761 gsl_vector v = to_vector (subs[0]);
2766 matrix_expr_evaluate_d_m (const char *name UNUSED,
2767 gsl_matrix *subs[], size_t n_subs,
2768 matrix_proto_d_m *f)
2770 assert (n_subs == 1);
2772 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2773 gsl_matrix_set (m, 0, 0, f (subs[0]));
2778 matrix_expr_evaluate_m_any (const char *name UNUSED,
2779 gsl_matrix *subs[], size_t n_subs,
2780 matrix_proto_m_any *f)
2782 return f (subs, n_subs);
2786 matrix_expr_evaluate_IDENT (const char *name,
2787 gsl_matrix *subs[], size_t n_subs,
2788 matrix_proto_IDENT *f)
2790 assert (n_subs <= 2);
2793 if (!to_scalar_args (name, subs, n_subs, d))
2796 return f (d[0], d[n_subs - 1]);
2800 matrix_expr_evaluate (const struct matrix_expr *e)
2802 if (e->op == MOP_NUMBER)
2804 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2805 gsl_matrix_set (m, 0, 0, e->number);
2808 else if (e->op == MOP_VARIABLE)
2810 const gsl_matrix *src = e->variable->value;
2813 msg (SE, _("Uninitialized variable %s used in expression."),
2818 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2819 gsl_matrix_memcpy (dst, src);
2822 else if (e->op == MOP_EOF)
2824 struct dfm_reader *reader = read_file_open (e->eof);
2825 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2826 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2830 enum { N_LOCAL = 3 };
2831 gsl_matrix *local_subs[N_LOCAL];
2832 gsl_matrix **subs = (e->n_subs < N_LOCAL
2834 : xmalloc (e->n_subs * sizeof *subs));
2836 for (size_t i = 0; i < e->n_subs; i++)
2838 subs[i] = matrix_expr_evaluate (e->subs[i]);
2841 for (size_t j = 0; j < i; j++)
2842 gsl_matrix_free (subs[j]);
2843 if (subs != local_subs)
2849 gsl_matrix *result = NULL;
2852 #define F(NAME, PROTOTYPE) \
2853 case MOP_F_##NAME: \
2854 result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \
2856 matrix_eval_##NAME); \
2862 gsl_matrix_scale (subs[0], -1.0);
2880 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2884 result = matrix_expr_evaluate_not (subs[0]);
2888 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2892 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2896 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2900 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2903 case MOP_PASTE_HORZ:
2904 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2907 case MOP_PASTE_VERT:
2908 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2912 result = gsl_matrix_alloc (0, 0);
2916 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
2920 result = matrix_expr_evaluate_vec_all (subs[0]);
2924 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
2928 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
2932 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
2941 for (size_t i = 0; i < e->n_subs; i++)
2942 if (subs[i] != result)
2943 gsl_matrix_free (subs[i]);
2944 if (subs != local_subs)
2950 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
2953 gsl_matrix *m = matrix_expr_evaluate (e);
2959 msg (SE, _("Expression for %s must evaluate to scalar, "
2960 "not a %zu×%zu matrix."),
2961 context, m->size1, m->size2);
2962 gsl_matrix_free (m);
2967 gsl_matrix_free (m);
2972 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
2976 if (!matrix_expr_evaluate_scalar (e, context, &d))
2980 if (d < LONG_MIN || d > LONG_MAX)
2982 msg (SE, _("Expression for %s is outside the integer range."), context);
2989 struct matrix_lvalue
2991 struct matrix_var *var;
2992 struct matrix_expr *indexes[2];
2997 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
3001 for (size_t i = 0; i < lvalue->n_indexes; i++)
3002 matrix_expr_destroy (lvalue->indexes[i]);
3007 static struct matrix_lvalue *
3008 matrix_lvalue_parse (struct matrix_state *s)
3010 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3011 if (!lex_force_id (s->lexer))
3013 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3014 if (lex_next_token (s->lexer, 1) == T_LPAREN)
3018 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3023 lex_get_n (s->lexer, 2);
3025 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3027 lvalue->n_indexes++;
3029 if (lex_match (s->lexer, T_COMMA))
3031 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3033 lvalue->n_indexes++;
3035 if (!lex_force_match (s->lexer, T_RPAREN))
3041 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3047 matrix_lvalue_destroy (lvalue);
3052 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3053 enum index_type index_type, size_t other_size,
3054 struct index_vector *iv)
3059 m = matrix_expr_evaluate (e);
3066 bool ok = matrix_normalize_index_vector (m, size, index_type,
3068 gsl_matrix_free (m);
3073 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3074 struct index_vector *iv, gsl_matrix *sm)
3076 gsl_vector dv = to_vector (lvalue->var->value);
3078 /* Convert source matrix 'sm' to source vector 'sv'. */
3079 if (!is_vector (sm))
3081 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
3082 sm->size1, sm->size2);
3085 gsl_vector sv = to_vector (sm);
3086 if (iv->n != sv.size)
3088 msg (SE, _("Can't assign %zu-element vector "
3089 "to %zu-element subvector."), sv.size, iv->n);
3093 /* Assign elements. */
3094 for (size_t x = 0; x < iv->n; x++)
3095 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3100 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3101 struct index_vector *iv0,
3102 struct index_vector *iv1, gsl_matrix *sm)
3104 gsl_matrix *dm = lvalue->var->value;
3106 /* Convert source matrix 'sm' to source vector 'sv'. */
3107 if (iv0->n != sm->size1)
3109 msg (SE, _("Row index vector for assignment to %s has %zu elements "
3110 "but source matrix has %zu rows."),
3111 lvalue->var->name, iv0->n, sm->size1);
3114 if (iv1->n != sm->size2)
3116 msg (SE, _("Column index vector for assignment to %s has %zu elements "
3117 "but source matrix has %zu columns."),
3118 lvalue->var->name, iv1->n, sm->size2);
3122 /* Assign elements. */
3123 for (size_t y = 0; y < iv0->n; y++)
3125 size_t f0 = iv0->indexes[y];
3126 for (size_t x = 0; x < iv1->n; x++)
3128 size_t f1 = iv1->indexes[x];
3129 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3135 /* Takes ownership of M. */
3137 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3138 struct index_vector *iv0, struct index_vector *iv1,
3141 if (!lvalue->n_indexes)
3143 gsl_matrix_free (lvalue->var->value);
3144 lvalue->var->value = sm;
3149 bool ok = (lvalue->n_indexes == 1
3150 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3151 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3152 free (iv0->indexes);
3153 free (iv1->indexes);
3154 gsl_matrix_free (sm);
3160 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3161 struct index_vector *iv0,
3162 struct index_vector *iv1)
3164 *iv0 = INDEX_VECTOR_INIT;
3165 *iv1 = INDEX_VECTOR_INIT;
3166 if (!lvalue->n_indexes)
3169 /* Validate destination matrix exists and has the right shape. */
3170 gsl_matrix *dm = lvalue->var->value;
3173 msg (SE, _("Undefined variable %s."), lvalue->var->name);
3176 else if (dm->size1 == 0 || dm->size2 == 0)
3178 msg (SE, _("Cannot index %zu×%zu matrix."),
3179 dm->size1, dm->size2);
3182 else if (lvalue->n_indexes == 1)
3184 if (!is_vector (dm))
3186 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
3187 dm->size1, dm->size2, lvalue->var->name);
3190 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3191 MAX (dm->size1, dm->size2),
3196 assert (lvalue->n_indexes == 2);
3197 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3198 IV_ROW, dm->size2, iv0))
3201 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3202 IV_COLUMN, dm->size1, iv1))
3204 free (iv0->indexes);
3211 /* Takes ownership of M. */
3213 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3215 struct index_vector iv0, iv1;
3216 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3218 gsl_matrix_free (sm);
3222 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3226 /* Matrix command. */
3230 struct matrix_cmd **commands;
3234 static void matrix_cmds_uninit (struct matrix_cmds *);
3238 enum matrix_cmd_type
3261 struct compute_command
3263 struct matrix_lvalue *lvalue;
3264 struct matrix_expr *rvalue;
3268 struct print_command
3270 struct matrix_expr *expression;
3271 bool use_default_format;
3272 struct fmt_spec format;
3274 int space; /* -1 means NEWPAGE. */
3278 struct string_array literals; /* CLABELS/RLABELS. */
3279 struct matrix_expr *expr; /* CNAMES/RNAMES. */
3285 struct do_if_command
3289 struct matrix_expr *condition;
3290 struct matrix_cmds commands;
3300 struct matrix_var *var;
3301 struct matrix_expr *start;
3302 struct matrix_expr *end;
3303 struct matrix_expr *increment;
3305 /* Loop conditions. */
3306 struct matrix_expr *top_condition;
3307 struct matrix_expr *bottom_condition;
3310 struct matrix_cmds commands;
3314 struct display_command
3316 struct matrix_state *state;
3322 struct release_command
3324 struct matrix_var **vars;
3331 struct matrix_expr *expression;
3332 struct file_handle *outfile;
3333 struct string_array *variables;
3334 struct matrix_expr *names;
3335 struct stringi_set strings;
3341 struct read_file *rf;
3342 struct matrix_lvalue *dst;
3343 struct matrix_expr *size;
3345 enum fmt_type format;
3353 struct write_command
3355 struct matrix_expr *expression;
3356 struct file_handle *outfile;
3359 enum fmt_type format;
3363 bool hold; /* XXX */
3369 struct matrix_lvalue *dst;
3370 struct file_handle *file;
3372 struct string_array variables;
3373 struct matrix_var *names;
3375 /* Treatment of missing values. */
3380 MGET_ERROR, /* Flag error on command. */
3381 MGET_ACCEPT, /* Accept without change, user-missing only. */
3382 MGET_OMIT, /* Drop this case. */
3383 MGET_RECODE /* Recode to 'substitute'. */
3386 double substitute; /* MGET_RECODE only. */
3392 struct msave_command
3394 struct msave_common *common;
3396 struct matrix_expr *expr;
3397 const char *rowtype;
3398 struct matrix_expr *factors;
3399 struct matrix_expr *splits;
3405 struct matrix_state *state;
3406 struct file_handle *file;
3407 struct stringi_set rowtypes;
3411 struct eigen_command
3413 struct matrix_expr *expr;
3414 struct matrix_var *evec;
3415 struct matrix_var *eval;
3419 struct setdiag_command
3421 struct matrix_var *dst;
3422 struct matrix_expr *expr;
3428 struct matrix_expr *expr;
3429 struct matrix_var *u;
3430 struct matrix_var *s;
3431 struct matrix_var *v;
3437 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3438 static bool matrix_cmd_execute (struct matrix_cmd *);
3439 static void matrix_cmd_destroy (struct matrix_cmd *);
3442 static struct matrix_cmd *
3443 matrix_parse_compute (struct matrix_state *s)
3445 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3446 *cmd = (struct matrix_cmd) {
3447 .type = MCMD_COMPUTE,
3448 .compute = { .lvalue = NULL },
3451 struct compute_command *compute = &cmd->compute;
3452 compute->lvalue = matrix_lvalue_parse (s);
3453 if (!compute->lvalue)
3456 if (!lex_force_match (s->lexer, T_EQUALS))
3459 compute->rvalue = matrix_parse_expr (s);
3460 if (!compute->rvalue)
3466 matrix_cmd_destroy (cmd);
3471 matrix_cmd_execute_compute (struct compute_command *compute)
3473 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3477 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3481 print_labels_destroy (struct print_labels *labels)
3485 string_array_destroy (&labels->literals);
3486 matrix_expr_destroy (labels->expr);
3492 parse_literal_print_labels (struct matrix_state *s,
3493 struct print_labels **labelsp)
3495 lex_match (s->lexer, T_EQUALS);
3496 print_labels_destroy (*labelsp);
3497 *labelsp = xzalloc (sizeof **labelsp);
3498 while (lex_token (s->lexer) != T_SLASH
3499 && lex_token (s->lexer) != T_ENDCMD
3500 && lex_token (s->lexer) != T_STOP)
3502 struct string label = DS_EMPTY_INITIALIZER;
3503 while (lex_token (s->lexer) != T_COMMA
3504 && lex_token (s->lexer) != T_SLASH
3505 && lex_token (s->lexer) != T_ENDCMD
3506 && lex_token (s->lexer) != T_STOP)
3508 if (!ds_is_empty (&label))
3509 ds_put_byte (&label, ' ');
3511 if (lex_token (s->lexer) == T_STRING)
3512 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3515 char *rep = lex_next_representation (s->lexer, 0, 0);
3516 ds_put_cstr (&label, rep);
3521 string_array_append_nocopy (&(*labelsp)->literals,
3522 ds_steal_cstr (&label));
3524 if (!lex_match (s->lexer, T_COMMA))
3530 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3532 lex_match (s->lexer, T_EQUALS);
3533 struct matrix_expr *e = matrix_parse_exp (s);
3537 print_labels_destroy (*labelsp);
3538 *labelsp = xzalloc (sizeof **labelsp);
3539 (*labelsp)->expr = e;
3543 static struct matrix_cmd *
3544 matrix_parse_print (struct matrix_state *s)
3546 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3547 *cmd = (struct matrix_cmd) {
3550 .use_default_format = true,
3554 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3557 for (size_t i = 0; ; i++)
3559 enum token_type t = lex_next_token (s->lexer, i);
3560 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3562 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3564 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3567 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3572 cmd->print.expression = matrix_parse_exp (s);
3573 if (!cmd->print.expression)
3577 while (lex_match (s->lexer, T_SLASH))
3579 if (lex_match_id (s->lexer, "FORMAT"))
3581 lex_match (s->lexer, T_EQUALS);
3582 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3584 cmd->print.use_default_format = false;
3586 else if (lex_match_id (s->lexer, "TITLE"))
3588 lex_match (s->lexer, T_EQUALS);
3589 if (!lex_force_string (s->lexer))
3591 free (cmd->print.title);
3592 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3595 else if (lex_match_id (s->lexer, "SPACE"))
3597 lex_match (s->lexer, T_EQUALS);
3598 if (lex_match_id (s->lexer, "NEWPAGE"))
3599 cmd->print.space = -1;
3600 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3602 cmd->print.space = lex_integer (s->lexer);
3608 else if (lex_match_id (s->lexer, "RLABELS"))
3609 parse_literal_print_labels (s, &cmd->print.rlabels);
3610 else if (lex_match_id (s->lexer, "CLABELS"))
3611 parse_literal_print_labels (s, &cmd->print.clabels);
3612 else if (lex_match_id (s->lexer, "RNAMES"))
3614 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3617 else if (lex_match_id (s->lexer, "CNAMES"))
3619 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3624 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3625 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3633 matrix_cmd_destroy (cmd);
3638 matrix_is_integer (const gsl_matrix *m)
3640 for (size_t y = 0; y < m->size1; y++)
3641 for (size_t x = 0; x < m->size2; x++)
3643 double d = gsl_matrix_get (m, y, x);
3651 matrix_max_magnitude (const gsl_matrix *m)
3654 for (size_t y = 0; y < m->size1; y++)
3655 for (size_t x = 0; x < m->size2; x++)
3657 double d = fabs (gsl_matrix_get (m, y, x));
3665 format_fits (struct fmt_spec format, double x)
3667 char *s = data_out (&(union value) { .f = x }, NULL,
3668 &format, settings_get_fmt_settings ());
3669 bool fits = *s != '*' && !strchr (s, 'E');
3674 static struct fmt_spec
3675 default_format (const gsl_matrix *m, int *log_scale)
3679 double max = matrix_max_magnitude (m);
3681 if (matrix_is_integer (m))
3682 for (int w = 1; w <= 12; w++)
3684 struct fmt_spec format = { .type = FMT_F, .w = w };
3685 if (format_fits (format, -max))
3689 if (max >= 1e9 || max <= 1e-4)
3692 snprintf (s, sizeof s, "%.1e", max);
3694 const char *e = strchr (s, 'e');
3696 *log_scale = atoi (e + 1);
3699 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3703 trimmed_string (double d)
3705 struct substring s = ss_buffer ((char *) &d, sizeof d);
3706 ss_rtrim (&s, ss_cstr (" "));
3707 return ss_xstrdup (s);
3710 static struct string_array *
3711 print_labels_get (const struct print_labels *labels, size_t n,
3712 const char *prefix, bool truncate)
3717 struct string_array *out = xzalloc (sizeof *out);
3718 if (labels->literals.n)
3719 string_array_clone (out, &labels->literals);
3720 else if (labels->expr)
3722 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3723 if (m && is_vector (m))
3725 gsl_vector v = to_vector (m);
3726 for (size_t i = 0; i < v.size; i++)
3727 string_array_append_nocopy (out, trimmed_string (
3728 gsl_vector_get (&v, i)));
3730 gsl_matrix_free (m);
3736 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3739 string_array_append (out, "");
3742 string_array_delete (out, out->n - 1);
3745 for (size_t i = 0; i < out->n; i++)
3747 char *s = out->strings[i];
3748 s[strnlen (s, 8)] = '\0';
3755 matrix_cmd_print_space (int space)
3758 output_item_submit (page_break_item_create ());
3759 for (int i = 0; i < space; i++)
3760 output_log ("%s", "");
3764 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3765 struct fmt_spec format, int log_scale)
3767 matrix_cmd_print_space (print->space);
3769 output_log ("%s", print->title);
3771 output_log (" 10 ** %d X", log_scale);
3773 struct string_array *clabels = print_labels_get (print->clabels,
3774 m->size2, "col", true);
3775 if (clabels && format.w < 8)
3777 struct string_array *rlabels = print_labels_get (print->rlabels,
3778 m->size1, "row", true);
3782 struct string line = DS_EMPTY_INITIALIZER;
3784 ds_put_byte_multiple (&line, ' ', 8);
3785 for (size_t x = 0; x < m->size2; x++)
3786 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3787 output_log_nocopy (ds_steal_cstr (&line));
3790 double scale = pow (10.0, log_scale);
3791 bool numeric = fmt_is_numeric (format.type);
3792 for (size_t y = 0; y < m->size1; y++)
3794 struct string line = DS_EMPTY_INITIALIZER;
3796 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3798 for (size_t x = 0; x < m->size2; x++)
3800 double f = gsl_matrix_get (m, y, x);
3802 ? data_out (&(union value) { .f = f / scale}, NULL,
3803 &format, settings_get_fmt_settings ())
3804 : trimmed_string (f));
3805 ds_put_format (&line, " %s", s);
3808 output_log_nocopy (ds_steal_cstr (&line));
3811 string_array_destroy (rlabels);
3813 string_array_destroy (clabels);
3818 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3819 const struct print_labels *print_labels, size_t n,
3820 const char *name, const char *prefix)
3822 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3824 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3825 for (size_t i = 0; i < n; i++)
3826 pivot_category_create_leaf (
3828 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3829 : pivot_value_new_integer (i + 1)));
3831 d->hide_all_labels = true;
3832 string_array_destroy (labels);
3837 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3838 struct fmt_spec format, int log_scale)
3840 struct pivot_table *table = pivot_table_create__ (
3841 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3844 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3846 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3847 N_("Columns"), "col");
3849 struct pivot_footnote *footnote = NULL;
3852 char *s = xasprintf ("× 10**%d", log_scale);
3853 footnote = pivot_table_create_footnote (
3854 table, pivot_value_new_user_text_nocopy (s));
3857 double scale = pow (10.0, log_scale);
3858 bool numeric = fmt_is_numeric (format.type);
3859 for (size_t y = 0; y < m->size1; y++)
3860 for (size_t x = 0; x < m->size2; x++)
3862 double f = gsl_matrix_get (m, y, x);
3863 struct pivot_value *v;
3866 v = pivot_value_new_number (f / scale);
3867 v->numeric.format = format;
3870 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3872 pivot_value_add_footnote (v, footnote);
3873 pivot_table_put2 (table, y, x, v);
3876 pivot_table_submit (table);
3880 matrix_cmd_execute_print (const struct print_command *print)
3882 if (print->expression)
3884 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3889 struct fmt_spec format = (print->use_default_format
3890 ? default_format (m, &log_scale)
3893 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3894 matrix_cmd_print_text (print, m, format, log_scale);
3896 matrix_cmd_print_tables (print, m, format, log_scale);
3898 gsl_matrix_free (m);
3902 matrix_cmd_print_space (print->space);
3904 output_log ("%s", print->title);
3911 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3912 const char *command_name,
3913 const char *stop1, const char *stop2)
3915 lex_end_of_command (s->lexer);
3916 lex_discard_rest_of_command (s->lexer);
3918 size_t allocated = 0;
3921 while (lex_token (s->lexer) == T_ENDCMD)
3924 if (lex_at_phrase (s->lexer, stop1)
3925 || (stop2 && lex_at_phrase (s->lexer, stop2)))
3928 if (lex_at_phrase (s->lexer, "END MATRIX"))
3930 msg (SE, _("Premature END MATRIX within %s."), command_name);
3934 struct matrix_cmd *cmd = matrix_parse_command (s);
3938 if (c->n >= allocated)
3939 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
3940 c->commands[c->n++] = cmd;
3945 matrix_parse_do_if_clause (struct matrix_state *s,
3946 struct do_if_command *ifc,
3947 bool parse_condition,
3948 size_t *allocated_clauses)
3950 if (ifc->n_clauses >= *allocated_clauses)
3951 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
3952 sizeof *ifc->clauses);
3953 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
3954 *c = (struct do_if_clause) { .condition = NULL };
3956 if (parse_condition)
3958 c->condition = matrix_parse_expr (s);
3963 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
3966 static struct matrix_cmd *
3967 matrix_parse_do_if (struct matrix_state *s)
3969 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3970 *cmd = (struct matrix_cmd) {
3972 .do_if = { .n_clauses = 0 }
3975 size_t allocated_clauses = 0;
3978 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
3981 while (lex_match_phrase (s->lexer, "ELSE IF"));
3983 if (lex_match_id (s->lexer, "ELSE")
3984 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
3987 if (!lex_match_phrase (s->lexer, "END IF"))
3992 matrix_cmd_destroy (cmd);
3997 matrix_cmd_execute_do_if (struct do_if_command *cmd)
3999 for (size_t i = 0; i < cmd->n_clauses; i++)
4001 struct do_if_clause *c = &cmd->clauses[i];
4005 if (!matrix_expr_evaluate_scalar (c->condition,
4006 i ? "ELSE IF" : "DO IF",
4011 for (size_t j = 0; j < c->commands.n; j++)
4012 if (!matrix_cmd_execute (c->commands.commands[j]))
4019 static struct matrix_cmd *
4020 matrix_parse_loop (struct matrix_state *s)
4022 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4023 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4025 struct loop_command *loop = &cmd->loop;
4026 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4028 struct substring name = lex_tokss (s->lexer);
4029 loop->var = matrix_var_lookup (s, name);
4031 loop->var = matrix_var_create (s, name);
4036 loop->start = matrix_parse_expr (s);
4037 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4040 loop->end = matrix_parse_expr (s);
4044 if (lex_match (s->lexer, T_BY))
4046 loop->increment = matrix_parse_expr (s);
4047 if (!loop->increment)
4052 if (lex_match_id (s->lexer, "IF"))
4054 loop->top_condition = matrix_parse_expr (s);
4055 if (!loop->top_condition)
4059 bool was_in_loop = s->in_loop;
4061 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4063 s->in_loop = was_in_loop;
4067 if (!lex_match_phrase (s->lexer, "END LOOP"))
4070 if (lex_match_id (s->lexer, "IF"))
4072 loop->bottom_condition = matrix_parse_expr (s);
4073 if (!loop->bottom_condition)
4080 matrix_cmd_destroy (cmd);
4085 matrix_cmd_execute_loop (struct loop_command *cmd)
4089 long int increment = 1;
4092 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4093 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4095 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4099 if (increment > 0 ? value > end
4100 : increment < 0 ? value < end
4105 int mxloops = settings_get_mxloops ();
4106 for (int i = 0; i < mxloops; i++)
4111 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4113 gsl_matrix_free (cmd->var->value);
4114 cmd->var->value = NULL;
4116 if (!cmd->var->value)
4117 cmd->var->value = gsl_matrix_alloc (1, 1);
4119 gsl_matrix_set (cmd->var->value, 0, 0, value);
4122 if (cmd->top_condition)
4125 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4130 for (size_t j = 0; j < cmd->commands.n; j++)
4131 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4134 if (cmd->bottom_condition)
4137 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4138 "END LOOP IF", &d) || d > 0)
4145 if (increment > 0 ? value > end : value < end)
4151 static struct matrix_cmd *
4152 matrix_parse_break (struct matrix_state *s)
4156 msg (SE, _("BREAK not inside LOOP."));
4160 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4161 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4165 static struct matrix_cmd *
4166 matrix_parse_display (struct matrix_state *s)
4168 bool dictionary = false;
4169 bool status = false;
4170 while (lex_token (s->lexer) == T_ID)
4172 if (lex_match_id (s->lexer, "DICTIONARY"))
4174 else if (lex_match_id (s->lexer, "STATUS"))
4178 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4182 if (!dictionary && !status)
4183 dictionary = status = true;
4185 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4186 *cmd = (struct matrix_cmd) {
4187 .type = MCMD_DISPLAY,
4190 .dictionary = dictionary,
4198 compare_matrix_var_pointers (const void *a_, const void *b_)
4200 const struct matrix_var *const *ap = a_;
4201 const struct matrix_var *const *bp = b_;
4202 const struct matrix_var *a = *ap;
4203 const struct matrix_var *b = *bp;
4204 return strcmp (a->name, b->name);
4208 matrix_cmd_execute_display (struct display_command *cmd)
4210 const struct matrix_state *s = cmd->state;
4212 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4213 pivot_dimension_create (
4214 table, PIVOT_AXIS_COLUMN, N_("Property"),
4215 N_("Rows"), N_("Columns"), N_("Size (kB)"));
4217 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4220 const struct matrix_var *var;
4221 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4223 vars[n_vars++] = var;
4224 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4226 struct pivot_dimension *rows = pivot_dimension_create (
4227 table, PIVOT_AXIS_ROW, N_("Variable"));
4228 for (size_t i = 0; i < n_vars; i++)
4230 const struct matrix_var *var = vars[i];
4231 pivot_category_create_leaf (
4232 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4234 size_t r = var->value->size1;
4235 size_t c = var->value->size2;
4236 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4237 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4238 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4240 pivot_table_submit (table);
4243 static struct matrix_cmd *
4244 matrix_parse_release (struct matrix_state *s)
4246 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4247 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4249 size_t allocated_vars = 0;
4250 while (lex_token (s->lexer) == T_ID)
4252 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4255 if (cmd->release.n_vars >= allocated_vars)
4256 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4257 sizeof *cmd->release.vars);
4258 cmd->release.vars[cmd->release.n_vars++] = var;
4261 lex_error (s->lexer, _("Variable name expected."));
4264 if (!lex_match (s->lexer, T_COMMA))
4272 matrix_cmd_execute_release (struct release_command *cmd)
4274 for (size_t i = 0; i < cmd->n_vars; i++)
4276 struct matrix_var *var = cmd->vars[i];
4277 gsl_matrix_free (var->value);
4282 static struct matrix_cmd *
4283 matrix_parse_save (struct matrix_state *s)
4285 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4286 *cmd = (struct matrix_cmd) {
4289 .strings = STRINGI_SET_INITIALIZER (cmd->save.strings)
4293 struct save_command *save = &cmd->save;
4294 save->expression = matrix_parse_exp (s);
4295 if (!save->expression)
4298 while (lex_match (s->lexer, T_SLASH))
4300 if (lex_match_id (s->lexer, "OUTFILE"))
4302 lex_match (s->lexer, T_EQUALS);
4304 fh_unref (save->outfile);
4305 save->outfile = (lex_match (s->lexer, T_ASTERISK)
4307 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4311 else if (lex_match_id (s->lexer, "VARIABLES"))
4313 lex_match (s->lexer, T_EQUALS);
4317 struct dictionary *d = dict_create (get_default_encoding ());
4318 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4319 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4324 string_array_destroy (save->variables);
4325 if (!save->variables)
4326 save->variables = xmalloc (sizeof *save->variables);
4327 *save->variables = (struct string_array) {
4333 else if (lex_match_id (s->lexer, "NAMES"))
4335 lex_match (s->lexer, T_EQUALS);
4336 matrix_expr_destroy (save->names);
4337 save->names = matrix_parse_exp (s);
4341 else if (lex_match_id (s->lexer, "STRINGS"))
4343 lex_match (s->lexer, T_EQUALS);
4344 while (lex_token (s->lexer) == T_ID)
4346 stringi_set_insert (&save->strings, lex_tokcstr (s->lexer));
4348 if (!lex_match (s->lexer, T_COMMA))
4354 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4362 if (s->prev_save_outfile)
4363 save->outfile = fh_ref (s->prev_save_outfile);
4366 lex_sbc_missing ("OUTFILE");
4370 fh_unref (s->prev_save_outfile);
4371 s->prev_save_outfile = fh_ref (save->outfile);
4373 if (save->variables && save->names)
4375 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4376 matrix_expr_destroy (save->names);
4383 matrix_cmd_destroy (cmd);
4388 matrix_cmd_execute_save (const struct save_command *save)
4390 assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */
4392 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4397 struct dictionary *dict = dict_create (get_default_encoding ());
4398 struct string_array names = { .n = 0 };
4399 if (save->variables)
4400 string_array_clone (&names, save->variables);
4401 else if (save->names)
4403 gsl_matrix *nm = matrix_expr_evaluate (save->names);
4404 if (nm && is_vector (nm))
4406 gsl_vector nv = to_vector (nm);
4407 for (size_t i = 0; i < nv.size; i++)
4409 char *name = trimmed_string (gsl_vector_get (&nv, i));
4410 if (dict_id_is_valid (dict, name, true))
4411 string_array_append_nocopy (&names, name);
4418 struct stringi_set strings;
4419 stringi_set_clone (&strings, &save->strings);
4421 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4426 name = names.strings[i];
4429 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4433 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4434 struct variable *var = dict_create_var (dict, name, width);
4437 msg (SE, _("Duplicate variable name %s in SAVE statement."),
4443 if (!stringi_set_is_empty (&strings))
4445 const char *example = stringi_set_node_get_string (
4446 stringi_set_first (&strings));
4447 msg (SE, ngettext ("STRINGS specified a variable %s, but no variable "
4448 "with that name was found on SAVE.",
4449 "STRINGS specified %2$zu variables, including %1$s, "
4450 "whose names were not found on SAVE.",
4451 stringi_set_count (&strings)),
4452 example, stringi_set_count (&strings));
4455 stringi_set_destroy (&strings);
4456 string_array_destroy (&names);
4460 gsl_matrix_free (m);
4465 struct casewriter *writer = any_writer_open (save->outfile, dict);
4468 gsl_matrix_free (m);
4473 const struct caseproto *proto = dict_get_proto (dict);
4474 for (size_t y = 0; y < m->size1; y++)
4476 struct ccase *c = case_create (proto);
4477 for (size_t x = 0; x < m->size2; x++)
4479 double d = gsl_matrix_get (m, y, x);
4480 int width = caseproto_get_width (proto, x);
4481 union value *value = case_data_rw_idx (c, x);
4485 memcpy (value->s, &d, width);
4487 casewriter_write (writer, c);
4489 casewriter_destroy (writer);
4491 gsl_matrix_free (m);
4495 static struct matrix_cmd *
4496 matrix_parse_read (struct matrix_state *s)
4498 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4499 *cmd = (struct matrix_cmd) {
4501 .read = { .format = FMT_F },
4504 struct file_handle *fh = NULL;
4505 char *encoding = NULL;
4506 struct read_command *read = &cmd->read;
4507 read->dst = matrix_lvalue_parse (s);
4512 int repetitions = 0;
4513 int record_width = 0;
4514 bool seen_format = false;
4515 while (lex_match (s->lexer, T_SLASH))
4517 if (lex_match_id (s->lexer, "FILE"))
4519 lex_match (s->lexer, T_EQUALS);
4522 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4526 else if (lex_match_id (s->lexer, "ENCODING"))
4528 lex_match (s->lexer, T_EQUALS);
4529 if (!lex_force_string (s->lexer))
4533 encoding = ss_xstrdup (lex_tokss (s->lexer));
4537 else if (lex_match_id (s->lexer, "FIELD"))
4539 lex_match (s->lexer, T_EQUALS);
4541 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4543 read->c1 = lex_integer (s->lexer);
4545 if (!lex_force_match (s->lexer, T_TO)
4546 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4548 read->c2 = lex_integer (s->lexer) + 1;
4551 record_width = read->c2 - read->c1;
4552 if (lex_match (s->lexer, T_BY))
4554 if (!lex_force_int_range (s->lexer, "BY", 1,
4555 read->c2 - read->c1))
4557 by = lex_integer (s->lexer);
4560 if (record_width % by)
4562 msg (SE, _("BY %d does not evenly divide record width %d."),
4570 else if (lex_match_id (s->lexer, "SIZE"))
4572 lex_match (s->lexer, T_EQUALS);
4573 read->size = matrix_parse_exp (s);
4577 else if (lex_match_id (s->lexer, "MODE"))
4579 lex_match (s->lexer, T_EQUALS);
4580 if (lex_match_id (s->lexer, "RECTANGULAR"))
4581 read->symmetric = false;
4582 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4583 read->symmetric = true;
4586 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4590 else if (lex_match_id (s->lexer, "REREAD"))
4591 read->reread = true;
4592 else if (lex_match_id (s->lexer, "FORMAT"))
4596 lex_sbc_only_once ("FORMAT");
4601 lex_match (s->lexer, T_EQUALS);
4603 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4606 const char *p = lex_tokcstr (s->lexer);
4607 if (c_isdigit (p[0]))
4609 repetitions = atoi (p);
4610 p += strspn (p, "0123456789");
4611 if (!fmt_from_name (p, &read->format))
4613 lex_error (s->lexer, _("Unknown format %s."), p);
4618 else if (!fmt_from_name (p, &read->format))
4620 struct fmt_spec format;
4621 if (!parse_format_specifier (s->lexer, &format))
4623 read->format = format.type;
4630 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4631 "REREAD", "FORMAT");
4638 lex_sbc_missing ("FIELD");
4642 if (!read->dst->n_indexes && !read->size)
4644 msg (SE, _("SIZE is required for reading data into a full matrix "
4645 "(as opposed to a submatrix)."));
4651 if (s->prev_read_file)
4652 fh = fh_ref (s->prev_read_file);
4655 lex_sbc_missing ("FILE");
4659 fh_unref (s->prev_read_file);
4660 s->prev_read_file = fh_ref (fh);
4662 read->rf = read_file_create (s, fh);
4665 free (read->rf->encoding);
4666 read->rf->encoding = encoding;
4670 /* Field width may be specified in multiple ways:
4673 2. The format on FORMAT.
4674 3. The repetition factor on FORMAT.
4676 (2) and (3) are mutually exclusive.
4678 If more than one of these is present, they must agree. If none of them is
4679 present, then free-field format is used.
4681 if (repetitions > record_width)
4683 msg (SE, _("%d repetitions cannot fit in record width %d."),
4684 repetitions, record_width);
4687 int w = (repetitions ? record_width / repetitions
4693 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4694 "which implies field width %d, "
4695 "but BY specifies field width %d."),
4696 repetitions, record_width, w, by);
4698 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4707 matrix_cmd_destroy (cmd);
4713 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4714 gsl_matrix *m, struct substring p, size_t y, size_t x)
4716 const char *input_encoding = dfm_reader_get_encoding (reader);
4718 char *error = data_in (p, input_encoding, read->format,
4719 settings_get_fmt_settings (), &v, 0, NULL);
4720 /* XXX report error if value is missing */
4722 msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
4725 gsl_matrix_set (m, y, x, v.f);
4726 if (read->symmetric && x != y)
4727 gsl_matrix_set (m, x, y, v.f);
4732 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4733 struct substring *line)
4735 if (dfm_eof (reader))
4737 msg (SE, _("Unexpected end of file reading matrix data."));
4740 dfm_expand_tabs (reader);
4741 *line = ss_substr (dfm_get_record (reader),
4742 read->c1 - 1, read->c2 - read->c1);
4747 matrix_read (struct read_command *read, struct dfm_reader *reader,
4750 for (size_t y = 0; y < m->size1; y++)
4752 size_t nx = read->symmetric ? y + 1 : m->size2;
4754 struct substring line = ss_empty ();
4755 for (size_t x = 0; x < nx; x++)
4762 ss_ltrim (&line, ss_cstr (" ,"));
4763 if (!ss_is_empty (line))
4765 if (!matrix_read_line (read, reader, &line))
4767 dfm_forward_record (reader);
4770 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
4774 if (!matrix_read_line (read, reader, &line))
4776 size_t fields_per_line = (read->c2 - read->c1) / read->w;
4777 int f = x % fields_per_line;
4778 if (f == fields_per_line - 1)
4779 dfm_forward_record (reader);
4781 p = ss_substr (line, read->w * f, read->w);
4784 matrix_read_set_field (read, reader, m, p, y, x);
4788 dfm_forward_record (reader);
4791 ss_ltrim (&line, ss_cstr (" ,"));
4792 if (!ss_is_empty (line))
4793 msg (SW, _("Trailing garbage on line \"%.*s\""),
4794 (int) line.length, line.string);
4800 matrix_cmd_execute_read (struct read_command *read)
4802 struct index_vector iv0, iv1;
4803 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
4806 size_t size[2] = { SIZE_MAX, SIZE_MAX };
4809 gsl_matrix *m = matrix_expr_evaluate (read->size);
4815 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4816 "not a %zu×%zu matrix."), m->size1, m->size2);
4817 gsl_matrix_free (m);
4823 gsl_vector v = to_vector (m);
4827 d[0] = gsl_vector_get (&v, 0);
4830 else if (v.size == 2)
4832 d[0] = gsl_vector_get (&v, 0);
4833 d[1] = gsl_vector_get (&v, 1);
4837 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4838 "not a %zu×%zu matrix."),
4839 m->size1, m->size2),
4840 gsl_matrix_free (m);
4845 gsl_matrix_free (m);
4847 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
4849 msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]);
4859 if (read->dst->n_indexes)
4861 size_t submatrix_size[2];
4862 if (read->dst->n_indexes == 2)
4864 submatrix_size[0] = iv0.n;
4865 submatrix_size[1] = iv1.n;
4867 else if (read->dst->var->value->size1 == 1)
4869 submatrix_size[0] = 1;
4870 submatrix_size[1] = iv0.n;
4874 submatrix_size[0] = iv0.n;
4875 submatrix_size[1] = 1;
4880 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
4882 msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
4885 submatrix_size[0], submatrix_size[1]);
4893 size[0] = submatrix_size[0];
4894 size[1] = submatrix_size[1];
4898 struct dfm_reader *reader = read_file_open (read->rf);
4900 dfm_reread_record (reader, 1);
4902 if (read->symmetric && size[0] != size[1])
4904 msg (SE, _("Cannot read non-square %zu×%zu matrix "
4905 "using READ with MODE=SYMMETRIC."),
4911 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
4912 matrix_read (read, reader, tmp);
4913 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
4916 static struct matrix_cmd *
4917 matrix_parse_write (struct matrix_state *s)
4919 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4920 *cmd = (struct matrix_cmd) {
4922 .write = { .format = FMT_F },
4925 struct write_command *write = &cmd->write;
4926 write->expression = matrix_parse_exp (s);
4927 if (!write->expression)
4931 int repetitions = 0;
4932 int record_width = 0;
4933 bool seen_format = false;
4934 while (lex_match (s->lexer, T_SLASH))
4936 if (lex_match_id (s->lexer, "OUTFILE"))
4938 lex_match (s->lexer, T_EQUALS);
4940 fh_unref (write->outfile);
4941 write->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
4942 if (!write->outfile)
4945 else if (lex_match_id (s->lexer, "ENCODING"))
4947 lex_match (s->lexer, T_EQUALS);
4948 if (!lex_force_string (s->lexer))
4951 free (write->encoding);
4952 write->encoding = ss_xstrdup (lex_tokss (s->lexer));
4956 else if (lex_match_id (s->lexer, "FIELD"))
4958 lex_match (s->lexer, T_EQUALS);
4960 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4962 write->c1 = lex_integer (s->lexer);
4964 if (!lex_force_match (s->lexer, T_TO)
4965 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
4967 write->c2 = lex_integer (s->lexer) + 1;
4970 record_width = write->c2 - write->c1;
4971 if (lex_match (s->lexer, T_BY))
4973 if (!lex_force_int_range (s->lexer, "BY", 1,
4974 write->c2 - write->c1))
4976 by = lex_integer (s->lexer);
4979 if (record_width % by)
4981 msg (SE, _("BY %d does not evenly divide record width %d."),
4989 else if (lex_match_id (s->lexer, "MODE"))
4991 lex_match (s->lexer, T_EQUALS);
4992 if (lex_match_id (s->lexer, "RECTANGULAR"))
4993 write->triangular = false;
4994 else if (lex_match_id (s->lexer, "TRIANGULAR"))
4995 write->triangular = true;
4998 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5002 else if (lex_match_id (s->lexer, "HOLD"))
5004 else if (lex_match_id (s->lexer, "FORMAT"))
5008 lex_sbc_only_once ("FORMAT");
5013 lex_match (s->lexer, T_EQUALS);
5015 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5018 const char *p = lex_tokcstr (s->lexer);
5019 if (c_isdigit (p[0]))
5021 repetitions = atoi (p);
5022 p += strspn (p, "0123456789");
5023 if (!fmt_from_name (p, &write->format))
5025 lex_error (s->lexer, _("Unknown format %s."), p);
5030 else if (!fmt_from_name (p, &write->format))
5032 struct fmt_spec format;
5033 if (!parse_format_specifier (s->lexer, &format))
5035 write->format = format.type;
5036 write->w = format.w;
5037 write->d = format.d;
5042 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5050 lex_sbc_missing ("FIELD");
5054 if (!write->outfile)
5056 if (s->prev_write_outfile)
5057 write->outfile = fh_ref (s->prev_write_outfile);
5060 lex_sbc_missing ("OUTFILE");
5064 fh_unref (s->prev_write_outfile);
5065 s->prev_write_outfile = fh_ref (write->outfile);
5067 /* Field width may be specified in multiple ways:
5070 2. The format on FORMAT.
5071 3. The repetition factor on FORMAT.
5073 (2) and (3) are mutually exclusive.
5075 If more than one of these is present, they must agree. If none of them is
5076 present, then free-field format is used.
5078 if (repetitions > record_width)
5080 msg (SE, _("%d repetitions cannot fit in record width %d."),
5081 repetitions, record_width);
5084 int w = (repetitions ? record_width / repetitions
5085 : write->w ? write->w
5090 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5091 "which implies field width %d, "
5092 "but BY specifies field width %d."),
5093 repetitions, record_width, w, by);
5095 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5103 matrix_cmd_destroy (cmd);
5108 matrix_cmd_execute_write (struct write_command *write)
5110 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5114 if (write->triangular && m->size1 != m->size2)
5116 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5117 "the matrix to be written has dimensions %zu×%zu."),
5118 m->size1, m->size2);
5119 gsl_matrix_free (m);
5123 struct dfm_writer *writer = dfm_open_writer (write->outfile, write->encoding);
5126 gsl_matrix_free (m);
5130 const struct fmt_settings *settings = settings_get_fmt_settings ();
5131 struct fmt_spec format = {
5132 .type = write->format,
5133 .w = write->w ? write->w : 40,
5136 struct u8_line line = U8_LINE_INITIALIZER;
5137 for (size_t y = 0; y < m->size1; y++)
5139 size_t nx = write->triangular ? y + 1 : m->size2;
5141 for (size_t x = 0; x < nx; x++)
5143 /* XXX string values */
5144 union value v = { .f = gsl_matrix_get (m, y, x) };
5146 ? data_out (&v, NULL, &format, settings)
5147 : data_out_stretchy (&v, NULL, &format, settings, NULL));
5148 size_t len = strlen (s);
5149 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5150 if (width + x0 > write->c2)
5152 dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5153 u8_line_clear (&line);
5156 u8_line_put (&line, x0, x0 + width, s, len);
5159 x0 += write->w ? write->w : width + 1;
5162 dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5163 u8_line_clear (&line);
5165 u8_line_destroy (&line);
5166 dfm_close_writer (writer);
5168 gsl_matrix_free (m);
5171 static struct matrix_cmd *
5172 matrix_parse_get (struct matrix_state *s)
5174 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5175 *cmd = (struct matrix_cmd) {
5178 .user = { .treatment = MGET_ERROR },
5179 .system = { .treatment = MGET_ERROR },
5183 struct get_command *get = &cmd->get;
5184 get->dst = matrix_lvalue_parse (s);
5188 while (lex_match (s->lexer, T_SLASH))
5190 if (lex_match_id (s->lexer, "FILE"))
5192 if (get->variables.n)
5194 lex_error (s->lexer, _("FILE must precede VARIABLES"));
5197 lex_match (s->lexer, T_EQUALS);
5199 fh_unref (get->file);
5200 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5204 else if (lex_match_id (s->lexer, "ENCODING"))
5206 if (get->variables.n)
5208 lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
5211 lex_match (s->lexer, T_EQUALS);
5212 if (!lex_force_string (s->lexer))
5215 free (get->encoding);
5216 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5220 else if (lex_match_id (s->lexer, "VARIABLES"))
5222 lex_match (s->lexer, T_EQUALS);
5224 struct dictionary *dict = NULL;
5227 dict = dataset_dict (s->dataset);
5228 if (dict_get_var_cnt (dict) == 0)
5230 lex_error (s->lexer, _("GET cannot read empty active file."));
5236 struct casereader *reader = any_reader_open_and_decode (
5237 get->file, get->encoding, &dict, NULL);
5240 casereader_destroy (reader);
5243 struct variable **vars;
5245 bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
5246 PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
5253 string_array_clear (&get->variables);
5254 for (size_t i = 0; i < n_vars; i++)
5255 string_array_append (&get->variables, var_get_name (vars[i]));
5259 else if (lex_match_id (s->lexer, "NAMES"))
5261 lex_match (s->lexer, T_EQUALS);
5262 if (!lex_force_id (s->lexer))
5265 struct substring name = lex_tokss (s->lexer);
5266 get->names = matrix_var_lookup (s, name);
5268 get->names = matrix_var_create (s, name);
5271 else if (lex_match_id (s->lexer, "MISSING"))
5273 lex_match (s->lexer, T_EQUALS);
5274 if (lex_match_id (s->lexer, "ACCEPT"))
5275 get->user.treatment = MGET_ACCEPT;
5276 else if (lex_match_id (s->lexer, "OMIT"))
5277 get->user.treatment = MGET_OMIT;
5278 else if (lex_is_number (s->lexer))
5280 get->user.treatment = MGET_RECODE;
5281 get->user.substitute = lex_number (s->lexer);
5286 lex_error (s->lexer, NULL);
5290 else if (lex_match_id (s->lexer, "SYSMIS"))
5292 lex_match (s->lexer, T_EQUALS);
5293 if (lex_match_id (s->lexer, "OMIT"))
5294 get->user.treatment = MGET_OMIT;
5295 else if (lex_is_number (s->lexer))
5297 get->user.treatment = MGET_RECODE;
5298 get->user.substitute = lex_number (s->lexer);
5303 lex_error (s->lexer, NULL);
5309 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5310 "MISSING", "SYSMIS");
5317 matrix_cmd_destroy (cmd);
5322 matrix_cmd_execute_get (struct get_command *get)
5324 assert (get->file); /* XXX */
5326 struct dictionary *dict;
5327 struct casereader *reader = any_reader_open_and_decode (
5328 get->file, get->encoding, &dict, NULL);
5332 const struct variable **vars = xnmalloc (
5333 get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
5337 if (get->variables.n)
5339 for (size_t i = 0; i < get->variables.n; i++)
5341 const char *name = get->variables.strings[i];
5342 const struct variable *var = dict_lookup_var (dict, name);
5345 msg (SE, _("GET: Data file does not contain variable %s."),
5351 if (!var_is_numeric (var))
5353 msg (SE, _("GET: Variable %s is not numeric."), name);
5358 vars[n_vars++] = var;
5363 for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
5365 const struct variable *var = dict_get_var (dict, i);
5366 if (!var_is_numeric (var))
5368 msg (SE, _("GET: Variable %s is not numeric."),
5369 var_get_name (var));
5374 vars[n_vars++] = var;
5379 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5380 long long int casenum = 1;
5382 for (struct ccase *c = casereader_read (reader); c;
5383 c = casereader_read (reader), casenum++)
5385 if (n_rows >= m->size1)
5387 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5388 for (size_t y = 0; y < n_rows; y++)
5389 for (size_t x = 0; x < n_vars; x++)
5390 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5391 gsl_matrix_free (m);
5396 for (size_t x = 0; x < n_vars; x++)
5398 const struct variable *var = vars[x];
5399 double d = case_num (c, var);
5402 if (get->system.treatment == MGET_RECODE)
5403 d = get->system.substitute;
5404 else if (get->system.treatment == MGET_OMIT)
5408 msg (SE, _("GET: Variable %s in case %lld "
5409 "is system-missing."),
5410 var_get_name (var), casenum);
5414 else if (var_is_num_missing (var, d, MV_USER))
5416 if (get->user.treatment == MGET_RECODE)
5417 d = get->user.substitute;
5418 else if (get->user.treatment == MGET_OMIT)
5420 else if (get->user.treatment != MGET_ACCEPT)
5422 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5424 var_get_name (var), casenum, d);
5428 gsl_matrix_set (m, n_rows, x, d);
5436 casereader_destroy (reader);
5440 matrix_lvalue_evaluate_and_assign (get->dst, m);
5443 gsl_matrix_free (m);
5449 match_rowtype (struct lexer *lexer)
5451 static const char *rowtypes[] = {
5452 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5454 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5456 for (size_t i = 0; i < n_rowtypes; i++)
5457 if (lex_match_id (lexer, rowtypes[i]))
5460 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5465 parse_var_names (struct lexer *lexer, struct string_array *sa)
5467 lex_match (lexer, T_EQUALS);
5469 string_array_clear (sa);
5471 struct dictionary *dict = dict_create (get_default_encoding ());
5472 dict_create_var_assert (dict, "ROWTYPE_", 8);
5473 dict_create_var_assert (dict, "VARNAME_", 8);
5476 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5477 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5482 string_array_clear (sa);
5483 sa->strings = names;
5484 sa->n = sa->allocated = n_names;
5490 msave_common_uninit (struct msave_common *common)
5494 fh_unref (common->outfile);
5495 string_array_destroy (&common->variables);
5496 string_array_destroy (&common->fnames);
5497 string_array_destroy (&common->snames);
5502 compare_variables (const char *keyword,
5503 const struct string_array *new,
5504 const struct string_array *old)
5510 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5511 "on the first MSAVE within MATRIX."), keyword);
5514 else if (!string_array_equal_case (old, new))
5516 msg (SE, _("%s must specify the same variables each time within "
5517 "a given MATRIX."), keyword);
5524 static struct matrix_cmd *
5525 matrix_parse_msave (struct matrix_state *s)
5527 struct msave_common common = { .outfile = NULL };
5528 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5529 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5531 struct msave_command *msave = &cmd->msave;
5532 if (lex_token (s->lexer) == T_ID)
5533 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5534 msave->expr = matrix_parse_exp (s);
5538 while (lex_match (s->lexer, T_SLASH))
5540 if (lex_match_id (s->lexer, "TYPE"))
5542 lex_match (s->lexer, T_EQUALS);
5544 msave->rowtype = match_rowtype (s->lexer);
5545 if (!msave->rowtype)
5548 else if (lex_match_id (s->lexer, "OUTFILE"))
5550 lex_match (s->lexer, T_EQUALS);
5552 fh_unref (common.outfile);
5553 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5554 if (!common.outfile)
5557 else if (lex_match_id (s->lexer, "VARIABLES"))
5559 if (!parse_var_names (s->lexer, &common.variables))
5562 else if (lex_match_id (s->lexer, "FNAMES"))
5564 if (!parse_var_names (s->lexer, &common.fnames))
5567 else if (lex_match_id (s->lexer, "SNAMES"))
5569 if (!parse_var_names (s->lexer, &common.snames))
5572 else if (lex_match_id (s->lexer, "SPLIT"))
5574 lex_match (s->lexer, T_EQUALS);
5576 matrix_expr_destroy (msave->splits);
5577 msave->splits = matrix_parse_exp (s);
5581 else if (lex_match_id (s->lexer, "FACTOR"))
5583 lex_match (s->lexer, T_EQUALS);
5585 matrix_expr_destroy (msave->factors);
5586 msave->factors = matrix_parse_exp (s);
5587 if (!msave->factors)
5592 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5593 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5597 if (!msave->rowtype)
5599 lex_sbc_missing ("TYPE");
5602 common.has_splits = msave->splits || common.snames.n;
5603 common.has_factors = msave->factors || common.fnames.n;
5605 struct msave_common *c = s->common ? s->common : &common;
5606 if (c->fnames.n && !msave->factors)
5608 msg (SE, _("FNAMES requires FACTOR."));
5611 if (c->snames.n && !msave->splits)
5613 msg (SE, _("SNAMES requires SPLIT."));
5616 if (c->has_factors && !common.has_factors)
5618 msg (SE, _("%s is required because it was present on the first "
5619 "MSAVE in this MATRIX command."), "FACTOR");
5622 if (c->has_splits && !common.has_splits)
5624 msg (SE, _("%s is required because it was present on the first "
5625 "MSAVE in this MATRIX command."), "SPLIT");
5631 if (!common.outfile)
5633 lex_sbc_missing ("OUTFILE");
5636 s->common = xmemdup (&common, sizeof common);
5642 bool same = common.outfile == s->common->outfile;
5643 fh_unref (common.outfile);
5646 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5647 "within a single MATRIX command."));
5651 if (!compare_variables ("VARIABLES",
5652 &common.variables, &s->common->variables)
5653 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5654 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5656 msave_common_uninit (&common);
5658 msave->common = s->common;
5659 if (!msave->varname_)
5660 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5664 msave_common_uninit (&common);
5665 matrix_cmd_destroy (cmd);
5670 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5672 gsl_matrix *m = matrix_expr_evaluate (e);
5678 msg (SE, _("%s expression must evaluate to vector, "
5679 "not a %zu×%zu matrix."),
5680 name, m->size1, m->size2);
5681 gsl_matrix_free (m);
5685 return matrix_to_vector (m);
5689 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5691 for (size_t i = 0; i < vars->n; i++)
5692 if (!dict_create_var (d, vars->strings[i], 0))
5693 return vars->strings[i];
5697 static struct dictionary *
5698 msave_create_dict (const struct msave_common *common)
5700 struct dictionary *dict = dict_create (get_default_encoding ());
5702 const char *dup_split = msave_add_vars (dict, &common->snames);
5705 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5709 dict_create_var_assert (dict, "ROWTYPE_", 8);
5711 const char *dup_factor = msave_add_vars (dict, &common->fnames);
5714 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5718 dict_create_var_assert (dict, "VARNAME_", 8);
5720 const char *dup_var = msave_add_vars (dict, &common->variables);
5723 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
5735 matrix_cmd_execute_msave (struct msave_command *msave)
5737 struct msave_common *common = msave->common;
5738 gsl_matrix *m = NULL;
5739 gsl_vector *factors = NULL;
5740 gsl_vector *splits = NULL;
5742 m = matrix_expr_evaluate (msave->expr);
5746 if (!common->variables.n)
5747 for (size_t i = 0; i < m->size2; i++)
5748 string_array_append_nocopy (&common->variables,
5749 xasprintf ("COL%zu", i + 1));
5751 if (m->size2 != common->variables.n)
5753 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
5754 m->size2, common->variables.n);
5760 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
5764 if (!common->fnames.n)
5765 for (size_t i = 0; i < factors->size; i++)
5766 string_array_append_nocopy (&common->fnames,
5767 xasprintf ("FAC%zu", i + 1));
5769 if (factors->size != common->fnames.n)
5771 msg (SE, _("There are %zu factor variables, "
5772 "but %zu split values were supplied."),
5773 common->fnames.n, factors->size);
5780 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
5784 if (!common->fnames.n)
5785 for (size_t i = 0; i < splits->size; i++)
5786 string_array_append_nocopy (&common->fnames,
5787 xasprintf ("SPL%zu", i + 1));
5789 if (splits->size != common->fnames.n)
5791 msg (SE, _("There are %zu split variables, "
5792 "but %zu split values were supplied."),
5793 common->fnames.n, splits->size);
5798 if (!common->writer)
5800 struct dictionary *dict = msave_create_dict (common);
5804 common->writer = any_writer_open (common->outfile, dict);
5805 if (!common->writer)
5811 common->dict = dict;
5814 for (size_t y = 0; y < m->size1; y++)
5816 struct ccase *c = case_create (dict_get_proto (common->dict));
5819 /* Split variables */
5821 for (size_t i = 0; i < splits->size; i++)
5822 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
5825 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5826 msave->rowtype, ' ');
5830 for (size_t i = 0; i < factors->size; i++)
5831 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
5834 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5835 msave->varname_, ' ');
5837 /* Continuous variables. */
5838 for (size_t x = 0; x < m->size2; x++)
5839 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
5840 casewriter_write (common->writer, c);
5844 gsl_matrix_free (m);
5845 gsl_vector_free (factors);
5846 gsl_vector_free (splits);
5849 static struct matrix_cmd *
5850 matrix_parse_mget (struct matrix_state *s)
5852 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5853 *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
5855 struct mget_command *mget = &cmd->mget;
5859 if (lex_match_id (s->lexer, "FILE"))
5861 lex_match (s->lexer, T_EQUALS);
5863 fh_unref (mget->file);
5864 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5868 else if (lex_match_id (s->lexer, "TYPE"))
5870 lex_match (s->lexer, T_EQUALS);
5871 while (lex_token (s->lexer) != T_SLASH
5872 && lex_token (s->lexer) != T_ENDCMD)
5874 const char *rowtype = match_rowtype (s->lexer);
5878 stringi_set_insert (&mget->rowtypes, rowtype);
5883 lex_error_expecting (s->lexer, "FILE", "TYPE");
5886 if (lex_token (s->lexer) == T_ENDCMD)
5889 if (!lex_force_match (s->lexer, T_SLASH))
5895 matrix_cmd_destroy (cmd);
5899 static const struct variable *
5900 get_a8_var (const struct dictionary *d, const char *name)
5902 const struct variable *v = dict_lookup_var (d, name);
5905 msg (SE, _("Matrix data file lacks %s variable."), name);
5908 if (var_get_width (v) != 8)
5910 msg (SE, _("%s variable in matrix data file must be "
5911 "8-byte string, but it has width %d."),
5912 name, var_get_width (v));
5919 is_rowtype (const union value *v, const char *rowtype)
5921 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
5922 ss_rtrim (&vs, ss_cstr (" "));
5923 return ss_equals_case (vs, ss_cstr (rowtype));
5927 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
5928 const struct dictionary *d,
5929 const struct variable *rowtype_var,
5930 struct matrix_state *s, size_t si, size_t fi,
5931 size_t cs, size_t cn)
5936 const union value *rowtype_ = case_data (rows[0], rowtype_var);
5937 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
5938 : is_rowtype (rowtype_, "CORR") ? "CR"
5939 : is_rowtype (rowtype_, "MEAN") ? "MN"
5940 : is_rowtype (rowtype_, "STDDEV") ? "SD"
5941 : is_rowtype (rowtype_, "N") ? "NC"
5944 struct string name = DS_EMPTY_INITIALIZER;
5945 ds_put_cstr (&name, name_prefix);
5947 ds_put_format (&name, "F%zu", fi);
5949 ds_put_format (&name, "S%zu", si);
5951 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
5953 mv = matrix_var_create (s, ds_ss (&name));
5956 msg (SW, _("Matrix data file contains variable with existing name %s."),
5961 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
5962 size_t n_missing = 0;
5963 for (size_t y = 0; y < n_rows; y++)
5965 for (size_t x = 0; x < cn; x++)
5967 struct variable *var = dict_get_var (d, cs + x);
5968 double value = case_num (rows[y], var);
5969 if (var_is_num_missing (var, value, MV_ANY))
5974 gsl_matrix_set (m, y, x, value);
5979 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
5980 "value, which was treated as zero.",
5981 "Matrix data file variable %s contains %zu missing "
5982 "values, which were treated as zero.", n_missing),
5983 ds_cstr (&name), n_missing);
5988 for (size_t y = 0; y < n_rows; y++)
5989 case_unref (rows[y]);
5993 var_changed (const struct ccase *ca, const struct ccase *cb,
5994 const struct variable *var)
5996 return !value_equal (case_data (ca, var), case_data (cb, var),
5997 var_get_width (var));
6001 vars_changed (const struct ccase *ca, const struct ccase *cb,
6002 const struct dictionary *d,
6003 size_t first_var, size_t n_vars)
6005 for (size_t i = 0; i < n_vars; i++)
6007 const struct variable *v = dict_get_var (d, first_var + i);
6008 if (var_changed (ca, cb, v))
6015 matrix_cmd_execute_mget (struct mget_command *mget)
6017 struct dictionary *d;
6018 struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6023 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6024 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6025 if (!rowtype_ || !varname_)
6028 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6030 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6033 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6035 msg (SE, _("Matrix data file contains no continuous variables."));
6039 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6041 const struct variable *v = dict_get_var (d, i);
6042 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6045 _("Matrix data file contains unexpected string variable %s."),
6051 /* SPLIT variables. */
6053 size_t sn = var_get_dict_index (rowtype_);
6054 struct ccase *sc = NULL;
6057 /* FACTOR variables. */
6058 size_t fs = var_get_dict_index (rowtype_) + 1;
6059 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6060 struct ccase *fc = NULL;
6063 /* Continuous variables. */
6064 size_t cs = var_get_dict_index (varname_) + 1;
6065 size_t cn = dict_get_var_cnt (d) - cs;
6066 struct ccase *cc = NULL;
6069 struct ccase **rows = NULL;
6070 size_t allocated_rows = 0;
6074 while ((c = casereader_read (r)) != NULL)
6076 bool sd = vars_changed (sc, c, d, ss, sn);
6077 bool fd = sd || vars_changed (fc, c, d, fs, fn);
6078 bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6093 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6094 mget->state, si, fi, cs, cn);
6100 if (n_rows >= allocated_rows)
6101 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6104 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6105 mget->state, si, fi, cs, cn);
6109 casereader_destroy (r);
6113 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6115 if (!lex_force_id (s->lexer))
6118 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6120 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6125 static struct matrix_cmd *
6126 matrix_parse_eigen (struct matrix_state *s)
6128 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6129 *cmd = (struct matrix_cmd) {
6131 .eigen = { .expr = NULL }
6134 struct eigen_command *eigen = &cmd->eigen;
6135 if (!lex_force_match (s->lexer, T_LPAREN))
6137 eigen->expr = matrix_parse_expr (s);
6139 || !lex_force_match (s->lexer, T_COMMA)
6140 || !matrix_parse_dst_var (s, &eigen->evec)
6141 || !lex_force_match (s->lexer, T_COMMA)
6142 || !matrix_parse_dst_var (s, &eigen->eval)
6143 || !lex_force_match (s->lexer, T_RPAREN))
6149 matrix_cmd_destroy (cmd);
6154 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6156 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6159 if (!is_symmetric (A))
6161 msg (SE, _("Argument of EIGEN must be symmetric."));
6162 gsl_matrix_free (A);
6166 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6167 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6168 gsl_vector v_eval = to_vector (eval);
6169 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6170 gsl_eigen_symmv (A, &v_eval, evec, w);
6171 gsl_eigen_symmv_free (w);
6173 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6175 gsl_matrix_free (A);
6177 gsl_matrix_free (eigen->eval->value);
6178 eigen->eval->value = eval;
6180 gsl_matrix_free (eigen->evec->value);
6181 eigen->evec->value = evec;
6184 static struct matrix_cmd *
6185 matrix_parse_setdiag (struct matrix_state *s)
6187 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6188 *cmd = (struct matrix_cmd) {
6189 .type = MCMD_SETDIAG,
6190 .setdiag = { .dst = NULL }
6193 struct setdiag_command *setdiag = &cmd->setdiag;
6194 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6197 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6200 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6205 if (!lex_force_match (s->lexer, T_COMMA))
6208 setdiag->expr = matrix_parse_expr (s);
6212 if (!lex_force_match (s->lexer, T_RPAREN))
6218 matrix_cmd_destroy (cmd);
6223 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6225 gsl_matrix *dst = setdiag->dst->value;
6228 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6229 setdiag->dst->name);
6233 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6237 size_t n = MIN (dst->size1, dst->size2);
6238 if (is_scalar (src))
6240 double d = to_scalar (src);
6241 for (size_t i = 0; i < n; i++)
6242 gsl_matrix_set (dst, i, i, d);
6244 else if (is_vector (src))
6246 gsl_vector v = to_vector (src);
6247 for (size_t i = 0; i < n && i < v.size; i++)
6248 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6251 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6252 "not a %zu×%zu matrix."),
6253 src->size1, src->size2);
6254 gsl_matrix_free (src);
6257 static struct matrix_cmd *
6258 matrix_parse_svd (struct matrix_state *s)
6260 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6261 *cmd = (struct matrix_cmd) {
6263 .svd = { .expr = NULL }
6266 struct svd_command *svd = &cmd->svd;
6267 if (!lex_force_match (s->lexer, T_LPAREN))
6269 svd->expr = matrix_parse_expr (s);
6271 || !lex_force_match (s->lexer, T_COMMA)
6272 || !matrix_parse_dst_var (s, &svd->u)
6273 || !lex_force_match (s->lexer, T_COMMA)
6274 || !matrix_parse_dst_var (s, &svd->s)
6275 || !lex_force_match (s->lexer, T_COMMA)
6276 || !matrix_parse_dst_var (s, &svd->v)
6277 || !lex_force_match (s->lexer, T_RPAREN))
6283 matrix_cmd_destroy (cmd);
6288 matrix_cmd_execute_svd (struct svd_command *svd)
6290 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6294 if (m->size1 >= m->size2)
6297 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6298 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6299 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6300 gsl_vector *work = gsl_vector_alloc (A->size2);
6301 gsl_linalg_SV_decomp (A, V, &Sv, work);
6302 gsl_vector_free (work);
6304 matrix_var_set (svd->u, A);
6305 matrix_var_set (svd->s, S);
6306 matrix_var_set (svd->v, V);
6310 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6311 gsl_matrix_transpose_memcpy (At, m);
6312 gsl_matrix_free (m);
6314 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6315 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6316 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6317 gsl_vector *work = gsl_vector_alloc (At->size2);
6318 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6319 gsl_vector_free (work);
6321 matrix_var_set (svd->v, At);
6322 matrix_var_set (svd->s, St);
6323 matrix_var_set (svd->u, Vt);
6328 matrix_cmd_execute (struct matrix_cmd *cmd)
6333 matrix_cmd_execute_compute (&cmd->compute);
6337 matrix_cmd_execute_print (&cmd->print);
6341 return matrix_cmd_execute_do_if (&cmd->do_if);
6344 matrix_cmd_execute_loop (&cmd->loop);
6351 matrix_cmd_execute_display (&cmd->display);
6355 matrix_cmd_execute_release (&cmd->release);
6359 matrix_cmd_execute_save (&cmd->save);
6363 matrix_cmd_execute_read (&cmd->read);
6367 matrix_cmd_execute_write (&cmd->write);
6371 matrix_cmd_execute_get (&cmd->get);
6375 matrix_cmd_execute_msave (&cmd->msave);
6379 matrix_cmd_execute_mget (&cmd->mget);
6383 matrix_cmd_execute_eigen (&cmd->eigen);
6387 matrix_cmd_execute_setdiag (&cmd->setdiag);
6391 matrix_cmd_execute_svd (&cmd->svd);
6399 matrix_cmds_uninit (struct matrix_cmds *cmds)
6401 for (size_t i = 0; i < cmds->n; i++)
6402 matrix_cmd_destroy (cmds->commands[i]);
6403 free (cmds->commands);
6407 matrix_cmd_destroy (struct matrix_cmd *cmd)
6415 matrix_lvalue_destroy (cmd->compute.lvalue);
6416 matrix_expr_destroy (cmd->compute.rvalue);
6420 matrix_expr_destroy (cmd->print.expression);
6421 free (cmd->print.title);
6422 print_labels_destroy (cmd->print.rlabels);
6423 print_labels_destroy (cmd->print.clabels);
6427 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6429 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6430 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6432 free (cmd->do_if.clauses);
6436 matrix_expr_destroy (cmd->loop.start);
6437 matrix_expr_destroy (cmd->loop.end);
6438 matrix_expr_destroy (cmd->loop.increment);
6439 matrix_expr_destroy (cmd->loop.top_condition);
6440 matrix_expr_destroy (cmd->loop.bottom_condition);
6441 matrix_cmds_uninit (&cmd->loop.commands);
6451 free (cmd->release.vars);
6455 matrix_expr_destroy (cmd->save.expression);
6456 fh_unref (cmd->save.outfile);
6457 string_array_destroy (cmd->save.variables);
6458 matrix_expr_destroy (cmd->save.names);
6459 stringi_set_destroy (&cmd->save.strings);
6463 matrix_lvalue_destroy (cmd->read.dst);
6464 matrix_expr_destroy (cmd->read.size);
6468 matrix_expr_destroy (cmd->write.expression);
6469 free (cmd->write.encoding);
6473 matrix_lvalue_destroy (cmd->get.dst);
6474 fh_unref (cmd->get.file);
6475 free (cmd->get.encoding);
6476 string_array_destroy (&cmd->get.variables);
6480 free (cmd->msave.varname_);
6481 matrix_expr_destroy (cmd->msave.expr);
6482 matrix_expr_destroy (cmd->msave.factors);
6483 matrix_expr_destroy (cmd->msave.splits);
6487 fh_unref (cmd->mget.file);
6488 stringi_set_destroy (&cmd->mget.rowtypes);
6492 matrix_expr_destroy (cmd->eigen.expr);
6496 matrix_expr_destroy (cmd->setdiag.expr);
6500 matrix_expr_destroy (cmd->svd.expr);
6506 struct matrix_command_name
6509 struct matrix_cmd *(*parse) (struct matrix_state *);
6512 static const struct matrix_command_name *
6513 matrix_parse_command_name (struct lexer *lexer)
6515 static const struct matrix_command_name commands[] = {
6516 { "COMPUTE", matrix_parse_compute },
6517 { "CALL EIGEN", matrix_parse_eigen },
6518 { "CALL SETDIAG", matrix_parse_setdiag },
6519 { "CALL SVD", matrix_parse_svd },
6520 { "PRINT", matrix_parse_print },
6521 { "DO IF", matrix_parse_do_if },
6522 { "LOOP", matrix_parse_loop },
6523 { "BREAK", matrix_parse_break },
6524 { "READ", matrix_parse_read },
6525 { "WRITE", matrix_parse_write },
6526 { "GET", matrix_parse_get },
6527 { "SAVE", matrix_parse_save },
6528 { "MGET", matrix_parse_mget },
6529 { "MSAVE", matrix_parse_msave },
6530 { "DISPLAY", matrix_parse_display },
6531 { "RELEASE", matrix_parse_release },
6533 static size_t n = sizeof commands / sizeof *commands;
6535 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6536 if (lex_match_phrase (lexer, c->name))
6541 static struct matrix_cmd *
6542 matrix_parse_command (struct matrix_state *s)
6544 size_t nesting_level = SIZE_MAX;
6546 struct matrix_cmd *c = NULL;
6547 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6549 lex_error (s->lexer, _("Unknown matrix command."));
6550 else if (!cmd->parse)
6551 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6555 nesting_level = output_open_group (
6556 group_item_create_nocopy (utf8_to_title (cmd->name),
6557 utf8_to_title (cmd->name)));
6562 lex_end_of_command (s->lexer);
6563 lex_discard_rest_of_command (s->lexer);
6564 if (nesting_level != SIZE_MAX)
6565 output_close_groups (nesting_level);
6571 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6573 if (!lex_force_match (lexer, T_ENDCMD))
6576 struct matrix_state state = {
6577 .session = dataset_session (ds),
6579 .vars = HMAP_INITIALIZER (state.vars),
6584 while (lex_match (lexer, T_ENDCMD))
6586 if (lex_token (lexer) == T_STOP)
6588 msg (SE, _("Unexpected end of input expecting matrix command."));
6592 if (lex_match_phrase (lexer, "END MATRIX"))
6595 struct matrix_cmd *c = matrix_parse_command (&state);
6598 matrix_cmd_execute (c);
6599 matrix_cmd_destroy (c);
6603 struct matrix_var *var, *next;
6604 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6607 gsl_matrix_free (var->value);
6608 hmap_delete (&state.vars, &var->hmap_node);
6611 hmap_destroy (&state.vars);
6612 fh_unref (state.prev_save_outfile);
6613 fh_unref (state.prev_write_outfile);
6616 dict_unref (state.common->dict);
6617 casewriter_destroy (state.common->writer);
6618 free (state.common);
6620 fh_unref (state.prev_read_file);
6621 for (size_t i = 0; i < state.n_read_files; i++)
6622 read_file_destroy (state.read_files[i]);
6623 free (state.read_files);