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 dataset *dataset;
111 struct session *session;
115 struct file_handle *prev_save_outfile;
116 struct msave_common *common;
118 struct file_handle *prev_read_file;
119 struct read_file **read_files;
122 struct file_handle *prev_write_file;
123 struct write_file **write_files;
124 size_t n_write_files;
127 static struct matrix_var *
128 matrix_var_lookup (struct matrix_state *s, struct substring name)
130 struct matrix_var *var;
132 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
133 utf8_hash_case_substring (name, 0), &s->vars)
134 if (!utf8_sscasecmp (ss_cstr (var->name), name))
140 static struct matrix_var *
141 matrix_var_create (struct matrix_state *s, struct substring name)
143 struct matrix_var *var = xmalloc (sizeof *var);
144 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
145 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
150 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
152 gsl_matrix_free (var->value);
156 #define MATRIX_FUNCTIONS \
212 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
213 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
214 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
215 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
216 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
217 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
218 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
219 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
220 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
221 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
222 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
224 typedef gsl_matrix *matrix_proto_m_d (double);
225 typedef gsl_matrix *matrix_proto_m_dd (double, double);
226 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
227 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
228 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
229 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
230 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
231 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
232 typedef double matrix_proto_d_m (gsl_matrix *);
233 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
234 typedef gsl_matrix *matrix_proto_IDENT (double, double);
236 #define F(NAME, PROTOTYPE) static matrix_proto_##PROTOTYPE matrix_eval_##NAME;
240 /* Matrix expressions. */
247 #define F(NAME, PROTOTYPE) MOP_F_##NAME,
251 /* Elementwise and scalar arithmetic. */
252 MOP_NEGATE, /* unary - */
253 MOP_ADD_ELEMS, /* + */
254 MOP_SUB_ELEMS, /* - */
255 MOP_MUL_ELEMS, /* &* */
256 MOP_DIV_ELEMS, /* / and &/ */
257 MOP_EXP_ELEMS, /* &** */
259 MOP_SEQ_BY, /* a:b:c */
261 /* Matrix arithmetic. */
263 MOP_EXP_MAT, /* ** */
280 MOP_PASTE_HORZ, /* a, b, c, ... */
281 MOP_PASTE_VERT, /* a; b; c; ... */
285 MOP_VEC_INDEX, /* x(y) */
286 MOP_VEC_ALL, /* x(:) */
287 MOP_MAT_INDEX, /* x(y,z) */
288 MOP_ROW_INDEX, /* x(y,:) */
289 MOP_COL_INDEX, /* x(:,z) */
296 MOP_EOF, /* EOF('file') */
304 struct matrix_expr **subs;
309 struct matrix_var *variable;
310 struct read_file *eof;
315 matrix_expr_destroy (struct matrix_expr *e)
322 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
353 for (size_t i = 0; i < e->n_subs; i++)
354 matrix_expr_destroy (e->subs[i]);
366 static struct matrix_expr *
367 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
370 struct matrix_expr *e = xmalloc (sizeof *e);
371 *e = (struct matrix_expr) {
373 .subs = xmemdup (subs, n_subs * sizeof *subs),
379 static struct matrix_expr *
380 matrix_expr_create_0 (enum matrix_op op)
382 struct matrix_expr *sub;
383 return matrix_expr_create_subs (op, &sub, 0);
386 static struct matrix_expr *
387 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
389 return matrix_expr_create_subs (op, &sub, 1);
392 static struct matrix_expr *
393 matrix_expr_create_2 (enum matrix_op op,
394 struct matrix_expr *sub0, struct matrix_expr *sub1)
396 struct matrix_expr *subs[] = { sub0, sub1 };
397 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
400 static struct matrix_expr *
401 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
402 struct matrix_expr *sub1, struct matrix_expr *sub2)
404 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
405 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
408 static struct matrix_expr *
409 matrix_expr_create_number (double number)
411 struct matrix_expr *e = xmalloc (sizeof *e);
412 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
416 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
418 static struct matrix_expr *
419 matrix_parse_curly_comma (struct matrix_state *s)
421 struct matrix_expr *lhs = matrix_parse_or_xor (s);
425 while (lex_match (s->lexer, T_COMMA))
427 struct matrix_expr *rhs = matrix_parse_or_xor (s);
430 matrix_expr_destroy (lhs);
433 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
438 static struct matrix_expr *
439 matrix_parse_curly_semi (struct matrix_state *s)
441 if (lex_token (s->lexer) == T_RCURLY)
442 return matrix_expr_create_0 (MOP_EMPTY);
444 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
448 while (lex_match (s->lexer, T_SEMICOLON))
450 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
453 matrix_expr_destroy (lhs);
456 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
461 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
462 for (size_t Y = 0; Y < (M)->size1; Y++) \
463 for (size_t X = 0; X < (M)->size2; X++) \
464 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
467 is_vector (const gsl_matrix *m)
469 return m->size1 <= 1 || m->size2 <= 1;
473 to_vector (gsl_matrix *m)
475 return (m->size1 == 1
476 ? gsl_matrix_row (m, 0).vector
477 : gsl_matrix_column (m, 0).vector);
482 matrix_eval_ABS (gsl_matrix *m)
484 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
490 matrix_eval_ALL (gsl_matrix *m)
492 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
499 matrix_eval_ANY (gsl_matrix *m)
501 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
508 matrix_eval_ARSIN (gsl_matrix *m)
510 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
516 matrix_eval_ARTAN (gsl_matrix *m)
518 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
524 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
528 for (size_t i = 0; i < n; i++)
533 gsl_matrix *block = gsl_matrix_calloc (r, c);
535 for (size_t i = 0; i < n; i++)
537 for (size_t y = 0; y < m[i]->size1; y++)
538 for (size_t x = 0; x < m[i]->size2; x++)
539 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
547 matrix_eval_CHOL (gsl_matrix *m)
549 if (!gsl_linalg_cholesky_decomp1 (m))
551 for (size_t y = 0; y < m->size1; y++)
552 for (size_t x = y + 1; x < m->size2; x++)
553 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
555 for (size_t y = 0; y < m->size1; y++)
556 for (size_t x = 0; x < y; x++)
557 gsl_matrix_set (m, y, x, 0);
562 msg (SE, _("Input to CHOL function is not positive-definite."));
568 matrix_eval_col_extremum (gsl_matrix *m, bool min)
573 return gsl_matrix_alloc (1, 0);
575 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
576 for (size_t x = 0; x < m->size2; x++)
578 double ext = gsl_matrix_get (m, 0, x);
579 for (size_t y = 1; y < m->size1; y++)
581 double value = gsl_matrix_get (m, y, x);
582 if (min ? value < ext : value > ext)
585 gsl_matrix_set (cext, 0, x, ext);
591 matrix_eval_CMAX (gsl_matrix *m)
593 return matrix_eval_col_extremum (m, false);
597 matrix_eval_CMIN (gsl_matrix *m)
599 return matrix_eval_col_extremum (m, true);
603 matrix_eval_COS (gsl_matrix *m)
605 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
611 matrix_eval_col_sum (gsl_matrix *m, bool square)
616 return gsl_matrix_alloc (1, 0);
618 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
619 for (size_t x = 0; x < m->size2; x++)
622 for (size_t y = 0; y < m->size1; y++)
624 double d = gsl_matrix_get (m, y, x);
625 sum += square ? pow2 (d) : d;
627 gsl_matrix_set (result, 0, x, sum);
633 matrix_eval_CSSQ (gsl_matrix *m)
635 return matrix_eval_col_sum (m, true);
639 matrix_eval_CSUM (gsl_matrix *m)
641 return matrix_eval_col_sum (m, false);
645 compare_double_3way (const void *a_, const void *b_)
647 const double *a = a_;
648 const double *b = b_;
649 return *a < *b ? -1 : *a > *b;
653 matrix_eval_DESIGN (gsl_matrix *m)
655 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
656 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
657 gsl_matrix_transpose_memcpy (&m2, m);
659 for (size_t y = 0; y < m2.size1; y++)
660 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
662 size_t *n = xcalloc (m2.size1, sizeof *n);
664 for (size_t i = 0; i < m2.size1; i++)
666 double *row = tmp + m2.size2 * i;
667 for (size_t j = 0; j < m2.size2; )
670 for (k = j + 1; k < m2.size2; k++)
671 if (row[j] != row[k])
673 row[n[i]++] = row[j];
678 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
683 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
685 for (size_t i = 0; i < m->size2; i++)
690 const double *unique = tmp + m2.size2 * i;
691 for (size_t j = 0; j < n[i]; j++, x++)
693 double value = unique[j];
694 for (size_t y = 0; y < m->size1; y++)
695 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
706 matrix_eval_DET (gsl_matrix *m)
708 gsl_permutation *p = gsl_permutation_alloc (m->size1);
710 gsl_linalg_LU_decomp (m, p, &signum);
711 gsl_permutation_free (p);
712 return gsl_linalg_LU_det (m, signum);
716 matrix_eval_DIAG (gsl_matrix *m)
718 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
719 for (size_t i = 0; i < diag->size1; i++)
720 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
725 is_symmetric (const gsl_matrix *m)
727 if (m->size1 != m->size2)
730 for (size_t y = 0; y < m->size1; y++)
731 for (size_t x = 0; x < y; x++)
732 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
739 compare_double_desc (const void *a_, const void *b_)
741 const double *a = a_;
742 const double *b = b_;
743 return *a > *b ? -1 : *a < *b;
747 matrix_eval_EVAL (gsl_matrix *m)
749 if (!is_symmetric (m))
751 msg (SE, _("Argument of EVAL must be symmetric."));
755 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
756 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
757 gsl_vector v_eval = to_vector (eval);
758 gsl_eigen_symm (m, &v_eval, w);
759 gsl_eigen_symm_free (w);
761 assert (v_eval.stride == 1);
762 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
768 matrix_eval_EXP (gsl_matrix *m)
770 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
775 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
778 Charl Linssen <charl@itfromb.it>
782 matrix_eval_GINV (gsl_matrix *A)
787 gsl_matrix *tmp_mat = NULL;
790 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
791 tmp_mat = gsl_matrix_alloc (m, n);
792 gsl_matrix_transpose_memcpy (tmp_mat, A);
800 gsl_matrix *V = gsl_matrix_alloc (m, m);
801 gsl_vector *u = gsl_vector_alloc (m);
803 gsl_vector *tmp_vec = gsl_vector_alloc (m);
804 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
805 gsl_vector_free (tmp_vec);
808 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
809 gsl_matrix_set_zero (Sigma_pinv);
810 double cutoff = 1e-15 * gsl_vector_max (u);
812 for (size_t i = 0; i < m; ++i)
814 double x = gsl_vector_get (u, i);
815 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
818 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
819 gsl_matrix *U = gsl_matrix_calloc (n, n);
820 for (size_t i = 0; i < n; i++)
821 for (size_t j = 0; j < m; j++)
822 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
824 /* two dot products to obtain pseudoinverse */
825 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
826 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
831 A_pinv = gsl_matrix_alloc (n, m);
832 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
836 A_pinv = gsl_matrix_alloc (m, n);
837 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
840 gsl_matrix_free (tmp_mat);
841 gsl_matrix_free (tmp_mat2);
843 gsl_matrix_free (Sigma_pinv);
857 grade_compare_3way (const void *a_, const void *b_)
859 const struct grade *a = a_;
860 const struct grade *b = b_;
862 return (a->value < b->value ? -1
863 : a->value > b->value ? 1
871 matrix_eval_GRADE (gsl_matrix *m)
873 size_t n = m->size1 * m->size2;
874 struct grade *grades = xmalloc (n * sizeof *grades);
877 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
878 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
879 qsort (grades, n, sizeof *grades, grade_compare_3way);
881 for (size_t i = 0; i < n; i++)
882 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
890 dot (gsl_vector *a, gsl_vector *b)
893 for (size_t i = 0; i < a->size; i++)
894 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
899 norm2 (gsl_vector *v)
902 for (size_t i = 0; i < v->size; i++)
903 result += pow2 (gsl_vector_get (v, i));
910 return sqrt (norm2 (v));
914 matrix_eval_GSCH (gsl_matrix *v)
916 if (v->size2 < v->size1)
918 msg (SE, _("GSCH requires its argument to have at least as many columns "
919 "as rows, but it has dimensions %zu×%zu."),
923 if (!v->size1 || !v->size2)
926 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
928 for (size_t vx = 0; vx < v->size2; vx++)
930 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
931 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
933 gsl_vector_memcpy (&u_i, &v_i);
934 for (size_t j = 0; j < ux; j++)
936 gsl_vector u_j = gsl_matrix_column (u, j).vector;
937 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
938 for (size_t k = 0; k < u_i.size; k++)
939 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
940 - scale * gsl_vector_get (&u_j, k)));
943 double len = norm (&u_i);
946 gsl_vector_scale (&u_i, 1.0 / len);
947 if (++ux >= v->size1)
954 msg (SE, _("%zu×%zu argument to GSCH contains only "
955 "%zu linearly independent columns."),
956 v->size1, v->size2, ux);
966 matrix_eval_IDENT (double s1, double s2)
968 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
970 msg (SE, _("Arguments to IDENT must be non-negative integers."));
974 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
975 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
981 invert_matrix (gsl_matrix *x)
983 gsl_permutation *p = gsl_permutation_alloc (x->size1);
985 gsl_linalg_LU_decomp (x, p, &signum);
986 gsl_linalg_LU_invx (x, p);
987 gsl_permutation_free (p);
991 matrix_eval_INV (gsl_matrix *m)
998 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1000 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1001 a->size2 * b->size2);
1003 for (size_t ar = 0; ar < a->size1; ar++)
1004 for (size_t br = 0; br < b->size1; br++, y++)
1007 for (size_t ac = 0; ac < a->size2; ac++)
1008 for (size_t bc = 0; bc < b->size2; bc++, x++)
1010 double av = gsl_matrix_get (a, ar, ac);
1011 double bv = gsl_matrix_get (b, br, bc);
1012 gsl_matrix_set (k, y, x, av * bv);
1019 matrix_eval_LG10 (gsl_matrix *m)
1021 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1027 matrix_eval_LN (gsl_matrix *m)
1029 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1035 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1037 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1040 for (size_t i = 1; i <= n * n; i++)
1042 gsl_matrix_set (m, y, x, i);
1044 size_t y1 = !y ? n - 1 : y - 1;
1045 size_t x1 = x + 1 >= n ? 0 : x + 1;
1046 if (gsl_matrix_get (m, y1, x1) == 0)
1052 y = y + 1 >= n ? 0 : y + 1;
1057 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1059 double a = gsl_matrix_get (m, y1, x1);
1060 double b = gsl_matrix_get (m, y2, x2);
1061 gsl_matrix_set (m, y1, x1, b);
1062 gsl_matrix_set (m, y2, x2, a);
1066 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1070 /* A. Umar, "On the Construction of Even Order Magic Squares",
1071 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1073 for (size_t i = 1; i <= n * n / 2; i++)
1075 gsl_matrix_set (m, y, x, i);
1085 for (size_t i = n * n; i > n * n / 2; i--)
1087 gsl_matrix_set (m, y, x, i);
1095 for (size_t y = 0; y < n; y++)
1096 for (size_t x = 0; x < n / 2; x++)
1098 unsigned int d = gsl_matrix_get (m, y, x);
1099 if (d % 2 != (y < n / 2))
1100 magic_exchange (m, y, x, y, n - x - 1);
1105 size_t x1 = n / 2 - 1;
1107 magic_exchange (m, y1, x1, y2, x1);
1108 magic_exchange (m, y1, x2, y2, x2);
1112 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1114 /* A. Umar, "On the Construction of Even Order Magic Squares",
1115 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1119 for (size_t i = 1; ; i++)
1121 gsl_matrix_set (m, y, x, i);
1122 if (++y == n / 2 - 1)
1134 for (size_t i = n * n; ; i--)
1136 gsl_matrix_set (m, y, x, i);
1137 if (++y == n / 2 - 1)
1146 for (size_t y = 0; y < n; y++)
1147 if (y != n / 2 - 1 && y != n / 2)
1148 for (size_t x = 0; x < n / 2; x++)
1150 unsigned int d = gsl_matrix_get (m, y, x);
1151 if (d % 2 != (y < n / 2))
1152 magic_exchange (m, y, x, y, n - x - 1);
1155 size_t a0 = (n * n - 2 * n) / 2 + 1;
1156 for (size_t i = 0; i < n / 2; i++)
1159 gsl_matrix_set (m, n / 2, i, a);
1160 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1162 for (size_t i = 0; i < n / 2; i++)
1164 size_t a = a0 + i + n / 2;
1165 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1166 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1168 for (size_t x = 1; x < n / 2; x += 2)
1169 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1170 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1171 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1172 size_t x1 = n / 2 - 2;
1173 size_t x2 = n / 2 + 1;
1174 size_t y1 = n / 2 - 2;
1175 size_t y2 = n / 2 + 1;
1176 magic_exchange (m, y1, x1, y2, x1);
1177 magic_exchange (m, y1, x2, y2, x2);
1181 matrix_eval_MAGIC (double n_)
1183 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1185 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1190 gsl_matrix *m = gsl_matrix_calloc (n, n);
1192 matrix_eval_MAGIC_odd (m, n);
1194 matrix_eval_MAGIC_singly_even (m, n);
1196 matrix_eval_MAGIC_doubly_even (m, n);
1201 matrix_eval_MAKE (double r, double c, double value)
1203 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1205 msg (SE, _("Size arguments to MAKE must be integers."));
1209 gsl_matrix *m = gsl_matrix_alloc (r, c);
1210 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1216 matrix_eval_MDIAG (gsl_vector *v)
1218 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1219 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1220 gsl_vector_memcpy (&diagonal, v);
1225 matrix_eval_MMAX (gsl_matrix *m)
1227 return gsl_matrix_max (m);
1231 matrix_eval_MMIN (gsl_matrix *m)
1233 return gsl_matrix_min (m);
1237 matrix_eval_MOD (gsl_matrix *m, double divisor)
1241 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1245 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1246 *d = fmod (*d, divisor);
1251 matrix_eval_MSSQ (gsl_matrix *m)
1254 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1260 matrix_eval_MSUM (gsl_matrix *m)
1263 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1269 matrix_eval_NCOL (gsl_matrix *m)
1275 matrix_eval_NROW (gsl_matrix *m)
1281 matrix_eval_RANK (gsl_matrix *m)
1283 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1284 gsl_linalg_QR_decomp (m, tau);
1285 gsl_vector_free (tau);
1287 return gsl_linalg_QRPT_rank (m, -1);
1291 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1293 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1295 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1300 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1302 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1303 "product of matrix dimensions."));
1307 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1310 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1312 gsl_matrix_set (dst, y1, x1, *d);
1323 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1328 return gsl_matrix_alloc (0, 1);
1330 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1331 for (size_t y = 0; y < m->size1; y++)
1333 double ext = gsl_matrix_get (m, y, 0);
1334 for (size_t x = 1; x < m->size2; x++)
1336 double value = gsl_matrix_get (m, y, x);
1337 if (min ? value < ext : value > ext)
1340 gsl_matrix_set (rext, y, 0, ext);
1346 matrix_eval_RMAX (gsl_matrix *m)
1348 return matrix_eval_row_extremum (m, false);
1352 matrix_eval_RMIN (gsl_matrix *m)
1354 return matrix_eval_row_extremum (m, true);
1358 matrix_eval_RND (gsl_matrix *m)
1360 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1372 rank_compare_3way (const void *a_, const void *b_)
1374 const struct rank *a = a_;
1375 const struct rank *b = b_;
1377 return a->value < b->value ? -1 : a->value > b->value;
1381 matrix_eval_RNKORDER (gsl_matrix *m)
1383 size_t n = m->size1 * m->size2;
1384 struct rank *ranks = xmalloc (n * sizeof *ranks);
1386 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1387 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1388 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1390 for (size_t i = 0; i < n; )
1393 for (j = i + 1; j < n; j++)
1394 if (ranks[i].value != ranks[j].value)
1397 double rank = (i + j + 1.0) / 2.0;
1398 for (size_t k = i; k < j; k++)
1399 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1410 matrix_eval_row_sum (gsl_matrix *m, bool square)
1415 return gsl_matrix_alloc (0, 1);
1417 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1418 for (size_t y = 0; y < m->size1; y++)
1421 for (size_t x = 0; x < m->size2; x++)
1423 double d = gsl_matrix_get (m, y, x);
1424 sum += square ? pow2 (d) : d;
1426 gsl_matrix_set (result, y, 0, sum);
1432 matrix_eval_RSSQ (gsl_matrix *m)
1434 return matrix_eval_row_sum (m, true);
1438 matrix_eval_RSUM (gsl_matrix *m)
1440 return matrix_eval_row_sum (m, false);
1444 matrix_eval_SIN (gsl_matrix *m)
1446 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1452 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1454 if (m1->size1 != m2->size1)
1456 msg (SE, _("SOLVE requires its arguments to have the same number of "
1457 "rows, but the first argument has dimensions %zu×%zu and "
1458 "the second %zu×%zu."),
1459 m1->size1, m1->size2,
1460 m2->size1, m2->size2);
1464 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1465 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1467 gsl_linalg_LU_decomp (m1, p, &signum);
1468 for (size_t i = 0; i < m2->size2; i++)
1470 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1471 gsl_vector xi = gsl_matrix_column (x, i).vector;
1472 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1474 gsl_permutation_free (p);
1479 matrix_eval_SQRT (gsl_matrix *m)
1481 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1485 msg (SE, _("Argument to SQRT must be nonnegative."));
1494 matrix_eval_SSCP (gsl_matrix *m)
1496 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1497 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1502 matrix_eval_SVAL (gsl_matrix *m)
1504 gsl_matrix *tmp_mat = NULL;
1505 if (m->size2 > m->size1)
1507 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1508 gsl_matrix_transpose_memcpy (tmp_mat, m);
1513 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1514 gsl_vector *S = gsl_vector_alloc (m->size2);
1515 gsl_vector *work = gsl_vector_alloc (m->size2);
1516 gsl_linalg_SV_decomp (m, V, S, work);
1518 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1519 for (size_t i = 0; i < m->size2; i++)
1520 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1522 gsl_matrix_free (V);
1523 gsl_vector_free (S);
1524 gsl_vector_free (work);
1525 gsl_matrix_free (tmp_mat);
1531 matrix_eval_SWEEP (gsl_matrix *m, double d)
1533 if (d < 1 || d > SIZE_MAX)
1535 msg (SE, _("Scalar argument to SWEEP must be integer."));
1539 if (k >= MIN (m->size1, m->size2))
1541 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1542 "equal to the smaller of the matrix argument's rows and "
1547 double m_kk = gsl_matrix_get (m, k, k);
1548 if (fabs (m_kk) > 1e-19)
1550 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1551 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1553 double m_ij = gsl_matrix_get (m, i, j);
1554 double m_ik = gsl_matrix_get (m, i, k);
1555 double m_kj = gsl_matrix_get (m, k, j);
1556 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1565 for (size_t i = 0; i < m->size1; i++)
1567 gsl_matrix_set (m, i, k, 0);
1568 gsl_matrix_set (m, k, i, 0);
1575 matrix_eval_TRACE (gsl_matrix *m)
1578 size_t n = MIN (m->size1, m->size2);
1579 for (size_t i = 0; i < n; i++)
1580 sum += gsl_matrix_get (m, i, i);
1585 matrix_eval_T (gsl_matrix *m)
1587 return matrix_eval_TRANSPOS (m);
1591 matrix_eval_TRANSPOS (gsl_matrix *m)
1593 if (m->size1 == m->size2)
1595 gsl_matrix_transpose (m);
1600 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1601 gsl_matrix_transpose_memcpy (t, m);
1607 matrix_eval_TRUNC (gsl_matrix *m)
1609 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1615 matrix_eval_UNIFORM (double r_, double c_)
1617 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1619 msg (SE, _("Arguments to UNIFORM must be integers."));
1624 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1626 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1630 gsl_matrix *m = gsl_matrix_alloc (r, c);
1631 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1632 *d = gsl_ran_flat (get_rng (), 0, 1);
1636 struct matrix_function
1640 size_t min_args, max_args;
1643 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
1645 static const struct matrix_function *
1646 matrix_parse_function_name (const char *token)
1648 static const struct matrix_function functions[] =
1650 #define F(NAME, PROTO) { #NAME, MOP_F_##NAME, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
1654 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
1656 for (size_t i = 0; i < N_FUNCTIONS; i++)
1658 if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3))
1659 return &functions[i];
1664 static struct read_file *
1665 read_file_create (struct matrix_state *s, struct file_handle *fh)
1667 for (size_t i = 0; i < s->n_read_files; i++)
1669 struct read_file *rf = s->read_files[i];
1677 struct read_file *rf = xmalloc (sizeof *rf);
1678 *rf = (struct read_file) { .file = fh };
1680 s->read_files = xrealloc (s->read_files,
1681 (s->n_read_files + 1) * sizeof *s->read_files);
1682 s->read_files[s->n_read_files++] = rf;
1686 static struct dfm_reader *
1687 read_file_open (struct read_file *rf)
1690 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
1695 read_file_destroy (struct read_file *rf)
1699 fh_unref (rf->file);
1700 dfm_close_reader (rf->reader);
1701 free (rf->encoding);
1706 static struct write_file *
1707 write_file_create (struct matrix_state *s, struct file_handle *fh)
1709 for (size_t i = 0; i < s->n_write_files; i++)
1711 struct write_file *wf = s->write_files[i];
1719 struct write_file *wf = xmalloc (sizeof *wf);
1720 *wf = (struct write_file) { .file = fh };
1722 s->write_files = xrealloc (s->write_files,
1723 (s->n_write_files + 1) * sizeof *s->write_files);
1724 s->write_files[s->n_write_files++] = wf;
1728 static struct dfm_writer *
1729 write_file_open (struct write_file *wf)
1732 wf->writer = dfm_open_writer (wf->file, wf->encoding);
1737 write_file_destroy (struct write_file *wf)
1743 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
1744 wf->held->s.ss.length);
1745 u8_line_destroy (wf->held);
1749 fh_unref (wf->file);
1750 dfm_close_writer (wf->writer);
1751 free (wf->encoding);
1757 matrix_parse_function (struct matrix_state *s, const char *token,
1758 struct matrix_expr **exprp)
1761 if (lex_next_token (s->lexer, 1) != T_LPAREN)
1764 if (lex_match_id (s->lexer, "EOF"))
1767 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
1771 if (!lex_force_match (s->lexer, T_RPAREN))
1777 struct read_file *rf = read_file_create (s, fh);
1779 struct matrix_expr *e = xmalloc (sizeof *e);
1780 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
1785 const struct matrix_function *f = matrix_parse_function_name (token);
1789 lex_get_n (s->lexer, 2);
1791 struct matrix_expr *e = xmalloc (sizeof *e);
1792 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
1794 size_t allocated_subs = 0;
1797 struct matrix_expr *sub = matrix_parse_expr (s);
1801 if (e->n_subs >= allocated_subs)
1802 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
1803 e->subs[e->n_subs++] = sub;
1805 while (lex_match (s->lexer, T_COMMA));
1806 if (!lex_force_match (s->lexer, T_RPAREN))
1813 matrix_expr_destroy (e);
1817 static struct matrix_expr *
1818 matrix_parse_primary (struct matrix_state *s)
1820 if (lex_is_number (s->lexer))
1822 double number = lex_number (s->lexer);
1824 return matrix_expr_create_number (number);
1826 else if (lex_is_string (s->lexer))
1828 char string[sizeof (double)];
1829 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1833 memcpy (&number, string, sizeof number);
1834 return matrix_expr_create_number (number);
1836 else if (lex_match (s->lexer, T_LPAREN))
1838 struct matrix_expr *e = matrix_parse_or_xor (s);
1839 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1841 matrix_expr_destroy (e);
1846 else if (lex_match (s->lexer, T_LCURLY))
1848 struct matrix_expr *e = matrix_parse_curly_semi (s);
1849 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1851 matrix_expr_destroy (e);
1856 else if (lex_token (s->lexer) == T_ID)
1858 struct matrix_expr *retval;
1859 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1862 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1865 lex_error (s->lexer, _("Unknown variable %s."),
1866 lex_tokcstr (s->lexer));
1871 struct matrix_expr *e = xmalloc (sizeof *e);
1872 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1875 else if (lex_token (s->lexer) == T_ALL)
1877 struct matrix_expr *retval;
1878 if (matrix_parse_function (s, "ALL", &retval))
1882 lex_error (s->lexer, NULL);
1886 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1889 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
1891 if (lex_match (s->lexer, T_COLON))
1898 *indexp = matrix_parse_expr (s);
1899 return *indexp != NULL;
1903 static struct matrix_expr *
1904 matrix_parse_postfix (struct matrix_state *s)
1906 struct matrix_expr *lhs = matrix_parse_primary (s);
1907 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1910 struct matrix_expr *i0;
1911 if (!matrix_parse_index_expr (s, &i0))
1913 matrix_expr_destroy (lhs);
1916 if (lex_match (s->lexer, T_RPAREN))
1918 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1919 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1920 else if (lex_match (s->lexer, T_COMMA))
1922 struct matrix_expr *i1;
1923 if (!matrix_parse_index_expr (s, &i1)
1924 || !lex_force_match (s->lexer, T_RPAREN))
1926 matrix_expr_destroy (lhs);
1927 matrix_expr_destroy (i0);
1928 matrix_expr_destroy (i1);
1931 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1932 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1933 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1938 lex_error_expecting (s->lexer, "`)'", "`,'");
1943 static struct matrix_expr *
1944 matrix_parse_unary (struct matrix_state *s)
1946 if (lex_match (s->lexer, T_DASH))
1948 struct matrix_expr *lhs = matrix_parse_unary (s);
1949 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
1951 else if (lex_match (s->lexer, T_PLUS))
1952 return matrix_parse_unary (s);
1954 return matrix_parse_postfix (s);
1957 static struct matrix_expr *
1958 matrix_parse_seq (struct matrix_state *s)
1960 struct matrix_expr *start = matrix_parse_unary (s);
1961 if (!start || !lex_match (s->lexer, T_COLON))
1964 struct matrix_expr *end = matrix_parse_unary (s);
1967 matrix_expr_destroy (start);
1971 if (lex_match (s->lexer, T_COLON))
1973 struct matrix_expr *increment = matrix_parse_unary (s);
1976 matrix_expr_destroy (start);
1977 matrix_expr_destroy (end);
1980 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
1983 return matrix_expr_create_2 (MOP_SEQ, start, end);
1986 static struct matrix_expr *
1987 matrix_parse_exp (struct matrix_state *s)
1989 struct matrix_expr *lhs = matrix_parse_seq (s);
1996 if (lex_match (s->lexer, T_EXP))
1998 else if (lex_match_phrase (s->lexer, "&**"))
2003 struct matrix_expr *rhs = matrix_parse_seq (s);
2006 matrix_expr_destroy (lhs);
2009 lhs = matrix_expr_create_2 (op, lhs, rhs);
2013 static struct matrix_expr *
2014 matrix_parse_mul_div (struct matrix_state *s)
2016 struct matrix_expr *lhs = matrix_parse_exp (s);
2023 if (lex_match (s->lexer, T_ASTERISK))
2025 else if (lex_match (s->lexer, T_SLASH))
2027 else if (lex_match_phrase (s->lexer, "&*"))
2029 else if (lex_match_phrase (s->lexer, "&/"))
2034 struct matrix_expr *rhs = matrix_parse_exp (s);
2037 matrix_expr_destroy (lhs);
2040 lhs = matrix_expr_create_2 (op, lhs, rhs);
2044 static struct matrix_expr *
2045 matrix_parse_add_sub (struct matrix_state *s)
2047 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2054 if (lex_match (s->lexer, T_PLUS))
2056 else if (lex_match (s->lexer, T_DASH))
2058 else if (lex_token (s->lexer) == T_NEG_NUM)
2063 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2066 matrix_expr_destroy (lhs);
2069 lhs = matrix_expr_create_2 (op, lhs, rhs);
2073 static struct matrix_expr *
2074 matrix_parse_relational (struct matrix_state *s)
2076 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2083 if (lex_match (s->lexer, T_GT))
2085 else if (lex_match (s->lexer, T_GE))
2087 else if (lex_match (s->lexer, T_LT))
2089 else if (lex_match (s->lexer, T_LE))
2091 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2093 else if (lex_match (s->lexer, T_NE))
2098 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2101 matrix_expr_destroy (lhs);
2104 lhs = matrix_expr_create_2 (op, lhs, rhs);
2108 static struct matrix_expr *
2109 matrix_parse_not (struct matrix_state *s)
2111 if (lex_match (s->lexer, T_NOT))
2113 struct matrix_expr *lhs = matrix_parse_not (s);
2114 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2117 return matrix_parse_relational (s);
2120 static struct matrix_expr *
2121 matrix_parse_and (struct matrix_state *s)
2123 struct matrix_expr *lhs = matrix_parse_not (s);
2127 while (lex_match (s->lexer, T_AND))
2129 struct matrix_expr *rhs = matrix_parse_not (s);
2132 matrix_expr_destroy (lhs);
2135 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2140 static struct matrix_expr *
2141 matrix_parse_or_xor (struct matrix_state *s)
2143 struct matrix_expr *lhs = matrix_parse_and (s);
2150 if (lex_match (s->lexer, T_OR))
2152 else if (lex_match_id (s->lexer, "XOR"))
2157 struct matrix_expr *rhs = matrix_parse_and (s);
2160 matrix_expr_destroy (lhs);
2163 lhs = matrix_expr_create_2 (op, lhs, rhs);
2167 static struct matrix_expr *
2168 matrix_parse_expr (struct matrix_state *s)
2170 return matrix_parse_or_xor (s);
2173 /* Expression evaluation. */
2176 matrix_op_eval (enum matrix_op op, double a, double b)
2180 case MOP_ADD_ELEMS: return a + b;
2181 case MOP_SUB_ELEMS: return a - b;
2182 case MOP_MUL_ELEMS: return a * b;
2183 case MOP_DIV_ELEMS: return a / b;
2184 case MOP_EXP_ELEMS: return pow (a, b);
2185 case MOP_GT: return a > b;
2186 case MOP_GE: return a >= b;
2187 case MOP_LT: return a < b;
2188 case MOP_LE: return a <= b;
2189 case MOP_EQ: return a == b;
2190 case MOP_NE: return a != b;
2191 case MOP_AND: return (a > 0) && (b > 0);
2192 case MOP_OR: return (a > 0) || (b > 0);
2193 case MOP_XOR: return (a > 0) != (b > 0);
2195 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2204 case MOP_PASTE_HORZ:
2205 case MOP_PASTE_VERT:
2221 matrix_op_name (enum matrix_op op)
2225 case MOP_ADD_ELEMS: return "+";
2226 case MOP_SUB_ELEMS: return "-";
2227 case MOP_MUL_ELEMS: return "&*";
2228 case MOP_DIV_ELEMS: return "&/";
2229 case MOP_EXP_ELEMS: return "&**";
2230 case MOP_GT: return ">";
2231 case MOP_GE: return ">=";
2232 case MOP_LT: return "<";
2233 case MOP_LE: return "<=";
2234 case MOP_EQ: return "=";
2235 case MOP_NE: return "<>";
2236 case MOP_AND: return "AND";
2237 case MOP_OR: return "OR";
2238 case MOP_XOR: return "XOR";
2240 #define F(NAME, PROTOTYPE) case MOP_F_##NAME:
2249 case MOP_PASTE_HORZ:
2250 case MOP_PASTE_VERT:
2266 is_scalar (const gsl_matrix *m)
2268 return m->size1 == 1 && m->size2 == 1;
2272 to_scalar (const gsl_matrix *m)
2274 assert (is_scalar (m));
2275 return gsl_matrix_get (m, 0, 0);
2279 matrix_expr_evaluate_elementwise (enum matrix_op op,
2280 gsl_matrix *a, gsl_matrix *b)
2284 double be = to_scalar (b);
2285 for (size_t r = 0; r < a->size1; r++)
2286 for (size_t c = 0; c < a->size2; c++)
2288 double *ae = gsl_matrix_ptr (a, r, c);
2289 *ae = matrix_op_eval (op, *ae, be);
2293 else if (is_scalar (a))
2295 double ae = to_scalar (a);
2296 for (size_t r = 0; r < b->size1; r++)
2297 for (size_t c = 0; c < b->size2; c++)
2299 double *be = gsl_matrix_ptr (b, r, c);
2300 *be = matrix_op_eval (op, ae, *be);
2304 else if (a->size1 == b->size1 && a->size2 == b->size2)
2306 for (size_t r = 0; r < a->size1; r++)
2307 for (size_t c = 0; c < a->size2; c++)
2309 double *ae = gsl_matrix_ptr (a, r, c);
2310 double be = gsl_matrix_get (b, r, c);
2311 *ae = matrix_op_eval (op, *ae, be);
2317 msg (SE, _("Operands to %s must have the same dimensions or one "
2318 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
2319 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
2325 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
2327 if (is_scalar (a) || is_scalar (b))
2328 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
2330 if (a->size2 != b->size1)
2332 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
2333 "not conformable for multiplication."),
2334 a->size1, a->size2, b->size1, b->size2);
2338 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
2339 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
2344 swap_matrix (gsl_matrix **a, gsl_matrix **b)
2346 gsl_matrix *tmp = *a;
2352 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
2355 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
2356 swap_matrix (z, tmp);
2360 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
2362 mul_matrix (x, *x, *x, tmp);
2366 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
2369 if (x->size1 != x->size2)
2371 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
2372 "the left-hand size, not one with dimensions %zu×%zu."),
2373 x->size1, x->size2);
2378 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
2379 "right-hand side, not a matrix with dimensions %zu×%zu."),
2380 b->size1, b->size2);
2383 double bf = to_scalar (b);
2384 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
2386 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
2387 "or outside the valid range."), bf);
2392 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
2394 gsl_matrix_set_identity (y);
2398 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
2400 for (unsigned long int n = labs (bl); n > 1; n /= 2)
2403 mul_matrix (&y, x, y, &t);
2404 square_matrix (&x, &t);
2407 square_matrix (&x, &t);
2409 mul_matrix (&y, x, y, &t);
2413 /* Garbage collection.
2415 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
2416 a permutation of them. We are returning one of them; that one must not be
2417 destroyed. We must not destroy 'x_' because the caller owns it. */
2419 gsl_matrix_free (y_);
2421 gsl_matrix_free (t_);
2427 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
2430 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
2432 msg (SE, _("All operands of : operator must be scalars."));
2436 long int start = to_scalar (start_);
2437 long int end = to_scalar (end_);
2438 long int by = by_ ? to_scalar (by_) : 1;
2442 msg (SE, _("The increment operand to : must be nonzero."));
2446 long int n = (end >= start && by > 0 ? (end - start + by) / by
2447 : end <= start && by < 0 ? (start - end - by) / -by
2449 gsl_matrix *m = gsl_matrix_alloc (1, n);
2450 for (long int i = 0; i < n; i++)
2451 gsl_matrix_set (m, 0, i, start + i * by);
2456 matrix_expr_evaluate_not (gsl_matrix *a)
2458 for (size_t r = 0; r < a->size1; r++)
2459 for (size_t c = 0; c < a->size2; c++)
2461 double *ae = gsl_matrix_ptr (a, r, c);
2468 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
2470 if (a->size1 != b->size1)
2472 if (!a->size1 || !a->size2)
2474 else if (!b->size1 || !b->size2)
2477 msg (SE, _("All columns in a matrix must have the same number of rows, "
2478 "but this tries to paste matrices with %zu and %zu rows."),
2479 a->size1, b->size1);
2483 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
2484 for (size_t y = 0; y < a->size1; y++)
2486 for (size_t x = 0; x < a->size2; x++)
2487 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2488 for (size_t x = 0; x < b->size2; x++)
2489 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
2495 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
2497 if (a->size2 != b->size2)
2499 if (!a->size1 || !a->size2)
2501 else if (!b->size1 || !b->size2)
2504 msg (SE, _("All rows in a matrix must have the same number of columns, "
2505 "but this tries to stack matrices with %zu and %zu columns."),
2506 a->size2, b->size2);
2510 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
2511 for (size_t x = 0; x < a->size2; x++)
2513 for (size_t y = 0; y < a->size1; y++)
2514 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
2515 for (size_t y = 0; y < b->size1; y++)
2516 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
2522 matrix_to_vector (gsl_matrix *m)
2525 gsl_vector v = to_vector (m);
2526 assert (v.block == m->block || !v.block);
2530 return xmemdup (&v, sizeof v);
2544 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
2547 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
2548 enum index_type index_type, size_t other_size,
2549 struct index_vector *iv)
2558 msg (SE, _("Vector index must be scalar or vector, not a "
2560 m->size1, m->size2);
2564 msg (SE, _("Matrix row index must be scalar or vector, not a "
2566 m->size1, m->size2);
2570 msg (SE, _("Matrix column index must be scalar or vector, not a "
2572 m->size1, m->size2);
2578 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
2579 *iv = (struct index_vector) {
2580 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
2583 for (size_t i = 0; i < v.size; i++)
2585 double index = gsl_vector_get (&v, i);
2586 if (index < 1 || index >= size + 1)
2591 msg (SE, _("Index %g is out of range for vector "
2592 "with %zu elements."), index, size);
2596 msg (SE, _("%g is not a valid row index for "
2597 "a %zu×%zu matrix."),
2598 index, size, other_size);
2602 msg (SE, _("%g is not a valid column index for "
2603 "a %zu×%zu matrix."),
2604 index, other_size, size);
2611 iv->indexes[i] = index - 1;
2617 *iv = (struct index_vector) {
2618 .indexes = xnmalloc (size, sizeof *iv->indexes),
2621 for (size_t i = 0; i < size; i++)
2628 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
2630 if (!is_vector (sm))
2632 msg (SE, _("Vector index operator may not be applied to "
2633 "a %zu×%zu matrix."),
2634 sm->size1, sm->size2);
2642 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
2644 if (!matrix_expr_evaluate_vec_all (sm))
2647 gsl_vector sv = to_vector (sm);
2648 struct index_vector iv;
2649 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
2652 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
2653 sm->size1 == 1 ? iv.n : 1);
2654 gsl_vector dv = to_vector (dm);
2655 for (size_t dx = 0; dx < iv.n; dx++)
2657 size_t sx = iv.indexes[dx];
2658 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
2666 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
2669 struct index_vector iv0;
2670 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
2673 struct index_vector iv1;
2674 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
2681 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
2682 for (size_t dy = 0; dy < iv0.n; dy++)
2684 size_t sy = iv0.indexes[dy];
2686 for (size_t dx = 0; dx < iv1.n; dx++)
2688 size_t sx = iv1.indexes[dx];
2689 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
2697 #define F(NAME, PROTOTYPE) \
2698 static gsl_matrix *matrix_expr_evaluate_##PROTOTYPE ( \
2699 const char *name, gsl_matrix *[], size_t, \
2700 matrix_proto_##PROTOTYPE *);
2705 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
2707 if (!is_scalar (subs[index]))
2709 msg (SE, _("Function %s argument %zu must be a scalar, "
2710 "not a %zu×%zu matrix."),
2711 name, index + 1, subs[index]->size1, subs[index]->size2);
2718 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
2720 if (!is_vector (subs[index]))
2722 msg (SE, _("Function %s argument %zu must be a vector, "
2723 "not a %zu×%zu matrix."),
2724 name, index + 1, subs[index]->size1, subs[index]->size2);
2731 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
2733 for (size_t i = 0; i < n_subs; i++)
2735 if (!check_scalar_arg (name, subs, i))
2737 d[i] = to_scalar (subs[i]);
2743 matrix_expr_evaluate_m_d (const char *name,
2744 gsl_matrix *subs[], size_t n_subs,
2745 matrix_proto_m_d *f)
2747 assert (n_subs == 1);
2750 return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL;
2754 matrix_expr_evaluate_m_dd (const char *name,
2755 gsl_matrix *subs[], size_t n_subs,
2756 matrix_proto_m_dd *f)
2758 assert (n_subs == 2);
2761 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
2765 matrix_expr_evaluate_m_ddd (const char *name,
2766 gsl_matrix *subs[], size_t n_subs,
2767 matrix_proto_m_ddd *f)
2769 assert (n_subs == 3);
2772 return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
2776 matrix_expr_evaluate_m_m (const char *name UNUSED,
2777 gsl_matrix *subs[], size_t n_subs,
2778 matrix_proto_m_m *f)
2780 assert (n_subs == 1);
2785 matrix_expr_evaluate_m_md (const char *name UNUSED,
2786 gsl_matrix *subs[], size_t n_subs,
2787 matrix_proto_m_md *f)
2789 assert (n_subs == 2);
2790 if (!check_scalar_arg (name, subs, 1))
2792 return f (subs[0], to_scalar (subs[1]));
2796 matrix_expr_evaluate_m_mdd (const char *name UNUSED,
2797 gsl_matrix *subs[], size_t n_subs,
2798 matrix_proto_m_mdd *f)
2800 assert (n_subs == 3);
2801 if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2))
2803 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
2807 matrix_expr_evaluate_m_mm (const char *name UNUSED,
2808 gsl_matrix *subs[], size_t n_subs,
2809 matrix_proto_m_mm *f)
2811 assert (n_subs == 2);
2812 return f (subs[0], subs[1]);
2816 matrix_expr_evaluate_m_v (const char *name,
2817 gsl_matrix *subs[], size_t n_subs,
2818 matrix_proto_m_v *f)
2820 assert (n_subs == 1);
2821 if (!check_vector_arg (name, subs, 0))
2823 gsl_vector v = to_vector (subs[0]);
2828 matrix_expr_evaluate_d_m (const char *name UNUSED,
2829 gsl_matrix *subs[], size_t n_subs,
2830 matrix_proto_d_m *f)
2832 assert (n_subs == 1);
2834 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2835 gsl_matrix_set (m, 0, 0, f (subs[0]));
2840 matrix_expr_evaluate_m_any (const char *name UNUSED,
2841 gsl_matrix *subs[], size_t n_subs,
2842 matrix_proto_m_any *f)
2844 return f (subs, n_subs);
2848 matrix_expr_evaluate_IDENT (const char *name,
2849 gsl_matrix *subs[], size_t n_subs,
2850 matrix_proto_IDENT *f)
2852 assert (n_subs <= 2);
2855 if (!to_scalar_args (name, subs, n_subs, d))
2858 return f (d[0], d[n_subs - 1]);
2862 matrix_expr_evaluate (const struct matrix_expr *e)
2864 if (e->op == MOP_NUMBER)
2866 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2867 gsl_matrix_set (m, 0, 0, e->number);
2870 else if (e->op == MOP_VARIABLE)
2872 const gsl_matrix *src = e->variable->value;
2875 msg (SE, _("Uninitialized variable %s used in expression."),
2880 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
2881 gsl_matrix_memcpy (dst, src);
2884 else if (e->op == MOP_EOF)
2886 struct dfm_reader *reader = read_file_open (e->eof);
2887 gsl_matrix *m = gsl_matrix_alloc (1, 1);
2888 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
2892 enum { N_LOCAL = 3 };
2893 gsl_matrix *local_subs[N_LOCAL];
2894 gsl_matrix **subs = (e->n_subs < N_LOCAL
2896 : xmalloc (e->n_subs * sizeof *subs));
2898 for (size_t i = 0; i < e->n_subs; i++)
2900 subs[i] = matrix_expr_evaluate (e->subs[i]);
2903 for (size_t j = 0; j < i; j++)
2904 gsl_matrix_free (subs[j]);
2905 if (subs != local_subs)
2911 gsl_matrix *result = NULL;
2914 #define F(NAME, PROTOTYPE) \
2915 case MOP_F_##NAME: \
2916 result = matrix_expr_evaluate_##PROTOTYPE (#NAME, \
2918 matrix_eval_##NAME); \
2924 gsl_matrix_scale (subs[0], -1.0);
2942 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
2946 result = matrix_expr_evaluate_not (subs[0]);
2950 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
2954 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
2958 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
2962 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
2965 case MOP_PASTE_HORZ:
2966 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
2969 case MOP_PASTE_VERT:
2970 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
2974 result = gsl_matrix_alloc (0, 0);
2978 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
2982 result = matrix_expr_evaluate_vec_all (subs[0]);
2986 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
2990 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
2994 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
3003 for (size_t i = 0; i < e->n_subs; i++)
3004 if (subs[i] != result)
3005 gsl_matrix_free (subs[i]);
3006 if (subs != local_subs)
3012 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
3015 gsl_matrix *m = matrix_expr_evaluate (e);
3021 msg (SE, _("Expression for %s must evaluate to scalar, "
3022 "not a %zu×%zu matrix."),
3023 context, m->size1, m->size2);
3024 gsl_matrix_free (m);
3029 gsl_matrix_free (m);
3034 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
3038 if (!matrix_expr_evaluate_scalar (e, context, &d))
3042 if (d < LONG_MIN || d > LONG_MAX)
3044 msg (SE, _("Expression for %s is outside the integer range."), context);
3051 struct matrix_lvalue
3053 struct matrix_var *var;
3054 struct matrix_expr *indexes[2];
3059 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
3063 for (size_t i = 0; i < lvalue->n_indexes; i++)
3064 matrix_expr_destroy (lvalue->indexes[i]);
3069 static struct matrix_lvalue *
3070 matrix_lvalue_parse (struct matrix_state *s)
3072 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
3073 if (!lex_force_id (s->lexer))
3075 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
3076 if (lex_next_token (s->lexer, 1) == T_LPAREN)
3080 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
3085 lex_get_n (s->lexer, 2);
3087 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
3089 lvalue->n_indexes++;
3091 if (lex_match (s->lexer, T_COMMA))
3093 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
3095 lvalue->n_indexes++;
3097 if (!lex_force_match (s->lexer, T_RPAREN))
3103 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
3109 matrix_lvalue_destroy (lvalue);
3114 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
3115 enum index_type index_type, size_t other_size,
3116 struct index_vector *iv)
3121 m = matrix_expr_evaluate (e);
3128 bool ok = matrix_normalize_index_vector (m, size, index_type,
3130 gsl_matrix_free (m);
3135 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
3136 struct index_vector *iv, gsl_matrix *sm)
3138 gsl_vector dv = to_vector (lvalue->var->value);
3140 /* Convert source matrix 'sm' to source vector 'sv'. */
3141 if (!is_vector (sm))
3143 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
3144 sm->size1, sm->size2);
3147 gsl_vector sv = to_vector (sm);
3148 if (iv->n != sv.size)
3150 msg (SE, _("Can't assign %zu-element vector "
3151 "to %zu-element subvector."), sv.size, iv->n);
3155 /* Assign elements. */
3156 for (size_t x = 0; x < iv->n; x++)
3157 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
3162 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
3163 struct index_vector *iv0,
3164 struct index_vector *iv1, gsl_matrix *sm)
3166 gsl_matrix *dm = lvalue->var->value;
3168 /* Convert source matrix 'sm' to source vector 'sv'. */
3169 if (iv0->n != sm->size1)
3171 msg (SE, _("Row index vector for assignment to %s has %zu elements "
3172 "but source matrix has %zu rows."),
3173 lvalue->var->name, iv0->n, sm->size1);
3176 if (iv1->n != sm->size2)
3178 msg (SE, _("Column index vector for assignment to %s has %zu elements "
3179 "but source matrix has %zu columns."),
3180 lvalue->var->name, iv1->n, sm->size2);
3184 /* Assign elements. */
3185 for (size_t y = 0; y < iv0->n; y++)
3187 size_t f0 = iv0->indexes[y];
3188 for (size_t x = 0; x < iv1->n; x++)
3190 size_t f1 = iv1->indexes[x];
3191 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
3197 /* Takes ownership of M. */
3199 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
3200 struct index_vector *iv0, struct index_vector *iv1,
3203 if (!lvalue->n_indexes)
3205 gsl_matrix_free (lvalue->var->value);
3206 lvalue->var->value = sm;
3211 bool ok = (lvalue->n_indexes == 1
3212 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
3213 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
3214 free (iv0->indexes);
3215 free (iv1->indexes);
3216 gsl_matrix_free (sm);
3222 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
3223 struct index_vector *iv0,
3224 struct index_vector *iv1)
3226 *iv0 = INDEX_VECTOR_INIT;
3227 *iv1 = INDEX_VECTOR_INIT;
3228 if (!lvalue->n_indexes)
3231 /* Validate destination matrix exists and has the right shape. */
3232 gsl_matrix *dm = lvalue->var->value;
3235 msg (SE, _("Undefined variable %s."), lvalue->var->name);
3238 else if (dm->size1 == 0 || dm->size2 == 0)
3240 msg (SE, _("Cannot index %zu×%zu matrix."),
3241 dm->size1, dm->size2);
3244 else if (lvalue->n_indexes == 1)
3246 if (!is_vector (dm))
3248 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
3249 dm->size1, dm->size2, lvalue->var->name);
3252 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
3253 MAX (dm->size1, dm->size2),
3258 assert (lvalue->n_indexes == 2);
3259 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
3260 IV_ROW, dm->size2, iv0))
3263 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
3264 IV_COLUMN, dm->size1, iv1))
3266 free (iv0->indexes);
3273 /* Takes ownership of M. */
3275 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
3277 struct index_vector iv0, iv1;
3278 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
3280 gsl_matrix_free (sm);
3284 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
3288 /* Matrix command. */
3292 struct matrix_cmd **commands;
3296 static void matrix_cmds_uninit (struct matrix_cmds *);
3300 enum matrix_cmd_type
3323 struct compute_command
3325 struct matrix_lvalue *lvalue;
3326 struct matrix_expr *rvalue;
3330 struct print_command
3332 struct matrix_expr *expression;
3333 bool use_default_format;
3334 struct fmt_spec format;
3336 int space; /* -1 means NEWPAGE. */
3340 struct string_array literals; /* CLABELS/RLABELS. */
3341 struct matrix_expr *expr; /* CNAMES/RNAMES. */
3347 struct do_if_command
3351 struct matrix_expr *condition;
3352 struct matrix_cmds commands;
3362 struct matrix_var *var;
3363 struct matrix_expr *start;
3364 struct matrix_expr *end;
3365 struct matrix_expr *increment;
3367 /* Loop conditions. */
3368 struct matrix_expr *top_condition;
3369 struct matrix_expr *bottom_condition;
3372 struct matrix_cmds commands;
3376 struct display_command
3378 struct matrix_state *state;
3384 struct release_command
3386 struct matrix_var **vars;
3393 struct matrix_expr *expression;
3394 struct file_handle *outfile;
3395 struct string_array *variables;
3396 struct matrix_expr *names;
3397 struct stringi_set strings;
3403 struct read_file *rf;
3404 struct matrix_lvalue *dst;
3405 struct matrix_expr *size;
3407 enum fmt_type format;
3415 struct write_command
3417 struct write_file *wf;
3418 struct matrix_expr *expression;
3421 /* If this is nonnull, WRITE uses this format.
3423 If this is NULL, WRITE uses free-field format with as many
3424 digits of precision as needed. */
3425 struct fmt_spec *format;
3434 struct matrix_lvalue *dst;
3435 struct dataset *dataset;
3436 struct file_handle *file;
3438 struct string_array variables;
3439 struct matrix_var *names;
3441 /* Treatment of missing values. */
3446 MGET_ERROR, /* Flag error on command. */
3447 MGET_ACCEPT, /* Accept without change, user-missing only. */
3448 MGET_OMIT, /* Drop this case. */
3449 MGET_RECODE /* Recode to 'substitute'. */
3452 double substitute; /* MGET_RECODE only. */
3458 struct msave_command
3460 struct msave_common *common;
3462 struct matrix_expr *expr;
3463 const char *rowtype;
3464 struct matrix_expr *factors;
3465 struct matrix_expr *splits;
3471 struct matrix_state *state;
3472 struct file_handle *file;
3473 struct stringi_set rowtypes;
3477 struct eigen_command
3479 struct matrix_expr *expr;
3480 struct matrix_var *evec;
3481 struct matrix_var *eval;
3485 struct setdiag_command
3487 struct matrix_var *dst;
3488 struct matrix_expr *expr;
3494 struct matrix_expr *expr;
3495 struct matrix_var *u;
3496 struct matrix_var *s;
3497 struct matrix_var *v;
3503 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3504 static bool matrix_cmd_execute (struct matrix_cmd *);
3505 static void matrix_cmd_destroy (struct matrix_cmd *);
3508 static struct matrix_cmd *
3509 matrix_parse_compute (struct matrix_state *s)
3511 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3512 *cmd = (struct matrix_cmd) {
3513 .type = MCMD_COMPUTE,
3514 .compute = { .lvalue = NULL },
3517 struct compute_command *compute = &cmd->compute;
3518 compute->lvalue = matrix_lvalue_parse (s);
3519 if (!compute->lvalue)
3522 if (!lex_force_match (s->lexer, T_EQUALS))
3525 compute->rvalue = matrix_parse_expr (s);
3526 if (!compute->rvalue)
3532 matrix_cmd_destroy (cmd);
3537 matrix_cmd_execute_compute (struct compute_command *compute)
3539 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3543 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3547 print_labels_destroy (struct print_labels *labels)
3551 string_array_destroy (&labels->literals);
3552 matrix_expr_destroy (labels->expr);
3558 parse_literal_print_labels (struct matrix_state *s,
3559 struct print_labels **labelsp)
3561 lex_match (s->lexer, T_EQUALS);
3562 print_labels_destroy (*labelsp);
3563 *labelsp = xzalloc (sizeof **labelsp);
3564 while (lex_token (s->lexer) != T_SLASH
3565 && lex_token (s->lexer) != T_ENDCMD
3566 && lex_token (s->lexer) != T_STOP)
3568 struct string label = DS_EMPTY_INITIALIZER;
3569 while (lex_token (s->lexer) != T_COMMA
3570 && lex_token (s->lexer) != T_SLASH
3571 && lex_token (s->lexer) != T_ENDCMD
3572 && lex_token (s->lexer) != T_STOP)
3574 if (!ds_is_empty (&label))
3575 ds_put_byte (&label, ' ');
3577 if (lex_token (s->lexer) == T_STRING)
3578 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3581 char *rep = lex_next_representation (s->lexer, 0, 0);
3582 ds_put_cstr (&label, rep);
3587 string_array_append_nocopy (&(*labelsp)->literals,
3588 ds_steal_cstr (&label));
3590 if (!lex_match (s->lexer, T_COMMA))
3596 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3598 lex_match (s->lexer, T_EQUALS);
3599 struct matrix_expr *e = matrix_parse_exp (s);
3603 print_labels_destroy (*labelsp);
3604 *labelsp = xzalloc (sizeof **labelsp);
3605 (*labelsp)->expr = e;
3609 static struct matrix_cmd *
3610 matrix_parse_print (struct matrix_state *s)
3612 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3613 *cmd = (struct matrix_cmd) {
3616 .use_default_format = true,
3620 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3623 for (size_t i = 0; ; i++)
3625 enum token_type t = lex_next_token (s->lexer, i);
3626 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3628 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3630 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3633 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3638 cmd->print.expression = matrix_parse_exp (s);
3639 if (!cmd->print.expression)
3643 while (lex_match (s->lexer, T_SLASH))
3645 if (lex_match_id (s->lexer, "FORMAT"))
3647 lex_match (s->lexer, T_EQUALS);
3648 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3650 cmd->print.use_default_format = false;
3652 else if (lex_match_id (s->lexer, "TITLE"))
3654 lex_match (s->lexer, T_EQUALS);
3655 if (!lex_force_string (s->lexer))
3657 free (cmd->print.title);
3658 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3661 else if (lex_match_id (s->lexer, "SPACE"))
3663 lex_match (s->lexer, T_EQUALS);
3664 if (lex_match_id (s->lexer, "NEWPAGE"))
3665 cmd->print.space = -1;
3666 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3668 cmd->print.space = lex_integer (s->lexer);
3674 else if (lex_match_id (s->lexer, "RLABELS"))
3675 parse_literal_print_labels (s, &cmd->print.rlabels);
3676 else if (lex_match_id (s->lexer, "CLABELS"))
3677 parse_literal_print_labels (s, &cmd->print.clabels);
3678 else if (lex_match_id (s->lexer, "RNAMES"))
3680 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3683 else if (lex_match_id (s->lexer, "CNAMES"))
3685 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3690 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3691 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3699 matrix_cmd_destroy (cmd);
3704 matrix_is_integer (const gsl_matrix *m)
3706 for (size_t y = 0; y < m->size1; y++)
3707 for (size_t x = 0; x < m->size2; x++)
3709 double d = gsl_matrix_get (m, y, x);
3717 matrix_max_magnitude (const gsl_matrix *m)
3720 for (size_t y = 0; y < m->size1; y++)
3721 for (size_t x = 0; x < m->size2; x++)
3723 double d = fabs (gsl_matrix_get (m, y, x));
3731 format_fits (struct fmt_spec format, double x)
3733 char *s = data_out (&(union value) { .f = x }, NULL,
3734 &format, settings_get_fmt_settings ());
3735 bool fits = *s != '*' && !strchr (s, 'E');
3740 static struct fmt_spec
3741 default_format (const gsl_matrix *m, int *log_scale)
3745 double max = matrix_max_magnitude (m);
3747 if (matrix_is_integer (m))
3748 for (int w = 1; w <= 12; w++)
3750 struct fmt_spec format = { .type = FMT_F, .w = w };
3751 if (format_fits (format, -max))
3755 if (max >= 1e9 || max <= 1e-4)
3758 snprintf (s, sizeof s, "%.1e", max);
3760 const char *e = strchr (s, 'e');
3762 *log_scale = atoi (e + 1);
3765 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3769 trimmed_string (double d)
3771 struct substring s = ss_buffer ((char *) &d, sizeof d);
3772 ss_rtrim (&s, ss_cstr (" "));
3773 return ss_xstrdup (s);
3776 static struct string_array *
3777 print_labels_get (const struct print_labels *labels, size_t n,
3778 const char *prefix, bool truncate)
3783 struct string_array *out = xzalloc (sizeof *out);
3784 if (labels->literals.n)
3785 string_array_clone (out, &labels->literals);
3786 else if (labels->expr)
3788 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3789 if (m && is_vector (m))
3791 gsl_vector v = to_vector (m);
3792 for (size_t i = 0; i < v.size; i++)
3793 string_array_append_nocopy (out, trimmed_string (
3794 gsl_vector_get (&v, i)));
3796 gsl_matrix_free (m);
3802 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3805 string_array_append (out, "");
3808 string_array_delete (out, out->n - 1);
3811 for (size_t i = 0; i < out->n; i++)
3813 char *s = out->strings[i];
3814 s[strnlen (s, 8)] = '\0';
3821 matrix_cmd_print_space (int space)
3824 output_item_submit (page_break_item_create ());
3825 for (int i = 0; i < space; i++)
3826 output_log ("%s", "");
3830 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3831 struct fmt_spec format, int log_scale)
3833 matrix_cmd_print_space (print->space);
3835 output_log ("%s", print->title);
3837 output_log (" 10 ** %d X", log_scale);
3839 struct string_array *clabels = print_labels_get (print->clabels,
3840 m->size2, "col", true);
3841 if (clabels && format.w < 8)
3843 struct string_array *rlabels = print_labels_get (print->rlabels,
3844 m->size1, "row", true);
3848 struct string line = DS_EMPTY_INITIALIZER;
3850 ds_put_byte_multiple (&line, ' ', 8);
3851 for (size_t x = 0; x < m->size2; x++)
3852 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3853 output_log_nocopy (ds_steal_cstr (&line));
3856 double scale = pow (10.0, log_scale);
3857 bool numeric = fmt_is_numeric (format.type);
3858 for (size_t y = 0; y < m->size1; y++)
3860 struct string line = DS_EMPTY_INITIALIZER;
3862 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3864 for (size_t x = 0; x < m->size2; x++)
3866 double f = gsl_matrix_get (m, y, x);
3868 ? data_out (&(union value) { .f = f / scale}, NULL,
3869 &format, settings_get_fmt_settings ())
3870 : trimmed_string (f));
3871 ds_put_format (&line, " %s", s);
3874 output_log_nocopy (ds_steal_cstr (&line));
3877 string_array_destroy (rlabels);
3879 string_array_destroy (clabels);
3884 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3885 const struct print_labels *print_labels, size_t n,
3886 const char *name, const char *prefix)
3888 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3890 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3891 for (size_t i = 0; i < n; i++)
3892 pivot_category_create_leaf (
3894 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3895 : pivot_value_new_integer (i + 1)));
3897 d->hide_all_labels = true;
3898 string_array_destroy (labels);
3903 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3904 struct fmt_spec format, int log_scale)
3906 struct pivot_table *table = pivot_table_create__ (
3907 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3910 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3912 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3913 N_("Columns"), "col");
3915 struct pivot_footnote *footnote = NULL;
3918 char *s = xasprintf ("× 10**%d", log_scale);
3919 footnote = pivot_table_create_footnote (
3920 table, pivot_value_new_user_text_nocopy (s));
3923 double scale = pow (10.0, log_scale);
3924 bool numeric = fmt_is_numeric (format.type);
3925 for (size_t y = 0; y < m->size1; y++)
3926 for (size_t x = 0; x < m->size2; x++)
3928 double f = gsl_matrix_get (m, y, x);
3929 struct pivot_value *v;
3932 v = pivot_value_new_number (f / scale);
3933 v->numeric.format = format;
3936 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3938 pivot_value_add_footnote (v, footnote);
3939 pivot_table_put2 (table, y, x, v);
3942 pivot_table_submit (table);
3946 matrix_cmd_execute_print (const struct print_command *print)
3948 if (print->expression)
3950 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3955 struct fmt_spec format = (print->use_default_format
3956 ? default_format (m, &log_scale)
3959 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3960 matrix_cmd_print_text (print, m, format, log_scale);
3962 matrix_cmd_print_tables (print, m, format, log_scale);
3964 gsl_matrix_free (m);
3968 matrix_cmd_print_space (print->space);
3970 output_log ("%s", print->title);
3977 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3978 const char *command_name,
3979 const char *stop1, const char *stop2)
3981 lex_end_of_command (s->lexer);
3982 lex_discard_rest_of_command (s->lexer);
3984 size_t allocated = 0;
3987 while (lex_token (s->lexer) == T_ENDCMD)
3990 if (lex_at_phrase (s->lexer, stop1)
3991 || (stop2 && lex_at_phrase (s->lexer, stop2)))
3994 if (lex_at_phrase (s->lexer, "END MATRIX"))
3996 msg (SE, _("Premature END MATRIX within %s."), command_name);
4000 struct matrix_cmd *cmd = matrix_parse_command (s);
4004 if (c->n >= allocated)
4005 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
4006 c->commands[c->n++] = cmd;
4011 matrix_parse_do_if_clause (struct matrix_state *s,
4012 struct do_if_command *ifc,
4013 bool parse_condition,
4014 size_t *allocated_clauses)
4016 if (ifc->n_clauses >= *allocated_clauses)
4017 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
4018 sizeof *ifc->clauses);
4019 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
4020 *c = (struct do_if_clause) { .condition = NULL };
4022 if (parse_condition)
4024 c->condition = matrix_parse_expr (s);
4029 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
4032 static struct matrix_cmd *
4033 matrix_parse_do_if (struct matrix_state *s)
4035 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4036 *cmd = (struct matrix_cmd) {
4038 .do_if = { .n_clauses = 0 }
4041 size_t allocated_clauses = 0;
4044 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
4047 while (lex_match_phrase (s->lexer, "ELSE IF"));
4049 if (lex_match_id (s->lexer, "ELSE")
4050 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
4053 if (!lex_match_phrase (s->lexer, "END IF"))
4058 matrix_cmd_destroy (cmd);
4063 matrix_cmd_execute_do_if (struct do_if_command *cmd)
4065 for (size_t i = 0; i < cmd->n_clauses; i++)
4067 struct do_if_clause *c = &cmd->clauses[i];
4071 if (!matrix_expr_evaluate_scalar (c->condition,
4072 i ? "ELSE IF" : "DO IF",
4077 for (size_t j = 0; j < c->commands.n; j++)
4078 if (!matrix_cmd_execute (c->commands.commands[j]))
4085 static struct matrix_cmd *
4086 matrix_parse_loop (struct matrix_state *s)
4088 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4089 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4091 struct loop_command *loop = &cmd->loop;
4092 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4094 struct substring name = lex_tokss (s->lexer);
4095 loop->var = matrix_var_lookup (s, name);
4097 loop->var = matrix_var_create (s, name);
4102 loop->start = matrix_parse_expr (s);
4103 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4106 loop->end = matrix_parse_expr (s);
4110 if (lex_match (s->lexer, T_BY))
4112 loop->increment = matrix_parse_expr (s);
4113 if (!loop->increment)
4118 if (lex_match_id (s->lexer, "IF"))
4120 loop->top_condition = matrix_parse_expr (s);
4121 if (!loop->top_condition)
4125 bool was_in_loop = s->in_loop;
4127 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4129 s->in_loop = was_in_loop;
4133 if (!lex_match_phrase (s->lexer, "END LOOP"))
4136 if (lex_match_id (s->lexer, "IF"))
4138 loop->bottom_condition = matrix_parse_expr (s);
4139 if (!loop->bottom_condition)
4146 matrix_cmd_destroy (cmd);
4151 matrix_cmd_execute_loop (struct loop_command *cmd)
4155 long int increment = 1;
4158 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4159 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4161 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4165 if (increment > 0 ? value > end
4166 : increment < 0 ? value < end
4171 int mxloops = settings_get_mxloops ();
4172 for (int i = 0; i < mxloops; i++)
4177 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4179 gsl_matrix_free (cmd->var->value);
4180 cmd->var->value = NULL;
4182 if (!cmd->var->value)
4183 cmd->var->value = gsl_matrix_alloc (1, 1);
4185 gsl_matrix_set (cmd->var->value, 0, 0, value);
4188 if (cmd->top_condition)
4191 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4196 for (size_t j = 0; j < cmd->commands.n; j++)
4197 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4200 if (cmd->bottom_condition)
4203 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4204 "END LOOP IF", &d) || d > 0)
4211 if (increment > 0 ? value > end : value < end)
4217 static struct matrix_cmd *
4218 matrix_parse_break (struct matrix_state *s)
4222 msg (SE, _("BREAK not inside LOOP."));
4226 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4227 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4231 static struct matrix_cmd *
4232 matrix_parse_display (struct matrix_state *s)
4234 bool dictionary = false;
4235 bool status = false;
4236 while (lex_token (s->lexer) == T_ID)
4238 if (lex_match_id (s->lexer, "DICTIONARY"))
4240 else if (lex_match_id (s->lexer, "STATUS"))
4244 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4248 if (!dictionary && !status)
4249 dictionary = status = true;
4251 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4252 *cmd = (struct matrix_cmd) {
4253 .type = MCMD_DISPLAY,
4256 .dictionary = dictionary,
4264 compare_matrix_var_pointers (const void *a_, const void *b_)
4266 const struct matrix_var *const *ap = a_;
4267 const struct matrix_var *const *bp = b_;
4268 const struct matrix_var *a = *ap;
4269 const struct matrix_var *b = *bp;
4270 return strcmp (a->name, b->name);
4274 matrix_cmd_execute_display (struct display_command *cmd)
4276 const struct matrix_state *s = cmd->state;
4278 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4279 pivot_dimension_create (
4280 table, PIVOT_AXIS_COLUMN, N_("Property"),
4281 N_("Rows"), N_("Columns"), N_("Size (kB)"));
4283 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4286 const struct matrix_var *var;
4287 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4289 vars[n_vars++] = var;
4290 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4292 struct pivot_dimension *rows = pivot_dimension_create (
4293 table, PIVOT_AXIS_ROW, N_("Variable"));
4294 for (size_t i = 0; i < n_vars; i++)
4296 const struct matrix_var *var = vars[i];
4297 pivot_category_create_leaf (
4298 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4300 size_t r = var->value->size1;
4301 size_t c = var->value->size2;
4302 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4303 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4304 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4306 pivot_table_submit (table);
4309 static struct matrix_cmd *
4310 matrix_parse_release (struct matrix_state *s)
4312 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4313 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4315 size_t allocated_vars = 0;
4316 while (lex_token (s->lexer) == T_ID)
4318 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4321 if (cmd->release.n_vars >= allocated_vars)
4322 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4323 sizeof *cmd->release.vars);
4324 cmd->release.vars[cmd->release.n_vars++] = var;
4327 lex_error (s->lexer, _("Variable name expected."));
4330 if (!lex_match (s->lexer, T_COMMA))
4338 matrix_cmd_execute_release (struct release_command *cmd)
4340 for (size_t i = 0; i < cmd->n_vars; i++)
4342 struct matrix_var *var = cmd->vars[i];
4343 gsl_matrix_free (var->value);
4348 static struct matrix_cmd *
4349 matrix_parse_save (struct matrix_state *s)
4351 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4352 *cmd = (struct matrix_cmd) {
4355 .strings = STRINGI_SET_INITIALIZER (cmd->save.strings)
4359 struct save_command *save = &cmd->save;
4360 save->expression = matrix_parse_exp (s);
4361 if (!save->expression)
4364 while (lex_match (s->lexer, T_SLASH))
4366 if (lex_match_id (s->lexer, "OUTFILE"))
4368 lex_match (s->lexer, T_EQUALS);
4370 fh_unref (save->outfile);
4371 save->outfile = (lex_match (s->lexer, T_ASTERISK)
4373 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4377 else if (lex_match_id (s->lexer, "VARIABLES"))
4379 lex_match (s->lexer, T_EQUALS);
4383 struct dictionary *d = dict_create (get_default_encoding ());
4384 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4385 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4390 string_array_destroy (save->variables);
4391 if (!save->variables)
4392 save->variables = xmalloc (sizeof *save->variables);
4393 *save->variables = (struct string_array) {
4399 else if (lex_match_id (s->lexer, "NAMES"))
4401 lex_match (s->lexer, T_EQUALS);
4402 matrix_expr_destroy (save->names);
4403 save->names = matrix_parse_exp (s);
4407 else if (lex_match_id (s->lexer, "STRINGS"))
4409 lex_match (s->lexer, T_EQUALS);
4410 while (lex_token (s->lexer) == T_ID)
4412 stringi_set_insert (&save->strings, lex_tokcstr (s->lexer));
4414 if (!lex_match (s->lexer, T_COMMA))
4420 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4428 if (s->prev_save_outfile)
4429 save->outfile = fh_ref (s->prev_save_outfile);
4432 lex_sbc_missing ("OUTFILE");
4436 fh_unref (s->prev_save_outfile);
4437 s->prev_save_outfile = fh_ref (save->outfile);
4439 if (save->variables && save->names)
4441 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4442 matrix_expr_destroy (save->names);
4449 matrix_cmd_destroy (cmd);
4454 matrix_cmd_execute_save (const struct save_command *save)
4456 assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */
4458 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4463 struct dictionary *dict = dict_create (get_default_encoding ());
4464 struct string_array names = { .n = 0 };
4465 if (save->variables)
4466 string_array_clone (&names, save->variables);
4467 else if (save->names)
4469 gsl_matrix *nm = matrix_expr_evaluate (save->names);
4470 if (nm && is_vector (nm))
4472 gsl_vector nv = to_vector (nm);
4473 for (size_t i = 0; i < nv.size; i++)
4475 char *name = trimmed_string (gsl_vector_get (&nv, i));
4476 if (dict_id_is_valid (dict, name, true))
4477 string_array_append_nocopy (&names, name);
4484 struct stringi_set strings;
4485 stringi_set_clone (&strings, &save->strings);
4487 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4492 name = names.strings[i];
4495 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4499 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4500 struct variable *var = dict_create_var (dict, name, width);
4503 msg (SE, _("Duplicate variable name %s in SAVE statement."),
4509 if (!stringi_set_is_empty (&strings))
4511 const char *example = stringi_set_node_get_string (
4512 stringi_set_first (&strings));
4513 msg (SE, ngettext ("STRINGS specified a variable %s, but no variable "
4514 "with that name was found on SAVE.",
4515 "STRINGS specified %2$zu variables, including %1$s, "
4516 "whose names were not found on SAVE.",
4517 stringi_set_count (&strings)),
4518 example, stringi_set_count (&strings));
4521 stringi_set_destroy (&strings);
4522 string_array_destroy (&names);
4526 gsl_matrix_free (m);
4531 struct casewriter *writer = any_writer_open (save->outfile, dict);
4534 gsl_matrix_free (m);
4539 const struct caseproto *proto = dict_get_proto (dict);
4540 for (size_t y = 0; y < m->size1; y++)
4542 struct ccase *c = case_create (proto);
4543 for (size_t x = 0; x < m->size2; x++)
4545 double d = gsl_matrix_get (m, y, x);
4546 int width = caseproto_get_width (proto, x);
4547 union value *value = case_data_rw_idx (c, x);
4551 memcpy (value->s, &d, width);
4553 casewriter_write (writer, c);
4555 casewriter_destroy (writer);
4557 gsl_matrix_free (m);
4561 static struct matrix_cmd *
4562 matrix_parse_read (struct matrix_state *s)
4564 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4565 *cmd = (struct matrix_cmd) {
4567 .read = { .format = FMT_F },
4570 struct file_handle *fh = NULL;
4571 char *encoding = NULL;
4572 struct read_command *read = &cmd->read;
4573 read->dst = matrix_lvalue_parse (s);
4578 int repetitions = 0;
4579 int record_width = 0;
4580 bool seen_format = false;
4581 while (lex_match (s->lexer, T_SLASH))
4583 if (lex_match_id (s->lexer, "FILE"))
4585 lex_match (s->lexer, T_EQUALS);
4588 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4592 else if (lex_match_id (s->lexer, "ENCODING"))
4594 lex_match (s->lexer, T_EQUALS);
4595 if (!lex_force_string (s->lexer))
4599 encoding = ss_xstrdup (lex_tokss (s->lexer));
4603 else if (lex_match_id (s->lexer, "FIELD"))
4605 lex_match (s->lexer, T_EQUALS);
4607 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4609 read->c1 = lex_integer (s->lexer);
4611 if (!lex_force_match (s->lexer, T_TO)
4612 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4614 read->c2 = lex_integer (s->lexer) + 1;
4617 record_width = read->c2 - read->c1;
4618 if (lex_match (s->lexer, T_BY))
4620 if (!lex_force_int_range (s->lexer, "BY", 1,
4621 read->c2 - read->c1))
4623 by = lex_integer (s->lexer);
4626 if (record_width % by)
4628 msg (SE, _("BY %d does not evenly divide record width %d."),
4636 else if (lex_match_id (s->lexer, "SIZE"))
4638 lex_match (s->lexer, T_EQUALS);
4639 matrix_expr_destroy (read->size);
4640 read->size = matrix_parse_exp (s);
4644 else if (lex_match_id (s->lexer, "MODE"))
4646 lex_match (s->lexer, T_EQUALS);
4647 if (lex_match_id (s->lexer, "RECTANGULAR"))
4648 read->symmetric = false;
4649 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4650 read->symmetric = true;
4653 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4657 else if (lex_match_id (s->lexer, "REREAD"))
4658 read->reread = true;
4659 else if (lex_match_id (s->lexer, "FORMAT"))
4663 lex_sbc_only_once ("FORMAT");
4668 lex_match (s->lexer, T_EQUALS);
4670 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4673 const char *p = lex_tokcstr (s->lexer);
4674 if (c_isdigit (p[0]))
4676 repetitions = atoi (p);
4677 p += strspn (p, "0123456789");
4678 if (!fmt_from_name (p, &read->format))
4680 lex_error (s->lexer, _("Unknown format %s."), p);
4685 else if (fmt_from_name (p, &read->format))
4689 struct fmt_spec format;
4690 if (!parse_format_specifier (s->lexer, &format))
4692 read->format = format.type;
4698 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4699 "REREAD", "FORMAT");
4706 lex_sbc_missing ("FIELD");
4710 if (!read->dst->n_indexes && !read->size)
4712 msg (SE, _("SIZE is required for reading data into a full matrix "
4713 "(as opposed to a submatrix)."));
4719 if (s->prev_read_file)
4720 fh = fh_ref (s->prev_read_file);
4723 lex_sbc_missing ("FILE");
4727 fh_unref (s->prev_read_file);
4728 s->prev_read_file = fh_ref (fh);
4730 read->rf = read_file_create (s, fh);
4734 free (read->rf->encoding);
4735 read->rf->encoding = encoding;
4739 /* Field width may be specified in multiple ways:
4742 2. The format on FORMAT.
4743 3. The repetition factor on FORMAT.
4745 (2) and (3) are mutually exclusive.
4747 If more than one of these is present, they must agree. If none of them is
4748 present, then free-field format is used.
4750 if (repetitions > record_width)
4752 msg (SE, _("%d repetitions cannot fit in record width %d."),
4753 repetitions, record_width);
4756 int w = (repetitions ? record_width / repetitions
4762 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4763 "which implies field width %d, "
4764 "but BY specifies field width %d."),
4765 repetitions, record_width, w, by);
4767 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4776 matrix_cmd_destroy (cmd);
4782 parse_error (const struct dfm_reader *reader, enum fmt_type format,
4783 struct substring data, size_t y, size_t x,
4784 int first_column, int last_column, char *error)
4786 int line_number = dfm_get_line_number (reader);
4787 struct msg_location *location = xmalloc (sizeof *location);
4788 *location = (struct msg_location) {
4789 .file_name = xstrdup (dfm_get_file_name (reader)),
4790 .first_line = line_number,
4791 .last_line = line_number + 1,
4792 .first_column = first_column,
4793 .last_column = last_column,
4795 struct msg *m = xmalloc (sizeof *m);
4797 .category = MSG_C_DATA,
4798 .severity = MSG_S_WARNING,
4799 .location = location,
4800 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
4801 "for matrix row %zu, column %zu: %s"),
4802 (int) data.length, data.string, fmt_name (format),
4803 y + 1, x + 1, error),
4811 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4812 gsl_matrix *m, struct substring p, size_t y, size_t x,
4813 const char *line_start)
4815 const char *input_encoding = dfm_reader_get_encoding (reader);
4818 if (fmt_is_numeric (read->format))
4821 error = data_in (p, input_encoding, read->format,
4822 settings_get_fmt_settings (), &v, 0, NULL);
4823 if (!error && v.f == SYSMIS)
4824 error = xstrdup (_("Matrix data may not contain missing value."));
4829 uint8_t s[sizeof (double)];
4830 union value v = { .s = s };
4831 error = data_in (p, input_encoding, read->format,
4832 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
4833 memcpy (&f, s, sizeof f);
4838 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
4839 int c2 = c1 + ss_utf8_count_columns (p) - 1;
4840 parse_error (reader, read->format, p, y, x, c1, c2, error);
4844 gsl_matrix_set (m, y, x, f);
4845 if (read->symmetric && x != y)
4846 gsl_matrix_set (m, x, y, f);
4851 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4852 struct substring *line, const char **startp)
4854 if (dfm_eof (reader))
4856 msg (SE, _("Unexpected end of file reading matrix data."));
4859 dfm_expand_tabs (reader);
4860 struct substring record = dfm_get_record (reader);
4861 /* XXX need to recode record into UTF-8 */
4862 *startp = record.string;
4863 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
4868 matrix_read (struct read_command *read, struct dfm_reader *reader,
4871 for (size_t y = 0; y < m->size1; y++)
4873 size_t nx = read->symmetric ? y + 1 : m->size2;
4875 struct substring line = ss_empty ();
4876 const char *line_start = line.string;
4877 for (size_t x = 0; x < nx; x++)
4884 ss_ltrim (&line, ss_cstr (" ,"));
4885 if (!ss_is_empty (line))
4887 if (!matrix_read_line (read, reader, &line, &line_start))
4889 dfm_forward_record (reader);
4892 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
4896 if (!matrix_read_line (read, reader, &line, &line_start))
4898 size_t fields_per_line = (read->c2 - read->c1) / read->w;
4899 int f = x % fields_per_line;
4900 if (f == fields_per_line - 1)
4901 dfm_forward_record (reader);
4903 p = ss_substr (line, read->w * f, read->w);
4906 matrix_read_set_field (read, reader, m, p, y, x, line_start);
4910 dfm_forward_record (reader);
4913 ss_ltrim (&line, ss_cstr (" ,"));
4914 if (!ss_is_empty (line))
4917 msg (SW, _("Trailing garbage on line \"%.*s\""),
4918 (int) line.length, line.string);
4925 matrix_cmd_execute_read (struct read_command *read)
4927 struct index_vector iv0, iv1;
4928 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
4931 size_t size[2] = { SIZE_MAX, SIZE_MAX };
4934 gsl_matrix *m = matrix_expr_evaluate (read->size);
4940 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4941 "not a %zu×%zu matrix."), m->size1, m->size2);
4942 gsl_matrix_free (m);
4948 gsl_vector v = to_vector (m);
4952 d[0] = gsl_vector_get (&v, 0);
4955 else if (v.size == 2)
4957 d[0] = gsl_vector_get (&v, 0);
4958 d[1] = gsl_vector_get (&v, 1);
4962 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4963 "not a %zu×%zu matrix."),
4964 m->size1, m->size2),
4965 gsl_matrix_free (m);
4970 gsl_matrix_free (m);
4972 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
4974 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
4975 "are outside valid range."),
4986 if (read->dst->n_indexes)
4988 size_t submatrix_size[2];
4989 if (read->dst->n_indexes == 2)
4991 submatrix_size[0] = iv0.n;
4992 submatrix_size[1] = iv1.n;
4994 else if (read->dst->var->value->size1 == 1)
4996 submatrix_size[0] = 1;
4997 submatrix_size[1] = iv0.n;
5001 submatrix_size[0] = iv0.n;
5002 submatrix_size[1] = 1;
5007 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
5009 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
5010 "differ from submatrix dimensions %zu×%zu."),
5012 submatrix_size[0], submatrix_size[1]);
5020 size[0] = submatrix_size[0];
5021 size[1] = submatrix_size[1];
5025 struct dfm_reader *reader = read_file_open (read->rf);
5027 dfm_reread_record (reader, 1);
5029 if (read->symmetric && size[0] != size[1])
5031 msg (SE, _("Cannot read non-square %zu×%zu matrix "
5032 "using READ with MODE=SYMMETRIC."),
5038 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
5039 matrix_read (read, reader, tmp);
5040 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
5043 static struct matrix_cmd *
5044 matrix_parse_write (struct matrix_state *s)
5046 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5047 *cmd = (struct matrix_cmd) {
5051 struct file_handle *fh = NULL;
5052 char *encoding = NULL;
5053 struct write_command *write = &cmd->write;
5054 write->expression = matrix_parse_exp (s);
5055 if (!write->expression)
5059 int repetitions = 0;
5060 int record_width = 0;
5061 enum fmt_type format = FMT_F;
5062 bool has_format = false;
5063 while (lex_match (s->lexer, T_SLASH))
5065 if (lex_match_id (s->lexer, "OUTFILE"))
5067 lex_match (s->lexer, T_EQUALS);
5070 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5074 else if (lex_match_id (s->lexer, "ENCODING"))
5076 lex_match (s->lexer, T_EQUALS);
5077 if (!lex_force_string (s->lexer))
5081 encoding = ss_xstrdup (lex_tokss (s->lexer));
5085 else if (lex_match_id (s->lexer, "FIELD"))
5087 lex_match (s->lexer, T_EQUALS);
5089 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5091 write->c1 = lex_integer (s->lexer);
5093 if (!lex_force_match (s->lexer, T_TO)
5094 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5096 write->c2 = lex_integer (s->lexer) + 1;
5099 record_width = write->c2 - write->c1;
5100 if (lex_match (s->lexer, T_BY))
5102 if (!lex_force_int_range (s->lexer, "BY", 1,
5103 write->c2 - write->c1))
5105 by = lex_integer (s->lexer);
5108 if (record_width % by)
5110 msg (SE, _("BY %d does not evenly divide record width %d."),
5118 else if (lex_match_id (s->lexer, "MODE"))
5120 lex_match (s->lexer, T_EQUALS);
5121 if (lex_match_id (s->lexer, "RECTANGULAR"))
5122 write->triangular = false;
5123 else if (lex_match_id (s->lexer, "TRIANGULAR"))
5124 write->triangular = true;
5127 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5131 else if (lex_match_id (s->lexer, "HOLD"))
5133 else if (lex_match_id (s->lexer, "FORMAT"))
5135 if (has_format || write->format)
5137 lex_sbc_only_once ("FORMAT");
5141 lex_match (s->lexer, T_EQUALS);
5143 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5146 const char *p = lex_tokcstr (s->lexer);
5147 if (c_isdigit (p[0]))
5149 repetitions = atoi (p);
5150 p += strspn (p, "0123456789");
5151 if (!fmt_from_name (p, &format))
5153 lex_error (s->lexer, _("Unknown format %s."), p);
5159 else if (fmt_from_name (p, &format))
5166 struct fmt_spec spec;
5167 if (!parse_format_specifier (s->lexer, &spec))
5169 write->format = xmemdup (&spec, sizeof spec);
5174 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5182 lex_sbc_missing ("FIELD");
5188 if (s->prev_write_file)
5189 fh = fh_ref (s->prev_write_file);
5192 lex_sbc_missing ("OUTFILE");
5196 fh_unref (s->prev_write_file);
5197 s->prev_write_file = fh_ref (fh);
5199 write->wf = write_file_create (s, fh);
5203 free (write->wf->encoding);
5204 write->wf->encoding = encoding;
5208 /* Field width may be specified in multiple ways:
5211 2. The format on FORMAT.
5212 3. The repetition factor on FORMAT.
5214 (2) and (3) are mutually exclusive.
5216 If more than one of these is present, they must agree. If none of them is
5217 present, then free-field format is used.
5219 if (repetitions > record_width)
5221 msg (SE, _("%d repetitions cannot fit in record width %d."),
5222 repetitions, record_width);
5225 int w = (repetitions ? record_width / repetitions
5226 : write->format ? write->format->w
5231 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5232 "which implies field width %d, "
5233 "but BY specifies field width %d."),
5234 repetitions, record_width, w, by);
5236 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5240 if (w && !write->format)
5242 write->format = xmalloc (sizeof *write->format);
5243 *write->format = (struct fmt_spec) { .type = format, .w = w };
5245 if (!fmt_check_output (write->format))
5249 if (write->format && fmt_var_width (write->format) > sizeof (double))
5251 char s[FMT_STRING_LEN_MAX + 1];
5252 fmt_to_string (write->format, s);
5253 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
5254 s, sizeof (double));
5262 matrix_cmd_destroy (cmd);
5267 matrix_cmd_execute_write (struct write_command *write)
5269 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5273 if (write->triangular && m->size1 != m->size2)
5275 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5276 "the matrix to be written has dimensions %zu×%zu."),
5277 m->size1, m->size2);
5278 gsl_matrix_free (m);
5282 struct dfm_writer *writer = write_file_open (write->wf);
5283 if (!writer || !m->size1)
5285 gsl_matrix_free (m);
5289 const struct fmt_settings *settings = settings_get_fmt_settings ();
5290 struct u8_line *line = write->wf->held;
5291 for (size_t y = 0; y < m->size1; y++)
5295 line = xmalloc (sizeof *line);
5296 u8_line_init (line);
5298 size_t nx = write->triangular ? y + 1 : m->size2;
5300 for (size_t x = 0; x < nx; x++)
5303 double f = gsl_matrix_get (m, y, x);
5307 if (fmt_is_numeric (write->format->type))
5310 v.s = (uint8_t *) &f;
5311 s = data_out (&v, NULL, write->format, settings);
5315 s = xmalloc (DBL_BUFSIZE_BOUND);
5316 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
5317 >= DBL_BUFSIZE_BOUND)
5320 size_t len = strlen (s);
5321 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5322 if (width + x0 > write->c2)
5324 dfm_put_record_utf8 (writer, line->s.ss.string,
5326 u8_line_clear (line);
5329 u8_line_put (line, x0, x0 + width, s, len);
5332 x0 += write->format ? write->format->w : width + 1;
5335 if (y + 1 >= m->size1 && write->hold)
5337 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5338 u8_line_clear (line);
5342 u8_line_destroy (line);
5346 write->wf->held = line;
5348 gsl_matrix_free (m);
5351 static struct matrix_cmd *
5352 matrix_parse_get (struct matrix_state *s)
5354 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5355 *cmd = (struct matrix_cmd) {
5358 .dataset = s->dataset,
5359 .user = { .treatment = MGET_ERROR },
5360 .system = { .treatment = MGET_ERROR },
5364 struct get_command *get = &cmd->get;
5365 get->dst = matrix_lvalue_parse (s);
5369 while (lex_match (s->lexer, T_SLASH))
5371 if (lex_match_id (s->lexer, "FILE"))
5373 if (get->variables.n)
5375 lex_error (s->lexer, _("FILE must precede VARIABLES"));
5378 lex_match (s->lexer, T_EQUALS);
5380 fh_unref (get->file);
5381 if (lex_match (s->lexer, T_ASTERISK))
5385 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5390 else if (lex_match_id (s->lexer, "ENCODING"))
5392 if (get->variables.n)
5394 lex_error (s->lexer, _("ENCODING must precede VARIABLES"));
5397 lex_match (s->lexer, T_EQUALS);
5398 if (!lex_force_string (s->lexer))
5401 free (get->encoding);
5402 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5406 else if (lex_match_id (s->lexer, "VARIABLES"))
5408 lex_match (s->lexer, T_EQUALS);
5410 struct dictionary *dict = NULL;
5413 dict = dict_ref (dataset_dict (s->dataset));
5414 if (dict_get_var_cnt (dict) == 0)
5416 lex_error (s->lexer, _("GET cannot read empty active file."));
5422 struct casereader *reader = any_reader_open_and_decode (
5423 get->file, get->encoding, &dict, NULL);
5426 casereader_destroy (reader);
5429 struct variable **vars;
5431 bool ok = parse_variables (s->lexer, dict, &vars, &n_vars,
5432 PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH);
5439 string_array_clear (&get->variables);
5440 for (size_t i = 0; i < n_vars; i++)
5441 string_array_append (&get->variables, var_get_name (vars[i]));
5445 else if (lex_match_id (s->lexer, "NAMES"))
5447 lex_match (s->lexer, T_EQUALS);
5448 if (!lex_force_id (s->lexer))
5451 struct substring name = lex_tokss (s->lexer);
5452 get->names = matrix_var_lookup (s, name);
5454 get->names = matrix_var_create (s, name);
5457 else if (lex_match_id (s->lexer, "MISSING"))
5459 lex_match (s->lexer, T_EQUALS);
5460 if (lex_match_id (s->lexer, "ACCEPT"))
5461 get->user.treatment = MGET_ACCEPT;
5462 else if (lex_match_id (s->lexer, "OMIT"))
5463 get->user.treatment = MGET_OMIT;
5464 else if (lex_is_number (s->lexer))
5466 get->user.treatment = MGET_RECODE;
5467 get->user.substitute = lex_number (s->lexer);
5472 lex_error (s->lexer, NULL);
5476 else if (lex_match_id (s->lexer, "SYSMIS"))
5478 lex_match (s->lexer, T_EQUALS);
5479 if (lex_match_id (s->lexer, "OMIT"))
5480 get->system.treatment = MGET_OMIT;
5481 else if (lex_is_number (s->lexer))
5483 get->system.treatment = MGET_RECODE;
5484 get->system.substitute = lex_number (s->lexer);
5489 lex_error (s->lexer, NULL);
5495 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5496 "MISSING", "SYSMIS");
5501 if (get->user.treatment != MGET_ACCEPT)
5502 get->system.treatment = MGET_ERROR;
5507 matrix_cmd_destroy (cmd);
5512 matrix_cmd_execute_get__ (struct get_command *get,
5513 const struct dictionary *dict,
5514 struct casereader *reader)
5516 const struct variable **vars = xnmalloc (
5517 get->variables.n ? get->variables.n : dict_get_var_cnt (dict),
5521 if (get->variables.n)
5523 for (size_t i = 0; i < get->variables.n; i++)
5525 const char *name = get->variables.strings[i];
5526 const struct variable *var = dict_lookup_var (dict, name);
5529 msg (SE, _("GET: Data file does not contain variable %s."),
5534 if (!var_is_numeric (var))
5536 msg (SE, _("GET: Variable %s is not numeric."), name);
5540 vars[n_vars++] = var;
5545 for (size_t i = 0; i < dict_get_var_cnt (dict); i++)
5547 const struct variable *var = dict_get_var (dict, i);
5548 if (!var_is_numeric (var))
5550 msg (SE, _("GET: Variable %s is not numeric."),
5551 var_get_name (var));
5555 vars[n_vars++] = var;
5561 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
5562 for (size_t i = 0; i < n_vars; i++)
5564 char s[sizeof (double)];
5566 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
5567 memcpy (&f, s, sizeof f);
5568 gsl_matrix_set (names, i, 0, f);
5571 gsl_matrix_free (get->names->value);
5572 get->names->value = names;
5576 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5577 long long int casenum = 1;
5579 for (struct ccase *c = casereader_read (reader); c;
5580 c = casereader_read (reader), casenum++)
5582 if (n_rows >= m->size1)
5584 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5585 for (size_t y = 0; y < n_rows; y++)
5586 for (size_t x = 0; x < n_vars; x++)
5587 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5588 gsl_matrix_free (m);
5593 for (size_t x = 0; x < n_vars; x++)
5595 const struct variable *var = vars[x];
5596 double d = case_num (c, var);
5599 if (get->system.treatment == MGET_RECODE)
5600 d = get->system.substitute;
5601 else if (get->system.treatment == MGET_OMIT)
5605 msg (SE, _("GET: Variable %s in case %lld "
5606 "is system-missing."),
5607 var_get_name (var), casenum);
5611 else if (var_is_num_missing (var, d, MV_USER))
5613 if (get->user.treatment == MGET_RECODE)
5614 d = get->user.substitute;
5615 else if (get->user.treatment == MGET_OMIT)
5617 else if (get->user.treatment != MGET_ACCEPT)
5619 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5621 var_get_name (var), casenum, d);
5625 gsl_matrix_set (m, n_rows, x, d);
5636 matrix_lvalue_evaluate_and_assign (get->dst, m);
5639 gsl_matrix_free (m);
5644 matrix_cmd_execute_get (struct get_command *get)
5646 struct dictionary *dict;
5647 struct casereader *reader;
5650 reader = any_reader_open_and_decode (get->file, get->encoding,
5657 reader = proc_open (get->dataset);
5658 dict = dict_ref (dataset_dict (get->dataset));
5661 matrix_cmd_execute_get__ (get, dict, reader);
5664 casereader_destroy (reader);
5666 proc_commit (get->dataset);
5670 match_rowtype (struct lexer *lexer)
5672 static const char *rowtypes[] = {
5673 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5675 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5677 for (size_t i = 0; i < n_rowtypes; i++)
5678 if (lex_match_id (lexer, rowtypes[i]))
5681 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5686 parse_var_names (struct lexer *lexer, struct string_array *sa)
5688 lex_match (lexer, T_EQUALS);
5690 string_array_clear (sa);
5692 struct dictionary *dict = dict_create (get_default_encoding ());
5693 dict_create_var_assert (dict, "ROWTYPE_", 8);
5694 dict_create_var_assert (dict, "VARNAME_", 8);
5697 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5698 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5703 string_array_clear (sa);
5704 sa->strings = names;
5705 sa->n = sa->allocated = n_names;
5711 msave_common_uninit (struct msave_common *common)
5715 fh_unref (common->outfile);
5716 string_array_destroy (&common->variables);
5717 string_array_destroy (&common->fnames);
5718 string_array_destroy (&common->snames);
5723 compare_variables (const char *keyword,
5724 const struct string_array *new,
5725 const struct string_array *old)
5731 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5732 "on the first MSAVE within MATRIX."), keyword);
5735 else if (!string_array_equal_case (old, new))
5737 msg (SE, _("%s must specify the same variables each time within "
5738 "a given MATRIX."), keyword);
5745 static struct matrix_cmd *
5746 matrix_parse_msave (struct matrix_state *s)
5748 struct msave_common common = { .outfile = NULL };
5749 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5750 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5752 struct msave_command *msave = &cmd->msave;
5753 if (lex_token (s->lexer) == T_ID)
5754 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5755 msave->expr = matrix_parse_exp (s);
5759 while (lex_match (s->lexer, T_SLASH))
5761 if (lex_match_id (s->lexer, "TYPE"))
5763 lex_match (s->lexer, T_EQUALS);
5765 msave->rowtype = match_rowtype (s->lexer);
5766 if (!msave->rowtype)
5769 else if (lex_match_id (s->lexer, "OUTFILE"))
5771 lex_match (s->lexer, T_EQUALS);
5773 fh_unref (common.outfile);
5774 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5775 if (!common.outfile)
5778 else if (lex_match_id (s->lexer, "VARIABLES"))
5780 if (!parse_var_names (s->lexer, &common.variables))
5783 else if (lex_match_id (s->lexer, "FNAMES"))
5785 if (!parse_var_names (s->lexer, &common.fnames))
5788 else if (lex_match_id (s->lexer, "SNAMES"))
5790 if (!parse_var_names (s->lexer, &common.snames))
5793 else if (lex_match_id (s->lexer, "SPLIT"))
5795 lex_match (s->lexer, T_EQUALS);
5797 matrix_expr_destroy (msave->splits);
5798 msave->splits = matrix_parse_exp (s);
5802 else if (lex_match_id (s->lexer, "FACTOR"))
5804 lex_match (s->lexer, T_EQUALS);
5806 matrix_expr_destroy (msave->factors);
5807 msave->factors = matrix_parse_exp (s);
5808 if (!msave->factors)
5813 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5814 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5818 if (!msave->rowtype)
5820 lex_sbc_missing ("TYPE");
5823 common.has_splits = msave->splits || common.snames.n;
5824 common.has_factors = msave->factors || common.fnames.n;
5826 struct msave_common *c = s->common ? s->common : &common;
5827 if (c->fnames.n && !msave->factors)
5829 msg (SE, _("FNAMES requires FACTOR."));
5832 if (c->snames.n && !msave->splits)
5834 msg (SE, _("SNAMES requires SPLIT."));
5837 if (c->has_factors && !common.has_factors)
5839 msg (SE, _("%s is required because it was present on the first "
5840 "MSAVE in this MATRIX command."), "FACTOR");
5843 if (c->has_splits && !common.has_splits)
5845 msg (SE, _("%s is required because it was present on the first "
5846 "MSAVE in this MATRIX command."), "SPLIT");
5852 if (!common.outfile)
5854 lex_sbc_missing ("OUTFILE");
5857 s->common = xmemdup (&common, sizeof common);
5863 bool same = common.outfile == s->common->outfile;
5864 fh_unref (common.outfile);
5867 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5868 "within a single MATRIX command."));
5872 if (!compare_variables ("VARIABLES",
5873 &common.variables, &s->common->variables)
5874 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5875 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5877 msave_common_uninit (&common);
5879 msave->common = s->common;
5880 if (!msave->varname_)
5881 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5885 msave_common_uninit (&common);
5886 matrix_cmd_destroy (cmd);
5891 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5893 gsl_matrix *m = matrix_expr_evaluate (e);
5899 msg (SE, _("%s expression must evaluate to vector, "
5900 "not a %zu×%zu matrix."),
5901 name, m->size1, m->size2);
5902 gsl_matrix_free (m);
5906 return matrix_to_vector (m);
5910 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5912 for (size_t i = 0; i < vars->n; i++)
5913 if (!dict_create_var (d, vars->strings[i], 0))
5914 return vars->strings[i];
5918 static struct dictionary *
5919 msave_create_dict (const struct msave_common *common)
5921 struct dictionary *dict = dict_create (get_default_encoding ());
5923 const char *dup_split = msave_add_vars (dict, &common->snames);
5926 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5930 dict_create_var_assert (dict, "ROWTYPE_", 8);
5932 const char *dup_factor = msave_add_vars (dict, &common->fnames);
5935 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5939 dict_create_var_assert (dict, "VARNAME_", 8);
5941 const char *dup_var = msave_add_vars (dict, &common->variables);
5944 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
5956 matrix_cmd_execute_msave (struct msave_command *msave)
5958 struct msave_common *common = msave->common;
5959 gsl_matrix *m = NULL;
5960 gsl_vector *factors = NULL;
5961 gsl_vector *splits = NULL;
5963 m = matrix_expr_evaluate (msave->expr);
5967 if (!common->variables.n)
5968 for (size_t i = 0; i < m->size2; i++)
5969 string_array_append_nocopy (&common->variables,
5970 xasprintf ("COL%zu", i + 1));
5972 if (m->size2 != common->variables.n)
5974 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
5975 m->size2, common->variables.n);
5981 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
5985 if (!common->fnames.n)
5986 for (size_t i = 0; i < factors->size; i++)
5987 string_array_append_nocopy (&common->fnames,
5988 xasprintf ("FAC%zu", i + 1));
5990 if (factors->size != common->fnames.n)
5992 msg (SE, _("There are %zu factor variables, "
5993 "but %zu split values were supplied."),
5994 common->fnames.n, factors->size);
6001 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
6005 if (!common->fnames.n)
6006 for (size_t i = 0; i < splits->size; i++)
6007 string_array_append_nocopy (&common->fnames,
6008 xasprintf ("SPL%zu", i + 1));
6010 if (splits->size != common->fnames.n)
6012 msg (SE, _("There are %zu split variables, "
6013 "but %zu split values were supplied."),
6014 common->fnames.n, splits->size);
6019 if (!common->writer)
6021 struct dictionary *dict = msave_create_dict (common);
6025 common->writer = any_writer_open (common->outfile, dict);
6026 if (!common->writer)
6032 common->dict = dict;
6035 for (size_t y = 0; y < m->size1; y++)
6037 struct ccase *c = case_create (dict_get_proto (common->dict));
6040 /* Split variables */
6042 for (size_t i = 0; i < splits->size; i++)
6043 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
6046 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6047 msave->rowtype, ' ');
6051 for (size_t i = 0; i < factors->size; i++)
6052 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
6055 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6056 msave->varname_, ' ');
6058 /* Continuous variables. */
6059 for (size_t x = 0; x < m->size2; x++)
6060 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
6061 casewriter_write (common->writer, c);
6065 gsl_matrix_free (m);
6066 gsl_vector_free (factors);
6067 gsl_vector_free (splits);
6070 static struct matrix_cmd *
6071 matrix_parse_mget (struct matrix_state *s)
6073 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6074 *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
6076 struct mget_command *mget = &cmd->mget;
6080 if (lex_match_id (s->lexer, "FILE"))
6082 lex_match (s->lexer, T_EQUALS);
6084 fh_unref (mget->file);
6085 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6089 else if (lex_match_id (s->lexer, "TYPE"))
6091 lex_match (s->lexer, T_EQUALS);
6092 while (lex_token (s->lexer) != T_SLASH
6093 && lex_token (s->lexer) != T_ENDCMD)
6095 const char *rowtype = match_rowtype (s->lexer);
6099 stringi_set_insert (&mget->rowtypes, rowtype);
6104 lex_error_expecting (s->lexer, "FILE", "TYPE");
6107 if (lex_token (s->lexer) == T_ENDCMD)
6110 if (!lex_force_match (s->lexer, T_SLASH))
6116 matrix_cmd_destroy (cmd);
6120 static const struct variable *
6121 get_a8_var (const struct dictionary *d, const char *name)
6123 const struct variable *v = dict_lookup_var (d, name);
6126 msg (SE, _("Matrix data file lacks %s variable."), name);
6129 if (var_get_width (v) != 8)
6131 msg (SE, _("%s variable in matrix data file must be "
6132 "8-byte string, but it has width %d."),
6133 name, var_get_width (v));
6140 is_rowtype (const union value *v, const char *rowtype)
6142 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
6143 ss_rtrim (&vs, ss_cstr (" "));
6144 return ss_equals_case (vs, ss_cstr (rowtype));
6148 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6149 const struct dictionary *d,
6150 const struct variable *rowtype_var,
6151 struct matrix_state *s, size_t si, size_t fi,
6152 size_t cs, size_t cn)
6157 const union value *rowtype_ = case_data (rows[0], rowtype_var);
6158 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
6159 : is_rowtype (rowtype_, "CORR") ? "CR"
6160 : is_rowtype (rowtype_, "MEAN") ? "MN"
6161 : is_rowtype (rowtype_, "STDDEV") ? "SD"
6162 : is_rowtype (rowtype_, "N") ? "NC"
6165 struct string name = DS_EMPTY_INITIALIZER;
6166 ds_put_cstr (&name, name_prefix);
6168 ds_put_format (&name, "F%zu", fi);
6170 ds_put_format (&name, "S%zu", si);
6172 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6174 mv = matrix_var_create (s, ds_ss (&name));
6177 msg (SW, _("Matrix data file contains variable with existing name %s."),
6182 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6183 size_t n_missing = 0;
6184 for (size_t y = 0; y < n_rows; y++)
6186 for (size_t x = 0; x < cn; x++)
6188 struct variable *var = dict_get_var (d, cs + x);
6189 double value = case_num (rows[y], var);
6190 if (var_is_num_missing (var, value, MV_ANY))
6195 gsl_matrix_set (m, y, x, value);
6200 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6201 "value, which was treated as zero.",
6202 "Matrix data file variable %s contains %zu missing "
6203 "values, which were treated as zero.", n_missing),
6204 ds_cstr (&name), n_missing);
6209 for (size_t y = 0; y < n_rows; y++)
6210 case_unref (rows[y]);
6214 var_changed (const struct ccase *ca, const struct ccase *cb,
6215 const struct variable *var)
6217 return !value_equal (case_data (ca, var), case_data (cb, var),
6218 var_get_width (var));
6222 vars_changed (const struct ccase *ca, const struct ccase *cb,
6223 const struct dictionary *d,
6224 size_t first_var, size_t n_vars)
6226 for (size_t i = 0; i < n_vars; i++)
6228 const struct variable *v = dict_get_var (d, first_var + i);
6229 if (var_changed (ca, cb, v))
6236 matrix_cmd_execute_mget (struct mget_command *mget)
6238 struct dictionary *d;
6239 struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6244 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6245 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6246 if (!rowtype_ || !varname_)
6249 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6251 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6254 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6256 msg (SE, _("Matrix data file contains no continuous variables."));
6260 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6262 const struct variable *v = dict_get_var (d, i);
6263 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6266 _("Matrix data file contains unexpected string variable %s."),
6272 /* SPLIT variables. */
6274 size_t sn = var_get_dict_index (rowtype_);
6275 struct ccase *sc = NULL;
6278 /* FACTOR variables. */
6279 size_t fs = var_get_dict_index (rowtype_) + 1;
6280 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6281 struct ccase *fc = NULL;
6284 /* Continuous variables. */
6285 size_t cs = var_get_dict_index (varname_) + 1;
6286 size_t cn = dict_get_var_cnt (d) - cs;
6287 struct ccase *cc = NULL;
6290 struct ccase **rows = NULL;
6291 size_t allocated_rows = 0;
6295 while ((c = casereader_read (r)) != NULL)
6297 bool sd = vars_changed (sc, c, d, ss, sn);
6298 bool fd = sd || vars_changed (fc, c, d, fs, fn);
6299 bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6314 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6315 mget->state, si, fi, cs, cn);
6321 if (n_rows >= allocated_rows)
6322 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6325 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6326 mget->state, si, fi, cs, cn);
6330 casereader_destroy (r);
6334 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6336 if (!lex_force_id (s->lexer))
6339 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6341 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6346 static struct matrix_cmd *
6347 matrix_parse_eigen (struct matrix_state *s)
6349 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6350 *cmd = (struct matrix_cmd) {
6352 .eigen = { .expr = NULL }
6355 struct eigen_command *eigen = &cmd->eigen;
6356 if (!lex_force_match (s->lexer, T_LPAREN))
6358 eigen->expr = matrix_parse_expr (s);
6360 || !lex_force_match (s->lexer, T_COMMA)
6361 || !matrix_parse_dst_var (s, &eigen->evec)
6362 || !lex_force_match (s->lexer, T_COMMA)
6363 || !matrix_parse_dst_var (s, &eigen->eval)
6364 || !lex_force_match (s->lexer, T_RPAREN))
6370 matrix_cmd_destroy (cmd);
6375 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6377 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6380 if (!is_symmetric (A))
6382 msg (SE, _("Argument of EIGEN must be symmetric."));
6383 gsl_matrix_free (A);
6387 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6388 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6389 gsl_vector v_eval = to_vector (eval);
6390 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6391 gsl_eigen_symmv (A, &v_eval, evec, w);
6392 gsl_eigen_symmv_free (w);
6394 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6396 gsl_matrix_free (A);
6398 gsl_matrix_free (eigen->eval->value);
6399 eigen->eval->value = eval;
6401 gsl_matrix_free (eigen->evec->value);
6402 eigen->evec->value = evec;
6405 static struct matrix_cmd *
6406 matrix_parse_setdiag (struct matrix_state *s)
6408 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6409 *cmd = (struct matrix_cmd) {
6410 .type = MCMD_SETDIAG,
6411 .setdiag = { .dst = NULL }
6414 struct setdiag_command *setdiag = &cmd->setdiag;
6415 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6418 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6421 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6426 if (!lex_force_match (s->lexer, T_COMMA))
6429 setdiag->expr = matrix_parse_expr (s);
6433 if (!lex_force_match (s->lexer, T_RPAREN))
6439 matrix_cmd_destroy (cmd);
6444 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6446 gsl_matrix *dst = setdiag->dst->value;
6449 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6450 setdiag->dst->name);
6454 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6458 size_t n = MIN (dst->size1, dst->size2);
6459 if (is_scalar (src))
6461 double d = to_scalar (src);
6462 for (size_t i = 0; i < n; i++)
6463 gsl_matrix_set (dst, i, i, d);
6465 else if (is_vector (src))
6467 gsl_vector v = to_vector (src);
6468 for (size_t i = 0; i < n && i < v.size; i++)
6469 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6472 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6473 "not a %zu×%zu matrix."),
6474 src->size1, src->size2);
6475 gsl_matrix_free (src);
6478 static struct matrix_cmd *
6479 matrix_parse_svd (struct matrix_state *s)
6481 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6482 *cmd = (struct matrix_cmd) {
6484 .svd = { .expr = NULL }
6487 struct svd_command *svd = &cmd->svd;
6488 if (!lex_force_match (s->lexer, T_LPAREN))
6490 svd->expr = matrix_parse_expr (s);
6492 || !lex_force_match (s->lexer, T_COMMA)
6493 || !matrix_parse_dst_var (s, &svd->u)
6494 || !lex_force_match (s->lexer, T_COMMA)
6495 || !matrix_parse_dst_var (s, &svd->s)
6496 || !lex_force_match (s->lexer, T_COMMA)
6497 || !matrix_parse_dst_var (s, &svd->v)
6498 || !lex_force_match (s->lexer, T_RPAREN))
6504 matrix_cmd_destroy (cmd);
6509 matrix_cmd_execute_svd (struct svd_command *svd)
6511 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6515 if (m->size1 >= m->size2)
6518 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6519 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6520 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6521 gsl_vector *work = gsl_vector_alloc (A->size2);
6522 gsl_linalg_SV_decomp (A, V, &Sv, work);
6523 gsl_vector_free (work);
6525 matrix_var_set (svd->u, A);
6526 matrix_var_set (svd->s, S);
6527 matrix_var_set (svd->v, V);
6531 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6532 gsl_matrix_transpose_memcpy (At, m);
6533 gsl_matrix_free (m);
6535 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6536 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6537 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6538 gsl_vector *work = gsl_vector_alloc (At->size2);
6539 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6540 gsl_vector_free (work);
6542 matrix_var_set (svd->v, At);
6543 matrix_var_set (svd->s, St);
6544 matrix_var_set (svd->u, Vt);
6549 matrix_cmd_execute (struct matrix_cmd *cmd)
6554 matrix_cmd_execute_compute (&cmd->compute);
6558 matrix_cmd_execute_print (&cmd->print);
6562 return matrix_cmd_execute_do_if (&cmd->do_if);
6565 matrix_cmd_execute_loop (&cmd->loop);
6572 matrix_cmd_execute_display (&cmd->display);
6576 matrix_cmd_execute_release (&cmd->release);
6580 matrix_cmd_execute_save (&cmd->save);
6584 matrix_cmd_execute_read (&cmd->read);
6588 matrix_cmd_execute_write (&cmd->write);
6592 matrix_cmd_execute_get (&cmd->get);
6596 matrix_cmd_execute_msave (&cmd->msave);
6600 matrix_cmd_execute_mget (&cmd->mget);
6604 matrix_cmd_execute_eigen (&cmd->eigen);
6608 matrix_cmd_execute_setdiag (&cmd->setdiag);
6612 matrix_cmd_execute_svd (&cmd->svd);
6620 matrix_cmds_uninit (struct matrix_cmds *cmds)
6622 for (size_t i = 0; i < cmds->n; i++)
6623 matrix_cmd_destroy (cmds->commands[i]);
6624 free (cmds->commands);
6628 matrix_cmd_destroy (struct matrix_cmd *cmd)
6636 matrix_lvalue_destroy (cmd->compute.lvalue);
6637 matrix_expr_destroy (cmd->compute.rvalue);
6641 matrix_expr_destroy (cmd->print.expression);
6642 free (cmd->print.title);
6643 print_labels_destroy (cmd->print.rlabels);
6644 print_labels_destroy (cmd->print.clabels);
6648 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6650 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6651 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6653 free (cmd->do_if.clauses);
6657 matrix_expr_destroy (cmd->loop.start);
6658 matrix_expr_destroy (cmd->loop.end);
6659 matrix_expr_destroy (cmd->loop.increment);
6660 matrix_expr_destroy (cmd->loop.top_condition);
6661 matrix_expr_destroy (cmd->loop.bottom_condition);
6662 matrix_cmds_uninit (&cmd->loop.commands);
6672 free (cmd->release.vars);
6676 matrix_expr_destroy (cmd->save.expression);
6677 fh_unref (cmd->save.outfile);
6678 string_array_destroy (cmd->save.variables);
6679 matrix_expr_destroy (cmd->save.names);
6680 stringi_set_destroy (&cmd->save.strings);
6684 matrix_lvalue_destroy (cmd->read.dst);
6685 matrix_expr_destroy (cmd->read.size);
6689 matrix_expr_destroy (cmd->write.expression);
6690 free (cmd->write.format);
6694 matrix_lvalue_destroy (cmd->get.dst);
6695 fh_unref (cmd->get.file);
6696 free (cmd->get.encoding);
6697 string_array_destroy (&cmd->get.variables);
6701 free (cmd->msave.varname_);
6702 matrix_expr_destroy (cmd->msave.expr);
6703 matrix_expr_destroy (cmd->msave.factors);
6704 matrix_expr_destroy (cmd->msave.splits);
6708 fh_unref (cmd->mget.file);
6709 stringi_set_destroy (&cmd->mget.rowtypes);
6713 matrix_expr_destroy (cmd->eigen.expr);
6717 matrix_expr_destroy (cmd->setdiag.expr);
6721 matrix_expr_destroy (cmd->svd.expr);
6727 struct matrix_command_name
6730 struct matrix_cmd *(*parse) (struct matrix_state *);
6733 static const struct matrix_command_name *
6734 matrix_parse_command_name (struct lexer *lexer)
6736 static const struct matrix_command_name commands[] = {
6737 { "COMPUTE", matrix_parse_compute },
6738 { "CALL EIGEN", matrix_parse_eigen },
6739 { "CALL SETDIAG", matrix_parse_setdiag },
6740 { "CALL SVD", matrix_parse_svd },
6741 { "PRINT", matrix_parse_print },
6742 { "DO IF", matrix_parse_do_if },
6743 { "LOOP", matrix_parse_loop },
6744 { "BREAK", matrix_parse_break },
6745 { "READ", matrix_parse_read },
6746 { "WRITE", matrix_parse_write },
6747 { "GET", matrix_parse_get },
6748 { "SAVE", matrix_parse_save },
6749 { "MGET", matrix_parse_mget },
6750 { "MSAVE", matrix_parse_msave },
6751 { "DISPLAY", matrix_parse_display },
6752 { "RELEASE", matrix_parse_release },
6754 static size_t n = sizeof commands / sizeof *commands;
6756 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6757 if (lex_match_phrase (lexer, c->name))
6762 static struct matrix_cmd *
6763 matrix_parse_command (struct matrix_state *s)
6765 size_t nesting_level = SIZE_MAX;
6767 struct matrix_cmd *c = NULL;
6768 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6770 lex_error (s->lexer, _("Unknown matrix command."));
6771 else if (!cmd->parse)
6772 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6776 nesting_level = output_open_group (
6777 group_item_create_nocopy (utf8_to_title (cmd->name),
6778 utf8_to_title (cmd->name)));
6783 lex_end_of_command (s->lexer);
6784 lex_discard_rest_of_command (s->lexer);
6785 if (nesting_level != SIZE_MAX)
6786 output_close_groups (nesting_level);
6792 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6794 if (!lex_force_match (lexer, T_ENDCMD))
6797 struct matrix_state state = {
6799 .session = dataset_session (ds),
6801 .vars = HMAP_INITIALIZER (state.vars),
6806 while (lex_match (lexer, T_ENDCMD))
6808 if (lex_token (lexer) == T_STOP)
6810 msg (SE, _("Unexpected end of input expecting matrix command."));
6814 if (lex_match_phrase (lexer, "END MATRIX"))
6817 struct matrix_cmd *c = matrix_parse_command (&state);
6820 matrix_cmd_execute (c);
6821 matrix_cmd_destroy (c);
6825 struct matrix_var *var, *next;
6826 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6829 gsl_matrix_free (var->value);
6830 hmap_delete (&state.vars, &var->hmap_node);
6833 hmap_destroy (&state.vars);
6834 fh_unref (state.prev_save_outfile);
6837 dict_unref (state.common->dict);
6838 casewriter_destroy (state.common->writer);
6839 free (state.common);
6841 fh_unref (state.prev_read_file);
6842 for (size_t i = 0; i < state.n_read_files; i++)
6843 read_file_destroy (state.read_files[i]);
6844 free (state.read_files);
6845 fh_unref (state.prev_write_file);
6846 for (size_t i = 0; i < state.n_write_files; i++)
6847 write_file_destroy (state.write_files[i]);
6848 free (state.write_files);