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;
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."),
4706 matrix_cmd_destroy (cmd);
4712 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4713 gsl_matrix *m, struct substring p, size_t y, size_t x)
4715 const char *input_encoding = dfm_reader_get_encoding (reader);
4717 char *error = data_in (p, input_encoding, read->format,
4718 settings_get_fmt_settings (), &v, 0, NULL);
4719 /* XXX report error if value is missing */
4721 msg (SW, _("GET parse error (%.*s): %s"), (int) p.length, p.string, error);
4724 gsl_matrix_set (m, y, x, v.f);
4725 if (read->symmetric && x != y)
4726 gsl_matrix_set (m, x, y, v.f);
4731 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4732 struct substring *line)
4734 if (dfm_eof (reader))
4736 msg (SE, _("Unexpected end of file reading matrix data."));
4739 dfm_expand_tabs (reader);
4740 *line = ss_substr (dfm_get_record (reader),
4741 read->c1 - 1, read->c2 - read->c1);
4746 matrix_read (struct read_command *read, struct dfm_reader *reader,
4749 for (size_t y = 0; y < m->size1; y++)
4751 size_t nx = read->symmetric ? y + 1 : m->size2;
4753 struct substring line = ss_empty ();
4754 for (size_t x = 0; x < nx; x++)
4761 ss_ltrim (&line, ss_cstr (" ,"));
4762 if (!ss_is_empty (line))
4764 if (!matrix_read_line (read, reader, &line))
4766 dfm_forward_record (reader);
4769 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
4773 if (!matrix_read_line (read, reader, &line))
4775 size_t fields_per_line = (read->c2 - read->c1) / read->w;
4776 int f = x % fields_per_line;
4777 if (f == fields_per_line - 1)
4778 dfm_forward_record (reader);
4780 p = ss_substr (line, read->w * f, read->w);
4783 matrix_read_set_field (read, reader, m, p, y, x);
4787 dfm_forward_record (reader);
4790 ss_ltrim (&line, ss_cstr (" ,"));
4791 if (!ss_is_empty (line))
4792 msg (SW, _("Trailing garbage on line \"%.*s\""),
4793 (int) line.length, line.string);
4799 matrix_cmd_execute_read (struct read_command *read)
4801 struct index_vector iv0, iv1;
4802 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
4805 size_t size[2] = { SIZE_MAX, SIZE_MAX };
4808 gsl_matrix *m = matrix_expr_evaluate (read->size);
4814 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
4815 gsl_matrix_free (m);
4821 gsl_vector v = to_vector (m);
4825 d[0] = gsl_vector_get (&v, 0);
4828 else if (v.size == 2)
4830 d[0] = gsl_vector_get (&v, 0);
4831 d[1] = gsl_vector_get (&v, 1);
4835 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector"));
4836 gsl_matrix_free (m);
4841 gsl_matrix_free (m);
4843 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
4845 msg (SE, _("SIZE (%g,%g) is outside valid range."), d[0], d[1]);
4855 if (read->dst->n_indexes)
4857 size_t submatrix_size[2];
4858 if (read->dst->n_indexes == 2)
4860 submatrix_size[0] = iv0.n;
4861 submatrix_size[1] = iv1.n;
4863 else if (read->dst->var->value->size1 == 1)
4865 submatrix_size[0] = 1;
4866 submatrix_size[1] = iv0.n;
4870 submatrix_size[0] = iv0.n;
4871 submatrix_size[1] = 1;
4876 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
4878 msg (SE, _("SIZE (%zu,%zu) differs from submatrix dimensions "
4881 submatrix_size[0], submatrix_size[1]);
4889 size[0] = submatrix_size[0];
4890 size[1] = submatrix_size[1];
4894 struct dfm_reader *reader = read_file_open (read->rf);
4896 dfm_reread_record (reader, 1);
4898 if (read->symmetric && size[0] != size[1])
4900 msg (SE, _("Cannot read non-square %zu×%zu matrix "
4901 "using READ with MODE=SYMMETRIC."),
4907 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
4908 matrix_read (read, reader, tmp);
4909 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
4912 static struct matrix_cmd *
4913 matrix_parse_write (struct matrix_state *s)
4915 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4916 *cmd = (struct matrix_cmd) {
4918 .write = { .format = FMT_F },
4921 struct write_command *write = &cmd->write;
4922 write->expression = matrix_parse_exp (s);
4923 if (!write->expression)
4927 int repetitions = 0;
4928 int record_width = 0;
4929 bool seen_format = false;
4930 while (lex_match (s->lexer, T_SLASH))
4932 if (lex_match_id (s->lexer, "OUTFILE"))
4934 lex_match (s->lexer, T_EQUALS);
4936 fh_unref (write->outfile);
4937 write->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
4938 if (!write->outfile)
4941 else if (lex_match_id (s->lexer, "ENCODING"))
4943 lex_match (s->lexer, T_EQUALS);
4944 if (!lex_force_string (s->lexer))
4947 free (write->encoding);
4948 write->encoding = ss_xstrdup (lex_tokss (s->lexer));
4952 else if (lex_match_id (s->lexer, "FIELD"))
4954 lex_match (s->lexer, T_EQUALS);
4956 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4958 write->c1 = lex_integer (s->lexer);
4960 if (!lex_force_match (s->lexer, T_TO)
4961 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
4963 write->c2 = lex_integer (s->lexer) + 1;
4966 record_width = write->c2 - write->c1;
4967 if (lex_match (s->lexer, T_BY))
4969 if (!lex_force_int_range (s->lexer, "BY", 1,
4970 write->c2 - write->c1))
4972 by = lex_integer (s->lexer);
4975 if (record_width % by)
4977 msg (SE, _("BY %d does not evenly divide record width %d."),
4985 else if (lex_match_id (s->lexer, "MODE"))
4987 lex_match (s->lexer, T_EQUALS);
4988 if (lex_match_id (s->lexer, "RECTANGULAR"))
4989 write->triangular = false;
4990 else if (lex_match_id (s->lexer, "TRIANGULAR"))
4991 write->triangular = true;
4994 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
4998 else if (lex_match_id (s->lexer, "HOLD"))
5000 else if (lex_match_id (s->lexer, "FORMAT"))
5004 lex_sbc_only_once ("FORMAT");
5009 lex_match (s->lexer, T_EQUALS);
5011 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5014 const char *p = lex_tokcstr (s->lexer);
5015 if (c_isdigit (p[0]))
5017 repetitions = atoi (p);
5018 p += strspn (p, "0123456789");
5019 if (!fmt_from_name (p, &write->format))
5021 lex_error (s->lexer, _("Unknown format %s."), p);
5026 else if (!fmt_from_name (p, &write->format))
5028 struct fmt_spec format;
5029 if (!parse_format_specifier (s->lexer, &format))
5031 write->format = format.type;
5032 write->w = format.w;
5033 write->d = format.d;
5038 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5046 lex_sbc_missing ("FIELD");
5050 if (!write->outfile)
5052 if (s->prev_write_outfile)
5053 write->outfile = fh_ref (s->prev_write_outfile);
5056 lex_sbc_missing ("OUTFILE");
5060 fh_unref (s->prev_write_outfile);
5061 s->prev_write_outfile = fh_ref (write->outfile);
5063 /* Field width may be specified in multiple ways:
5066 2. The format on FORMAT.
5067 3. The repetition factor on FORMAT.
5069 (2) and (3) are mutually exclusive.
5071 If more than one of these is present, they must agree. If none of them is
5072 present, then free-field format is used.
5074 if (repetitions > record_width)
5076 msg (SE, _("%d repetitions cannot fit in record width %d."),
5077 repetitions, record_width);
5080 int w = (repetitions ? record_width / repetitions
5081 : write->w ? write->w
5086 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5087 "which implies field width %d, "
5088 "but BY specifies field width %d."),
5089 repetitions, record_width, w, by);
5091 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5099 matrix_cmd_destroy (cmd);
5104 matrix_cmd_execute_write (struct write_command *write)
5106 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5110 if (write->triangular && m->size1 != m->size2)
5112 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5113 "the matrix to be written has dimensions %zu×%zu."),
5114 m->size1, m->size2);
5115 gsl_matrix_free (m);
5119 struct dfm_writer *writer = dfm_open_writer (write->outfile, write->encoding);
5122 gsl_matrix_free (m);
5126 const struct fmt_settings *settings = settings_get_fmt_settings ();
5127 struct fmt_spec format = {
5128 .type = write->format,
5129 .w = write->w ? write->w : 40,
5132 struct u8_line line = U8_LINE_INITIALIZER;
5133 for (size_t y = 0; y < m->size1; y++)
5135 size_t nx = write->triangular ? y + 1 : m->size2;
5137 for (size_t x = 0; x < nx; x++)
5139 /* XXX string values */
5140 union value v = { .f = gsl_matrix_get (m, y, x) };
5142 ? data_out (&v, NULL, &format, settings)
5143 : data_out_stretchy (&v, NULL, &format, settings, NULL));
5144 size_t len = strlen (s);
5145 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5146 if (width + x0 > write->c2)
5148 dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5149 u8_line_clear (&line);
5152 u8_line_put (&line, x0, x0 + width, s, len);
5155 x0 += write->w ? write->w : width + 1;
5158 dfm_put_record_utf8 (writer, line.s.ss.string, line.s.ss.length);
5159 u8_line_clear (&line);
5161 u8_line_destroy (&line);
5162 dfm_close_writer (writer);
5164 gsl_matrix_free (m);
5167 static struct matrix_cmd *
5168 matrix_parse_get (struct matrix_state *s)
5170 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5171 *cmd = (struct matrix_cmd) {
5174 .user = { .treatment = MGET_ERROR },
5175 .system = { .treatment = MGET_ERROR },
5179 struct get_command *get = &cmd->get;
5180 get->dst = matrix_lvalue_parse (s);
5184 while (lex_match (s->lexer, T_SLASH))
5186 if (lex_match_id (s->lexer, "FILE"))
5188 if (get->variables.n)
5190 lex_error (s->lexer, _("FILE must precede VARIABLES"));
5193 lex_match (s->lexer, T_EQUALS);
5195 fh_unref (get->file);
5196 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5200 else if (lex_match_id (s->lexer, "ENCODING"))
5202 if (get->variables.n)
5204 lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
5207 lex_match (s->lexer, T_EQUALS);
5208 if (!lex_force_string (s->lexer))
5211 free (get->encoding);
5212 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5216 else if (lex_match_id (s->lexer, "VARIABLES"))
5218 lex_match (s->lexer, T_EQUALS);
5220 struct dictionary *dict = NULL;
5223 dict = dataset_dict (s->dataset);
5224 if (dict_get_var_cnt (dict) == 0)
5226 lex_error (s->lexer, _("GET cannot read empty active file."));
5232 struct casereader *reader = any_reader_open_and_decode (
5233 get->file, get->encoding, &dict, NULL);
5236 casereader_destroy (reader);
5239 struct variable **vars;
5241 bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
5242 PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
5249 string_array_clear (&get->variables);
5250 for (size_t i = 0; i < n_vars; i++)
5251 string_array_append (&get->variables, var_get_name (vars[i]));
5255 else if (lex_match_id (s->lexer, "NAMES"))
5257 lex_match (s->lexer, T_EQUALS);
5258 if (!lex_force_id (s->lexer))
5261 struct substring name = lex_tokss (s->lexer);
5262 get->names = matrix_var_lookup (s, name);
5264 get->names = matrix_var_create (s, name);
5267 else if (lex_match_id (s->lexer, "MISSING"))
5269 lex_match (s->lexer, T_EQUALS);
5270 if (lex_match_id (s->lexer, "ACCEPT"))
5271 get->user.treatment = MGET_ACCEPT;
5272 else if (lex_match_id (s->lexer, "OMIT"))
5273 get->user.treatment = MGET_OMIT;
5274 else if (lex_is_number (s->lexer))
5276 get->user.treatment = MGET_RECODE;
5277 get->user.substitute = lex_number (s->lexer);
5282 lex_error (s->lexer, NULL);
5286 else if (lex_match_id (s->lexer, "SYSMIS"))
5288 lex_match (s->lexer, T_EQUALS);
5289 if (lex_match_id (s->lexer, "OMIT"))
5290 get->user.treatment = MGET_OMIT;
5291 else if (lex_is_number (s->lexer))
5293 get->user.treatment = MGET_RECODE;
5294 get->user.substitute = lex_number (s->lexer);
5299 lex_error (s->lexer, NULL);
5305 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5306 "MISSING", "SYSMIS");
5313 matrix_cmd_destroy (cmd);
5318 matrix_cmd_execute_get (struct get_command *get)
5320 assert (get->file); /* XXX */
5322 struct dictionary *dict;
5323 struct casereader *reader = any_reader_open_and_decode (
5324 get->file, get->encoding, &dict, NULL);
5328 const struct variable **vars = xnmalloc (
5329 get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
5333 if (get->variables.n)
5335 for (size_t i = 0; i < get->variables.n; i++)
5337 const char *name = get->variables.strings[i];
5338 const struct variable *var = dict_lookup_var (dict, name);
5341 msg (SE, _("GET: Data file does not contain variable %s."),
5347 if (!var_is_numeric (var))
5349 msg (SE, _("GET: Variable %s is not numeric."), name);
5354 vars[n_vars++] = var;
5359 for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
5361 const struct variable *var = dict_get_var (dict, i);
5362 if (!var_is_numeric (var))
5364 msg (SE, _("GET: Variable %s is not numeric."),
5365 var_get_name (var));
5370 vars[n_vars++] = var;
5375 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5376 long long int casenum = 1;
5378 for (struct ccase *c = casereader_read (reader); c;
5379 c = casereader_read (reader), casenum++)
5381 if (n_rows >= m->size1)
5383 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5384 for (size_t y = 0; y < n_rows; y++)
5385 for (size_t x = 0; x < n_vars; x++)
5386 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5387 gsl_matrix_free (m);
5392 for (size_t x = 0; x < n_vars; x++)
5394 const struct variable *var = vars[x];
5395 double d = case_num (c, var);
5398 if (get->system.treatment == MGET_RECODE)
5399 d = get->system.substitute;
5400 else if (get->system.treatment == MGET_OMIT)
5404 msg (SE, _("GET: Variable %s in case %lld "
5405 "is system-missing."),
5406 var_get_name (var), casenum);
5410 else if (var_is_num_missing (var, d, MV_USER))
5412 if (get->user.treatment == MGET_RECODE)
5413 d = get->user.substitute;
5414 else if (get->user.treatment == MGET_OMIT)
5416 else if (get->user.treatment != MGET_ACCEPT)
5418 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5420 var_get_name (var), casenum, d);
5424 gsl_matrix_set (m, n_rows, x, d);
5432 casereader_destroy (reader);
5436 matrix_lvalue_evaluate_and_assign (get->dst, m);
5439 gsl_matrix_free (m);
5445 match_rowtype (struct lexer *lexer)
5447 static const char *rowtypes[] = {
5448 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5450 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5452 for (size_t i = 0; i < n_rowtypes; i++)
5453 if (lex_match_id (lexer, rowtypes[i]))
5456 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5461 parse_var_names (struct lexer *lexer, struct string_array *sa)
5463 lex_match (lexer, T_EQUALS);
5465 string_array_clear (sa);
5467 struct dictionary *dict = dict_create (get_default_encoding ());
5468 dict_create_var_assert (dict, "ROWTYPE_", 8);
5469 dict_create_var_assert (dict, "VARNAME_", 8);
5472 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5473 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5478 string_array_clear (sa);
5479 sa->strings = names;
5480 sa->n = sa->allocated = n_names;
5486 msave_common_uninit (struct msave_common *common)
5490 fh_unref (common->outfile);
5491 string_array_destroy (&common->variables);
5492 string_array_destroy (&common->fnames);
5493 string_array_destroy (&common->snames);
5498 compare_variables (const char *keyword,
5499 const struct string_array *new,
5500 const struct string_array *old)
5506 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5507 "on the first MSAVE within MATRIX."), keyword);
5510 else if (!string_array_equal_case (old, new))
5512 msg (SE, _("%s must specify the same variables each time within "
5513 "a given MATRIX."), keyword);
5520 static struct matrix_cmd *
5521 matrix_parse_msave (struct matrix_state *s)
5523 struct msave_common common = { .outfile = NULL };
5524 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5525 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5527 struct msave_command *msave = &cmd->msave;
5528 if (lex_token (s->lexer) == T_ID)
5529 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5530 msave->expr = matrix_parse_exp (s);
5534 while (lex_match (s->lexer, T_SLASH))
5536 if (lex_match_id (s->lexer, "TYPE"))
5538 lex_match (s->lexer, T_EQUALS);
5540 msave->rowtype = match_rowtype (s->lexer);
5541 if (!msave->rowtype)
5544 else if (lex_match_id (s->lexer, "OUTFILE"))
5546 lex_match (s->lexer, T_EQUALS);
5548 fh_unref (common.outfile);
5549 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5550 if (!common.outfile)
5553 else if (lex_match_id (s->lexer, "VARIABLES"))
5555 if (!parse_var_names (s->lexer, &common.variables))
5558 else if (lex_match_id (s->lexer, "FNAMES"))
5560 if (!parse_var_names (s->lexer, &common.fnames))
5563 else if (lex_match_id (s->lexer, "SNAMES"))
5565 if (!parse_var_names (s->lexer, &common.snames))
5568 else if (lex_match_id (s->lexer, "SPLIT"))
5570 lex_match (s->lexer, T_EQUALS);
5572 matrix_expr_destroy (msave->splits);
5573 msave->splits = matrix_parse_exp (s);
5577 else if (lex_match_id (s->lexer, "FACTOR"))
5579 lex_match (s->lexer, T_EQUALS);
5581 matrix_expr_destroy (msave->factors);
5582 msave->factors = matrix_parse_exp (s);
5583 if (!msave->factors)
5588 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5589 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5593 if (!msave->rowtype)
5595 lex_sbc_missing ("TYPE");
5598 common.has_splits = msave->splits || common.snames.n;
5599 common.has_factors = msave->factors || common.fnames.n;
5601 struct msave_common *c = s->common ? s->common : &common;
5602 if (c->fnames.n && !msave->factors)
5604 msg (SE, _("FNAMES requires FACTOR."));
5607 if (c->snames.n && !msave->splits)
5609 msg (SE, _("SNAMES requires SPLIT."));
5612 if (c->has_factors && !common.has_factors)
5614 msg (SE, _("%s is required because it was present on the first "
5615 "MSAVE in this MATRIX command."), "FACTOR");
5618 if (c->has_splits && !common.has_splits)
5620 msg (SE, _("%s is required because it was present on the first "
5621 "MSAVE in this MATRIX command."), "SPLIT");
5627 if (!common.outfile)
5629 lex_sbc_missing ("OUTFILE");
5632 s->common = xmemdup (&common, sizeof common);
5638 bool same = common.outfile == s->common->outfile;
5639 fh_unref (common.outfile);
5642 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5643 "within a single MATRIX command."));
5647 if (!compare_variables ("VARIABLES",
5648 &common.variables, &s->common->variables)
5649 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5650 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5652 msave_common_uninit (&common);
5654 msave->common = s->common;
5655 if (!msave->varname_)
5656 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5660 msave_common_uninit (&common);
5661 matrix_cmd_destroy (cmd);
5666 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5668 gsl_matrix *m = matrix_expr_evaluate (e);
5674 msg (SE, _("%s expression must evaluate to vector, "
5675 "not a %zu×%zu matrix."),
5676 name, m->size1, m->size2);
5677 gsl_matrix_free (m);
5681 return matrix_to_vector (m);
5685 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5687 for (size_t i = 0; i < vars->n; i++)
5688 if (!dict_create_var (d, vars->strings[i], 0))
5689 return vars->strings[i];
5693 static struct dictionary *
5694 msave_create_dict (const struct msave_common *common)
5696 struct dictionary *dict = dict_create (get_default_encoding ());
5698 const char *dup_split = msave_add_vars (dict, &common->snames);
5701 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5705 dict_create_var_assert (dict, "ROWTYPE_", 8);
5707 const char *dup_factor = msave_add_vars (dict, &common->fnames);
5710 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5714 dict_create_var_assert (dict, "VARNAME_", 8);
5716 const char *dup_var = msave_add_vars (dict, &common->variables);
5719 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
5731 matrix_cmd_execute_msave (struct msave_command *msave)
5733 struct msave_common *common = msave->common;
5734 gsl_matrix *m = NULL;
5735 gsl_vector *factors = NULL;
5736 gsl_vector *splits = NULL;
5738 m = matrix_expr_evaluate (msave->expr);
5742 if (!common->variables.n)
5743 for (size_t i = 0; i < m->size2; i++)
5744 string_array_append_nocopy (&common->variables,
5745 xasprintf ("COL%zu", i + 1));
5747 if (m->size2 != common->variables.n)
5749 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
5750 m->size2, common->variables.n);
5756 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
5760 if (!common->fnames.n)
5761 for (size_t i = 0; i < factors->size; i++)
5762 string_array_append_nocopy (&common->fnames,
5763 xasprintf ("FAC%zu", i + 1));
5765 if (factors->size != common->fnames.n)
5767 msg (SE, _("There are %zu factor variables, "
5768 "but %zu split values were supplied."),
5769 common->fnames.n, factors->size);
5776 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
5780 if (!common->fnames.n)
5781 for (size_t i = 0; i < splits->size; i++)
5782 string_array_append_nocopy (&common->fnames,
5783 xasprintf ("SPL%zu", i + 1));
5785 if (splits->size != common->fnames.n)
5787 msg (SE, _("There are %zu split variables, "
5788 "but %zu split values were supplied."),
5789 common->fnames.n, splits->size);
5794 if (!common->writer)
5796 struct dictionary *dict = msave_create_dict (common);
5800 common->writer = any_writer_open (common->outfile, dict);
5801 if (!common->writer)
5807 common->dict = dict;
5810 for (size_t y = 0; y < m->size1; y++)
5812 struct ccase *c = case_create (dict_get_proto (common->dict));
5815 /* Split variables */
5817 for (size_t i = 0; i < splits->size; i++)
5818 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
5821 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5822 msave->rowtype, ' ');
5826 for (size_t i = 0; i < factors->size; i++)
5827 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
5830 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
5831 msave->varname_, ' ');
5833 /* Continuous variables. */
5834 for (size_t x = 0; x < m->size2; x++)
5835 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
5836 casewriter_write (common->writer, c);
5840 gsl_matrix_free (m);
5841 gsl_vector_free (factors);
5842 gsl_vector_free (splits);
5845 static struct matrix_cmd *
5846 matrix_parse_mget (struct matrix_state *s)
5848 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5849 *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
5851 struct mget_command *mget = &cmd->mget;
5855 if (lex_match_id (s->lexer, "FILE"))
5857 lex_match (s->lexer, T_EQUALS);
5859 fh_unref (mget->file);
5860 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5864 else if (lex_match_id (s->lexer, "TYPE"))
5866 lex_match (s->lexer, T_EQUALS);
5867 while (lex_token (s->lexer) != T_SLASH
5868 && lex_token (s->lexer) != T_ENDCMD)
5870 const char *rowtype = match_rowtype (s->lexer);
5874 stringi_set_insert (&mget->rowtypes, rowtype);
5879 lex_error_expecting (s->lexer, "FILE", "TYPE");
5882 if (lex_token (s->lexer) == T_ENDCMD)
5885 if (!lex_force_match (s->lexer, T_SLASH))
5891 matrix_cmd_destroy (cmd);
5895 static const struct variable *
5896 get_a8_var (const struct dictionary *d, const char *name)
5898 const struct variable *v = dict_lookup_var (d, name);
5901 msg (SE, _("Matrix data file lacks %s variable."), name);
5904 if (var_get_width (v) != 8)
5906 msg (SE, _("%s variable in matrix data file must be "
5907 "8-byte string, but it has width %d."),
5908 name, var_get_width (v));
5915 is_rowtype (const union value *v, const char *rowtype)
5917 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
5918 ss_rtrim (&vs, ss_cstr (" "));
5919 return ss_equals_case (vs, ss_cstr (rowtype));
5923 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
5924 const struct dictionary *d,
5925 const struct variable *rowtype_var,
5926 struct matrix_state *s, size_t si, size_t fi,
5927 size_t cs, size_t cn)
5932 const union value *rowtype_ = case_data (rows[0], rowtype_var);
5933 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
5934 : is_rowtype (rowtype_, "CORR") ? "CR"
5935 : is_rowtype (rowtype_, "MEAN") ? "MN"
5936 : is_rowtype (rowtype_, "STDDEV") ? "SD"
5937 : is_rowtype (rowtype_, "N") ? "NC"
5940 struct string name = DS_EMPTY_INITIALIZER;
5941 ds_put_cstr (&name, name_prefix);
5943 ds_put_format (&name, "F%zu", fi);
5945 ds_put_format (&name, "S%zu", si);
5947 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
5949 mv = matrix_var_create (s, ds_ss (&name));
5952 msg (SW, _("Matrix data file contains variable with existing name %s."),
5957 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
5958 size_t n_missing = 0;
5959 for (size_t y = 0; y < n_rows; y++)
5961 for (size_t x = 0; x < cn; x++)
5963 struct variable *var = dict_get_var (d, cs + x);
5964 double value = case_num (rows[y], var);
5965 if (var_is_num_missing (var, value, MV_ANY))
5970 gsl_matrix_set (m, y, x, value);
5975 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
5976 "value, which was treated as zero.",
5977 "Matrix data file variable %s contains %zu missing "
5978 "values, which were treated as zero.", n_missing),
5979 ds_cstr (&name), n_missing);
5984 for (size_t y = 0; y < n_rows; y++)
5985 case_unref (rows[y]);
5989 var_changed (const struct ccase *ca, const struct ccase *cb,
5990 const struct variable *var)
5992 return !value_equal (case_data (ca, var), case_data (cb, var),
5993 var_get_width (var));
5997 vars_changed (const struct ccase *ca, const struct ccase *cb,
5998 const struct dictionary *d,
5999 size_t first_var, size_t n_vars)
6001 for (size_t i = 0; i < n_vars; i++)
6003 const struct variable *v = dict_get_var (d, first_var + i);
6004 if (var_changed (ca, cb, v))
6011 matrix_cmd_execute_mget (struct mget_command *mget)
6013 struct dictionary *d;
6014 struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6019 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6020 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6021 if (!rowtype_ || !varname_)
6024 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6026 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6029 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6031 msg (SE, _("Matrix data file contains no continuous variables."));
6035 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6037 const struct variable *v = dict_get_var (d, i);
6038 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6041 _("Matrix data file contains unexpected string variable %s."),
6047 /* SPLIT variables. */
6049 size_t sn = var_get_dict_index (rowtype_);
6050 struct ccase *sc = NULL;
6053 /* FACTOR variables. */
6054 size_t fs = var_get_dict_index (rowtype_) + 1;
6055 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6056 struct ccase *fc = NULL;
6059 /* Continuous variables. */
6060 size_t cs = var_get_dict_index (varname_) + 1;
6061 size_t cn = dict_get_var_cnt (d) - cs;
6062 struct ccase *cc = NULL;
6065 struct ccase **rows = NULL;
6066 size_t allocated_rows = 0;
6070 while ((c = casereader_read (r)) != NULL)
6072 bool sd = vars_changed (sc, c, d, ss, sn);
6073 bool fd = sd || vars_changed (fc, c, d, fs, fn);
6074 bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6089 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6090 mget->state, si, fi, cs, cn);
6096 if (n_rows >= allocated_rows)
6097 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6100 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6101 mget->state, si, fi, cs, cn);
6105 casereader_destroy (r);
6109 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6111 if (!lex_force_id (s->lexer))
6114 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6116 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6121 static struct matrix_cmd *
6122 matrix_parse_eigen (struct matrix_state *s)
6124 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6125 *cmd = (struct matrix_cmd) {
6127 .eigen = { .expr = NULL }
6130 struct eigen_command *eigen = &cmd->eigen;
6131 if (!lex_force_match (s->lexer, T_LPAREN))
6133 eigen->expr = matrix_parse_expr (s);
6135 || !lex_force_match (s->lexer, T_COMMA)
6136 || !matrix_parse_dst_var (s, &eigen->evec)
6137 || !lex_force_match (s->lexer, T_COMMA)
6138 || !matrix_parse_dst_var (s, &eigen->eval)
6139 || !lex_force_match (s->lexer, T_RPAREN))
6145 matrix_cmd_destroy (cmd);
6150 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6152 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6155 if (!is_symmetric (A))
6157 msg (SE, _("Argument of EIGEN must be symmetric."));
6158 gsl_matrix_free (A);
6162 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6163 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6164 gsl_vector v_eval = to_vector (eval);
6165 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6166 gsl_eigen_symmv (A, &v_eval, evec, w);
6167 gsl_eigen_symmv_free (w);
6169 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6171 gsl_matrix_free (A);
6173 gsl_matrix_free (eigen->eval->value);
6174 eigen->eval->value = eval;
6176 gsl_matrix_free (eigen->evec->value);
6177 eigen->evec->value = evec;
6180 static struct matrix_cmd *
6181 matrix_parse_setdiag (struct matrix_state *s)
6183 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6184 *cmd = (struct matrix_cmd) {
6185 .type = MCMD_SETDIAG,
6186 .setdiag = { .dst = NULL }
6189 struct setdiag_command *setdiag = &cmd->setdiag;
6190 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6193 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6196 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6201 if (!lex_force_match (s->lexer, T_COMMA))
6204 setdiag->expr = matrix_parse_expr (s);
6208 if (!lex_force_match (s->lexer, T_RPAREN))
6214 matrix_cmd_destroy (cmd);
6219 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6221 gsl_matrix *dst = setdiag->dst->value;
6224 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6225 setdiag->dst->name);
6229 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6233 size_t n = MIN (dst->size1, dst->size2);
6234 if (is_scalar (src))
6236 double d = to_scalar (src);
6237 for (size_t i = 0; i < n; i++)
6238 gsl_matrix_set (dst, i, i, d);
6240 else if (is_vector (src))
6242 gsl_vector v = to_vector (src);
6243 for (size_t i = 0; i < n && i < v.size; i++)
6244 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6247 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6248 "not a %zu×%zu matrix."),
6249 src->size1, src->size2);
6250 gsl_matrix_free (src);
6253 static struct matrix_cmd *
6254 matrix_parse_svd (struct matrix_state *s)
6256 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6257 *cmd = (struct matrix_cmd) {
6259 .svd = { .expr = NULL }
6262 struct svd_command *svd = &cmd->svd;
6263 if (!lex_force_match (s->lexer, T_LPAREN))
6265 svd->expr = matrix_parse_expr (s);
6267 || !lex_force_match (s->lexer, T_COMMA)
6268 || !matrix_parse_dst_var (s, &svd->u)
6269 || !lex_force_match (s->lexer, T_COMMA)
6270 || !matrix_parse_dst_var (s, &svd->s)
6271 || !lex_force_match (s->lexer, T_COMMA)
6272 || !matrix_parse_dst_var (s, &svd->v)
6273 || !lex_force_match (s->lexer, T_RPAREN))
6279 matrix_cmd_destroy (cmd);
6284 matrix_cmd_execute_svd (struct svd_command *svd)
6286 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6290 if (m->size1 >= m->size2)
6293 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6294 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6295 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6296 gsl_vector *work = gsl_vector_alloc (A->size2);
6297 gsl_linalg_SV_decomp (A, V, &Sv, work);
6298 gsl_vector_free (work);
6300 matrix_var_set (svd->u, A);
6301 matrix_var_set (svd->s, S);
6302 matrix_var_set (svd->v, V);
6306 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6307 gsl_matrix_transpose_memcpy (At, m);
6308 gsl_matrix_free (m);
6310 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6311 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6312 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6313 gsl_vector *work = gsl_vector_alloc (At->size2);
6314 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6315 gsl_vector_free (work);
6317 matrix_var_set (svd->v, At);
6318 matrix_var_set (svd->s, St);
6319 matrix_var_set (svd->u, Vt);
6324 matrix_cmd_execute (struct matrix_cmd *cmd)
6329 matrix_cmd_execute_compute (&cmd->compute);
6333 matrix_cmd_execute_print (&cmd->print);
6337 return matrix_cmd_execute_do_if (&cmd->do_if);
6340 matrix_cmd_execute_loop (&cmd->loop);
6347 matrix_cmd_execute_display (&cmd->display);
6351 matrix_cmd_execute_release (&cmd->release);
6355 matrix_cmd_execute_save (&cmd->save);
6359 matrix_cmd_execute_read (&cmd->read);
6363 matrix_cmd_execute_write (&cmd->write);
6367 matrix_cmd_execute_get (&cmd->get);
6371 matrix_cmd_execute_msave (&cmd->msave);
6375 matrix_cmd_execute_mget (&cmd->mget);
6379 matrix_cmd_execute_eigen (&cmd->eigen);
6383 matrix_cmd_execute_setdiag (&cmd->setdiag);
6387 matrix_cmd_execute_svd (&cmd->svd);
6395 matrix_cmds_uninit (struct matrix_cmds *cmds)
6397 for (size_t i = 0; i < cmds->n; i++)
6398 matrix_cmd_destroy (cmds->commands[i]);
6399 free (cmds->commands);
6403 matrix_cmd_destroy (struct matrix_cmd *cmd)
6411 matrix_lvalue_destroy (cmd->compute.lvalue);
6412 matrix_expr_destroy (cmd->compute.rvalue);
6416 matrix_expr_destroy (cmd->print.expression);
6417 free (cmd->print.title);
6418 print_labels_destroy (cmd->print.rlabels);
6419 print_labels_destroy (cmd->print.clabels);
6423 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6425 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6426 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6428 free (cmd->do_if.clauses);
6432 matrix_expr_destroy (cmd->loop.start);
6433 matrix_expr_destroy (cmd->loop.end);
6434 matrix_expr_destroy (cmd->loop.increment);
6435 matrix_expr_destroy (cmd->loop.top_condition);
6436 matrix_expr_destroy (cmd->loop.bottom_condition);
6437 matrix_cmds_uninit (&cmd->loop.commands);
6447 free (cmd->release.vars);
6451 matrix_expr_destroy (cmd->save.expression);
6452 fh_unref (cmd->save.outfile);
6453 string_array_destroy (cmd->save.variables);
6454 matrix_expr_destroy (cmd->save.names);
6455 stringi_set_destroy (&cmd->save.strings);
6459 matrix_lvalue_destroy (cmd->read.dst);
6460 matrix_expr_destroy (cmd->read.size);
6464 matrix_expr_destroy (cmd->write.expression);
6465 free (cmd->write.encoding);
6469 matrix_lvalue_destroy (cmd->get.dst);
6470 fh_unref (cmd->get.file);
6471 free (cmd->get.encoding);
6472 string_array_destroy (&cmd->get.variables);
6476 free (cmd->msave.varname_);
6477 matrix_expr_destroy (cmd->msave.expr);
6478 matrix_expr_destroy (cmd->msave.factors);
6479 matrix_expr_destroy (cmd->msave.splits);
6483 fh_unref (cmd->mget.file);
6484 stringi_set_destroy (&cmd->mget.rowtypes);
6488 matrix_expr_destroy (cmd->eigen.expr);
6492 matrix_expr_destroy (cmd->setdiag.expr);
6496 matrix_expr_destroy (cmd->svd.expr);
6502 struct matrix_command_name
6505 struct matrix_cmd *(*parse) (struct matrix_state *);
6508 static const struct matrix_command_name *
6509 matrix_parse_command_name (struct lexer *lexer)
6511 static const struct matrix_command_name commands[] = {
6512 { "COMPUTE", matrix_parse_compute },
6513 { "CALL EIGEN", matrix_parse_eigen },
6514 { "CALL SETDIAG", matrix_parse_setdiag },
6515 { "CALL SVD", matrix_parse_svd },
6516 { "PRINT", matrix_parse_print },
6517 { "DO IF", matrix_parse_do_if },
6518 { "LOOP", matrix_parse_loop },
6519 { "BREAK", matrix_parse_break },
6520 { "READ", matrix_parse_read },
6521 { "WRITE", matrix_parse_write },
6522 { "GET", matrix_parse_get },
6523 { "SAVE", matrix_parse_save },
6524 { "MGET", matrix_parse_mget },
6525 { "MSAVE", matrix_parse_msave },
6526 { "DISPLAY", matrix_parse_display },
6527 { "RELEASE", matrix_parse_release },
6529 static size_t n = sizeof commands / sizeof *commands;
6531 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6532 if (lex_match_phrase (lexer, c->name))
6537 static struct matrix_cmd *
6538 matrix_parse_command (struct matrix_state *s)
6540 size_t nesting_level = SIZE_MAX;
6542 struct matrix_cmd *c = NULL;
6543 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6545 lex_error (s->lexer, _("Unknown matrix command."));
6546 else if (!cmd->parse)
6547 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6551 nesting_level = output_open_group (
6552 group_item_create_nocopy (utf8_to_title (cmd->name),
6553 utf8_to_title (cmd->name)));
6558 lex_end_of_command (s->lexer);
6559 lex_discard_rest_of_command (s->lexer);
6560 if (nesting_level != SIZE_MAX)
6561 output_close_groups (nesting_level);
6567 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6569 if (!lex_force_match (lexer, T_ENDCMD))
6572 struct matrix_state state = {
6573 .session = dataset_session (ds),
6575 .vars = HMAP_INITIALIZER (state.vars),
6580 while (lex_match (lexer, T_ENDCMD))
6582 if (lex_token (lexer) == T_STOP)
6584 msg (SE, _("Unexpected end of input expecting matrix command."));
6588 if (lex_match_phrase (lexer, "END MATRIX"))
6591 struct matrix_cmd *c = matrix_parse_command (&state);
6594 matrix_cmd_execute (c);
6595 matrix_cmd_destroy (c);
6599 struct matrix_var *var, *next;
6600 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6603 gsl_matrix_free (var->value);
6604 hmap_delete (&state.vars, &var->hmap_node);
6607 hmap_destroy (&state.vars);
6608 fh_unref (state.prev_save_outfile);
6609 fh_unref (state.prev_write_outfile);
6612 dict_unref (state.common->dict);
6613 casewriter_destroy (state.common->writer);
6614 free (state.common);
6616 fh_unref (state.prev_read_file);
6617 for (size_t i = 0; i < state.n_read_files; i++)
6618 read_file_destroy (state.read_files[i]);
6619 free (state.read_files);