1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2021 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include <gsl/gsl_blas.h>
20 #include <gsl/gsl_eigen.h>
21 #include <gsl/gsl_linalg.h>
22 #include <gsl/gsl_matrix.h>
23 #include <gsl/gsl_permutation.h>
24 #include <gsl/gsl_randist.h>
25 #include <gsl/gsl_vector.h>
30 #include "data/any-reader.h"
31 #include "data/any-writer.h"
32 #include "data/casereader.h"
33 #include "data/casewriter.h"
34 #include "data/data-in.h"
35 #include "data/data-out.h"
36 #include "data/dataset.h"
37 #include "data/dictionary.h"
38 #include "data/file-handle-def.h"
39 #include "language/command.h"
40 #include "language/data-io/data-reader.h"
41 #include "language/data-io/data-writer.h"
42 #include "language/data-io/file-handle.h"
43 #include "language/lexer/format-parser.h"
44 #include "language/lexer/lexer.h"
45 #include "language/lexer/variable-parser.h"
46 #include "libpspp/array.h"
47 #include "libpspp/assertion.h"
48 #include "libpspp/compiler.h"
49 #include "libpspp/hmap.h"
50 #include "libpspp/i18n.h"
51 #include "libpspp/misc.h"
52 #include "libpspp/str.h"
53 #include "libpspp/string-array.h"
54 #include "libpspp/stringi-set.h"
55 #include "libpspp/u8-line.h"
56 #include "math/random.h"
57 #include "output/driver.h"
58 #include "output/output-item.h"
59 #include "output/pivot-table.h"
61 #include "gl/c-ctype.h"
62 #include "gl/ftoastr.h"
63 #include "gl/minmax.h"
67 #define _(msgid) gettext (msgid)
68 #define N_(msgid) (msgid)
72 struct hmap_node hmap_node;
80 struct file_handle *outfile;
81 struct string_array variables;
82 struct string_array fnames;
83 struct string_array snames;
88 /* Execution state. */
89 struct dictionary *dict;
90 struct casewriter *writer;
95 struct file_handle *file;
96 struct dfm_reader *reader;
102 struct file_handle *file;
103 struct dfm_writer *writer;
105 struct u8_line *held;
110 struct file_handle *file;
111 struct dataset *dataset;
113 /* Parameters from parsing the first SAVE command for 'file'. */
114 struct string_array variables;
115 struct matrix_expr *names;
116 struct stringi_set strings;
118 /* Results from the first (only) attempt to open this save_file. */
120 struct casewriter *writer;
121 struct dictionary *dict;
126 struct dataset *dataset;
127 struct session *session;
131 struct msave_common *common;
133 struct file_handle *prev_read_file;
134 struct read_file **read_files;
137 struct file_handle *prev_write_file;
138 struct write_file **write_files;
139 size_t n_write_files;
141 struct file_handle *prev_save_file;
142 struct save_file **save_files;
146 static struct matrix_var *
147 matrix_var_lookup (struct matrix_state *s, struct substring name)
149 struct matrix_var *var;
151 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
152 utf8_hash_case_substring (name, 0), &s->vars)
153 if (!utf8_sscasecmp (ss_cstr (var->name), name))
159 static struct matrix_var *
160 matrix_var_create (struct matrix_state *s, struct substring name)
162 struct matrix_var *var = xmalloc (sizeof *var);
163 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
164 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
169 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
171 gsl_matrix_free (var->value);
175 #define MATRIX_FUNCTIONS \
231 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
232 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
233 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
234 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
235 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
236 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
237 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
238 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
239 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
240 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
241 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
243 typedef gsl_matrix *matrix_proto_m_d (double);
244 typedef gsl_matrix *matrix_proto_m_dd (double, double);
245 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
246 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
247 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
248 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
249 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
250 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
251 typedef double matrix_proto_d_m (gsl_matrix *);
252 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
253 typedef gsl_matrix *matrix_proto_IDENT (double, double);
255 #define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
259 /* Matrix expressions. */
266 #define F(NAME, PROTOTYPE) MOP_F_##NAME,
270 /* Elementwise and scalar arithmetic. */
271 MOP_NEGATE, /* unary - */
272 MOP_ADD_ELEMS, /* + */
273 MOP_SUB_ELEMS, /* - */
274 MOP_MUL_ELEMS, /* &* */
275 MOP_DIV_ELEMS, /* / and &/ */
276 MOP_EXP_ELEMS, /* &** */
278 MOP_SEQ_BY, /* a:b:c */
280 /* Matrix arithmetic. */
282 MOP_EXP_MAT, /* ** */
299 MOP_PASTE_HORZ, /* a, b, c, ... */
300 MOP_PASTE_VERT, /* a; b; c; ... */
304 MOP_VEC_INDEX, /* x(y) */
305 MOP_VEC_ALL, /* x(:) */
306 MOP_MAT_INDEX, /* x(y,z) */
307 MOP_ROW_INDEX, /* x(y,:) */
308 MOP_COL_INDEX, /* x(:,z) */
315 MOP_EOF, /* EOF('file') */
323 struct matrix_expr **subs;
328 struct matrix_var *variable;
329 struct read_file *eof;
334 matrix_expr_destroy (struct matrix_expr *e)
341 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
372 for (size_t i = 0; i < e->n_subs; i++)
373 matrix_expr_destroy (e->subs[i]);
385 static struct matrix_expr *
386 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
389 struct matrix_expr *e = xmalloc (sizeof *e);
390 *e = (struct matrix_expr) {
392 .subs = xmemdup (subs, n_subs * sizeof *subs),
398 static struct matrix_expr *
399 matrix_expr_create_0 (enum matrix_op op)
401 struct matrix_expr *sub;
402 return matrix_expr_create_subs (op, &sub, 0);
405 static struct matrix_expr *
406 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
408 return matrix_expr_create_subs (op, &sub, 1);
411 static struct matrix_expr *
412 matrix_expr_create_2 (enum matrix_op op,
413 struct matrix_expr *sub0, struct matrix_expr *sub1)
415 struct matrix_expr *subs[] = { sub0, sub1 };
416 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
419 static struct matrix_expr *
420 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
421 struct matrix_expr *sub1, struct matrix_expr *sub2)
423 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
424 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
427 static struct matrix_expr *
428 matrix_expr_create_number (double number)
430 struct matrix_expr *e = xmalloc (sizeof *e);
431 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
435 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
437 static struct matrix_expr *
438 matrix_parse_curly_comma (struct matrix_state *s)
440 struct matrix_expr *lhs = matrix_parse_or_xor (s);
444 while (lex_match (s->lexer, T_COMMA))
446 struct matrix_expr *rhs = matrix_parse_or_xor (s);
449 matrix_expr_destroy (lhs);
452 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
457 static struct matrix_expr *
458 matrix_parse_curly_semi (struct matrix_state *s)
460 if (lex_token (s->lexer) == T_RCURLY)
461 return matrix_expr_create_0 (MOP_EMPTY);
463 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
467 while (lex_match (s->lexer, T_SEMICOLON))
469 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
472 matrix_expr_destroy (lhs);
475 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
480 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
481 for (size_t Y = 0; Y < (M)->size1; Y++) \
482 for (size_t X = 0; X < (M)->size2; X++) \
483 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
486 is_vector (const gsl_matrix *m)
488 return m->size1 <= 1 || m->size2 <= 1;
492 to_vector (gsl_matrix *m)
494 return (m->size1 == 1
495 ? gsl_matrix_row (m, 0).vector
496 : gsl_matrix_column (m, 0).vector);
501 matrix_eval_ABS (gsl_matrix *m)
503 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
509 matrix_eval_ALL (gsl_matrix *m)
511 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
518 matrix_eval_ANY (gsl_matrix *m)
520 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
527 matrix_eval_ARSIN (gsl_matrix *m)
529 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
535 matrix_eval_ARTAN (gsl_matrix *m)
537 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
543 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
547 for (size_t i = 0; i < n; i++)
552 gsl_matrix *block = gsl_matrix_calloc (r, c);
554 for (size_t i = 0; i < n; i++)
556 for (size_t y = 0; y < m[i]->size1; y++)
557 for (size_t x = 0; x < m[i]->size2; x++)
558 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
566 matrix_eval_CHOL (gsl_matrix *m)
568 if (!gsl_linalg_cholesky_decomp1 (m))
570 for (size_t y = 0; y < m->size1; y++)
571 for (size_t x = y + 1; x < m->size2; x++)
572 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
574 for (size_t y = 0; y < m->size1; y++)
575 for (size_t x = 0; x < y; x++)
576 gsl_matrix_set (m, y, x, 0);
581 msg (SE, _("Input to CHOL function is not positive-definite."));
587 matrix_eval_col_extremum (gsl_matrix *m, bool min)
592 return gsl_matrix_alloc (1, 0);
594 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
595 for (size_t x = 0; x < m->size2; x++)
597 double ext = gsl_matrix_get (m, 0, x);
598 for (size_t y = 1; y < m->size1; y++)
600 double value = gsl_matrix_get (m, y, x);
601 if (min ? value < ext : value > ext)
604 gsl_matrix_set (cext, 0, x, ext);
610 matrix_eval_CMAX (gsl_matrix *m)
612 return matrix_eval_col_extremum (m, false);
616 matrix_eval_CMIN (gsl_matrix *m)
618 return matrix_eval_col_extremum (m, true);
622 matrix_eval_COS (gsl_matrix *m)
624 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
630 matrix_eval_col_sum (gsl_matrix *m, bool square)
635 return gsl_matrix_alloc (1, 0);
637 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
638 for (size_t x = 0; x < m->size2; x++)
641 for (size_t y = 0; y < m->size1; y++)
643 double d = gsl_matrix_get (m, y, x);
644 sum += square ? pow2 (d) : d;
646 gsl_matrix_set (result, 0, x, sum);
652 matrix_eval_CSSQ (gsl_matrix *m)
654 return matrix_eval_col_sum (m, true);
658 matrix_eval_CSUM (gsl_matrix *m)
660 return matrix_eval_col_sum (m, false);
664 compare_double_3way (const void *a_, const void *b_)
666 const double *a = a_;
667 const double *b = b_;
668 return *a < *b ? -1 : *a > *b;
672 matrix_eval_DESIGN (gsl_matrix *m)
674 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
675 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
676 gsl_matrix_transpose_memcpy (&m2, m);
678 for (size_t y = 0; y < m2.size1; y++)
679 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
681 size_t *n = xcalloc (m2.size1, sizeof *n);
683 for (size_t i = 0; i < m2.size1; i++)
685 double *row = tmp + m2.size2 * i;
686 for (size_t j = 0; j < m2.size2; )
689 for (k = j + 1; k < m2.size2; k++)
690 if (row[j] != row[k])
692 row[n[i]++] = row[j];
697 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
702 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
704 for (size_t i = 0; i < m->size2; i++)
709 const double *unique = tmp + m2.size2 * i;
710 for (size_t j = 0; j < n[i]; j++, x++)
712 double value = unique[j];
713 for (size_t y = 0; y < m->size1; y++)
714 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
725 matrix_eval_DET (gsl_matrix *m)
727 gsl_permutation *p = gsl_permutation_alloc (m->size1);
729 gsl_linalg_LU_decomp (m, p, &signum);
730 gsl_permutation_free (p);
731 return gsl_linalg_LU_det (m, signum);
735 matrix_eval_DIAG (gsl_matrix *m)
737 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
738 for (size_t i = 0; i < diag->size1; i++)
739 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
744 is_symmetric (const gsl_matrix *m)
746 if (m->size1 != m->size2)
749 for (size_t y = 0; y < m->size1; y++)
750 for (size_t x = 0; x < y; x++)
751 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
758 compare_double_desc (const void *a_, const void *b_)
760 const double *a = a_;
761 const double *b = b_;
762 return *a > *b ? -1 : *a < *b;
766 matrix_eval_EVAL (gsl_matrix *m)
768 if (!is_symmetric (m))
770 msg (SE, _("Argument of EVAL must be symmetric."));
774 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
775 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
776 gsl_vector v_eval = to_vector (eval);
777 gsl_eigen_symm (m, &v_eval, w);
778 gsl_eigen_symm_free (w);
780 assert (v_eval.stride == 1);
781 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
787 matrix_eval_EXP (gsl_matrix *m)
789 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
794 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
797 Charl Linssen <charl@itfromb.it>
801 matrix_eval_GINV (gsl_matrix *A)
806 gsl_matrix *tmp_mat = NULL;
809 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
810 tmp_mat = gsl_matrix_alloc (m, n);
811 gsl_matrix_transpose_memcpy (tmp_mat, A);
819 gsl_matrix *V = gsl_matrix_alloc (m, m);
820 gsl_vector *u = gsl_vector_alloc (m);
822 gsl_vector *tmp_vec = gsl_vector_alloc (m);
823 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
824 gsl_vector_free (tmp_vec);
827 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
828 gsl_matrix_set_zero (Sigma_pinv);
829 double cutoff = 1e-15 * gsl_vector_max (u);
831 for (size_t i = 0; i < m; ++i)
833 double x = gsl_vector_get (u, i);
834 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
837 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
838 gsl_matrix *U = gsl_matrix_calloc (n, n);
839 for (size_t i = 0; i < n; i++)
840 for (size_t j = 0; j < m; j++)
841 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
843 /* two dot products to obtain pseudoinverse */
844 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
845 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
850 A_pinv = gsl_matrix_alloc (n, m);
851 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
855 A_pinv = gsl_matrix_alloc (m, n);
856 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
859 gsl_matrix_free (tmp_mat);
860 gsl_matrix_free (tmp_mat2);
862 gsl_matrix_free (Sigma_pinv);
876 grade_compare_3way (const void *a_, const void *b_)
878 const struct grade *a = a_;
879 const struct grade *b = b_;
881 return (a->value < b->value ? -1
882 : a->value > b->value ? 1
890 matrix_eval_GRADE (gsl_matrix *m)
892 size_t n = m->size1 * m->size2;
893 struct grade *grades = xmalloc (n * sizeof *grades);
896 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
897 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
898 qsort (grades, n, sizeof *grades, grade_compare_3way);
900 for (size_t i = 0; i < n; i++)
901 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
909 dot (gsl_vector *a, gsl_vector *b)
912 for (size_t i = 0; i < a->size; i++)
913 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
918 norm2 (gsl_vector *v)
921 for (size_t i = 0; i < v->size; i++)
922 result += pow2 (gsl_vector_get (v, i));
929 return sqrt (norm2 (v));
933 matrix_eval_GSCH (gsl_matrix *v)
935 if (v->size2 < v->size1)
937 msg (SE, _("GSCH requires its argument to have at least as many columns "
938 "as rows, but it has dimensions %zu×%zu."),
942 if (!v->size1 || !v->size2)
945 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
947 for (size_t vx = 0; vx < v->size2; vx++)
949 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
950 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
952 gsl_vector_memcpy (&u_i, &v_i);
953 for (size_t j = 0; j < ux; j++)
955 gsl_vector u_j = gsl_matrix_column (u, j).vector;
956 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
957 for (size_t k = 0; k < u_i.size; k++)
958 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
959 - scale * gsl_vector_get (&u_j, k)));
962 double len = norm (&u_i);
965 gsl_vector_scale (&u_i, 1.0 / len);
966 if (++ux >= v->size1)
973 msg (SE, _("%zu×%zu argument to GSCH contains only "
974 "%zu linearly independent columns."),
975 v->size1, v->size2, ux);
985 matrix_eval_IDENT (double s1, double s2)
987 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
989 msg (SE, _("Arguments to IDENT must be non-negative integers."));
993 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
994 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1000 invert_matrix (gsl_matrix *x)
1002 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1004 gsl_linalg_LU_decomp (x, p, &signum);
1005 gsl_linalg_LU_invx (x, p);
1006 gsl_permutation_free (p);
1010 matrix_eval_INV (gsl_matrix *m)
1017 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1019 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1020 a->size2 * b->size2);
1022 for (size_t ar = 0; ar < a->size1; ar++)
1023 for (size_t br = 0; br < b->size1; br++, y++)
1026 for (size_t ac = 0; ac < a->size2; ac++)
1027 for (size_t bc = 0; bc < b->size2; bc++, x++)
1029 double av = gsl_matrix_get (a, ar, ac);
1030 double bv = gsl_matrix_get (b, br, bc);
1031 gsl_matrix_set (k, y, x, av * bv);
1038 matrix_eval_LG10 (gsl_matrix *m)
1040 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1046 matrix_eval_LN (gsl_matrix *m)
1048 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1054 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1056 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1059 for (size_t i = 1; i <= n * n; i++)
1061 gsl_matrix_set (m, y, x, i);
1063 size_t y1 = !y ? n - 1 : y - 1;
1064 size_t x1 = x + 1 >= n ? 0 : x + 1;
1065 if (gsl_matrix_get (m, y1, x1) == 0)
1071 y = y + 1 >= n ? 0 : y + 1;
1076 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1078 double a = gsl_matrix_get (m, y1, x1);
1079 double b = gsl_matrix_get (m, y2, x2);
1080 gsl_matrix_set (m, y1, x1, b);
1081 gsl_matrix_set (m, y2, x2, a);
1085 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1089 /* A. Umar, "On the Construction of Even Order Magic Squares",
1090 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1092 for (size_t i = 1; i <= n * n / 2; i++)
1094 gsl_matrix_set (m, y, x, i);
1104 for (size_t i = n * n; i > n * n / 2; i--)
1106 gsl_matrix_set (m, y, x, i);
1114 for (size_t y = 0; y < n; y++)
1115 for (size_t x = 0; x < n / 2; x++)
1117 unsigned int d = gsl_matrix_get (m, y, x);
1118 if (d % 2 != (y < n / 2))
1119 magic_exchange (m, y, x, y, n - x - 1);
1124 size_t x1 = n / 2 - 1;
1126 magic_exchange (m, y1, x1, y2, x1);
1127 magic_exchange (m, y1, x2, y2, x2);
1131 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1133 /* A. Umar, "On the Construction of Even Order Magic Squares",
1134 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1138 for (size_t i = 1; ; i++)
1140 gsl_matrix_set (m, y, x, i);
1141 if (++y == n / 2 - 1)
1153 for (size_t i = n * n; ; i--)
1155 gsl_matrix_set (m, y, x, i);
1156 if (++y == n / 2 - 1)
1165 for (size_t y = 0; y < n; y++)
1166 if (y != n / 2 - 1 && y != n / 2)
1167 for (size_t x = 0; x < n / 2; x++)
1169 unsigned int d = gsl_matrix_get (m, y, x);
1170 if (d % 2 != (y < n / 2))
1171 magic_exchange (m, y, x, y, n - x - 1);
1174 size_t a0 = (n * n - 2 * n) / 2 + 1;
1175 for (size_t i = 0; i < n / 2; i++)
1178 gsl_matrix_set (m, n / 2, i, a);
1179 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1181 for (size_t i = 0; i < n / 2; i++)
1183 size_t a = a0 + i + n / 2;
1184 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1185 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1187 for (size_t x = 1; x < n / 2; x += 2)
1188 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1189 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1190 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1191 size_t x1 = n / 2 - 2;
1192 size_t x2 = n / 2 + 1;
1193 size_t y1 = n / 2 - 2;
1194 size_t y2 = n / 2 + 1;
1195 magic_exchange (m, y1, x1, y2, x1);
1196 magic_exchange (m, y1, x2, y2, x2);
1200 matrix_eval_MAGIC (double n_)
1202 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1204 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1209 gsl_matrix *m = gsl_matrix_calloc (n, n);
1211 matrix_eval_MAGIC_odd (m, n);
1213 matrix_eval_MAGIC_singly_even (m, n);
1215 matrix_eval_MAGIC_doubly_even (m, n);
1220 matrix_eval_MAKE (double r, double c, double value)
1222 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1224 msg (SE, _("Size arguments to MAKE must be integers."));
1228 gsl_matrix *m = gsl_matrix_alloc (r, c);
1229 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1235 matrix_eval_MDIAG (gsl_vector *v)
1237 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1238 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1239 gsl_vector_memcpy (&diagonal, v);
1244 matrix_eval_MMAX (gsl_matrix *m)
1246 return gsl_matrix_max (m);
1250 matrix_eval_MMIN (gsl_matrix *m)
1252 return gsl_matrix_min (m);
1256 matrix_eval_MOD (gsl_matrix *m, double divisor)
1260 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1264 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1265 *d = fmod (*d, divisor);
1270 matrix_eval_MSSQ (gsl_matrix *m)
1273 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1279 matrix_eval_MSUM (gsl_matrix *m)
1282 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1288 matrix_eval_NCOL (gsl_matrix *m)
1294 matrix_eval_NROW (gsl_matrix *m)
1300 matrix_eval_RANK (gsl_matrix *m)
1302 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1303 gsl_linalg_QR_decomp (m, tau);
1304 gsl_vector_free (tau);
1306 return gsl_linalg_QRPT_rank (m, -1);
1310 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1312 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1314 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1319 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1321 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1322 "product of matrix dimensions."));
1326 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1329 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1331 gsl_matrix_set (dst, y1, x1, *d);
1342 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1347 return gsl_matrix_alloc (0, 1);
1349 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1350 for (size_t y = 0; y < m->size1; y++)
1352 double ext = gsl_matrix_get (m, y, 0);
1353 for (size_t x = 1; x < m->size2; x++)
1355 double value = gsl_matrix_get (m, y, x);
1356 if (min ? value < ext : value > ext)
1359 gsl_matrix_set (rext, y, 0, ext);
1365 matrix_eval_RMAX (gsl_matrix *m)
1367 return matrix_eval_row_extremum (m, false);
1371 matrix_eval_RMIN (gsl_matrix *m)
1373 return matrix_eval_row_extremum (m, true);
1377 matrix_eval_RND (gsl_matrix *m)
1379 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1391 rank_compare_3way (const void *a_, const void *b_)
1393 const struct rank *a = a_;
1394 const struct rank *b = b_;
1396 return a->value < b->value ? -1 : a->value > b->value;
1400 matrix_eval_RNKORDER (gsl_matrix *m)
1402 size_t n = m->size1 * m->size2;
1403 struct rank *ranks = xmalloc (n * sizeof *ranks);
1405 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1406 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1407 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1409 for (size_t i = 0; i < n; )
1412 for (j = i + 1; j < n; j++)
1413 if (ranks[i].value != ranks[j].value)
1416 double rank = (i + j + 1.0) / 2.0;
1417 for (size_t k = i; k < j; k++)
1418 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1429 matrix_eval_row_sum (gsl_matrix *m, bool square)
1434 return gsl_matrix_alloc (0, 1);
1436 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1437 for (size_t y = 0; y < m->size1; y++)
1440 for (size_t x = 0; x < m->size2; x++)
1442 double d = gsl_matrix_get (m, y, x);
1443 sum += square ? pow2 (d) : d;
1445 gsl_matrix_set (result, y, 0, sum);
1451 matrix_eval_RSSQ (gsl_matrix *m)
1453 return matrix_eval_row_sum (m, true);
1457 matrix_eval_RSUM (gsl_matrix *m)
1459 return matrix_eval_row_sum (m, false);
1463 matrix_eval_SIN (gsl_matrix *m)
1465 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1471 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1473 if (m1->size1 != m2->size1)
1475 msg (SE, _("SOLVE requires its arguments to have the same number of "
1476 "rows, but the first argument has dimensions %zu×%zu and "
1477 "the second %zu×%zu."),
1478 m1->size1, m1->size2,
1479 m2->size1, m2->size2);
1483 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1484 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1486 gsl_linalg_LU_decomp (m1, p, &signum);
1487 for (size_t i = 0; i < m2->size2; i++)
1489 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1490 gsl_vector xi = gsl_matrix_column (x, i).vector;
1491 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1493 gsl_permutation_free (p);
1498 matrix_eval_SQRT (gsl_matrix *m)
1500 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1504 msg (SE, _("Argument to SQRT must be nonnegative."));
1513 matrix_eval_SSCP (gsl_matrix *m)
1515 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1516 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1521 matrix_eval_SVAL (gsl_matrix *m)
1523 gsl_matrix *tmp_mat = NULL;
1524 if (m->size2 > m->size1)
1526 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1527 gsl_matrix_transpose_memcpy (tmp_mat, m);
1532 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1533 gsl_vector *S = gsl_vector_alloc (m->size2);
1534 gsl_vector *work = gsl_vector_alloc (m->size2);
1535 gsl_linalg_SV_decomp (m, V, S, work);
1537 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1538 for (size_t i = 0; i < m->size2; i++)
1539 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1541 gsl_matrix_free (V);
1542 gsl_vector_free (S);
1543 gsl_vector_free (work);
1544 gsl_matrix_free (tmp_mat);
1550 matrix_eval_SWEEP (gsl_matrix *m, double d)
1552 if (d < 1 || d > SIZE_MAX)
1554 msg (SE, _("Scalar argument to SWEEP must be integer."));
1558 if (k >= MIN (m->size1, m->size2))
1560 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1561 "equal to the smaller of the matrix argument's rows and "
1566 double m_kk = gsl_matrix_get (m, k, k);
1567 if (fabs (m_kk) > 1e-19)
1569 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1570 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1572 double m_ij = gsl_matrix_get (m, i, j);
1573 double m_ik = gsl_matrix_get (m, i, k);
1574 double m_kj = gsl_matrix_get (m, k, j);
1575 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1584 for (size_t i = 0; i < m->size1; i++)
1586 gsl_matrix_set (m, i, k, 0);
1587 gsl_matrix_set (m, k, i, 0);
1594 matrix_eval_TRACE (gsl_matrix *m)
1597 size_t n = MIN (m->size1, m->size2);
1598 for (size_t i = 0; i < n; i++)
1599 sum += gsl_matrix_get (m, i, i);
1604 matrix_eval_T (gsl_matrix *m)
1606 return matrix_eval_TRANSPOS (m);
1610 matrix_eval_TRANSPOS (gsl_matrix *m)
1612 if (m->size1 == m->size2)
1614 gsl_matrix_transpose (m);
1619 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1620 gsl_matrix_transpose_memcpy (t, m);
1626 matrix_eval_TRUNC (gsl_matrix *m)
1628 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1634 matrix_eval_UNIFORM (double r_, double c_)
1636 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1638 msg (SE, _("Arguments to UNIFORM must be integers."));
1643 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1645 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1649 gsl_matrix *m = gsl_matrix_alloc (r, c);
1650 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1651 *d = gsl_ran_flat (get_rng (), 0, 1);
1655 struct matrix_function
1659 size_t min_args, max_args;
1662 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1664 static const struct matrix_function *
1665 matrix_parse_function_name (const char *token)
1667 static const struct matrix_function functions[] =
1669 #define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1673 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1675 for (size_t i = 0; i < N_FUNCTIONS; i++)
1677 if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1678 return &functions[i];
1683 static struct read_file *
1684 read_file_create (struct matrix_state *s, struct file_handle *fh)
1686 for (size_t i = 0; i < s->n_read_files; i++)
1688 struct read_file *rf = s->read_files[i];
1696 struct read_file *rf = xmalloc (sizeof *rf);
1697 *rf = (struct read_file) { .file = fh };
1699 s->read_files = xrealloc (s->read_files,
1700 (s->n_read_files + 1) * sizeof *s->read_files);
1701 s->read_files[s->n_read_files++] = rf;
1705 static struct dfm_reader *
1706 read_file_open (struct read_file *rf)
1709 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1714 read_file_destroy (struct read_file *rf)
1718 fh_unref (rf->file);
1719 dfm_close_reader (rf->reader);
1720 free (rf->encoding);
1725 static struct write_file *
1726 write_file_create (struct matrix_state *s, struct file_handle *fh)
1728 for (size_t i = 0; i < s->n_write_files; i++)
1730 struct write_file *wf = s->write_files[i];
1738 struct write_file *wf = xmalloc (sizeof *wf);
1739 *wf = (struct write_file) { .file = fh };
1741 s->write_files = xrealloc (s->write_files,
1742 (s->n_write_files + 1) * sizeof *s->write_files);
1743 s->write_files[s->n_write_files++] = wf;
1747 static struct dfm_writer *
1748 write_file_open (struct write_file *wf)
1751 wf->writer = dfm_open_writer (wf->file, wf->encoding);
1756 write_file_destroy (struct write_file *wf)
1762 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
1763 wf->held->s.ss.length);
1764 u8_line_destroy (wf->held);
1768 fh_unref (wf->file);
1769 dfm_close_writer (wf->writer);
1770 free (wf->encoding);
1776 matrix_parse_function (struct matrix_state *s, const char *token,
1777 struct matrix_expr **exprp)
1780 if (lex_next_token (s->lexer, 1) != T_LPAREN)
1783 if (lex_match_id (s->lexer, "EOF"))
1786 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1790 if (!lex_force_match (s->lexer, T_RPAREN))
1796 struct read_file *rf = read_file_create (s, fh);
1798 struct matrix_expr *e = xmalloc (sizeof *e);
1799 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1804 const struct matrix_function *f = matrix_parse_function_name (token);
1808 lex_get_n (s->lexer, 2);
1810 struct matrix_expr *e = xmalloc (sizeof *e);
1811 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1813 size_t allocated_subs = 0;
1816 struct matrix_expr *sub = matrix_parse_expr (s);
1820 if (e->n_subs >= allocated_subs)
1821 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1822 e->subs[e->n_subs++] = sub;
1824 while (lex_match (s->lexer, T_COMMA));
1825 if (!lex_force_match (s->lexer, T_RPAREN))
1832 matrix_expr_destroy (e);
1836 static struct matrix_expr *
1837 matrix_parse_primary (struct matrix_state *s)
1839 if (lex_is_number (s->lexer))
1841 double number = lex_number (s->lexer);
1843 return matrix_expr_create_number (number);
1845 else if (lex_is_string (s->lexer))
1847 char string[sizeof (double)];
1848 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1852 memcpy (&number, string, sizeof number);
1853 return matrix_expr_create_number (number);
1855 else if (lex_match (s->lexer, T_LPAREN))
1857 struct matrix_expr *e = matrix_parse_or_xor (s);
1858 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1860 matrix_expr_destroy (e);
1865 else if (lex_match (s->lexer, T_LCURLY))
1867 struct matrix_expr *e = matrix_parse_curly_semi (s);
1868 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1870 matrix_expr_destroy (e);
1875 else if (lex_token (s->lexer) == T_ID)
1877 struct matrix_expr *retval;
1878 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1881 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1884 lex_error (s->lexer, _("Unknown variable %s."),
1885 lex_tokcstr (s->lexer));
1890 struct matrix_expr *e = xmalloc (sizeof *e);
1891 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1894 else if (lex_token (s->lexer) == T_ALL)
1896 struct matrix_expr *retval;
1897 if (matrix_parse_function (s, "ALL", &retval))
1901 lex_error (s->lexer, NULL);
1905 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1908 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1910 if (lex_match (s->lexer, T_COLON))
1917 *indexp = matrix_parse_expr (s);
1918 return *indexp != NULL;
1922 static struct matrix_expr *
1923 matrix_parse_postfix (struct matrix_state *s)
1925 struct matrix_expr *lhs = matrix_parse_primary (s);
1926 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1929 struct matrix_expr *i0;
1930 if (!matrix_parse_index_expr (s, &i0))
1932 matrix_expr_destroy (lhs);
1935 if (lex_match (s->lexer, T_RPAREN))
1937 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1938 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1939 else if (lex_match (s->lexer, T_COMMA))
1941 struct matrix_expr *i1;
1942 if (!matrix_parse_index_expr (s, &i1)
1943 || !lex_force_match (s->lexer, T_RPAREN))
1945 matrix_expr_destroy (lhs);
1946 matrix_expr_destroy (i0);
1947 matrix_expr_destroy (i1);
1950 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1951 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1952 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1957 lex_error_expecting (s->lexer, "`)'", "`,'");
1962 static struct matrix_expr *
1963 matrix_parse_unary (struct matrix_state *s)
1965 if (lex_match (s->lexer, T_DASH))
1967 struct matrix_expr *lhs = matrix_parse_unary (s);
1968 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1970 else if (lex_match (s->lexer, T_PLUS))
1971 return matrix_parse_unary (s);
1973 return matrix_parse_postfix (s);
1976 static struct matrix_expr *
1977 matrix_parse_seq (struct matrix_state *s)
1979 struct matrix_expr *start = matrix_parse_unary (s);
1980 if (!start || !lex_match (s->lexer, T_COLON))
1983 struct matrix_expr *end = matrix_parse_unary (s);
1986 matrix_expr_destroy (start);
1990 if (lex_match (s->lexer, T_COLON))
1992 struct matrix_expr *increment = matrix_parse_unary (s);
1995 matrix_expr_destroy (start);
1996 matrix_expr_destroy (end);
1999 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2002 return matrix_expr_create_2 (MOP_SEQ, start, end);
2005 static struct matrix_expr *
2006 matrix_parse_exp (struct matrix_state *s)
2008 struct matrix_expr *lhs = matrix_parse_seq (s);
2015 if (lex_match (s->lexer, T_EXP))
2017 else if (lex_match_phrase (s->lexer, "&**"))
2022 struct matrix_expr *rhs = matrix_parse_seq (s);
2025 matrix_expr_destroy (lhs);
2028 lhs = matrix_expr_create_2 (op, lhs, rhs);
2032 static struct matrix_expr *
2033 matrix_parse_mul_div (struct matrix_state *s)
2035 struct matrix_expr *lhs = matrix_parse_exp (s);
2042 if (lex_match (s->lexer, T_ASTERISK))
2044 else if (lex_match (s->lexer, T_SLASH))
2046 else if (lex_match_phrase (s->lexer, "&*"))
2048 else if (lex_match_phrase (s->lexer, "&/"))
2053 struct matrix_expr *rhs = matrix_parse_exp (s);
2056 matrix_expr_destroy (lhs);
2059 lhs = matrix_expr_create_2 (op, lhs, rhs);
2063 static struct matrix_expr *
2064 matrix_parse_add_sub (struct matrix_state *s)
2066 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2073 if (lex_match (s->lexer, T_PLUS))
2075 else if (lex_match (s->lexer, T_DASH))
2077 else if (lex_token (s->lexer) == T_NEG_NUM)
2082 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2085 matrix_expr_destroy (lhs);
2088 lhs = matrix_expr_create_2 (op, lhs, rhs);
2092 static struct matrix_expr *
2093 matrix_parse_relational (struct matrix_state *s)
2095 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2102 if (lex_match (s->lexer, T_GT))
2104 else if (lex_match (s->lexer, T_GE))
2106 else if (lex_match (s->lexer, T_LT))
2108 else if (lex_match (s->lexer, T_LE))
2110 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2112 else if (lex_match (s->lexer, T_NE))
2117 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2120 matrix_expr_destroy (lhs);
2123 lhs = matrix_expr_create_2 (op, lhs, rhs);
2127 static struct matrix_expr *
2128 matrix_parse_not (struct matrix_state *s)
2130 if (lex_match (s->lexer, T_NOT))
2132 struct matrix_expr *lhs = matrix_parse_not (s);
2133 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2136 return matrix_parse_relational (s);
2139 static struct matrix_expr *
2140 matrix_parse_and (struct matrix_state *s)
2142 struct matrix_expr *lhs = matrix_parse_not (s);
2146 while (lex_match (s->lexer, T_AND))
2148 struct matrix_expr *rhs = matrix_parse_not (s);
2151 matrix_expr_destroy (lhs);
2154 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2159 static struct matrix_expr *
2160 matrix_parse_or_xor (struct matrix_state *s)
2162 struct matrix_expr *lhs = matrix_parse_and (s);
2169 if (lex_match (s->lexer, T_OR))
2171 else if (lex_match_id (s->lexer, "XOR"))
2176 struct matrix_expr *rhs = matrix_parse_and (s);
2179 matrix_expr_destroy (lhs);
2182 lhs = matrix_expr_create_2 (op, lhs, rhs);
2186 static struct matrix_expr *
2187 matrix_parse_expr (struct matrix_state *s)
2189 return matrix_parse_or_xor (s);
2192 /* Expression evaluation. */
2195 matrix_op_eval (enum matrix_op op, double a, double b)
2199 case MOP_ADD_ELEMS: return a + b;
2200 case MOP_SUB_ELEMS: return a - b;
2201 case MOP_MUL_ELEMS: return a * b;
2202 case MOP_DIV_ELEMS: return a / b;
2203 case MOP_EXP_ELEMS: return pow (a, b);
2204 case MOP_GT: return a > b;
2205 case MOP_GE: return a >= b;
2206 case MOP_LT: return a < b;
2207 case MOP_LE: return a <= b;
2208 case MOP_EQ: return a == b;
2209 case MOP_NE: return a != b;
2210 case MOP_AND: return (a > 0) && (b > 0);
2211 case MOP_OR: return (a > 0) || (b > 0);
2212 case MOP_XOR: return (a > 0) != (b > 0);
2214 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2223 case MOP_PASTE_HORZ:
2224 case MOP_PASTE_VERT:
2240 matrix_op_name (enum matrix_op op)
2244 case MOP_ADD_ELEMS: return "+";
2245 case MOP_SUB_ELEMS: return "-";
2246 case MOP_MUL_ELEMS: return "&*";
2247 case MOP_DIV_ELEMS: return "&/";
2248 case MOP_EXP_ELEMS: return "&**";
2249 case MOP_GT: return ">";
2250 case MOP_GE: return ">=";
2251 case MOP_LT: return "<";
2252 case MOP_LE: return "<=";
2253 case MOP_EQ: return "=";
2254 case MOP_NE: return "<>";
2255 case MOP_AND: return "AND";
2256 case MOP_OR: return "OR";
2257 case MOP_XOR: return "XOR";
2259 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2268 case MOP_PASTE_HORZ:
2269 case MOP_PASTE_VERT:
2285 is_scalar (const gsl_matrix *m)
2287 return m->size1 == 1 && m->size2 == 1;
2291 to_scalar (const gsl_matrix *m)
2293 assert (is_scalar (m));
2294 return gsl_matrix_get (m, 0, 0);
2298 matrix_expr_evaluate_elementwise (enum matrix_op op,
2299 gsl_matrix *a, gsl_matrix *b)
2303 double be = to_scalar (b);
2304 for (size_t r = 0; r < a->size1; r++)
2305 for (size_t c = 0; c < a->size2; c++)
2307 double *ae = gsl_matrix_ptr (a, r, c);
2308 *ae = matrix_op_eval (op, *ae, be);
2312 else if (is_scalar (a))
2314 double ae = to_scalar (a);
2315 for (size_t r = 0; r < b->size1; r++)
2316 for (size_t c = 0; c < b->size2; c++)
2318 double *be = gsl_matrix_ptr (b, r, c);
2319 *be = matrix_op_eval (op, ae, *be);
2323 else if (a->size1 == b->size1 && a->size2 == b->size2)
2325 for (size_t r = 0; r < a->size1; r++)
2326 for (size_t c = 0; c < a->size2; c++)
2328 double *ae = gsl_matrix_ptr (a, r, c);
2329 double be = gsl_matrix_get (b, r, c);
2330 *ae = matrix_op_eval (op, *ae, be);
2336 msg (SE, _("Operands to %s must have the same dimensions or one "
2337 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
2338 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2344 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2346 if (is_scalar (a) || is_scalar (b))
2347 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2349 if (a->size2 != b->size1)
2351 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
2352 "not conformable for multiplication."),
2353 a->size1, a->size2, b->size1, b->size2);
2357 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2358 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2363 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2365 gsl_matrix *tmp = *a;
2371 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2374 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2375 swap_matrix (z, tmp);
2379 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2381 mul_matrix (x, *x, *x, tmp);
2385 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2388 if (x->size1 != x->size2)
2390 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2391 "the left-hand size, not one with dimensions %zu×%zu."),
2392 x->size1, x->size2);
2397 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2398 "right-hand side, not a matrix with dimensions %zu×%zu."),
2399 b->size1, b->size2);
2402 double bf = to_scalar (b);
2403 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2405 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2406 "or outside the valid range."), bf);
2411 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2413 gsl_matrix_set_identity (y);
2417 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2419 for (unsigned long int n = labs (bl); n > 1; n /= 2)
2422 mul_matrix (&y, x, y, &t);
2423 square_matrix (&x, &t);
2426 square_matrix (&x, &t);
2428 mul_matrix (&y, x, y, &t);
2432 /* Garbage collection.
2434 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2435 a permutation of them. We are returning one of them; that one must not be
2436 destroyed. We must not destroy 'x_' because the caller owns it. */
2438 gsl_matrix_free (y_);
2440 gsl_matrix_free (t_);
2446 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2449 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2451 msg (SE, _("All operands of : operator must be scalars."));
2455 long int start = to_scalar (start_);
2456 long int end = to_scalar (end_);
2457 long int by = by_ ? to_scalar (by_) : 1;
2461 msg (SE, _("The increment operand to : must be nonzero."));
2465 long int n = (end >= start && by > 0 ? (end - start + by) / by
2466 : end <= start && by < 0 ? (start - end - by) / -by
2468 gsl_matrix *m = gsl_matrix_alloc (1, n);
2469 for (long int i = 0; i < n; i++)
2470 gsl_matrix_set (m, 0, i, start + i * by);
2475 matrix_expr_evaluate_not (gsl_matrix *a)
2477 for (size_t r = 0; r < a->size1; r++)
2478 for (size_t c = 0; c < a->size2; c++)
2480 double *ae = gsl_matrix_ptr (a, r, c);
2487 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2489 if (a->size1 != b->size1)
2491 if (!a->size1 || !a->size2)
2493 else if (!b->size1 || !b->size2)
2496 msg (SE, _("All columns in a matrix must have the same number of rows, "
2497 "but this tries to paste matrices with %zu and %zu rows."),
2498 a->size1, b->size1);
2502 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2503 for (size_t y = 0; y < a->size1; y++)
2505 for (size_t x = 0; x < a->size2; x++)
2506 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2507 for (size_t x = 0; x < b->size2; x++)
2508 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2514 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2516 if (a->size2 != b->size2)
2518 if (!a->size1 || !a->size2)
2520 else if (!b->size1 || !b->size2)
2523 msg (SE, _("All rows in a matrix must have the same number of columns, "
2524 "but this tries to stack matrices with %zu and %zu columns."),
2525 a->size2, b->size2);
2529 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2530 for (size_t x = 0; x < a->size2; x++)
2532 for (size_t y = 0; y < a->size1; y++)
2533 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2534 for (size_t y = 0; y < b->size1; y++)
2535 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2541 matrix_to_vector (gsl_matrix *m)
2544 gsl_vector v = to_vector (m);
2545 assert (v.block == m->block || !v.block);
2549 return xmemdup (&v, sizeof v);
2563 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2566 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2567 enum index_type index_type, size_t other_size,
2568 struct index_vector *iv)
2577 msg (SE, _("Vector index must be scalar or vector, not a "
2579 m->size1, m->size2);
2583 msg (SE, _("Matrix row index must be scalar or vector, not a "
2585 m->size1, m->size2);
2589 msg (SE, _("Matrix column index must be scalar or vector, not a "
2591 m->size1, m->size2);
2597 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2598 *iv = (struct index_vector) {
2599 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2602 for (size_t i = 0; i < v.size; i++)
2604 double index = gsl_vector_get (&v, i);
2605 if (index < 1 || index >= size + 1)
2610 msg (SE, _("Index %g is out of range for vector "
2611 "with %zu elements."), index, size);
2615 msg (SE, _("%g is not a valid row index for "
2616 "a %zu×%zu matrix."),
2617 index, size, other_size);
2621 msg (SE, _("%g is not a valid column index for "
2622 "a %zu×%zu matrix."),
2623 index, other_size, size);
2630 iv->indexes[i] = index - 1;
2636 *iv = (struct index_vector) {
2637 .indexes = xnmalloc (size, sizeof *iv->indexes),
2640 for (size_t i = 0; i < size; i++)
2647 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2649 if (!is_vector (sm))
2651 msg (SE, _("Vector index operator may not be applied to "
2652 "a %zu×%zu matrix."),
2653 sm->size1, sm->size2);
2661 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2663 if (!matrix_expr_evaluate_vec_all (sm))
2666 gsl_vector sv = to_vector (sm);
2667 struct index_vector iv;
2668 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2671 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2672 sm->size1 == 1 ? iv.n : 1);
2673 gsl_vector dv = to_vector (dm);
2674 for (size_t dx = 0; dx < iv.n; dx++)
2676 size_t sx = iv.indexes[dx];
2677 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2685 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2688 struct index_vector iv0;
2689 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2692 struct index_vector iv1;
2693 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2700 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2701 for (size_t dy = 0; dy < iv0.n; dy++)
2703 size_t sy = iv0.indexes[dy];
2705 for (size_t dx = 0; dx < iv1.n; dx++)
2707 size_t sx = iv1.indexes[dx];
2708 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2716 #define F(NAME, PROTOTYPE) \
2717 static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2718 const char *name, gsl_matrix *[], size_t, \
2719 matrix_proto_##PROTOTYPE *);
2724 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2726 if (!is_scalar (subs[index]))
2728 msg (SE, _("Function %s argument %zu must be a scalar, "
2729 "not a %zu×%zu matrix."),
2730 name, index + 1, subs[index]->size1, subs[index]->size2);
2737 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2739 if (!is_vector (subs[index]))
2741 msg (SE, _("Function %s argument %zu must be a vector, "
2742 "not a %zu×%zu matrix."),
2743 name, index + 1, subs[index]->size1, subs[index]->size2);
2750 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2752 for (size_t i = 0; i < n_subs; i++)
2754 if (!check_scalar_arg (name, subs, i))
2756 d[i] = to_scalar (subs[i]);
2762 matrix_expr_evaluate_m_d (const char *name,
2763 gsl_matrix *subs[], size_t n_subs,
2764 matrix_proto_m_d *f)
2766 assert (n_subs == 1);
2769 return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2773 matrix_expr_evaluate_m_dd (const char *name,
2774 gsl_matrix *subs[], size_t n_subs,
2775 matrix_proto_m_dd *f)
2777 assert (n_subs == 2);
2780 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2784 matrix_expr_evaluate_m_ddd (const char *name,
2785 gsl_matrix *subs[], size_t n_subs,
2786 matrix_proto_m_ddd *f)
2788 assert (n_subs == 3);
2791 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2795 matrix_expr_evaluate_m_m (const char *name UNUSED,
2796 gsl_matrix *subs[], size_t n_subs,
2797 matrix_proto_m_m *f)
2799 assert (n_subs == 1);
2804 matrix_expr_evaluate_m_md (const char *name UNUSED,
2805 gsl_matrix *subs[], size_t n_subs,
2806 matrix_proto_m_md *f)
2808 assert (n_subs == 2);
2809 if (!check_scalar_arg (name, subs, 1))
2811 return f (subs[0], to_scalar (subs[1]));
2815 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2816 gsl_matrix *subs[], size_t n_subs,
2817 matrix_proto_m_mdd *f)
2819 assert (n_subs == 3);
2820 if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2822 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2826 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2827 gsl_matrix *subs[], size_t n_subs,
2828 matrix_proto_m_mm *f)
2830 assert (n_subs == 2);
2831 return f (subs[0], subs[1]);
2835 matrix_expr_evaluate_m_v (const char *name,
2836 gsl_matrix *subs[], size_t n_subs,
2837 matrix_proto_m_v *f)
2839 assert (n_subs == 1);
2840 if (!check_vector_arg (name, subs, 0))
2842 gsl_vector v = to_vector (subs[0]);
2847 matrix_expr_evaluate_d_m (const char *name UNUSED,
2848 gsl_matrix *subs[], size_t n_subs,
2849 matrix_proto_d_m *f)
2851 assert (n_subs == 1);
2853 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2854 gsl_matrix_set (m, 0, 0, f (subs[0]));
2859 matrix_expr_evaluate_m_any (const char *name UNUSED,
2860 gsl_matrix *subs[], size_t n_subs,
2861 matrix_proto_m_any *f)
2863 return f (subs, n_subs);
2867 matrix_expr_evaluate_IDENT (const char *name,
2868 gsl_matrix *subs[], size_t n_subs,
2869 matrix_proto_IDENT *f)
2871 assert (n_subs <= 2);
2874 if (!to_scalar_args (name, subs, n_subs, d))
2877 return f (d[0], d[n_subs - 1]);
2881 matrix_expr_evaluate (const struct matrix_expr *e)
2883 if (e->op == MOP_NUMBER)
2885 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2886 gsl_matrix_set (m, 0, 0, e->number);
2889 else if (e->op == MOP_VARIABLE)
2891 const gsl_matrix *src = e->variable->value;
2894 msg (SE, _("Uninitialized variable %s used in expression."),
2899 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2900 gsl_matrix_memcpy (dst, src);
2903 else if (e->op == MOP_EOF)
2905 struct dfm_reader *reader = read_file_open (e->eof);
2906 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2907 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2911 enum { N_LOCAL = 3 };
2912 gsl_matrix *local_subs[N_LOCAL];
2913 gsl_matrix **subs = (e->n_subs < N_LOCAL
2915 : xmalloc (e->n_subs * sizeof *subs));
2917 for (size_t i = 0; i < e->n_subs; i++)
2919 subs[i] = matrix_expr_evaluate (e->subs[i]);
2922 for (size_t j = 0; j < i; j++)
2923 gsl_matrix_free (subs[j]);
2924 if (subs != local_subs)
2930 gsl_matrix *result = NULL;
2933 #define F(NAME, PROTOTYPE) \
2934 case MOP_F_##NAME: \
2935 result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \
2937 matrix_eval_##NAME); \
2943 gsl_matrix_scale (subs[0], -1.0);
2961 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2965 result = matrix_expr_evaluate_not (subs[0]);
2969 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2973 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2977 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2981 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2984 case MOP_PASTE_HORZ:
2985 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2988 case MOP_PASTE_VERT:
2989 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2993 result = gsl_matrix_alloc (0, 0);
2997 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
3001 result = matrix_expr_evaluate_vec_all (subs[0]);
3005 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
3009 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
3013 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
3022 for (size_t i = 0; i < e->n_subs; i++)
3023 if (subs[i] != result)
3024 gsl_matrix_free (subs[i]);
3025 if (subs != local_subs)
3031 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
3034 gsl_matrix *m = matrix_expr_evaluate (e);
3040 msg (SE, _("Expression for %s must evaluate to scalar, "
3041 "not a %zu×%zu matrix."),
3042 context, m->size1, m->size2);
3043 gsl_matrix_free (m);
3048 gsl_matrix_free (m);
3053 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
3057 if (!matrix_expr_evaluate_scalar (e, context, &d))
3061 if (d < LONG_MIN || d > LONG_MAX)
3063 msg (SE, _("Expression for %s is outside the integer range."), context);
3070 struct matrix_lvalue
3072 struct matrix_var *var;
3073 struct matrix_expr *indexes[2];
3078 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
3082 for (size_t i = 0; i < lvalue->n_indexes; i++)
3083 matrix_expr_destroy (lvalue->indexes[i]);
3088 static struct matrix_lvalue *
3089 matrix_lvalue_parse (struct matrix_state *s)
3091 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3092 if (!lex_force_id (s->lexer))
3094 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3095 if (lex_next_token (s->lexer, 1) == T_LPAREN)
3099 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3104 lex_get_n (s->lexer, 2);
3106 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3108 lvalue->n_indexes++;
3110 if (lex_match (s->lexer, T_COMMA))
3112 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3114 lvalue->n_indexes++;
3116 if (!lex_force_match (s->lexer, T_RPAREN))
3122 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3128 matrix_lvalue_destroy (lvalue);
3133 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3134 enum index_type index_type, size_t other_size,
3135 struct index_vector *iv)
3140 m = matrix_expr_evaluate (e);
3147 bool ok = matrix_normalize_index_vector (m, size, index_type,
3149 gsl_matrix_free (m);
3154 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3155 struct index_vector *iv, gsl_matrix *sm)
3157 gsl_vector dv = to_vector (lvalue->var->value);
3159 /* Convert source matrix 'sm' to source vector 'sv'. */
3160 if (!is_vector (sm))
3162 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
3163 sm->size1, sm->size2);
3166 gsl_vector sv = to_vector (sm);
3167 if (iv->n != sv.size)
3169 msg (SE, _("Can't assign %zu-element vector "
3170 "to %zu-element subvector."), sv.size, iv->n);
3174 /* Assign elements. */
3175 for (size_t x = 0; x < iv->n; x++)
3176 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3181 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3182 struct index_vector *iv0,
3183 struct index_vector *iv1, gsl_matrix *sm)
3185 gsl_matrix *dm = lvalue->var->value;
3187 /* Convert source matrix 'sm' to source vector 'sv'. */
3188 if (iv0->n != sm->size1)
3190 msg (SE, _("Row index vector for assignment to %s has %zu elements "
3191 "but source matrix has %zu rows."),
3192 lvalue->var->name, iv0->n, sm->size1);
3195 if (iv1->n != sm->size2)
3197 msg (SE, _("Column index vector for assignment to %s has %zu elements "
3198 "but source matrix has %zu columns."),
3199 lvalue->var->name, iv1->n, sm->size2);
3203 /* Assign elements. */
3204 for (size_t y = 0; y < iv0->n; y++)
3206 size_t f0 = iv0->indexes[y];
3207 for (size_t x = 0; x < iv1->n; x++)
3209 size_t f1 = iv1->indexes[x];
3210 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3216 /* Takes ownership of M. */
3218 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3219 struct index_vector *iv0, struct index_vector *iv1,
3222 if (!lvalue->n_indexes)
3224 gsl_matrix_free (lvalue->var->value);
3225 lvalue->var->value = sm;
3230 bool ok = (lvalue->n_indexes == 1
3231 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3232 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3233 free (iv0->indexes);
3234 free (iv1->indexes);
3235 gsl_matrix_free (sm);
3241 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3242 struct index_vector *iv0,
3243 struct index_vector *iv1)
3245 *iv0 = INDEX_VECTOR_INIT;
3246 *iv1 = INDEX_VECTOR_INIT;
3247 if (!lvalue->n_indexes)
3250 /* Validate destination matrix exists and has the right shape. */
3251 gsl_matrix *dm = lvalue->var->value;
3254 msg (SE, _("Undefined variable %s."), lvalue->var->name);
3257 else if (dm->size1 == 0 || dm->size2 == 0)
3259 msg (SE, _("Cannot index %zu×%zu matrix."),
3260 dm->size1, dm->size2);
3263 else if (lvalue->n_indexes == 1)
3265 if (!is_vector (dm))
3267 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
3268 dm->size1, dm->size2, lvalue->var->name);
3271 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3272 MAX (dm->size1, dm->size2),
3277 assert (lvalue->n_indexes == 2);
3278 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3279 IV_ROW, dm->size2, iv0))
3282 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3283 IV_COLUMN, dm->size1, iv1))
3285 free (iv0->indexes);
3292 /* Takes ownership of M. */
3294 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3296 struct index_vector iv0, iv1;
3297 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3299 gsl_matrix_free (sm);
3303 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3307 /* Matrix command. */
3311 struct matrix_cmd **commands;
3315 static void matrix_cmds_uninit (struct matrix_cmds *);
3319 enum matrix_cmd_type
3342 struct compute_command
3344 struct matrix_lvalue *lvalue;
3345 struct matrix_expr *rvalue;
3349 struct print_command
3351 struct matrix_expr *expression;
3352 bool use_default_format;
3353 struct fmt_spec format;
3355 int space; /* -1 means NEWPAGE. */
3359 struct string_array literals; /* CLABELS/RLABELS. */
3360 struct matrix_expr *expr; /* CNAMES/RNAMES. */
3366 struct do_if_command
3370 struct matrix_expr *condition;
3371 struct matrix_cmds commands;
3381 struct matrix_var *var;
3382 struct matrix_expr *start;
3383 struct matrix_expr *end;
3384 struct matrix_expr *increment;
3386 /* Loop conditions. */
3387 struct matrix_expr *top_condition;
3388 struct matrix_expr *bottom_condition;
3391 struct matrix_cmds commands;
3395 struct display_command
3397 struct matrix_state *state;
3403 struct release_command
3405 struct matrix_var **vars;
3412 struct matrix_expr *expression;
3413 struct save_file *sf;
3419 struct read_file *rf;
3420 struct matrix_lvalue *dst;
3421 struct matrix_expr *size;
3423 enum fmt_type format;
3430 struct write_command
3432 struct write_file *wf;
3433 struct matrix_expr *expression;
3436 /* If this is nonnull, WRITE uses this format.
3438 If this is NULL, WRITE uses free-field format with as many
3439 digits of precision as needed. */
3440 struct fmt_spec *format;
3449 struct matrix_lvalue *dst;
3450 struct dataset *dataset;
3451 struct file_handle *file;
3453 struct var_syntax *vars;
3455 struct matrix_var *names;
3457 /* Treatment of missing values. */
3462 MGET_ERROR, /* Flag error on command. */
3463 MGET_ACCEPT, /* Accept without change, user-missing only. */
3464 MGET_OMIT, /* Drop this case. */
3465 MGET_RECODE /* Recode to 'substitute'. */
3468 double substitute; /* MGET_RECODE only. */
3474 struct msave_command
3476 struct msave_common *common;
3478 struct matrix_expr *expr;
3479 const char *rowtype;
3480 struct matrix_expr *factors;
3481 struct matrix_expr *splits;
3487 struct matrix_state *state;
3488 struct file_handle *file;
3489 struct stringi_set rowtypes;
3493 struct eigen_command
3495 struct matrix_expr *expr;
3496 struct matrix_var *evec;
3497 struct matrix_var *eval;
3501 struct setdiag_command
3503 struct matrix_var *dst;
3504 struct matrix_expr *expr;
3510 struct matrix_expr *expr;
3511 struct matrix_var *u;
3512 struct matrix_var *s;
3513 struct matrix_var *v;
3519 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3520 static bool matrix_cmd_execute (struct matrix_cmd *);
3521 static void matrix_cmd_destroy (struct matrix_cmd *);
3524 static struct matrix_cmd *
3525 matrix_parse_compute (struct matrix_state *s)
3527 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3528 *cmd = (struct matrix_cmd) {
3529 .type = MCMD_COMPUTE,
3530 .compute = { .lvalue = NULL },
3533 struct compute_command *compute = &cmd->compute;
3534 compute->lvalue = matrix_lvalue_parse (s);
3535 if (!compute->lvalue)
3538 if (!lex_force_match (s->lexer, T_EQUALS))
3541 compute->rvalue = matrix_parse_expr (s);
3542 if (!compute->rvalue)
3548 matrix_cmd_destroy (cmd);
3553 matrix_cmd_execute_compute (struct compute_command *compute)
3555 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3559 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3563 print_labels_destroy (struct print_labels *labels)
3567 string_array_destroy (&labels->literals);
3568 matrix_expr_destroy (labels->expr);
3574 parse_literal_print_labels (struct matrix_state *s,
3575 struct print_labels **labelsp)
3577 lex_match (s->lexer, T_EQUALS);
3578 print_labels_destroy (*labelsp);
3579 *labelsp = xzalloc (sizeof **labelsp);
3580 while (lex_token (s->lexer) != T_SLASH
3581 && lex_token (s->lexer) != T_ENDCMD
3582 && lex_token (s->lexer) != T_STOP)
3584 struct string label = DS_EMPTY_INITIALIZER;
3585 while (lex_token (s->lexer) != T_COMMA
3586 && lex_token (s->lexer) != T_SLASH
3587 && lex_token (s->lexer) != T_ENDCMD
3588 && lex_token (s->lexer) != T_STOP)
3590 if (!ds_is_empty (&label))
3591 ds_put_byte (&label, ' ');
3593 if (lex_token (s->lexer) == T_STRING)
3594 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3597 char *rep = lex_next_representation (s->lexer, 0, 0);
3598 ds_put_cstr (&label, rep);
3603 string_array_append_nocopy (&(*labelsp)->literals,
3604 ds_steal_cstr (&label));
3606 if (!lex_match (s->lexer, T_COMMA))
3612 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3614 lex_match (s->lexer, T_EQUALS);
3615 struct matrix_expr *e = matrix_parse_exp (s);
3619 print_labels_destroy (*labelsp);
3620 *labelsp = xzalloc (sizeof **labelsp);
3621 (*labelsp)->expr = e;
3625 static struct matrix_cmd *
3626 matrix_parse_print (struct matrix_state *s)
3628 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3629 *cmd = (struct matrix_cmd) {
3632 .use_default_format = true,
3636 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3639 for (size_t i = 0; ; i++)
3641 enum token_type t = lex_next_token (s->lexer, i);
3642 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3644 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3646 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3649 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3654 cmd->print.expression = matrix_parse_exp (s);
3655 if (!cmd->print.expression)
3659 while (lex_match (s->lexer, T_SLASH))
3661 if (lex_match_id (s->lexer, "FORMAT"))
3663 lex_match (s->lexer, T_EQUALS);
3664 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3666 cmd->print.use_default_format = false;
3668 else if (lex_match_id (s->lexer, "TITLE"))
3670 lex_match (s->lexer, T_EQUALS);
3671 if (!lex_force_string (s->lexer))
3673 free (cmd->print.title);
3674 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3677 else if (lex_match_id (s->lexer, "SPACE"))
3679 lex_match (s->lexer, T_EQUALS);
3680 if (lex_match_id (s->lexer, "NEWPAGE"))
3681 cmd->print.space = -1;
3682 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3684 cmd->print.space = lex_integer (s->lexer);
3690 else if (lex_match_id (s->lexer, "RLABELS"))
3691 parse_literal_print_labels (s, &cmd->print.rlabels);
3692 else if (lex_match_id (s->lexer, "CLABELS"))
3693 parse_literal_print_labels (s, &cmd->print.clabels);
3694 else if (lex_match_id (s->lexer, "RNAMES"))
3696 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3699 else if (lex_match_id (s->lexer, "CNAMES"))
3701 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3706 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3707 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3715 matrix_cmd_destroy (cmd);
3720 matrix_is_integer (const gsl_matrix *m)
3722 for (size_t y = 0; y < m->size1; y++)
3723 for (size_t x = 0; x < m->size2; x++)
3725 double d = gsl_matrix_get (m, y, x);
3733 matrix_max_magnitude (const gsl_matrix *m)
3736 for (size_t y = 0; y < m->size1; y++)
3737 for (size_t x = 0; x < m->size2; x++)
3739 double d = fabs (gsl_matrix_get (m, y, x));
3747 format_fits (struct fmt_spec format, double x)
3749 char *s = data_out (&(union value) { .f = x }, NULL,
3750 &format, settings_get_fmt_settings ());
3751 bool fits = *s != '*' && !strchr (s, 'E');
3756 static struct fmt_spec
3757 default_format (const gsl_matrix *m, int *log_scale)
3761 double max = matrix_max_magnitude (m);
3763 if (matrix_is_integer (m))
3764 for (int w = 1; w <= 12; w++)
3766 struct fmt_spec format = { .type = FMT_F, .w = w };
3767 if (format_fits (format, -max))
3771 if (max >= 1e9 || max <= 1e-4)
3774 snprintf (s, sizeof s, "%.1e", max);
3776 const char *e = strchr (s, 'e');
3778 *log_scale = atoi (e + 1);
3781 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3785 trimmed_string (double d)
3787 struct substring s = ss_buffer ((char *) &d, sizeof d);
3788 ss_rtrim (&s, ss_cstr (" "));
3789 return ss_xstrdup (s);
3792 static struct string_array *
3793 print_labels_get (const struct print_labels *labels, size_t n,
3794 const char *prefix, bool truncate)
3799 struct string_array *out = xzalloc (sizeof *out);
3800 if (labels->literals.n)
3801 string_array_clone (out, &labels->literals);
3802 else if (labels->expr)
3804 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3805 if (m && is_vector (m))
3807 gsl_vector v = to_vector (m);
3808 for (size_t i = 0; i < v.size; i++)
3809 string_array_append_nocopy (out, trimmed_string (
3810 gsl_vector_get (&v, i)));
3812 gsl_matrix_free (m);
3818 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3821 string_array_append (out, "");
3824 string_array_delete (out, out->n - 1);
3827 for (size_t i = 0; i < out->n; i++)
3829 char *s = out->strings[i];
3830 s[strnlen (s, 8)] = '\0';
3837 matrix_cmd_print_space (int space)
3840 output_item_submit (page_break_item_create ());
3841 for (int i = 0; i < space; i++)
3842 output_log ("%s", "");
3846 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3847 struct fmt_spec format, int log_scale)
3849 matrix_cmd_print_space (print->space);
3851 output_log ("%s", print->title);
3853 output_log (" 10 ** %d X", log_scale);
3855 struct string_array *clabels = print_labels_get (print->clabels,
3856 m->size2, "col", true);
3857 if (clabels && format.w < 8)
3859 struct string_array *rlabels = print_labels_get (print->rlabels,
3860 m->size1, "row", true);
3864 struct string line = DS_EMPTY_INITIALIZER;
3866 ds_put_byte_multiple (&line, ' ', 8);
3867 for (size_t x = 0; x < m->size2; x++)
3868 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3869 output_log_nocopy (ds_steal_cstr (&line));
3872 double scale = pow (10.0, log_scale);
3873 bool numeric = fmt_is_numeric (format.type);
3874 for (size_t y = 0; y < m->size1; y++)
3876 struct string line = DS_EMPTY_INITIALIZER;
3878 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3880 for (size_t x = 0; x < m->size2; x++)
3882 double f = gsl_matrix_get (m, y, x);
3884 ? data_out (&(union value) { .f = f / scale}, NULL,
3885 &format, settings_get_fmt_settings ())
3886 : trimmed_string (f));
3887 ds_put_format (&line, " %s", s);
3890 output_log_nocopy (ds_steal_cstr (&line));
3893 string_array_destroy (rlabels);
3895 string_array_destroy (clabels);
3900 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3901 const struct print_labels *print_labels, size_t n,
3902 const char *name, const char *prefix)
3904 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3906 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3907 for (size_t i = 0; i < n; i++)
3908 pivot_category_create_leaf (
3910 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3911 : pivot_value_new_integer (i + 1)));
3913 d->hide_all_labels = true;
3914 string_array_destroy (labels);
3919 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3920 struct fmt_spec format, int log_scale)
3922 struct pivot_table *table = pivot_table_create__ (
3923 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3926 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3928 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3929 N_("Columns"), "col");
3931 struct pivot_footnote *footnote = NULL;
3934 char *s = xasprintf ("× 10**%d", log_scale);
3935 footnote = pivot_table_create_footnote (
3936 table, pivot_value_new_user_text_nocopy (s));
3939 double scale = pow (10.0, log_scale);
3940 bool numeric = fmt_is_numeric (format.type);
3941 for (size_t y = 0; y < m->size1; y++)
3942 for (size_t x = 0; x < m->size2; x++)
3944 double f = gsl_matrix_get (m, y, x);
3945 struct pivot_value *v;
3948 v = pivot_value_new_number (f / scale);
3949 v->numeric.format = format;
3952 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3954 pivot_value_add_footnote (v, footnote);
3955 pivot_table_put2 (table, y, x, v);
3958 pivot_table_submit (table);
3962 matrix_cmd_execute_print (const struct print_command *print)
3964 if (print->expression)
3966 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3971 struct fmt_spec format = (print->use_default_format
3972 ? default_format (m, &log_scale)
3975 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3976 matrix_cmd_print_text (print, m, format, log_scale);
3978 matrix_cmd_print_tables (print, m, format, log_scale);
3980 gsl_matrix_free (m);
3984 matrix_cmd_print_space (print->space);
3986 output_log ("%s", print->title);
3993 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3994 const char *command_name,
3995 const char *stop1, const char *stop2)
3997 lex_end_of_command (s->lexer);
3998 lex_discard_rest_of_command (s->lexer);
4000 size_t allocated = 0;
4003 while (lex_token (s->lexer) == T_ENDCMD)
4006 if (lex_at_phrase (s->lexer, stop1)
4007 || (stop2 && lex_at_phrase (s->lexer, stop2)))
4010 if (lex_at_phrase (s->lexer, "END MATRIX"))
4012 msg (SE, _("Premature END MATRIX within %s."), command_name);
4016 struct matrix_cmd *cmd = matrix_parse_command (s);
4020 if (c->n >= allocated)
4021 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
4022 c->commands[c->n++] = cmd;
4027 matrix_parse_do_if_clause (struct matrix_state *s,
4028 struct do_if_command *ifc,
4029 bool parse_condition,
4030 size_t *allocated_clauses)
4032 if (ifc->n_clauses >= *allocated_clauses)
4033 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
4034 sizeof *ifc->clauses);
4035 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
4036 *c = (struct do_if_clause) { .condition = NULL };
4038 if (parse_condition)
4040 c->condition = matrix_parse_expr (s);
4045 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
4048 static struct matrix_cmd *
4049 matrix_parse_do_if (struct matrix_state *s)
4051 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4052 *cmd = (struct matrix_cmd) {
4054 .do_if = { .n_clauses = 0 }
4057 size_t allocated_clauses = 0;
4060 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
4063 while (lex_match_phrase (s->lexer, "ELSE IF"));
4065 if (lex_match_id (s->lexer, "ELSE")
4066 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
4069 if (!lex_match_phrase (s->lexer, "END IF"))
4074 matrix_cmd_destroy (cmd);
4079 matrix_cmd_execute_do_if (struct do_if_command *cmd)
4081 for (size_t i = 0; i < cmd->n_clauses; i++)
4083 struct do_if_clause *c = &cmd->clauses[i];
4087 if (!matrix_expr_evaluate_scalar (c->condition,
4088 i ? "ELSE IF" : "DO IF",
4093 for (size_t j = 0; j < c->commands.n; j++)
4094 if (!matrix_cmd_execute (c->commands.commands[j]))
4101 static struct matrix_cmd *
4102 matrix_parse_loop (struct matrix_state *s)
4104 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4105 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4107 struct loop_command *loop = &cmd->loop;
4108 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4110 struct substring name = lex_tokss (s->lexer);
4111 loop->var = matrix_var_lookup (s, name);
4113 loop->var = matrix_var_create (s, name);
4118 loop->start = matrix_parse_expr (s);
4119 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4122 loop->end = matrix_parse_expr (s);
4126 if (lex_match (s->lexer, T_BY))
4128 loop->increment = matrix_parse_expr (s);
4129 if (!loop->increment)
4134 if (lex_match_id (s->lexer, "IF"))
4136 loop->top_condition = matrix_parse_expr (s);
4137 if (!loop->top_condition)
4141 bool was_in_loop = s->in_loop;
4143 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4145 s->in_loop = was_in_loop;
4149 if (!lex_match_phrase (s->lexer, "END LOOP"))
4152 if (lex_match_id (s->lexer, "IF"))
4154 loop->bottom_condition = matrix_parse_expr (s);
4155 if (!loop->bottom_condition)
4162 matrix_cmd_destroy (cmd);
4167 matrix_cmd_execute_loop (struct loop_command *cmd)
4171 long int increment = 1;
4174 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4175 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4177 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4181 if (increment > 0 ? value > end
4182 : increment < 0 ? value < end
4187 int mxloops = settings_get_mxloops ();
4188 for (int i = 0; i < mxloops; i++)
4193 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4195 gsl_matrix_free (cmd->var->value);
4196 cmd->var->value = NULL;
4198 if (!cmd->var->value)
4199 cmd->var->value = gsl_matrix_alloc (1, 1);
4201 gsl_matrix_set (cmd->var->value, 0, 0, value);
4204 if (cmd->top_condition)
4207 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4212 for (size_t j = 0; j < cmd->commands.n; j++)
4213 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4216 if (cmd->bottom_condition)
4219 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4220 "END LOOP IF", &d) || d > 0)
4227 if (increment > 0 ? value > end : value < end)
4233 static struct matrix_cmd *
4234 matrix_parse_break (struct matrix_state *s)
4238 msg (SE, _("BREAK not inside LOOP."));
4242 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4243 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4247 static struct matrix_cmd *
4248 matrix_parse_display (struct matrix_state *s)
4250 bool dictionary = false;
4251 bool status = false;
4252 while (lex_token (s->lexer) == T_ID)
4254 if (lex_match_id (s->lexer, "DICTIONARY"))
4256 else if (lex_match_id (s->lexer, "STATUS"))
4260 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4264 if (!dictionary && !status)
4265 dictionary = status = true;
4267 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4268 *cmd = (struct matrix_cmd) {
4269 .type = MCMD_DISPLAY,
4272 .dictionary = dictionary,
4280 compare_matrix_var_pointers (const void *a_, const void *b_)
4282 const struct matrix_var *const *ap = a_;
4283 const struct matrix_var *const *bp = b_;
4284 const struct matrix_var *a = *ap;
4285 const struct matrix_var *b = *bp;
4286 return strcmp (a->name, b->name);
4290 matrix_cmd_execute_display (struct display_command *cmd)
4292 const struct matrix_state *s = cmd->state;
4294 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4295 pivot_dimension_create (
4296 table, PIVOT_AXIS_COLUMN, N_("Property"),
4297 N_("Rows"), N_("Columns"), N_("Size (kB)"));
4299 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4302 const struct matrix_var *var;
4303 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4305 vars[n_vars++] = var;
4306 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4308 struct pivot_dimension *rows = pivot_dimension_create (
4309 table, PIVOT_AXIS_ROW, N_("Variable"));
4310 for (size_t i = 0; i < n_vars; i++)
4312 const struct matrix_var *var = vars[i];
4313 pivot_category_create_leaf (
4314 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4316 size_t r = var->value->size1;
4317 size_t c = var->value->size2;
4318 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4319 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4320 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4322 pivot_table_submit (table);
4325 static struct matrix_cmd *
4326 matrix_parse_release (struct matrix_state *s)
4328 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4329 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4331 size_t allocated_vars = 0;
4332 while (lex_token (s->lexer) == T_ID)
4334 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4337 if (cmd->release.n_vars >= allocated_vars)
4338 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4339 sizeof *cmd->release.vars);
4340 cmd->release.vars[cmd->release.n_vars++] = var;
4343 lex_error (s->lexer, _("Variable name expected."));
4346 if (!lex_match (s->lexer, T_COMMA))
4354 matrix_cmd_execute_release (struct release_command *cmd)
4356 for (size_t i = 0; i < cmd->n_vars; i++)
4358 struct matrix_var *var = cmd->vars[i];
4359 gsl_matrix_free (var->value);
4364 static struct save_file *
4365 save_file_create (struct matrix_state *s, struct file_handle *fh,
4366 struct string_array *variables,
4367 struct matrix_expr *names,
4368 struct stringi_set *strings)
4370 for (size_t i = 0; i < s->n_save_files; i++)
4372 struct save_file *sf = s->save_files[i];
4377 string_array_destroy (variables);
4378 matrix_expr_destroy (names);
4379 stringi_set_destroy (strings);
4385 struct save_file *sf = xmalloc (sizeof *sf);
4386 *sf = (struct save_file) {
4388 .dataset = s->dataset,
4389 .variables = *variables,
4391 .strings = STRINGI_SET_INITIALIZER (sf->strings),
4394 stringi_set_swap (&sf->strings, strings);
4395 stringi_set_destroy (strings);
4397 s->save_files = xrealloc (s->save_files,
4398 (s->n_save_files + 1) * sizeof *s->save_files);
4399 s->save_files[s->n_save_files++] = sf;
4403 static struct casewriter *
4404 save_file_open (struct save_file *sf, gsl_matrix *m)
4406 if (sf->writer || sf->error)
4410 size_t n_variables = caseproto_get_n_widths (
4411 casewriter_get_proto (sf->writer));
4412 if (m->size2 != n_variables)
4414 msg (ME, _("The first SAVE to %s within this matrix program "
4415 "had %zu columns, so a %zu×%zu matrix cannot be "
4417 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
4418 n_variables, m->size1, m->size2);
4426 struct dictionary *dict = dict_create (get_default_encoding ());
4428 /* Fill 'names' with user-specified names if there were any, either from
4429 sf->variables or sf->names. */
4430 struct string_array names = { .n = 0 };
4431 if (sf->variables.n)
4432 string_array_clone (&names, &sf->variables);
4435 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
4436 if (nm && is_vector (nm))
4438 gsl_vector nv = to_vector (nm);
4439 for (size_t i = 0; i < nv.size; i++)
4441 char *name = trimmed_string (gsl_vector_get (&nv, i));
4442 if (dict_id_is_valid (dict, name, true))
4443 string_array_append_nocopy (&names, name);
4448 gsl_matrix_free (nm);
4451 struct stringi_set strings;
4452 stringi_set_clone (&strings, &sf->strings);
4454 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4459 name = names.strings[i];
4462 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4466 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4467 struct variable *var = dict_create_var (dict, name, width);
4470 msg (ME, _("Duplicate variable name %s in SAVE statement."),
4476 if (!stringi_set_is_empty (&strings))
4478 const char *example = stringi_set_node_get_string (
4479 stringi_set_first (&strings));
4480 msg (ME, ngettext ("STRINGS specified a variable %s, but no variable "
4481 "with that name was found on SAVE.",
4482 "STRINGS specified %2$zu variables, including %1$s, "
4483 "whose names were not found on SAVE.",
4484 stringi_set_count (&strings)),
4485 example, stringi_set_count (&strings));
4488 stringi_set_destroy (&strings);
4489 string_array_destroy (&names);
4498 if (sf->file == fh_inline_file ())
4499 sf->writer = autopaging_writer_create (dict_get_proto (dict));
4501 sf->writer = any_writer_open (sf->file, dict);
4513 save_file_destroy (struct save_file *sf)
4517 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
4519 dataset_set_dict (sf->dataset, sf->dict);
4520 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
4525 fh_unref (sf->file);
4526 string_array_destroy (&sf->variables);
4527 matrix_expr_destroy (sf->names);
4528 stringi_set_destroy (&sf->strings);
4529 casewriter_destroy (sf->writer);
4530 dict_unref (sf->dict);
4535 static struct matrix_cmd *
4536 matrix_parse_save (struct matrix_state *s)
4538 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4539 *cmd = (struct matrix_cmd) {
4541 .save = { .expression = NULL },
4544 struct file_handle *fh = NULL;
4545 struct save_command *save = &cmd->save;
4547 struct string_array variables = STRING_ARRAY_INITIALIZER;
4548 struct matrix_expr *names = NULL;
4549 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
4551 save->expression = matrix_parse_exp (s);
4552 if (!save->expression)
4555 while (lex_match (s->lexer, T_SLASH))
4557 if (lex_match_id (s->lexer, "OUTFILE"))
4559 lex_match (s->lexer, T_EQUALS);
4562 fh = (lex_match (s->lexer, T_ASTERISK)
4564 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4568 else if (lex_match_id (s->lexer, "VARIABLES"))
4570 lex_match (s->lexer, T_EQUALS);
4574 struct dictionary *d = dict_create (get_default_encoding ());
4575 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4576 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4581 string_array_clear (&variables);
4582 variables = (struct string_array) {
4588 else if (lex_match_id (s->lexer, "NAMES"))
4590 lex_match (s->lexer, T_EQUALS);
4591 matrix_expr_destroy (names);
4592 names = matrix_parse_exp (s);
4596 else if (lex_match_id (s->lexer, "STRINGS"))
4598 lex_match (s->lexer, T_EQUALS);
4599 while (lex_token (s->lexer) == T_ID)
4601 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
4603 if (!lex_match (s->lexer, T_COMMA))
4609 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4617 if (s->prev_save_file)
4618 fh = fh_ref (s->prev_save_file);
4621 lex_sbc_missing ("OUTFILE");
4625 fh_unref (s->prev_save_file);
4626 s->prev_save_file = fh_ref (fh);
4628 if (variables.n && names)
4630 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4631 matrix_expr_destroy (names);
4635 save->sf = save_file_create (s, fh, &variables, names, &strings);
4641 string_array_destroy (&variables);
4642 matrix_expr_destroy (names);
4643 stringi_set_destroy (&strings);
4645 matrix_cmd_destroy (cmd);
4650 matrix_cmd_execute_save (const struct save_command *save)
4652 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4656 struct casewriter *writer = save_file_open (save->sf, m);
4659 gsl_matrix_free (m);
4663 const struct caseproto *proto = casewriter_get_proto (writer);
4664 for (size_t y = 0; y < m->size1; y++)
4666 struct ccase *c = case_create (proto);
4667 for (size_t x = 0; x < m->size2; x++)
4669 double d = gsl_matrix_get (m, y, x);
4670 int width = caseproto_get_width (proto, x);
4671 union value *value = case_data_rw_idx (c, x);
4675 memcpy (value->s, &d, width);
4677 casewriter_write (writer, c);
4679 gsl_matrix_free (m);
4682 static struct matrix_cmd *
4683 matrix_parse_read (struct matrix_state *s)
4685 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4686 *cmd = (struct matrix_cmd) {
4688 .read = { .format = FMT_F },
4691 struct file_handle *fh = NULL;
4692 char *encoding = NULL;
4693 struct read_command *read = &cmd->read;
4694 read->dst = matrix_lvalue_parse (s);
4699 int repetitions = 0;
4700 int record_width = 0;
4701 bool seen_format = false;
4702 while (lex_match (s->lexer, T_SLASH))
4704 if (lex_match_id (s->lexer, "FILE"))
4706 lex_match (s->lexer, T_EQUALS);
4709 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4713 else if (lex_match_id (s->lexer, "ENCODING"))
4715 lex_match (s->lexer, T_EQUALS);
4716 if (!lex_force_string (s->lexer))
4720 encoding = ss_xstrdup (lex_tokss (s->lexer));
4724 else if (lex_match_id (s->lexer, "FIELD"))
4726 lex_match (s->lexer, T_EQUALS);
4728 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4730 read->c1 = lex_integer (s->lexer);
4732 if (!lex_force_match (s->lexer, T_TO)
4733 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4735 read->c2 = lex_integer (s->lexer) + 1;
4738 record_width = read->c2 - read->c1;
4739 if (lex_match (s->lexer, T_BY))
4741 if (!lex_force_int_range (s->lexer, "BY", 1,
4742 read->c2 - read->c1))
4744 by = lex_integer (s->lexer);
4747 if (record_width % by)
4749 msg (SE, _("BY %d does not evenly divide record width %d."),
4757 else if (lex_match_id (s->lexer, "SIZE"))
4759 lex_match (s->lexer, T_EQUALS);
4760 matrix_expr_destroy (read->size);
4761 read->size = matrix_parse_exp (s);
4765 else if (lex_match_id (s->lexer, "MODE"))
4767 lex_match (s->lexer, T_EQUALS);
4768 if (lex_match_id (s->lexer, "RECTANGULAR"))
4769 read->symmetric = false;
4770 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4771 read->symmetric = true;
4774 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4778 else if (lex_match_id (s->lexer, "REREAD"))
4779 read->reread = true;
4780 else if (lex_match_id (s->lexer, "FORMAT"))
4784 lex_sbc_only_once ("FORMAT");
4789 lex_match (s->lexer, T_EQUALS);
4791 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4794 const char *p = lex_tokcstr (s->lexer);
4795 if (c_isdigit (p[0]))
4797 repetitions = atoi (p);
4798 p += strspn (p, "0123456789");
4799 if (!fmt_from_name (p, &read->format))
4801 lex_error (s->lexer, _("Unknown format %s."), p);
4806 else if (fmt_from_name (p, &read->format))
4810 struct fmt_spec format;
4811 if (!parse_format_specifier (s->lexer, &format))
4813 read->format = format.type;
4819 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4820 "REREAD", "FORMAT");
4827 lex_sbc_missing ("FIELD");
4831 if (!read->dst->n_indexes && !read->size)
4833 msg (SE, _("SIZE is required for reading data into a full matrix "
4834 "(as opposed to a submatrix)."));
4840 if (s->prev_read_file)
4841 fh = fh_ref (s->prev_read_file);
4844 lex_sbc_missing ("FILE");
4848 fh_unref (s->prev_read_file);
4849 s->prev_read_file = fh_ref (fh);
4851 read->rf = read_file_create (s, fh);
4855 free (read->rf->encoding);
4856 read->rf->encoding = encoding;
4860 /* Field width may be specified in multiple ways:
4863 2. The format on FORMAT.
4864 3. The repetition factor on FORMAT.
4866 (2) and (3) are mutually exclusive.
4868 If more than one of these is present, they must agree. If none of them is
4869 present, then free-field format is used.
4871 if (repetitions > record_width)
4873 msg (SE, _("%d repetitions cannot fit in record width %d."),
4874 repetitions, record_width);
4877 int w = (repetitions ? record_width / repetitions
4883 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4884 "which implies field width %d, "
4885 "but BY specifies field width %d."),
4886 repetitions, record_width, w, by);
4888 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4897 matrix_cmd_destroy (cmd);
4903 parse_error (const struct dfm_reader *reader, enum fmt_type format,
4904 struct substring data, size_t y, size_t x,
4905 int first_column, int last_column, char *error)
4907 int line_number = dfm_get_line_number (reader);
4908 struct msg_location *location = xmalloc (sizeof *location);
4909 *location = (struct msg_location) {
4910 .file_name = xstrdup (dfm_get_file_name (reader)),
4911 .first_line = line_number,
4912 .last_line = line_number + 1,
4913 .first_column = first_column,
4914 .last_column = last_column,
4916 struct msg *m = xmalloc (sizeof *m);
4918 .category = MSG_C_DATA,
4919 .severity = MSG_S_WARNING,
4920 .location = location,
4921 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
4922 "for matrix row %zu, column %zu: %s"),
4923 (int) data.length, data.string, fmt_name (format),
4924 y + 1, x + 1, error),
4932 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4933 gsl_matrix *m, struct substring p, size_t y, size_t x,
4934 const char *line_start)
4936 const char *input_encoding = dfm_reader_get_encoding (reader);
4939 if (fmt_is_numeric (read->format))
4942 error = data_in (p, input_encoding, read->format,
4943 settings_get_fmt_settings (), &v, 0, NULL);
4944 if (!error && v.f == SYSMIS)
4945 error = xstrdup (_("Matrix data may not contain missing value."));
4950 uint8_t s[sizeof (double)];
4951 union value v = { .s = s };
4952 error = data_in (p, input_encoding, read->format,
4953 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
4954 memcpy (&f, s, sizeof f);
4959 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
4960 int c2 = c1 + ss_utf8_count_columns (p) - 1;
4961 parse_error (reader, read->format, p, y, x, c1, c2, error);
4965 gsl_matrix_set (m, y, x, f);
4966 if (read->symmetric && x != y)
4967 gsl_matrix_set (m, x, y, f);
4972 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4973 struct substring *line, const char **startp)
4975 if (dfm_eof (reader))
4977 msg (SE, _("Unexpected end of file reading matrix data."));
4980 dfm_expand_tabs (reader);
4981 struct substring record = dfm_get_record (reader);
4982 /* XXX need to recode record into UTF-8 */
4983 *startp = record.string;
4984 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
4989 matrix_read (struct read_command *read, struct dfm_reader *reader,
4992 for (size_t y = 0; y < m->size1; y++)
4994 size_t nx = read->symmetric ? y + 1 : m->size2;
4996 struct substring line = ss_empty ();
4997 const char *line_start = line.string;
4998 for (size_t x = 0; x < nx; x++)
5005 ss_ltrim (&line, ss_cstr (" ,"));
5006 if (!ss_is_empty (line))
5008 if (!matrix_read_line (read, reader, &line, &line_start))
5010 dfm_forward_record (reader);
5013 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
5017 if (!matrix_read_line (read, reader, &line, &line_start))
5019 size_t fields_per_line = (read->c2 - read->c1) / read->w;
5020 int f = x % fields_per_line;
5021 if (f == fields_per_line - 1)
5022 dfm_forward_record (reader);
5024 p = ss_substr (line, read->w * f, read->w);
5027 matrix_read_set_field (read, reader, m, p, y, x, line_start);
5031 dfm_forward_record (reader);
5034 ss_ltrim (&line, ss_cstr (" ,"));
5035 if (!ss_is_empty (line))
5038 msg (SW, _("Trailing garbage on line \"%.*s\""),
5039 (int) line.length, line.string);
5046 matrix_cmd_execute_read (struct read_command *read)
5048 struct index_vector iv0, iv1;
5049 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
5052 size_t size[2] = { SIZE_MAX, SIZE_MAX };
5055 gsl_matrix *m = matrix_expr_evaluate (read->size);
5061 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5062 "not a %zu×%zu matrix."), m->size1, m->size2);
5063 gsl_matrix_free (m);
5069 gsl_vector v = to_vector (m);
5073 d[0] = gsl_vector_get (&v, 0);
5076 else if (v.size == 2)
5078 d[0] = gsl_vector_get (&v, 0);
5079 d[1] = gsl_vector_get (&v, 1);
5083 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5084 "not a %zu×%zu matrix."),
5085 m->size1, m->size2),
5086 gsl_matrix_free (m);
5091 gsl_matrix_free (m);
5093 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
5095 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
5096 "are outside valid range."),
5107 if (read->dst->n_indexes)
5109 size_t submatrix_size[2];
5110 if (read->dst->n_indexes == 2)
5112 submatrix_size[0] = iv0.n;
5113 submatrix_size[1] = iv1.n;
5115 else if (read->dst->var->value->size1 == 1)
5117 submatrix_size[0] = 1;
5118 submatrix_size[1] = iv0.n;
5122 submatrix_size[0] = iv0.n;
5123 submatrix_size[1] = 1;
5128 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
5130 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
5131 "differ from submatrix dimensions %zu×%zu."),
5133 submatrix_size[0], submatrix_size[1]);
5141 size[0] = submatrix_size[0];
5142 size[1] = submatrix_size[1];
5146 struct dfm_reader *reader = read_file_open (read->rf);
5148 dfm_reread_record (reader, 1);
5150 if (read->symmetric && size[0] != size[1])
5152 msg (SE, _("Cannot read non-square %zu×%zu matrix "
5153 "using READ with MODE=SYMMETRIC."),
5159 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
5160 matrix_read (read, reader, tmp);
5161 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
5164 static struct matrix_cmd *
5165 matrix_parse_write (struct matrix_state *s)
5167 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5168 *cmd = (struct matrix_cmd) {
5172 struct file_handle *fh = NULL;
5173 char *encoding = NULL;
5174 struct write_command *write = &cmd->write;
5175 write->expression = matrix_parse_exp (s);
5176 if (!write->expression)
5180 int repetitions = 0;
5181 int record_width = 0;
5182 enum fmt_type format = FMT_F;
5183 bool has_format = false;
5184 while (lex_match (s->lexer, T_SLASH))
5186 if (lex_match_id (s->lexer, "OUTFILE"))
5188 lex_match (s->lexer, T_EQUALS);
5191 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5195 else if (lex_match_id (s->lexer, "ENCODING"))
5197 lex_match (s->lexer, T_EQUALS);
5198 if (!lex_force_string (s->lexer))
5202 encoding = ss_xstrdup (lex_tokss (s->lexer));
5206 else if (lex_match_id (s->lexer, "FIELD"))
5208 lex_match (s->lexer, T_EQUALS);
5210 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5212 write->c1 = lex_integer (s->lexer);
5214 if (!lex_force_match (s->lexer, T_TO)
5215 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5217 write->c2 = lex_integer (s->lexer) + 1;
5220 record_width = write->c2 - write->c1;
5221 if (lex_match (s->lexer, T_BY))
5223 if (!lex_force_int_range (s->lexer, "BY", 1,
5224 write->c2 - write->c1))
5226 by = lex_integer (s->lexer);
5229 if (record_width % by)
5231 msg (SE, _("BY %d does not evenly divide record width %d."),
5239 else if (lex_match_id (s->lexer, "MODE"))
5241 lex_match (s->lexer, T_EQUALS);
5242 if (lex_match_id (s->lexer, "RECTANGULAR"))
5243 write->triangular = false;
5244 else if (lex_match_id (s->lexer, "TRIANGULAR"))
5245 write->triangular = true;
5248 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5252 else if (lex_match_id (s->lexer, "HOLD"))
5254 else if (lex_match_id (s->lexer, "FORMAT"))
5256 if (has_format || write->format)
5258 lex_sbc_only_once ("FORMAT");
5262 lex_match (s->lexer, T_EQUALS);
5264 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5267 const char *p = lex_tokcstr (s->lexer);
5268 if (c_isdigit (p[0]))
5270 repetitions = atoi (p);
5271 p += strspn (p, "0123456789");
5272 if (!fmt_from_name (p, &format))
5274 lex_error (s->lexer, _("Unknown format %s."), p);
5280 else if (fmt_from_name (p, &format))
5287 struct fmt_spec spec;
5288 if (!parse_format_specifier (s->lexer, &spec))
5290 write->format = xmemdup (&spec, sizeof spec);
5295 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5303 lex_sbc_missing ("FIELD");
5309 if (s->prev_write_file)
5310 fh = fh_ref (s->prev_write_file);
5313 lex_sbc_missing ("OUTFILE");
5317 fh_unref (s->prev_write_file);
5318 s->prev_write_file = fh_ref (fh);
5320 write->wf = write_file_create (s, fh);
5324 free (write->wf->encoding);
5325 write->wf->encoding = encoding;
5329 /* Field width may be specified in multiple ways:
5332 2. The format on FORMAT.
5333 3. The repetition factor on FORMAT.
5335 (2) and (3) are mutually exclusive.
5337 If more than one of these is present, they must agree. If none of them is
5338 present, then free-field format is used.
5340 if (repetitions > record_width)
5342 msg (SE, _("%d repetitions cannot fit in record width %d."),
5343 repetitions, record_width);
5346 int w = (repetitions ? record_width / repetitions
5347 : write->format ? write->format->w
5352 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5353 "which implies field width %d, "
5354 "but BY specifies field width %d."),
5355 repetitions, record_width, w, by);
5357 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5361 if (w && !write->format)
5363 write->format = xmalloc (sizeof *write->format);
5364 *write->format = (struct fmt_spec) { .type = format, .w = w };
5366 if (!fmt_check_output (write->format))
5370 if (write->format && fmt_var_width (write->format) > sizeof (double))
5372 char s[FMT_STRING_LEN_MAX + 1];
5373 fmt_to_string (write->format, s);
5374 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
5375 s, sizeof (double));
5383 matrix_cmd_destroy (cmd);
5388 matrix_cmd_execute_write (struct write_command *write)
5390 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5394 if (write->triangular && m->size1 != m->size2)
5396 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5397 "the matrix to be written has dimensions %zu×%zu."),
5398 m->size1, m->size2);
5399 gsl_matrix_free (m);
5403 struct dfm_writer *writer = write_file_open (write->wf);
5404 if (!writer || !m->size1)
5406 gsl_matrix_free (m);
5410 const struct fmt_settings *settings = settings_get_fmt_settings ();
5411 struct u8_line *line = write->wf->held;
5412 for (size_t y = 0; y < m->size1; y++)
5416 line = xmalloc (sizeof *line);
5417 u8_line_init (line);
5419 size_t nx = write->triangular ? y + 1 : m->size2;
5421 for (size_t x = 0; x < nx; x++)
5424 double f = gsl_matrix_get (m, y, x);
5428 if (fmt_is_numeric (write->format->type))
5431 v.s = (uint8_t *) &f;
5432 s = data_out (&v, NULL, write->format, settings);
5436 s = xmalloc (DBL_BUFSIZE_BOUND);
5437 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
5438 >= DBL_BUFSIZE_BOUND)
5441 size_t len = strlen (s);
5442 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5443 if (width + x0 > write->c2)
5445 dfm_put_record_utf8 (writer, line->s.ss.string,
5447 u8_line_clear (line);
5450 u8_line_put (line, x0, x0 + width, s, len);
5453 x0 += write->format ? write->format->w : width + 1;
5456 if (y + 1 >= m->size1 && write->hold)
5458 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5459 u8_line_clear (line);
5463 u8_line_destroy (line);
5467 write->wf->held = line;
5469 gsl_matrix_free (m);
5472 static struct matrix_cmd *
5473 matrix_parse_get (struct matrix_state *s)
5475 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5476 *cmd = (struct matrix_cmd) {
5479 .dataset = s->dataset,
5480 .user = { .treatment = MGET_ERROR },
5481 .system = { .treatment = MGET_ERROR },
5485 struct get_command *get = &cmd->get;
5486 get->dst = matrix_lvalue_parse (s);
5490 while (lex_match (s->lexer, T_SLASH))
5492 if (lex_match_id (s->lexer, "FILE"))
5494 lex_match (s->lexer, T_EQUALS);
5496 fh_unref (get->file);
5497 if (lex_match (s->lexer, T_ASTERISK))
5501 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5506 else if (lex_match_id (s->lexer, "ENCODING"))
5508 lex_match (s->lexer, T_EQUALS);
5509 if (!lex_force_string (s->lexer))
5512 free (get->encoding);
5513 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5517 else if (lex_match_id (s->lexer, "VARIABLES"))
5519 lex_match (s->lexer, T_EQUALS);
5523 lex_sbc_only_once ("VARIABLES");
5527 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
5530 else if (lex_match_id (s->lexer, "NAMES"))
5532 lex_match (s->lexer, T_EQUALS);
5533 if (!lex_force_id (s->lexer))
5536 struct substring name = lex_tokss (s->lexer);
5537 get->names = matrix_var_lookup (s, name);
5539 get->names = matrix_var_create (s, name);
5542 else if (lex_match_id (s->lexer, "MISSING"))
5544 lex_match (s->lexer, T_EQUALS);
5545 if (lex_match_id (s->lexer, "ACCEPT"))
5546 get->user.treatment = MGET_ACCEPT;
5547 else if (lex_match_id (s->lexer, "OMIT"))
5548 get->user.treatment = MGET_OMIT;
5549 else if (lex_is_number (s->lexer))
5551 get->user.treatment = MGET_RECODE;
5552 get->user.substitute = lex_number (s->lexer);
5557 lex_error (s->lexer, NULL);
5561 else if (lex_match_id (s->lexer, "SYSMIS"))
5563 lex_match (s->lexer, T_EQUALS);
5564 if (lex_match_id (s->lexer, "OMIT"))
5565 get->system.treatment = MGET_OMIT;
5566 else if (lex_is_number (s->lexer))
5568 get->system.treatment = MGET_RECODE;
5569 get->system.substitute = lex_number (s->lexer);
5574 lex_error (s->lexer, NULL);
5580 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5581 "MISSING", "SYSMIS");
5586 if (get->user.treatment != MGET_ACCEPT)
5587 get->system.treatment = MGET_ERROR;
5592 matrix_cmd_destroy (cmd);
5597 matrix_cmd_execute_get__ (struct get_command *get,
5598 const struct dictionary *dict,
5599 struct casereader *reader)
5601 struct variable **vars;
5606 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
5607 &vars, &n_vars, PV_NUMERIC))
5612 n_vars = dict_get_var_cnt (dict);
5613 vars = xnmalloc (n_vars, sizeof *vars);
5614 for (size_t i = 0; i < n_vars; i++)
5616 struct variable *var = dict_get_var (dict, i);
5617 if (!var_is_numeric (var))
5619 msg (SE, _("GET: Variable %s is not numeric."),
5620 var_get_name (var));
5630 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
5631 for (size_t i = 0; i < n_vars; i++)
5633 char s[sizeof (double)];
5635 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
5636 memcpy (&f, s, sizeof f);
5637 gsl_matrix_set (names, i, 0, f);
5640 gsl_matrix_free (get->names->value);
5641 get->names->value = names;
5645 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5646 long long int casenum = 1;
5648 for (struct ccase *c = casereader_read (reader); c;
5649 c = casereader_read (reader), casenum++)
5651 if (n_rows >= m->size1)
5653 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5654 for (size_t y = 0; y < n_rows; y++)
5655 for (size_t x = 0; x < n_vars; x++)
5656 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5657 gsl_matrix_free (m);
5662 for (size_t x = 0; x < n_vars; x++)
5664 const struct variable *var = vars[x];
5665 double d = case_num (c, var);
5668 if (get->system.treatment == MGET_RECODE)
5669 d = get->system.substitute;
5670 else if (get->system.treatment == MGET_OMIT)
5674 msg (SE, _("GET: Variable %s in case %lld "
5675 "is system-missing."),
5676 var_get_name (var), casenum);
5680 else if (var_is_num_missing (var, d, MV_USER))
5682 if (get->user.treatment == MGET_RECODE)
5683 d = get->user.substitute;
5684 else if (get->user.treatment == MGET_OMIT)
5686 else if (get->user.treatment != MGET_ACCEPT)
5688 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5690 var_get_name (var), casenum, d);
5694 gsl_matrix_set (m, n_rows, x, d);
5705 matrix_lvalue_evaluate_and_assign (get->dst, m);
5708 gsl_matrix_free (m);
5713 matrix_cmd_execute_get (struct get_command *get)
5715 struct dictionary *dict;
5716 struct casereader *reader;
5719 reader = any_reader_open_and_decode (get->file, get->encoding,
5726 if (dict_get_var_cnt (dataset_dict (get->dataset)) == 0)
5728 msg (ME, _("GET cannot read empty active file."));
5731 reader = proc_open (get->dataset);
5732 dict = dict_ref (dataset_dict (get->dataset));
5735 matrix_cmd_execute_get__ (get, dict, reader);
5738 casereader_destroy (reader);
5740 proc_commit (get->dataset);
5744 match_rowtype (struct lexer *lexer)
5746 static const char *rowtypes[] = {
5747 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5749 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5751 for (size_t i = 0; i < n_rowtypes; i++)
5752 if (lex_match_id (lexer, rowtypes[i]))
5755 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5760 parse_var_names (struct lexer *lexer, struct string_array *sa)
5762 lex_match (lexer, T_EQUALS);
5764 string_array_clear (sa);
5766 struct dictionary *dict = dict_create (get_default_encoding ());
5767 dict_create_var_assert (dict, "ROWTYPE_", 8);
5768 dict_create_var_assert (dict, "VARNAME_", 8);
5771 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5772 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5777 string_array_clear (sa);
5778 sa->strings = names;
5779 sa->n = sa->allocated = n_names;
5785 msave_common_uninit (struct msave_common *common)
5789 fh_unref (common->outfile);
5790 string_array_destroy (&common->variables);
5791 string_array_destroy (&common->fnames);
5792 string_array_destroy (&common->snames);
5797 compare_variables (const char *keyword,
5798 const struct string_array *new,
5799 const struct string_array *old)
5805 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5806 "on the first MSAVE within MATRIX."), keyword);
5809 else if (!string_array_equal_case (old, new))
5811 msg (SE, _("%s must specify the same variables each time within "
5812 "a given MATRIX."), keyword);
5819 static struct matrix_cmd *
5820 matrix_parse_msave (struct matrix_state *s)
5822 struct msave_common common = { .outfile = NULL };
5823 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5824 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5826 struct msave_command *msave = &cmd->msave;
5827 if (lex_token (s->lexer) == T_ID)
5828 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5829 msave->expr = matrix_parse_exp (s);
5833 while (lex_match (s->lexer, T_SLASH))
5835 if (lex_match_id (s->lexer, "TYPE"))
5837 lex_match (s->lexer, T_EQUALS);
5839 msave->rowtype = match_rowtype (s->lexer);
5840 if (!msave->rowtype)
5843 else if (lex_match_id (s->lexer, "OUTFILE"))
5845 lex_match (s->lexer, T_EQUALS);
5847 fh_unref (common.outfile);
5848 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5849 if (!common.outfile)
5852 else if (lex_match_id (s->lexer, "VARIABLES"))
5854 if (!parse_var_names (s->lexer, &common.variables))
5857 else if (lex_match_id (s->lexer, "FNAMES"))
5859 if (!parse_var_names (s->lexer, &common.fnames))
5862 else if (lex_match_id (s->lexer, "SNAMES"))
5864 if (!parse_var_names (s->lexer, &common.snames))
5867 else if (lex_match_id (s->lexer, "SPLIT"))
5869 lex_match (s->lexer, T_EQUALS);
5871 matrix_expr_destroy (msave->splits);
5872 msave->splits = matrix_parse_exp (s);
5876 else if (lex_match_id (s->lexer, "FACTOR"))
5878 lex_match (s->lexer, T_EQUALS);
5880 matrix_expr_destroy (msave->factors);
5881 msave->factors = matrix_parse_exp (s);
5882 if (!msave->factors)
5887 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5888 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5892 if (!msave->rowtype)
5894 lex_sbc_missing ("TYPE");
5897 common.has_splits = msave->splits || common.snames.n;
5898 common.has_factors = msave->factors || common.fnames.n;
5900 struct msave_common *c = s->common ? s->common : &common;
5901 if (c->fnames.n && !msave->factors)
5903 msg (SE, _("FNAMES requires FACTOR."));
5906 if (c->snames.n && !msave->splits)
5908 msg (SE, _("SNAMES requires SPLIT."));
5911 if (c->has_factors && !common.has_factors)
5913 msg (SE, _("%s is required because it was present on the first "
5914 "MSAVE in this MATRIX command."), "FACTOR");
5917 if (c->has_splits && !common.has_splits)
5919 msg (SE, _("%s is required because it was present on the first "
5920 "MSAVE in this MATRIX command."), "SPLIT");
5926 if (!common.outfile)
5928 lex_sbc_missing ("OUTFILE");
5931 s->common = xmemdup (&common, sizeof common);
5937 bool same = common.outfile == s->common->outfile;
5938 fh_unref (common.outfile);
5941 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5942 "within a single MATRIX command."));
5946 if (!compare_variables ("VARIABLES",
5947 &common.variables, &s->common->variables)
5948 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5949 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5951 msave_common_uninit (&common);
5953 msave->common = s->common;
5954 if (!msave->varname_)
5955 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5959 msave_common_uninit (&common);
5960 matrix_cmd_destroy (cmd);
5965 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5967 gsl_matrix *m = matrix_expr_evaluate (e);
5973 msg (SE, _("%s expression must evaluate to vector, "
5974 "not a %zu×%zu matrix."),
5975 name, m->size1, m->size2);
5976 gsl_matrix_free (m);
5980 return matrix_to_vector (m);
5984 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5986 for (size_t i = 0; i < vars->n; i++)
5987 if (!dict_create_var (d, vars->strings[i], 0))
5988 return vars->strings[i];
5992 static struct dictionary *
5993 msave_create_dict (const struct msave_common *common)
5995 struct dictionary *dict = dict_create (get_default_encoding ());
5997 const char *dup_split = msave_add_vars (dict, &common->snames);
6000 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
6004 dict_create_var_assert (dict, "ROWTYPE_", 8);
6006 const char *dup_factor = msave_add_vars (dict, &common->fnames);
6009 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
6013 dict_create_var_assert (dict, "VARNAME_", 8);
6015 const char *dup_var = msave_add_vars (dict, &common->variables);
6018 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
6030 matrix_cmd_execute_msave (struct msave_command *msave)
6032 struct msave_common *common = msave->common;
6033 gsl_matrix *m = NULL;
6034 gsl_vector *factors = NULL;
6035 gsl_vector *splits = NULL;
6037 m = matrix_expr_evaluate (msave->expr);
6041 if (!common->variables.n)
6042 for (size_t i = 0; i < m->size2; i++)
6043 string_array_append_nocopy (&common->variables,
6044 xasprintf ("COL%zu", i + 1));
6046 if (m->size2 != common->variables.n)
6048 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
6049 m->size2, common->variables.n);
6055 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
6059 if (!common->fnames.n)
6060 for (size_t i = 0; i < factors->size; i++)
6061 string_array_append_nocopy (&common->fnames,
6062 xasprintf ("FAC%zu", i + 1));
6064 if (factors->size != common->fnames.n)
6066 msg (SE, _("There are %zu factor variables, "
6067 "but %zu split values were supplied."),
6068 common->fnames.n, factors->size);
6075 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
6079 if (!common->fnames.n)
6080 for (size_t i = 0; i < splits->size; i++)
6081 string_array_append_nocopy (&common->fnames,
6082 xasprintf ("SPL%zu", i + 1));
6084 if (splits->size != common->fnames.n)
6086 msg (SE, _("There are %zu split variables, "
6087 "but %zu split values were supplied."),
6088 common->fnames.n, splits->size);
6093 if (!common->writer)
6095 struct dictionary *dict = msave_create_dict (common);
6099 common->writer = any_writer_open (common->outfile, dict);
6100 if (!common->writer)
6106 common->dict = dict;
6109 for (size_t y = 0; y < m->size1; y++)
6111 struct ccase *c = case_create (dict_get_proto (common->dict));
6114 /* Split variables */
6116 for (size_t i = 0; i < splits->size; i++)
6117 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
6120 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6121 msave->rowtype, ' ');
6125 for (size_t i = 0; i < factors->size; i++)
6126 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
6129 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6130 msave->varname_, ' ');
6132 /* Continuous variables. */
6133 for (size_t x = 0; x < m->size2; x++)
6134 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
6135 casewriter_write (common->writer, c);
6139 gsl_matrix_free (m);
6140 gsl_vector_free (factors);
6141 gsl_vector_free (splits);
6144 static struct matrix_cmd *
6145 matrix_parse_mget (struct matrix_state *s)
6147 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6148 *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
6150 struct mget_command *mget = &cmd->mget;
6154 if (lex_match_id (s->lexer, "FILE"))
6156 lex_match (s->lexer, T_EQUALS);
6158 fh_unref (mget->file);
6159 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6163 else if (lex_match_id (s->lexer, "TYPE"))
6165 lex_match (s->lexer, T_EQUALS);
6166 while (lex_token (s->lexer) != T_SLASH
6167 && lex_token (s->lexer) != T_ENDCMD)
6169 const char *rowtype = match_rowtype (s->lexer);
6173 stringi_set_insert (&mget->rowtypes, rowtype);
6178 lex_error_expecting (s->lexer, "FILE", "TYPE");
6181 if (lex_token (s->lexer) == T_ENDCMD)
6184 if (!lex_force_match (s->lexer, T_SLASH))
6190 matrix_cmd_destroy (cmd);
6194 static const struct variable *
6195 get_a8_var (const struct dictionary *d, const char *name)
6197 const struct variable *v = dict_lookup_var (d, name);
6200 msg (SE, _("Matrix data file lacks %s variable."), name);
6203 if (var_get_width (v) != 8)
6205 msg (SE, _("%s variable in matrix data file must be "
6206 "8-byte string, but it has width %d."),
6207 name, var_get_width (v));
6214 is_rowtype (const union value *v, const char *rowtype)
6216 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
6217 ss_rtrim (&vs, ss_cstr (" "));
6218 return ss_equals_case (vs, ss_cstr (rowtype));
6222 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6223 const struct dictionary *d,
6224 const struct variable *rowtype_var,
6225 struct matrix_state *s, size_t si, size_t fi,
6226 size_t cs, size_t cn)
6231 const union value *rowtype_ = case_data (rows[0], rowtype_var);
6232 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
6233 : is_rowtype (rowtype_, "CORR") ? "CR"
6234 : is_rowtype (rowtype_, "MEAN") ? "MN"
6235 : is_rowtype (rowtype_, "STDDEV") ? "SD"
6236 : is_rowtype (rowtype_, "N") ? "NC"
6239 struct string name = DS_EMPTY_INITIALIZER;
6240 ds_put_cstr (&name, name_prefix);
6242 ds_put_format (&name, "F%zu", fi);
6244 ds_put_format (&name, "S%zu", si);
6246 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6248 mv = matrix_var_create (s, ds_ss (&name));
6251 msg (SW, _("Matrix data file contains variable with existing name %s."),
6256 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6257 size_t n_missing = 0;
6258 for (size_t y = 0; y < n_rows; y++)
6260 for (size_t x = 0; x < cn; x++)
6262 struct variable *var = dict_get_var (d, cs + x);
6263 double value = case_num (rows[y], var);
6264 if (var_is_num_missing (var, value, MV_ANY))
6269 gsl_matrix_set (m, y, x, value);
6274 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6275 "value, which was treated as zero.",
6276 "Matrix data file variable %s contains %zu missing "
6277 "values, which were treated as zero.", n_missing),
6278 ds_cstr (&name), n_missing);
6283 for (size_t y = 0; y < n_rows; y++)
6284 case_unref (rows[y]);
6288 var_changed (const struct ccase *ca, const struct ccase *cb,
6289 const struct variable *var)
6291 return !value_equal (case_data (ca, var), case_data (cb, var),
6292 var_get_width (var));
6296 vars_changed (const struct ccase *ca, const struct ccase *cb,
6297 const struct dictionary *d,
6298 size_t first_var, size_t n_vars)
6300 for (size_t i = 0; i < n_vars; i++)
6302 const struct variable *v = dict_get_var (d, first_var + i);
6303 if (var_changed (ca, cb, v))
6310 matrix_cmd_execute_mget (struct mget_command *mget)
6312 struct dictionary *d;
6313 struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6318 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6319 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6320 if (!rowtype_ || !varname_)
6323 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6325 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6328 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6330 msg (SE, _("Matrix data file contains no continuous variables."));
6334 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6336 const struct variable *v = dict_get_var (d, i);
6337 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6340 _("Matrix data file contains unexpected string variable %s."),
6346 /* SPLIT variables. */
6348 size_t sn = var_get_dict_index (rowtype_);
6349 struct ccase *sc = NULL;
6352 /* FACTOR variables. */
6353 size_t fs = var_get_dict_index (rowtype_) + 1;
6354 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6355 struct ccase *fc = NULL;
6358 /* Continuous variables. */
6359 size_t cs = var_get_dict_index (varname_) + 1;
6360 size_t cn = dict_get_var_cnt (d) - cs;
6361 struct ccase *cc = NULL;
6364 struct ccase **rows = NULL;
6365 size_t allocated_rows = 0;
6369 while ((c = casereader_read (r)) != NULL)
6371 bool sd = vars_changed (sc, c, d, ss, sn);
6372 bool fd = sd || vars_changed (fc, c, d, fs, fn);
6373 bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6388 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6389 mget->state, si, fi, cs, cn);
6395 if (n_rows >= allocated_rows)
6396 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6399 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6400 mget->state, si, fi, cs, cn);
6404 casereader_destroy (r);
6408 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6410 if (!lex_force_id (s->lexer))
6413 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6415 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6420 static struct matrix_cmd *
6421 matrix_parse_eigen (struct matrix_state *s)
6423 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6424 *cmd = (struct matrix_cmd) {
6426 .eigen = { .expr = NULL }
6429 struct eigen_command *eigen = &cmd->eigen;
6430 if (!lex_force_match (s->lexer, T_LPAREN))
6432 eigen->expr = matrix_parse_expr (s);
6434 || !lex_force_match (s->lexer, T_COMMA)
6435 || !matrix_parse_dst_var (s, &eigen->evec)
6436 || !lex_force_match (s->lexer, T_COMMA)
6437 || !matrix_parse_dst_var (s, &eigen->eval)
6438 || !lex_force_match (s->lexer, T_RPAREN))
6444 matrix_cmd_destroy (cmd);
6449 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6451 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6454 if (!is_symmetric (A))
6456 msg (SE, _("Argument of EIGEN must be symmetric."));
6457 gsl_matrix_free (A);
6461 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6462 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6463 gsl_vector v_eval = to_vector (eval);
6464 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6465 gsl_eigen_symmv (A, &v_eval, evec, w);
6466 gsl_eigen_symmv_free (w);
6468 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6470 gsl_matrix_free (A);
6472 gsl_matrix_free (eigen->eval->value);
6473 eigen->eval->value = eval;
6475 gsl_matrix_free (eigen->evec->value);
6476 eigen->evec->value = evec;
6479 static struct matrix_cmd *
6480 matrix_parse_setdiag (struct matrix_state *s)
6482 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6483 *cmd = (struct matrix_cmd) {
6484 .type = MCMD_SETDIAG,
6485 .setdiag = { .dst = NULL }
6488 struct setdiag_command *setdiag = &cmd->setdiag;
6489 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6492 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6495 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6500 if (!lex_force_match (s->lexer, T_COMMA))
6503 setdiag->expr = matrix_parse_expr (s);
6507 if (!lex_force_match (s->lexer, T_RPAREN))
6513 matrix_cmd_destroy (cmd);
6518 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6520 gsl_matrix *dst = setdiag->dst->value;
6523 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6524 setdiag->dst->name);
6528 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6532 size_t n = MIN (dst->size1, dst->size2);
6533 if (is_scalar (src))
6535 double d = to_scalar (src);
6536 for (size_t i = 0; i < n; i++)
6537 gsl_matrix_set (dst, i, i, d);
6539 else if (is_vector (src))
6541 gsl_vector v = to_vector (src);
6542 for (size_t i = 0; i < n && i < v.size; i++)
6543 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6546 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6547 "not a %zu×%zu matrix."),
6548 src->size1, src->size2);
6549 gsl_matrix_free (src);
6552 static struct matrix_cmd *
6553 matrix_parse_svd (struct matrix_state *s)
6555 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6556 *cmd = (struct matrix_cmd) {
6558 .svd = { .expr = NULL }
6561 struct svd_command *svd = &cmd->svd;
6562 if (!lex_force_match (s->lexer, T_LPAREN))
6564 svd->expr = matrix_parse_expr (s);
6566 || !lex_force_match (s->lexer, T_COMMA)
6567 || !matrix_parse_dst_var (s, &svd->u)
6568 || !lex_force_match (s->lexer, T_COMMA)
6569 || !matrix_parse_dst_var (s, &svd->s)
6570 || !lex_force_match (s->lexer, T_COMMA)
6571 || !matrix_parse_dst_var (s, &svd->v)
6572 || !lex_force_match (s->lexer, T_RPAREN))
6578 matrix_cmd_destroy (cmd);
6583 matrix_cmd_execute_svd (struct svd_command *svd)
6585 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6589 if (m->size1 >= m->size2)
6592 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6593 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6594 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6595 gsl_vector *work = gsl_vector_alloc (A->size2);
6596 gsl_linalg_SV_decomp (A, V, &Sv, work);
6597 gsl_vector_free (work);
6599 matrix_var_set (svd->u, A);
6600 matrix_var_set (svd->s, S);
6601 matrix_var_set (svd->v, V);
6605 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6606 gsl_matrix_transpose_memcpy (At, m);
6607 gsl_matrix_free (m);
6609 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6610 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6611 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6612 gsl_vector *work = gsl_vector_alloc (At->size2);
6613 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6614 gsl_vector_free (work);
6616 matrix_var_set (svd->v, At);
6617 matrix_var_set (svd->s, St);
6618 matrix_var_set (svd->u, Vt);
6623 matrix_cmd_execute (struct matrix_cmd *cmd)
6628 matrix_cmd_execute_compute (&cmd->compute);
6632 matrix_cmd_execute_print (&cmd->print);
6636 return matrix_cmd_execute_do_if (&cmd->do_if);
6639 matrix_cmd_execute_loop (&cmd->loop);
6646 matrix_cmd_execute_display (&cmd->display);
6650 matrix_cmd_execute_release (&cmd->release);
6654 matrix_cmd_execute_save (&cmd->save);
6658 matrix_cmd_execute_read (&cmd->read);
6662 matrix_cmd_execute_write (&cmd->write);
6666 matrix_cmd_execute_get (&cmd->get);
6670 matrix_cmd_execute_msave (&cmd->msave);
6674 matrix_cmd_execute_mget (&cmd->mget);
6678 matrix_cmd_execute_eigen (&cmd->eigen);
6682 matrix_cmd_execute_setdiag (&cmd->setdiag);
6686 matrix_cmd_execute_svd (&cmd->svd);
6694 matrix_cmds_uninit (struct matrix_cmds *cmds)
6696 for (size_t i = 0; i < cmds->n; i++)
6697 matrix_cmd_destroy (cmds->commands[i]);
6698 free (cmds->commands);
6702 matrix_cmd_destroy (struct matrix_cmd *cmd)
6710 matrix_lvalue_destroy (cmd->compute.lvalue);
6711 matrix_expr_destroy (cmd->compute.rvalue);
6715 matrix_expr_destroy (cmd->print.expression);
6716 free (cmd->print.title);
6717 print_labels_destroy (cmd->print.rlabels);
6718 print_labels_destroy (cmd->print.clabels);
6722 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6724 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6725 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6727 free (cmd->do_if.clauses);
6731 matrix_expr_destroy (cmd->loop.start);
6732 matrix_expr_destroy (cmd->loop.end);
6733 matrix_expr_destroy (cmd->loop.increment);
6734 matrix_expr_destroy (cmd->loop.top_condition);
6735 matrix_expr_destroy (cmd->loop.bottom_condition);
6736 matrix_cmds_uninit (&cmd->loop.commands);
6746 free (cmd->release.vars);
6750 matrix_expr_destroy (cmd->save.expression);
6754 matrix_lvalue_destroy (cmd->read.dst);
6755 matrix_expr_destroy (cmd->read.size);
6759 matrix_expr_destroy (cmd->write.expression);
6760 free (cmd->write.format);
6764 matrix_lvalue_destroy (cmd->get.dst);
6765 fh_unref (cmd->get.file);
6766 free (cmd->get.encoding);
6767 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
6771 free (cmd->msave.varname_);
6772 matrix_expr_destroy (cmd->msave.expr);
6773 matrix_expr_destroy (cmd->msave.factors);
6774 matrix_expr_destroy (cmd->msave.splits);
6778 fh_unref (cmd->mget.file);
6779 stringi_set_destroy (&cmd->mget.rowtypes);
6783 matrix_expr_destroy (cmd->eigen.expr);
6787 matrix_expr_destroy (cmd->setdiag.expr);
6791 matrix_expr_destroy (cmd->svd.expr);
6797 struct matrix_command_name
6800 struct matrix_cmd *(*parse) (struct matrix_state *);
6803 static const struct matrix_command_name *
6804 matrix_parse_command_name (struct lexer *lexer)
6806 static const struct matrix_command_name commands[] = {
6807 { "COMPUTE", matrix_parse_compute },
6808 { "CALL EIGEN", matrix_parse_eigen },
6809 { "CALL SETDIAG", matrix_parse_setdiag },
6810 { "CALL SVD", matrix_parse_svd },
6811 { "PRINT", matrix_parse_print },
6812 { "DO IF", matrix_parse_do_if },
6813 { "LOOP", matrix_parse_loop },
6814 { "BREAK", matrix_parse_break },
6815 { "READ", matrix_parse_read },
6816 { "WRITE", matrix_parse_write },
6817 { "GET", matrix_parse_get },
6818 { "SAVE", matrix_parse_save },
6819 { "MGET", matrix_parse_mget },
6820 { "MSAVE", matrix_parse_msave },
6821 { "DISPLAY", matrix_parse_display },
6822 { "RELEASE", matrix_parse_release },
6824 static size_t n = sizeof commands / sizeof *commands;
6826 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6827 if (lex_match_phrase (lexer, c->name))
6832 static struct matrix_cmd *
6833 matrix_parse_command (struct matrix_state *s)
6835 size_t nesting_level = SIZE_MAX;
6837 struct matrix_cmd *c = NULL;
6838 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6840 lex_error (s->lexer, _("Unknown matrix command."));
6841 else if (!cmd->parse)
6842 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6846 nesting_level = output_open_group (
6847 group_item_create_nocopy (utf8_to_title (cmd->name),
6848 utf8_to_title (cmd->name)));
6853 lex_end_of_command (s->lexer);
6854 lex_discard_rest_of_command (s->lexer);
6855 if (nesting_level != SIZE_MAX)
6856 output_close_groups (nesting_level);
6862 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6864 if (!lex_force_match (lexer, T_ENDCMD))
6867 struct matrix_state state = {
6869 .session = dataset_session (ds),
6871 .vars = HMAP_INITIALIZER (state.vars),
6876 while (lex_match (lexer, T_ENDCMD))
6878 if (lex_token (lexer) == T_STOP)
6880 msg (SE, _("Unexpected end of input expecting matrix command."));
6884 if (lex_match_phrase (lexer, "END MATRIX"))
6887 struct matrix_cmd *c = matrix_parse_command (&state);
6890 matrix_cmd_execute (c);
6891 matrix_cmd_destroy (c);
6895 struct matrix_var *var, *next;
6896 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6899 gsl_matrix_free (var->value);
6900 hmap_delete (&state.vars, &var->hmap_node);
6903 hmap_destroy (&state.vars);
6906 dict_unref (state.common->dict);
6907 casewriter_destroy (state.common->writer);
6908 free (state.common);
6910 fh_unref (state.prev_read_file);
6911 for (size_t i = 0; i < state.n_read_files; i++)
6912 read_file_destroy (state.read_files[i]);
6913 free (state.read_files);
6914 fh_unref (state.prev_write_file);
6915 for (size_t i = 0; i < state.n_write_files; i++)
6916 write_file_destroy (state.write_files[i]);
6917 free (state.write_files);
6918 fh_unref (state.prev_save_file);
6919 for (size_t i = 0; i < state.n_save_files; i++)
6920 save_file_destroy (state.save_files[i]);
6921 free (state.save_files);