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 var_syntax *vars;
3440 struct matrix_var *names;
3442 /* Treatment of missing values. */
3447 MGET_ERROR, /* Flag error on command. */
3448 MGET_ACCEPT, /* Accept without change, user-missing only. */
3449 MGET_OMIT, /* Drop this case. */
3450 MGET_RECODE /* Recode to 'substitute'. */
3453 double substitute; /* MGET_RECODE only. */
3459 struct msave_command
3461 struct msave_common *common;
3463 struct matrix_expr *expr;
3464 const char *rowtype;
3465 struct matrix_expr *factors;
3466 struct matrix_expr *splits;
3472 struct matrix_state *state;
3473 struct file_handle *file;
3474 struct stringi_set rowtypes;
3478 struct eigen_command
3480 struct matrix_expr *expr;
3481 struct matrix_var *evec;
3482 struct matrix_var *eval;
3486 struct setdiag_command
3488 struct matrix_var *dst;
3489 struct matrix_expr *expr;
3495 struct matrix_expr *expr;
3496 struct matrix_var *u;
3497 struct matrix_var *s;
3498 struct matrix_var *v;
3504 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
3505 static bool matrix_cmd_execute (struct matrix_cmd *);
3506 static void matrix_cmd_destroy (struct matrix_cmd *);
3509 static struct matrix_cmd *
3510 matrix_parse_compute (struct matrix_state *s)
3512 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3513 *cmd = (struct matrix_cmd) {
3514 .type = MCMD_COMPUTE,
3515 .compute = { .lvalue = NULL },
3518 struct compute_command *compute = &cmd->compute;
3519 compute->lvalue = matrix_lvalue_parse (s);
3520 if (!compute->lvalue)
3523 if (!lex_force_match (s->lexer, T_EQUALS))
3526 compute->rvalue = matrix_parse_expr (s);
3527 if (!compute->rvalue)
3533 matrix_cmd_destroy (cmd);
3538 matrix_cmd_execute_compute (struct compute_command *compute)
3540 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
3544 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
3548 print_labels_destroy (struct print_labels *labels)
3552 string_array_destroy (&labels->literals);
3553 matrix_expr_destroy (labels->expr);
3559 parse_literal_print_labels (struct matrix_state *s,
3560 struct print_labels **labelsp)
3562 lex_match (s->lexer, T_EQUALS);
3563 print_labels_destroy (*labelsp);
3564 *labelsp = xzalloc (sizeof **labelsp);
3565 while (lex_token (s->lexer) != T_SLASH
3566 && lex_token (s->lexer) != T_ENDCMD
3567 && lex_token (s->lexer) != T_STOP)
3569 struct string label = DS_EMPTY_INITIALIZER;
3570 while (lex_token (s->lexer) != T_COMMA
3571 && lex_token (s->lexer) != T_SLASH
3572 && lex_token (s->lexer) != T_ENDCMD
3573 && lex_token (s->lexer) != T_STOP)
3575 if (!ds_is_empty (&label))
3576 ds_put_byte (&label, ' ');
3578 if (lex_token (s->lexer) == T_STRING)
3579 ds_put_cstr (&label, lex_tokcstr (s->lexer));
3582 char *rep = lex_next_representation (s->lexer, 0, 0);
3583 ds_put_cstr (&label, rep);
3588 string_array_append_nocopy (&(*labelsp)->literals,
3589 ds_steal_cstr (&label));
3591 if (!lex_match (s->lexer, T_COMMA))
3597 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
3599 lex_match (s->lexer, T_EQUALS);
3600 struct matrix_expr *e = matrix_parse_exp (s);
3604 print_labels_destroy (*labelsp);
3605 *labelsp = xzalloc (sizeof **labelsp);
3606 (*labelsp)->expr = e;
3610 static struct matrix_cmd *
3611 matrix_parse_print (struct matrix_state *s)
3613 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
3614 *cmd = (struct matrix_cmd) {
3617 .use_default_format = true,
3621 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
3624 for (size_t i = 0; ; i++)
3626 enum token_type t = lex_next_token (s->lexer, i);
3627 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
3629 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
3631 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
3634 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
3639 cmd->print.expression = matrix_parse_exp (s);
3640 if (!cmd->print.expression)
3644 while (lex_match (s->lexer, T_SLASH))
3646 if (lex_match_id (s->lexer, "FORMAT"))
3648 lex_match (s->lexer, T_EQUALS);
3649 if (!parse_format_specifier (s->lexer, &cmd->print.format))
3651 cmd->print.use_default_format = false;
3653 else if (lex_match_id (s->lexer, "TITLE"))
3655 lex_match (s->lexer, T_EQUALS);
3656 if (!lex_force_string (s->lexer))
3658 free (cmd->print.title);
3659 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
3662 else if (lex_match_id (s->lexer, "SPACE"))
3664 lex_match (s->lexer, T_EQUALS);
3665 if (lex_match_id (s->lexer, "NEWPAGE"))
3666 cmd->print.space = -1;
3667 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
3669 cmd->print.space = lex_integer (s->lexer);
3675 else if (lex_match_id (s->lexer, "RLABELS"))
3676 parse_literal_print_labels (s, &cmd->print.rlabels);
3677 else if (lex_match_id (s->lexer, "CLABELS"))
3678 parse_literal_print_labels (s, &cmd->print.clabels);
3679 else if (lex_match_id (s->lexer, "RNAMES"))
3681 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
3684 else if (lex_match_id (s->lexer, "CNAMES"))
3686 if (!parse_expr_print_labels (s, &cmd->print.clabels))
3691 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
3692 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
3700 matrix_cmd_destroy (cmd);
3705 matrix_is_integer (const gsl_matrix *m)
3707 for (size_t y = 0; y < m->size1; y++)
3708 for (size_t x = 0; x < m->size2; x++)
3710 double d = gsl_matrix_get (m, y, x);
3718 matrix_max_magnitude (const gsl_matrix *m)
3721 for (size_t y = 0; y < m->size1; y++)
3722 for (size_t x = 0; x < m->size2; x++)
3724 double d = fabs (gsl_matrix_get (m, y, x));
3732 format_fits (struct fmt_spec format, double x)
3734 char *s = data_out (&(union value) { .f = x }, NULL,
3735 &format, settings_get_fmt_settings ());
3736 bool fits = *s != '*' && !strchr (s, 'E');
3741 static struct fmt_spec
3742 default_format (const gsl_matrix *m, int *log_scale)
3746 double max = matrix_max_magnitude (m);
3748 if (matrix_is_integer (m))
3749 for (int w = 1; w <= 12; w++)
3751 struct fmt_spec format = { .type = FMT_F, .w = w };
3752 if (format_fits (format, -max))
3756 if (max >= 1e9 || max <= 1e-4)
3759 snprintf (s, sizeof s, "%.1e", max);
3761 const char *e = strchr (s, 'e');
3763 *log_scale = atoi (e + 1);
3766 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
3770 trimmed_string (double d)
3772 struct substring s = ss_buffer ((char *) &d, sizeof d);
3773 ss_rtrim (&s, ss_cstr (" "));
3774 return ss_xstrdup (s);
3777 static struct string_array *
3778 print_labels_get (const struct print_labels *labels, size_t n,
3779 const char *prefix, bool truncate)
3784 struct string_array *out = xzalloc (sizeof *out);
3785 if (labels->literals.n)
3786 string_array_clone (out, &labels->literals);
3787 else if (labels->expr)
3789 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
3790 if (m && is_vector (m))
3792 gsl_vector v = to_vector (m);
3793 for (size_t i = 0; i < v.size; i++)
3794 string_array_append_nocopy (out, trimmed_string (
3795 gsl_vector_get (&v, i)));
3797 gsl_matrix_free (m);
3803 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
3806 string_array_append (out, "");
3809 string_array_delete (out, out->n - 1);
3812 for (size_t i = 0; i < out->n; i++)
3814 char *s = out->strings[i];
3815 s[strnlen (s, 8)] = '\0';
3822 matrix_cmd_print_space (int space)
3825 output_item_submit (page_break_item_create ());
3826 for (int i = 0; i < space; i++)
3827 output_log ("%s", "");
3831 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
3832 struct fmt_spec format, int log_scale)
3834 matrix_cmd_print_space (print->space);
3836 output_log ("%s", print->title);
3838 output_log (" 10 ** %d X", log_scale);
3840 struct string_array *clabels = print_labels_get (print->clabels,
3841 m->size2, "col", true);
3842 if (clabels && format.w < 8)
3844 struct string_array *rlabels = print_labels_get (print->rlabels,
3845 m->size1, "row", true);
3849 struct string line = DS_EMPTY_INITIALIZER;
3851 ds_put_byte_multiple (&line, ' ', 8);
3852 for (size_t x = 0; x < m->size2; x++)
3853 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
3854 output_log_nocopy (ds_steal_cstr (&line));
3857 double scale = pow (10.0, log_scale);
3858 bool numeric = fmt_is_numeric (format.type);
3859 for (size_t y = 0; y < m->size1; y++)
3861 struct string line = DS_EMPTY_INITIALIZER;
3863 ds_put_format (&line, "%-8s", rlabels->strings[y]);
3865 for (size_t x = 0; x < m->size2; x++)
3867 double f = gsl_matrix_get (m, y, x);
3869 ? data_out (&(union value) { .f = f / scale}, NULL,
3870 &format, settings_get_fmt_settings ())
3871 : trimmed_string (f));
3872 ds_put_format (&line, " %s", s);
3875 output_log_nocopy (ds_steal_cstr (&line));
3878 string_array_destroy (rlabels);
3880 string_array_destroy (clabels);
3885 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
3886 const struct print_labels *print_labels, size_t n,
3887 const char *name, const char *prefix)
3889 struct string_array *labels = print_labels_get (print_labels, n, prefix,
3891 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
3892 for (size_t i = 0; i < n; i++)
3893 pivot_category_create_leaf (
3895 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
3896 : pivot_value_new_integer (i + 1)));
3898 d->hide_all_labels = true;
3899 string_array_destroy (labels);
3904 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
3905 struct fmt_spec format, int log_scale)
3907 struct pivot_table *table = pivot_table_create__ (
3908 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
3911 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
3913 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
3914 N_("Columns"), "col");
3916 struct pivot_footnote *footnote = NULL;
3919 char *s = xasprintf ("× 10**%d", log_scale);
3920 footnote = pivot_table_create_footnote (
3921 table, pivot_value_new_user_text_nocopy (s));
3924 double scale = pow (10.0, log_scale);
3925 bool numeric = fmt_is_numeric (format.type);
3926 for (size_t y = 0; y < m->size1; y++)
3927 for (size_t x = 0; x < m->size2; x++)
3929 double f = gsl_matrix_get (m, y, x);
3930 struct pivot_value *v;
3933 v = pivot_value_new_number (f / scale);
3934 v->numeric.format = format;
3937 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
3939 pivot_value_add_footnote (v, footnote);
3940 pivot_table_put2 (table, y, x, v);
3943 pivot_table_submit (table);
3947 matrix_cmd_execute_print (const struct print_command *print)
3949 if (print->expression)
3951 gsl_matrix *m = matrix_expr_evaluate (print->expression);
3956 struct fmt_spec format = (print->use_default_format
3957 ? default_format (m, &log_scale)
3960 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
3961 matrix_cmd_print_text (print, m, format, log_scale);
3963 matrix_cmd_print_tables (print, m, format, log_scale);
3965 gsl_matrix_free (m);
3969 matrix_cmd_print_space (print->space);
3971 output_log ("%s", print->title);
3978 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
3979 const char *command_name,
3980 const char *stop1, const char *stop2)
3982 lex_end_of_command (s->lexer);
3983 lex_discard_rest_of_command (s->lexer);
3985 size_t allocated = 0;
3988 while (lex_token (s->lexer) == T_ENDCMD)
3991 if (lex_at_phrase (s->lexer, stop1)
3992 || (stop2 && lex_at_phrase (s->lexer, stop2)))
3995 if (lex_at_phrase (s->lexer, "END MATRIX"))
3997 msg (SE, _("Premature END MATRIX within %s."), command_name);
4001 struct matrix_cmd *cmd = matrix_parse_command (s);
4005 if (c->n >= allocated)
4006 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
4007 c->commands[c->n++] = cmd;
4012 matrix_parse_do_if_clause (struct matrix_state *s,
4013 struct do_if_command *ifc,
4014 bool parse_condition,
4015 size_t *allocated_clauses)
4017 if (ifc->n_clauses >= *allocated_clauses)
4018 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
4019 sizeof *ifc->clauses);
4020 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
4021 *c = (struct do_if_clause) { .condition = NULL };
4023 if (parse_condition)
4025 c->condition = matrix_parse_expr (s);
4030 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
4033 static struct matrix_cmd *
4034 matrix_parse_do_if (struct matrix_state *s)
4036 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4037 *cmd = (struct matrix_cmd) {
4039 .do_if = { .n_clauses = 0 }
4042 size_t allocated_clauses = 0;
4045 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
4048 while (lex_match_phrase (s->lexer, "ELSE IF"));
4050 if (lex_match_id (s->lexer, "ELSE")
4051 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
4054 if (!lex_match_phrase (s->lexer, "END IF"))
4059 matrix_cmd_destroy (cmd);
4064 matrix_cmd_execute_do_if (struct do_if_command *cmd)
4066 for (size_t i = 0; i < cmd->n_clauses; i++)
4068 struct do_if_clause *c = &cmd->clauses[i];
4072 if (!matrix_expr_evaluate_scalar (c->condition,
4073 i ? "ELSE IF" : "DO IF",
4078 for (size_t j = 0; j < c->commands.n; j++)
4079 if (!matrix_cmd_execute (c->commands.commands[j]))
4086 static struct matrix_cmd *
4087 matrix_parse_loop (struct matrix_state *s)
4089 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4090 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
4092 struct loop_command *loop = &cmd->loop;
4093 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
4095 struct substring name = lex_tokss (s->lexer);
4096 loop->var = matrix_var_lookup (s, name);
4098 loop->var = matrix_var_create (s, name);
4103 loop->start = matrix_parse_expr (s);
4104 if (!loop->start || !lex_force_match (s->lexer, T_TO))
4107 loop->end = matrix_parse_expr (s);
4111 if (lex_match (s->lexer, T_BY))
4113 loop->increment = matrix_parse_expr (s);
4114 if (!loop->increment)
4119 if (lex_match_id (s->lexer, "IF"))
4121 loop->top_condition = matrix_parse_expr (s);
4122 if (!loop->top_condition)
4126 bool was_in_loop = s->in_loop;
4128 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
4130 s->in_loop = was_in_loop;
4134 if (!lex_match_phrase (s->lexer, "END LOOP"))
4137 if (lex_match_id (s->lexer, "IF"))
4139 loop->bottom_condition = matrix_parse_expr (s);
4140 if (!loop->bottom_condition)
4147 matrix_cmd_destroy (cmd);
4152 matrix_cmd_execute_loop (struct loop_command *cmd)
4156 long int increment = 1;
4159 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
4160 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
4162 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
4166 if (increment > 0 ? value > end
4167 : increment < 0 ? value < end
4172 int mxloops = settings_get_mxloops ();
4173 for (int i = 0; i < mxloops; i++)
4178 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
4180 gsl_matrix_free (cmd->var->value);
4181 cmd->var->value = NULL;
4183 if (!cmd->var->value)
4184 cmd->var->value = gsl_matrix_alloc (1, 1);
4186 gsl_matrix_set (cmd->var->value, 0, 0, value);
4189 if (cmd->top_condition)
4192 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
4197 for (size_t j = 0; j < cmd->commands.n; j++)
4198 if (!matrix_cmd_execute (cmd->commands.commands[j]))
4201 if (cmd->bottom_condition)
4204 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
4205 "END LOOP IF", &d) || d > 0)
4212 if (increment > 0 ? value > end : value < end)
4218 static struct matrix_cmd *
4219 matrix_parse_break (struct matrix_state *s)
4223 msg (SE, _("BREAK not inside LOOP."));
4227 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4228 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
4232 static struct matrix_cmd *
4233 matrix_parse_display (struct matrix_state *s)
4235 bool dictionary = false;
4236 bool status = false;
4237 while (lex_token (s->lexer) == T_ID)
4239 if (lex_match_id (s->lexer, "DICTIONARY"))
4241 else if (lex_match_id (s->lexer, "STATUS"))
4245 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
4249 if (!dictionary && !status)
4250 dictionary = status = true;
4252 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4253 *cmd = (struct matrix_cmd) {
4254 .type = MCMD_DISPLAY,
4257 .dictionary = dictionary,
4265 compare_matrix_var_pointers (const void *a_, const void *b_)
4267 const struct matrix_var *const *ap = a_;
4268 const struct matrix_var *const *bp = b_;
4269 const struct matrix_var *a = *ap;
4270 const struct matrix_var *b = *bp;
4271 return strcmp (a->name, b->name);
4275 matrix_cmd_execute_display (struct display_command *cmd)
4277 const struct matrix_state *s = cmd->state;
4279 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
4280 pivot_dimension_create (
4281 table, PIVOT_AXIS_COLUMN, N_("Property"),
4282 N_("Rows"), N_("Columns"), N_("Size (kB)"));
4284 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
4287 const struct matrix_var *var;
4288 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
4290 vars[n_vars++] = var;
4291 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
4293 struct pivot_dimension *rows = pivot_dimension_create (
4294 table, PIVOT_AXIS_ROW, N_("Variable"));
4295 for (size_t i = 0; i < n_vars; i++)
4297 const struct matrix_var *var = vars[i];
4298 pivot_category_create_leaf (
4299 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
4301 size_t r = var->value->size1;
4302 size_t c = var->value->size2;
4303 double values[] = { r, c, r * c * sizeof (double) / 1024 };
4304 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
4305 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
4307 pivot_table_submit (table);
4310 static struct matrix_cmd *
4311 matrix_parse_release (struct matrix_state *s)
4313 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4314 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
4316 size_t allocated_vars = 0;
4317 while (lex_token (s->lexer) == T_ID)
4319 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
4322 if (cmd->release.n_vars >= allocated_vars)
4323 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
4324 sizeof *cmd->release.vars);
4325 cmd->release.vars[cmd->release.n_vars++] = var;
4328 lex_error (s->lexer, _("Variable name expected."));
4331 if (!lex_match (s->lexer, T_COMMA))
4339 matrix_cmd_execute_release (struct release_command *cmd)
4341 for (size_t i = 0; i < cmd->n_vars; i++)
4343 struct matrix_var *var = cmd->vars[i];
4344 gsl_matrix_free (var->value);
4349 static struct matrix_cmd *
4350 matrix_parse_save (struct matrix_state *s)
4352 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4353 *cmd = (struct matrix_cmd) {
4356 .strings = STRINGI_SET_INITIALIZER (cmd->save.strings)
4360 struct save_command *save = &cmd->save;
4361 save->expression = matrix_parse_exp (s);
4362 if (!save->expression)
4365 while (lex_match (s->lexer, T_SLASH))
4367 if (lex_match_id (s->lexer, "OUTFILE"))
4369 lex_match (s->lexer, T_EQUALS);
4371 fh_unref (save->outfile);
4372 save->outfile = (lex_match (s->lexer, T_ASTERISK)
4374 : fh_parse (s->lexer, FH_REF_FILE, s->session));
4378 else if (lex_match_id (s->lexer, "VARIABLES"))
4380 lex_match (s->lexer, T_EQUALS);
4384 struct dictionary *d = dict_create (get_default_encoding ());
4385 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
4386 PV_NO_SCRATCH | PV_NO_DUPLICATE);
4391 string_array_destroy (save->variables);
4392 if (!save->variables)
4393 save->variables = xmalloc (sizeof *save->variables);
4394 *save->variables = (struct string_array) {
4400 else if (lex_match_id (s->lexer, "NAMES"))
4402 lex_match (s->lexer, T_EQUALS);
4403 matrix_expr_destroy (save->names);
4404 save->names = matrix_parse_exp (s);
4408 else if (lex_match_id (s->lexer, "STRINGS"))
4410 lex_match (s->lexer, T_EQUALS);
4411 while (lex_token (s->lexer) == T_ID)
4413 stringi_set_insert (&save->strings, lex_tokcstr (s->lexer));
4415 if (!lex_match (s->lexer, T_COMMA))
4421 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
4429 if (s->prev_save_outfile)
4430 save->outfile = fh_ref (s->prev_save_outfile);
4433 lex_sbc_missing ("OUTFILE");
4437 fh_unref (s->prev_save_outfile);
4438 s->prev_save_outfile = fh_ref (save->outfile);
4440 if (save->variables && save->names)
4442 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
4443 matrix_expr_destroy (save->names);
4450 matrix_cmd_destroy (cmd);
4455 matrix_cmd_execute_save (const struct save_command *save)
4457 assert (save->outfile != fh_inline_file ()); /* XXX not yet implemented */
4459 gsl_matrix *m = matrix_expr_evaluate (save->expression);
4464 struct dictionary *dict = dict_create (get_default_encoding ());
4465 struct string_array names = { .n = 0 };
4466 if (save->variables)
4467 string_array_clone (&names, save->variables);
4468 else if (save->names)
4470 gsl_matrix *nm = matrix_expr_evaluate (save->names);
4471 if (nm && is_vector (nm))
4473 gsl_vector nv = to_vector (nm);
4474 for (size_t i = 0; i < nv.size; i++)
4476 char *name = trimmed_string (gsl_vector_get (&nv, i));
4477 if (dict_id_is_valid (dict, name, true))
4478 string_array_append_nocopy (&names, name);
4485 struct stringi_set strings;
4486 stringi_set_clone (&strings, &save->strings);
4488 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
4493 name = names.strings[i];
4496 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
4500 int width = stringi_set_delete (&strings, name) ? 8 : 0;
4501 struct variable *var = dict_create_var (dict, name, width);
4504 msg (SE, _("Duplicate variable name %s in SAVE statement."),
4510 if (!stringi_set_is_empty (&strings))
4512 const char *example = stringi_set_node_get_string (
4513 stringi_set_first (&strings));
4514 msg (SE, ngettext ("STRINGS specified a variable %s, but no variable "
4515 "with that name was found on SAVE.",
4516 "STRINGS specified %2$zu variables, including %1$s, "
4517 "whose names were not found on SAVE.",
4518 stringi_set_count (&strings)),
4519 example, stringi_set_count (&strings));
4522 stringi_set_destroy (&strings);
4523 string_array_destroy (&names);
4527 gsl_matrix_free (m);
4532 struct casewriter *writer = any_writer_open (save->outfile, dict);
4535 gsl_matrix_free (m);
4540 const struct caseproto *proto = dict_get_proto (dict);
4541 for (size_t y = 0; y < m->size1; y++)
4543 struct ccase *c = case_create (proto);
4544 for (size_t x = 0; x < m->size2; x++)
4546 double d = gsl_matrix_get (m, y, x);
4547 int width = caseproto_get_width (proto, x);
4548 union value *value = case_data_rw_idx (c, x);
4552 memcpy (value->s, &d, width);
4554 casewriter_write (writer, c);
4556 casewriter_destroy (writer);
4558 gsl_matrix_free (m);
4562 static struct matrix_cmd *
4563 matrix_parse_read (struct matrix_state *s)
4565 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4566 *cmd = (struct matrix_cmd) {
4568 .read = { .format = FMT_F },
4571 struct file_handle *fh = NULL;
4572 char *encoding = NULL;
4573 struct read_command *read = &cmd->read;
4574 read->dst = matrix_lvalue_parse (s);
4579 int repetitions = 0;
4580 int record_width = 0;
4581 bool seen_format = false;
4582 while (lex_match (s->lexer, T_SLASH))
4584 if (lex_match_id (s->lexer, "FILE"))
4586 lex_match (s->lexer, T_EQUALS);
4589 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
4593 else if (lex_match_id (s->lexer, "ENCODING"))
4595 lex_match (s->lexer, T_EQUALS);
4596 if (!lex_force_string (s->lexer))
4600 encoding = ss_xstrdup (lex_tokss (s->lexer));
4604 else if (lex_match_id (s->lexer, "FIELD"))
4606 lex_match (s->lexer, T_EQUALS);
4608 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
4610 read->c1 = lex_integer (s->lexer);
4612 if (!lex_force_match (s->lexer, T_TO)
4613 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
4615 read->c2 = lex_integer (s->lexer) + 1;
4618 record_width = read->c2 - read->c1;
4619 if (lex_match (s->lexer, T_BY))
4621 if (!lex_force_int_range (s->lexer, "BY", 1,
4622 read->c2 - read->c1))
4624 by = lex_integer (s->lexer);
4627 if (record_width % by)
4629 msg (SE, _("BY %d does not evenly divide record width %d."),
4637 else if (lex_match_id (s->lexer, "SIZE"))
4639 lex_match (s->lexer, T_EQUALS);
4640 matrix_expr_destroy (read->size);
4641 read->size = matrix_parse_exp (s);
4645 else if (lex_match_id (s->lexer, "MODE"))
4647 lex_match (s->lexer, T_EQUALS);
4648 if (lex_match_id (s->lexer, "RECTANGULAR"))
4649 read->symmetric = false;
4650 else if (lex_match_id (s->lexer, "SYMMETRIC"))
4651 read->symmetric = true;
4654 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
4658 else if (lex_match_id (s->lexer, "REREAD"))
4659 read->reread = true;
4660 else if (lex_match_id (s->lexer, "FORMAT"))
4664 lex_sbc_only_once ("FORMAT");
4669 lex_match (s->lexer, T_EQUALS);
4671 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
4674 const char *p = lex_tokcstr (s->lexer);
4675 if (c_isdigit (p[0]))
4677 repetitions = atoi (p);
4678 p += strspn (p, "0123456789");
4679 if (!fmt_from_name (p, &read->format))
4681 lex_error (s->lexer, _("Unknown format %s."), p);
4686 else if (fmt_from_name (p, &read->format))
4690 struct fmt_spec format;
4691 if (!parse_format_specifier (s->lexer, &format))
4693 read->format = format.type;
4699 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
4700 "REREAD", "FORMAT");
4707 lex_sbc_missing ("FIELD");
4711 if (!read->dst->n_indexes && !read->size)
4713 msg (SE, _("SIZE is required for reading data into a full matrix "
4714 "(as opposed to a submatrix)."));
4720 if (s->prev_read_file)
4721 fh = fh_ref (s->prev_read_file);
4724 lex_sbc_missing ("FILE");
4728 fh_unref (s->prev_read_file);
4729 s->prev_read_file = fh_ref (fh);
4731 read->rf = read_file_create (s, fh);
4735 free (read->rf->encoding);
4736 read->rf->encoding = encoding;
4740 /* Field width may be specified in multiple ways:
4743 2. The format on FORMAT.
4744 3. The repetition factor on FORMAT.
4746 (2) and (3) are mutually exclusive.
4748 If more than one of these is present, they must agree. If none of them is
4749 present, then free-field format is used.
4751 if (repetitions > record_width)
4753 msg (SE, _("%d repetitions cannot fit in record width %d."),
4754 repetitions, record_width);
4757 int w = (repetitions ? record_width / repetitions
4763 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
4764 "which implies field width %d, "
4765 "but BY specifies field width %d."),
4766 repetitions, record_width, w, by);
4768 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
4777 matrix_cmd_destroy (cmd);
4783 parse_error (const struct dfm_reader *reader, enum fmt_type format,
4784 struct substring data, size_t y, size_t x,
4785 int first_column, int last_column, char *error)
4787 int line_number = dfm_get_line_number (reader);
4788 struct msg_location *location = xmalloc (sizeof *location);
4789 *location = (struct msg_location) {
4790 .file_name = xstrdup (dfm_get_file_name (reader)),
4791 .first_line = line_number,
4792 .last_line = line_number + 1,
4793 .first_column = first_column,
4794 .last_column = last_column,
4796 struct msg *m = xmalloc (sizeof *m);
4798 .category = MSG_C_DATA,
4799 .severity = MSG_S_WARNING,
4800 .location = location,
4801 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
4802 "for matrix row %zu, column %zu: %s"),
4803 (int) data.length, data.string, fmt_name (format),
4804 y + 1, x + 1, error),
4812 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
4813 gsl_matrix *m, struct substring p, size_t y, size_t x,
4814 const char *line_start)
4816 const char *input_encoding = dfm_reader_get_encoding (reader);
4819 if (fmt_is_numeric (read->format))
4822 error = data_in (p, input_encoding, read->format,
4823 settings_get_fmt_settings (), &v, 0, NULL);
4824 if (!error && v.f == SYSMIS)
4825 error = xstrdup (_("Matrix data may not contain missing value."));
4830 uint8_t s[sizeof (double)];
4831 union value v = { .s = s };
4832 error = data_in (p, input_encoding, read->format,
4833 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
4834 memcpy (&f, s, sizeof f);
4839 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
4840 int c2 = c1 + ss_utf8_count_columns (p) - 1;
4841 parse_error (reader, read->format, p, y, x, c1, c2, error);
4845 gsl_matrix_set (m, y, x, f);
4846 if (read->symmetric && x != y)
4847 gsl_matrix_set (m, x, y, f);
4852 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
4853 struct substring *line, const char **startp)
4855 if (dfm_eof (reader))
4857 msg (SE, _("Unexpected end of file reading matrix data."));
4860 dfm_expand_tabs (reader);
4861 struct substring record = dfm_get_record (reader);
4862 /* XXX need to recode record into UTF-8 */
4863 *startp = record.string;
4864 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
4869 matrix_read (struct read_command *read, struct dfm_reader *reader,
4872 for (size_t y = 0; y < m->size1; y++)
4874 size_t nx = read->symmetric ? y + 1 : m->size2;
4876 struct substring line = ss_empty ();
4877 const char *line_start = line.string;
4878 for (size_t x = 0; x < nx; x++)
4885 ss_ltrim (&line, ss_cstr (" ,"));
4886 if (!ss_is_empty (line))
4888 if (!matrix_read_line (read, reader, &line, &line_start))
4890 dfm_forward_record (reader);
4893 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
4897 if (!matrix_read_line (read, reader, &line, &line_start))
4899 size_t fields_per_line = (read->c2 - read->c1) / read->w;
4900 int f = x % fields_per_line;
4901 if (f == fields_per_line - 1)
4902 dfm_forward_record (reader);
4904 p = ss_substr (line, read->w * f, read->w);
4907 matrix_read_set_field (read, reader, m, p, y, x, line_start);
4911 dfm_forward_record (reader);
4914 ss_ltrim (&line, ss_cstr (" ,"));
4915 if (!ss_is_empty (line))
4918 msg (SW, _("Trailing garbage on line \"%.*s\""),
4919 (int) line.length, line.string);
4926 matrix_cmd_execute_read (struct read_command *read)
4928 struct index_vector iv0, iv1;
4929 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
4932 size_t size[2] = { SIZE_MAX, SIZE_MAX };
4935 gsl_matrix *m = matrix_expr_evaluate (read->size);
4941 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4942 "not a %zu×%zu matrix."), m->size1, m->size2);
4943 gsl_matrix_free (m);
4949 gsl_vector v = to_vector (m);
4953 d[0] = gsl_vector_get (&v, 0);
4956 else if (v.size == 2)
4958 d[0] = gsl_vector_get (&v, 0);
4959 d[1] = gsl_vector_get (&v, 1);
4963 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
4964 "not a %zu×%zu matrix."),
4965 m->size1, m->size2),
4966 gsl_matrix_free (m);
4971 gsl_matrix_free (m);
4973 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
4975 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
4976 "are outside valid range."),
4987 if (read->dst->n_indexes)
4989 size_t submatrix_size[2];
4990 if (read->dst->n_indexes == 2)
4992 submatrix_size[0] = iv0.n;
4993 submatrix_size[1] = iv1.n;
4995 else if (read->dst->var->value->size1 == 1)
4997 submatrix_size[0] = 1;
4998 submatrix_size[1] = iv0.n;
5002 submatrix_size[0] = iv0.n;
5003 submatrix_size[1] = 1;
5008 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
5010 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
5011 "differ from submatrix dimensions %zu×%zu."),
5013 submatrix_size[0], submatrix_size[1]);
5021 size[0] = submatrix_size[0];
5022 size[1] = submatrix_size[1];
5026 struct dfm_reader *reader = read_file_open (read->rf);
5028 dfm_reread_record (reader, 1);
5030 if (read->symmetric && size[0] != size[1])
5032 msg (SE, _("Cannot read non-square %zu×%zu matrix "
5033 "using READ with MODE=SYMMETRIC."),
5039 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
5040 matrix_read (read, reader, tmp);
5041 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
5044 static struct matrix_cmd *
5045 matrix_parse_write (struct matrix_state *s)
5047 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5048 *cmd = (struct matrix_cmd) {
5052 struct file_handle *fh = NULL;
5053 char *encoding = NULL;
5054 struct write_command *write = &cmd->write;
5055 write->expression = matrix_parse_exp (s);
5056 if (!write->expression)
5060 int repetitions = 0;
5061 int record_width = 0;
5062 enum fmt_type format = FMT_F;
5063 bool has_format = false;
5064 while (lex_match (s->lexer, T_SLASH))
5066 if (lex_match_id (s->lexer, "OUTFILE"))
5068 lex_match (s->lexer, T_EQUALS);
5071 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5075 else if (lex_match_id (s->lexer, "ENCODING"))
5077 lex_match (s->lexer, T_EQUALS);
5078 if (!lex_force_string (s->lexer))
5082 encoding = ss_xstrdup (lex_tokss (s->lexer));
5086 else if (lex_match_id (s->lexer, "FIELD"))
5088 lex_match (s->lexer, T_EQUALS);
5090 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5092 write->c1 = lex_integer (s->lexer);
5094 if (!lex_force_match (s->lexer, T_TO)
5095 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
5097 write->c2 = lex_integer (s->lexer) + 1;
5100 record_width = write->c2 - write->c1;
5101 if (lex_match (s->lexer, T_BY))
5103 if (!lex_force_int_range (s->lexer, "BY", 1,
5104 write->c2 - write->c1))
5106 by = lex_integer (s->lexer);
5109 if (record_width % by)
5111 msg (SE, _("BY %d does not evenly divide record width %d."),
5119 else if (lex_match_id (s->lexer, "MODE"))
5121 lex_match (s->lexer, T_EQUALS);
5122 if (lex_match_id (s->lexer, "RECTANGULAR"))
5123 write->triangular = false;
5124 else if (lex_match_id (s->lexer, "TRIANGULAR"))
5125 write->triangular = true;
5128 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
5132 else if (lex_match_id (s->lexer, "HOLD"))
5134 else if (lex_match_id (s->lexer, "FORMAT"))
5136 if (has_format || write->format)
5138 lex_sbc_only_once ("FORMAT");
5142 lex_match (s->lexer, T_EQUALS);
5144 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
5147 const char *p = lex_tokcstr (s->lexer);
5148 if (c_isdigit (p[0]))
5150 repetitions = atoi (p);
5151 p += strspn (p, "0123456789");
5152 if (!fmt_from_name (p, &format))
5154 lex_error (s->lexer, _("Unknown format %s."), p);
5160 else if (fmt_from_name (p, &format))
5167 struct fmt_spec spec;
5168 if (!parse_format_specifier (s->lexer, &spec))
5170 write->format = xmemdup (&spec, sizeof spec);
5175 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
5183 lex_sbc_missing ("FIELD");
5189 if (s->prev_write_file)
5190 fh = fh_ref (s->prev_write_file);
5193 lex_sbc_missing ("OUTFILE");
5197 fh_unref (s->prev_write_file);
5198 s->prev_write_file = fh_ref (fh);
5200 write->wf = write_file_create (s, fh);
5204 free (write->wf->encoding);
5205 write->wf->encoding = encoding;
5209 /* Field width may be specified in multiple ways:
5212 2. The format on FORMAT.
5213 3. The repetition factor on FORMAT.
5215 (2) and (3) are mutually exclusive.
5217 If more than one of these is present, they must agree. If none of them is
5218 present, then free-field format is used.
5220 if (repetitions > record_width)
5222 msg (SE, _("%d repetitions cannot fit in record width %d."),
5223 repetitions, record_width);
5226 int w = (repetitions ? record_width / repetitions
5227 : write->format ? write->format->w
5232 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
5233 "which implies field width %d, "
5234 "but BY specifies field width %d."),
5235 repetitions, record_width, w, by);
5237 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
5241 if (w && !write->format)
5243 write->format = xmalloc (sizeof *write->format);
5244 *write->format = (struct fmt_spec) { .type = format, .w = w };
5246 if (!fmt_check_output (write->format))
5250 if (write->format && fmt_var_width (write->format) > sizeof (double))
5252 char s[FMT_STRING_LEN_MAX + 1];
5253 fmt_to_string (write->format, s);
5254 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
5255 s, sizeof (double));
5263 matrix_cmd_destroy (cmd);
5268 matrix_cmd_execute_write (struct write_command *write)
5270 gsl_matrix *m = matrix_expr_evaluate (write->expression);
5274 if (write->triangular && m->size1 != m->size2)
5276 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
5277 "the matrix to be written has dimensions %zu×%zu."),
5278 m->size1, m->size2);
5279 gsl_matrix_free (m);
5283 struct dfm_writer *writer = write_file_open (write->wf);
5284 if (!writer || !m->size1)
5286 gsl_matrix_free (m);
5290 const struct fmt_settings *settings = settings_get_fmt_settings ();
5291 struct u8_line *line = write->wf->held;
5292 for (size_t y = 0; y < m->size1; y++)
5296 line = xmalloc (sizeof *line);
5297 u8_line_init (line);
5299 size_t nx = write->triangular ? y + 1 : m->size2;
5301 for (size_t x = 0; x < nx; x++)
5304 double f = gsl_matrix_get (m, y, x);
5308 if (fmt_is_numeric (write->format->type))
5311 v.s = (uint8_t *) &f;
5312 s = data_out (&v, NULL, write->format, settings);
5316 s = xmalloc (DBL_BUFSIZE_BOUND);
5317 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
5318 >= DBL_BUFSIZE_BOUND)
5321 size_t len = strlen (s);
5322 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
5323 if (width + x0 > write->c2)
5325 dfm_put_record_utf8 (writer, line->s.ss.string,
5327 u8_line_clear (line);
5330 u8_line_put (line, x0, x0 + width, s, len);
5333 x0 += write->format ? write->format->w : width + 1;
5336 if (y + 1 >= m->size1 && write->hold)
5338 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
5339 u8_line_clear (line);
5343 u8_line_destroy (line);
5347 write->wf->held = line;
5349 gsl_matrix_free (m);
5352 static struct matrix_cmd *
5353 matrix_parse_get (struct matrix_state *s)
5355 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5356 *cmd = (struct matrix_cmd) {
5359 .dataset = s->dataset,
5360 .user = { .treatment = MGET_ERROR },
5361 .system = { .treatment = MGET_ERROR },
5365 struct get_command *get = &cmd->get;
5366 get->dst = matrix_lvalue_parse (s);
5370 while (lex_match (s->lexer, T_SLASH))
5372 if (lex_match_id (s->lexer, "FILE"))
5374 lex_match (s->lexer, T_EQUALS);
5376 fh_unref (get->file);
5377 if (lex_match (s->lexer, T_ASTERISK))
5381 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
5386 else if (lex_match_id (s->lexer, "ENCODING"))
5388 lex_match (s->lexer, T_EQUALS);
5389 if (!lex_force_string (s->lexer))
5392 free (get->encoding);
5393 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
5397 else if (lex_match_id (s->lexer, "VARIABLES"))
5399 lex_match (s->lexer, T_EQUALS);
5403 lex_sbc_only_once ("VARIABLES");
5407 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
5410 else if (lex_match_id (s->lexer, "NAMES"))
5412 lex_match (s->lexer, T_EQUALS);
5413 if (!lex_force_id (s->lexer))
5416 struct substring name = lex_tokss (s->lexer);
5417 get->names = matrix_var_lookup (s, name);
5419 get->names = matrix_var_create (s, name);
5422 else if (lex_match_id (s->lexer, "MISSING"))
5424 lex_match (s->lexer, T_EQUALS);
5425 if (lex_match_id (s->lexer, "ACCEPT"))
5426 get->user.treatment = MGET_ACCEPT;
5427 else if (lex_match_id (s->lexer, "OMIT"))
5428 get->user.treatment = MGET_OMIT;
5429 else if (lex_is_number (s->lexer))
5431 get->user.treatment = MGET_RECODE;
5432 get->user.substitute = lex_number (s->lexer);
5437 lex_error (s->lexer, NULL);
5441 else if (lex_match_id (s->lexer, "SYSMIS"))
5443 lex_match (s->lexer, T_EQUALS);
5444 if (lex_match_id (s->lexer, "OMIT"))
5445 get->system.treatment = MGET_OMIT;
5446 else if (lex_is_number (s->lexer))
5448 get->system.treatment = MGET_RECODE;
5449 get->system.substitute = lex_number (s->lexer);
5454 lex_error (s->lexer, NULL);
5460 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
5461 "MISSING", "SYSMIS");
5466 if (get->user.treatment != MGET_ACCEPT)
5467 get->system.treatment = MGET_ERROR;
5472 matrix_cmd_destroy (cmd);
5477 matrix_cmd_execute_get__ (struct get_command *get,
5478 const struct dictionary *dict,
5479 struct casereader *reader)
5481 struct variable **vars;
5486 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
5487 &vars, &n_vars, PV_NUMERIC))
5492 n_vars = dict_get_var_cnt (dict);
5493 vars = xnmalloc (n_vars, sizeof *vars);
5494 for (size_t i = 0; i < n_vars; i++)
5496 struct variable *var = dict_get_var (dict, i);
5497 if (!var_is_numeric (var))
5499 msg (SE, _("GET: Variable %s is not numeric."),
5500 var_get_name (var));
5510 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
5511 for (size_t i = 0; i < n_vars; i++)
5513 char s[sizeof (double)];
5515 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
5516 memcpy (&f, s, sizeof f);
5517 gsl_matrix_set (names, i, 0, f);
5520 gsl_matrix_free (get->names->value);
5521 get->names->value = names;
5525 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
5526 long long int casenum = 1;
5528 for (struct ccase *c = casereader_read (reader); c;
5529 c = casereader_read (reader), casenum++)
5531 if (n_rows >= m->size1)
5533 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
5534 for (size_t y = 0; y < n_rows; y++)
5535 for (size_t x = 0; x < n_vars; x++)
5536 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
5537 gsl_matrix_free (m);
5542 for (size_t x = 0; x < n_vars; x++)
5544 const struct variable *var = vars[x];
5545 double d = case_num (c, var);
5548 if (get->system.treatment == MGET_RECODE)
5549 d = get->system.substitute;
5550 else if (get->system.treatment == MGET_OMIT)
5554 msg (SE, _("GET: Variable %s in case %lld "
5555 "is system-missing."),
5556 var_get_name (var), casenum);
5560 else if (var_is_num_missing (var, d, MV_USER))
5562 if (get->user.treatment == MGET_RECODE)
5563 d = get->user.substitute;
5564 else if (get->user.treatment == MGET_OMIT)
5566 else if (get->user.treatment != MGET_ACCEPT)
5568 msg (SE, _("GET: Variable %s in case %lld has user-missing "
5570 var_get_name (var), casenum, d);
5574 gsl_matrix_set (m, n_rows, x, d);
5585 matrix_lvalue_evaluate_and_assign (get->dst, m);
5588 gsl_matrix_free (m);
5593 matrix_cmd_execute_get (struct get_command *get)
5595 struct dictionary *dict;
5596 struct casereader *reader;
5599 reader = any_reader_open_and_decode (get->file, get->encoding,
5606 if (dict_get_var_cnt (dataset_dict (get->dataset)) == 0)
5608 msg (ME, _("GET cannot read empty active file."));
5611 reader = proc_open (get->dataset);
5612 dict = dict_ref (dataset_dict (get->dataset));
5615 matrix_cmd_execute_get__ (get, dict, reader);
5618 casereader_destroy (reader);
5620 proc_commit (get->dataset);
5624 match_rowtype (struct lexer *lexer)
5626 static const char *rowtypes[] = {
5627 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
5629 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
5631 for (size_t i = 0; i < n_rowtypes; i++)
5632 if (lex_match_id (lexer, rowtypes[i]))
5635 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
5640 parse_var_names (struct lexer *lexer, struct string_array *sa)
5642 lex_match (lexer, T_EQUALS);
5644 string_array_clear (sa);
5646 struct dictionary *dict = dict_create (get_default_encoding ());
5647 dict_create_var_assert (dict, "ROWTYPE_", 8);
5648 dict_create_var_assert (dict, "VARNAME_", 8);
5651 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
5652 PV_NO_DUPLICATE | PV_NO_SCRATCH);
5657 string_array_clear (sa);
5658 sa->strings = names;
5659 sa->n = sa->allocated = n_names;
5665 msave_common_uninit (struct msave_common *common)
5669 fh_unref (common->outfile);
5670 string_array_destroy (&common->variables);
5671 string_array_destroy (&common->fnames);
5672 string_array_destroy (&common->snames);
5677 compare_variables (const char *keyword,
5678 const struct string_array *new,
5679 const struct string_array *old)
5685 msg (SE, _("%s may only be specified on MSAVE if it was specified "
5686 "on the first MSAVE within MATRIX."), keyword);
5689 else if (!string_array_equal_case (old, new))
5691 msg (SE, _("%s must specify the same variables each time within "
5692 "a given MATRIX."), keyword);
5699 static struct matrix_cmd *
5700 matrix_parse_msave (struct matrix_state *s)
5702 struct msave_common common = { .outfile = NULL };
5703 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5704 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
5706 struct msave_command *msave = &cmd->msave;
5707 if (lex_token (s->lexer) == T_ID)
5708 msave->varname_ = ss_xstrdup (lex_tokss (s->lexer));
5709 msave->expr = matrix_parse_exp (s);
5713 while (lex_match (s->lexer, T_SLASH))
5715 if (lex_match_id (s->lexer, "TYPE"))
5717 lex_match (s->lexer, T_EQUALS);
5719 msave->rowtype = match_rowtype (s->lexer);
5720 if (!msave->rowtype)
5723 else if (lex_match_id (s->lexer, "OUTFILE"))
5725 lex_match (s->lexer, T_EQUALS);
5727 fh_unref (common.outfile);
5728 common.outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
5729 if (!common.outfile)
5732 else if (lex_match_id (s->lexer, "VARIABLES"))
5734 if (!parse_var_names (s->lexer, &common.variables))
5737 else if (lex_match_id (s->lexer, "FNAMES"))
5739 if (!parse_var_names (s->lexer, &common.fnames))
5742 else if (lex_match_id (s->lexer, "SNAMES"))
5744 if (!parse_var_names (s->lexer, &common.snames))
5747 else if (lex_match_id (s->lexer, "SPLIT"))
5749 lex_match (s->lexer, T_EQUALS);
5751 matrix_expr_destroy (msave->splits);
5752 msave->splits = matrix_parse_exp (s);
5756 else if (lex_match_id (s->lexer, "FACTOR"))
5758 lex_match (s->lexer, T_EQUALS);
5760 matrix_expr_destroy (msave->factors);
5761 msave->factors = matrix_parse_exp (s);
5762 if (!msave->factors)
5767 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
5768 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
5772 if (!msave->rowtype)
5774 lex_sbc_missing ("TYPE");
5777 common.has_splits = msave->splits || common.snames.n;
5778 common.has_factors = msave->factors || common.fnames.n;
5780 struct msave_common *c = s->common ? s->common : &common;
5781 if (c->fnames.n && !msave->factors)
5783 msg (SE, _("FNAMES requires FACTOR."));
5786 if (c->snames.n && !msave->splits)
5788 msg (SE, _("SNAMES requires SPLIT."));
5791 if (c->has_factors && !common.has_factors)
5793 msg (SE, _("%s is required because it was present on the first "
5794 "MSAVE in this MATRIX command."), "FACTOR");
5797 if (c->has_splits && !common.has_splits)
5799 msg (SE, _("%s is required because it was present on the first "
5800 "MSAVE in this MATRIX command."), "SPLIT");
5806 if (!common.outfile)
5808 lex_sbc_missing ("OUTFILE");
5811 s->common = xmemdup (&common, sizeof common);
5817 bool same = common.outfile == s->common->outfile;
5818 fh_unref (common.outfile);
5821 msg (SE, _("OUTFILE must name the same file on each MSAVE "
5822 "within a single MATRIX command."));
5826 if (!compare_variables ("VARIABLES",
5827 &common.variables, &s->common->variables)
5828 || !compare_variables ("FNAMES", &common.fnames, &s->common->fnames)
5829 || !compare_variables ("SNAMES", &common.snames, &s->common->snames))
5831 msave_common_uninit (&common);
5833 msave->common = s->common;
5834 if (!msave->varname_)
5835 msave->varname_ = xasprintf ("MAT%zu", ++s->common->n_varnames);
5839 msave_common_uninit (&common);
5840 matrix_cmd_destroy (cmd);
5845 matrix_expr_evaluate_vector (struct matrix_expr *e, const char *name)
5847 gsl_matrix *m = matrix_expr_evaluate (e);
5853 msg (SE, _("%s expression must evaluate to vector, "
5854 "not a %zu×%zu matrix."),
5855 name, m->size1, m->size2);
5856 gsl_matrix_free (m);
5860 return matrix_to_vector (m);
5864 msave_add_vars (struct dictionary *d, const struct string_array *vars)
5866 for (size_t i = 0; i < vars->n; i++)
5867 if (!dict_create_var (d, vars->strings[i], 0))
5868 return vars->strings[i];
5872 static struct dictionary *
5873 msave_create_dict (const struct msave_common *common)
5875 struct dictionary *dict = dict_create (get_default_encoding ());
5877 const char *dup_split = msave_add_vars (dict, &common->snames);
5880 msg (SE, _("Duplicate SPLIT variable name %s."), dup_split);
5884 dict_create_var_assert (dict, "ROWTYPE_", 8);
5886 const char *dup_factor = msave_add_vars (dict, &common->fnames);
5889 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
5893 dict_create_var_assert (dict, "VARNAME_", 8);
5895 const char *dup_var = msave_add_vars (dict, &common->variables);
5898 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
5910 matrix_cmd_execute_msave (struct msave_command *msave)
5912 struct msave_common *common = msave->common;
5913 gsl_matrix *m = NULL;
5914 gsl_vector *factors = NULL;
5915 gsl_vector *splits = NULL;
5917 m = matrix_expr_evaluate (msave->expr);
5921 if (!common->variables.n)
5922 for (size_t i = 0; i < m->size2; i++)
5923 string_array_append_nocopy (&common->variables,
5924 xasprintf ("COL%zu", i + 1));
5926 if (m->size2 != common->variables.n)
5928 msg (SE, _("Matrix on MSAVE has %zu columns instead of required %zu."),
5929 m->size2, common->variables.n);
5935 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
5939 if (!common->fnames.n)
5940 for (size_t i = 0; i < factors->size; i++)
5941 string_array_append_nocopy (&common->fnames,
5942 xasprintf ("FAC%zu", i + 1));
5944 if (factors->size != common->fnames.n)
5946 msg (SE, _("There are %zu factor variables, "
5947 "but %zu split values were supplied."),
5948 common->fnames.n, factors->size);
5955 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
5959 if (!common->fnames.n)
5960 for (size_t i = 0; i < splits->size; i++)
5961 string_array_append_nocopy (&common->fnames,
5962 xasprintf ("SPL%zu", i + 1));
5964 if (splits->size != common->fnames.n)
5966 msg (SE, _("There are %zu split variables, "
5967 "but %zu split values were supplied."),
5968 common->fnames.n, splits->size);
5973 if (!common->writer)
5975 struct dictionary *dict = msave_create_dict (common);
5979 common->writer = any_writer_open (common->outfile, dict);
5980 if (!common->writer)
5986 common->dict = dict;
5989 for (size_t y = 0; y < m->size1; y++)
5991 struct ccase *c = case_create (dict_get_proto (common->dict));
5994 /* Split variables */
5996 for (size_t i = 0; i < splits->size; i++)
5997 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
6000 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6001 msave->rowtype, ' ');
6005 for (size_t i = 0; i < factors->size; i++)
6006 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
6009 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
6010 msave->varname_, ' ');
6012 /* Continuous variables. */
6013 for (size_t x = 0; x < m->size2; x++)
6014 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
6015 casewriter_write (common->writer, c);
6019 gsl_matrix_free (m);
6020 gsl_vector_free (factors);
6021 gsl_vector_free (splits);
6024 static struct matrix_cmd *
6025 matrix_parse_mget (struct matrix_state *s)
6027 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6028 *cmd = (struct matrix_cmd) { .type = MCMD_MGET, .mget = { .state = s } };
6030 struct mget_command *mget = &cmd->mget;
6034 if (lex_match_id (s->lexer, "FILE"))
6036 lex_match (s->lexer, T_EQUALS);
6038 fh_unref (mget->file);
6039 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6043 else if (lex_match_id (s->lexer, "TYPE"))
6045 lex_match (s->lexer, T_EQUALS);
6046 while (lex_token (s->lexer) != T_SLASH
6047 && lex_token (s->lexer) != T_ENDCMD)
6049 const char *rowtype = match_rowtype (s->lexer);
6053 stringi_set_insert (&mget->rowtypes, rowtype);
6058 lex_error_expecting (s->lexer, "FILE", "TYPE");
6061 if (lex_token (s->lexer) == T_ENDCMD)
6064 if (!lex_force_match (s->lexer, T_SLASH))
6070 matrix_cmd_destroy (cmd);
6074 static const struct variable *
6075 get_a8_var (const struct dictionary *d, const char *name)
6077 const struct variable *v = dict_lookup_var (d, name);
6080 msg (SE, _("Matrix data file lacks %s variable."), name);
6083 if (var_get_width (v) != 8)
6085 msg (SE, _("%s variable in matrix data file must be "
6086 "8-byte string, but it has width %d."),
6087 name, var_get_width (v));
6094 is_rowtype (const union value *v, const char *rowtype)
6096 struct substring vs = ss_buffer (CHAR_CAST (char *, v->s), 8);
6097 ss_rtrim (&vs, ss_cstr (" "));
6098 return ss_equals_case (vs, ss_cstr (rowtype));
6102 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
6103 const struct dictionary *d,
6104 const struct variable *rowtype_var,
6105 struct matrix_state *s, size_t si, size_t fi,
6106 size_t cs, size_t cn)
6111 const union value *rowtype_ = case_data (rows[0], rowtype_var);
6112 const char *name_prefix = (is_rowtype (rowtype_, "COV") ? "CV"
6113 : is_rowtype (rowtype_, "CORR") ? "CR"
6114 : is_rowtype (rowtype_, "MEAN") ? "MN"
6115 : is_rowtype (rowtype_, "STDDEV") ? "SD"
6116 : is_rowtype (rowtype_, "N") ? "NC"
6119 struct string name = DS_EMPTY_INITIALIZER;
6120 ds_put_cstr (&name, name_prefix);
6122 ds_put_format (&name, "F%zu", fi);
6124 ds_put_format (&name, "S%zu", si);
6126 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
6128 mv = matrix_var_create (s, ds_ss (&name));
6131 msg (SW, _("Matrix data file contains variable with existing name %s."),
6136 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
6137 size_t n_missing = 0;
6138 for (size_t y = 0; y < n_rows; y++)
6140 for (size_t x = 0; x < cn; x++)
6142 struct variable *var = dict_get_var (d, cs + x);
6143 double value = case_num (rows[y], var);
6144 if (var_is_num_missing (var, value, MV_ANY))
6149 gsl_matrix_set (m, y, x, value);
6154 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
6155 "value, which was treated as zero.",
6156 "Matrix data file variable %s contains %zu missing "
6157 "values, which were treated as zero.", n_missing),
6158 ds_cstr (&name), n_missing);
6163 for (size_t y = 0; y < n_rows; y++)
6164 case_unref (rows[y]);
6168 var_changed (const struct ccase *ca, const struct ccase *cb,
6169 const struct variable *var)
6171 return !value_equal (case_data (ca, var), case_data (cb, var),
6172 var_get_width (var));
6176 vars_changed (const struct ccase *ca, const struct ccase *cb,
6177 const struct dictionary *d,
6178 size_t first_var, size_t n_vars)
6180 for (size_t i = 0; i < n_vars; i++)
6182 const struct variable *v = dict_get_var (d, first_var + i);
6183 if (var_changed (ca, cb, v))
6190 matrix_cmd_execute_mget (struct mget_command *mget)
6192 struct dictionary *d;
6193 struct casereader *r = any_reader_open_and_decode (mget->file, "UTF-8",
6198 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
6199 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
6200 if (!rowtype_ || !varname_)
6203 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
6205 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
6208 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
6210 msg (SE, _("Matrix data file contains no continuous variables."));
6214 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
6216 const struct variable *v = dict_get_var (d, i);
6217 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
6220 _("Matrix data file contains unexpected string variable %s."),
6226 /* SPLIT variables. */
6228 size_t sn = var_get_dict_index (rowtype_);
6229 struct ccase *sc = NULL;
6232 /* FACTOR variables. */
6233 size_t fs = var_get_dict_index (rowtype_) + 1;
6234 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
6235 struct ccase *fc = NULL;
6238 /* Continuous variables. */
6239 size_t cs = var_get_dict_index (varname_) + 1;
6240 size_t cn = dict_get_var_cnt (d) - cs;
6241 struct ccase *cc = NULL;
6244 struct ccase **rows = NULL;
6245 size_t allocated_rows = 0;
6249 while ((c = casereader_read (r)) != NULL)
6251 bool sd = vars_changed (sc, c, d, ss, sn);
6252 bool fd = sd || vars_changed (fc, c, d, fs, fn);
6253 bool md = fd || !cc || var_changed (cc, c, rowtype_) || var_changed (cc, c, varname_);
6268 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6269 mget->state, si, fi, cs, cn);
6275 if (n_rows >= allocated_rows)
6276 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
6279 matrix_mget_commit_var (rows, n_rows, d, rowtype_,
6280 mget->state, si, fi, cs, cn);
6284 casereader_destroy (r);
6288 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
6290 if (!lex_force_id (s->lexer))
6293 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
6295 *varp = matrix_var_create (s, lex_tokss (s->lexer));
6300 static struct matrix_cmd *
6301 matrix_parse_eigen (struct matrix_state *s)
6303 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6304 *cmd = (struct matrix_cmd) {
6306 .eigen = { .expr = NULL }
6309 struct eigen_command *eigen = &cmd->eigen;
6310 if (!lex_force_match (s->lexer, T_LPAREN))
6312 eigen->expr = matrix_parse_expr (s);
6314 || !lex_force_match (s->lexer, T_COMMA)
6315 || !matrix_parse_dst_var (s, &eigen->evec)
6316 || !lex_force_match (s->lexer, T_COMMA)
6317 || !matrix_parse_dst_var (s, &eigen->eval)
6318 || !lex_force_match (s->lexer, T_RPAREN))
6324 matrix_cmd_destroy (cmd);
6329 matrix_cmd_execute_eigen (struct eigen_command *eigen)
6331 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
6334 if (!is_symmetric (A))
6336 msg (SE, _("Argument of EIGEN must be symmetric."));
6337 gsl_matrix_free (A);
6341 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
6342 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
6343 gsl_vector v_eval = to_vector (eval);
6344 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
6345 gsl_eigen_symmv (A, &v_eval, evec, w);
6346 gsl_eigen_symmv_free (w);
6348 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
6350 gsl_matrix_free (A);
6352 gsl_matrix_free (eigen->eval->value);
6353 eigen->eval->value = eval;
6355 gsl_matrix_free (eigen->evec->value);
6356 eigen->evec->value = evec;
6359 static struct matrix_cmd *
6360 matrix_parse_setdiag (struct matrix_state *s)
6362 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6363 *cmd = (struct matrix_cmd) {
6364 .type = MCMD_SETDIAG,
6365 .setdiag = { .dst = NULL }
6368 struct setdiag_command *setdiag = &cmd->setdiag;
6369 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
6372 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
6375 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
6380 if (!lex_force_match (s->lexer, T_COMMA))
6383 setdiag->expr = matrix_parse_expr (s);
6387 if (!lex_force_match (s->lexer, T_RPAREN))
6393 matrix_cmd_destroy (cmd);
6398 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
6400 gsl_matrix *dst = setdiag->dst->value;
6403 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
6404 setdiag->dst->name);
6408 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
6412 size_t n = MIN (dst->size1, dst->size2);
6413 if (is_scalar (src))
6415 double d = to_scalar (src);
6416 for (size_t i = 0; i < n; i++)
6417 gsl_matrix_set (dst, i, i, d);
6419 else if (is_vector (src))
6421 gsl_vector v = to_vector (src);
6422 for (size_t i = 0; i < n && i < v.size; i++)
6423 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
6426 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
6427 "not a %zu×%zu matrix."),
6428 src->size1, src->size2);
6429 gsl_matrix_free (src);
6432 static struct matrix_cmd *
6433 matrix_parse_svd (struct matrix_state *s)
6435 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6436 *cmd = (struct matrix_cmd) {
6438 .svd = { .expr = NULL }
6441 struct svd_command *svd = &cmd->svd;
6442 if (!lex_force_match (s->lexer, T_LPAREN))
6444 svd->expr = matrix_parse_expr (s);
6446 || !lex_force_match (s->lexer, T_COMMA)
6447 || !matrix_parse_dst_var (s, &svd->u)
6448 || !lex_force_match (s->lexer, T_COMMA)
6449 || !matrix_parse_dst_var (s, &svd->s)
6450 || !lex_force_match (s->lexer, T_COMMA)
6451 || !matrix_parse_dst_var (s, &svd->v)
6452 || !lex_force_match (s->lexer, T_RPAREN))
6458 matrix_cmd_destroy (cmd);
6463 matrix_cmd_execute_svd (struct svd_command *svd)
6465 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
6469 if (m->size1 >= m->size2)
6472 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
6473 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
6474 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
6475 gsl_vector *work = gsl_vector_alloc (A->size2);
6476 gsl_linalg_SV_decomp (A, V, &Sv, work);
6477 gsl_vector_free (work);
6479 matrix_var_set (svd->u, A);
6480 matrix_var_set (svd->s, S);
6481 matrix_var_set (svd->v, V);
6485 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
6486 gsl_matrix_transpose_memcpy (At, m);
6487 gsl_matrix_free (m);
6489 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
6490 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
6491 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
6492 gsl_vector *work = gsl_vector_alloc (At->size2);
6493 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
6494 gsl_vector_free (work);
6496 matrix_var_set (svd->v, At);
6497 matrix_var_set (svd->s, St);
6498 matrix_var_set (svd->u, Vt);
6503 matrix_cmd_execute (struct matrix_cmd *cmd)
6508 matrix_cmd_execute_compute (&cmd->compute);
6512 matrix_cmd_execute_print (&cmd->print);
6516 return matrix_cmd_execute_do_if (&cmd->do_if);
6519 matrix_cmd_execute_loop (&cmd->loop);
6526 matrix_cmd_execute_display (&cmd->display);
6530 matrix_cmd_execute_release (&cmd->release);
6534 matrix_cmd_execute_save (&cmd->save);
6538 matrix_cmd_execute_read (&cmd->read);
6542 matrix_cmd_execute_write (&cmd->write);
6546 matrix_cmd_execute_get (&cmd->get);
6550 matrix_cmd_execute_msave (&cmd->msave);
6554 matrix_cmd_execute_mget (&cmd->mget);
6558 matrix_cmd_execute_eigen (&cmd->eigen);
6562 matrix_cmd_execute_setdiag (&cmd->setdiag);
6566 matrix_cmd_execute_svd (&cmd->svd);
6574 matrix_cmds_uninit (struct matrix_cmds *cmds)
6576 for (size_t i = 0; i < cmds->n; i++)
6577 matrix_cmd_destroy (cmds->commands[i]);
6578 free (cmds->commands);
6582 matrix_cmd_destroy (struct matrix_cmd *cmd)
6590 matrix_lvalue_destroy (cmd->compute.lvalue);
6591 matrix_expr_destroy (cmd->compute.rvalue);
6595 matrix_expr_destroy (cmd->print.expression);
6596 free (cmd->print.title);
6597 print_labels_destroy (cmd->print.rlabels);
6598 print_labels_destroy (cmd->print.clabels);
6602 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
6604 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
6605 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
6607 free (cmd->do_if.clauses);
6611 matrix_expr_destroy (cmd->loop.start);
6612 matrix_expr_destroy (cmd->loop.end);
6613 matrix_expr_destroy (cmd->loop.increment);
6614 matrix_expr_destroy (cmd->loop.top_condition);
6615 matrix_expr_destroy (cmd->loop.bottom_condition);
6616 matrix_cmds_uninit (&cmd->loop.commands);
6626 free (cmd->release.vars);
6630 matrix_expr_destroy (cmd->save.expression);
6631 fh_unref (cmd->save.outfile);
6632 string_array_destroy (cmd->save.variables);
6633 matrix_expr_destroy (cmd->save.names);
6634 stringi_set_destroy (&cmd->save.strings);
6638 matrix_lvalue_destroy (cmd->read.dst);
6639 matrix_expr_destroy (cmd->read.size);
6643 matrix_expr_destroy (cmd->write.expression);
6644 free (cmd->write.format);
6648 matrix_lvalue_destroy (cmd->get.dst);
6649 fh_unref (cmd->get.file);
6650 free (cmd->get.encoding);
6651 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
6655 free (cmd->msave.varname_);
6656 matrix_expr_destroy (cmd->msave.expr);
6657 matrix_expr_destroy (cmd->msave.factors);
6658 matrix_expr_destroy (cmd->msave.splits);
6662 fh_unref (cmd->mget.file);
6663 stringi_set_destroy (&cmd->mget.rowtypes);
6667 matrix_expr_destroy (cmd->eigen.expr);
6671 matrix_expr_destroy (cmd->setdiag.expr);
6675 matrix_expr_destroy (cmd->svd.expr);
6681 struct matrix_command_name
6684 struct matrix_cmd *(*parse) (struct matrix_state *);
6687 static const struct matrix_command_name *
6688 matrix_parse_command_name (struct lexer *lexer)
6690 static const struct matrix_command_name commands[] = {
6691 { "COMPUTE", matrix_parse_compute },
6692 { "CALL EIGEN", matrix_parse_eigen },
6693 { "CALL SETDIAG", matrix_parse_setdiag },
6694 { "CALL SVD", matrix_parse_svd },
6695 { "PRINT", matrix_parse_print },
6696 { "DO IF", matrix_parse_do_if },
6697 { "LOOP", matrix_parse_loop },
6698 { "BREAK", matrix_parse_break },
6699 { "READ", matrix_parse_read },
6700 { "WRITE", matrix_parse_write },
6701 { "GET", matrix_parse_get },
6702 { "SAVE", matrix_parse_save },
6703 { "MGET", matrix_parse_mget },
6704 { "MSAVE", matrix_parse_msave },
6705 { "DISPLAY", matrix_parse_display },
6706 { "RELEASE", matrix_parse_release },
6708 static size_t n = sizeof commands / sizeof *commands;
6710 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
6711 if (lex_match_phrase (lexer, c->name))
6716 static struct matrix_cmd *
6717 matrix_parse_command (struct matrix_state *s)
6719 size_t nesting_level = SIZE_MAX;
6721 struct matrix_cmd *c = NULL;
6722 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
6724 lex_error (s->lexer, _("Unknown matrix command."));
6725 else if (!cmd->parse)
6726 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
6730 nesting_level = output_open_group (
6731 group_item_create_nocopy (utf8_to_title (cmd->name),
6732 utf8_to_title (cmd->name)));
6737 lex_end_of_command (s->lexer);
6738 lex_discard_rest_of_command (s->lexer);
6739 if (nesting_level != SIZE_MAX)
6740 output_close_groups (nesting_level);
6746 cmd_matrix (struct lexer *lexer, struct dataset *ds)
6748 if (!lex_force_match (lexer, T_ENDCMD))
6751 struct matrix_state state = {
6753 .session = dataset_session (ds),
6755 .vars = HMAP_INITIALIZER (state.vars),
6760 while (lex_match (lexer, T_ENDCMD))
6762 if (lex_token (lexer) == T_STOP)
6764 msg (SE, _("Unexpected end of input expecting matrix command."));
6768 if (lex_match_phrase (lexer, "END MATRIX"))
6771 struct matrix_cmd *c = matrix_parse_command (&state);
6774 matrix_cmd_execute (c);
6775 matrix_cmd_destroy (c);
6779 struct matrix_var *var, *next;
6780 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
6783 gsl_matrix_free (var->value);
6784 hmap_delete (&state.vars, &var->hmap_node);
6787 hmap_destroy (&state.vars);
6788 fh_unref (state.prev_save_outfile);
6791 dict_unref (state.common->dict);
6792 casewriter_destroy (state.common->writer);
6793 free (state.common);
6795 fh_unref (state.prev_read_file);
6796 for (size_t i = 0; i < state.n_read_files; i++)
6797 read_file_destroy (state.read_files[i]);
6798 free (state.read_files);
6799 fh_unref (state.prev_write_file);
6800 for (size_t i = 0; i < state.n_write_files; i++)
6801 write_file_destroy (state.write_files[i]);
6802 free (state.write_files);