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_cdf.h>
21 #include <gsl/gsl_eigen.h>
22 #include <gsl/gsl_linalg.h>
23 #include <gsl/gsl_matrix.h>
24 #include <gsl/gsl_permutation.h>
25 #include <gsl/gsl_randist.h>
26 #include <gsl/gsl_vector.h>
31 #include "data/any-reader.h"
32 #include "data/any-writer.h"
33 #include "data/casereader.h"
34 #include "data/casewriter.h"
35 #include "data/data-in.h"
36 #include "data/data-out.h"
37 #include "data/dataset.h"
38 #include "data/dictionary.h"
39 #include "data/file-handle-def.h"
40 #include "language/command.h"
41 #include "language/data-io/data-reader.h"
42 #include "language/data-io/data-writer.h"
43 #include "language/data-io/file-handle.h"
44 #include "language/lexer/format-parser.h"
45 #include "language/lexer/lexer.h"
46 #include "language/lexer/variable-parser.h"
47 #include "libpspp/array.h"
48 #include "libpspp/assertion.h"
49 #include "libpspp/compiler.h"
50 #include "libpspp/hmap.h"
51 #include "libpspp/i18n.h"
52 #include "libpspp/intern.h"
53 #include "libpspp/misc.h"
54 #include "libpspp/str.h"
55 #include "libpspp/string-array.h"
56 #include "libpspp/stringi-set.h"
57 #include "libpspp/u8-line.h"
58 #include "math/distributions.h"
59 #include "math/random.h"
60 #include "output/driver.h"
61 #include "output/output-item.h"
62 #include "output/pivot-table.h"
64 #include "gl/c-ctype.h"
65 #include "gl/c-strcase.h"
66 #include "gl/ftoastr.h"
67 #include "gl/minmax.h"
71 #define _(msgid) gettext (msgid)
72 #define N_(msgid) (msgid)
76 struct hmap_node hmap_node;
84 struct file_handle *outfile;
85 struct string_array variables;
86 struct string_array fnames;
87 struct string_array snames;
90 /* Collects and owns factors and splits. The individual msave_command
91 structs point to these but do not own them. */
92 struct matrix_expr **factors;
93 size_t n_factors, allocated_factors;
94 struct matrix_expr **splits;
95 size_t n_splits, allocated_splits;
97 /* Execution state. */
98 struct dictionary *dict;
99 struct casewriter *writer;
104 struct file_handle *file;
105 struct dfm_reader *reader;
111 struct file_handle *file;
112 struct dfm_writer *writer;
114 struct u8_line *held;
119 struct file_handle *file;
120 struct dataset *dataset;
122 /* Parameters from parsing the first SAVE command for 'file'. */
123 struct string_array variables;
124 struct matrix_expr *names;
125 struct stringi_set strings;
127 /* Results from the first (only) attempt to open this save_file. */
129 struct casewriter *writer;
130 struct dictionary *dict;
135 struct dataset *dataset;
136 struct session *session;
140 struct msave_common *common;
142 struct file_handle *prev_read_file;
143 struct read_file **read_files;
146 struct file_handle *prev_write_file;
147 struct write_file **write_files;
148 size_t n_write_files;
150 struct file_handle *prev_save_file;
151 struct save_file **save_files;
155 static struct matrix_var *
156 matrix_var_lookup (struct matrix_state *s, struct substring name)
158 struct matrix_var *var;
160 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
161 utf8_hash_case_substring (name, 0), &s->vars)
162 if (!utf8_sscasecmp (ss_cstr (var->name), name))
168 static struct matrix_var *
169 matrix_var_create (struct matrix_state *s, struct substring name)
171 struct matrix_var *var = xmalloc (sizeof *var);
172 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
173 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
178 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
180 gsl_matrix_free (var->value);
184 /* The third argument to F() is a "prototype". For most prototypes, the first
185 letter (before the _) represents the return type and each other letter
186 (after the _) is an argument type. The types are:
188 - "m": A matrix of unrestricted dimensions.
192 - "v": A row or column vector.
194 - "e": Primarily for the first argument, this is a matrix with
195 unrestricted dimensions treated elementwise. Each element in the matrix
196 is passed to the implementation function separately.
198 The fourth argument is an optional constraints string. For this purpose the
199 first argument is named "a", the second "b", and so on. The following kinds
200 of constraints are supported. For matrix arguments, the constraints are
201 applied to each value in the matrix separately:
203 - "a(0,1)" or "a[0,1]": 0 < a < 1 or 0 <= a <= 1, respectively. Any
204 integer may substitute for 0 and 1. Half-open constraints (] and [) are
207 - "ai": Restrict a to integer values.
209 - "a>0", "a<0", "a>=0", "a<=0".
211 - "a<b", "a>b", "a<=b", "a>=b".
213 #define MATRIX_FUNCTIONS \
214 F(ABS, "ABS", m_e, NULL) \
215 F(ALL, "ALL", d_m, NULL) \
216 F(ANY, "ANY", d_m, NULL) \
217 F(ARSIN, "ARSIN", m_e, "a[-1,1]") \
218 F(ARTAN, "ARTAN", m_e, NULL) \
219 F(BLOCK, "BLOCK", m_any, NULL) \
220 F(CHOL, "CHOL", m_m, NULL) \
221 F(CMIN, "CMIN", m_m, NULL) \
222 F(CMAX, "CMAX", m_m, NULL) \
223 F(COS, "COS", m_e, NULL) \
224 F(CSSQ, "CSSQ", m_m, NULL) \
225 F(CSUM, "CSUM", m_m, NULL) \
226 F(DESIGN, "DESIGN", m_m, NULL) \
227 F(DET, "DET", d_m, NULL) \
228 F(DIAG, "DIAG", m_m, NULL) \
229 F(EVAL, "EVAL", m_m, NULL) \
230 F(EXP, "EXP", m_e, NULL) \
231 F(GINV, "GINV", m_m, NULL) \
232 F(GRADE, "GRADE", m_m, NULL) \
233 F(GSCH, "GSCH", m_m, NULL) \
234 F(IDENT, "IDENT", IDENT, NULL) \
235 F(INV, "INV", m_m, NULL) \
236 F(KRONEKER, "KRONEKER", m_mm, NULL) \
237 F(LG10, "LG10", m_e, "a>0") \
238 F(LN, "LN", m_e, "a>0") \
239 F(MAGIC, "MAGIC", m_d, "ai>=3") \
240 F(MAKE, "MAKE", m_ddd, NULL) \
241 F(MDIAG, "MDIAG", m_v, NULL) \
242 F(MMAX, "MMAX", d_m, NULL) \
243 F(MMIN, "MMIN", d_m, NULL) \
244 F(MOD, "MOD", m_md, NULL) \
245 F(MSSQ, "MSSQ", d_m, NULL) \
246 F(MSUM, "MSUM", d_m, NULL) \
247 F(NCOL, "NCOL", d_m, NULL) \
248 F(NROW, "NROW", d_m, NULL) \
249 F(RANK, "RANK", d_m, NULL) \
250 F(RESHAPE, "RESHAPE", m_mdd, NULL) \
251 F(RMAX, "RMAX", m_m, NULL) \
252 F(RMIN, "RMIN", m_m, NULL) \
253 F(RND, "RND", m_e, NULL) \
254 F(RNKORDER, "RNKORDER", m_m, NULL) \
255 F(RSSQ, "RSSQ", m_m, NULL) \
256 F(RSUM, "RSUM", m_m, NULL) \
257 F(SIN, "SIN", m_e, NULL) \
258 F(SOLVE, "SOLVE", m_mm, NULL) \
259 F(SQRT, "SQRT", m_e, "a>=0") \
260 F(SSCP, "SSCP", m_m, NULL) \
261 F(SVAL, "SVAL", m_m, NULL) \
262 F(SWEEP, "SWEEP", m_md, NULL) \
263 F(T, "T", m_m, NULL) \
264 F(TRACE, "TRACE", d_m, NULL) \
265 F(TRANSPOS, "TRANSPOS", m_m, NULL) \
266 F(TRUNC, "TRUNC", m_e, NULL) \
267 F(UNIFORM, "UNIFORM", m_dd, "ai>=0 bi>=0") \
268 F(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0") \
269 F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0") \
270 F(IDF_BETA, "IDF.BETA", m_edd, "a[0,1] b>0 c>0") \
271 F(RV_BETA, "RV.BETA", d_dd, "a>0 b>0") \
272 F(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
273 F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
274 F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]") \
275 F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]") \
276 F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0") \
277 F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0") \
278 F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "c>0") \
279 F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \
280 F(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0") \
281 F(CHICDF, "CHICDF", m_ed, "a>=0 b>0") \
282 F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0") \
283 F(PDF_CHISQ, "PDF.CHISQ", m_ed, "a>=0 b>0") \
284 F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \
285 F(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0") \
286 F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0") \
287 F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0") \
288 F(PDF_EXP, "PDF.EXP", m_ed, "a>=0 b>0") \
289 F(RV_EXP, "RV.EXP", d_d, "a>0") \
290 F(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0") \
291 F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \
292 F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0") \
293 F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0") \
294 F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0") \
295 F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0") \
296 F(RV_F, "RV.F", d_dd, "a>0 b>0") \
297 F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0") \
298 F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
299 F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0") \
300 F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
301 F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0") \
302 F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL) \
303 F(RV_LANDAU, "RV.LANDAU", d_none, NULL) \
304 F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0") \
305 F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0") \
306 F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0") \
307 F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0") \
308 F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]") \
309 F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]") \
310 F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0") \
311 F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0") \
312 F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0") \
313 F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0") \
314 F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
315 F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0") \
316 F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
317 F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0") \
318 F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0") \
319 F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0") \
320 F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0") \
321 F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0") \
322 F(CDFNORM, "CDFNORM", m_e, NULL) \
323 F(PROBIT, "PROBIT", m_e, "a(0,1)") \
324 F(NORMAL, "NORMAL", m_e, "a>0") \
325 F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0") \
326 F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0") \
327 F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0") \
328 F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0") \
329 F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0") \
330 F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0") \
331 F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0") \
332 F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0") \
333 F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0") \
334 F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0") \
335 F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL) \
336 F(RV_RTAIL, "RV.RTAIL", d_dd, NULL) \
337 F(CDF_T, "CDF.T", m_ed, "b>0") \
338 F(TCDF, "TCDF", m_ed, "b>0") \
339 F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0") \
340 F(PDF_T, "PDF.T", m_ed, "b>0") \
341 F(RV_T, "RV.T", d_d, "a>0") \
342 F(CDF_T1G, "CDF.T1G", m_edd, NULL) \
343 F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)") \
344 F(PDF_T1G, "PDF.T1G", m_edd, NULL) \
345 F(RV_T1G, "RV.T1G", d_dd, NULL) \
346 F(CDF_T2G, "CDF.T2G", m_edd, NULL) \
347 F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)") \
348 F(PDF_T2G, "PDF.T2G", m_edd, NULL) \
349 F(RV_T2G, "RV.T2G", d_dd, NULL) \
350 F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c") \
351 F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c") \
352 F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c") \
353 F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b") \
354 F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
355 F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0") \
356 F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
357 F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0") \
358 F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
359 F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
360 F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]") \
361 F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]") \
362 F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]") \
363 F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]") \
364 F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]") \
365 F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]") \
366 F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]") \
367 F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
368 F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
369 F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a") \
370 F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]") \
371 F(RV_LOG, "RV.LOG", d_d, "a(0,1]") \
372 F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
373 F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
374 F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]") \
375 F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0") \
376 F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0") \
377 F(RV_POISSON, "RV.POISSON", d_d, "a>0")
379 struct matrix_function_properties
382 const char *constraints;
385 enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 };
386 enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
387 enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
388 enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 };
389 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
390 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
391 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
392 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
393 enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 };
394 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
395 enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 };
396 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
397 enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
398 enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
399 enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
400 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
401 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
402 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
403 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
404 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
406 typedef double matrix_proto_d_none (void);
407 typedef double matrix_proto_d_d (double);
408 typedef double matrix_proto_d_dd (double, double);
409 typedef double matrix_proto_d_dd (double, double);
410 typedef double matrix_proto_d_ddd (double, double, double);
411 typedef gsl_matrix *matrix_proto_m_d (double);
412 typedef gsl_matrix *matrix_proto_m_dd (double, double);
413 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
414 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
415 typedef double matrix_proto_m_e (double);
416 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
417 typedef double matrix_proto_m_ed (double, double);
418 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
419 typedef double matrix_proto_m_edd (double, double, double);
420 typedef double matrix_proto_m_eddd (double, double, double, double);
421 typedef double matrix_proto_m_eed (double, double, double);
422 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
423 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
424 typedef double matrix_proto_d_m (gsl_matrix *);
425 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
426 typedef gsl_matrix *matrix_proto_IDENT (double, double);
428 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
429 static matrix_proto_##PROTO matrix_eval_##ENUM;
433 /* Matrix expressions. */
440 #define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
444 /* Elementwise and scalar arithmetic. */
445 MOP_NEGATE, /* unary - */
446 MOP_ADD_ELEMS, /* + */
447 MOP_SUB_ELEMS, /* - */
448 MOP_MUL_ELEMS, /* &* */
449 MOP_DIV_ELEMS, /* / and &/ */
450 MOP_EXP_ELEMS, /* &** */
452 MOP_SEQ_BY, /* a:b:c */
454 /* Matrix arithmetic. */
456 MOP_EXP_MAT, /* ** */
473 MOP_PASTE_HORZ, /* a, b, c, ... */
474 MOP_PASTE_VERT, /* a; b; c; ... */
478 MOP_VEC_INDEX, /* x(y) */
479 MOP_VEC_ALL, /* x(:) */
480 MOP_MAT_INDEX, /* x(y,z) */
481 MOP_ROW_INDEX, /* x(y,:) */
482 MOP_COL_INDEX, /* x(:,z) */
489 MOP_EOF, /* EOF('file') */
497 struct matrix_expr **subs;
502 struct matrix_var *variable;
503 struct read_file *eof;
508 matrix_expr_destroy (struct matrix_expr *e)
515 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
546 for (size_t i = 0; i < e->n_subs; i++)
547 matrix_expr_destroy (e->subs[i]);
559 static struct matrix_expr *
560 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
563 struct matrix_expr *e = xmalloc (sizeof *e);
564 *e = (struct matrix_expr) {
566 .subs = xmemdup (subs, n_subs * sizeof *subs),
572 static struct matrix_expr *
573 matrix_expr_create_0 (enum matrix_op op)
575 struct matrix_expr *sub;
576 return matrix_expr_create_subs (op, &sub, 0);
579 static struct matrix_expr *
580 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
582 return matrix_expr_create_subs (op, &sub, 1);
585 static struct matrix_expr *
586 matrix_expr_create_2 (enum matrix_op op,
587 struct matrix_expr *sub0, struct matrix_expr *sub1)
589 struct matrix_expr *subs[] = { sub0, sub1 };
590 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
593 static struct matrix_expr *
594 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
595 struct matrix_expr *sub1, struct matrix_expr *sub2)
597 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
598 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
601 static struct matrix_expr *
602 matrix_expr_create_number (double number)
604 struct matrix_expr *e = xmalloc (sizeof *e);
605 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
609 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
611 static struct matrix_expr *
612 matrix_parse_curly_comma (struct matrix_state *s)
614 struct matrix_expr *lhs = matrix_parse_or_xor (s);
618 while (lex_match (s->lexer, T_COMMA))
620 struct matrix_expr *rhs = matrix_parse_or_xor (s);
623 matrix_expr_destroy (lhs);
626 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
631 static struct matrix_expr *
632 matrix_parse_curly_semi (struct matrix_state *s)
634 if (lex_token (s->lexer) == T_RCURLY)
635 return matrix_expr_create_0 (MOP_EMPTY);
637 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
641 while (lex_match (s->lexer, T_SEMICOLON))
643 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
646 matrix_expr_destroy (lhs);
649 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
654 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
655 for (size_t Y = 0; Y < (M)->size1; Y++) \
656 for (size_t X = 0; X < (M)->size2; X++) \
657 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
660 is_vector (const gsl_matrix *m)
662 return m->size1 <= 1 || m->size2 <= 1;
666 to_vector (gsl_matrix *m)
668 return (m->size1 == 1
669 ? gsl_matrix_row (m, 0).vector
670 : gsl_matrix_column (m, 0).vector);
675 matrix_eval_ABS (double d)
681 matrix_eval_ALL (gsl_matrix *m)
683 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
690 matrix_eval_ANY (gsl_matrix *m)
692 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
699 matrix_eval_ARSIN (double d)
705 matrix_eval_ARTAN (double d)
711 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
715 for (size_t i = 0; i < n; i++)
720 gsl_matrix *block = gsl_matrix_calloc (r, c);
722 for (size_t i = 0; i < n; i++)
724 for (size_t y = 0; y < m[i]->size1; y++)
725 for (size_t x = 0; x < m[i]->size2; x++)
726 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
734 matrix_eval_CHOL (gsl_matrix *m)
736 if (!gsl_linalg_cholesky_decomp1 (m))
738 for (size_t y = 0; y < m->size1; y++)
739 for (size_t x = y + 1; x < m->size2; x++)
740 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
742 for (size_t y = 0; y < m->size1; y++)
743 for (size_t x = 0; x < y; x++)
744 gsl_matrix_set (m, y, x, 0);
749 msg (SE, _("Input to CHOL function is not positive-definite."));
755 matrix_eval_col_extremum (gsl_matrix *m, bool min)
760 return gsl_matrix_alloc (1, 0);
762 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
763 for (size_t x = 0; x < m->size2; x++)
765 double ext = gsl_matrix_get (m, 0, x);
766 for (size_t y = 1; y < m->size1; y++)
768 double value = gsl_matrix_get (m, y, x);
769 if (min ? value < ext : value > ext)
772 gsl_matrix_set (cext, 0, x, ext);
778 matrix_eval_CMAX (gsl_matrix *m)
780 return matrix_eval_col_extremum (m, false);
784 matrix_eval_CMIN (gsl_matrix *m)
786 return matrix_eval_col_extremum (m, true);
790 matrix_eval_COS (double d)
796 matrix_eval_col_sum (gsl_matrix *m, bool square)
801 return gsl_matrix_alloc (1, 0);
803 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
804 for (size_t x = 0; x < m->size2; x++)
807 for (size_t y = 0; y < m->size1; y++)
809 double d = gsl_matrix_get (m, y, x);
810 sum += square ? pow2 (d) : d;
812 gsl_matrix_set (result, 0, x, sum);
818 matrix_eval_CSSQ (gsl_matrix *m)
820 return matrix_eval_col_sum (m, true);
824 matrix_eval_CSUM (gsl_matrix *m)
826 return matrix_eval_col_sum (m, false);
830 compare_double_3way (const void *a_, const void *b_)
832 const double *a = a_;
833 const double *b = b_;
834 return *a < *b ? -1 : *a > *b;
838 matrix_eval_DESIGN (gsl_matrix *m)
840 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
841 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
842 gsl_matrix_transpose_memcpy (&m2, m);
844 for (size_t y = 0; y < m2.size1; y++)
845 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
847 size_t *n = xcalloc (m2.size1, sizeof *n);
849 for (size_t i = 0; i < m2.size1; i++)
851 double *row = tmp + m2.size2 * i;
852 for (size_t j = 0; j < m2.size2; )
855 for (k = j + 1; k < m2.size2; k++)
856 if (row[j] != row[k])
858 row[n[i]++] = row[j];
863 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
868 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
870 for (size_t i = 0; i < m->size2; i++)
875 const double *unique = tmp + m2.size2 * i;
876 for (size_t j = 0; j < n[i]; j++, x++)
878 double value = unique[j];
879 for (size_t y = 0; y < m->size1; y++)
880 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
891 matrix_eval_DET (gsl_matrix *m)
893 gsl_permutation *p = gsl_permutation_alloc (m->size1);
895 gsl_linalg_LU_decomp (m, p, &signum);
896 gsl_permutation_free (p);
897 return gsl_linalg_LU_det (m, signum);
901 matrix_eval_DIAG (gsl_matrix *m)
903 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
904 for (size_t i = 0; i < diag->size1; i++)
905 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
910 is_symmetric (const gsl_matrix *m)
912 if (m->size1 != m->size2)
915 for (size_t y = 0; y < m->size1; y++)
916 for (size_t x = 0; x < y; x++)
917 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
924 compare_double_desc (const void *a_, const void *b_)
926 const double *a = a_;
927 const double *b = b_;
928 return *a > *b ? -1 : *a < *b;
932 matrix_eval_EVAL (gsl_matrix *m)
934 if (!is_symmetric (m))
936 msg (SE, _("Argument of EVAL must be symmetric."));
940 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
941 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
942 gsl_vector v_eval = to_vector (eval);
943 gsl_eigen_symm (m, &v_eval, w);
944 gsl_eigen_symm_free (w);
946 assert (v_eval.stride == 1);
947 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
953 matrix_eval_EXP (double d)
958 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
961 Charl Linssen <charl@itfromb.it>
965 matrix_eval_GINV (gsl_matrix *A)
970 gsl_matrix *tmp_mat = NULL;
973 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
974 tmp_mat = gsl_matrix_alloc (m, n);
975 gsl_matrix_transpose_memcpy (tmp_mat, A);
983 gsl_matrix *V = gsl_matrix_alloc (m, m);
984 gsl_vector *u = gsl_vector_alloc (m);
986 gsl_vector *tmp_vec = gsl_vector_alloc (m);
987 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
988 gsl_vector_free (tmp_vec);
991 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
992 gsl_matrix_set_zero (Sigma_pinv);
993 double cutoff = 1e-15 * gsl_vector_max (u);
995 for (size_t i = 0; i < m; ++i)
997 double x = gsl_vector_get (u, i);
998 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
1001 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
1002 gsl_matrix *U = gsl_matrix_calloc (n, n);
1003 for (size_t i = 0; i < n; i++)
1004 for (size_t j = 0; j < m; j++)
1005 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
1007 /* two dot products to obtain pseudoinverse */
1008 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
1009 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
1014 A_pinv = gsl_matrix_alloc (n, m);
1015 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
1019 A_pinv = gsl_matrix_alloc (m, n);
1020 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
1023 gsl_matrix_free (tmp_mat);
1024 gsl_matrix_free (tmp_mat2);
1025 gsl_matrix_free (U);
1026 gsl_matrix_free (Sigma_pinv);
1027 gsl_vector_free (u);
1028 gsl_matrix_free (V);
1040 grade_compare_3way (const void *a_, const void *b_)
1042 const struct grade *a = a_;
1043 const struct grade *b = b_;
1045 return (a->value < b->value ? -1
1046 : a->value > b->value ? 1
1054 matrix_eval_GRADE (gsl_matrix *m)
1056 size_t n = m->size1 * m->size2;
1057 struct grade *grades = xmalloc (n * sizeof *grades);
1060 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1061 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
1062 qsort (grades, n, sizeof *grades, grade_compare_3way);
1064 for (size_t i = 0; i < n; i++)
1065 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
1073 dot (gsl_vector *a, gsl_vector *b)
1075 double result = 0.0;
1076 for (size_t i = 0; i < a->size; i++)
1077 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
1082 norm2 (gsl_vector *v)
1084 double result = 0.0;
1085 for (size_t i = 0; i < v->size; i++)
1086 result += pow2 (gsl_vector_get (v, i));
1091 norm (gsl_vector *v)
1093 return sqrt (norm2 (v));
1097 matrix_eval_GSCH (gsl_matrix *v)
1099 if (v->size2 < v->size1)
1101 msg (SE, _("GSCH requires its argument to have at least as many columns "
1102 "as rows, but it has dimensions %zu×%zu."),
1103 v->size1, v->size2);
1106 if (!v->size1 || !v->size2)
1109 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
1111 for (size_t vx = 0; vx < v->size2; vx++)
1113 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
1114 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
1116 gsl_vector_memcpy (&u_i, &v_i);
1117 for (size_t j = 0; j < ux; j++)
1119 gsl_vector u_j = gsl_matrix_column (u, j).vector;
1120 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
1121 for (size_t k = 0; k < u_i.size; k++)
1122 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
1123 - scale * gsl_vector_get (&u_j, k)));
1126 double len = norm (&u_i);
1129 gsl_vector_scale (&u_i, 1.0 / len);
1130 if (++ux >= v->size1)
1137 msg (SE, _("%zu×%zu argument to GSCH contains only "
1138 "%zu linearly independent columns."),
1139 v->size1, v->size2, ux);
1140 gsl_matrix_free (u);
1144 u->size2 = v->size1;
1149 matrix_eval_IDENT (double s1, double s2)
1151 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
1153 msg (SE, _("Arguments to IDENT must be non-negative integers."));
1157 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
1158 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1164 invert_matrix (gsl_matrix *x)
1166 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1168 gsl_linalg_LU_decomp (x, p, &signum);
1169 gsl_linalg_LU_invx (x, p);
1170 gsl_permutation_free (p);
1174 matrix_eval_INV (gsl_matrix *m)
1181 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1183 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1184 a->size2 * b->size2);
1186 for (size_t ar = 0; ar < a->size1; ar++)
1187 for (size_t br = 0; br < b->size1; br++, y++)
1190 for (size_t ac = 0; ac < a->size2; ac++)
1191 for (size_t bc = 0; bc < b->size2; bc++, x++)
1193 double av = gsl_matrix_get (a, ar, ac);
1194 double bv = gsl_matrix_get (b, br, bc);
1195 gsl_matrix_set (k, y, x, av * bv);
1202 matrix_eval_LG10 (double d)
1208 matrix_eval_LN (double d)
1214 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1216 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1219 for (size_t i = 1; i <= n * n; i++)
1221 gsl_matrix_set (m, y, x, i);
1223 size_t y1 = !y ? n - 1 : y - 1;
1224 size_t x1 = x + 1 >= n ? 0 : x + 1;
1225 if (gsl_matrix_get (m, y1, x1) == 0)
1231 y = y + 1 >= n ? 0 : y + 1;
1236 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1238 double a = gsl_matrix_get (m, y1, x1);
1239 double b = gsl_matrix_get (m, y2, x2);
1240 gsl_matrix_set (m, y1, x1, b);
1241 gsl_matrix_set (m, y2, x2, a);
1245 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1249 /* A. Umar, "On the Construction of Even Order Magic Squares",
1250 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1252 for (size_t i = 1; i <= n * n / 2; i++)
1254 gsl_matrix_set (m, y, x, i);
1264 for (size_t i = n * n; i > n * n / 2; i--)
1266 gsl_matrix_set (m, y, x, i);
1274 for (size_t y = 0; y < n; y++)
1275 for (size_t x = 0; x < n / 2; x++)
1277 unsigned int d = gsl_matrix_get (m, y, x);
1278 if (d % 2 != (y < n / 2))
1279 magic_exchange (m, y, x, y, n - x - 1);
1284 size_t x1 = n / 2 - 1;
1286 magic_exchange (m, y1, x1, y2, x1);
1287 magic_exchange (m, y1, x2, y2, x2);
1291 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1293 /* A. Umar, "On the Construction of Even Order Magic Squares",
1294 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1298 for (size_t i = 1; ; i++)
1300 gsl_matrix_set (m, y, x, i);
1301 if (++y == n / 2 - 1)
1313 for (size_t i = n * n; ; i--)
1315 gsl_matrix_set (m, y, x, i);
1316 if (++y == n / 2 - 1)
1325 for (size_t y = 0; y < n; y++)
1326 if (y != n / 2 - 1 && y != n / 2)
1327 for (size_t x = 0; x < n / 2; x++)
1329 unsigned int d = gsl_matrix_get (m, y, x);
1330 if (d % 2 != (y < n / 2))
1331 magic_exchange (m, y, x, y, n - x - 1);
1334 size_t a0 = (n * n - 2 * n) / 2 + 1;
1335 for (size_t i = 0; i < n / 2; i++)
1338 gsl_matrix_set (m, n / 2, i, a);
1339 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1341 for (size_t i = 0; i < n / 2; i++)
1343 size_t a = a0 + i + n / 2;
1344 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1345 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1347 for (size_t x = 1; x < n / 2; x += 2)
1348 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1349 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1350 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1351 size_t x1 = n / 2 - 2;
1352 size_t x2 = n / 2 + 1;
1353 size_t y1 = n / 2 - 2;
1354 size_t y2 = n / 2 + 1;
1355 magic_exchange (m, y1, x1, y2, x1);
1356 magic_exchange (m, y1, x2, y2, x2);
1360 matrix_eval_MAGIC (double n_)
1364 gsl_matrix *m = gsl_matrix_calloc (n, n);
1366 matrix_eval_MAGIC_odd (m, n);
1368 matrix_eval_MAGIC_singly_even (m, n);
1370 matrix_eval_MAGIC_doubly_even (m, n);
1375 matrix_eval_MAKE (double r, double c, double value)
1377 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1379 msg (SE, _("Size arguments to MAKE must be integers."));
1383 gsl_matrix *m = gsl_matrix_alloc (r, c);
1384 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1390 matrix_eval_MDIAG (gsl_vector *v)
1392 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1393 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1394 gsl_vector_memcpy (&diagonal, v);
1399 matrix_eval_MMAX (gsl_matrix *m)
1401 return gsl_matrix_max (m);
1405 matrix_eval_MMIN (gsl_matrix *m)
1407 return gsl_matrix_min (m);
1411 matrix_eval_MOD (gsl_matrix *m, double divisor)
1415 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1419 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1420 *d = fmod (*d, divisor);
1425 matrix_eval_MSSQ (gsl_matrix *m)
1428 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1434 matrix_eval_MSUM (gsl_matrix *m)
1437 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1443 matrix_eval_NCOL (gsl_matrix *m)
1449 matrix_eval_NROW (gsl_matrix *m)
1455 matrix_eval_RANK (gsl_matrix *m)
1457 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1458 gsl_linalg_QR_decomp (m, tau);
1459 gsl_vector_free (tau);
1461 return gsl_linalg_QRPT_rank (m, -1);
1465 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1467 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1469 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1474 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1476 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1477 "product of matrix dimensions."));
1481 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1484 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1486 gsl_matrix_set (dst, y1, x1, *d);
1497 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1502 return gsl_matrix_alloc (0, 1);
1504 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1505 for (size_t y = 0; y < m->size1; y++)
1507 double ext = gsl_matrix_get (m, y, 0);
1508 for (size_t x = 1; x < m->size2; x++)
1510 double value = gsl_matrix_get (m, y, x);
1511 if (min ? value < ext : value > ext)
1514 gsl_matrix_set (rext, y, 0, ext);
1520 matrix_eval_RMAX (gsl_matrix *m)
1522 return matrix_eval_row_extremum (m, false);
1526 matrix_eval_RMIN (gsl_matrix *m)
1528 return matrix_eval_row_extremum (m, true);
1532 matrix_eval_RND (double d)
1544 rank_compare_3way (const void *a_, const void *b_)
1546 const struct rank *a = a_;
1547 const struct rank *b = b_;
1549 return a->value < b->value ? -1 : a->value > b->value;
1553 matrix_eval_RNKORDER (gsl_matrix *m)
1555 size_t n = m->size1 * m->size2;
1556 struct rank *ranks = xmalloc (n * sizeof *ranks);
1558 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1559 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1560 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1562 for (size_t i = 0; i < n; )
1565 for (j = i + 1; j < n; j++)
1566 if (ranks[i].value != ranks[j].value)
1569 double rank = (i + j + 1.0) / 2.0;
1570 for (size_t k = i; k < j; k++)
1571 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1582 matrix_eval_row_sum (gsl_matrix *m, bool square)
1587 return gsl_matrix_alloc (0, 1);
1589 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1590 for (size_t y = 0; y < m->size1; y++)
1593 for (size_t x = 0; x < m->size2; x++)
1595 double d = gsl_matrix_get (m, y, x);
1596 sum += square ? pow2 (d) : d;
1598 gsl_matrix_set (result, y, 0, sum);
1604 matrix_eval_RSSQ (gsl_matrix *m)
1606 return matrix_eval_row_sum (m, true);
1610 matrix_eval_RSUM (gsl_matrix *m)
1612 return matrix_eval_row_sum (m, false);
1616 matrix_eval_SIN (double d)
1622 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1624 if (m1->size1 != m2->size1)
1626 msg (SE, _("SOLVE requires its arguments to have the same number of "
1627 "rows, but the first argument has dimensions %zu×%zu and "
1628 "the second %zu×%zu."),
1629 m1->size1, m1->size2,
1630 m2->size1, m2->size2);
1634 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1635 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1637 gsl_linalg_LU_decomp (m1, p, &signum);
1638 for (size_t i = 0; i < m2->size2; i++)
1640 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1641 gsl_vector xi = gsl_matrix_column (x, i).vector;
1642 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1644 gsl_permutation_free (p);
1649 matrix_eval_SQRT (double d)
1655 matrix_eval_SSCP (gsl_matrix *m)
1657 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1658 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1663 matrix_eval_SVAL (gsl_matrix *m)
1665 gsl_matrix *tmp_mat = NULL;
1666 if (m->size2 > m->size1)
1668 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1669 gsl_matrix_transpose_memcpy (tmp_mat, m);
1674 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1675 gsl_vector *S = gsl_vector_alloc (m->size2);
1676 gsl_vector *work = gsl_vector_alloc (m->size2);
1677 gsl_linalg_SV_decomp (m, V, S, work);
1679 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1680 for (size_t i = 0; i < m->size2; i++)
1681 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1683 gsl_matrix_free (V);
1684 gsl_vector_free (S);
1685 gsl_vector_free (work);
1686 gsl_matrix_free (tmp_mat);
1692 matrix_eval_SWEEP (gsl_matrix *m, double d)
1694 if (d < 1 || d > SIZE_MAX)
1696 msg (SE, _("Scalar argument to SWEEP must be integer."));
1700 if (k >= MIN (m->size1, m->size2))
1702 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1703 "equal to the smaller of the matrix argument's rows and "
1708 double m_kk = gsl_matrix_get (m, k, k);
1709 if (fabs (m_kk) > 1e-19)
1711 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1712 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1714 double m_ij = gsl_matrix_get (m, i, j);
1715 double m_ik = gsl_matrix_get (m, i, k);
1716 double m_kj = gsl_matrix_get (m, k, j);
1717 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1726 for (size_t i = 0; i < m->size1; i++)
1728 gsl_matrix_set (m, i, k, 0);
1729 gsl_matrix_set (m, k, i, 0);
1736 matrix_eval_TRACE (gsl_matrix *m)
1739 size_t n = MIN (m->size1, m->size2);
1740 for (size_t i = 0; i < n; i++)
1741 sum += gsl_matrix_get (m, i, i);
1746 matrix_eval_T (gsl_matrix *m)
1748 return matrix_eval_TRANSPOS (m);
1752 matrix_eval_TRANSPOS (gsl_matrix *m)
1754 if (m->size1 == m->size2)
1756 gsl_matrix_transpose (m);
1761 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1762 gsl_matrix_transpose_memcpy (t, m);
1768 matrix_eval_TRUNC (double d)
1774 matrix_eval_UNIFORM (double r_, double c_)
1778 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1780 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1784 gsl_matrix *m = gsl_matrix_alloc (r, c);
1785 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1786 *d = gsl_ran_flat (get_rng (), 0, 1);
1791 matrix_eval_PDF_BETA (double x, double a, double b)
1793 return gsl_ran_beta_pdf (x, a, b);
1797 matrix_eval_CDF_BETA (double x, double a, double b)
1799 return gsl_cdf_beta_P (x, a, b);
1803 matrix_eval_IDF_BETA (double P, double a, double b)
1805 return gsl_cdf_beta_Pinv (P, a, b);
1809 matrix_eval_RV_BETA (double a, double b)
1811 return gsl_ran_beta (get_rng (), a, b);
1815 matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
1817 return ncdf_beta (x, a, b, lambda);
1821 matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
1823 return npdf_beta (x, a, b, lambda);
1827 matrix_eval_CDF_BVNOR (double x0, double x1, double r)
1829 return cdf_bvnor (x0, x1, r);
1833 matrix_eval_PDF_BVNOR (double x0, double x1, double r)
1835 return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
1839 matrix_eval_CDF_CAUCHY (double x, double a, double b)
1841 return gsl_cdf_cauchy_P ((x - a) / b, 1);
1845 matrix_eval_IDF_CAUCHY (double P, double a, double b)
1847 return a + b * gsl_cdf_cauchy_Pinv (P, 1);
1851 matrix_eval_PDF_CAUCHY (double x, double a, double b)
1853 return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
1857 matrix_eval_RV_CAUCHY (double a, double b)
1859 return a + b * gsl_ran_cauchy (get_rng (), 1);
1863 matrix_eval_CDF_CHISQ (double x, double df)
1865 return gsl_cdf_chisq_P (x, df);
1869 matrix_eval_CHICDF (double x, double df)
1871 return matrix_eval_CDF_CHISQ (x, df);
1875 matrix_eval_IDF_CHISQ (double P, double df)
1877 return gsl_cdf_chisq_Pinv (P, df);
1881 matrix_eval_PDF_CHISQ (double x, double df)
1883 return gsl_ran_chisq_pdf (x, df);
1887 matrix_eval_RV_CHISQ (double df)
1889 return gsl_ran_chisq (get_rng (), df);
1893 matrix_eval_SIG_CHISQ (double x, double df)
1895 return gsl_cdf_chisq_Q (x, df);
1899 matrix_eval_CDF_EXP (double x, double a)
1901 return gsl_cdf_exponential_P (x, 1. / a);
1905 matrix_eval_IDF_EXP (double P, double a)
1907 return gsl_cdf_exponential_Pinv (P, 1. / a);
1911 matrix_eval_PDF_EXP (double x, double a)
1913 return gsl_ran_exponential_pdf (x, 1. / a);
1917 matrix_eval_RV_EXP (double a)
1919 return gsl_ran_exponential (get_rng (), 1. / a);
1923 matrix_eval_PDF_XPOWER (double x, double a, double b)
1925 return gsl_ran_exppow_pdf (x, a, b);
1929 matrix_eval_RV_XPOWER (double a, double b)
1931 return gsl_ran_exppow (get_rng (), a, b);
1935 matrix_eval_CDF_F (double x, double df1, double df2)
1937 return gsl_cdf_fdist_P (x, df1, df2);
1941 matrix_eval_FCDF (double x, double df1, double df2)
1943 return matrix_eval_CDF_F (x, df1, df2);
1947 matrix_eval_IDF_F (double P, double df1, double df2)
1949 return idf_fdist (P, df1, df2);
1953 matrix_eval_RV_F (double df1, double df2)
1955 return gsl_ran_fdist (get_rng (), df1, df2);
1959 matrix_eval_PDF_F (double x, double df1, double df2)
1961 return gsl_ran_fdist_pdf (x, df1, df2);
1965 matrix_eval_SIG_F (double x, double df1, double df2)
1967 return gsl_cdf_fdist_Q (x, df1, df2);
1971 matrix_eval_CDF_GAMMA (double x, double a, double b)
1973 return gsl_cdf_gamma_P (x, a, 1. / b);
1977 matrix_eval_IDF_GAMMA (double P, double a, double b)
1979 return gsl_cdf_gamma_Pinv (P, a, 1. / b);
1983 matrix_eval_PDF_GAMMA (double x, double a, double b)
1985 return gsl_ran_gamma_pdf (x, a, 1. / b);
1989 matrix_eval_RV_GAMMA (double a, double b)
1991 return gsl_ran_gamma (get_rng (), a, 1. / b);
1995 matrix_eval_PDF_LANDAU (double x)
1997 return gsl_ran_landau_pdf (x);
2001 matrix_eval_RV_LANDAU (void)
2003 return gsl_ran_landau (get_rng ());
2007 matrix_eval_CDF_LAPLACE (double x, double a, double b)
2009 return gsl_cdf_laplace_P ((x - a) / b, 1);
2013 matrix_eval_IDF_LAPLACE (double P, double a, double b)
2015 return a + b * gsl_cdf_laplace_Pinv (P, 1);
2019 matrix_eval_PDF_LAPLACE (double x, double a, double b)
2021 return gsl_ran_laplace_pdf ((x - a) / b, 1);
2025 matrix_eval_RV_LAPLACE (double a, double b)
2027 return a + b * gsl_ran_laplace (get_rng (), 1);
2031 matrix_eval_RV_LEVY (double c, double alpha)
2033 return gsl_ran_levy (get_rng (), c, alpha);
2037 matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
2039 return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
2043 matrix_eval_CDF_LOGISTIC (double x, double a, double b)
2045 return gsl_cdf_logistic_P ((x - a) / b, 1);
2049 matrix_eval_IDF_LOGISTIC (double P, double a, double b)
2051 return a + b * gsl_cdf_logistic_Pinv (P, 1);
2055 matrix_eval_PDF_LOGISTIC (double x, double a, double b)
2057 return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
2061 matrix_eval_RV_LOGISTIC (double a, double b)
2063 return a + b * gsl_ran_logistic (get_rng (), 1);
2067 matrix_eval_CDF_LNORMAL (double x, double m, double s)
2069 return gsl_cdf_lognormal_P (x, log (m), s);
2073 matrix_eval_IDF_LNORMAL (double P, double m, double s)
2075 return gsl_cdf_lognormal_Pinv (P, log (m), s);;
2079 matrix_eval_PDF_LNORMAL (double x, double m, double s)
2081 return gsl_ran_lognormal_pdf (x, log (m), s);
2085 matrix_eval_RV_LNORMAL (double m, double s)
2087 return gsl_ran_lognormal (get_rng (), log (m), s);
2091 matrix_eval_CDF_NORMAL (double x, double u, double s)
2093 return gsl_cdf_gaussian_P (x - u, s);
2097 matrix_eval_IDF_NORMAL (double P, double u, double s)
2099 return u + gsl_cdf_gaussian_Pinv (P, s);
2103 matrix_eval_PDF_NORMAL (double x, double u, double s)
2105 return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
2109 matrix_eval_RV_NORMAL (double u, double s)
2111 return u + gsl_ran_gaussian (get_rng (), s);
2115 matrix_eval_CDFNORM (double x)
2117 return gsl_cdf_ugaussian_P (x);
2121 matrix_eval_PROBIT (double P)
2123 return gsl_cdf_ugaussian_Pinv (P);
2127 matrix_eval_NORMAL (double s)
2129 return gsl_ran_gaussian (get_rng (), s);
2133 matrix_eval_PDF_NTAIL (double x, double a, double sigma)
2135 return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
2139 matrix_eval_RV_NTAIL (double a, double sigma)
2141 return gsl_ran_gaussian_tail (get_rng (), a, sigma);
2145 matrix_eval_CDF_PARETO (double x, double a, double b)
2147 return gsl_cdf_pareto_P (x, b, a);
2151 matrix_eval_IDF_PARETO (double P, double a, double b)
2153 return gsl_cdf_pareto_Pinv (P, b, a);
2157 matrix_eval_PDF_PARETO (double x, double a, double b)
2159 return gsl_ran_pareto_pdf (x, b, a);
2163 matrix_eval_RV_PARETO (double a, double b)
2165 return gsl_ran_pareto (get_rng (), b, a);
2169 matrix_eval_CDF_RAYLEIGH (double x, double sigma)
2171 return gsl_cdf_rayleigh_P (x, sigma);
2175 matrix_eval_IDF_RAYLEIGH (double P, double sigma)
2177 return gsl_cdf_rayleigh_Pinv (P, sigma);
2181 matrix_eval_PDF_RAYLEIGH (double x, double sigma)
2183 return gsl_ran_rayleigh_pdf (x, sigma);
2187 matrix_eval_RV_RAYLEIGH (double sigma)
2189 return gsl_ran_rayleigh (get_rng (), sigma);
2193 matrix_eval_PDF_RTAIL (double x, double a, double sigma)
2195 return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
2199 matrix_eval_RV_RTAIL (double a, double sigma)
2201 return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
2205 matrix_eval_CDF_T (double x, double df)
2207 return gsl_cdf_tdist_P (x, df);
2211 matrix_eval_TCDF (double x, double df)
2213 return matrix_eval_CDF_T (x, df);
2217 matrix_eval_IDF_T (double P, double df)
2219 return gsl_cdf_tdist_Pinv (P, df);
2223 matrix_eval_PDF_T (double x, double df)
2225 return gsl_ran_tdist_pdf (x, df);
2229 matrix_eval_RV_T (double df)
2231 return gsl_ran_tdist (get_rng (), df);
2235 matrix_eval_CDF_T1G (double x, double a, double b)
2237 return gsl_cdf_gumbel1_P (x, a, b);
2241 matrix_eval_IDF_T1G (double P, double a, double b)
2243 return gsl_cdf_gumbel1_Pinv (P, a, b);
2247 matrix_eval_PDF_T1G (double x, double a, double b)
2249 return gsl_ran_gumbel1_pdf (x, a, b);
2253 matrix_eval_RV_T1G (double a, double b)
2255 return gsl_ran_gumbel1 (get_rng (), a, b);
2259 matrix_eval_CDF_T2G (double x, double a, double b)
2261 return gsl_cdf_gumbel1_P (x, a, b);
2265 matrix_eval_IDF_T2G (double P, double a, double b)
2267 return gsl_cdf_gumbel1_Pinv (P, a, b);
2271 matrix_eval_PDF_T2G (double x, double a, double b)
2273 return gsl_ran_gumbel1_pdf (x, a, b);
2277 matrix_eval_RV_T2G (double a, double b)
2279 return gsl_ran_gumbel1 (get_rng (), a, b);
2283 matrix_eval_CDF_UNIFORM (double x, double a, double b)
2285 return gsl_cdf_flat_P (x, a, b);
2289 matrix_eval_IDF_UNIFORM (double P, double a, double b)
2291 return gsl_cdf_flat_Pinv (P, a, b);
2295 matrix_eval_PDF_UNIFORM (double x, double a, double b)
2297 return gsl_ran_flat_pdf (x, a, b);
2301 matrix_eval_RV_UNIFORM (double a, double b)
2303 return gsl_ran_flat (get_rng (), a, b);
2307 matrix_eval_CDF_WEIBULL (double x, double a, double b)
2309 return gsl_cdf_weibull_P (x, a, b);
2313 matrix_eval_IDF_WEIBULL (double P, double a, double b)
2315 return gsl_cdf_weibull_Pinv (P, a, b);
2319 matrix_eval_PDF_WEIBULL (double x, double a, double b)
2321 return gsl_ran_weibull_pdf (x, a, b);
2325 matrix_eval_RV_WEIBULL (double a, double b)
2327 return gsl_ran_weibull (get_rng (), a, b);
2331 matrix_eval_CDF_BERNOULLI (double k, double p)
2333 return k ? 1 : 1 - p;
2337 matrix_eval_PDF_BERNOULLI (double k, double p)
2339 return gsl_ran_bernoulli_pdf (k, p);
2343 matrix_eval_RV_BERNOULLI (double p)
2345 return gsl_ran_bernoulli (get_rng (), p);
2349 matrix_eval_CDF_BINOM (double k, double n, double p)
2351 return gsl_cdf_binomial_P (k, p, n);
2355 matrix_eval_PDF_BINOM (double k, double n, double p)
2357 return gsl_ran_binomial_pdf (k, p, n);
2361 matrix_eval_RV_BINOM (double n, double p)
2363 return gsl_ran_binomial (get_rng (), p, n);
2367 matrix_eval_CDF_GEOM (double k, double p)
2369 return gsl_cdf_geometric_P (k, p);
2373 matrix_eval_PDF_GEOM (double k, double p)
2375 return gsl_ran_geometric_pdf (k, p);
2379 matrix_eval_RV_GEOM (double p)
2381 return gsl_ran_geometric (get_rng (), p);
2385 matrix_eval_CDF_HYPER (double k, double a, double b, double c)
2387 return gsl_cdf_hypergeometric_P (k, c, a - c, b);
2391 matrix_eval_PDF_HYPER (double k, double a, double b, double c)
2393 return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
2397 matrix_eval_RV_HYPER (double a, double b, double c)
2399 return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
2403 matrix_eval_PDF_LOG (double k, double p)
2405 return gsl_ran_logarithmic_pdf (k, p);
2409 matrix_eval_RV_LOG (double p)
2411 return gsl_ran_logarithmic (get_rng (), p);
2415 matrix_eval_CDF_NEGBIN (double k, double n, double p)
2417 return gsl_cdf_negative_binomial_P (k, p, n);
2421 matrix_eval_PDF_NEGBIN (double k, double n, double p)
2423 return gsl_ran_negative_binomial_pdf (k, p, n);
2427 matrix_eval_RV_NEGBIN (double n, double p)
2429 return gsl_ran_negative_binomial (get_rng (), p, n);
2433 matrix_eval_CDF_POISSON (double k, double mu)
2435 return gsl_cdf_poisson_P (k, mu);
2439 matrix_eval_PDF_POISSON (double k, double mu)
2441 return gsl_ran_poisson_pdf (k, mu);
2445 matrix_eval_RV_POISSON (double mu)
2447 return gsl_ran_poisson (get_rng (), mu);
2450 struct matrix_function
2454 size_t min_args, max_args;
2457 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
2460 word_matches (const char **test, const char **name)
2462 size_t test_len = strcspn (*test, ".");
2463 size_t name_len = strcspn (*name, ".");
2464 if (test_len == name_len)
2466 if (buf_compare_case (*test, *name, test_len))
2469 else if (test_len < 3 || test_len > name_len)
2473 if (buf_compare_case (*test, *name, test_len))
2479 if (**test != **name)
2490 /* Returns 0 if TOKEN and FUNC do not match,
2491 1 if TOKEN is an acceptable abbreviation for FUNC,
2492 2 if TOKEN equals FUNC. */
2494 compare_function_names (const char *token_, const char *func_)
2496 const char *token = token_;
2497 const char *func = func_;
2498 while (*token || *func)
2499 if (!word_matches (&token, &func))
2501 return !c_strcasecmp (token_, func_) ? 2 : 1;
2504 static const struct matrix_function *
2505 matrix_parse_function_name (const char *token)
2507 static const struct matrix_function functions[] =
2509 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
2510 { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
2514 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
2516 for (size_t i = 0; i < N_FUNCTIONS; i++)
2518 if (compare_function_names (token, functions[i].name) > 0)
2519 return &functions[i];
2524 static struct read_file *
2525 read_file_create (struct matrix_state *s, struct file_handle *fh)
2527 for (size_t i = 0; i < s->n_read_files; i++)
2529 struct read_file *rf = s->read_files[i];
2537 struct read_file *rf = xmalloc (sizeof *rf);
2538 *rf = (struct read_file) { .file = fh };
2540 s->read_files = xrealloc (s->read_files,
2541 (s->n_read_files + 1) * sizeof *s->read_files);
2542 s->read_files[s->n_read_files++] = rf;
2546 static struct dfm_reader *
2547 read_file_open (struct read_file *rf)
2550 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
2555 read_file_destroy (struct read_file *rf)
2559 fh_unref (rf->file);
2560 dfm_close_reader (rf->reader);
2561 free (rf->encoding);
2566 static struct write_file *
2567 write_file_create (struct matrix_state *s, struct file_handle *fh)
2569 for (size_t i = 0; i < s->n_write_files; i++)
2571 struct write_file *wf = s->write_files[i];
2579 struct write_file *wf = xmalloc (sizeof *wf);
2580 *wf = (struct write_file) { .file = fh };
2582 s->write_files = xrealloc (s->write_files,
2583 (s->n_write_files + 1) * sizeof *s->write_files);
2584 s->write_files[s->n_write_files++] = wf;
2588 static struct dfm_writer *
2589 write_file_open (struct write_file *wf)
2592 wf->writer = dfm_open_writer (wf->file, wf->encoding);
2597 write_file_destroy (struct write_file *wf)
2603 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
2604 wf->held->s.ss.length);
2605 u8_line_destroy (wf->held);
2609 fh_unref (wf->file);
2610 dfm_close_writer (wf->writer);
2611 free (wf->encoding);
2617 matrix_parse_function (struct matrix_state *s, const char *token,
2618 struct matrix_expr **exprp)
2621 if (lex_next_token (s->lexer, 1) != T_LPAREN)
2624 if (lex_match_id (s->lexer, "EOF"))
2627 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
2631 if (!lex_force_match (s->lexer, T_RPAREN))
2637 struct read_file *rf = read_file_create (s, fh);
2639 struct matrix_expr *e = xmalloc (sizeof *e);
2640 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
2645 const struct matrix_function *f = matrix_parse_function_name (token);
2649 lex_get_n (s->lexer, 2);
2651 struct matrix_expr *e = xmalloc (sizeof *e);
2652 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
2654 if (lex_token (s->lexer) != T_RPAREN)
2656 size_t allocated_subs = 0;
2659 struct matrix_expr *sub = matrix_parse_expr (s);
2663 if (e->n_subs >= allocated_subs)
2664 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
2665 e->subs[e->n_subs++] = sub;
2667 while (lex_match (s->lexer, T_COMMA));
2669 if (!lex_force_match (s->lexer, T_RPAREN))
2672 if (e->n_subs < f->min_args || e->n_subs > f->max_args)
2674 if (f->min_args == f->max_args)
2675 msg (SE, ngettext ("Matrix function %s requires %zu argument.",
2676 "Matrix function %s requires %zu arguments.",
2678 f->name, f->min_args);
2679 else if (f->min_args == 1 && f->max_args == 2)
2680 msg (SE, ngettext ("Matrix function %s requires 1 or 2 arguments, "
2681 "but %zu was provided.",
2682 "Matrix function %s requires 1 or 2 arguments, "
2683 "but %zu were provided.",
2685 f->name, e->n_subs);
2686 else if (f->min_args == 1 && f->max_args == INT_MAX)
2687 msg (SE, _("Matrix function %s requires at least one argument."),
2699 matrix_expr_destroy (e);
2703 static struct matrix_expr *
2704 matrix_parse_primary (struct matrix_state *s)
2706 if (lex_is_number (s->lexer))
2708 double number = lex_number (s->lexer);
2710 return matrix_expr_create_number (number);
2712 else if (lex_is_string (s->lexer))
2714 char string[sizeof (double)];
2715 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
2719 memcpy (&number, string, sizeof number);
2720 return matrix_expr_create_number (number);
2722 else if (lex_match (s->lexer, T_LPAREN))
2724 struct matrix_expr *e = matrix_parse_or_xor (s);
2725 if (!e || !lex_force_match (s->lexer, T_RPAREN))
2727 matrix_expr_destroy (e);
2732 else if (lex_match (s->lexer, T_LCURLY))
2734 struct matrix_expr *e = matrix_parse_curly_semi (s);
2735 if (!e || !lex_force_match (s->lexer, T_RCURLY))
2737 matrix_expr_destroy (e);
2742 else if (lex_token (s->lexer) == T_ID)
2744 struct matrix_expr *retval;
2745 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
2748 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
2751 lex_error (s->lexer, _("Unknown variable %s."),
2752 lex_tokcstr (s->lexer));
2757 struct matrix_expr *e = xmalloc (sizeof *e);
2758 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
2761 else if (lex_token (s->lexer) == T_ALL)
2763 struct matrix_expr *retval;
2764 if (matrix_parse_function (s, "ALL", &retval))
2768 lex_error (s->lexer, NULL);
2772 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
2775 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
2777 if (lex_match (s->lexer, T_COLON))
2784 *indexp = matrix_parse_expr (s);
2785 return *indexp != NULL;
2789 static struct matrix_expr *
2790 matrix_parse_postfix (struct matrix_state *s)
2792 struct matrix_expr *lhs = matrix_parse_primary (s);
2793 if (!lhs || !lex_match (s->lexer, T_LPAREN))
2796 struct matrix_expr *i0;
2797 if (!matrix_parse_index_expr (s, &i0))
2799 matrix_expr_destroy (lhs);
2802 if (lex_match (s->lexer, T_RPAREN))
2804 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
2805 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
2806 else if (lex_match (s->lexer, T_COMMA))
2808 struct matrix_expr *i1;
2809 if (!matrix_parse_index_expr (s, &i1)
2810 || !lex_force_match (s->lexer, T_RPAREN))
2812 matrix_expr_destroy (lhs);
2813 matrix_expr_destroy (i0);
2814 matrix_expr_destroy (i1);
2817 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
2818 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
2819 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
2824 lex_error_expecting (s->lexer, "`)'", "`,'");
2829 static struct matrix_expr *
2830 matrix_parse_unary (struct matrix_state *s)
2832 if (lex_match (s->lexer, T_DASH))
2834 struct matrix_expr *lhs = matrix_parse_unary (s);
2835 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
2837 else if (lex_match (s->lexer, T_PLUS))
2838 return matrix_parse_unary (s);
2840 return matrix_parse_postfix (s);
2843 static struct matrix_expr *
2844 matrix_parse_seq (struct matrix_state *s)
2846 struct matrix_expr *start = matrix_parse_unary (s);
2847 if (!start || !lex_match (s->lexer, T_COLON))
2850 struct matrix_expr *end = matrix_parse_unary (s);
2853 matrix_expr_destroy (start);
2857 if (lex_match (s->lexer, T_COLON))
2859 struct matrix_expr *increment = matrix_parse_unary (s);
2862 matrix_expr_destroy (start);
2863 matrix_expr_destroy (end);
2866 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2869 return matrix_expr_create_2 (MOP_SEQ, start, end);
2872 static struct matrix_expr *
2873 matrix_parse_exp (struct matrix_state *s)
2875 struct matrix_expr *lhs = matrix_parse_seq (s);
2882 if (lex_match (s->lexer, T_EXP))
2884 else if (lex_match_phrase (s->lexer, "&**"))
2889 struct matrix_expr *rhs = matrix_parse_seq (s);
2892 matrix_expr_destroy (lhs);
2895 lhs = matrix_expr_create_2 (op, lhs, rhs);
2899 static struct matrix_expr *
2900 matrix_parse_mul_div (struct matrix_state *s)
2902 struct matrix_expr *lhs = matrix_parse_exp (s);
2909 if (lex_match (s->lexer, T_ASTERISK))
2911 else if (lex_match (s->lexer, T_SLASH))
2913 else if (lex_match_phrase (s->lexer, "&*"))
2915 else if (lex_match_phrase (s->lexer, "&/"))
2920 struct matrix_expr *rhs = matrix_parse_exp (s);
2923 matrix_expr_destroy (lhs);
2926 lhs = matrix_expr_create_2 (op, lhs, rhs);
2930 static struct matrix_expr *
2931 matrix_parse_add_sub (struct matrix_state *s)
2933 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2940 if (lex_match (s->lexer, T_PLUS))
2942 else if (lex_match (s->lexer, T_DASH))
2944 else if (lex_token (s->lexer) == T_NEG_NUM)
2949 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2952 matrix_expr_destroy (lhs);
2955 lhs = matrix_expr_create_2 (op, lhs, rhs);
2959 static struct matrix_expr *
2960 matrix_parse_relational (struct matrix_state *s)
2962 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2969 if (lex_match (s->lexer, T_GT))
2971 else if (lex_match (s->lexer, T_GE))
2973 else if (lex_match (s->lexer, T_LT))
2975 else if (lex_match (s->lexer, T_LE))
2977 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2979 else if (lex_match (s->lexer, T_NE))
2984 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2987 matrix_expr_destroy (lhs);
2990 lhs = matrix_expr_create_2 (op, lhs, rhs);
2994 static struct matrix_expr *
2995 matrix_parse_not (struct matrix_state *s)
2997 if (lex_match (s->lexer, T_NOT))
2999 struct matrix_expr *lhs = matrix_parse_not (s);
3000 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
3003 return matrix_parse_relational (s);
3006 static struct matrix_expr *
3007 matrix_parse_and (struct matrix_state *s)
3009 struct matrix_expr *lhs = matrix_parse_not (s);
3013 while (lex_match (s->lexer, T_AND))
3015 struct matrix_expr *rhs = matrix_parse_not (s);
3018 matrix_expr_destroy (lhs);
3021 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
3026 static struct matrix_expr *
3027 matrix_parse_or_xor (struct matrix_state *s)
3029 struct matrix_expr *lhs = matrix_parse_and (s);
3036 if (lex_match (s->lexer, T_OR))
3038 else if (lex_match_id (s->lexer, "XOR"))
3043 struct matrix_expr *rhs = matrix_parse_and (s);
3046 matrix_expr_destroy (lhs);
3049 lhs = matrix_expr_create_2 (op, lhs, rhs);
3053 static struct matrix_expr *
3054 matrix_parse_expr (struct matrix_state *s)
3056 return matrix_parse_or_xor (s);
3059 /* Expression evaluation. */
3062 matrix_op_eval (enum matrix_op op, double a, double b)
3066 case MOP_ADD_ELEMS: return a + b;
3067 case MOP_SUB_ELEMS: return a - b;
3068 case MOP_MUL_ELEMS: return a * b;
3069 case MOP_DIV_ELEMS: return a / b;
3070 case MOP_EXP_ELEMS: return pow (a, b);
3071 case MOP_GT: return a > b;
3072 case MOP_GE: return a >= b;
3073 case MOP_LT: return a < b;
3074 case MOP_LE: return a <= b;
3075 case MOP_EQ: return a == b;
3076 case MOP_NE: return a != b;
3077 case MOP_AND: return (a > 0) && (b > 0);
3078 case MOP_OR: return (a > 0) || (b > 0);
3079 case MOP_XOR: return (a > 0) != (b > 0);
3081 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3090 case MOP_PASTE_HORZ:
3091 case MOP_PASTE_VERT:
3107 matrix_op_name (enum matrix_op op)
3111 case MOP_ADD_ELEMS: return "+";
3112 case MOP_SUB_ELEMS: return "-";
3113 case MOP_MUL_ELEMS: return "&*";
3114 case MOP_DIV_ELEMS: return "&/";
3115 case MOP_EXP_ELEMS: return "&**";
3116 case MOP_GT: return ">";
3117 case MOP_GE: return ">=";
3118 case MOP_LT: return "<";
3119 case MOP_LE: return "<=";
3120 case MOP_EQ: return "=";
3121 case MOP_NE: return "<>";
3122 case MOP_AND: return "AND";
3123 case MOP_OR: return "OR";
3124 case MOP_XOR: return "XOR";
3126 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3135 case MOP_PASTE_HORZ:
3136 case MOP_PASTE_VERT:
3152 is_scalar (const gsl_matrix *m)
3154 return m->size1 == 1 && m->size2 == 1;
3158 to_scalar (const gsl_matrix *m)
3160 assert (is_scalar (m));
3161 return gsl_matrix_get (m, 0, 0);
3165 matrix_expr_evaluate_elementwise (enum matrix_op op,
3166 gsl_matrix *a, gsl_matrix *b)
3170 double be = to_scalar (b);
3171 for (size_t r = 0; r < a->size1; r++)
3172 for (size_t c = 0; c < a->size2; c++)
3174 double *ae = gsl_matrix_ptr (a, r, c);
3175 *ae = matrix_op_eval (op, *ae, be);
3179 else if (is_scalar (a))
3181 double ae = to_scalar (a);
3182 for (size_t r = 0; r < b->size1; r++)
3183 for (size_t c = 0; c < b->size2; c++)
3185 double *be = gsl_matrix_ptr (b, r, c);
3186 *be = matrix_op_eval (op, ae, *be);
3190 else if (a->size1 == b->size1 && a->size2 == b->size2)
3192 for (size_t r = 0; r < a->size1; r++)
3193 for (size_t c = 0; c < a->size2; c++)
3195 double *ae = gsl_matrix_ptr (a, r, c);
3196 double be = gsl_matrix_get (b, r, c);
3197 *ae = matrix_op_eval (op, *ae, be);
3203 msg (SE, _("Operands to %s must have the same dimensions or one "
3204 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
3205 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
3211 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
3213 if (is_scalar (a) || is_scalar (b))
3214 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
3216 if (a->size2 != b->size1)
3218 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
3219 "not conformable for multiplication."),
3220 a->size1, a->size2, b->size1, b->size2);
3224 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3225 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3230 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3232 gsl_matrix *tmp = *a;
3238 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3241 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3242 swap_matrix (z, tmp);
3246 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3248 mul_matrix (x, *x, *x, tmp);
3252 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
3255 if (x->size1 != x->size2)
3257 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
3258 "the left-hand size, not one with dimensions %zu×%zu."),
3259 x->size1, x->size2);
3264 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
3265 "right-hand side, not a matrix with dimensions %zu×%zu."),
3266 b->size1, b->size2);
3269 double bf = to_scalar (b);
3270 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3272 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
3273 "or outside the valid range."), bf);
3278 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3280 gsl_matrix_set_identity (y);
3284 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3286 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3289 mul_matrix (&y, x, y, &t);
3290 square_matrix (&x, &t);
3293 square_matrix (&x, &t);
3295 mul_matrix (&y, x, y, &t);
3299 /* Garbage collection.
3301 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3302 a permutation of them. We are returning one of them; that one must not be
3303 destroyed. We must not destroy 'x_' because the caller owns it. */
3305 gsl_matrix_free (y_);
3307 gsl_matrix_free (t_);
3313 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
3316 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3318 msg (SE, _("All operands of : operator must be scalars."));
3322 long int start = to_scalar (start_);
3323 long int end = to_scalar (end_);
3324 long int by = by_ ? to_scalar (by_) : 1;
3328 msg (SE, _("The increment operand to : must be nonzero."));
3332 long int n = (end >= start && by > 0 ? (end - start + by) / by
3333 : end <= start && by < 0 ? (start - end - by) / -by
3335 gsl_matrix *m = gsl_matrix_alloc (1, n);
3336 for (long int i = 0; i < n; i++)
3337 gsl_matrix_set (m, 0, i, start + i * by);
3342 matrix_expr_evaluate_not (gsl_matrix *a)
3344 for (size_t r = 0; r < a->size1; r++)
3345 for (size_t c = 0; c < a->size2; c++)
3347 double *ae = gsl_matrix_ptr (a, r, c);
3354 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
3356 if (a->size1 != b->size1)
3358 if (!a->size1 || !a->size2)
3360 else if (!b->size1 || !b->size2)
3363 msg (SE, _("All columns in a matrix must have the same number of rows, "
3364 "but this tries to paste matrices with %zu and %zu rows."),
3365 a->size1, b->size1);
3369 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3370 for (size_t y = 0; y < a->size1; y++)
3372 for (size_t x = 0; x < a->size2; x++)
3373 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3374 for (size_t x = 0; x < b->size2; x++)
3375 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3381 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
3383 if (a->size2 != b->size2)
3385 if (!a->size1 || !a->size2)
3387 else if (!b->size1 || !b->size2)
3390 msg (SE, _("All rows in a matrix must have the same number of columns, "
3391 "but this tries to stack matrices with %zu and %zu columns."),
3392 a->size2, b->size2);
3396 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3397 for (size_t x = 0; x < a->size2; x++)
3399 for (size_t y = 0; y < a->size1; y++)
3400 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3401 for (size_t y = 0; y < b->size1; y++)
3402 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3408 matrix_to_vector (gsl_matrix *m)
3411 gsl_vector v = to_vector (m);
3412 assert (v.block == m->block || !v.block);
3416 gsl_matrix_free (m);
3417 return xmemdup (&v, sizeof v);
3431 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3434 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
3435 enum index_type index_type, size_t other_size,
3436 struct index_vector *iv)
3445 msg (SE, _("Vector index must be scalar or vector, not a "
3447 m->size1, m->size2);
3451 msg (SE, _("Matrix row index must be scalar or vector, not a "
3453 m->size1, m->size2);
3457 msg (SE, _("Matrix column index must be scalar or vector, not a "
3459 m->size1, m->size2);
3465 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3466 *iv = (struct index_vector) {
3467 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3470 for (size_t i = 0; i < v.size; i++)
3472 double index = gsl_vector_get (&v, i);
3473 if (index < 1 || index >= size + 1)
3478 msg (SE, _("Index %g is out of range for vector "
3479 "with %zu elements."), index, size);
3483 msg (SE, _("%g is not a valid row index for "
3484 "a %zu×%zu matrix."),
3485 index, size, other_size);
3489 msg (SE, _("%g is not a valid column index for "
3490 "a %zu×%zu matrix."),
3491 index, other_size, size);
3498 iv->indexes[i] = index - 1;
3504 *iv = (struct index_vector) {
3505 .indexes = xnmalloc (size, sizeof *iv->indexes),
3508 for (size_t i = 0; i < size; i++)
3515 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
3517 if (!is_vector (sm))
3519 msg (SE, _("Vector index operator may not be applied to "
3520 "a %zu×%zu matrix."),
3521 sm->size1, sm->size2);
3529 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
3531 if (!matrix_expr_evaluate_vec_all (sm))
3534 gsl_vector sv = to_vector (sm);
3535 struct index_vector iv;
3536 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
3539 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3540 sm->size1 == 1 ? iv.n : 1);
3541 gsl_vector dv = to_vector (dm);
3542 for (size_t dx = 0; dx < iv.n; dx++)
3544 size_t sx = iv.indexes[dx];
3545 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3553 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
3556 struct index_vector iv0;
3557 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
3560 struct index_vector iv1;
3561 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
3568 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3569 for (size_t dy = 0; dy < iv0.n; dy++)
3571 size_t sy = iv0.indexes[dy];
3573 for (size_t dx = 0; dx < iv1.n; dx++)
3575 size_t sx = iv1.indexes[dx];
3576 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3584 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3585 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3586 const struct matrix_function_properties *, gsl_matrix *[], size_t, \
3587 matrix_proto_##PROTO *);
3592 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
3594 if (!is_scalar (subs[index]))
3596 msg (SE, _("Function %s argument %zu must be a scalar, "
3597 "not a %zu×%zu matrix."),
3598 name, index + 1, subs[index]->size1, subs[index]->size2);
3605 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
3607 if (!is_vector (subs[index]))
3609 msg (SE, _("Function %s argument %zu must be a vector, "
3610 "not a %zu×%zu matrix."),
3611 name, index + 1, subs[index]->size1, subs[index]->size2);
3618 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
3620 for (size_t i = 0; i < n_subs; i++)
3622 if (!check_scalar_arg (name, subs, i))
3624 d[i] = to_scalar (subs[i]);
3630 parse_constraint_value (const char **constraintsp)
3633 long retval = strtol (*constraintsp, &tail, 10);
3634 assert (tail > *constraintsp);
3635 *constraintsp = tail;
3640 argument_invalid (const struct matrix_function_properties *props,
3641 size_t arg_index, gsl_matrix *a, size_t y, size_t x,
3645 ds_put_format (out, _("Argument %zu to matrix function %s "
3646 "has invalid value %g."),
3647 arg_index, props->name, gsl_matrix_get (a, y, x));
3649 ds_put_format (out, _("Row %zu, column %zu of argument %zu "
3650 "to matrix function %s has "
3651 "invalid value %g."),
3652 y + 1, x + 1, arg_index, props->name,
3653 gsl_matrix_get (a, y, x));
3657 argument_inequality_error (const struct matrix_function_properties *props,
3658 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3659 size_t b_index, double b,
3660 bool greater, bool equal)
3662 struct string s = DS_EMPTY_INITIALIZER;
3664 argument_invalid (props, a_index, a, y, x, &s);
3665 ds_put_cstr (&s, " ");
3666 if (greater && equal)
3667 ds_put_format (&s, _("This argument must be greater than or "
3668 "equal to argument %zu, but that has value %g."),
3670 else if (greater && !equal)
3671 ds_put_format (&s, _("This argument must be greater than argument %zu, "
3672 "but that has value %g."),
3675 ds_put_format (&s, _("This argument must be less than or "
3676 "equal to argument %zu, but that has value %g."),
3679 ds_put_format (&s, _("This argument must be less than argument %zu, "
3680 "but that has value %g."),
3682 msg (ME, ds_cstr (&s));
3687 argument_value_error (const struct matrix_function_properties *props,
3688 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3690 bool greater, bool equal)
3692 struct string s = DS_EMPTY_INITIALIZER;
3693 argument_invalid (props, a_index, a, y, x, &s);
3694 ds_put_cstr (&s, " ");
3695 if (greater && equal)
3696 ds_put_format (&s, _("This argument must be greater than or equal to %g."),
3698 else if (greater && !equal)
3699 ds_put_format (&s, _("This argument must be greater than %g."), b);
3701 ds_put_format (&s, _("This argument must be less than or equal to %g."), b);
3703 ds_put_format (&s, _("This argument must be less than %g."), b);
3704 msg (ME, ds_cstr (&s));
3709 check_constraints (const struct matrix_function_properties *props,
3710 gsl_matrix *args[], size_t n_args)
3712 const char *constraints = props->constraints;
3716 size_t arg_index = SIZE_MAX;
3717 while (*constraints)
3719 if (*constraints >= 'a' && *constraints <= 'd')
3721 arg_index = *constraints++ - 'a';
3722 assert (arg_index < n_args);
3724 else if (*constraints == '[' || *constraints == '(')
3726 assert (arg_index < n_args);
3727 bool open_lower = *constraints++ == '(';
3728 int minimum = parse_constraint_value (&constraints);
3729 assert (*constraints == ',');
3731 int maximum = parse_constraint_value (&constraints);
3732 assert (*constraints == ']' || *constraints == ')');
3733 bool open_upper = *constraints++ == ')';
3735 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3736 if ((open_lower ? *d <= minimum : *d < minimum)
3737 || (open_upper ? *d >= maximum : *d > maximum))
3739 if (!is_scalar (args[arg_index]))
3740 msg (ME, _("Row %zu, column %zu of argument %zu to matrix "
3741 "function %s has value %g, which is outside "
3742 "the valid range %c%d,%d%c."),
3743 y + 1, x + 1, arg_index + 1, props->name, *d,
3744 open_lower ? '(' : '[',
3746 open_upper ? ')' : ']');
3748 msg (ME, _("Argument %zu to matrix function %s has value %g, "
3749 "which is outside the valid range %c%d,%d%c."),
3750 arg_index + 1, props->name, *d,
3751 open_lower ? '(' : '[',
3753 open_upper ? ')' : ']');
3757 else if (*constraints == '>' || *constraints == '<')
3759 bool greater = *constraints++ == '>';
3761 if (*constraints == '=')
3769 if (*constraints >= 'a' && *constraints <= 'd')
3771 size_t a_index = arg_index;
3772 size_t b_index = *constraints - 'a';
3773 assert (a_index < n_args);
3774 assert (b_index < n_args);
3776 /* We only support one of the two arguments being non-scalar.
3777 It's easier to support only the first one being non-scalar, so
3778 flip things around if it's the other way. */
3779 if (!is_scalar (args[b_index]))
3781 assert (is_scalar (args[a_index]));
3782 size_t tmp_index = a_index;
3784 b_index = tmp_index;
3789 double b = to_scalar (args[b_index]);
3790 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
3792 ? (equal ? !(*a >= b) : !(*a > b))
3793 : (equal ? !(*a <= b) : !(*a < b)))
3795 argument_inequality_error (props,
3796 a_index + 1, args[a_index], y, x,
3804 int comparand = parse_constraint_value (&constraints);
3806 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3808 ? (equal ? !(*d >= comparand) : !(*d > comparand))
3809 : (equal ? !(*d <= comparand) : !(*d < comparand)))
3811 argument_value_error (props,
3812 arg_index + 1, args[arg_index], y, x,
3821 assert (*constraints == ' ');
3823 arg_index = SIZE_MAX;
3830 matrix_expr_evaluate_d_none (
3831 const struct matrix_function_properties *props UNUSED,
3832 gsl_matrix *subs[] UNUSED, size_t n_subs,
3833 matrix_proto_d_none *f)
3835 assert (n_subs == 0);
3837 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3838 gsl_matrix_set (m, 0, 0, f ());
3843 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
3844 gsl_matrix *subs[], size_t n_subs,
3845 matrix_proto_d_d *f)
3847 assert (n_subs == 1);
3850 if (!to_scalar_args (props->name, subs, n_subs, &d))
3853 if (!check_constraints (props, subs, n_subs))
3856 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3857 gsl_matrix_set (m, 0, 0, f (d));
3862 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
3863 gsl_matrix *subs[], size_t n_subs,
3864 matrix_proto_d_dd *f)
3866 assert (n_subs == 2);
3869 if (!to_scalar_args (props->name, subs, n_subs, d))
3872 if (!check_constraints (props, subs, n_subs))
3875 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3876 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
3881 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
3882 gsl_matrix *subs[], size_t n_subs,
3883 matrix_proto_d_ddd *f)
3885 assert (n_subs == 3);
3888 if (!to_scalar_args (props->name, subs, n_subs, d))
3891 if (!check_constraints (props, subs, n_subs))
3894 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3895 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
3900 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
3901 gsl_matrix *subs[], size_t n_subs,
3902 matrix_proto_m_d *f)
3904 assert (n_subs == 1);
3907 return to_scalar_args (props->name, subs, n_subs, &d) ? f (d) : NULL;
3911 matrix_expr_evaluate_m_dd (const struct matrix_function_properties *props,
3912 gsl_matrix *subs[], size_t n_subs,
3913 matrix_proto_m_dd *f)
3915 assert (n_subs == 2);
3918 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
3922 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
3923 gsl_matrix *subs[], size_t n_subs,
3924 matrix_proto_m_ddd *f)
3926 assert (n_subs == 3);
3929 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
3933 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED,
3934 gsl_matrix *subs[], size_t n_subs,
3935 matrix_proto_m_m *f)
3937 assert (n_subs == 1);
3942 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
3943 gsl_matrix *subs[], size_t n_subs,
3944 matrix_proto_m_e *f)
3946 assert (n_subs == 1);
3948 if (!check_constraints (props, subs, n_subs))
3951 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3957 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED,
3958 gsl_matrix *subs[], size_t n_subs,
3959 matrix_proto_m_md *f)
3961 assert (n_subs == 2);
3962 if (!check_scalar_arg (props->name, subs, 1))
3964 return f (subs[0], to_scalar (subs[1]));
3968 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
3969 gsl_matrix *subs[], size_t n_subs,
3970 matrix_proto_m_ed *f)
3972 assert (n_subs == 2);
3973 if (!check_scalar_arg (props->name, subs, 1))
3976 if (!check_constraints (props, subs, n_subs))
3979 double b = to_scalar (subs[1]);
3980 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3986 matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props UNUSED,
3987 gsl_matrix *subs[], size_t n_subs,
3988 matrix_proto_m_mdd *f)
3990 assert (n_subs == 3);
3991 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
3993 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
3997 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
3998 gsl_matrix *subs[], size_t n_subs,
3999 matrix_proto_m_edd *f)
4001 assert (n_subs == 3);
4002 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
4005 if (!check_constraints (props, subs, n_subs))
4008 double b = to_scalar (subs[1]);
4009 double c = to_scalar (subs[2]);
4010 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4016 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4017 gsl_matrix *subs[], size_t n_subs,
4018 matrix_proto_m_eddd *f)
4020 assert (n_subs == 4);
4021 for (size_t i = 1; i < 4; i++)
4022 if (!check_scalar_arg (props->name, subs, i))
4025 if (!check_constraints (props, subs, n_subs))
4028 double b = to_scalar (subs[1]);
4029 double c = to_scalar (subs[2]);
4030 double d = to_scalar (subs[3]);
4031 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4032 *a = f (*a, b, c, d);
4037 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4038 gsl_matrix *subs[], size_t n_subs,
4039 matrix_proto_m_eed *f)
4041 assert (n_subs == 3);
4042 if (!check_scalar_arg (props->name, subs, 2))
4045 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4046 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4048 msg (ME, _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4049 "%zu×%zu, but %s requires these arguments either to have "
4050 "the same dimensions or for one of them to be a scalar."),
4052 subs[0]->size1, subs[0]->size2,
4053 subs[1]->size1, subs[1]->size2,
4058 if (!check_constraints (props, subs, n_subs))
4061 double c = to_scalar (subs[2]);
4063 if (is_scalar (subs[0]))
4065 double a = to_scalar (subs[0]);
4066 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4072 double b = to_scalar (subs[1]);
4073 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4080 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED,
4081 gsl_matrix *subs[], size_t n_subs,
4082 matrix_proto_m_mm *f)
4084 assert (n_subs == 2);
4085 return f (subs[0], subs[1]);
4089 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4090 gsl_matrix *subs[], size_t n_subs,
4091 matrix_proto_m_v *f)
4093 assert (n_subs == 1);
4094 if (!check_vector_arg (props->name, subs, 0))
4096 gsl_vector v = to_vector (subs[0]);
4101 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED,
4102 gsl_matrix *subs[], size_t n_subs,
4103 matrix_proto_d_m *f)
4105 assert (n_subs == 1);
4107 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4108 gsl_matrix_set (m, 0, 0, f (subs[0]));
4113 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED,
4114 gsl_matrix *subs[], size_t n_subs,
4115 matrix_proto_m_any *f)
4117 return f (subs, n_subs);
4121 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props,
4122 gsl_matrix *subs[], size_t n_subs,
4123 matrix_proto_IDENT *f)
4125 assert (n_subs <= 2);
4128 if (!to_scalar_args (props->name, subs, n_subs, d))
4131 return f (d[0], d[n_subs - 1]);
4135 matrix_expr_evaluate (const struct matrix_expr *e)
4137 if (e->op == MOP_NUMBER)
4139 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4140 gsl_matrix_set (m, 0, 0, e->number);
4143 else if (e->op == MOP_VARIABLE)
4145 const gsl_matrix *src = e->variable->value;
4148 msg (SE, _("Uninitialized variable %s used in expression."),
4153 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4154 gsl_matrix_memcpy (dst, src);
4157 else if (e->op == MOP_EOF)
4159 struct dfm_reader *reader = read_file_open (e->eof);
4160 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4161 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4165 enum { N_LOCAL = 3 };
4166 gsl_matrix *local_subs[N_LOCAL];
4167 gsl_matrix **subs = (e->n_subs < N_LOCAL
4169 : xmalloc (e->n_subs * sizeof *subs));
4171 for (size_t i = 0; i < e->n_subs; i++)
4173 subs[i] = matrix_expr_evaluate (e->subs[i]);
4176 for (size_t j = 0; j < i; j++)
4177 gsl_matrix_free (subs[j]);
4178 if (subs != local_subs)
4184 gsl_matrix *result = NULL;
4187 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4188 case MOP_F_##ENUM: \
4190 static const struct matrix_function_properties props = { \
4192 .constraints = CONSTRAINTS, \
4194 result = matrix_expr_evaluate_##PROTO (&props, subs, e->n_subs, \
4195 matrix_eval_##ENUM); \
4202 gsl_matrix_scale (subs[0], -1.0);
4220 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
4224 result = matrix_expr_evaluate_not (subs[0]);
4228 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
4232 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
4236 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
4240 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
4243 case MOP_PASTE_HORZ:
4244 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
4247 case MOP_PASTE_VERT:
4248 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
4252 result = gsl_matrix_alloc (0, 0);
4256 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
4260 result = matrix_expr_evaluate_vec_all (subs[0]);
4264 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
4268 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
4272 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
4281 for (size_t i = 0; i < e->n_subs; i++)
4282 if (subs[i] != result)
4283 gsl_matrix_free (subs[i]);
4284 if (subs != local_subs)
4290 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4293 gsl_matrix *m = matrix_expr_evaluate (e);
4299 msg (SE, _("Expression for %s must evaluate to scalar, "
4300 "not a %zu×%zu matrix."),
4301 context, m->size1, m->size2);
4302 gsl_matrix_free (m);
4307 gsl_matrix_free (m);
4312 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4316 if (!matrix_expr_evaluate_scalar (e, context, &d))
4320 if (d < LONG_MIN || d > LONG_MAX)
4322 msg (SE, _("Expression for %s is outside the integer range."), context);
4329 struct matrix_lvalue
4331 struct matrix_var *var;
4332 struct matrix_expr *indexes[2];
4337 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4341 for (size_t i = 0; i < lvalue->n_indexes; i++)
4342 matrix_expr_destroy (lvalue->indexes[i]);
4347 static struct matrix_lvalue *
4348 matrix_lvalue_parse (struct matrix_state *s)
4350 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4351 if (!lex_force_id (s->lexer))
4353 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4354 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4358 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
4363 lex_get_n (s->lexer, 2);
4365 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
4367 lvalue->n_indexes++;
4369 if (lex_match (s->lexer, T_COMMA))
4371 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
4373 lvalue->n_indexes++;
4375 if (!lex_force_match (s->lexer, T_RPAREN))
4381 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4387 matrix_lvalue_destroy (lvalue);
4392 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4393 enum index_type index_type, size_t other_size,
4394 struct index_vector *iv)
4399 m = matrix_expr_evaluate (e);
4406 bool ok = matrix_normalize_index_vector (m, size, index_type,
4408 gsl_matrix_free (m);
4413 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4414 struct index_vector *iv, gsl_matrix *sm)
4416 gsl_vector dv = to_vector (lvalue->var->value);
4418 /* Convert source matrix 'sm' to source vector 'sv'. */
4419 if (!is_vector (sm))
4421 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
4422 sm->size1, sm->size2);
4425 gsl_vector sv = to_vector (sm);
4426 if (iv->n != sv.size)
4428 msg (SE, _("Can't assign %zu-element vector "
4429 "to %zu-element subvector."), sv.size, iv->n);
4433 /* Assign elements. */
4434 for (size_t x = 0; x < iv->n; x++)
4435 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4440 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4441 struct index_vector *iv0,
4442 struct index_vector *iv1, gsl_matrix *sm)
4444 gsl_matrix *dm = lvalue->var->value;
4446 /* Convert source matrix 'sm' to source vector 'sv'. */
4447 if (iv0->n != sm->size1)
4449 msg (SE, _("Row index vector for assignment to %s has %zu elements "
4450 "but source matrix has %zu rows."),
4451 lvalue->var->name, iv0->n, sm->size1);
4454 if (iv1->n != sm->size2)
4456 msg (SE, _("Column index vector for assignment to %s has %zu elements "
4457 "but source matrix has %zu columns."),
4458 lvalue->var->name, iv1->n, sm->size2);
4462 /* Assign elements. */
4463 for (size_t y = 0; y < iv0->n; y++)
4465 size_t f0 = iv0->indexes[y];
4466 for (size_t x = 0; x < iv1->n; x++)
4468 size_t f1 = iv1->indexes[x];
4469 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4475 /* Takes ownership of M. */
4477 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4478 struct index_vector *iv0, struct index_vector *iv1,
4481 if (!lvalue->n_indexes)
4483 gsl_matrix_free (lvalue->var->value);
4484 lvalue->var->value = sm;
4489 bool ok = (lvalue->n_indexes == 1
4490 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
4491 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
4492 free (iv0->indexes);
4493 free (iv1->indexes);
4494 gsl_matrix_free (sm);
4500 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4501 struct index_vector *iv0,
4502 struct index_vector *iv1)
4504 *iv0 = INDEX_VECTOR_INIT;
4505 *iv1 = INDEX_VECTOR_INIT;
4506 if (!lvalue->n_indexes)
4509 /* Validate destination matrix exists and has the right shape. */
4510 gsl_matrix *dm = lvalue->var->value;
4513 msg (SE, _("Undefined variable %s."), lvalue->var->name);
4516 else if (dm->size1 == 0 || dm->size2 == 0)
4518 msg (SE, _("Cannot index %zu×%zu matrix."),
4519 dm->size1, dm->size2);
4522 else if (lvalue->n_indexes == 1)
4524 if (!is_vector (dm))
4526 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
4527 dm->size1, dm->size2, lvalue->var->name);
4530 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4531 MAX (dm->size1, dm->size2),
4536 assert (lvalue->n_indexes == 2);
4537 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4538 IV_ROW, dm->size2, iv0))
4541 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4542 IV_COLUMN, dm->size1, iv1))
4544 free (iv0->indexes);
4551 /* Takes ownership of M. */
4553 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
4555 struct index_vector iv0, iv1;
4556 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
4558 gsl_matrix_free (sm);
4562 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
4566 /* Matrix command. */
4570 struct matrix_cmd **commands;
4574 static void matrix_cmds_uninit (struct matrix_cmds *);
4578 enum matrix_cmd_type
4601 struct compute_command
4603 struct matrix_lvalue *lvalue;
4604 struct matrix_expr *rvalue;
4608 struct print_command
4610 struct matrix_expr *expression;
4611 bool use_default_format;
4612 struct fmt_spec format;
4614 int space; /* -1 means NEWPAGE. */
4618 struct string_array literals; /* CLABELS/RLABELS. */
4619 struct matrix_expr *expr; /* CNAMES/RNAMES. */
4625 struct do_if_command
4629 struct matrix_expr *condition;
4630 struct matrix_cmds commands;
4640 struct matrix_var *var;
4641 struct matrix_expr *start;
4642 struct matrix_expr *end;
4643 struct matrix_expr *increment;
4645 /* Loop conditions. */
4646 struct matrix_expr *top_condition;
4647 struct matrix_expr *bottom_condition;
4650 struct matrix_cmds commands;
4654 struct display_command
4656 struct matrix_state *state;
4660 struct release_command
4662 struct matrix_var **vars;
4669 struct matrix_expr *expression;
4670 struct save_file *sf;
4676 struct read_file *rf;
4677 struct matrix_lvalue *dst;
4678 struct matrix_expr *size;
4680 enum fmt_type format;
4687 struct write_command
4689 struct write_file *wf;
4690 struct matrix_expr *expression;
4693 /* If this is nonnull, WRITE uses this format.
4695 If this is NULL, WRITE uses free-field format with as many
4696 digits of precision as needed. */
4697 struct fmt_spec *format;
4706 struct matrix_lvalue *dst;
4707 struct dataset *dataset;
4708 struct file_handle *file;
4710 struct var_syntax *vars;
4712 struct matrix_var *names;
4714 /* Treatment of missing values. */
4719 MGET_ERROR, /* Flag error on command. */
4720 MGET_ACCEPT, /* Accept without change, user-missing only. */
4721 MGET_OMIT, /* Drop this case. */
4722 MGET_RECODE /* Recode to 'substitute'. */
4725 double substitute; /* MGET_RECODE only. */
4731 struct msave_command
4733 struct msave_common *common;
4734 struct matrix_expr *expr;
4735 const char *rowtype;
4736 const struct matrix_expr *factors;
4737 const struct matrix_expr *splits;
4743 struct matrix_state *state;
4744 struct file_handle *file;
4746 struct stringi_set rowtypes;
4750 struct eigen_command
4752 struct matrix_expr *expr;
4753 struct matrix_var *evec;
4754 struct matrix_var *eval;
4758 struct setdiag_command
4760 struct matrix_var *dst;
4761 struct matrix_expr *expr;
4767 struct matrix_expr *expr;
4768 struct matrix_var *u;
4769 struct matrix_var *s;
4770 struct matrix_var *v;
4776 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
4777 static bool matrix_cmd_execute (struct matrix_cmd *);
4778 static void matrix_cmd_destroy (struct matrix_cmd *);
4781 static struct matrix_cmd *
4782 matrix_parse_compute (struct matrix_state *s)
4784 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4785 *cmd = (struct matrix_cmd) {
4786 .type = MCMD_COMPUTE,
4787 .compute = { .lvalue = NULL },
4790 struct compute_command *compute = &cmd->compute;
4791 compute->lvalue = matrix_lvalue_parse (s);
4792 if (!compute->lvalue)
4795 if (!lex_force_match (s->lexer, T_EQUALS))
4798 compute->rvalue = matrix_parse_expr (s);
4799 if (!compute->rvalue)
4805 matrix_cmd_destroy (cmd);
4810 matrix_cmd_execute_compute (struct compute_command *compute)
4812 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
4816 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
4820 print_labels_destroy (struct print_labels *labels)
4824 string_array_destroy (&labels->literals);
4825 matrix_expr_destroy (labels->expr);
4831 parse_literal_print_labels (struct matrix_state *s,
4832 struct print_labels **labelsp)
4834 lex_match (s->lexer, T_EQUALS);
4835 print_labels_destroy (*labelsp);
4836 *labelsp = xzalloc (sizeof **labelsp);
4837 while (lex_token (s->lexer) != T_SLASH
4838 && lex_token (s->lexer) != T_ENDCMD
4839 && lex_token (s->lexer) != T_STOP)
4841 struct string label = DS_EMPTY_INITIALIZER;
4842 while (lex_token (s->lexer) != T_COMMA
4843 && lex_token (s->lexer) != T_SLASH
4844 && lex_token (s->lexer) != T_ENDCMD
4845 && lex_token (s->lexer) != T_STOP)
4847 if (!ds_is_empty (&label))
4848 ds_put_byte (&label, ' ');
4850 if (lex_token (s->lexer) == T_STRING)
4851 ds_put_cstr (&label, lex_tokcstr (s->lexer));
4854 char *rep = lex_next_representation (s->lexer, 0, 0);
4855 ds_put_cstr (&label, rep);
4860 string_array_append_nocopy (&(*labelsp)->literals,
4861 ds_steal_cstr (&label));
4863 if (!lex_match (s->lexer, T_COMMA))
4869 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
4871 lex_match (s->lexer, T_EQUALS);
4872 struct matrix_expr *e = matrix_parse_exp (s);
4876 print_labels_destroy (*labelsp);
4877 *labelsp = xzalloc (sizeof **labelsp);
4878 (*labelsp)->expr = e;
4882 static struct matrix_cmd *
4883 matrix_parse_print (struct matrix_state *s)
4885 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4886 *cmd = (struct matrix_cmd) {
4889 .use_default_format = true,
4893 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
4896 for (size_t i = 0; ; i++)
4898 enum token_type t = lex_next_token (s->lexer, i);
4899 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
4901 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
4903 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
4906 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
4911 cmd->print.expression = matrix_parse_exp (s);
4912 if (!cmd->print.expression)
4916 while (lex_match (s->lexer, T_SLASH))
4918 if (lex_match_id (s->lexer, "FORMAT"))
4920 lex_match (s->lexer, T_EQUALS);
4921 if (!parse_format_specifier (s->lexer, &cmd->print.format))
4923 cmd->print.use_default_format = false;
4925 else if (lex_match_id (s->lexer, "TITLE"))
4927 lex_match (s->lexer, T_EQUALS);
4928 if (!lex_force_string (s->lexer))
4930 free (cmd->print.title);
4931 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
4934 else if (lex_match_id (s->lexer, "SPACE"))
4936 lex_match (s->lexer, T_EQUALS);
4937 if (lex_match_id (s->lexer, "NEWPAGE"))
4938 cmd->print.space = -1;
4939 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
4941 cmd->print.space = lex_integer (s->lexer);
4947 else if (lex_match_id (s->lexer, "RLABELS"))
4948 parse_literal_print_labels (s, &cmd->print.rlabels);
4949 else if (lex_match_id (s->lexer, "CLABELS"))
4950 parse_literal_print_labels (s, &cmd->print.clabels);
4951 else if (lex_match_id (s->lexer, "RNAMES"))
4953 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
4956 else if (lex_match_id (s->lexer, "CNAMES"))
4958 if (!parse_expr_print_labels (s, &cmd->print.clabels))
4963 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
4964 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
4972 matrix_cmd_destroy (cmd);
4977 matrix_is_integer (const gsl_matrix *m)
4979 for (size_t y = 0; y < m->size1; y++)
4980 for (size_t x = 0; x < m->size2; x++)
4982 double d = gsl_matrix_get (m, y, x);
4990 matrix_max_magnitude (const gsl_matrix *m)
4993 for (size_t y = 0; y < m->size1; y++)
4994 for (size_t x = 0; x < m->size2; x++)
4996 double d = fabs (gsl_matrix_get (m, y, x));
5004 format_fits (struct fmt_spec format, double x)
5006 char *s = data_out (&(union value) { .f = x }, NULL,
5007 &format, settings_get_fmt_settings ());
5008 bool fits = *s != '*' && !strchr (s, 'E');
5013 static struct fmt_spec
5014 default_format (const gsl_matrix *m, int *log_scale)
5018 double max = matrix_max_magnitude (m);
5020 if (matrix_is_integer (m))
5021 for (int w = 1; w <= 12; w++)
5023 struct fmt_spec format = { .type = FMT_F, .w = w };
5024 if (format_fits (format, -max))
5028 if (max >= 1e9 || max <= 1e-4)
5031 snprintf (s, sizeof s, "%.1e", max);
5033 const char *e = strchr (s, 'e');
5035 *log_scale = atoi (e + 1);
5038 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5042 trimmed_string (double d)
5044 struct substring s = ss_buffer ((char *) &d, sizeof d);
5045 ss_rtrim (&s, ss_cstr (" "));
5046 return ss_xstrdup (s);
5049 static struct string_array *
5050 print_labels_get (const struct print_labels *labels, size_t n,
5051 const char *prefix, bool truncate)
5056 struct string_array *out = xzalloc (sizeof *out);
5057 if (labels->literals.n)
5058 string_array_clone (out, &labels->literals);
5059 else if (labels->expr)
5061 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5062 if (m && is_vector (m))
5064 gsl_vector v = to_vector (m);
5065 for (size_t i = 0; i < v.size; i++)
5066 string_array_append_nocopy (out, trimmed_string (
5067 gsl_vector_get (&v, i)));
5069 gsl_matrix_free (m);
5075 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5078 string_array_append (out, "");
5081 string_array_delete (out, out->n - 1);
5084 for (size_t i = 0; i < out->n; i++)
5086 char *s = out->strings[i];
5087 s[strnlen (s, 8)] = '\0';
5094 matrix_cmd_print_space (int space)
5097 output_item_submit (page_break_item_create ());
5098 for (int i = 0; i < space; i++)
5099 output_log ("%s", "");
5103 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
5104 struct fmt_spec format, int log_scale)
5106 matrix_cmd_print_space (print->space);
5108 output_log ("%s", print->title);
5110 output_log (" 10 ** %d X", log_scale);
5112 struct string_array *clabels = print_labels_get (print->clabels,
5113 m->size2, "col", true);
5114 if (clabels && format.w < 8)
5116 struct string_array *rlabels = print_labels_get (print->rlabels,
5117 m->size1, "row", true);
5121 struct string line = DS_EMPTY_INITIALIZER;
5123 ds_put_byte_multiple (&line, ' ', 8);
5124 for (size_t x = 0; x < m->size2; x++)
5125 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5126 output_log_nocopy (ds_steal_cstr (&line));
5129 double scale = pow (10.0, log_scale);
5130 bool numeric = fmt_is_numeric (format.type);
5131 for (size_t y = 0; y < m->size1; y++)
5133 struct string line = DS_EMPTY_INITIALIZER;
5135 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5137 for (size_t x = 0; x < m->size2; x++)
5139 double f = gsl_matrix_get (m, y, x);
5141 ? data_out (&(union value) { .f = f / scale}, NULL,
5142 &format, settings_get_fmt_settings ())
5143 : trimmed_string (f));
5144 ds_put_format (&line, " %s", s);
5147 output_log_nocopy (ds_steal_cstr (&line));
5150 string_array_destroy (rlabels);
5152 string_array_destroy (clabels);
5157 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5158 const struct print_labels *print_labels, size_t n,
5159 const char *name, const char *prefix)
5161 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5163 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5164 for (size_t i = 0; i < n; i++)
5165 pivot_category_create_leaf (
5167 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5168 : pivot_value_new_integer (i + 1)));
5170 d->hide_all_labels = true;
5171 string_array_destroy (labels);
5176 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
5177 struct fmt_spec format, int log_scale)
5179 struct pivot_table *table = pivot_table_create__ (
5180 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5183 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5185 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5186 N_("Columns"), "col");
5188 struct pivot_footnote *footnote = NULL;
5191 char *s = xasprintf ("× 10**%d", log_scale);
5192 footnote = pivot_table_create_footnote (
5193 table, pivot_value_new_user_text_nocopy (s));
5196 double scale = pow (10.0, log_scale);
5197 bool numeric = fmt_is_numeric (format.type);
5198 for (size_t y = 0; y < m->size1; y++)
5199 for (size_t x = 0; x < m->size2; x++)
5201 double f = gsl_matrix_get (m, y, x);
5202 struct pivot_value *v;
5205 v = pivot_value_new_number (f / scale);
5206 v->numeric.format = format;
5209 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5211 pivot_value_add_footnote (v, footnote);
5212 pivot_table_put2 (table, y, x, v);
5215 pivot_table_submit (table);
5219 matrix_cmd_execute_print (const struct print_command *print)
5221 if (print->expression)
5223 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5228 struct fmt_spec format = (print->use_default_format
5229 ? default_format (m, &log_scale)
5232 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5233 matrix_cmd_print_text (print, m, format, log_scale);
5235 matrix_cmd_print_tables (print, m, format, log_scale);
5237 gsl_matrix_free (m);
5241 matrix_cmd_print_space (print->space);
5243 output_log ("%s", print->title);
5250 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
5251 const char *command_name,
5252 const char *stop1, const char *stop2)
5254 lex_end_of_command (s->lexer);
5255 lex_discard_rest_of_command (s->lexer);
5257 size_t allocated = 0;
5260 while (lex_token (s->lexer) == T_ENDCMD)
5263 if (lex_at_phrase (s->lexer, stop1)
5264 || (stop2 && lex_at_phrase (s->lexer, stop2)))
5267 if (lex_at_phrase (s->lexer, "END MATRIX"))
5269 msg (SE, _("Premature END MATRIX within %s."), command_name);
5273 struct matrix_cmd *cmd = matrix_parse_command (s);
5277 if (c->n >= allocated)
5278 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
5279 c->commands[c->n++] = cmd;
5284 matrix_parse_do_if_clause (struct matrix_state *s,
5285 struct do_if_command *ifc,
5286 bool parse_condition,
5287 size_t *allocated_clauses)
5289 if (ifc->n_clauses >= *allocated_clauses)
5290 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5291 sizeof *ifc->clauses);
5292 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5293 *c = (struct do_if_clause) { .condition = NULL };
5295 if (parse_condition)
5297 c->condition = matrix_parse_expr (s);
5302 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
5305 static struct matrix_cmd *
5306 matrix_parse_do_if (struct matrix_state *s)
5308 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5309 *cmd = (struct matrix_cmd) {
5311 .do_if = { .n_clauses = 0 }
5314 size_t allocated_clauses = 0;
5317 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
5320 while (lex_match_phrase (s->lexer, "ELSE IF"));
5322 if (lex_match_id (s->lexer, "ELSE")
5323 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
5326 if (!lex_match_phrase (s->lexer, "END IF"))
5331 matrix_cmd_destroy (cmd);
5336 matrix_cmd_execute_do_if (struct do_if_command *cmd)
5338 for (size_t i = 0; i < cmd->n_clauses; i++)
5340 struct do_if_clause *c = &cmd->clauses[i];
5344 if (!matrix_expr_evaluate_scalar (c->condition,
5345 i ? "ELSE IF" : "DO IF",
5350 for (size_t j = 0; j < c->commands.n; j++)
5351 if (!matrix_cmd_execute (c->commands.commands[j]))
5358 static struct matrix_cmd *
5359 matrix_parse_loop (struct matrix_state *s)
5361 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5362 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5364 struct loop_command *loop = &cmd->loop;
5365 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5367 struct substring name = lex_tokss (s->lexer);
5368 loop->var = matrix_var_lookup (s, name);
5370 loop->var = matrix_var_create (s, name);
5375 loop->start = matrix_parse_expr (s);
5376 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5379 loop->end = matrix_parse_expr (s);
5383 if (lex_match (s->lexer, T_BY))
5385 loop->increment = matrix_parse_expr (s);
5386 if (!loop->increment)
5391 if (lex_match_id (s->lexer, "IF"))
5393 loop->top_condition = matrix_parse_expr (s);
5394 if (!loop->top_condition)
5398 bool was_in_loop = s->in_loop;
5400 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
5402 s->in_loop = was_in_loop;
5406 if (!lex_match_phrase (s->lexer, "END LOOP"))
5409 if (lex_match_id (s->lexer, "IF"))
5411 loop->bottom_condition = matrix_parse_expr (s);
5412 if (!loop->bottom_condition)
5419 matrix_cmd_destroy (cmd);
5424 matrix_cmd_execute_loop (struct loop_command *cmd)
5428 long int increment = 1;
5431 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5432 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5434 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5438 if (increment > 0 ? value > end
5439 : increment < 0 ? value < end
5444 int mxloops = settings_get_mxloops ();
5445 for (int i = 0; i < mxloops; i++)
5450 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5452 gsl_matrix_free (cmd->var->value);
5453 cmd->var->value = NULL;
5455 if (!cmd->var->value)
5456 cmd->var->value = gsl_matrix_alloc (1, 1);
5458 gsl_matrix_set (cmd->var->value, 0, 0, value);
5461 if (cmd->top_condition)
5464 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5469 for (size_t j = 0; j < cmd->commands.n; j++)
5470 if (!matrix_cmd_execute (cmd->commands.commands[j]))
5473 if (cmd->bottom_condition)
5476 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5477 "END LOOP IF", &d) || d > 0)
5484 if (increment > 0 ? value > end : value < end)
5490 static struct matrix_cmd *
5491 matrix_parse_break (struct matrix_state *s)
5495 msg (SE, _("BREAK not inside LOOP."));
5499 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5500 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
5504 static struct matrix_cmd *
5505 matrix_parse_display (struct matrix_state *s)
5507 while (lex_token (s->lexer) != T_ENDCMD)
5509 if (!lex_match_id (s->lexer, "DICTIONARY")
5510 && !lex_match_id (s->lexer, "STATUS"))
5512 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5517 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5518 *cmd = (struct matrix_cmd) { .type = MCMD_DISPLAY, .display = { s } };
5523 compare_matrix_var_pointers (const void *a_, const void *b_)
5525 const struct matrix_var *const *ap = a_;
5526 const struct matrix_var *const *bp = b_;
5527 const struct matrix_var *a = *ap;
5528 const struct matrix_var *b = *bp;
5529 return strcmp (a->name, b->name);
5533 matrix_cmd_execute_display (struct display_command *cmd)
5535 const struct matrix_state *s = cmd->state;
5537 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5538 struct pivot_dimension *attr_dimension
5539 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5540 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5541 N_("Rows"), N_("Columns"));
5542 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5544 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5547 const struct matrix_var *var;
5548 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5550 vars[n_vars++] = var;
5551 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5553 struct pivot_dimension *rows = pivot_dimension_create (
5554 table, PIVOT_AXIS_ROW, N_("Variable"));
5555 for (size_t i = 0; i < n_vars; i++)
5557 const struct matrix_var *var = vars[i];
5558 pivot_category_create_leaf (
5559 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
5561 size_t r = var->value->size1;
5562 size_t c = var->value->size2;
5563 double values[] = { r, c, r * c * sizeof (double) / 1024 };
5564 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
5565 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
5568 pivot_table_submit (table);
5571 static struct matrix_cmd *
5572 matrix_parse_release (struct matrix_state *s)
5574 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5575 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
5577 size_t allocated_vars = 0;
5578 while (lex_token (s->lexer) == T_ID)
5580 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
5583 if (cmd->release.n_vars >= allocated_vars)
5584 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
5585 sizeof *cmd->release.vars);
5586 cmd->release.vars[cmd->release.n_vars++] = var;
5589 lex_error (s->lexer, _("Variable name expected."));
5592 if (!lex_match (s->lexer, T_COMMA))
5600 matrix_cmd_execute_release (struct release_command *cmd)
5602 for (size_t i = 0; i < cmd->n_vars; i++)
5604 struct matrix_var *var = cmd->vars[i];
5605 gsl_matrix_free (var->value);
5610 static struct save_file *
5611 save_file_create (struct matrix_state *s, struct file_handle *fh,
5612 struct string_array *variables,
5613 struct matrix_expr *names,
5614 struct stringi_set *strings)
5616 for (size_t i = 0; i < s->n_save_files; i++)
5618 struct save_file *sf = s->save_files[i];
5619 if (fh_equal (sf->file, fh))
5623 string_array_destroy (variables);
5624 matrix_expr_destroy (names);
5625 stringi_set_destroy (strings);
5631 struct save_file *sf = xmalloc (sizeof *sf);
5632 *sf = (struct save_file) {
5634 .dataset = s->dataset,
5635 .variables = *variables,
5637 .strings = STRINGI_SET_INITIALIZER (sf->strings),
5640 stringi_set_swap (&sf->strings, strings);
5641 stringi_set_destroy (strings);
5643 s->save_files = xrealloc (s->save_files,
5644 (s->n_save_files + 1) * sizeof *s->save_files);
5645 s->save_files[s->n_save_files++] = sf;
5649 static struct casewriter *
5650 save_file_open (struct save_file *sf, gsl_matrix *m)
5652 if (sf->writer || sf->error)
5656 size_t n_variables = caseproto_get_n_widths (
5657 casewriter_get_proto (sf->writer));
5658 if (m->size2 != n_variables)
5660 msg (ME, _("The first SAVE to %s within this matrix program "
5661 "had %zu columns, so a %zu×%zu matrix cannot be "
5663 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
5664 n_variables, m->size1, m->size2);
5672 struct dictionary *dict = dict_create (get_default_encoding ());
5674 /* Fill 'names' with user-specified names if there were any, either from
5675 sf->variables or sf->names. */
5676 struct string_array names = { .n = 0 };
5677 if (sf->variables.n)
5678 string_array_clone (&names, &sf->variables);
5681 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
5682 if (nm && is_vector (nm))
5684 gsl_vector nv = to_vector (nm);
5685 for (size_t i = 0; i < nv.size; i++)
5687 char *name = trimmed_string (gsl_vector_get (&nv, i));
5688 if (dict_id_is_valid (dict, name, true))
5689 string_array_append_nocopy (&names, name);
5694 gsl_matrix_free (nm);
5697 struct stringi_set strings;
5698 stringi_set_clone (&strings, &sf->strings);
5700 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
5705 name = names.strings[i];
5708 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
5712 int width = stringi_set_delete (&strings, name) ? 8 : 0;
5713 struct variable *var = dict_create_var (dict, name, width);
5716 msg (ME, _("Duplicate variable name %s in SAVE statement."),
5722 if (!stringi_set_is_empty (&strings))
5724 size_t count = stringi_set_count (&strings);
5725 const char *example = stringi_set_node_get_string (
5726 stringi_set_first (&strings));
5729 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
5730 "variable %s."), example);
5732 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
5733 "unknown variable: %s.",
5734 "The SAVE command STRINGS subcommand specifies %zu "
5735 "unknown variables, including %s.",
5740 stringi_set_destroy (&strings);
5741 string_array_destroy (&names);
5750 if (sf->file == fh_inline_file ())
5751 sf->writer = autopaging_writer_create (dict_get_proto (dict));
5753 sf->writer = any_writer_open (sf->file, dict);
5765 save_file_destroy (struct save_file *sf)
5769 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
5771 dataset_set_dict (sf->dataset, sf->dict);
5772 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
5776 casewriter_destroy (sf->writer);
5777 dict_unref (sf->dict);
5779 fh_unref (sf->file);
5780 string_array_destroy (&sf->variables);
5781 matrix_expr_destroy (sf->names);
5782 stringi_set_destroy (&sf->strings);
5787 static struct matrix_cmd *
5788 matrix_parse_save (struct matrix_state *s)
5790 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5791 *cmd = (struct matrix_cmd) {
5793 .save = { .expression = NULL },
5796 struct file_handle *fh = NULL;
5797 struct save_command *save = &cmd->save;
5799 struct string_array variables = STRING_ARRAY_INITIALIZER;
5800 struct matrix_expr *names = NULL;
5801 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
5803 save->expression = matrix_parse_exp (s);
5804 if (!save->expression)
5807 while (lex_match (s->lexer, T_SLASH))
5809 if (lex_match_id (s->lexer, "OUTFILE"))
5811 lex_match (s->lexer, T_EQUALS);
5814 fh = (lex_match (s->lexer, T_ASTERISK)
5816 : fh_parse (s->lexer, FH_REF_FILE, s->session));
5820 else if (lex_match_id (s->lexer, "VARIABLES"))
5822 lex_match (s->lexer, T_EQUALS);
5826 struct dictionary *d = dict_create (get_default_encoding ());
5827 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
5828 PV_NO_SCRATCH | PV_NO_DUPLICATE);
5833 string_array_clear (&variables);
5834 variables = (struct string_array) {
5840 else if (lex_match_id (s->lexer, "NAMES"))
5842 lex_match (s->lexer, T_EQUALS);
5843 matrix_expr_destroy (names);
5844 names = matrix_parse_exp (s);
5848 else if (lex_match_id (s->lexer, "STRINGS"))
5850 lex_match (s->lexer, T_EQUALS);
5851 while (lex_token (s->lexer) == T_ID)
5853 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
5855 if (!lex_match (s->lexer, T_COMMA))
5861 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
5869 if (s->prev_save_file)
5870 fh = fh_ref (s->prev_save_file);
5873 lex_sbc_missing ("OUTFILE");
5877 fh_unref (s->prev_save_file);
5878 s->prev_save_file = fh_ref (fh);
5880 if (variables.n && names)
5882 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
5883 matrix_expr_destroy (names);
5887 save->sf = save_file_create (s, fh, &variables, names, &strings);
5891 string_array_destroy (&variables);
5892 matrix_expr_destroy (names);
5893 stringi_set_destroy (&strings);
5895 matrix_cmd_destroy (cmd);
5900 matrix_cmd_execute_save (const struct save_command *save)
5902 gsl_matrix *m = matrix_expr_evaluate (save->expression);
5906 struct casewriter *writer = save_file_open (save->sf, m);
5909 gsl_matrix_free (m);
5913 const struct caseproto *proto = casewriter_get_proto (writer);
5914 for (size_t y = 0; y < m->size1; y++)
5916 struct ccase *c = case_create (proto);
5917 for (size_t x = 0; x < m->size2; x++)
5919 double d = gsl_matrix_get (m, y, x);
5920 int width = caseproto_get_width (proto, x);
5921 union value *value = case_data_rw_idx (c, x);
5925 memcpy (value->s, &d, width);
5927 casewriter_write (writer, c);
5929 gsl_matrix_free (m);
5932 static struct matrix_cmd *
5933 matrix_parse_read (struct matrix_state *s)
5935 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5936 *cmd = (struct matrix_cmd) {
5938 .read = { .format = FMT_F },
5941 struct file_handle *fh = NULL;
5942 char *encoding = NULL;
5943 struct read_command *read = &cmd->read;
5944 read->dst = matrix_lvalue_parse (s);
5949 int repetitions = 0;
5950 int record_width = 0;
5951 bool seen_format = false;
5952 while (lex_match (s->lexer, T_SLASH))
5954 if (lex_match_id (s->lexer, "FILE"))
5956 lex_match (s->lexer, T_EQUALS);
5959 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5963 else if (lex_match_id (s->lexer, "ENCODING"))
5965 lex_match (s->lexer, T_EQUALS);
5966 if (!lex_force_string (s->lexer))
5970 encoding = ss_xstrdup (lex_tokss (s->lexer));
5974 else if (lex_match_id (s->lexer, "FIELD"))
5976 lex_match (s->lexer, T_EQUALS);
5978 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5980 read->c1 = lex_integer (s->lexer);
5982 if (!lex_force_match (s->lexer, T_TO)
5983 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
5985 read->c2 = lex_integer (s->lexer) + 1;
5988 record_width = read->c2 - read->c1;
5989 if (lex_match (s->lexer, T_BY))
5991 if (!lex_force_int_range (s->lexer, "BY", 1,
5992 read->c2 - read->c1))
5994 by = lex_integer (s->lexer);
5997 if (record_width % by)
5999 msg (SE, _("BY %d does not evenly divide record width %d."),
6007 else if (lex_match_id (s->lexer, "SIZE"))
6009 lex_match (s->lexer, T_EQUALS);
6010 matrix_expr_destroy (read->size);
6011 read->size = matrix_parse_exp (s);
6015 else if (lex_match_id (s->lexer, "MODE"))
6017 lex_match (s->lexer, T_EQUALS);
6018 if (lex_match_id (s->lexer, "RECTANGULAR"))
6019 read->symmetric = false;
6020 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6021 read->symmetric = true;
6024 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6028 else if (lex_match_id (s->lexer, "REREAD"))
6029 read->reread = true;
6030 else if (lex_match_id (s->lexer, "FORMAT"))
6034 lex_sbc_only_once ("FORMAT");
6039 lex_match (s->lexer, T_EQUALS);
6041 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6044 const char *p = lex_tokcstr (s->lexer);
6045 if (c_isdigit (p[0]))
6047 repetitions = atoi (p);
6048 p += strspn (p, "0123456789");
6049 if (!fmt_from_name (p, &read->format))
6051 lex_error (s->lexer, _("Unknown format %s."), p);
6056 else if (fmt_from_name (p, &read->format))
6060 struct fmt_spec format;
6061 if (!parse_format_specifier (s->lexer, &format))
6063 read->format = format.type;
6069 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6070 "REREAD", "FORMAT");
6077 lex_sbc_missing ("FIELD");
6081 if (!read->dst->n_indexes && !read->size)
6083 msg (SE, _("SIZE is required for reading data into a full matrix "
6084 "(as opposed to a submatrix)."));
6090 if (s->prev_read_file)
6091 fh = fh_ref (s->prev_read_file);
6094 lex_sbc_missing ("FILE");
6098 fh_unref (s->prev_read_file);
6099 s->prev_read_file = fh_ref (fh);
6101 read->rf = read_file_create (s, fh);
6105 free (read->rf->encoding);
6106 read->rf->encoding = encoding;
6110 /* Field width may be specified in multiple ways:
6113 2. The format on FORMAT.
6114 3. The repetition factor on FORMAT.
6116 (2) and (3) are mutually exclusive.
6118 If more than one of these is present, they must agree. If none of them is
6119 present, then free-field format is used.
6121 if (repetitions > record_width)
6123 msg (SE, _("%d repetitions cannot fit in record width %d."),
6124 repetitions, record_width);
6127 int w = (repetitions ? record_width / repetitions
6133 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6134 "which implies field width %d, "
6135 "but BY specifies field width %d."),
6136 repetitions, record_width, w, by);
6138 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6147 matrix_cmd_destroy (cmd);
6153 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6154 struct substring data, size_t y, size_t x,
6155 int first_column, int last_column, char *error)
6157 int line_number = dfm_get_line_number (reader);
6158 struct msg_location *location = xmalloc (sizeof *location);
6159 *location = (struct msg_location) {
6160 .file_name = intern_new (dfm_get_file_name (reader)),
6161 .first_line = line_number,
6162 .last_line = line_number + 1,
6163 .first_column = first_column,
6164 .last_column = last_column,
6166 struct msg *m = xmalloc (sizeof *m);
6168 .category = MSG_C_DATA,
6169 .severity = MSG_S_WARNING,
6170 .location = location,
6171 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
6172 "for matrix row %zu, column %zu: %s"),
6173 (int) data.length, data.string, fmt_name (format),
6174 y + 1, x + 1, error),
6182 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
6183 gsl_matrix *m, struct substring p, size_t y, size_t x,
6184 const char *line_start)
6186 const char *input_encoding = dfm_reader_get_encoding (reader);
6189 if (fmt_is_numeric (read->format))
6192 error = data_in (p, input_encoding, read->format,
6193 settings_get_fmt_settings (), &v, 0, NULL);
6194 if (!error && v.f == SYSMIS)
6195 error = xstrdup (_("Matrix data may not contain missing value."));
6200 uint8_t s[sizeof (double)];
6201 union value v = { .s = s };
6202 error = data_in (p, input_encoding, read->format,
6203 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6204 memcpy (&f, s, sizeof f);
6209 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6210 int c2 = c1 + ss_utf8_count_columns (p) - 1;
6211 parse_error (reader, read->format, p, y, x, c1, c2, error);
6215 gsl_matrix_set (m, y, x, f);
6216 if (read->symmetric && x != y)
6217 gsl_matrix_set (m, x, y, f);
6222 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
6223 struct substring *line, const char **startp)
6225 if (dfm_eof (reader))
6227 msg (SE, _("Unexpected end of file reading matrix data."));
6230 dfm_expand_tabs (reader);
6231 struct substring record = dfm_get_record (reader);
6232 /* XXX need to recode record into UTF-8 */
6233 *startp = record.string;
6234 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6239 matrix_read (struct read_command *read, struct dfm_reader *reader,
6242 for (size_t y = 0; y < m->size1; y++)
6244 size_t nx = read->symmetric ? y + 1 : m->size2;
6246 struct substring line = ss_empty ();
6247 const char *line_start = line.string;
6248 for (size_t x = 0; x < nx; x++)
6255 ss_ltrim (&line, ss_cstr (" ,"));
6256 if (!ss_is_empty (line))
6258 if (!matrix_read_line (read, reader, &line, &line_start))
6260 dfm_forward_record (reader);
6263 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6267 if (!matrix_read_line (read, reader, &line, &line_start))
6269 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6270 int f = x % fields_per_line;
6271 if (f == fields_per_line - 1)
6272 dfm_forward_record (reader);
6274 p = ss_substr (line, read->w * f, read->w);
6277 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6281 dfm_forward_record (reader);
6284 ss_ltrim (&line, ss_cstr (" ,"));
6285 if (!ss_is_empty (line))
6288 msg (SW, _("Trailing garbage on line \"%.*s\""),
6289 (int) line.length, line.string);
6296 matrix_cmd_execute_read (struct read_command *read)
6298 struct index_vector iv0, iv1;
6299 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6302 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6305 gsl_matrix *m = matrix_expr_evaluate (read->size);
6311 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6312 "not a %zu×%zu matrix."), m->size1, m->size2);
6313 gsl_matrix_free (m);
6319 gsl_vector v = to_vector (m);
6323 d[0] = gsl_vector_get (&v, 0);
6326 else if (v.size == 2)
6328 d[0] = gsl_vector_get (&v, 0);
6329 d[1] = gsl_vector_get (&v, 1);
6333 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6334 "not a %zu×%zu matrix."),
6335 m->size1, m->size2),
6336 gsl_matrix_free (m);
6341 gsl_matrix_free (m);
6343 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6345 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
6346 "are outside valid range."),
6357 if (read->dst->n_indexes)
6359 size_t submatrix_size[2];
6360 if (read->dst->n_indexes == 2)
6362 submatrix_size[0] = iv0.n;
6363 submatrix_size[1] = iv1.n;
6365 else if (read->dst->var->value->size1 == 1)
6367 submatrix_size[0] = 1;
6368 submatrix_size[1] = iv0.n;
6372 submatrix_size[0] = iv0.n;
6373 submatrix_size[1] = 1;
6378 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6380 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
6381 "differ from submatrix dimensions %zu×%zu."),
6383 submatrix_size[0], submatrix_size[1]);
6391 size[0] = submatrix_size[0];
6392 size[1] = submatrix_size[1];
6396 struct dfm_reader *reader = read_file_open (read->rf);
6398 dfm_reread_record (reader, 1);
6400 if (read->symmetric && size[0] != size[1])
6402 msg (SE, _("Cannot read non-square %zu×%zu matrix "
6403 "using READ with MODE=SYMMETRIC."),
6409 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6410 matrix_read (read, reader, tmp);
6411 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
6414 static struct matrix_cmd *
6415 matrix_parse_write (struct matrix_state *s)
6417 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6418 *cmd = (struct matrix_cmd) {
6422 struct file_handle *fh = NULL;
6423 char *encoding = NULL;
6424 struct write_command *write = &cmd->write;
6425 write->expression = matrix_parse_exp (s);
6426 if (!write->expression)
6430 int repetitions = 0;
6431 int record_width = 0;
6432 enum fmt_type format = FMT_F;
6433 bool has_format = false;
6434 while (lex_match (s->lexer, T_SLASH))
6436 if (lex_match_id (s->lexer, "OUTFILE"))
6438 lex_match (s->lexer, T_EQUALS);
6441 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6445 else if (lex_match_id (s->lexer, "ENCODING"))
6447 lex_match (s->lexer, T_EQUALS);
6448 if (!lex_force_string (s->lexer))
6452 encoding = ss_xstrdup (lex_tokss (s->lexer));
6456 else if (lex_match_id (s->lexer, "FIELD"))
6458 lex_match (s->lexer, T_EQUALS);
6460 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6462 write->c1 = lex_integer (s->lexer);
6464 if (!lex_force_match (s->lexer, T_TO)
6465 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
6467 write->c2 = lex_integer (s->lexer) + 1;
6470 record_width = write->c2 - write->c1;
6471 if (lex_match (s->lexer, T_BY))
6473 if (!lex_force_int_range (s->lexer, "BY", 1,
6474 write->c2 - write->c1))
6476 by = lex_integer (s->lexer);
6479 if (record_width % by)
6481 msg (SE, _("BY %d does not evenly divide record width %d."),
6489 else if (lex_match_id (s->lexer, "MODE"))
6491 lex_match (s->lexer, T_EQUALS);
6492 if (lex_match_id (s->lexer, "RECTANGULAR"))
6493 write->triangular = false;
6494 else if (lex_match_id (s->lexer, "TRIANGULAR"))
6495 write->triangular = true;
6498 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
6502 else if (lex_match_id (s->lexer, "HOLD"))
6504 else if (lex_match_id (s->lexer, "FORMAT"))
6506 if (has_format || write->format)
6508 lex_sbc_only_once ("FORMAT");
6512 lex_match (s->lexer, T_EQUALS);
6514 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6517 const char *p = lex_tokcstr (s->lexer);
6518 if (c_isdigit (p[0]))
6520 repetitions = atoi (p);
6521 p += strspn (p, "0123456789");
6522 if (!fmt_from_name (p, &format))
6524 lex_error (s->lexer, _("Unknown format %s."), p);
6530 else if (fmt_from_name (p, &format))
6537 struct fmt_spec spec;
6538 if (!parse_format_specifier (s->lexer, &spec))
6540 write->format = xmemdup (&spec, sizeof spec);
6545 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
6553 lex_sbc_missing ("FIELD");
6559 if (s->prev_write_file)
6560 fh = fh_ref (s->prev_write_file);
6563 lex_sbc_missing ("OUTFILE");
6567 fh_unref (s->prev_write_file);
6568 s->prev_write_file = fh_ref (fh);
6570 write->wf = write_file_create (s, fh);
6574 free (write->wf->encoding);
6575 write->wf->encoding = encoding;
6579 /* Field width may be specified in multiple ways:
6582 2. The format on FORMAT.
6583 3. The repetition factor on FORMAT.
6585 (2) and (3) are mutually exclusive.
6587 If more than one of these is present, they must agree. If none of them is
6588 present, then free-field format is used.
6590 if (repetitions > record_width)
6592 msg (SE, _("%d repetitions cannot fit in record width %d."),
6593 repetitions, record_width);
6596 int w = (repetitions ? record_width / repetitions
6597 : write->format ? write->format->w
6602 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6603 "which implies field width %d, "
6604 "but BY specifies field width %d."),
6605 repetitions, record_width, w, by);
6607 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6611 if (w && !write->format)
6613 write->format = xmalloc (sizeof *write->format);
6614 *write->format = (struct fmt_spec) { .type = format, .w = w };
6616 if (!fmt_check_output (write->format))
6620 if (write->format && fmt_var_width (write->format) > sizeof (double))
6622 char s[FMT_STRING_LEN_MAX + 1];
6623 fmt_to_string (write->format, s);
6624 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
6625 s, sizeof (double));
6633 matrix_cmd_destroy (cmd);
6638 matrix_cmd_execute_write (struct write_command *write)
6640 gsl_matrix *m = matrix_expr_evaluate (write->expression);
6644 if (write->triangular && m->size1 != m->size2)
6646 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
6647 "the matrix to be written has dimensions %zu×%zu."),
6648 m->size1, m->size2);
6649 gsl_matrix_free (m);
6653 struct dfm_writer *writer = write_file_open (write->wf);
6654 if (!writer || !m->size1)
6656 gsl_matrix_free (m);
6660 const struct fmt_settings *settings = settings_get_fmt_settings ();
6661 struct u8_line *line = write->wf->held;
6662 for (size_t y = 0; y < m->size1; y++)
6666 line = xmalloc (sizeof *line);
6667 u8_line_init (line);
6669 size_t nx = write->triangular ? y + 1 : m->size2;
6671 for (size_t x = 0; x < nx; x++)
6674 double f = gsl_matrix_get (m, y, x);
6678 if (fmt_is_numeric (write->format->type))
6681 v.s = (uint8_t *) &f;
6682 s = data_out (&v, NULL, write->format, settings);
6686 s = xmalloc (DBL_BUFSIZE_BOUND);
6687 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
6688 >= DBL_BUFSIZE_BOUND)
6691 size_t len = strlen (s);
6692 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
6693 if (width + x0 > write->c2)
6695 dfm_put_record_utf8 (writer, line->s.ss.string,
6697 u8_line_clear (line);
6700 u8_line_put (line, x0, x0 + width, s, len);
6703 x0 += write->format ? write->format->w : width + 1;
6706 if (y + 1 >= m->size1 && write->hold)
6708 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
6709 u8_line_clear (line);
6713 u8_line_destroy (line);
6717 write->wf->held = line;
6719 gsl_matrix_free (m);
6722 static struct matrix_cmd *
6723 matrix_parse_get (struct matrix_state *s)
6725 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6726 *cmd = (struct matrix_cmd) {
6729 .dataset = s->dataset,
6730 .user = { .treatment = MGET_ERROR },
6731 .system = { .treatment = MGET_ERROR },
6735 struct get_command *get = &cmd->get;
6736 get->dst = matrix_lvalue_parse (s);
6740 while (lex_match (s->lexer, T_SLASH))
6742 if (lex_match_id (s->lexer, "FILE"))
6744 lex_match (s->lexer, T_EQUALS);
6746 fh_unref (get->file);
6747 if (lex_match (s->lexer, T_ASTERISK))
6751 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6756 else if (lex_match_id (s->lexer, "ENCODING"))
6758 lex_match (s->lexer, T_EQUALS);
6759 if (!lex_force_string (s->lexer))
6762 free (get->encoding);
6763 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
6767 else if (lex_match_id (s->lexer, "VARIABLES"))
6769 lex_match (s->lexer, T_EQUALS);
6773 lex_sbc_only_once ("VARIABLES");
6777 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
6780 else if (lex_match_id (s->lexer, "NAMES"))
6782 lex_match (s->lexer, T_EQUALS);
6783 if (!lex_force_id (s->lexer))
6786 struct substring name = lex_tokss (s->lexer);
6787 get->names = matrix_var_lookup (s, name);
6789 get->names = matrix_var_create (s, name);
6792 else if (lex_match_id (s->lexer, "MISSING"))
6794 lex_match (s->lexer, T_EQUALS);
6795 if (lex_match_id (s->lexer, "ACCEPT"))
6796 get->user.treatment = MGET_ACCEPT;
6797 else if (lex_match_id (s->lexer, "OMIT"))
6798 get->user.treatment = MGET_OMIT;
6799 else if (lex_is_number (s->lexer))
6801 get->user.treatment = MGET_RECODE;
6802 get->user.substitute = lex_number (s->lexer);
6807 lex_error (s->lexer, NULL);
6811 else if (lex_match_id (s->lexer, "SYSMIS"))
6813 lex_match (s->lexer, T_EQUALS);
6814 if (lex_match_id (s->lexer, "OMIT"))
6815 get->system.treatment = MGET_OMIT;
6816 else if (lex_is_number (s->lexer))
6818 get->system.treatment = MGET_RECODE;
6819 get->system.substitute = lex_number (s->lexer);
6824 lex_error (s->lexer, NULL);
6830 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
6831 "MISSING", "SYSMIS");
6836 if (get->user.treatment != MGET_ACCEPT)
6837 get->system.treatment = MGET_ERROR;
6842 matrix_cmd_destroy (cmd);
6847 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
6848 const struct dictionary *dict)
6850 struct variable **vars;
6855 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
6856 &vars, &n_vars, PV_NUMERIC))
6861 n_vars = dict_get_var_cnt (dict);
6862 vars = xnmalloc (n_vars, sizeof *vars);
6863 for (size_t i = 0; i < n_vars; i++)
6865 struct variable *var = dict_get_var (dict, i);
6866 if (!var_is_numeric (var))
6868 msg (SE, _("GET: Variable %s is not numeric."),
6869 var_get_name (var));
6879 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
6880 for (size_t i = 0; i < n_vars; i++)
6882 char s[sizeof (double)];
6884 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
6885 memcpy (&f, s, sizeof f);
6886 gsl_matrix_set (names, i, 0, f);
6889 gsl_matrix_free (get->names->value);
6890 get->names->value = names;
6894 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
6895 long long int casenum = 1;
6897 for (struct ccase *c = casereader_read (reader); c;
6898 c = casereader_read (reader), casenum++)
6900 if (n_rows >= m->size1)
6902 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
6903 for (size_t y = 0; y < n_rows; y++)
6904 for (size_t x = 0; x < n_vars; x++)
6905 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
6906 gsl_matrix_free (m);
6911 for (size_t x = 0; x < n_vars; x++)
6913 const struct variable *var = vars[x];
6914 double d = case_num (c, var);
6917 if (get->system.treatment == MGET_RECODE)
6918 d = get->system.substitute;
6919 else if (get->system.treatment == MGET_OMIT)
6923 msg (SE, _("GET: Variable %s in case %lld "
6924 "is system-missing."),
6925 var_get_name (var), casenum);
6929 else if (var_is_num_missing (var, d, MV_USER))
6931 if (get->user.treatment == MGET_RECODE)
6932 d = get->user.substitute;
6933 else if (get->user.treatment == MGET_OMIT)
6935 else if (get->user.treatment != MGET_ACCEPT)
6937 msg (SE, _("GET: Variable %s in case %lld has user-missing "
6939 var_get_name (var), casenum, d);
6943 gsl_matrix_set (m, n_rows, x, d);
6954 matrix_lvalue_evaluate_and_assign (get->dst, m);
6957 gsl_matrix_free (m);
6962 matrix_open_casereader (const char *command_name,
6963 struct file_handle *file, const char *encoding,
6964 struct dataset *dataset,
6965 struct casereader **readerp, struct dictionary **dictp)
6969 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
6970 return *readerp != NULL;
6974 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
6976 msg (ME, _("The %s command cannot read an empty active file."),
6980 *readerp = proc_open (dataset);
6981 *dictp = dict_ref (dataset_dict (dataset));
6987 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
6988 struct casereader *reader, struct dictionary *dict)
6991 casereader_destroy (reader);
6993 proc_commit (dataset);
6997 matrix_cmd_execute_get (struct get_command *get)
6999 struct casereader *r;
7000 struct dictionary *d;
7001 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
7004 matrix_cmd_execute_get__ (get, r, d);
7005 matrix_close_casereader (get->file, get->dataset, r, d);
7010 match_rowtype (struct lexer *lexer)
7012 static const char *rowtypes[] = {
7013 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7015 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7017 for (size_t i = 0; i < n_rowtypes; i++)
7018 if (lex_match_id (lexer, rowtypes[i]))
7021 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7026 parse_var_names (struct lexer *lexer, struct string_array *sa)
7028 lex_match (lexer, T_EQUALS);
7030 string_array_clear (sa);
7032 struct dictionary *dict = dict_create (get_default_encoding ());
7035 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7036 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7041 for (size_t i = 0; i < n_names; i++)
7042 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7043 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7045 msg (SE, _("Variable name %s is reserved."), names[i]);
7046 for (size_t j = 0; j < n_names; j++)
7052 string_array_clear (sa);
7053 sa->strings = names;
7054 sa->n = sa->allocated = n_names;
7060 msave_common_destroy (struct msave_common *common)
7064 fh_unref (common->outfile);
7065 string_array_destroy (&common->variables);
7066 string_array_destroy (&common->fnames);
7067 string_array_destroy (&common->snames);
7069 for (size_t i = 0; i < common->n_factors; i++)
7070 matrix_expr_destroy (common->factors[i]);
7071 free (common->factors);
7073 for (size_t i = 0; i < common->n_splits; i++)
7074 matrix_expr_destroy (common->splits[i]);
7075 free (common->splits);
7077 dict_unref (common->dict);
7078 casewriter_destroy (common->writer);
7085 compare_variables (const char *keyword,
7086 const struct string_array *new,
7087 const struct string_array *old)
7093 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7094 "on the first MSAVE within MATRIX."), keyword);
7097 else if (!string_array_equal_case (old, new))
7099 msg (SE, _("%s must specify the same variables each time within "
7100 "a given MATRIX."), keyword);
7107 static struct matrix_cmd *
7108 matrix_parse_msave (struct matrix_state *s)
7110 struct msave_common *common = xmalloc (sizeof *common);
7111 *common = (struct msave_common) { .outfile = NULL };
7113 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7114 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7116 struct matrix_expr *splits = NULL;
7117 struct matrix_expr *factors = NULL;
7119 struct msave_command *msave = &cmd->msave;
7120 msave->expr = matrix_parse_exp (s);
7124 while (lex_match (s->lexer, T_SLASH))
7126 if (lex_match_id (s->lexer, "TYPE"))
7128 lex_match (s->lexer, T_EQUALS);
7130 msave->rowtype = match_rowtype (s->lexer);
7131 if (!msave->rowtype)
7134 else if (lex_match_id (s->lexer, "OUTFILE"))
7136 lex_match (s->lexer, T_EQUALS);
7138 fh_unref (common->outfile);
7139 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7140 if (!common->outfile)
7143 else if (lex_match_id (s->lexer, "VARIABLES"))
7145 if (!parse_var_names (s->lexer, &common->variables))
7148 else if (lex_match_id (s->lexer, "FNAMES"))
7150 if (!parse_var_names (s->lexer, &common->fnames))
7153 else if (lex_match_id (s->lexer, "SNAMES"))
7155 if (!parse_var_names (s->lexer, &common->snames))
7158 else if (lex_match_id (s->lexer, "SPLIT"))
7160 lex_match (s->lexer, T_EQUALS);
7162 matrix_expr_destroy (splits);
7163 splits = matrix_parse_exp (s);
7167 else if (lex_match_id (s->lexer, "FACTOR"))
7169 lex_match (s->lexer, T_EQUALS);
7171 matrix_expr_destroy (factors);
7172 factors = matrix_parse_exp (s);
7178 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7179 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7183 if (!msave->rowtype)
7185 lex_sbc_missing ("TYPE");
7191 if (common->fnames.n && !factors)
7193 msg (SE, _("FNAMES requires FACTOR."));
7196 if (common->snames.n && !splits)
7198 msg (SE, _("SNAMES requires SPLIT."));
7201 if (!common->outfile)
7203 lex_sbc_missing ("OUTFILE");
7210 if (common->outfile)
7212 if (!fh_equal (common->outfile, s->common->outfile))
7214 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7215 "within a single MATRIX command."));
7219 if (!compare_variables ("VARIABLES",
7220 &common->variables, &s->common->variables)
7221 || !compare_variables ("FNAMES", &common->fnames, &s->common->fnames)
7222 || !compare_variables ("SNAMES", &common->snames, &s->common->snames))
7224 msave_common_destroy (common);
7226 msave->common = s->common;
7228 struct msave_common *c = s->common;
7231 if (c->n_factors >= c->allocated_factors)
7232 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7233 sizeof *c->factors);
7234 c->factors[c->n_factors++] = factors;
7236 if (c->n_factors > 0)
7237 msave->factors = c->factors[c->n_factors - 1];
7241 if (c->n_splits >= c->allocated_splits)
7242 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7244 c->splits[c->n_splits++] = splits;
7246 if (c->n_splits > 0)
7247 msave->splits = c->splits[c->n_splits - 1];
7252 matrix_expr_destroy (splits);
7253 matrix_expr_destroy (factors);
7254 msave_common_destroy (common);
7255 matrix_cmd_destroy (cmd);
7260 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7262 gsl_matrix *m = matrix_expr_evaluate (e);
7268 msg (SE, _("%s expression must evaluate to vector, "
7269 "not a %zu×%zu matrix."),
7270 name, m->size1, m->size2);
7271 gsl_matrix_free (m);
7275 return matrix_to_vector (m);
7279 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7281 for (size_t i = 0; i < vars->n; i++)
7282 if (!dict_create_var (d, vars->strings[i], 0))
7283 return vars->strings[i];
7287 static struct dictionary *
7288 msave_create_dict (const struct msave_common *common)
7290 struct dictionary *dict = dict_create (get_default_encoding ());
7292 const char *dup_split = msave_add_vars (dict, &common->snames);
7295 /* Should not be possible because the parser ensures that the names are
7300 dict_create_var_assert (dict, "ROWTYPE_", 8);
7302 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7305 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
7309 dict_create_var_assert (dict, "VARNAME_", 8);
7311 const char *dup_var = msave_add_vars (dict, &common->variables);
7314 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
7326 matrix_cmd_execute_msave (struct msave_command *msave)
7328 struct msave_common *common = msave->common;
7329 gsl_matrix *m = NULL;
7330 gsl_vector *factors = NULL;
7331 gsl_vector *splits = NULL;
7333 m = matrix_expr_evaluate (msave->expr);
7337 if (!common->variables.n)
7338 for (size_t i = 0; i < m->size2; i++)
7339 string_array_append_nocopy (&common->variables,
7340 xasprintf ("COL%zu", i + 1));
7341 else if (m->size2 != common->variables.n)
7344 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7345 m->size2, common->variables.n);
7351 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7355 if (!common->fnames.n)
7356 for (size_t i = 0; i < factors->size; i++)
7357 string_array_append_nocopy (&common->fnames,
7358 xasprintf ("FAC%zu", i + 1));
7359 else if (factors->size != common->fnames.n)
7361 msg (SE, _("There are %zu factor variables, "
7362 "but %zu factor values were supplied."),
7363 common->fnames.n, factors->size);
7370 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7374 if (!common->snames.n)
7375 for (size_t i = 0; i < splits->size; i++)
7376 string_array_append_nocopy (&common->snames,
7377 xasprintf ("SPL%zu", i + 1));
7378 else if (splits->size != common->snames.n)
7380 msg (SE, _("There are %zu split variables, "
7381 "but %zu split values were supplied."),
7382 common->snames.n, splits->size);
7387 if (!common->writer)
7389 struct dictionary *dict = msave_create_dict (common);
7393 common->writer = any_writer_open (common->outfile, dict);
7394 if (!common->writer)
7400 common->dict = dict;
7403 bool matrix = (!strcmp (msave->rowtype, "COV")
7404 || !strcmp (msave->rowtype, "CORR"));
7405 for (size_t y = 0; y < m->size1; y++)
7407 struct ccase *c = case_create (dict_get_proto (common->dict));
7410 /* Split variables */
7412 for (size_t i = 0; i < splits->size; i++)
7413 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
7416 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7417 msave->rowtype, ' ');
7421 for (size_t i = 0; i < factors->size; i++)
7422 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
7425 const char *varname_ = (matrix && y < common->variables.n
7426 ? common->variables.strings[y]
7428 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7431 /* Continuous variables. */
7432 for (size_t x = 0; x < m->size2; x++)
7433 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
7434 casewriter_write (common->writer, c);
7438 gsl_matrix_free (m);
7439 gsl_vector_free (factors);
7440 gsl_vector_free (splits);
7443 static struct matrix_cmd *
7444 matrix_parse_mget (struct matrix_state *s)
7446 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7447 *cmd = (struct matrix_cmd) {
7451 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
7455 struct mget_command *mget = &cmd->mget;
7457 lex_match (s->lexer, T_SLASH);
7458 while (lex_token (s->lexer) != T_ENDCMD)
7460 if (lex_match_id (s->lexer, "FILE"))
7462 lex_match (s->lexer, T_EQUALS);
7464 fh_unref (mget->file);
7465 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7469 else if (lex_match_id (s->lexer, "ENCODING"))
7471 lex_match (s->lexer, T_EQUALS);
7472 if (!lex_force_string (s->lexer))
7475 free (mget->encoding);
7476 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
7480 else if (lex_match_id (s->lexer, "TYPE"))
7482 lex_match (s->lexer, T_EQUALS);
7483 while (lex_token (s->lexer) != T_SLASH
7484 && lex_token (s->lexer) != T_ENDCMD)
7486 const char *rowtype = match_rowtype (s->lexer);
7490 stringi_set_insert (&mget->rowtypes, rowtype);
7495 lex_error_expecting (s->lexer, "FILE", "TYPE");
7498 lex_match (s->lexer, T_SLASH);
7503 matrix_cmd_destroy (cmd);
7507 static const struct variable *
7508 get_a8_var (const struct dictionary *d, const char *name)
7510 const struct variable *v = dict_lookup_var (d, name);
7513 msg (SE, _("Matrix data file lacks %s variable."), name);
7516 if (var_get_width (v) != 8)
7518 msg (SE, _("%s variable in matrix data file must be "
7519 "8-byte string, but it has width %d."),
7520 name, var_get_width (v));
7527 var_changed (const struct ccase *ca, const struct ccase *cb,
7528 const struct variable *var)
7531 ? !value_equal (case_data (ca, var), case_data (cb, var),
7532 var_get_width (var))
7537 vars_changed (const struct ccase *ca, const struct ccase *cb,
7538 const struct dictionary *d,
7539 size_t first_var, size_t n_vars)
7541 for (size_t i = 0; i < n_vars; i++)
7543 const struct variable *v = dict_get_var (d, first_var + i);
7544 if (var_changed (ca, cb, v))
7551 vars_all_missing (const struct ccase *c, const struct dictionary *d,
7552 size_t first_var, size_t n_vars)
7554 for (size_t i = 0; i < n_vars; i++)
7555 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
7561 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
7562 const struct dictionary *d,
7563 const struct variable *rowtype_var,
7564 const struct stringi_set *accepted_rowtypes,
7565 struct matrix_state *s,
7566 size_t ss, size_t sn, size_t si,
7567 size_t fs, size_t fn, size_t fi,
7568 size_t cs, size_t cn,
7569 struct pivot_table *pt,
7570 struct pivot_dimension *var_dimension)
7575 /* Is this a matrix for pooled data, either where there are no factor
7576 variables or the factor variables are missing? */
7577 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
7579 struct substring rowtype = case_ss (rows[0], rowtype_var);
7580 ss_rtrim (&rowtype, ss_cstr (" "));
7581 if (!stringi_set_is_empty (accepted_rowtypes)
7582 && !stringi_set_contains_len (accepted_rowtypes,
7583 rowtype.string, rowtype.length))
7586 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
7587 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
7588 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
7589 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
7590 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
7591 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
7595 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
7596 (int) rowtype.length, rowtype.string);
7600 struct string name = DS_EMPTY_INITIALIZER;
7601 ds_put_cstr (&name, prefix);
7603 ds_put_format (&name, "F%zu", fi);
7605 ds_put_format (&name, "S%zu", si);
7607 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
7609 mv = matrix_var_create (s, ds_ss (&name));
7612 msg (SW, _("Matrix data file contains variable with existing name %s."),
7614 goto exit_free_name;
7617 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
7618 size_t n_missing = 0;
7619 for (size_t y = 0; y < n_rows; y++)
7621 for (size_t x = 0; x < cn; x++)
7623 struct variable *var = dict_get_var (d, cs + x);
7624 double value = case_num (rows[y], var);
7625 if (var_is_num_missing (var, value, MV_ANY))
7630 gsl_matrix_set (m, y, x, value);
7634 int var_index = pivot_category_create_leaf (
7635 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
7636 double values[] = { n_rows, cn };
7637 for (size_t j = 0; j < sn; j++)
7639 struct variable *var = dict_get_var (d, ss + j);
7640 const union value *value = case_data (rows[0], var);
7641 pivot_table_put2 (pt, j, var_index,
7642 pivot_value_new_var_value (var, value));
7644 for (size_t j = 0; j < fn; j++)
7646 struct variable *var = dict_get_var (d, fs + j);
7647 const union value sysmis = { .f = SYSMIS };
7648 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
7649 pivot_table_put2 (pt, j + sn, var_index,
7650 pivot_value_new_var_value (var, value));
7652 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
7653 pivot_table_put2 (pt, j + sn + fn, var_index,
7654 pivot_value_new_integer (values[j]));
7657 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
7658 "value, which was treated as zero.",
7659 "Matrix data file variable %s contains %zu missing "
7660 "values, which were treated as zero.", n_missing),
7661 ds_cstr (&name), n_missing);
7668 for (size_t y = 0; y < n_rows; y++)
7669 case_unref (rows[y]);
7673 matrix_cmd_execute_mget__ (struct mget_command *mget,
7674 struct casereader *r, const struct dictionary *d)
7676 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
7677 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
7678 if (!rowtype_ || !varname_)
7681 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
7683 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
7686 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
7688 msg (SE, _("Matrix data file contains no continuous variables."));
7692 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
7694 const struct variable *v = dict_get_var (d, i);
7695 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
7698 _("Matrix data file contains unexpected string variable %s."),
7704 /* SPLIT variables. */
7706 size_t sn = var_get_dict_index (rowtype_);
7707 struct ccase *sc = NULL;
7710 /* FACTOR variables. */
7711 size_t fs = var_get_dict_index (rowtype_) + 1;
7712 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
7713 struct ccase *fc = NULL;
7716 /* Continuous variables. */
7717 size_t cs = var_get_dict_index (varname_) + 1;
7718 size_t cn = dict_get_var_cnt (d) - cs;
7719 struct ccase *cc = NULL;
7722 struct pivot_table *pt = pivot_table_create (
7723 N_("Matrix Variables Created by MGET"));
7724 struct pivot_dimension *attr_dimension = pivot_dimension_create (
7725 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
7726 struct pivot_dimension *var_dimension = pivot_dimension_create (
7727 pt, PIVOT_AXIS_ROW, N_("Variable"));
7730 struct pivot_category *splits = pivot_category_create_group (
7731 attr_dimension->root, N_("Split Values"));
7732 for (size_t i = 0; i < sn; i++)
7733 pivot_category_create_leaf (splits, pivot_value_new_variable (
7734 dict_get_var (d, ss + i)));
7738 struct pivot_category *factors = pivot_category_create_group (
7739 attr_dimension->root, N_("Factors"));
7740 for (size_t i = 0; i < fn; i++)
7741 pivot_category_create_leaf (factors, pivot_value_new_variable (
7742 dict_get_var (d, fs + i)));
7744 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
7745 N_("Rows"), N_("Columns"));
7748 struct ccase **rows = NULL;
7749 size_t allocated_rows = 0;
7753 while ((c = casereader_read (r)) != NULL)
7755 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
7765 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
7766 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
7767 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
7770 if (change != NOTHING_CHANGED)
7772 matrix_mget_commit_var (rows, n_rows, d,
7773 rowtype_, &mget->rowtypes,
7784 if (n_rows >= allocated_rows)
7785 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
7788 if (change == SPLITS_CHANGED)
7794 /* Reset the factor number, if there are factors. */
7798 if (row_has_factors)
7804 else if (change == FACTORS_CHANGED)
7806 if (row_has_factors)
7812 matrix_mget_commit_var (rows, n_rows, d,
7813 rowtype_, &mget->rowtypes,
7825 if (var_dimension->n_leaves)
7826 pivot_table_submit (pt);
7828 pivot_table_unref (pt);
7832 matrix_cmd_execute_mget (struct mget_command *mget)
7834 struct casereader *r;
7835 struct dictionary *d;
7836 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
7837 mget->state->dataset, &r, &d))
7839 matrix_cmd_execute_mget__ (mget, r, d);
7840 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
7845 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
7847 if (!lex_force_id (s->lexer))
7850 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
7852 *varp = matrix_var_create (s, lex_tokss (s->lexer));
7857 static struct matrix_cmd *
7858 matrix_parse_eigen (struct matrix_state *s)
7860 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7861 *cmd = (struct matrix_cmd) {
7863 .eigen = { .expr = NULL }
7866 struct eigen_command *eigen = &cmd->eigen;
7867 if (!lex_force_match (s->lexer, T_LPAREN))
7869 eigen->expr = matrix_parse_expr (s);
7871 || !lex_force_match (s->lexer, T_COMMA)
7872 || !matrix_parse_dst_var (s, &eigen->evec)
7873 || !lex_force_match (s->lexer, T_COMMA)
7874 || !matrix_parse_dst_var (s, &eigen->eval)
7875 || !lex_force_match (s->lexer, T_RPAREN))
7881 matrix_cmd_destroy (cmd);
7886 matrix_cmd_execute_eigen (struct eigen_command *eigen)
7888 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
7891 if (!is_symmetric (A))
7893 msg (SE, _("Argument of EIGEN must be symmetric."));
7894 gsl_matrix_free (A);
7898 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
7899 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
7900 gsl_vector v_eval = to_vector (eval);
7901 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
7902 gsl_eigen_symmv (A, &v_eval, evec, w);
7903 gsl_eigen_symmv_free (w);
7905 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
7907 gsl_matrix_free (A);
7909 gsl_matrix_free (eigen->eval->value);
7910 eigen->eval->value = eval;
7912 gsl_matrix_free (eigen->evec->value);
7913 eigen->evec->value = evec;
7916 static struct matrix_cmd *
7917 matrix_parse_setdiag (struct matrix_state *s)
7919 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7920 *cmd = (struct matrix_cmd) {
7921 .type = MCMD_SETDIAG,
7922 .setdiag = { .dst = NULL }
7925 struct setdiag_command *setdiag = &cmd->setdiag;
7926 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
7929 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
7932 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
7937 if (!lex_force_match (s->lexer, T_COMMA))
7940 setdiag->expr = matrix_parse_expr (s);
7944 if (!lex_force_match (s->lexer, T_RPAREN))
7950 matrix_cmd_destroy (cmd);
7955 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
7957 gsl_matrix *dst = setdiag->dst->value;
7960 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
7961 setdiag->dst->name);
7965 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
7969 size_t n = MIN (dst->size1, dst->size2);
7970 if (is_scalar (src))
7972 double d = to_scalar (src);
7973 for (size_t i = 0; i < n; i++)
7974 gsl_matrix_set (dst, i, i, d);
7976 else if (is_vector (src))
7978 gsl_vector v = to_vector (src);
7979 for (size_t i = 0; i < n && i < v.size; i++)
7980 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
7983 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
7984 "not a %zu×%zu matrix."),
7985 src->size1, src->size2);
7986 gsl_matrix_free (src);
7989 static struct matrix_cmd *
7990 matrix_parse_svd (struct matrix_state *s)
7992 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7993 *cmd = (struct matrix_cmd) {
7995 .svd = { .expr = NULL }
7998 struct svd_command *svd = &cmd->svd;
7999 if (!lex_force_match (s->lexer, T_LPAREN))
8001 svd->expr = matrix_parse_expr (s);
8003 || !lex_force_match (s->lexer, T_COMMA)
8004 || !matrix_parse_dst_var (s, &svd->u)
8005 || !lex_force_match (s->lexer, T_COMMA)
8006 || !matrix_parse_dst_var (s, &svd->s)
8007 || !lex_force_match (s->lexer, T_COMMA)
8008 || !matrix_parse_dst_var (s, &svd->v)
8009 || !lex_force_match (s->lexer, T_RPAREN))
8015 matrix_cmd_destroy (cmd);
8020 matrix_cmd_execute_svd (struct svd_command *svd)
8022 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8026 if (m->size1 >= m->size2)
8029 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8030 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8031 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8032 gsl_vector *work = gsl_vector_alloc (A->size2);
8033 gsl_linalg_SV_decomp (A, V, &Sv, work);
8034 gsl_vector_free (work);
8036 matrix_var_set (svd->u, A);
8037 matrix_var_set (svd->s, S);
8038 matrix_var_set (svd->v, V);
8042 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8043 gsl_matrix_transpose_memcpy (At, m);
8044 gsl_matrix_free (m);
8046 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8047 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8048 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8049 gsl_vector *work = gsl_vector_alloc (At->size2);
8050 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8051 gsl_vector_free (work);
8053 matrix_var_set (svd->v, At);
8054 matrix_var_set (svd->s, St);
8055 matrix_var_set (svd->u, Vt);
8060 matrix_cmd_execute (struct matrix_cmd *cmd)
8065 matrix_cmd_execute_compute (&cmd->compute);
8069 matrix_cmd_execute_print (&cmd->print);
8073 return matrix_cmd_execute_do_if (&cmd->do_if);
8076 matrix_cmd_execute_loop (&cmd->loop);
8083 matrix_cmd_execute_display (&cmd->display);
8087 matrix_cmd_execute_release (&cmd->release);
8091 matrix_cmd_execute_save (&cmd->save);
8095 matrix_cmd_execute_read (&cmd->read);
8099 matrix_cmd_execute_write (&cmd->write);
8103 matrix_cmd_execute_get (&cmd->get);
8107 matrix_cmd_execute_msave (&cmd->msave);
8111 matrix_cmd_execute_mget (&cmd->mget);
8115 matrix_cmd_execute_eigen (&cmd->eigen);
8119 matrix_cmd_execute_setdiag (&cmd->setdiag);
8123 matrix_cmd_execute_svd (&cmd->svd);
8131 matrix_cmds_uninit (struct matrix_cmds *cmds)
8133 for (size_t i = 0; i < cmds->n; i++)
8134 matrix_cmd_destroy (cmds->commands[i]);
8135 free (cmds->commands);
8139 matrix_cmd_destroy (struct matrix_cmd *cmd)
8147 matrix_lvalue_destroy (cmd->compute.lvalue);
8148 matrix_expr_destroy (cmd->compute.rvalue);
8152 matrix_expr_destroy (cmd->print.expression);
8153 free (cmd->print.title);
8154 print_labels_destroy (cmd->print.rlabels);
8155 print_labels_destroy (cmd->print.clabels);
8159 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8161 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8162 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
8164 free (cmd->do_if.clauses);
8168 matrix_expr_destroy (cmd->loop.start);
8169 matrix_expr_destroy (cmd->loop.end);
8170 matrix_expr_destroy (cmd->loop.increment);
8171 matrix_expr_destroy (cmd->loop.top_condition);
8172 matrix_expr_destroy (cmd->loop.bottom_condition);
8173 matrix_cmds_uninit (&cmd->loop.commands);
8183 free (cmd->release.vars);
8187 matrix_expr_destroy (cmd->save.expression);
8191 matrix_lvalue_destroy (cmd->read.dst);
8192 matrix_expr_destroy (cmd->read.size);
8196 matrix_expr_destroy (cmd->write.expression);
8197 free (cmd->write.format);
8201 matrix_lvalue_destroy (cmd->get.dst);
8202 fh_unref (cmd->get.file);
8203 free (cmd->get.encoding);
8204 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8208 matrix_expr_destroy (cmd->msave.expr);
8212 fh_unref (cmd->mget.file);
8213 stringi_set_destroy (&cmd->mget.rowtypes);
8217 matrix_expr_destroy (cmd->eigen.expr);
8221 matrix_expr_destroy (cmd->setdiag.expr);
8225 matrix_expr_destroy (cmd->svd.expr);
8231 struct matrix_command_name
8234 struct matrix_cmd *(*parse) (struct matrix_state *);
8237 static const struct matrix_command_name *
8238 matrix_parse_command_name (struct lexer *lexer)
8240 static const struct matrix_command_name commands[] = {
8241 { "COMPUTE", matrix_parse_compute },
8242 { "CALL EIGEN", matrix_parse_eigen },
8243 { "CALL SETDIAG", matrix_parse_setdiag },
8244 { "CALL SVD", matrix_parse_svd },
8245 { "PRINT", matrix_parse_print },
8246 { "DO IF", matrix_parse_do_if },
8247 { "LOOP", matrix_parse_loop },
8248 { "BREAK", matrix_parse_break },
8249 { "READ", matrix_parse_read },
8250 { "WRITE", matrix_parse_write },
8251 { "GET", matrix_parse_get },
8252 { "SAVE", matrix_parse_save },
8253 { "MGET", matrix_parse_mget },
8254 { "MSAVE", matrix_parse_msave },
8255 { "DISPLAY", matrix_parse_display },
8256 { "RELEASE", matrix_parse_release },
8258 static size_t n = sizeof commands / sizeof *commands;
8260 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8261 if (lex_match_phrase (lexer, c->name))
8266 static struct matrix_cmd *
8267 matrix_parse_command (struct matrix_state *s)
8269 size_t nesting_level = SIZE_MAX;
8271 struct matrix_cmd *c = NULL;
8272 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
8274 lex_error (s->lexer, _("Unknown matrix command."));
8275 else if (!cmd->parse)
8276 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8280 nesting_level = output_open_group (
8281 group_item_create_nocopy (utf8_to_title (cmd->name),
8282 utf8_to_title (cmd->name)));
8287 lex_end_of_command (s->lexer);
8288 lex_discard_rest_of_command (s->lexer);
8289 if (nesting_level != SIZE_MAX)
8290 output_close_groups (nesting_level);
8296 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8298 if (!lex_force_match (lexer, T_ENDCMD))
8301 struct matrix_state state = {
8303 .session = dataset_session (ds),
8305 .vars = HMAP_INITIALIZER (state.vars),
8310 while (lex_match (lexer, T_ENDCMD))
8312 if (lex_token (lexer) == T_STOP)
8314 msg (SE, _("Unexpected end of input expecting matrix command."));
8318 if (lex_match_phrase (lexer, "END MATRIX"))
8321 struct matrix_cmd *c = matrix_parse_command (&state);
8324 matrix_cmd_execute (c);
8325 matrix_cmd_destroy (c);
8329 struct matrix_var *var, *next;
8330 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8333 gsl_matrix_free (var->value);
8334 hmap_delete (&state.vars, &var->hmap_node);
8337 hmap_destroy (&state.vars);
8338 msave_common_destroy (state.common);
8339 fh_unref (state.prev_read_file);
8340 for (size_t i = 0; i < state.n_read_files; i++)
8341 read_file_destroy (state.read_files[i]);
8342 free (state.read_files);
8343 fh_unref (state.prev_write_file);
8344 for (size_t i = 0; i < state.n_write_files; i++)
8345 write_file_destroy (state.write_files[i]);
8346 free (state.write_files);
8347 fh_unref (state.prev_save_file);
8348 for (size_t i = 0; i < state.n_save_files; i++)
8349 save_file_destroy (state.save_files[i]);
8350 free (state.save_files);