1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2021 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include <gsl/gsl_blas.h>
20 #include <gsl/gsl_eigen.h>
21 #include <gsl/gsl_linalg.h>
22 #include <gsl/gsl_matrix.h>
23 #include <gsl/gsl_permutation.h>
24 #include <gsl/gsl_randist.h>
25 #include <gsl/gsl_vector.h>
30 #include "data/any-reader.h"
31 #include "data/any-writer.h"
32 #include "data/casereader.h"
33 #include "data/casewriter.h"
34 #include "data/data-in.h"
35 #include "data/data-out.h"
36 #include "data/dataset.h"
37 #include "data/dictionary.h"
38 #include "data/file-handle-def.h"
39 #include "language/command.h"
40 #include "language/data-io/data-reader.h"
41 #include "language/data-io/data-writer.h"
42 #include "language/data-io/file-handle.h"
43 #include "language/lexer/format-parser.h"
44 #include "language/lexer/lexer.h"
45 #include "language/lexer/variable-parser.h"
46 #include "libpspp/array.h"
47 #include "libpspp/assertion.h"
48 #include "libpspp/compiler.h"
49 #include "libpspp/hmap.h"
50 #include "libpspp/i18n.h"
51 #include "libpspp/misc.h"
52 #include "libpspp/str.h"
53 #include "libpspp/string-array.h"
54 #include "libpspp/stringi-set.h"
55 #include "libpspp/u8-line.h"
56 #include "math/random.h"
57 #include "output/driver.h"
58 #include "output/output-item.h"
59 #include "output/pivot-table.h"
61 #include "gl/c-ctype.h"
62 #include "gl/ftoastr.h"
63 #include "gl/minmax.h"
67 #define _(msgid) gettext (msgid)
68 #define N_(msgid) (msgid)
72 struct hmap_node hmap_node;
80 struct file_handle *outfile;
81 struct string_array variables;
82 struct string_array fnames;
83 struct string_array snames;
88 /* Execution state. */
89 struct dictionary *dict;
90 struct casewriter *writer;
95 struct file_handle *file;
96 struct dfm_reader *reader;
102 struct file_handle *file;
103 struct dfm_writer *writer;
105 struct u8_line *held;
110 struct file_handle *file;
111 struct dataset *dataset;
113 /* Parameters from parsing the first SAVE command for 'file'. */
114 struct string_array variables;
115 struct matrix_expr *names;
116 struct stringi_set strings;
118 /* Results from the first (only) attempt to open this save_file. */
120 struct casewriter *writer;
121 struct dictionary *dict;
126 struct dataset *dataset;
127 struct session *session;
131 struct msave_common *common;
133 struct file_handle *prev_read_file;
134 struct read_file **read_files;
137 struct file_handle *prev_write_file;
138 struct write_file **write_files;
139 size_t n_write_files;
141 struct file_handle *prev_save_file;
142 struct save_file **save_files;
146 static struct matrix_var *
147 matrix_var_lookup (struct matrix_state *s, struct substring name)
149 struct matrix_var *var;
151 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
152 utf8_hash_case_substring (name, 0), &s->vars)
153 if (!utf8_sscasecmp (ss_cstr (var->name), name))
159 static struct matrix_var *
160 matrix_var_create (struct matrix_state *s, struct substring name)
162 struct matrix_var *var = xmalloc (sizeof *var);
163 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
164 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
169 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
171 gsl_matrix_free (var->value);
175 #define MATRIX_FUNCTIONS \
231 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
232 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
233 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
234 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
235 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
236 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
237 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
238 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
239 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
240 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
241 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
243 typedef gsl_matrix *matrix_proto_m_d (double);
244 typedef gsl_matrix *matrix_proto_m_dd (double, double);
245 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
246 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
247 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
248 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
249 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
250 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
251 typedef double matrix_proto_d_m (gsl_matrix *);
252 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
253 typedef gsl_matrix *matrix_proto_IDENT (double, double);
255 #define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
259 /* Matrix expressions. */
266 #define F(NAME, PROTOTYPE) MOP_F_##NAME,
270 /* Elementwise and scalar arithmetic. */
271 MOP_NEGATE, /* unary - */
272 MOP_ADD_ELEMS, /* + */
273 MOP_SUB_ELEMS, /* - */
274 MOP_MUL_ELEMS, /* &* */
275 MOP_DIV_ELEMS, /* / and &/ */
276 MOP_EXP_ELEMS, /* &** */
278 MOP_SEQ_BY, /* a:b:c */
280 /* Matrix arithmetic. */
282 MOP_EXP_MAT, /* ** */
299 MOP_PASTE_HORZ, /* a, b, c, ... */
300 MOP_PASTE_VERT, /* a; b; c; ... */
304 MOP_VEC_INDEX, /* x(y) */
305 MOP_VEC_ALL, /* x(:) */
306 MOP_MAT_INDEX, /* x(y,z) */
307 MOP_ROW_INDEX, /* x(y,:) */
308 MOP_COL_INDEX, /* x(:,z) */
315 MOP_EOF, /* EOF('file') */
323 struct matrix_expr **subs;
328 struct matrix_var *variable;
329 struct read_file *eof;
334 matrix_expr_destroy (struct matrix_expr *e)
341 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
372 for (size_t i = 0; i < e->n_subs; i++)
373 matrix_expr_destroy (e->subs[i]);
385 static struct matrix_expr *
386 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
389 struct matrix_expr *e = xmalloc (sizeof *e);
390 *e = (struct matrix_expr) {
392 .subs = xmemdup (subs, n_subs * sizeof *subs),
398 static struct matrix_expr *
399 matrix_expr_create_0 (enum matrix_op op)
401 struct matrix_expr *sub;
402 return matrix_expr_create_subs (op, &sub, 0);
405 static struct matrix_expr *
406 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
408 return matrix_expr_create_subs (op, &sub, 1);
411 static struct matrix_expr *
412 matrix_expr_create_2 (enum matrix_op op,
413 struct matrix_expr *sub0, struct matrix_expr *sub1)
415 struct matrix_expr *subs[] = { sub0, sub1 };
416 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
419 static struct matrix_expr *
420 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
421 struct matrix_expr *sub1, struct matrix_expr *sub2)
423 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
424 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
427 static struct matrix_expr *
428 matrix_expr_create_number (double number)
430 struct matrix_expr *e = xmalloc (sizeof *e);
431 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
435 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
437 static struct matrix_expr *
438 matrix_parse_curly_comma (struct matrix_state *s)
440 struct matrix_expr *lhs = matrix_parse_or_xor (s);
444 while (lex_match (s->lexer, T_COMMA))
446 struct matrix_expr *rhs = matrix_parse_or_xor (s);
449 matrix_expr_destroy (lhs);
452 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
457 static struct matrix_expr *
458 matrix_parse_curly_semi (struct matrix_state *s)
460 if (lex_token (s->lexer) == T_RCURLY)
461 return matrix_expr_create_0 (MOP_EMPTY);
463 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
467 while (lex_match (s->lexer, T_SEMICOLON))
469 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
472 matrix_expr_destroy (lhs);
475 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
480 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
481 for (size_t Y = 0; Y < (M)->size1; Y++) \
482 for (size_t X = 0; X < (M)->size2; X++) \
483 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
486 is_vector (const gsl_matrix *m)
488 return m->size1 <= 1 || m->size2 <= 1;
492 to_vector (gsl_matrix *m)
494 return (m->size1 == 1
495 ? gsl_matrix_row (m, 0).vector
496 : gsl_matrix_column (m, 0).vector);
501 matrix_eval_ABS (gsl_matrix *m)
503 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
509 matrix_eval_ALL (gsl_matrix *m)
511 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
518 matrix_eval_ANY (gsl_matrix *m)
520 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
527 matrix_eval_ARSIN (gsl_matrix *m)
529 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
535 matrix_eval_ARTAN (gsl_matrix *m)
537 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
543 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
547 for (size_t i = 0; i < n; i++)
552 gsl_matrix *block = gsl_matrix_calloc (r, c);
554 for (size_t i = 0; i < n; i++)
556 for (size_t y = 0; y < m[i]->size1; y++)
557 for (size_t x = 0; x < m[i]->size2; x++)
558 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
566 matrix_eval_CHOL (gsl_matrix *m)
568 if (!gsl_linalg_cholesky_decomp1 (m))
570 for (size_t y = 0; y < m->size1; y++)
571 for (size_t x = y + 1; x < m->size2; x++)
572 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
574 for (size_t y = 0; y < m->size1; y++)
575 for (size_t x = 0; x < y; x++)
576 gsl_matrix_set (m, y, x, 0);
581 msg (SE, _("Input to CHOL function is not positive-definite."));
587 matrix_eval_col_extremum (gsl_matrix *m, bool min)
592 return gsl_matrix_alloc (1, 0);
594 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
595 for (size_t x = 0; x < m->size2; x++)
597 double ext = gsl_matrix_get (m, 0, x);
598 for (size_t y = 1; y < m->size1; y++)
600 double value = gsl_matrix_get (m, y, x);
601 if (min ? value < ext : value > ext)
604 gsl_matrix_set (cext, 0, x, ext);
610 matrix_eval_CMAX (gsl_matrix *m)
612 return matrix_eval_col_extremum (m, false);
616 matrix_eval_CMIN (gsl_matrix *m)
618 return matrix_eval_col_extremum (m, true);
622 matrix_eval_COS (gsl_matrix *m)
624 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
630 matrix_eval_col_sum (gsl_matrix *m, bool square)
635 return gsl_matrix_alloc (1, 0);
637 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
638 for (size_t x = 0; x < m->size2; x++)
641 for (size_t y = 0; y < m->size1; y++)
643 double d = gsl_matrix_get (m, y, x);
644 sum += square ? pow2 (d) : d;
646 gsl_matrix_set (result, 0, x, sum);
652 matrix_eval_CSSQ (gsl_matrix *m)
654 return matrix_eval_col_sum (m, true);
658 matrix_eval_CSUM (gsl_matrix *m)
660 return matrix_eval_col_sum (m, false);
664 compare_double_3way (const void *a_, const void *b_)
666 const double *a = a_;
667 const double *b = b_;
668 return *a < *b ? -1 : *a > *b;
672 matrix_eval_DESIGN (gsl_matrix *m)
674 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
675 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
676 gsl_matrix_transpose_memcpy (&m2, m);
678 for (size_t y = 0; y < m2.size1; y++)
679 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
681 size_t *n = xcalloc (m2.size1, sizeof *n);
683 for (size_t i = 0; i < m2.size1; i++)
685 double *row = tmp + m2.size2 * i;
686 for (size_t j = 0; j < m2.size2; )
689 for (k = j + 1; k < m2.size2; k++)
690 if (row[j] != row[k])
692 row[n[i]++] = row[j];
697 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
702 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
704 for (size_t i = 0; i < m->size2; i++)
709 const double *unique = tmp + m2.size2 * i;
710 for (size_t j = 0; j < n[i]; j++, x++)
712 double value = unique[j];
713 for (size_t y = 0; y < m->size1; y++)
714 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
725 matrix_eval_DET (gsl_matrix *m)
727 gsl_permutation *p = gsl_permutation_alloc (m->size1);
729 gsl_linalg_LU_decomp (m, p, &signum);
730 gsl_permutation_free (p);
731 return gsl_linalg_LU_det (m, signum);
735 matrix_eval_DIAG (gsl_matrix *m)
737 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
738 for (size_t i = 0; i < diag->size1; i++)
739 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
744 is_symmetric (const gsl_matrix *m)
746 if (m->size1 != m->size2)
749 for (size_t y = 0; y < m->size1; y++)
750 for (size_t x = 0; x < y; x++)
751 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
758 compare_double_desc (const void *a_, const void *b_)
760 const double *a = a_;
761 const double *b = b_;
762 return *a > *b ? -1 : *a < *b;
766 matrix_eval_EVAL (gsl_matrix *m)
768 if (!is_symmetric (m))
770 msg (SE, _("Argument of EVAL must be symmetric."));
774 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
775 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
776 gsl_vector v_eval = to_vector (eval);
777 gsl_eigen_symm (m, &v_eval, w);
778 gsl_eigen_symm_free (w);
780 assert (v_eval.stride == 1);
781 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
787 matrix_eval_EXP (gsl_matrix *m)
789 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
794 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
797 Charl Linssen <charl@itfromb.it>
801 matrix_eval_GINV (gsl_matrix *A)
806 gsl_matrix *tmp_mat = NULL;
809 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
810 tmp_mat = gsl_matrix_alloc (m, n);
811 gsl_matrix_transpose_memcpy (tmp_mat, A);
819 gsl_matrix *V = gsl_matrix_alloc (m, m);
820 gsl_vector *u = gsl_vector_alloc (m);
822 gsl_vector *tmp_vec = gsl_vector_alloc (m);
823 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
824 gsl_vector_free (tmp_vec);
827 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
828 gsl_matrix_set_zero (Sigma_pinv);
829 double cutoff = 1e-15 * gsl_vector_max (u);
831 for (size_t i = 0; i < m; ++i)
833 double x = gsl_vector_get (u, i);
834 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
837 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
838 gsl_matrix *U = gsl_matrix_calloc (n, n);
839 for (size_t i = 0; i < n; i++)
840 for (size_t j = 0; j < m; j++)
841 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
843 /* two dot products to obtain pseudoinverse */
844 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
845 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
850 A_pinv = gsl_matrix_alloc (n, m);
851 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
855 A_pinv = gsl_matrix_alloc (m, n);
856 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
859 gsl_matrix_free (tmp_mat);
860 gsl_matrix_free (tmp_mat2);
862 gsl_matrix_free (Sigma_pinv);
876 grade_compare_3way (const void *a_, const void *b_)
878 const struct grade *a = a_;
879 const struct grade *b = b_;
881 return (a->value < b->value ? -1
882 : a->value > b->value ? 1
890 matrix_eval_GRADE (gsl_matrix *m)
892 size_t n = m->size1 * m->size2;
893 struct grade *grades = xmalloc (n * sizeof *grades);
896 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
897 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
898 qsort (grades, n, sizeof *grades, grade_compare_3way);
900 for (size_t i = 0; i < n; i++)
901 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
909 dot (gsl_vector *a, gsl_vector *b)
912 for (size_t i = 0; i < a->size; i++)
913 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
918 norm2 (gsl_vector *v)
921 for (size_t i = 0; i < v->size; i++)
922 result += pow2 (gsl_vector_get (v, i));
929 return sqrt (norm2 (v));
933 matrix_eval_GSCH (gsl_matrix *v)
935 if (v->size2 < v->size1)
937 msg (SE, _("GSCH requires its argument to have at least as many columns "
938 "as rows, but it has dimensions %zu×%zu."),
942 if (!v->size1 || !v->size2)
945 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
947 for (size_t vx = 0; vx < v->size2; vx++)
949 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
950 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
952 gsl_vector_memcpy (&u_i, &v_i);
953 for (size_t j = 0; j < ux; j++)
955 gsl_vector u_j = gsl_matrix_column (u, j).vector;
956 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
957 for (size_t k = 0; k < u_i.size; k++)
958 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
959 - scale * gsl_vector_get (&u_j, k)));
962 double len = norm (&u_i);
965 gsl_vector_scale (&u_i, 1.0 / len);
966 if (++ux >= v->size1)
973 msg (SE, _("%zu×%zu argument to GSCH contains only "
974 "%zu linearly independent columns."),
975 v->size1, v->size2, ux);
985 matrix_eval_IDENT (double s1, double s2)
987 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
989 msg (SE, _("Arguments to IDENT must be non-negative integers."));
993 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
994 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1000 invert_matrix (gsl_matrix *x)
1002 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1004 gsl_linalg_LU_decomp (x, p, &signum);
1005 gsl_linalg_LU_invx (x, p);
1006 gsl_permutation_free (p);
1010 matrix_eval_INV (gsl_matrix *m)
1017 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1019 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1020 a->size2 * b->size2);
1022 for (size_t ar = 0; ar < a->size1; ar++)
1023 for (size_t br = 0; br < b->size1; br++, y++)
1026 for (size_t ac = 0; ac < a->size2; ac++)
1027 for (size_t bc = 0; bc < b->size2; bc++, x++)
1029 double av = gsl_matrix_get (a, ar, ac);
1030 double bv = gsl_matrix_get (b, br, bc);
1031 gsl_matrix_set (k, y, x, av * bv);
1038 matrix_eval_LG10 (gsl_matrix *m)
1040 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1046 matrix_eval_LN (gsl_matrix *m)
1048 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1054 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1056 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1059 for (size_t i = 1; i <= n * n; i++)
1061 gsl_matrix_set (m, y, x, i);
1063 size_t y1 = !y ? n - 1 : y - 1;
1064 size_t x1 = x + 1 >= n ? 0 : x + 1;
1065 if (gsl_matrix_get (m, y1, x1) == 0)
1071 y = y + 1 >= n ? 0 : y + 1;
1076 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1078 double a = gsl_matrix_get (m, y1, x1);
1079 double b = gsl_matrix_get (m, y2, x2);
1080 gsl_matrix_set (m, y1, x1, b);
1081 gsl_matrix_set (m, y2, x2, a);
1085 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1089 /* A. Umar, "On the Construction of Even Order Magic Squares",
1090 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1092 for (size_t i = 1; i <= n * n / 2; i++)
1094 gsl_matrix_set (m, y, x, i);
1104 for (size_t i = n * n; i > n * n / 2; i--)
1106 gsl_matrix_set (m, y, x, i);
1114 for (size_t y = 0; y < n; y++)
1115 for (size_t x = 0; x < n / 2; x++)
1117 unsigned int d = gsl_matrix_get (m, y, x);
1118 if (d % 2 != (y < n / 2))
1119 magic_exchange (m, y, x, y, n - x - 1);
1124 size_t x1 = n / 2 - 1;
1126 magic_exchange (m, y1, x1, y2, x1);
1127 magic_exchange (m, y1, x2, y2, x2);
1131 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1133 /* A. Umar, "On the Construction of Even Order Magic Squares",
1134 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1138 for (size_t i = 1; ; i++)
1140 gsl_matrix_set (m, y, x, i);
1141 if (++y == n / 2 - 1)
1153 for (size_t i = n * n; ; i--)
1155 gsl_matrix_set (m, y, x, i);
1156 if (++y == n / 2 - 1)
1165 for (size_t y = 0; y < n; y++)
1166 if (y != n / 2 - 1 && y != n / 2)
1167 for (size_t x = 0; x < n / 2; x++)
1169 unsigned int d = gsl_matrix_get (m, y, x);
1170 if (d % 2 != (y < n / 2))
1171 magic_exchange (m, y, x, y, n - x - 1);
1174 size_t a0 = (n * n - 2 * n) / 2 + 1;
1175 for (size_t i = 0; i < n / 2; i++)
1178 gsl_matrix_set (m, n / 2, i, a);
1179 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1181 for (size_t i = 0; i < n / 2; i++)
1183 size_t a = a0 + i + n / 2;
1184 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1185 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1187 for (size_t x = 1; x < n / 2; x += 2)
1188 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1189 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1190 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1191 size_t x1 = n / 2 - 2;
1192 size_t x2 = n / 2 + 1;
1193 size_t y1 = n / 2 - 2;
1194 size_t y2 = n / 2 + 1;
1195 magic_exchange (m, y1, x1, y2, x1);
1196 magic_exchange (m, y1, x2, y2, x2);
1200 matrix_eval_MAGIC (double n_)
1202 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1204 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1209 gsl_matrix *m = gsl_matrix_calloc (n, n);
1211 matrix_eval_MAGIC_odd (m, n);
1213 matrix_eval_MAGIC_singly_even (m, n);
1215 matrix_eval_MAGIC_doubly_even (m, n);
1220 matrix_eval_MAKE (double r, double c, double value)
1222 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1224 msg (SE, _("Size arguments to MAKE must be integers."));
1228 gsl_matrix *m = gsl_matrix_alloc (r, c);
1229 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1235 matrix_eval_MDIAG (gsl_vector *v)
1237 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1238 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1239 gsl_vector_memcpy (&diagonal, v);
1244 matrix_eval_MMAX (gsl_matrix *m)
1246 return gsl_matrix_max (m);
1250 matrix_eval_MMIN (gsl_matrix *m)
1252 return gsl_matrix_min (m);
1256 matrix_eval_MOD (gsl_matrix *m, double divisor)
1260 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1264 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1265 *d = fmod (*d, divisor);
1270 matrix_eval_MSSQ (gsl_matrix *m)
1273 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1279 matrix_eval_MSUM (gsl_matrix *m)
1282 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1288 matrix_eval_NCOL (gsl_matrix *m)
1294 matrix_eval_NROW (gsl_matrix *m)
1300 matrix_eval_RANK (gsl_matrix *m)
1302 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1303 gsl_linalg_QR_decomp (m, tau);
1304 gsl_vector_free (tau);
1306 return gsl_linalg_QRPT_rank (m, -1);
1310 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1312 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1314 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1319 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1321 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1322 "product of matrix dimensions."));
1326 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1329 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1331 gsl_matrix_set (dst, y1, x1, *d);
1342 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1347 return gsl_matrix_alloc (0, 1);
1349 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1350 for (size_t y = 0; y < m->size1; y++)
1352 double ext = gsl_matrix_get (m, y, 0);
1353 for (size_t x = 1; x < m->size2; x++)
1355 double value = gsl_matrix_get (m, y, x);
1356 if (min ? value < ext : value > ext)
1359 gsl_matrix_set (rext, y, 0, ext);
1365 matrix_eval_RMAX (gsl_matrix *m)
1367 return matrix_eval_row_extremum (m, false);
1371 matrix_eval_RMIN (gsl_matrix *m)
1373 return matrix_eval_row_extremum (m, true);
1377 matrix_eval_RND (gsl_matrix *m)
1379 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1391 rank_compare_3way (const void *a_, const void *b_)
1393 const struct rank *a = a_;
1394 const struct rank *b = b_;
1396 return a->value < b->value ? -1 : a->value > b->value;
1400 matrix_eval_RNKORDER (gsl_matrix *m)
1402 size_t n = m->size1 * m->size2;
1403 struct rank *ranks = xmalloc (n * sizeof *ranks);
1405 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1406 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1407 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1409 for (size_t i = 0; i < n; )
1412 for (j = i + 1; j < n; j++)
1413 if (ranks[i].value != ranks[j].value)
1416 double rank = (i + j + 1.0) / 2.0;
1417 for (size_t k = i; k < j; k++)
1418 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1429 matrix_eval_row_sum (gsl_matrix *m, bool square)
1434 return gsl_matrix_alloc (0, 1);
1436 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1437 for (size_t y = 0; y < m->size1; y++)
1440 for (size_t x = 0; x < m->size2; x++)
1442 double d = gsl_matrix_get (m, y, x);
1443 sum += square ? pow2 (d) : d;
1445 gsl_matrix_set (result, y, 0, sum);
1451 matrix_eval_RSSQ (gsl_matrix *m)
1453 return matrix_eval_row_sum (m, true);
1457 matrix_eval_RSUM (gsl_matrix *m)
1459 return matrix_eval_row_sum (m, false);
1463 matrix_eval_SIN (gsl_matrix *m)
1465 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1471 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1473 if (m1->size1 != m2->size1)
1475 msg (SE, _("SOLVE requires its arguments to have the same number of "
1476 "rows, but the first argument has dimensions %zu×%zu and "
1477 "the second %zu×%zu."),
1478 m1->size1, m1->size2,
1479 m2->size1, m2->size2);
1483 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1484 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1486 gsl_linalg_LU_decomp (m1, p, &signum);
1487 for (size_t i = 0; i < m2->size2; i++)
1489 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1490 gsl_vector xi = gsl_matrix_column (x, i).vector;
1491 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1493 gsl_permutation_free (p);
1498 matrix_eval_SQRT (gsl_matrix *m)
1500 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1504 msg (SE, _("Argument to SQRT must be nonnegative."));
1513 matrix_eval_SSCP (gsl_matrix *m)
1515 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1516 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1521 matrix_eval_SVAL (gsl_matrix *m)
1523 gsl_matrix *tmp_mat = NULL;
1524 if (m->size2 > m->size1)
1526 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1527 gsl_matrix_transpose_memcpy (tmp_mat, m);
1532 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1533 gsl_vector *S = gsl_vector_alloc (m->size2);
1534 gsl_vector *work = gsl_vector_alloc (m->size2);
1535 gsl_linalg_SV_decomp (m, V, S, work);
1537 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1538 for (size_t i = 0; i < m->size2; i++)
1539 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1541 gsl_matrix_free (V);
1542 gsl_vector_free (S);
1543 gsl_vector_free (work);
1544 gsl_matrix_free (tmp_mat);
1550 matrix_eval_SWEEP (gsl_matrix *m, double d)
1552 if (d < 1 || d > SIZE_MAX)
1554 msg (SE, _("Scalar argument to SWEEP must be integer."));
1558 if (k >= MIN (m->size1, m->size2))
1560 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1561 "equal to the smaller of the matrix argument's rows and "
1566 double m_kk = gsl_matrix_get (m, k, k);
1567 if (fabs (m_kk) > 1e-19)
1569 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1570 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1572 double m_ij = gsl_matrix_get (m, i, j);
1573 double m_ik = gsl_matrix_get (m, i, k);
1574 double m_kj = gsl_matrix_get (m, k, j);
1575 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1584 for (size_t i = 0; i < m->size1; i++)
1586 gsl_matrix_set (m, i, k, 0);
1587 gsl_matrix_set (m, k, i, 0);
1594 matrix_eval_TRACE (gsl_matrix *m)
1597 size_t n = MIN (m->size1, m->size2);
1598 for (size_t i = 0; i < n; i++)
1599 sum += gsl_matrix_get (m, i, i);
1604 matrix_eval_T (gsl_matrix *m)
1606 return matrix_eval_TRANSPOS (m);
1610 matrix_eval_TRANSPOS (gsl_matrix *m)
1612 if (m->size1 == m->size2)
1614 gsl_matrix_transpose (m);
1619 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1620 gsl_matrix_transpose_memcpy (t, m);
1626 matrix_eval_TRUNC (gsl_matrix *m)
1628 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1634 matrix_eval_UNIFORM (double r_, double c_)
1636 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1638 msg (SE, _("Arguments to UNIFORM must be integers."));
1643 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1645 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1649 gsl_matrix *m = gsl_matrix_alloc (r, c);
1650 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1651 *d = gsl_ran_flat (get_rng (), 0, 1);
1655 struct matrix_function
1659 size_t min_args, max_args;
1662 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1664 static const struct matrix_function *
1665 matrix_parse_function_name (const char *token)
1667 static const struct matrix_function functions[] =
1669 #define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1673 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1675 for (size_t i = 0; i < N_FUNCTIONS; i++)
1677 if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1678 return &functions[i];
1683 static struct read_file *
1684 read_file_create (struct matrix_state *s, struct file_handle *fh)
1686 for (size_t i = 0; i < s->n_read_files; i++)
1688 struct read_file *rf = s->read_files[i];
1696 struct read_file *rf = xmalloc (sizeof *rf);
1697 *rf = (struct read_file) { .file = fh };
1699 s->read_files = xrealloc (s->read_files,
1700 (s->n_read_files + 1) * sizeof *s->read_files);
1701 s->read_files[s->n_read_files++] = rf;
1705 static struct dfm_reader *
1706 read_file_open (struct read_file *rf)
1709 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1714 read_file_destroy (struct read_file *rf)
1718 fh_unref (rf->file);
1719 dfm_close_reader (rf->reader);
1720 free (rf->encoding);
1725 static struct write_file *
1726 write_file_create (struct matrix_state *s, struct file_handle *fh)
1728 for (size_t i = 0; i < s->n_write_files; i++)
1730 struct write_file *wf = s->write_files[i];
1738 struct write_file *wf = xmalloc (sizeof *wf);
1739 *wf = (struct write_file) { .file = fh };
1741 s->write_files = xrealloc (s->write_files,
1742 (s->n_write_files + 1) * sizeof *s->write_files);
1743 s->write_files[s->n_write_files++] = wf;
1747 static struct dfm_writer *
1748 write_file_open (struct write_file *wf)
1751 wf->writer = dfm_open_writer (wf->file, wf->encoding);
1756 write_file_destroy (struct write_file *wf)
1762 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
1763 wf->held->s.ss.length);
1764 u8_line_destroy (wf->held);
1768 fh_unref (wf->file);
1769 dfm_close_writer (wf->writer);
1770 free (wf->encoding);
1776 matrix_parse_function (struct matrix_state *s, const char *token,
1777 struct matrix_expr **exprp)
1780 if (lex_next_token (s->lexer, 1) != T_LPAREN)
1783 if (lex_match_id (s->lexer, "EOF"))
1786 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1790 if (!lex_force_match (s->lexer, T_RPAREN))
1796 struct read_file *rf = read_file_create (s, fh);
1798 struct matrix_expr *e = xmalloc (sizeof *e);
1799 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1804 const struct matrix_function *f = matrix_parse_function_name (token);
1808 lex_get_n (s->lexer, 2);
1810 struct matrix_expr *e = xmalloc (sizeof *e);
1811 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1813 size_t allocated_subs = 0;
1816 struct matrix_expr *sub = matrix_parse_expr (s);
1820 if (e->n_subs >= allocated_subs)
1821 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1822 e->subs[e->n_subs++] = sub;
1824 while (lex_match (s->lexer, T_COMMA));
1825 if (!lex_force_match (s->lexer, T_RPAREN))
1832 matrix_expr_destroy (e);
1836 static struct matrix_expr *
1837 matrix_parse_primary (struct matrix_state *s)
1839 if (lex_is_number (s->lexer))
1841 double number = lex_number (s->lexer);
1843 return matrix_expr_create_number (number);
1845 else if (lex_is_string (s->lexer))
1847 char string[sizeof (double)];
1848 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1852 memcpy (&number, string, sizeof number);
1853 return matrix_expr_create_number (number);
1855 else if (lex_match (s->lexer, T_LPAREN))
1857 struct matrix_expr *e = matrix_parse_or_xor (s);
1858 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1860 matrix_expr_destroy (e);
1865 else if (lex_match (s->lexer, T_LCURLY))
1867 struct matrix_expr *e = matrix_parse_curly_semi (s);
1868 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1870 matrix_expr_destroy (e);
1875 else if (lex_token (s->lexer) == T_ID)
1877 struct matrix_expr *retval;
1878 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1881 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1884 lex_error (s->lexer, _("Unknown variable %s."),
1885 lex_tokcstr (s->lexer));
1890 struct matrix_expr *e = xmalloc (sizeof *e);
1891 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1894 else if (lex_token (s->lexer) == T_ALL)
1896 struct matrix_expr *retval;
1897 if (matrix_parse_function (s, "ALL", &retval))
1901 lex_error (s->lexer, NULL);
1905 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1908 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1910 if (lex_match (s->lexer, T_COLON))
1917 *indexp = matrix_parse_expr (s);
1918 return *indexp != NULL;
1922 static struct matrix_expr *
1923 matrix_parse_postfix (struct matrix_state *s)
1925 struct matrix_expr *lhs = matrix_parse_primary (s);
1926 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1929 struct matrix_expr *i0;
1930 if (!matrix_parse_index_expr (s, &i0))
1932 matrix_expr_destroy (lhs);
1935 if (lex_match (s->lexer, T_RPAREN))
1937 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1938 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1939 else if (lex_match (s->lexer, T_COMMA))
1941 struct matrix_expr *i1;
1942 if (!matrix_parse_index_expr (s, &i1)
1943 || !lex_force_match (s->lexer, T_RPAREN))
1945 matrix_expr_destroy (lhs);
1946 matrix_expr_destroy (i0);
1947 matrix_expr_destroy (i1);
1950 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1951 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1952 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1957 lex_error_expecting (s->lexer, "`)'", "`,'");
1962 static struct matrix_expr *
1963 matrix_parse_unary (struct matrix_state *s)
1965 if (lex_match (s->lexer, T_DASH))
1967 struct matrix_expr *lhs = matrix_parse_unary (s);
1968 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1970 else if (lex_match (s->lexer, T_PLUS))
1971 return matrix_parse_unary (s);
1973 return matrix_parse_postfix (s);
1976 static struct matrix_expr *
1977 matrix_parse_seq (struct matrix_state *s)
1979 struct matrix_expr *start = matrix_parse_unary (s);
1980 if (!start || !lex_match (s->lexer, T_COLON))
1983 struct matrix_expr *end = matrix_parse_unary (s);
1986 matrix_expr_destroy (start);
1990 if (lex_match (s->lexer, T_COLON))
1992 struct matrix_expr *increment = matrix_parse_unary (s);
1995 matrix_expr_destroy (start);
1996 matrix_expr_destroy (end);
1999 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2002 return matrix_expr_create_2 (MOP_SEQ, start, end);
2005 static struct matrix_expr *
2006 matrix_parse_exp (struct matrix_state *s)
2008 struct matrix_expr *lhs = matrix_parse_seq (s);
2015 if (lex_match (s->lexer, T_EXP))
2017 else if (lex_match_phrase (s->lexer, "&**"))
2022 struct matrix_expr *rhs = matrix_parse_seq (s);
2025 matrix_expr_destroy (lhs);
2028 lhs = matrix_expr_create_2 (op, lhs, rhs);
2032 static struct matrix_expr *
2033 matrix_parse_mul_div (struct matrix_state *s)
2035 struct matrix_expr *lhs = matrix_parse_exp (s);
2042 if (lex_match (s->lexer, T_ASTERISK))
2044 else if (lex_match (s->lexer, T_SLASH))
2046 else if (lex_match_phrase (s->lexer, "&*"))
2048 else if (lex_match_phrase (s->lexer, "&/"))
2053 struct matrix_expr *rhs = matrix_parse_exp (s);
2056 matrix_expr_destroy (lhs);
2059 lhs = matrix_expr_create_2 (op, lhs, rhs);
2063 static struct matrix_expr *
2064 matrix_parse_add_sub (struct matrix_state *s)
2066 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2073 if (lex_match (s->lexer, T_PLUS))
2075 else if (lex_match (s->lexer, T_DASH))
2077 else if (lex_token (s->lexer) == T_NEG_NUM)
2082 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2085 matrix_expr_destroy (lhs);
2088 lhs = matrix_expr_create_2 (op, lhs, rhs);
2092 static struct matrix_expr *
2093 matrix_parse_relational (struct matrix_state *s)
2095 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2102 if (lex_match (s->lexer, T_GT))
2104 else if (lex_match (s->lexer, T_GE))
2106 else if (lex_match (s->lexer, T_LT))
2108 else if (lex_match (s->lexer, T_LE))
2110 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2112 else if (lex_match (s->lexer, T_NE))
2117 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2120 matrix_expr_destroy (lhs);
2123 lhs = matrix_expr_create_2 (op, lhs, rhs);
2127 static struct matrix_expr *
2128 matrix_parse_not (struct matrix_state *s)
2130 if (lex_match (s->lexer, T_NOT))
2132 struct matrix_expr *lhs = matrix_parse_not (s);
2133 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2136 return matrix_parse_relational (s);
2139 static struct matrix_expr *
2140 matrix_parse_and (struct matrix_state *s)
2142 struct matrix_expr *lhs = matrix_parse_not (s);
2146 while (lex_match (s->lexer, T_AND))
2148 struct matrix_expr *rhs = matrix_parse_not (s);
2151 matrix_expr_destroy (lhs);
2154 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2159 static struct matrix_expr *
2160 matrix_parse_or_xor (struct matrix_state *s)
2162 struct matrix_expr *lhs = matrix_parse_and (s);
2169 if (lex_match (s->lexer, T_OR))
2171 else if (lex_match_id (s->lexer, "XOR"))
2176 struct matrix_expr *rhs = matrix_parse_and (s);
2179 matrix_expr_destroy (lhs);
2182 lhs = matrix_expr_create_2 (op, lhs, rhs);
2186 static struct matrix_expr *
2187 matrix_parse_expr (struct matrix_state *s)
2189 return matrix_parse_or_xor (s);
2192 /* Expression evaluation. */
2195 matrix_op_eval (enum matrix_op op, double a, double b)
2199 case MOP_ADD_ELEMS: return a + b;
2200 case MOP_SUB_ELEMS: return a - b;
2201 case MOP_MUL_ELEMS: return a * b;
2202 case MOP_DIV_ELEMS: return a / b;
2203 case MOP_EXP_ELEMS: return pow (a, b);
2204 case MOP_GT: return a > b;
2205 case MOP_GE: return a >= b;
2206 case MOP_LT: return a < b;
2207 case MOP_LE: return a <= b;
2208 case MOP_EQ: return a == b;
2209 case MOP_NE: return a != b;
2210 case MOP_AND: return (a > 0) && (b > 0);
2211 case MOP_OR: return (a > 0) || (b > 0);
2212 case MOP_XOR: return (a > 0) != (b > 0);
2214 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2223 case MOP_PASTE_HORZ:
2224 case MOP_PASTE_VERT:
2240 matrix_op_name (enum matrix_op op)
2244 case MOP_ADD_ELEMS: return "+";
2245 case MOP_SUB_ELEMS: return "-";
2246 case MOP_MUL_ELEMS: return "&*";
2247 case MOP_DIV_ELEMS: return "&/";
2248 case MOP_EXP_ELEMS: return "&**";
2249 case MOP_GT: return ">";
2250 case MOP_GE: return ">=";
2251 case MOP_LT: return "<";
2252 case MOP_LE: return "<=";
2253 case MOP_EQ: return "=";
2254 case MOP_NE: return "<>";
2255 case MOP_AND: return "AND";
2256 case MOP_OR: return "OR";
2257 case MOP_XOR: return "XOR";
2259 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2268 case MOP_PASTE_HORZ:
2269 case MOP_PASTE_VERT:
2285 is_scalar (const gsl_matrix *m)
2287 return m->size1 == 1 && m->size2 == 1;
2291 to_scalar (const gsl_matrix *m)
2293 assert (is_scalar (m));
2294 return gsl_matrix_get (m, 0, 0);
2298 matrix_expr_evaluate_elementwise (enum matrix_op op,
2299 gsl_matrix *a, gsl_matrix *b)
2303 double be = to_scalar (b);
2304 for (size_t r = 0; r < a->size1; r++)
2305 for (size_t c = 0; c < a->size2; c++)
2307 double *ae = gsl_matrix_ptr (a, r, c);
2308 *ae = matrix_op_eval (op, *ae, be);
2312 else if (is_scalar (a))
2314 double ae = to_scalar (a);
2315 for (size_t r = 0; r < b->size1; r++)
2316 for (size_t c = 0; c < b->size2; c++)
2318 double *be = gsl_matrix_ptr (b, r, c);
2319 *be = matrix_op_eval (op, ae, *be);
2323 else if (a->size1 == b->size1 && a->size2 == b->size2)
2325 for (size_t r = 0; r < a->size1; r++)
2326 for (size_t c = 0; c < a->size2; c++)
2328 double *ae = gsl_matrix_ptr (a, r, c);
2329 double be = gsl_matrix_get (b, r, c);
2330 *ae = matrix_op_eval (op, *ae, be);
2336 msg (SE, _("Operands to %s must have the same dimensions or one "
2337 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
2338 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2344 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2346 if (is_scalar (a) || is_scalar (b))
2347 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2349 if (a->size2 != b->size1)
2351 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
2352 "not conformable for multiplication."),
2353 a->size1, a->size2, b->size1, b->size2);
2357 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2358 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2363 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2365 gsl_matrix *tmp = *a;
2371 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2374 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2375 swap_matrix (z, tmp);
2379 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2381 mul_matrix (x, *x, *x, tmp);
2385 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2388 if (x->size1 != x->size2)
2390 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2391 "the left-hand size, not one with dimensions %zu×%zu."),
2392 x->size1, x->size2);
2397 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2398 "right-hand side, not a matrix with dimensions %zu×%zu."),
2399 b->size1, b->size2);
2402 double bf = to_scalar (b);
2403 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2405 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2406 "or outside the valid range."), bf);
2411 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2413 gsl_matrix_set_identity (y);
2417 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2419 for (unsigned long int n = labs (bl); n > 1; n /= 2)
2422 mul_matrix (&y, x, y, &t);
2423 square_matrix (&x, &t);
2426 square_matrix (&x, &t);
2428 mul_matrix (&y, x, y, &t);
2432 /* Garbage collection.
2434 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2435 a permutation of them. We are returning one of them; that one must not be
2436 destroyed. We must not destroy 'x_' because the caller owns it. */
2438 gsl_matrix_free (y_);
2440 gsl_matrix_free (t_);
2446 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2449 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2451 msg (SE, _("All operands of : operator must be scalars."));
2455 long int start = to_scalar (start_);
2456 long int end = to_scalar (end_);
2457 long int by = by_ ? to_scalar (by_) : 1;
2461 msg (SE, _("The increment operand to : must be nonzero."));
2465 long int n = (end >= start && by > 0 ? (end - start + by) / by
2466 : end <= start && by < 0 ? (start - end - by) / -by
2468 gsl_matrix *m = gsl_matrix_alloc (1, n);
2469 for (long int i = 0; i < n; i++)
2470 gsl_matrix_set (m, 0, i, start + i * by);
2475 matrix_expr_evaluate_not (gsl_matrix *a)
2477 for (size_t r = 0; r < a->size1; r++)
2478 for (size_t c = 0; c < a->size2; c++)
2480 double *ae = gsl_matrix_ptr (a, r, c);
2487 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2489 if (a->size1 != b->size1)
2491 if (!a->size1 || !a->size2)
2493 else if (!b->size1 || !b->size2)
2496 msg (SE, _("All columns in a matrix must have the same number of rows, "
2497 "but this tries to paste matrices with %zu and %zu rows."),
2498 a->size1, b->size1);
2502 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2503 for (size_t y = 0; y < a->size1; y++)
2505 for (size_t x = 0; x < a->size2; x++)
2506 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2507 for (size_t x = 0; x < b->size2; x++)
2508 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2514 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2516 if (a->size2 != b->size2)
2518 if (!a->size1 || !a->size2)
2520 else if (!b->size1 || !b->size2)
2523 msg (SE, _("All rows in a matrix must have the same number of columns, "
2524 "but this tries to stack matrices with %zu and %zu columns."),
2525 a->size2, b->size2);
2529 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2530 for (size_t x = 0; x < a->size2; x++)
2532 for (size_t y = 0; y < a->size1; y++)
2533 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2534 for (size_t y = 0; y < b->size1; y++)
2535 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2541 matrix_to_vector (gsl_matrix *m)
2544 gsl_vector v = to_vector (m);
2545 assert (v.block == m->block || !v.block);
2549 return xmemdup (&v, sizeof v);
2563 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2566 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2567 enum index_type index_type, size_t other_size,
2568 struct index_vector *iv)
2577 msg (SE, _("Vector index must be scalar or vector, not a "
2579 m->size1, m->size2);
2583 msg (SE, _("Matrix row index must be scalar or vector, not a "
2585 m->size1, m->size2);
2589 msg (SE, _("Matrix column index must be scalar or vector, not a "
2591 m->size1, m->size2);
2597 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2598 *iv = (struct index_vector) {
2599 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2602 for (size_t i = 0; i < v.size; i++)
2604 double index = gsl_vector_get (&v, i);
2605 if (index < 1 || index >= size + 1)
2610 msg (SE, _("Index %g is out of range for vector "
2611 "with %zu elements."), index, size);
2615 msg (SE, _("%g is not a valid row index for "
2616 "a %zu×%zu matrix."),
2617 index, size, other_size);
2621 msg (SE, _("%g is not a valid column index for "
2622 "a %zu×%zu matrix."),
2623 index, other_size, size);
2630 iv->indexes[i] = index - 1;
2636 *iv = (struct index_vector) {
2637 .indexes = xnmalloc (size, sizeof *iv->indexes),
2640 for (size_t i = 0; i < size; i++)
2647 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2649 if (!is_vector (sm))
2651 msg (SE, _("Vector index operator may not be applied to "
2652 "a %zu×%zu matrix."),
2653 sm->size1, sm->size2);
2661 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2663 if (!matrix_expr_evaluate_vec_all (sm))
2666 gsl_vector sv = to_vector (sm);
2667 struct index_vector iv;
2668 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2671 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2672 sm->size1 == 1 ? iv.n : 1);
2673 gsl_vector dv = to_vector (dm);
2674 for (size_t dx = 0; dx < iv.n; dx++)
2676 size_t sx = iv.indexes[dx];
2677 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2685 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2688 struct index_vector iv0;
2689 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2692 struct index_vector iv1;
2693 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2700 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2701 for (size_t dy = 0; dy < iv0.n; dy++)
2703 size_t sy = iv0.indexes[dy];
2705 for (size_t dx = 0; dx < iv1.n; dx++)
2707 size_t sx = iv1.indexes[dx];
2708 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2716 #define F(NAME, PROTOTYPE) \
2717 static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2718 const char *name, gsl_matrix *[], size_t, \
2719 matrix_proto_##PROTOTYPE *);
2724 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2726 if (!is_scalar (subs[index]))
2728 msg (SE, _("Function %s argument %zu must be a scalar, "
2729 "not a %zu×%zu matrix."),
2730 name, index + 1, subs[index]->size1, subs[index]->size2);
2737 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2739 if (!is_vector (subs[index]))
2741 msg (SE, _("Function %s argument %zu must be a vector, "
2742 "not a %zu×%zu matrix."),
2743 name, index + 1, subs[index]->size1, subs[index]->size2);
2750 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2752 for (size_t i = 0; i < n_subs; i++)
2754 if (!check_scalar_arg (name, subs, i))
2756 d[i] = to_scalar (subs[i]);
2762 matrix_expr_evaluate_m_d (const char *name,
2763 gsl_matrix *subs[], size_t n_subs,
2764 matrix_proto_m_d *f)
2766 assert (n_subs == 1);
2769 return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2773 matrix_expr_evaluate_m_dd (const char *name,
2774 gsl_matrix *subs[], size_t n_subs,
2775 matrix_proto_m_dd *f)
2777 assert (n_subs == 2);
2780 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2784 matrix_expr_evaluate_m_ddd (const char *name,
2785 gsl_matrix *subs[], size_t n_subs,
2786 matrix_proto_m_ddd *f)
2788 assert (n_subs == 3);
2791 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2795 matrix_expr_evaluate_m_m (const char *name UNUSED,
2796 gsl_matrix *subs[], size_t n_subs,
2797 matrix_proto_m_m *f)
2799 assert (n_subs == 1);
2804 matrix_expr_evaluate_m_md (const char *name UNUSED,
2805 gsl_matrix *subs[], size_t n_subs,
2806 matrix_proto_m_md *f)
2808 assert (n_subs == 2);
2809 if (!check_scalar_arg (name, subs, 1))
2811 return f (subs[0], to_scalar (subs[1]));
2815 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2816 gsl_matrix *subs[], size_t n_subs,
2817 matrix_proto_m_mdd *f)
2819 assert (n_subs == 3);
2820 if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2822 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2826 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2827 gsl_matrix *subs[], size_t n_subs,
2828 matrix_proto_m_mm *f)
2830 assert (n_subs == 2);
2831 return f (subs[0], subs[1]);
2835 matrix_expr_evaluate_m_v (const char *name,
2836 gsl_matrix *subs[], size_t n_subs,
2837 matrix_proto_m_v *f)
2839 assert (n_subs == 1);
2840 if (!check_vector_arg (name, subs, 0))
2842 gsl_vector v = to_vector (subs[0]);
2847 matrix_expr_evaluate_d_m (const char *name UNUSED,
2848 gsl_matrix *subs[], size_t n_subs,
2849 matrix_proto_d_m *f)
2851 assert (n_subs == 1);
2853 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2854 gsl_matrix_set (m, 0, 0, f (subs[0]));
2859 matrix_expr_evaluate_m_any (const char *name UNUSED,
2860 gsl_matrix *subs[], size_t n_subs,
2861 matrix_proto_m_any *f)
2863 return f (subs, n_subs);
2867 matrix_expr_evaluate_IDENT (const char *name,
2868 gsl_matrix *subs[], size_t n_subs,
2869 matrix_proto_IDENT *f)
2871 assert (n_subs <= 2);
2874 if (!to_scalar_args (name, subs, n_subs, d))
2877 return f (d[0], d[n_subs - 1]);
2881 matrix_expr_evaluate (const struct matrix_expr *e)
2883 if (e->op == MOP_NUMBER)
2885 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2886 gsl_matrix_set (m, 0, 0, e->number);
2889 else if (e->op == MOP_VARIABLE)
2891 const gsl_matrix *src = e->variable->value;
2894 msg (SE, _("Uninitialized variable %s used in expression."),
2899 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2900 gsl_matrix_memcpy (dst, src);
2903 else if (e->op == MOP_EOF)
2905 struct dfm_reader *reader = read_file_open (e->eof);
2906 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2907 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2911 enum { N_LOCAL = 3 };
2912 gsl_matrix *local_subs[N_LOCAL];
2913 gsl_matrix **subs = (e->n_subs < N_LOCAL
2915 : xmalloc (e->n_subs * sizeof *subs));
2917 for (size_t i = 0; i < e->n_subs; i++)
2919 subs[i] = matrix_expr_evaluate (e->subs[i]);
2922 for (size_t j = 0; j < i; j++)
2923 gsl_matrix_free (subs[j]);
2924 if (subs != local_subs)
2930 gsl_matrix *result = NULL;
2933 #define F(NAME, PROTOTYPE) \
2934 case MOP_F_##NAME: \
2935 result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \
2937 matrix_eval_##NAME); \
2943 gsl_matrix_scale (subs[0], -1.0);
2961 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2965 result = matrix_expr_evaluate_not (subs[0]);
2969 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2973 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2977 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2981 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2984 case MOP_PASTE_HORZ:
2985 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2988 case MOP_PASTE_VERT:
2989 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2993 result = gsl_matrix_alloc (0, 0);
2997 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
3001 result = matrix_expr_evaluate_vec_all (subs[0]);
3005 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
3009 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
3013 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
3022 for (size_t i = 0; i < e->n_subs; i++)
3023 if (subs[i] != result)
3024 gsl_matrix_free (subs[i]);
3025 if (subs != local_subs)
3031 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
3034 gsl_matrix *m = matrix_expr_evaluate (e);
3040 msg (SE, _("Expression for %s must evaluate to scalar, "
3041 "not a %zu×%zu matrix."),
3042 context, m->size1, m->size2);
3043 gsl_matrix_free (m);
3048 gsl_matrix_free (m);
3053 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
3057 if (!matrix_expr_evaluate_scalar (e, context, &d))
3061 if (d < LONG_MIN || d > LONG_MAX)
3063 msg (SE, _("Expression for %s is outside the integer range."), context);
3070 struct matrix_lvalue
3072 struct matrix_var *var;
3073 struct matrix_expr *indexes[2];
3078 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
3082 for (size_t i = 0; i < lvalue->n_indexes; i++)
3083 matrix_expr_destroy (lvalue->indexes[i]);
3088 static struct matrix_lvalue *
3089 matrix_lvalue_parse (struct matrix_state *s)
3091 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3092 if (!lex_force_id (s->lexer))
3094 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3095 if (lex_next_token (s->lexer, 1) == T_LPAREN)
3099 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3104 lex_get_n (s->lexer, 2);
3106 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3108 lvalue->n_indexes++;
3110 if (lex_match (s->lexer, T_COMMA))
3112 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3114 lvalue->n_indexes++;
3116 if (!lex_force_match (s->lexer, T_RPAREN))
3122 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3128 matrix_lvalue_destroy (lvalue);
3133 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3134 enum index_type index_type, size_t other_size,
3135 struct index_vector *iv)
3140 m = matrix_expr_evaluate (e);
3147 bool ok = matrix_normalize_index_vector (m, size, index_type,
3149 gsl_matrix_free (m);
3154 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3155 struct index_vector *iv, gsl_matrix *sm)
3157 gsl_vector dv = to_vector (lvalue->var->value);
3159 /* Convert source matrix 'sm' to source vector 'sv'. */
3160 if (!is_vector (sm))
3162 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
3163 sm->size1, sm->size2);
3166 gsl_vector sv = to_vector (sm);
3167 if (iv->n != sv.size)
3169 msg (SE, _("Can't assign %zu-element vector "
3170 "to %zu-element subvector."), sv.size, iv->n);
3174 /* Assign elements. */
3175 for (size_t x = 0; x < iv->n; x++)
3176 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3181 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3182 struct index_vector *iv0,
3183 struct index_vector *iv1, gsl_matrix *sm)
3185 gsl_matrix *dm = lvalue->var->value;
3187 /* Convert source matrix 'sm' to source vector 'sv'. */
3188 if (iv0->n != sm->size1)
3190 msg (SE, _("Row index vector for assignment to %s has %zu elements "
3191 "but source matrix has %zu rows."),
3192 lvalue->var->name, iv0->n, sm->size1);
3195 if (iv1->n != sm->size2)
3197 msg (SE, _("Column index vector for assignment to %s has %zu elements "
3198 "but source matrix has %zu columns."),
3199 lvalue->var->name, iv1->n, sm->size2);
3203 /* Assign elements. */
3204 for (size_t y = 0; y < iv0->n; y++)
3206 size_t f0 = iv0->indexes[y];
3207 for (size_t x = 0; x < iv1->n; x++)
3209 size_t f1 = iv1->indexes[x];
3210 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3216 /* Takes ownership of M. */
3218 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3219 struct index_vector *iv0, struct index_vector *iv1,
3222 if (!lvalue->n_indexes)
3224 gsl_matrix_free (lvalue->var->value);
3225 lvalue->var->value = sm;
3230 bool ok = (lvalue->n_indexes == 1
3231 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3232 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3233 free (iv0->indexes);
3234 free (iv1->indexes);
3235 gsl_matrix_free (sm);
3241 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3242 struct index_vector *iv0,
3243 struct index_vector *iv1)
3245 *iv0 = INDEX_VECTOR_INIT;
3246 *iv1 = INDEX_VECTOR_INIT;
3247 if (!lvalue->n_indexes)
3250 /* Validate destination matrix exists and has the right shape. */
3251 gsl_matrix *dm = lvalue->var->value;
3254 msg (SE, _("Undefined variable %s."), lvalue->var->name);
3257 else if (dm->size1 == 0 || dm->size2 == 0)
3259 msg (SE, _("Cannot index %zu×%zu matrix."),
3260 dm->size1, dm->size2);
3263 else if (lvalue->n_indexes == 1)
3265 if (!is_vector (dm))
3267 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
3268 dm->size1, dm->size2, lvalue->var->name);
3271 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3272 MAX (dm->size1, dm->size2),
3277 assert (lvalue->n_indexes == 2);
3278 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3279 IV_ROW, dm->size2, iv0))
3282 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3283 IV_COLUMN, dm->size1, iv1))
3285 free (iv0->indexes);
3292 /* Takes ownership of M. */
3294 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3296 struct index_vector iv0, iv1;
3297 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3299 gsl_matrix_free (sm);
3303 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3307 /* Matrix command. */
3311 struct matrix_cmd **commands;
3315 static void matrix_cmds_uninit (struct matrix_cmds *);
3319 enum matrix_cmd_type
3342 struct compute_command
3344 struct matrix_lvalue *lvalue;
3345 struct matrix_expr *rvalue;
3349 struct print_command
3351 struct matrix_expr *expression;
3352 bool use_default_format;
3353 struct fmt_spec format;
3355 int space; /* -1 means NEWPAGE. */
3359 struct string_array literals; /* CLABELS/RLABELS. */
3360 struct matrix_expr *expr; /* CNAMES/RNAMES. */
3366 struct do_if_command
3370 struct matrix_expr *condition;
3371 struct matrix_cmds commands;
3381 struct matrix_var *var;
3382 struct matrix_expr *start;
3383 struct matrix_expr *end;
3384 struct matrix_expr *increment;
3386 /* Loop conditions. */
3387 struct matrix_expr *top_condition;
3388 struct matrix_expr *bottom_condition;
3391 struct matrix_cmds commands;
3395 struct display_command
3397 struct matrix_state *state;
3403 struct release_command
3405 struct matrix_var **vars;
3412 struct matrix_expr *expression;
3413 struct save_file *sf;
3419 struct read_file *rf;
3420 struct matrix_lvalue *dst;
3421 struct matrix_expr *size;
3423 enum fmt_type format;
3430 struct write_command
3432 struct write_file *wf;
3433 struct matrix_expr *expression;
3436 /* If this is nonnull, WRITE uses this format.
3438 If this is NULL, WRITE uses free-field format with as many
3439 digits of precision as needed. */
3440 struct fmt_spec *format;
3449 struct matrix_lvalue *dst;
3450 struct dataset *dataset;
3451 struct file_handle *file;
3453 struct var_syntax *vars;
3455 struct matrix_var *names;
3457 /* Treatment of missing values. */
3462 MGET_ERROR, /* Flag error on command. */
3463 MGET_ACCEPT, /* Accept without change, user-missing only. */
3464 MGET_OMIT, /* Drop this case. */
3465 MGET_RECODE /* Recode to 'substitute'. */
3468 double substitute; /* MGET_RECODE only. */
3474 struct msave_command
3476 struct msave_common *common;
3478 struct matrix_expr *expr;
3479 const char *rowtype;
3480 struct matrix_expr *factors;
3481 struct matrix_expr *splits;
3487 struct matrix_state *state;
3488 struct file_handle *file;
3490 struct stringi_set rowtypes;
3494 struct eigen_command
3496 struct matrix_expr *expr;
3497 struct matrix_var *evec;
3498 struct matrix_var *eval;
3502 struct setdiag_command
3504 struct matrix_var *dst;
3505 struct matrix_expr *expr;
3511 struct matrix_expr *expr;
3512 struct matrix_var *u;
3513 struct matrix_var *s;
3514 struct matrix_var *v;
3520 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3521 static bool matrix_cmd_execute (struct matrix_cmd *);
3522 static void matrix_cmd_destroy (struct matrix_cmd *);
3525 static struct matrix_cmd *
3526 matrix_parse_compute (struct matrix_state *s)
3528 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3529 *cmd = (struct matrix_cmd) {
3530 .type = MCMD_COMPUTE,
3531 .compute = { .lvalue = NULL },
3534 struct compute_command *compute = &cmd->compute;
3535 compute->lvalue = matrix_lvalue_parse (s);
3536 if (!compute->lvalue)
3539 if (!lex_force_match (s->lexer, T_EQUALS))
3542 compute->rvalue = matrix_parse_expr (s);
3543 if (!compute->rvalue)
3549 matrix_cmd_destroy (cmd);
3554 matrix_cmd_execute_compute (struct compute_command *compute)
3556 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3560 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3564 print_labels_destroy (struct print_labels *labels)
3568 string_array_destroy (&labels->literals);
3569 matrix_expr_destroy (labels->expr);
3575 parse_literal_print_labels (struct matrix_state *s,
3576 struct print_labels **labelsp)
3578 lex_match (s->lexer, T_EQUALS);
3579 print_labels_destroy (*labelsp);
3580 *labelsp = xzalloc (sizeof **labelsp);
3581 while (lex_token (s->lexer) != T_SLASH
3582 && lex_token (s->lexer) != T_ENDCMD
3583 && lex_token (s->lexer) != T_STOP)
3585 struct string label = DS_EMPTY_INITIALIZER;
3586 while (lex_token (s->lexer) != T_COMMA
3587 && lex_token (s->lexer) != T_SLASH
3588 && lex_token (s->lexer) != T_ENDCMD
3589 && lex_token (s->lexer) != T_STOP)
3591 if (!ds_is_empty (&label))
3592 ds_put_byte (&label, ' ');
3594 if (lex_token (s->lexer) == T_STRING)
3595 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3598 char *rep = lex_next_representation (s->lexer, 0, 0);
3599 ds_put_cstr (&label, rep);
3604 string_array_append_nocopy (&(*labelsp)->literals,
3605 ds_steal_cstr (&label));
3607 if (!lex_match (s->lexer, T_COMMA))
3613 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3615 lex_match (s->lexer, T_EQUALS);
3616 struct matrix_expr *e = matrix_parse_exp (s);
3620 print_labels_destroy (*labelsp);
3621 *labelsp = xzalloc (sizeof **labelsp);
3622 (*labelsp)->expr = e;
3626 static struct matrix_cmd *
3627 matrix_parse_print (struct matrix_state *s)
3629 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3630 *cmd = (struct matrix_cmd) {
3633 .use_default_format = true,
3637 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3640 for (size_t i = 0; ; i++)
3642 enum token_type t = lex_next_token (s->lexer, i);
3643 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3645 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3647 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3650 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3655 cmd->print.expression = matrix_parse_exp (s);
3656 if (!cmd->print.expression)
3660 while (lex_match (s->lexer, T_SLASH))
3662 if (lex_match_id (s->lexer, "FORMAT"))
3664 lex_match (s->lexer, T_EQUALS);
3665 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3667 cmd->print.use_default_format = false;
3669 else if (lex_match_id (s->lexer, "TITLE"))
3671 lex_match (s->lexer, T_EQUALS);
3672 if (!lex_force_string (s->lexer))
3674 free (cmd->print.title);
3675 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3678 else if (lex_match_id (s->lexer, "SPACE"))
3680 lex_match (s->lexer, T_EQUALS);
3681 if (lex_match_id (s->lexer, "NEWPAGE"))
3682 cmd->print.space = -1;
3683 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3685 cmd->print.space = lex_integer (s->lexer);
3691 else if (lex_match_id (s->lexer, "RLABELS"))
3692 parse_literal_print_labels (s, &cmd->print.rlabels);
3693 else if (lex_match_id (s->lexer, "CLABELS"))
3694 parse_literal_print_labels (s, &cmd->print.clabels);
3695 else if (lex_match_id (s->lexer, "RNAMES"))
3697 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3700 else if (lex_match_id (s->lexer, "CNAMES"))
3702 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3707 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3708 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3716 matrix_cmd_destroy (cmd);
3721 matrix_is_integer (const gsl_matrix *m)
3723 for (size_t y = 0; y < m->size1; y++)
3724 for (size_t x = 0; x < m->size2; x++)
3726 double d = gsl_matrix_get (m, y, x);
3734 matrix_max_magnitude (const gsl_matrix *m)
3737 for (size_t y = 0; y < m->size1; y++)
3738 for (size_t x = 0; x < m->size2; x++)
3740 double d = fabs (gsl_matrix_get (m, y, x));
3748 format_fits (struct fmt_spec format, double x)
3750 char *s = data_out (&(union value) { .f = x }, NULL,
3751 &format, settings_get_fmt_settings ());
3752 bool fits = *s != '*' && !strchr (s, 'E');
3757 static struct fmt_spec
3758 default_format (const gsl_matrix *m, int *log_scale)
3762 double max = matrix_max_magnitude (m);
3764 if (matrix_is_integer (m))
3765 for (int w = 1; w <= 12; w++)
3767 struct fmt_spec format = { .type = FMT_F, .w = w };
3768 if (format_fits (format, -max))
3772 if (max >= 1e9 || max <= 1e-4)
3775 snprintf (s, sizeof s, "%.1e", max);
3777 const char *e = strchr (s, 'e');
3779 *log_scale = atoi (e + 1);
3782 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3786 trimmed_string (double d)
3788 struct substring s = ss_buffer ((char *) &d, sizeof d);
3789 ss_rtrim (&s, ss_cstr (" "));
3790 return ss_xstrdup (s);
3793 static struct string_array *
3794 print_labels_get (const struct print_labels *labels, size_t n,
3795 const char *prefix, bool truncate)
3800 struct string_array *out = xzalloc (sizeof *out);
3801 if (labels->literals.n)
3802 string_array_clone (out, &labels->literals);
3803 else if (labels->expr)
3805 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3806 if (m && is_vector (m))
3808 gsl_vector v = to_vector (m);
3809 for (size_t i = 0; i < v.size; i++)
3810 string_array_append_nocopy (out, trimmed_string (
3811 gsl_vector_get (&v, i)));
3813 gsl_matrix_free (m);
3819 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3822 string_array_append (out, "");
3825 string_array_delete (out, out->n - 1);
3828 for (size_t i = 0; i < out->n; i++)
3830 char *s = out->strings[i];
3831 s[strnlen (s, 8)] = '\0';
3838 matrix_cmd_print_space (int space)
3841 output_item_submit (page_break_item_create ());
3842 for (int i = 0; i < space; i++)
3843 output_log ("%s", "");
3847 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3848 struct fmt_spec format, int log_scale)
3850 matrix_cmd_print_space (print->space);
3852 output_log ("%s", print->title);
3854 output_log (" 10 ** %d X", log_scale);
3856 struct string_array *clabels = print_labels_get (print->clabels,
3857 m->size2, "col", true);
3858 if (clabels && format.w < 8)
3860 struct string_array *rlabels = print_labels_get (print->rlabels,
3861 m->size1, "row", true);
3865 struct string line = DS_EMPTY_INITIALIZER;
3867 ds_put_byte_multiple (&line, ' ', 8);
3868 for (size_t x = 0; x < m->size2; x++)
3869 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3870 output_log_nocopy (ds_steal_cstr (&line));
3873 double scale = pow (10.0, log_scale);
3874 bool numeric = fmt_is_numeric (format.type);
3875 for (size_t y = 0; y < m->size1; y++)
3877 struct string line = DS_EMPTY_INITIALIZER;
3879 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3881 for (size_t x = 0; x < m->size2; x++)
3883 double f = gsl_matrix_get (m, y, x);
3885 ? data_out (&(union value) { .f = f / scale}, NULL,
3886 &format, settings_get_fmt_settings ())
3887 : trimmed_string (f));
3888 ds_put_format (&line, " %s", s);
3891 output_log_nocopy (ds_steal_cstr (&line));
3894 string_array_destroy (rlabels);
3896 string_array_destroy (clabels);
3901 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3902 const struct print_labels *print_labels, size_t n,
3903 const char *name, const char *prefix)
3905 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3907 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3908 for (size_t i = 0; i < n; i++)
3909 pivot_category_create_leaf (
3911 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3912 : pivot_value_new_integer (i + 1)));
3914 d->hide_all_labels = true;
3915 string_array_destroy (labels);
3920 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3921 struct fmt_spec format, int log_scale)
3923 struct pivot_table *table = pivot_table_create__ (
3924 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3927 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3929 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3930 N_("Columns"), "col");
3932 struct pivot_footnote *footnote = NULL;
3935 char *s = xasprintf ("× 10**%d", log_scale);
3936 footnote = pivot_table_create_footnote (
3937 table, pivot_value_new_user_text_nocopy (s));
3940 double scale = pow (10.0, log_scale);
3941 bool numeric = fmt_is_numeric (format.type);
3942 for (size_t y = 0; y < m->size1; y++)
3943 for (size_t x = 0; x < m->size2; x++)
3945 double f = gsl_matrix_get (m, y, x);
3946 struct pivot_value *v;
3949 v = pivot_value_new_number (f / scale);
3950 v->numeric.format = format;
3953 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3955 pivot_value_add_footnote (v, footnote);
3956 pivot_table_put2 (table, y, x, v);
3959 pivot_table_submit (table);
3963 matrix_cmd_execute_print (const struct print_command *print)
3965 if (print->expression)
3967 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3972 struct fmt_spec format = (print->use_default_format
3973 ? default_format (m, &log_scale)
3976 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3977 matrix_cmd_print_text (print, m, format, log_scale);
3979 matrix_cmd_print_tables (print, m, format, log_scale);
3981 gsl_matrix_free (m);
3985 matrix_cmd_print_space (print->space);
3987 output_log ("%s", print->title);
3994 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3995 const char *command_name,
3996 const char *stop1, const char *stop2)
3998 lex_end_of_command (s->lexer);
3999 lex_discard_rest_of_command (s->lexer);
4001 size_t allocated = 0;
4004 while (lex_token (s->lexer) == T_ENDCMD)
4007 if (lex_at_phrase (s->lexer, stop1)
4008 || (stop2 && lex_at_phrase (s->lexer, stop2)))
4011 if (lex_at_phrase (s->lexer, "END MATRIX"))
4013 msg (SE, _("Premature END MATRIX within %s."), command_name);
4017 struct matrix_cmd *cmd = matrix_parse_command (s);
4021 if (c->n >= allocated)
4022 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
4023 c->commands[c->n++] = cmd;
4028 matrix_parse_do_if_clause (struct matrix_state *s,
4029 struct do_if_command *ifc,
4030 bool parse_condition,
4031 size_t *allocated_clauses)
4033 if (ifc->n_clauses >= *allocated_clauses)
4034 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
4035 sizeof *ifc->clauses);
4036 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
4037 *c = (struct do_if_clause) { .condition = NULL };
4039 if (parse_condition)
4041 c->condition = matrix_parse_expr (s);
4046 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
4049 static struct matrix_cmd *
4050 matrix_parse_do_if (struct matrix_state *s)
4052 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4053 *cmd = (struct matrix_cmd) {
4055 .do_if = { .n_clauses = 0 }
4058 size_t allocated_clauses = 0;
4061 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
4064 while (lex_match_phrase (s->lexer, "ELSE IF"));
4066 if (lex_match_id (s->lexer, "ELSE")
4067 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
4070 if (!lex_match_phrase (s->lexer, "END IF"))
4075 matrix_cmd_destroy (cmd);
4080 matrix_cmd_execute_do_if (struct do_if_command *cmd)
4082 for (size_t i = 0; i < cmd->n_clauses; i++)
4084 struct do_if_clause *c = &cmd->clauses[i];
4088 if (!matrix_expr_evaluate_scalar (c->condition,
4089 i ? "ELSE IF" : "DO IF",
4094 for (size_t j = 0; j < c->commands.n; j++)
4095 if (!matrix_cmd_execute (c->commands.commands[j]))
4102 static struct matrix_cmd *
4103 matrix_parse_loop (struct matrix_state *s)
4105 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4106 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4108 struct loop_command *loop = &cmd->loop;
4109 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4111 struct substring name = lex_tokss (s->lexer);
4112 loop->var = matrix_var_lookup (s, name);
4114 loop->var = matrix_var_create (s, name);
4119 loop->start = matrix_parse_expr (s);
4120 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4123 loop->end = matrix_parse_expr (s);
4127 if (lex_match (s->lexer, T_BY))
4129 loop->increment = matrix_parse_expr (s);
4130 if (!loop->increment)
4135 if (lex_match_id (s->lexer, "IF"))
4137 loop->top_condition = matrix_parse_expr (s);
4138 if (!loop->top_condition)
4142 bool was_in_loop = s->in_loop;
4144 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4146 s->in_loop = was_in_loop;
4150 if (!lex_match_phrase (s->lexer, "END LOOP"))
4153 if (lex_match_id (s->lexer, "IF"))
4155 loop->bottom_condition = matrix_parse_expr (s);
4156 if (!loop->bottom_condition)
4163 matrix_cmd_destroy (cmd);
4168 matrix_cmd_execute_loop (struct loop_command *cmd)
4172 long int increment = 1;
4175 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4176 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4178 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4182 if (increment > 0 ? value > end
4183 : increment < 0 ? value < end
4188 int mxloops = settings_get_mxloops ();
4189 for (int i = 0; i < mxloops; i++)
4194 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4196 gsl_matrix_free (cmd->var->value);
4197 cmd->var->value = NULL;
4199 if (!cmd->var->value)
4200 cmd->var->value = gsl_matrix_alloc (1, 1);
4202 gsl_matrix_set (cmd->var->value, 0, 0, value);
4205 if (cmd->top_condition)
4208 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4213 for (size_t j = 0; j < cmd->commands.n; j++)
4214 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4217 if (cmd->bottom_condition)
4220 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4221 "END LOOP IF", &d) || d > 0)
4228 if (increment > 0 ? value > end : value < end)
4234 static struct matrix_cmd *
4235 matrix_parse_break (struct matrix_state *s)
4239 msg (SE, _("BREAK not inside LOOP."));
4243 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4244 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4248 static struct matrix_cmd *
4249 matrix_parse_display (struct matrix_state *s)
4251 bool dictionary = false;
4252 bool status = false;
4253 while (lex_token (s->lexer) == T_ID)
4255 if (lex_match_id (s->lexer, "DICTIONARY"))
4257 else if (lex_match_id (s->lexer, "STATUS"))
4261 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4265 if (!dictionary && !status)
4266 dictionary = status = true;
4268 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4269 *cmd = (struct matrix_cmd) {
4270 .type = MCMD_DISPLAY,
4273 .dictionary = dictionary,
4281 compare_matrix_var_pointers (const void *a_, const void *b_)
4283 const struct matrix_var *const *ap = a_;
4284 const struct matrix_var *const *bp = b_;
4285 const struct matrix_var *a = *ap;
4286 const struct matrix_var *b = *bp;
4287 return strcmp (a->name, b->name);
4291 matrix_cmd_execute_display (struct display_command *cmd)
4293 const struct matrix_state *s = cmd->state;
4295 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4296 struct pivot_dimension *attr_dimension
4297 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
4298 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
4299 N_("Rows"), N_("Columns"));
4300 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
4302 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4305 const struct matrix_var *var;
4306 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4308 vars[n_vars++] = var;
4309 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4311 struct pivot_dimension *rows = pivot_dimension_create (
4312 table, PIVOT_AXIS_ROW, N_("Variable"));
4313 for (size_t i = 0; i < n_vars; i++)
4315 const struct matrix_var *var = vars[i];
4316 pivot_category_create_leaf (
4317 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4319 size_t r = var->value->size1;
4320 size_t c = var->value->size2;
4321 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4322 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4323 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4326 pivot_table_submit (table);
4329 static struct matrix_cmd *
4330 matrix_parse_release (struct matrix_state *s)
4332 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4333 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4335 size_t allocated_vars = 0;
4336 while (lex_token (s->lexer) == T_ID)
4338 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4341 if (cmd->release.n_vars >= allocated_vars)
4342 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4343 sizeof *cmd->release.vars);
4344 cmd->release.vars[cmd->release.n_vars++] = var;
4347 lex_error (s->lexer, _("Variable name expected."));
4350 if (!lex_match (s->lexer, T_COMMA))
4358 matrix_cmd_execute_release (struct release_command *cmd)
4360 for (size_t i = 0; i < cmd->n_vars; i++)
4362 struct matrix_var *var = cmd->vars[i];
4363 gsl_matrix_free (var->value);
4368 static struct save_file *
4369 save_file_create (struct matrix_state *s, struct file_handle *fh,
4370 struct string_array *variables,
4371 struct matrix_expr *names,
4372 struct stringi_set *strings)
4374 for (size_t i = 0; i < s->n_save_files; i++)
4376 struct save_file *sf = s->save_files[i];
4377 if (fh_equal (sf->file, fh))
4381 string_array_destroy (variables);
4382 matrix_expr_destroy (names);
4383 stringi_set_destroy (strings);
4389 struct save_file *sf = xmalloc (sizeof *sf);
4390 *sf = (struct save_file) {
4392 .dataset = s->dataset,
4393 .variables = *variables,
4395 .strings = STRINGI_SET_INITIALIZER (sf->strings),
4398 stringi_set_swap (&sf->strings, strings);
4399 stringi_set_destroy (strings);
4401 s->save_files = xrealloc (s->save_files,
4402 (s->n_save_files + 1) * sizeof *s->save_files);
4403 s->save_files[s->n_save_files++] = sf;
4407 static struct casewriter *
4408 save_file_open (struct save_file *sf, gsl_matrix *m)
4410 if (sf->writer || sf->error)
4414 size_t n_variables = caseproto_get_n_widths (
4415 casewriter_get_proto (sf->writer));
4416 if (m->size2 != n_variables)
4418 msg (ME, _("The first SAVE to %s within this matrix program "
4419 "had %zu columns, so a %zu×%zu matrix cannot be "
4421 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
4422 n_variables, m->size1, m->size2);
4430 struct dictionary *dict = dict_create (get_default_encoding ());
4432 /* Fill 'names' with user-specified names if there were any, either from
4433 sf->variables or sf->names. */
4434 struct string_array names = { .n = 0 };
4435 if (sf->variables.n)
4436 string_array_clone (&names, &sf->variables);
4439 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
4440 if (nm && is_vector (nm))
4442 gsl_vector nv = to_vector (nm);
4443 for (size_t i = 0; i < nv.size; i++)
4445 char *name = trimmed_string (gsl_vector_get (&nv, i));
4446 if (dict_id_is_valid (dict, name, true))
4447 string_array_append_nocopy (&names, name);
4452 gsl_matrix_free (nm);
4455 struct stringi_set strings;
4456 stringi_set_clone (&strings, &sf->strings);
4458 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4463 name = names.strings[i];
4466 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4470 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4471 struct variable *var = dict_create_var (dict, name, width);
4474 msg (ME, _("Duplicate variable name %s in SAVE statement."),
4480 if (!stringi_set_is_empty (&strings))
4482 size_t count = stringi_set_count (&strings);
4483 const char *example = stringi_set_node_get_string (
4484 stringi_set_first (&strings));
4487 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
4488 "variable %s."), example);
4490 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
4491 "unknown variable: %s.",
4492 "The SAVE command STRINGS subcommand specifies %zu "
4493 "unknown variables, including %s.",
4498 stringi_set_destroy (&strings);
4499 string_array_destroy (&names);
4508 if (sf->file == fh_inline_file ())
4509 sf->writer = autopaging_writer_create (dict_get_proto (dict));
4511 sf->writer = any_writer_open (sf->file, dict);
4523 save_file_destroy (struct save_file *sf)
4527 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
4529 dataset_set_dict (sf->dataset, sf->dict);
4530 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
4534 casewriter_destroy (sf->writer);
4535 dict_unref (sf->dict);
4537 fh_unref (sf->file);
4538 string_array_destroy (&sf->variables);
4539 matrix_expr_destroy (sf->names);
4540 stringi_set_destroy (&sf->strings);
4545 static struct matrix_cmd *
4546 matrix_parse_save (struct matrix_state *s)
4548 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4549 *cmd = (struct matrix_cmd) {
4551 .save = { .expression = NULL },
4554 struct file_handle *fh = NULL;
4555 struct save_command *save = &cmd->save;
4557 struct string_array variables = STRING_ARRAY_INITIALIZER;
4558 struct matrix_expr *names = NULL;
4559 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
4561 save->expression = matrix_parse_exp (s);
4562 if (!save->expression)
4565 while (lex_match (s->lexer, T_SLASH))
4567 if (lex_match_id (s->lexer, "OUTFILE"))
4569 lex_match (s->lexer, T_EQUALS);
4572 fh = (lex_match (s->lexer, T_ASTERISK)
4574 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4578 else if (lex_match_id (s->lexer, "VARIABLES"))
4580 lex_match (s->lexer, T_EQUALS);
4584 struct dictionary *d = dict_create (get_default_encoding ());
4585 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4586 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4591 string_array_clear (&variables);
4592 variables = (struct string_array) {
4598 else if (lex_match_id (s->lexer, "NAMES"))
4600 lex_match (s->lexer, T_EQUALS);
4601 matrix_expr_destroy (names);
4602 names = matrix_parse_exp (s);
4606 else if (lex_match_id (s->lexer, "STRINGS"))
4608 lex_match (s->lexer, T_EQUALS);
4609 while (lex_token (s->lexer) == T_ID)
4611 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
4613 if (!lex_match (s->lexer, T_COMMA))
4619 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4627 if (s->prev_save_file)
4628 fh = fh_ref (s->prev_save_file);
4631 lex_sbc_missing ("OUTFILE");
4635 fh_unref (s->prev_save_file);
4636 s->prev_save_file = fh_ref (fh);
4638 if (variables.n && names)
4640 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4641 matrix_expr_destroy (names);
4645 save->sf = save_file_create (s, fh, &variables, names, &strings);
4649 string_array_destroy (&variables);
4650 matrix_expr_destroy (names);
4651 stringi_set_destroy (&strings);
4653 matrix_cmd_destroy (cmd);
4658 matrix_cmd_execute_save (const struct save_command *save)
4660 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4664 struct casewriter *writer = save_file_open (save->sf, m);
4667 gsl_matrix_free (m);
4671 const struct caseproto *proto = casewriter_get_proto (writer);
4672 for (size_t y = 0; y < m->size1; y++)
4674 struct ccase *c = case_create (proto);
4675 for (size_t x = 0; x < m->size2; x++)
4677 double d = gsl_matrix_get (m, y, x);
4678 int width = caseproto_get_width (proto, x);
4679 union value *value = case_data_rw_idx (c, x);
4683 memcpy (value->s, &d, width);
4685 casewriter_write (writer, c);
4687 gsl_matrix_free (m);
4690 static struct matrix_cmd *
4691 matrix_parse_read (struct matrix_state *s)
4693 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4694 *cmd = (struct matrix_cmd) {
4696 .read = { .format = FMT_F },
4699 struct file_handle *fh = NULL;
4700 char *encoding = NULL;
4701 struct read_command *read = &cmd->read;
4702 read->dst = matrix_lvalue_parse (s);
4707 int repetitions = 0;
4708 int record_width = 0;
4709 bool seen_format = false;
4710 while (lex_match (s->lexer, T_SLASH))
4712 if (lex_match_id (s->lexer, "FILE"))
4714 lex_match (s->lexer, T_EQUALS);
4717 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4721 else if (lex_match_id (s->lexer, "ENCODING"))
4723 lex_match (s->lexer, T_EQUALS);
4724 if (!lex_force_string (s->lexer))
4728 encoding = ss_xstrdup (lex_tokss (s->lexer));
4732 else if (lex_match_id (s->lexer, "FIELD"))
4734 lex_match (s->lexer, T_EQUALS);
4736 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4738 read->c1 = lex_integer (s->lexer);
4740 if (!lex_force_match (s->lexer, T_TO)
4741 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4743 read->c2 = lex_integer (s->lexer) + 1;
4746 record_width = read->c2 - read->c1;
4747 if (lex_match (s->lexer, T_BY))
4749 if (!lex_force_int_range (s->lexer, "BY", 1,
4750 read->c2 - read->c1))
4752 by = lex_integer (s->lexer);
4755 if (record_width % by)
4757 msg (SE, _("BY %d does not evenly divide record width %d."),
4765 else if (lex_match_id (s->lexer, "SIZE"))
4767 lex_match (s->lexer, T_EQUALS);
4768 matrix_expr_destroy (read->size);
4769 read->size = matrix_parse_exp (s);
4773 else if (lex_match_id (s->lexer, "MODE"))
4775 lex_match (s->lexer, T_EQUALS);
4776 if (lex_match_id (s->lexer, "RECTANGULAR"))
4777 read->symmetric = false;
4778 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4779 read->symmetric = true;
4782 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4786 else if (lex_match_id (s->lexer, "REREAD"))
4787 read->reread = true;
4788 else if (lex_match_id (s->lexer, "FORMAT"))
4792 lex_sbc_only_once ("FORMAT");
4797 lex_match (s->lexer, T_EQUALS);
4799 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4802 const char *p = lex_tokcstr (s->lexer);
4803 if (c_isdigit (p[0]))
4805 repetitions = atoi (p);
4806 p += strspn (p, "0123456789");
4807 if (!fmt_from_name (p, &read->format))
4809 lex_error (s->lexer, _("Unknown format %s."), p);
4814 else if (fmt_from_name (p, &read->format))
4818 struct fmt_spec format;
4819 if (!parse_format_specifier (s->lexer, &format))
4821 read->format = format.type;
4827 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4828 "REREAD", "FORMAT");
4835 lex_sbc_missing ("FIELD");
4839 if (!read->dst->n_indexes && !read->size)
4841 msg (SE, _("SIZE is required for reading data into a full matrix "
4842 "(as opposed to a submatrix)."));
4848 if (s->prev_read_file)
4849 fh = fh_ref (s->prev_read_file);
4852 lex_sbc_missing ("FILE");
4856 fh_unref (s->prev_read_file);
4857 s->prev_read_file = fh_ref (fh);
4859 read->rf = read_file_create (s, fh);
4863 free (read->rf->encoding);
4864 read->rf->encoding = encoding;
4868 /* Field width may be specified in multiple ways:
4871 2. The format on FORMAT.
4872 3. The repetition factor on FORMAT.
4874 (2) and (3) are mutually exclusive.
4876 If more than one of these is present, they must agree. If none of them is
4877 present, then free-field format is used.
4879 if (repetitions > record_width)
4881 msg (SE, _("%d repetitions cannot fit in record width %d."),
4882 repetitions, record_width);
4885 int w = (repetitions ? record_width / repetitions
4891 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4892 "which implies field width %d, "
4893 "but BY specifies field width %d."),
4894 repetitions, record_width, w, by);
4896 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4905 matrix_cmd_destroy (cmd);
4911 parse_error (const struct dfm_reader *reader, enum fmt_type format,
4912 struct substring data, size_t y, size_t x,
4913 int first_column, int last_column, char *error)
4915 int line_number = dfm_get_line_number (reader);
4916 struct msg_location *location = xmalloc (sizeof *location);
4917 *location = (struct msg_location) {
4918 .file_name = xstrdup (dfm_get_file_name (reader)),
4919 .first_line = line_number,
4920 .last_line = line_number + 1,
4921 .first_column = first_column,
4922 .last_column = last_column,
4924 struct msg *m = xmalloc (sizeof *m);
4926 .category = MSG_C_DATA,
4927 .severity = MSG_S_WARNING,
4928 .location = location,
4929 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
4930 "for matrix row %zu, column %zu: %s"),
4931 (int) data.length, data.string, fmt_name (format),
4932 y + 1, x + 1, error),
4940 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4941 gsl_matrix *m, struct substring p, size_t y, size_t x,
4942 const char *line_start)
4944 const char *input_encoding = dfm_reader_get_encoding (reader);
4947 if (fmt_is_numeric (read->format))
4950 error = data_in (p, input_encoding, read->format,
4951 settings_get_fmt_settings (), &v, 0, NULL);
4952 if (!error && v.f == SYSMIS)
4953 error = xstrdup (_("Matrix data may not contain missing value."));
4958 uint8_t s[sizeof (double)];
4959 union value v = { .s = s };
4960 error = data_in (p, input_encoding, read->format,
4961 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
4962 memcpy (&f, s, sizeof f);
4967 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
4968 int c2 = c1 + ss_utf8_count_columns (p) - 1;
4969 parse_error (reader, read->format, p, y, x, c1, c2, error);
4973 gsl_matrix_set (m, y, x, f);
4974 if (read->symmetric && x != y)
4975 gsl_matrix_set (m, x, y, f);
4980 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4981 struct substring *line, const char **startp)
4983 if (dfm_eof (reader))
4985 msg (SE, _("Unexpected end of file reading matrix data."));
4988 dfm_expand_tabs (reader);
4989 struct substring record = dfm_get_record (reader);
4990 /* XXX need to recode record into UTF-8 */
4991 *startp = record.string;
4992 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
4997 matrix_read (struct read_command *read, struct dfm_reader *reader,
5000 for (size_t y = 0; y < m->size1; y++)
5002 size_t nx = read->symmetric ? y + 1 : m->size2;
5004 struct substring line = ss_empty ();
5005 const char *line_start = line.string;
5006 for (size_t x = 0; x < nx; x++)
5013 ss_ltrim (&line, ss_cstr (" ,"));
5014 if (!ss_is_empty (line))
5016 if (!matrix_read_line (read, reader, &line, &line_start))
5018 dfm_forward_record (reader);
5021 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
5025 if (!matrix_read_line (read, reader, &line, &line_start))
5027 size_t fields_per_line = (read->c2 - read->c1) / read->w;
5028 int f = x % fields_per_line;
5029 if (f == fields_per_line - 1)
5030 dfm_forward_record (reader);
5032 p = ss_substr (line, read->w * f, read->w);
5035 matrix_read_set_field (read, reader, m, p, y, x, line_start);
5039 dfm_forward_record (reader);
5042 ss_ltrim (&line, ss_cstr (" ,"));
5043 if (!ss_is_empty (line))
5046 msg (SW, _("Trailing garbage on line \"%.*s\""),
5047 (int) line.length, line.string);
5054 matrix_cmd_execute_read (struct read_command *read)
5056 struct index_vector iv0, iv1;
5057 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
5060 size_t size[2] = { SIZE_MAX, SIZE_MAX };
5063 gsl_matrix *m = matrix_expr_evaluate (read->size);
5069 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5070 "not a %zu×%zu matrix."), m->size1, m->size2);
5071 gsl_matrix_free (m);
5077 gsl_vector v = to_vector (m);
5081 d[0] = gsl_vector_get (&v, 0);
5084 else if (v.size == 2)
5086 d[0] = gsl_vector_get (&v, 0);
5087 d[1] = gsl_vector_get (&v, 1);
5091 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
5092 "not a %zu×%zu matrix."),
5093 m->size1, m->size2),
5094 gsl_matrix_free (m);
5099 gsl_matrix_free (m);
5101 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
5103 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
5104 "are outside valid range."),
5115 if (read->dst->n_indexes)
5117 size_t submatrix_size[2];
5118 if (read->dst->n_indexes == 2)
5120 submatrix_size[0] = iv0.n;
5121 submatrix_size[1] = iv1.n;
5123 else if (read->dst->var->value->size1 == 1)
5125 submatrix_size[0] = 1;
5126 submatrix_size[1] = iv0.n;
5130 submatrix_size[0] = iv0.n;
5131 submatrix_size[1] = 1;
5136 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
5138 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
5139 "differ from submatrix dimensions %zu×%zu."),
5141 submatrix_size[0], submatrix_size[1]);
5149 size[0] = submatrix_size[0];
5150 size[1] = submatrix_size[1];
5154 struct dfm_reader *reader = read_file_open (read->rf);
5156 dfm_reread_record (reader, 1);
5158 if (read->symmetric && size[0] != size[1])
5160 msg (SE, _("Cannot read non-square %zu×%zu matrix "
5161 "using READ with MODE=SYMMETRIC."),
5167 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
5168 matrix_read (read, reader, tmp);
5169 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
5172 static struct matrix_cmd *
5173 matrix_parse_write (struct matrix_state *s)
5175 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5176 *cmd = (struct matrix_cmd) {
5180 struct file_handle *fh = NULL;
5181 char *encoding = NULL;
5182 struct write_command *write = &cmd->write;
5183 write->expression = matrix_parse_exp (s);
5184 if (!write->expression)
5188 int repetitions = 0;
5189 int record_width = 0;
5190 enum fmt_type format = FMT_F;
5191 bool has_format = false;
5192 while (lex_match (s->lexer, T_SLASH))
5194 if (lex_match_id (s->lexer, "OUTFILE"))
5196 lex_match (s->lexer, T_EQUALS);
5199 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5203 else if (lex_match_id (s->lexer, "ENCODING"))
5205 lex_match (s->lexer, T_EQUALS);
5206 if (!lex_force_string (s->lexer))
5210 encoding = ss_xstrdup (lex_tokss (s->lexer));
5214 else if (lex_match_id (s->lexer, "FIELD"))
5216 lex_match (s->lexer, T_EQUALS);
5218 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5220 write->c1 = lex_integer (s->lexer);
5222 if (!lex_force_match (s->lexer, T_TO)
5223 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5225 write->c2 = lex_integer (s->lexer) + 1;
5228 record_width = write->c2 - write->c1;
5229 if (lex_match (s->lexer, T_BY))
5231 if (!lex_force_int_range (s->lexer, "BY", 1,
5232 write->c2 - write->c1))
5234 by = lex_integer (s->lexer);
5237 if (record_width % by)
5239 msg (SE, _("BY %d does not evenly divide record width %d."),
5247 else if (lex_match_id (s->lexer, "MODE"))
5249 lex_match (s->lexer, T_EQUALS);
5250 if (lex_match_id (s->lexer, "RECTANGULAR"))
5251 write->triangular = false;
5252 else if (lex_match_id (s->lexer, "TRIANGULAR"))
5253 write->triangular = true;
5256 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5260 else if (lex_match_id (s->lexer, "HOLD"))
5262 else if (lex_match_id (s->lexer, "FORMAT"))
5264 if (has_format || write->format)
5266 lex_sbc_only_once ("FORMAT");
5270 lex_match (s->lexer, T_EQUALS);
5272 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5275 const char *p = lex_tokcstr (s->lexer);
5276 if (c_isdigit (p[0]))
5278 repetitions = atoi (p);
5279 p += strspn (p, "0123456789");
5280 if (!fmt_from_name (p, &format))
5282 lex_error (s->lexer, _("Unknown format %s."), p);
5288 else if (fmt_from_name (p, &format))
5295 struct fmt_spec spec;
5296 if (!parse_format_specifier (s->lexer, &spec))
5298 write->format = xmemdup (&spec, sizeof spec);
5303 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5311 lex_sbc_missing ("FIELD");
5317 if (s->prev_write_file)
5318 fh = fh_ref (s->prev_write_file);
5321 lex_sbc_missing ("OUTFILE");
5325 fh_unref (s->prev_write_file);
5326 s->prev_write_file = fh_ref (fh);
5328 write->wf = write_file_create (s, fh);
5332 free (write->wf->encoding);
5333 write->wf->encoding = encoding;
5337 /* Field width may be specified in multiple ways:
5340 2. The format on FORMAT.
5341 3. The repetition factor on FORMAT.
5343 (2) and (3) are mutually exclusive.
5345 If more than one of these is present, they must agree. If none of them is
5346 present, then free-field format is used.
5348 if (repetitions > record_width)
5350 msg (SE, _("%d repetitions cannot fit in record width %d."),
5351 repetitions, record_width);
5354 int w = (repetitions ? record_width / repetitions
5355 : write->format ? write->format->w
5360 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5361 "which implies field width %d, "
5362 "but BY specifies field width %d."),
5363 repetitions, record_width, w, by);
5365 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5369 if (w && !write->format)
5371 write->format = xmalloc (sizeof *write->format);
5372 *write->format = (struct fmt_spec) { .type = format, .w = w };
5374 if (!fmt_check_output (write->format))
5378 if (write->format && fmt_var_width (write->format) > sizeof (double))
5380 char s[FMT_STRING_LEN_MAX + 1];
5381 fmt_to_string (write->format, s);
5382 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
5383 s, sizeof (double));
5391 matrix_cmd_destroy (cmd);
5396 matrix_cmd_execute_write (struct write_command *write)
5398 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5402 if (write->triangular && m->size1 != m->size2)
5404 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5405 "the matrix to be written has dimensions %zu×%zu."),
5406 m->size1, m->size2);
5407 gsl_matrix_free (m);
5411 struct dfm_writer *writer = write_file_open (write->wf);
5412 if (!writer || !m->size1)
5414 gsl_matrix_free (m);
5418 const struct fmt_settings *settings = settings_get_fmt_settings ();
5419 struct u8_line *line = write->wf->held;
5420 for (size_t y = 0; y < m->size1; y++)
5424 line = xmalloc (sizeof *line);
5425 u8_line_init (line);
5427 size_t nx = write->triangular ? y + 1 : m->size2;
5429 for (size_t x = 0; x < nx; x++)
5432 double f = gsl_matrix_get (m, y, x);
5436 if (fmt_is_numeric (write->format->type))
5439 v.s = (uint8_t *) &f;
5440 s = data_out (&v, NULL, write->format, settings);
5444 s = xmalloc (DBL_BUFSIZE_BOUND);
5445 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
5446 >= DBL_BUFSIZE_BOUND)
5449 size_t len = strlen (s);
5450 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5451 if (width + x0 > write->c2)
5453 dfm_put_record_utf8 (writer, line->s.ss.string,
5455 u8_line_clear (line);
5458 u8_line_put (line, x0, x0 + width, s, len);
5461 x0 += write->format ? write->format->w : width + 1;
5464 if (y + 1 >= m->size1 && write->hold)
5466 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5467 u8_line_clear (line);
5471 u8_line_destroy (line);
5475 write->wf->held = line;
5477 gsl_matrix_free (m);
5480 static struct matrix_cmd *
5481 matrix_parse_get (struct matrix_state *s)
5483 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5484 *cmd = (struct matrix_cmd) {
5487 .dataset = s->dataset,
5488 .user = { .treatment = MGET_ERROR },
5489 .system = { .treatment = MGET_ERROR },
5493 struct get_command *get = &cmd->get;
5494 get->dst = matrix_lvalue_parse (s);
5498 while (lex_match (s->lexer, T_SLASH))
5500 if (lex_match_id (s->lexer, "FILE"))
5502 lex_match (s->lexer, T_EQUALS);
5504 fh_unref (get->file);
5505 if (lex_match (s->lexer, T_ASTERISK))
5509 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5514 else if (lex_match_id (s->lexer, "ENCODING"))
5516 lex_match (s->lexer, T_EQUALS);
5517 if (!lex_force_string (s->lexer))
5520 free (get->encoding);
5521 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5525 else if (lex_match_id (s->lexer, "VARIABLES"))
5527 lex_match (s->lexer, T_EQUALS);
5531 lex_sbc_only_once ("VARIABLES");
5535 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
5538 else if (lex_match_id (s->lexer, "NAMES"))
5540 lex_match (s->lexer, T_EQUALS);
5541 if (!lex_force_id (s->lexer))
5544 struct substring name = lex_tokss (s->lexer);
5545 get->names = matrix_var_lookup (s, name);
5547 get->names = matrix_var_create (s, name);
5550 else if (lex_match_id (s->lexer, "MISSING"))
5552 lex_match (s->lexer, T_EQUALS);
5553 if (lex_match_id (s->lexer, "ACCEPT"))
5554 get->user.treatment = MGET_ACCEPT;
5555 else if (lex_match_id (s->lexer, "OMIT"))
5556 get->user.treatment = MGET_OMIT;
5557 else if (lex_is_number (s->lexer))
5559 get->user.treatment = MGET_RECODE;
5560 get->user.substitute = lex_number (s->lexer);
5565 lex_error (s->lexer, NULL);
5569 else if (lex_match_id (s->lexer, "SYSMIS"))
5571 lex_match (s->lexer, T_EQUALS);
5572 if (lex_match_id (s->lexer, "OMIT"))
5573 get->system.treatment = MGET_OMIT;
5574 else if (lex_is_number (s->lexer))
5576 get->system.treatment = MGET_RECODE;
5577 get->system.substitute = lex_number (s->lexer);
5582 lex_error (s->lexer, NULL);
5588 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5589 "MISSING", "SYSMIS");
5594 if (get->user.treatment != MGET_ACCEPT)
5595 get->system.treatment = MGET_ERROR;
5600 matrix_cmd_destroy (cmd);
5605 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
5606 const struct dictionary *dict)
5608 struct variable **vars;
5613 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
5614 &vars, &n_vars, PV_NUMERIC))
5619 n_vars = dict_get_var_cnt (dict);
5620 vars = xnmalloc (n_vars, sizeof *vars);
5621 for (size_t i = 0; i < n_vars; i++)
5623 struct variable *var = dict_get_var (dict, i);
5624 if (!var_is_numeric (var))
5626 msg (SE, _("GET: Variable %s is not numeric."),
5627 var_get_name (var));
5637 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
5638 for (size_t i = 0; i < n_vars; i++)
5640 char s[sizeof (double)];
5642 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
5643 memcpy (&f, s, sizeof f);
5644 gsl_matrix_set (names, i, 0, f);
5647 gsl_matrix_free (get->names->value);
5648 get->names->value = names;
5652 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5653 long long int casenum = 1;
5655 for (struct ccase *c = casereader_read (reader); c;
5656 c = casereader_read (reader), casenum++)
5658 if (n_rows >= m->size1)
5660 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5661 for (size_t y = 0; y < n_rows; y++)
5662 for (size_t x = 0; x < n_vars; x++)
5663 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5664 gsl_matrix_free (m);
5669 for (size_t x = 0; x < n_vars; x++)
5671 const struct variable *var = vars[x];
5672 double d = case_num (c, var);
5675 if (get->system.treatment == MGET_RECODE)
5676 d = get->system.substitute;
5677 else if (get->system.treatment == MGET_OMIT)
5681 msg (SE, _("GET: Variable %s in case %lld "
5682 "is system-missing."),
5683 var_get_name (var), casenum);
5687 else if (var_is_num_missing (var, d, MV_USER))
5689 if (get->user.treatment == MGET_RECODE)
5690 d = get->user.substitute;
5691 else if (get->user.treatment == MGET_OMIT)
5693 else if (get->user.treatment != MGET_ACCEPT)
5695 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5697 var_get_name (var), casenum, d);
5701 gsl_matrix_set (m, n_rows, x, d);
5712 matrix_lvalue_evaluate_and_assign (get->dst, m);
5715 gsl_matrix_free (m);
5720 matrix_open_casereader (const char *command_name,
5721 struct file_handle *file, const char *encoding,
5722 struct dataset *dataset,
5723 struct casereader **readerp, struct dictionary **dictp)
5727 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
5728 return *readerp != NULL;
5732 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
5734 msg (ME, _("The %s command cannot read empty active file."),
5738 *readerp = proc_open (dataset);
5739 *dictp = dict_ref (dataset_dict (dataset));
5745 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
5746 struct casereader *reader, struct dictionary *dict)
5749 casereader_destroy (reader);
5751 proc_commit (dataset);
5755 matrix_cmd_execute_get (struct get_command *get)
5757 struct casereader *r;
5758 struct dictionary *d;
5759 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
5762 matrix_cmd_execute_get__ (get, r, d);
5763 matrix_close_casereader (get->file, get->dataset, r, d);
5768 match_rowtype (struct lexer *lexer)
5770 static const char *rowtypes[] = {
5771 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5773 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5775 for (size_t i = 0; i < n_rowtypes; i++)
5776 if (lex_match_id (lexer, rowtypes[i]))
5779 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5784 parse_var_names (struct lexer *lexer, struct string_array *sa)
5786 lex_match (lexer, T_EQUALS);
5788 string_array_clear (sa);
5790 struct dictionary *dict = dict_create (get_default_encoding ());
5791 dict_create_var_assert (dict, "ROWTYPE_", 8);
5792 dict_create_var_assert (dict, "VARNAME_", 8);
5795 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5796 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5801 string_array_clear (sa);
5802 sa->strings = names;
5803 sa->n = sa->allocated = n_names;
5809 msave_common_uninit (struct msave_common *common)
5813 fh_unref (common->outfile);
5814 string_array_destroy (&common->variables);
5815 string_array_destroy (&common->fnames);
5816 string_array_destroy (&common->snames);
5821 compare_variables (const char *keyword,
5822 const struct string_array *new,
5823 const struct string_array *old)
5829 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5830 "on the first MSAVE within MATRIX."), keyword);
5833 else if (!string_array_equal_case (old, new))
5835 msg (SE, _("%s must specify the same variables each time within "
5836 "a given MATRIX."), keyword);
5843 static struct matrix_cmd *
5844 matrix_parse_msave (struct matrix_state *s)
5846 struct msave_common common = { .outfile = NULL };
5847 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5848 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5850 struct msave_command *msave = &cmd->msave;
5851 if (lex_token (s->lexer) == T_ID)
5852 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5853 msave->expr = matrix_parse_exp (s);
5857 while (lex_match (s->lexer, T_SLASH))
5859 if (lex_match_id (s->lexer, "TYPE"))
5861 lex_match (s->lexer, T_EQUALS);
5863 msave->rowtype = match_rowtype (s->lexer);
5864 if (!msave->rowtype)
5867 else if (lex_match_id (s->lexer, "OUTFILE"))
5869 lex_match (s->lexer, T_EQUALS);
5871 fh_unref (common.outfile);
5872 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5873 if (!common.outfile)
5876 else if (lex_match_id (s->lexer, "VARIABLES"))
5878 if (!parse_var_names (s->lexer, &common.variables))
5881 else if (lex_match_id (s->lexer, "FNAMES"))
5883 if (!parse_var_names (s->lexer, &common.fnames))
5886 else if (lex_match_id (s->lexer, "SNAMES"))
5888 if (!parse_var_names (s->lexer, &common.snames))
5891 else if (lex_match_id (s->lexer, "SPLIT"))
5893 lex_match (s->lexer, T_EQUALS);
5895 matrix_expr_destroy (msave->splits);
5896 msave->splits = matrix_parse_exp (s);
5900 else if (lex_match_id (s->lexer, "FACTOR"))
5902 lex_match (s->lexer, T_EQUALS);
5904 matrix_expr_destroy (msave->factors);
5905 msave->factors = matrix_parse_exp (s);
5906 if (!msave->factors)
5911 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5912 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5916 if (!msave->rowtype)
5918 lex_sbc_missing ("TYPE");
5921 common.has_splits = msave->splits || common.snames.n;
5922 common.has_factors = msave->factors || common.fnames.n;
5924 struct msave_common *c = s->common ? s->common : &common;
5925 if (c->fnames.n && !msave->factors)
5927 msg (SE, _("FNAMES requires FACTOR."));
5930 if (c->snames.n && !msave->splits)
5932 msg (SE, _("SNAMES requires SPLIT."));
5935 if (c->has_factors && !common.has_factors)
5937 msg (SE, _("%s is required because it was present on the first "
5938 "MSAVE in this MATRIX command."), "FACTOR");
5941 if (c->has_splits && !common.has_splits)
5943 msg (SE, _("%s is required because it was present on the first "
5944 "MSAVE in this MATRIX command."), "SPLIT");
5950 if (!common.outfile)
5952 lex_sbc_missing ("OUTFILE");
5955 s->common = xmemdup (&common, sizeof common);
5961 bool same = common.outfile == s->common->outfile;
5962 fh_unref (common.outfile);
5965 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5966 "within a single MATRIX command."));
5970 if (!compare_variables ("VARIABLES",
5971 &common.variables, &s->common->variables)
5972 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5973 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5975 msave_common_uninit (&common);
5977 msave->common = s->common;
5978 if (!msave->varname_)
5979 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5983 msave_common_uninit (&common);
5984 matrix_cmd_destroy (cmd);
5989 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5991 gsl_matrix *m = matrix_expr_evaluate (e);
5997 msg (SE, _("%s expression must evaluate to vector, "
5998 "not a %zu×%zu matrix."),
5999 name, m->size1, m->size2);
6000 gsl_matrix_free (m);
6004 return matrix_to_vector (m);
6008 msave_add_vars (struct dictionary *d, const struct string_array *vars)
6010 for (size_t i = 0; i < vars->n; i++)
6011 if (!dict_create_var (d, vars->strings[i], 0))
6012 return vars->strings[i];
6016 static struct dictionary *
6017 msave_create_dict (const struct msave_common *common)
6019 struct dictionary *dict = dict_create (get_default_encoding ());
6021 const char *dup_split = msave_add_vars (dict, &common->snames);
6024 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
6028 dict_create_var_assert (dict, "ROWTYPE_", 8);
6030 const char *dup_factor = msave_add_vars (dict, &common->fnames);
6033 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
6037 dict_create_var_assert (dict, "VARNAME_", 8);
6039 const char *dup_var = msave_add_vars (dict, &common->variables);
6042 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
6054 matrix_cmd_execute_msave (struct msave_command *msave)
6056 struct msave_common *common = msave->common;
6057 gsl_matrix *m = NULL;
6058 gsl_vector *factors = NULL;
6059 gsl_vector *splits = NULL;
6061 m = matrix_expr_evaluate (msave->expr);
6065 if (!common->variables.n)
6066 for (size_t i = 0; i < m->size2; i++)
6067 string_array_append_nocopy (&common->variables,
6068 xasprintf ("COL%zu", i + 1));
6070 if (m->size2 != common->variables.n)
6072 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
6073 m->size2, common->variables.n);
6079 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
6083 if (!common->fnames.n)
6084 for (size_t i = 0; i < factors->size; i++)
6085 string_array_append_nocopy (&common->fnames,
6086 xasprintf ("FAC%zu", i + 1));
6088 if (factors->size != common->fnames.n)
6090 msg (SE, _("There are %zu factor variables, "
6091 "but %zu split values were supplied."),
6092 common->fnames.n, factors->size);
6099 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
6103 if (!common->fnames.n)
6104 for (size_t i = 0; i < splits->size; i++)
6105 string_array_append_nocopy (&common->fnames,
6106 xasprintf ("SPL%zu", i + 1));
6108 if (splits->size != common->fnames.n)
6110 msg (SE, _("There are %zu split variables, "
6111 "but %zu split values were supplied."),
6112 common->fnames.n, splits->size);
6117 if (!common->writer)
6119 struct dictionary *dict = msave_create_dict (common);
6123 common->writer = any_writer_open (common->outfile, dict);
6124 if (!common->writer)
6130 common->dict = dict;
6133 for (size_t y = 0; y < m->size1; y++)
6135 struct ccase *c = case_create (dict_get_proto (common->dict));
6138 /* Split variables */
6140 for (size_t i = 0; i < splits->size; i++)
6141 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
6144 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6145 msave->rowtype, ' ');
6149 for (size_t i = 0; i < factors->size; i++)
6150 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
6153 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6154 msave->varname_, ' ');
6156 /* Continuous variables. */
6157 for (size_t x = 0; x < m->size2; x++)
6158 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
6159 casewriter_write (common->writer, c);
6163 gsl_matrix_free (m);
6164 gsl_vector_free (factors);
6165 gsl_vector_free (splits);
6168 static struct matrix_cmd *
6169 matrix_parse_mget (struct matrix_state *s)
6171 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6172 *cmd = (struct matrix_cmd) {
6176 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
6180 struct mget_command *mget = &cmd->mget;
6182 while (lex_token (s->lexer) != T_ENDCMD)
6184 if (lex_match_id (s->lexer, "FILE"))
6186 lex_match (s->lexer, T_EQUALS);
6188 fh_unref (mget->file);
6189 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6193 else if (lex_match_id (s->lexer, "ENCODING"))
6195 lex_match (s->lexer, T_EQUALS);
6196 if (!lex_force_string (s->lexer))
6199 free (mget->encoding);
6200 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
6204 else if (lex_match_id (s->lexer, "TYPE"))
6206 lex_match (s->lexer, T_EQUALS);
6207 while (lex_token (s->lexer) != T_SLASH
6208 && lex_token (s->lexer) != T_ENDCMD)
6210 const char *rowtype = match_rowtype (s->lexer);
6214 stringi_set_insert (&mget->rowtypes, rowtype);
6219 lex_error_expecting (s->lexer, "FILE", "TYPE");
6222 lex_match (s->lexer, T_SLASH);
6227 matrix_cmd_destroy (cmd);
6231 static const struct variable *
6232 get_a8_var (const struct dictionary *d, const char *name)
6234 const struct variable *v = dict_lookup_var (d, name);
6237 msg (SE, _("Matrix data file lacks %s variable."), name);
6240 if (var_get_width (v) != 8)
6242 msg (SE, _("%s variable in matrix data file must be "
6243 "8-byte string, but it has width %d."),
6244 name, var_get_width (v));
6251 is_rowtype (const union value *v, const char *rowtype)
6253 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
6254 ss_rtrim (&vs, ss_cstr (" "));
6255 return ss_equals_case (vs, ss_cstr (rowtype));
6259 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6260 const struct dictionary *d,
6261 const struct variable *rowtype_var,
6262 struct matrix_state *s,
6263 size_t ss, size_t sn, size_t si,
6264 size_t fs, size_t fn, size_t fi,
6265 size_t cs, size_t cn,
6266 struct pivot_table *pt,
6267 struct pivot_dimension *var_dimension)
6272 const union value *rowtype_ = case_data (rows[0], rowtype_var);
6273 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
6274 : is_rowtype (rowtype_, "CORR") ? "CR"
6275 : is_rowtype (rowtype_, "MEAN") ? "MN"
6276 : is_rowtype (rowtype_, "STDDEV") ? "SD"
6277 : is_rowtype (rowtype_, "N") ? "NC"
6280 struct string name = DS_EMPTY_INITIALIZER;
6281 ds_put_cstr (&name, name_prefix);
6283 ds_put_format (&name, "F%zu", fi);
6285 ds_put_format (&name, "S%zu", si);
6287 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6289 mv = matrix_var_create (s, ds_ss (&name));
6292 msg (SW, _("Matrix data file contains variable with existing name %s."),
6297 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6298 size_t n_missing = 0;
6299 for (size_t y = 0; y < n_rows; y++)
6301 for (size_t x = 0; x < cn; x++)
6303 struct variable *var = dict_get_var (d, cs + x);
6304 double value = case_num (rows[y], var);
6305 if (var_is_num_missing (var, value, MV_ANY))
6310 gsl_matrix_set (m, y, x, value);
6314 int var_index = pivot_category_create_leaf (
6315 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
6316 double values[] = { n_rows, cn };
6317 for (size_t j = 0; j < sn; j++)
6319 struct variable *var = dict_get_var (d, ss + j);
6320 const union value *value = case_data (rows[0], var);
6321 pivot_table_put2 (pt, j, var_index,
6322 pivot_value_new_var_value (var, value));
6324 for (size_t j = 0; j < fn; j++)
6326 struct variable *var = dict_get_var (d, fs + j);
6327 const union value *value = case_data (rows[0], var);
6328 pivot_table_put2 (pt, j + sn, var_index,
6329 pivot_value_new_var_value (var, value));
6331 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
6332 pivot_table_put2 (pt, j + sn + fn, var_index,
6333 pivot_value_new_integer (values[j]));
6336 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6337 "value, which was treated as zero.",
6338 "Matrix data file variable %s contains %zu missing "
6339 "values, which were treated as zero.", n_missing),
6340 ds_cstr (&name), n_missing);
6345 for (size_t y = 0; y < n_rows; y++)
6346 case_unref (rows[y]);
6350 var_changed (const struct ccase *ca, const struct ccase *cb,
6351 const struct variable *var)
6354 ? !value_equal (case_data (ca, var), case_data (cb, var),
6355 var_get_width (var))
6360 vars_changed (const struct ccase *ca, const struct ccase *cb,
6361 const struct dictionary *d,
6362 size_t first_var, size_t n_vars)
6364 for (size_t i = 0; i < n_vars; i++)
6366 const struct variable *v = dict_get_var (d, first_var + i);
6367 if (var_changed (ca, cb, v))
6374 matrix_cmd_execute_mget__ (struct mget_command *mget,
6375 struct casereader *r, const struct dictionary *d)
6377 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6378 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6379 if (!rowtype_ || !varname_)
6382 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6384 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6387 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6389 msg (SE, _("Matrix data file contains no continuous variables."));
6393 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6395 const struct variable *v = dict_get_var (d, i);
6396 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6399 _("Matrix data file contains unexpected string variable %s."),
6405 /* SPLIT variables. */
6407 size_t sn = var_get_dict_index (rowtype_);
6408 struct ccase *sc = NULL;
6411 /* FACTOR variables. */
6412 size_t fs = var_get_dict_index (rowtype_) + 1;
6413 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6414 struct ccase *fc = NULL;
6417 /* Continuous variables. */
6418 size_t cs = var_get_dict_index (varname_) + 1;
6419 size_t cn = dict_get_var_cnt (d) - cs;
6420 struct ccase *cc = NULL;
6423 struct pivot_table *pt = pivot_table_create (
6424 N_("Matrix Variables Created by MGET"));
6425 struct pivot_dimension *attr_dimension = pivot_dimension_create (
6426 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
6427 struct pivot_dimension *var_dimension = pivot_dimension_create (
6428 pt, PIVOT_AXIS_ROW, N_("Variable"));
6431 struct pivot_category *splits = pivot_category_create_group (
6432 attr_dimension->root, N_("Split Values"));
6433 for (size_t i = 0; i < sn; i++)
6434 pivot_category_create_leaf (splits, pivot_value_new_variable (
6435 dict_get_var (d, ss + i)));
6439 struct pivot_category *factors = pivot_category_create_group (
6440 attr_dimension->root, N_("Factor Values"));
6441 for (size_t i = 0; i < fn; i++)
6442 pivot_category_create_leaf (factors, pivot_value_new_variable (
6443 dict_get_var (d, fs + i)));
6445 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
6446 N_("Rows"), N_("Columns"));
6449 struct ccase **rows = NULL;
6450 size_t allocated_rows = 0;
6454 while ((c = casereader_read (r)) != NULL)
6464 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
6465 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
6466 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
6469 if (change != NOTHING_CHANGED)
6471 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6482 if (n_rows >= allocated_rows)
6483 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6486 if (change == SPLITS_CHANGED)
6492 /* Reset the factor number, if there are factors. */
6500 else if (change == FACTORS_CHANGED)
6507 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6519 if (var_dimension->n_leaves)
6520 pivot_table_submit (pt);
6522 pivot_table_unref (pt);
6526 matrix_cmd_execute_mget (struct mget_command *mget)
6528 struct casereader *r;
6529 struct dictionary *d;
6530 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
6531 mget->state->dataset, &r, &d))
6533 matrix_cmd_execute_mget__ (mget, r, d);
6534 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
6539 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6541 if (!lex_force_id (s->lexer))
6544 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6546 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6551 static struct matrix_cmd *
6552 matrix_parse_eigen (struct matrix_state *s)
6554 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6555 *cmd = (struct matrix_cmd) {
6557 .eigen = { .expr = NULL }
6560 struct eigen_command *eigen = &cmd->eigen;
6561 if (!lex_force_match (s->lexer, T_LPAREN))
6563 eigen->expr = matrix_parse_expr (s);
6565 || !lex_force_match (s->lexer, T_COMMA)
6566 || !matrix_parse_dst_var (s, &eigen->evec)
6567 || !lex_force_match (s->lexer, T_COMMA)
6568 || !matrix_parse_dst_var (s, &eigen->eval)
6569 || !lex_force_match (s->lexer, T_RPAREN))
6575 matrix_cmd_destroy (cmd);
6580 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6582 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6585 if (!is_symmetric (A))
6587 msg (SE, _("Argument of EIGEN must be symmetric."));
6588 gsl_matrix_free (A);
6592 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6593 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6594 gsl_vector v_eval = to_vector (eval);
6595 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6596 gsl_eigen_symmv (A, &v_eval, evec, w);
6597 gsl_eigen_symmv_free (w);
6599 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6601 gsl_matrix_free (A);
6603 gsl_matrix_free (eigen->eval->value);
6604 eigen->eval->value = eval;
6606 gsl_matrix_free (eigen->evec->value);
6607 eigen->evec->value = evec;
6610 static struct matrix_cmd *
6611 matrix_parse_setdiag (struct matrix_state *s)
6613 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6614 *cmd = (struct matrix_cmd) {
6615 .type = MCMD_SETDIAG,
6616 .setdiag = { .dst = NULL }
6619 struct setdiag_command *setdiag = &cmd->setdiag;
6620 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6623 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6626 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6631 if (!lex_force_match (s->lexer, T_COMMA))
6634 setdiag->expr = matrix_parse_expr (s);
6638 if (!lex_force_match (s->lexer, T_RPAREN))
6644 matrix_cmd_destroy (cmd);
6649 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6651 gsl_matrix *dst = setdiag->dst->value;
6654 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6655 setdiag->dst->name);
6659 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6663 size_t n = MIN (dst->size1, dst->size2);
6664 if (is_scalar (src))
6666 double d = to_scalar (src);
6667 for (size_t i = 0; i < n; i++)
6668 gsl_matrix_set (dst, i, i, d);
6670 else if (is_vector (src))
6672 gsl_vector v = to_vector (src);
6673 for (size_t i = 0; i < n && i < v.size; i++)
6674 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6677 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6678 "not a %zu×%zu matrix."),
6679 src->size1, src->size2);
6680 gsl_matrix_free (src);
6683 static struct matrix_cmd *
6684 matrix_parse_svd (struct matrix_state *s)
6686 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6687 *cmd = (struct matrix_cmd) {
6689 .svd = { .expr = NULL }
6692 struct svd_command *svd = &cmd->svd;
6693 if (!lex_force_match (s->lexer, T_LPAREN))
6695 svd->expr = matrix_parse_expr (s);
6697 || !lex_force_match (s->lexer, T_COMMA)
6698 || !matrix_parse_dst_var (s, &svd->u)
6699 || !lex_force_match (s->lexer, T_COMMA)
6700 || !matrix_parse_dst_var (s, &svd->s)
6701 || !lex_force_match (s->lexer, T_COMMA)
6702 || !matrix_parse_dst_var (s, &svd->v)
6703 || !lex_force_match (s->lexer, T_RPAREN))
6709 matrix_cmd_destroy (cmd);
6714 matrix_cmd_execute_svd (struct svd_command *svd)
6716 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6720 if (m->size1 >= m->size2)
6723 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6724 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6725 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6726 gsl_vector *work = gsl_vector_alloc (A->size2);
6727 gsl_linalg_SV_decomp (A, V, &Sv, work);
6728 gsl_vector_free (work);
6730 matrix_var_set (svd->u, A);
6731 matrix_var_set (svd->s, S);
6732 matrix_var_set (svd->v, V);
6736 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6737 gsl_matrix_transpose_memcpy (At, m);
6738 gsl_matrix_free (m);
6740 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6741 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6742 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6743 gsl_vector *work = gsl_vector_alloc (At->size2);
6744 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6745 gsl_vector_free (work);
6747 matrix_var_set (svd->v, At);
6748 matrix_var_set (svd->s, St);
6749 matrix_var_set (svd->u, Vt);
6754 matrix_cmd_execute (struct matrix_cmd *cmd)
6759 matrix_cmd_execute_compute (&cmd->compute);
6763 matrix_cmd_execute_print (&cmd->print);
6767 return matrix_cmd_execute_do_if (&cmd->do_if);
6770 matrix_cmd_execute_loop (&cmd->loop);
6777 matrix_cmd_execute_display (&cmd->display);
6781 matrix_cmd_execute_release (&cmd->release);
6785 matrix_cmd_execute_save (&cmd->save);
6789 matrix_cmd_execute_read (&cmd->read);
6793 matrix_cmd_execute_write (&cmd->write);
6797 matrix_cmd_execute_get (&cmd->get);
6801 matrix_cmd_execute_msave (&cmd->msave);
6805 matrix_cmd_execute_mget (&cmd->mget);
6809 matrix_cmd_execute_eigen (&cmd->eigen);
6813 matrix_cmd_execute_setdiag (&cmd->setdiag);
6817 matrix_cmd_execute_svd (&cmd->svd);
6825 matrix_cmds_uninit (struct matrix_cmds *cmds)
6827 for (size_t i = 0; i < cmds->n; i++)
6828 matrix_cmd_destroy (cmds->commands[i]);
6829 free (cmds->commands);
6833 matrix_cmd_destroy (struct matrix_cmd *cmd)
6841 matrix_lvalue_destroy (cmd->compute.lvalue);
6842 matrix_expr_destroy (cmd->compute.rvalue);
6846 matrix_expr_destroy (cmd->print.expression);
6847 free (cmd->print.title);
6848 print_labels_destroy (cmd->print.rlabels);
6849 print_labels_destroy (cmd->print.clabels);
6853 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6855 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6856 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6858 free (cmd->do_if.clauses);
6862 matrix_expr_destroy (cmd->loop.start);
6863 matrix_expr_destroy (cmd->loop.end);
6864 matrix_expr_destroy (cmd->loop.increment);
6865 matrix_expr_destroy (cmd->loop.top_condition);
6866 matrix_expr_destroy (cmd->loop.bottom_condition);
6867 matrix_cmds_uninit (&cmd->loop.commands);
6877 free (cmd->release.vars);
6881 matrix_expr_destroy (cmd->save.expression);
6885 matrix_lvalue_destroy (cmd->read.dst);
6886 matrix_expr_destroy (cmd->read.size);
6890 matrix_expr_destroy (cmd->write.expression);
6891 free (cmd->write.format);
6895 matrix_lvalue_destroy (cmd->get.dst);
6896 fh_unref (cmd->get.file);
6897 free (cmd->get.encoding);
6898 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
6902 free (cmd->msave.varname_);
6903 matrix_expr_destroy (cmd->msave.expr);
6904 matrix_expr_destroy (cmd->msave.factors);
6905 matrix_expr_destroy (cmd->msave.splits);
6909 fh_unref (cmd->mget.file);
6910 stringi_set_destroy (&cmd->mget.rowtypes);
6914 matrix_expr_destroy (cmd->eigen.expr);
6918 matrix_expr_destroy (cmd->setdiag.expr);
6922 matrix_expr_destroy (cmd->svd.expr);
6928 struct matrix_command_name
6931 struct matrix_cmd *(*parse) (struct matrix_state *);
6934 static const struct matrix_command_name *
6935 matrix_parse_command_name (struct lexer *lexer)
6937 static const struct matrix_command_name commands[] = {
6938 { "COMPUTE", matrix_parse_compute },
6939 { "CALL EIGEN", matrix_parse_eigen },
6940 { "CALL SETDIAG", matrix_parse_setdiag },
6941 { "CALL SVD", matrix_parse_svd },
6942 { "PRINT", matrix_parse_print },
6943 { "DO IF", matrix_parse_do_if },
6944 { "LOOP", matrix_parse_loop },
6945 { "BREAK", matrix_parse_break },
6946 { "READ", matrix_parse_read },
6947 { "WRITE", matrix_parse_write },
6948 { "GET", matrix_parse_get },
6949 { "SAVE", matrix_parse_save },
6950 { "MGET", matrix_parse_mget },
6951 { "MSAVE", matrix_parse_msave },
6952 { "DISPLAY", matrix_parse_display },
6953 { "RELEASE", matrix_parse_release },
6955 static size_t n = sizeof commands / sizeof *commands;
6957 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6958 if (lex_match_phrase (lexer, c->name))
6963 static struct matrix_cmd *
6964 matrix_parse_command (struct matrix_state *s)
6966 size_t nesting_level = SIZE_MAX;
6968 struct matrix_cmd *c = NULL;
6969 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6971 lex_error (s->lexer, _("Unknown matrix command."));
6972 else if (!cmd->parse)
6973 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6977 nesting_level = output_open_group (
6978 group_item_create_nocopy (utf8_to_title (cmd->name),
6979 utf8_to_title (cmd->name)));
6984 lex_end_of_command (s->lexer);
6985 lex_discard_rest_of_command (s->lexer);
6986 if (nesting_level != SIZE_MAX)
6987 output_close_groups (nesting_level);
6993 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6995 if (!lex_force_match (lexer, T_ENDCMD))
6998 struct matrix_state state = {
7000 .session = dataset_session (ds),
7002 .vars = HMAP_INITIALIZER (state.vars),
7007 while (lex_match (lexer, T_ENDCMD))
7009 if (lex_token (lexer) == T_STOP)
7011 msg (SE, _("Unexpected end of input expecting matrix command."));
7015 if (lex_match_phrase (lexer, "END MATRIX"))
7018 struct matrix_cmd *c = matrix_parse_command (&state);
7021 matrix_cmd_execute (c);
7022 matrix_cmd_destroy (c);
7026 struct matrix_var *var, *next;
7027 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
7030 gsl_matrix_free (var->value);
7031 hmap_delete (&state.vars, &var->hmap_node);
7034 hmap_destroy (&state.vars);
7037 dict_unref (state.common->dict);
7038 casewriter_destroy (state.common->writer);
7039 free (state.common);
7041 fh_unref (state.prev_read_file);
7042 for (size_t i = 0; i < state.n_read_files; i++)
7043 read_file_destroy (state.read_files[i]);
7044 free (state.read_files);
7045 fh_unref (state.prev_write_file);
7046 for (size_t i = 0; i < state.n_write_files; i++)
7047 write_file_destroy (state.write_files[i]);
7048 free (state.write_files);
7049 fh_unref (state.prev_save_file);
7050 for (size_t i = 0; i < state.n_save_files; i++)
7051 save_file_destroy (state.save_files[i]);
7052 free (state.save_files);