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", "a!=0".
211 - "a<b", "a>b", "a<=b", "a>=b", "b!=0".
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_me, 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, "ai>=0 bi>=0") \
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, "ai>=0 bi>=0") \
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, "b!=0") \
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_me_MIN_ARGS = 1, m_me_MAX_ARGS = 1 };
394 enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 };
395 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
396 enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 };
397 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
398 enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
399 enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
400 enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
401 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
402 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
403 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
404 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
405 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
407 typedef double matrix_proto_d_none (void);
408 typedef double matrix_proto_d_d (double);
409 typedef double matrix_proto_d_dd (double, double);
410 typedef double matrix_proto_d_dd (double, double);
411 typedef double matrix_proto_d_ddd (double, double, double);
412 typedef gsl_matrix *matrix_proto_m_d (double);
413 typedef gsl_matrix *matrix_proto_m_dd (double, double);
414 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
415 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
416 typedef gsl_matrix *matrix_proto_m_me (gsl_matrix *, const struct matrix_expr *);
417 typedef double matrix_proto_m_e (double);
418 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
419 typedef double matrix_proto_m_ed (double, double);
420 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
421 typedef double matrix_proto_m_edd (double, double, double);
422 typedef double matrix_proto_m_eddd (double, double, double, double);
423 typedef double matrix_proto_m_eed (double, double, double);
424 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
425 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
426 typedef double matrix_proto_d_m (gsl_matrix *);
427 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
428 typedef gsl_matrix *matrix_proto_IDENT (double, double);
430 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
431 static matrix_proto_##PROTO matrix_eval_##ENUM;
435 /* Matrix expressions. */
442 #define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
446 /* Elementwise and scalar arithmetic. */
447 MOP_NEGATE, /* unary - */
448 MOP_ADD_ELEMS, /* + */
449 MOP_SUB_ELEMS, /* - */
450 MOP_MUL_ELEMS, /* &* */
451 MOP_DIV_ELEMS, /* / and &/ */
452 MOP_EXP_ELEMS, /* &** */
454 MOP_SEQ_BY, /* a:b:c */
456 /* Matrix arithmetic. */
458 MOP_EXP_MAT, /* ** */
475 MOP_PASTE_HORZ, /* a, b, c, ... */
476 MOP_PASTE_VERT, /* a; b; c; ... */
480 MOP_VEC_INDEX, /* x(y) */
481 MOP_VEC_ALL, /* x(:) */
482 MOP_MAT_INDEX, /* x(y,z) */
483 MOP_ROW_INDEX, /* x(y,:) */
484 MOP_COL_INDEX, /* x(:,z) */
491 MOP_EOF, /* EOF('file') */
499 struct matrix_expr **subs;
504 struct matrix_var *variable;
505 struct read_file *eof;
508 struct msg_location *location;
512 matrix_expr_destroy (struct matrix_expr *e)
519 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
550 for (size_t i = 0; i < e->n_subs; i++)
551 matrix_expr_destroy (e->subs[i]);
560 msg_location_destroy (e->location);
564 static struct matrix_expr *
565 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
568 struct matrix_expr *e = xmalloc (sizeof *e);
569 *e = (struct matrix_expr) {
571 .subs = xmemdup (subs, n_subs * sizeof *subs),
578 static struct matrix_expr *
579 matrix_expr_create_0 (enum matrix_op op)
581 struct matrix_expr *sub;
582 return matrix_expr_create_subs (op, &sub, 0);
585 static struct matrix_expr *
586 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
588 return matrix_expr_create_subs (op, &sub, 1);
591 static struct matrix_expr *
592 matrix_expr_create_2 (enum matrix_op op,
593 struct matrix_expr *sub0, struct matrix_expr *sub1)
595 struct matrix_expr *subs[] = { sub0, sub1 };
596 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
599 static struct matrix_expr *
600 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
601 struct matrix_expr *sub1, struct matrix_expr *sub2)
603 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
604 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
607 static struct matrix_expr *
608 matrix_expr_create_number (struct lexer *lexer, double number)
610 struct matrix_expr *e = xmalloc (sizeof *e);
611 *e = (struct matrix_expr) {
614 .location = lex_get_location (lexer, 0, 0),
620 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
622 static struct matrix_expr *
623 matrix_parse_curly_comma (struct matrix_state *s)
625 struct matrix_expr *lhs = matrix_parse_or_xor (s);
629 while (lex_match (s->lexer, T_COMMA))
631 struct matrix_expr *rhs = matrix_parse_or_xor (s);
634 matrix_expr_destroy (lhs);
637 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
642 static struct matrix_expr *
643 matrix_parse_curly_semi (struct matrix_state *s)
645 if (lex_token (s->lexer) == T_RCURLY)
647 struct matrix_expr *e = matrix_expr_create_0 (MOP_EMPTY);
648 e->location = lex_get_location (s->lexer, -1, 0);
652 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
656 while (lex_match (s->lexer, T_SEMICOLON))
658 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
661 matrix_expr_destroy (lhs);
664 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
669 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
670 for (size_t Y = 0; Y < (M)->size1; Y++) \
671 for (size_t X = 0; X < (M)->size2; X++) \
672 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
675 is_vector (const gsl_matrix *m)
677 return m->size1 <= 1 || m->size2 <= 1;
681 to_vector (gsl_matrix *m)
683 return (m->size1 == 1
684 ? gsl_matrix_row (m, 0).vector
685 : gsl_matrix_column (m, 0).vector);
690 matrix_eval_ABS (double d)
696 matrix_eval_ALL (gsl_matrix *m)
698 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
705 matrix_eval_ANY (gsl_matrix *m)
707 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
714 matrix_eval_ARSIN (double d)
720 matrix_eval_ARTAN (double d)
726 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
730 for (size_t i = 0; i < n; i++)
735 gsl_matrix *block = gsl_matrix_calloc (r, c);
737 for (size_t i = 0; i < n; i++)
739 for (size_t y = 0; y < m[i]->size1; y++)
740 for (size_t x = 0; x < m[i]->size2; x++)
741 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
749 matrix_eval_CHOL (gsl_matrix *m, const struct matrix_expr *e)
751 if (!gsl_linalg_cholesky_decomp1 (m))
753 for (size_t y = 0; y < m->size1; y++)
754 for (size_t x = y + 1; x < m->size2; x++)
755 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
757 for (size_t y = 0; y < m->size1; y++)
758 for (size_t x = 0; x < y; x++)
759 gsl_matrix_set (m, y, x, 0);
764 msg_at (SE, e->location,
765 _("Input to CHOL function is not positive-definite."));
771 matrix_eval_col_extremum (gsl_matrix *m, bool min)
776 return gsl_matrix_alloc (1, 0);
778 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
779 for (size_t x = 0; x < m->size2; x++)
781 double ext = gsl_matrix_get (m, 0, x);
782 for (size_t y = 1; y < m->size1; y++)
784 double value = gsl_matrix_get (m, y, x);
785 if (min ? value < ext : value > ext)
788 gsl_matrix_set (cext, 0, x, ext);
794 matrix_eval_CMAX (gsl_matrix *m)
796 return matrix_eval_col_extremum (m, false);
800 matrix_eval_CMIN (gsl_matrix *m)
802 return matrix_eval_col_extremum (m, true);
806 matrix_eval_COS (double d)
812 matrix_eval_col_sum (gsl_matrix *m, bool square)
817 return gsl_matrix_alloc (1, 0);
819 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
820 for (size_t x = 0; x < m->size2; x++)
823 for (size_t y = 0; y < m->size1; y++)
825 double d = gsl_matrix_get (m, y, x);
826 sum += square ? pow2 (d) : d;
828 gsl_matrix_set (result, 0, x, sum);
834 matrix_eval_CSSQ (gsl_matrix *m)
836 return matrix_eval_col_sum (m, true);
840 matrix_eval_CSUM (gsl_matrix *m)
842 return matrix_eval_col_sum (m, false);
846 compare_double_3way (const void *a_, const void *b_)
848 const double *a = a_;
849 const double *b = b_;
850 return *a < *b ? -1 : *a > *b;
854 matrix_eval_DESIGN (gsl_matrix *m)
856 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
857 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
858 gsl_matrix_transpose_memcpy (&m2, m);
860 for (size_t y = 0; y < m2.size1; y++)
861 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
863 size_t *n = xcalloc (m2.size1, sizeof *n);
865 for (size_t i = 0; i < m2.size1; i++)
867 double *row = tmp + m2.size2 * i;
868 for (size_t j = 0; j < m2.size2; )
871 for (k = j + 1; k < m2.size2; k++)
872 if (row[j] != row[k])
874 row[n[i]++] = row[j];
879 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
884 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
886 for (size_t i = 0; i < m->size2; i++)
891 const double *unique = tmp + m2.size2 * i;
892 for (size_t j = 0; j < n[i]; j++, x++)
894 double value = unique[j];
895 for (size_t y = 0; y < m->size1; y++)
896 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
907 matrix_eval_DET (gsl_matrix *m)
909 gsl_permutation *p = gsl_permutation_alloc (m->size1);
911 gsl_linalg_LU_decomp (m, p, &signum);
912 gsl_permutation_free (p);
913 return gsl_linalg_LU_det (m, signum);
917 matrix_eval_DIAG (gsl_matrix *m)
919 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
920 for (size_t i = 0; i < diag->size1; i++)
921 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
926 is_symmetric (const gsl_matrix *m)
928 if (m->size1 != m->size2)
931 for (size_t y = 0; y < m->size1; y++)
932 for (size_t x = 0; x < y; x++)
933 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
940 compare_double_desc (const void *a_, const void *b_)
942 const double *a = a_;
943 const double *b = b_;
944 return *a > *b ? -1 : *a < *b;
948 matrix_eval_EVAL (gsl_matrix *m)
950 if (!is_symmetric (m))
952 msg (SE, _("Argument of EVAL must be symmetric."));
956 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
957 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
958 gsl_vector v_eval = to_vector (eval);
959 gsl_eigen_symm (m, &v_eval, w);
960 gsl_eigen_symm_free (w);
962 assert (v_eval.stride == 1);
963 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
969 matrix_eval_EXP (double d)
974 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
977 Charl Linssen <charl@itfromb.it>
981 matrix_eval_GINV (gsl_matrix *A)
986 gsl_matrix *tmp_mat = NULL;
989 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
990 tmp_mat = gsl_matrix_alloc (m, n);
991 gsl_matrix_transpose_memcpy (tmp_mat, A);
999 gsl_matrix *V = gsl_matrix_alloc (m, m);
1000 gsl_vector *u = gsl_vector_alloc (m);
1002 gsl_vector *tmp_vec = gsl_vector_alloc (m);
1003 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
1004 gsl_vector_free (tmp_vec);
1007 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
1008 gsl_matrix_set_zero (Sigma_pinv);
1009 double cutoff = 1e-15 * gsl_vector_max (u);
1011 for (size_t i = 0; i < m; ++i)
1013 double x = gsl_vector_get (u, i);
1014 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
1017 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
1018 gsl_matrix *U = gsl_matrix_calloc (n, n);
1019 for (size_t i = 0; i < n; i++)
1020 for (size_t j = 0; j < m; j++)
1021 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
1023 /* two dot products to obtain pseudoinverse */
1024 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
1025 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
1030 A_pinv = gsl_matrix_alloc (n, m);
1031 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
1035 A_pinv = gsl_matrix_alloc (m, n);
1036 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
1039 gsl_matrix_free (tmp_mat);
1040 gsl_matrix_free (tmp_mat2);
1041 gsl_matrix_free (U);
1042 gsl_matrix_free (Sigma_pinv);
1043 gsl_vector_free (u);
1044 gsl_matrix_free (V);
1056 grade_compare_3way (const void *a_, const void *b_)
1058 const struct grade *a = a_;
1059 const struct grade *b = b_;
1061 return (a->value < b->value ? -1
1062 : a->value > b->value ? 1
1070 matrix_eval_GRADE (gsl_matrix *m)
1072 size_t n = m->size1 * m->size2;
1073 struct grade *grades = xmalloc (n * sizeof *grades);
1076 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1077 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
1078 qsort (grades, n, sizeof *grades, grade_compare_3way);
1080 for (size_t i = 0; i < n; i++)
1081 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
1089 dot (gsl_vector *a, gsl_vector *b)
1091 double result = 0.0;
1092 for (size_t i = 0; i < a->size; i++)
1093 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
1098 norm2 (gsl_vector *v)
1100 double result = 0.0;
1101 for (size_t i = 0; i < v->size; i++)
1102 result += pow2 (gsl_vector_get (v, i));
1107 norm (gsl_vector *v)
1109 return sqrt (norm2 (v));
1113 matrix_eval_GSCH (gsl_matrix *v)
1115 if (v->size2 < v->size1)
1117 msg (SE, _("GSCH requires its argument to have at least as many columns "
1118 "as rows, but it has dimensions %zu×%zu."),
1119 v->size1, v->size2);
1122 if (!v->size1 || !v->size2)
1125 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
1127 for (size_t vx = 0; vx < v->size2; vx++)
1129 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
1130 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
1132 gsl_vector_memcpy (&u_i, &v_i);
1133 for (size_t j = 0; j < ux; j++)
1135 gsl_vector u_j = gsl_matrix_column (u, j).vector;
1136 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
1137 for (size_t k = 0; k < u_i.size; k++)
1138 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
1139 - scale * gsl_vector_get (&u_j, k)));
1142 double len = norm (&u_i);
1145 gsl_vector_scale (&u_i, 1.0 / len);
1146 if (++ux >= v->size1)
1153 msg (SE, _("%zu×%zu argument to GSCH contains only "
1154 "%zu linearly independent columns."),
1155 v->size1, v->size2, ux);
1156 gsl_matrix_free (u);
1160 u->size2 = v->size1;
1165 matrix_eval_IDENT (double s1, double s2)
1167 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
1168 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1174 invert_matrix (gsl_matrix *x)
1176 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1178 gsl_linalg_LU_decomp (x, p, &signum);
1179 gsl_linalg_LU_invx (x, p);
1180 gsl_permutation_free (p);
1184 matrix_eval_INV (gsl_matrix *m)
1191 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1193 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1194 a->size2 * b->size2);
1196 for (size_t ar = 0; ar < a->size1; ar++)
1197 for (size_t br = 0; br < b->size1; br++, y++)
1200 for (size_t ac = 0; ac < a->size2; ac++)
1201 for (size_t bc = 0; bc < b->size2; bc++, x++)
1203 double av = gsl_matrix_get (a, ar, ac);
1204 double bv = gsl_matrix_get (b, br, bc);
1205 gsl_matrix_set (k, y, x, av * bv);
1212 matrix_eval_LG10 (double d)
1218 matrix_eval_LN (double d)
1224 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1226 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1229 for (size_t i = 1; i <= n * n; i++)
1231 gsl_matrix_set (m, y, x, i);
1233 size_t y1 = !y ? n - 1 : y - 1;
1234 size_t x1 = x + 1 >= n ? 0 : x + 1;
1235 if (gsl_matrix_get (m, y1, x1) == 0)
1241 y = y + 1 >= n ? 0 : y + 1;
1246 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1248 double a = gsl_matrix_get (m, y1, x1);
1249 double b = gsl_matrix_get (m, y2, x2);
1250 gsl_matrix_set (m, y1, x1, b);
1251 gsl_matrix_set (m, y2, x2, a);
1255 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1259 /* A. Umar, "On the Construction of Even Order Magic Squares",
1260 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1262 for (size_t i = 1; i <= n * n / 2; i++)
1264 gsl_matrix_set (m, y, x, i);
1274 for (size_t i = n * n; i > n * n / 2; i--)
1276 gsl_matrix_set (m, y, x, i);
1284 for (size_t y = 0; y < n; y++)
1285 for (size_t x = 0; x < n / 2; x++)
1287 unsigned int d = gsl_matrix_get (m, y, x);
1288 if (d % 2 != (y < n / 2))
1289 magic_exchange (m, y, x, y, n - x - 1);
1294 size_t x1 = n / 2 - 1;
1296 magic_exchange (m, y1, x1, y2, x1);
1297 magic_exchange (m, y1, x2, y2, x2);
1301 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1303 /* A. Umar, "On the Construction of Even Order Magic Squares",
1304 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1308 for (size_t i = 1; ; i++)
1310 gsl_matrix_set (m, y, x, i);
1311 if (++y == n / 2 - 1)
1323 for (size_t i = n * n; ; i--)
1325 gsl_matrix_set (m, y, x, i);
1326 if (++y == n / 2 - 1)
1335 for (size_t y = 0; y < n; y++)
1336 if (y != n / 2 - 1 && y != n / 2)
1337 for (size_t x = 0; x < n / 2; x++)
1339 unsigned int d = gsl_matrix_get (m, y, x);
1340 if (d % 2 != (y < n / 2))
1341 magic_exchange (m, y, x, y, n - x - 1);
1344 size_t a0 = (n * n - 2 * n) / 2 + 1;
1345 for (size_t i = 0; i < n / 2; i++)
1348 gsl_matrix_set (m, n / 2, i, a);
1349 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1351 for (size_t i = 0; i < n / 2; i++)
1353 size_t a = a0 + i + n / 2;
1354 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1355 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1357 for (size_t x = 1; x < n / 2; x += 2)
1358 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1359 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1360 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1361 size_t x1 = n / 2 - 2;
1362 size_t x2 = n / 2 + 1;
1363 size_t y1 = n / 2 - 2;
1364 size_t y2 = n / 2 + 1;
1365 magic_exchange (m, y1, x1, y2, x1);
1366 magic_exchange (m, y1, x2, y2, x2);
1370 matrix_eval_MAGIC (double n_)
1374 gsl_matrix *m = gsl_matrix_calloc (n, n);
1376 matrix_eval_MAGIC_odd (m, n);
1378 matrix_eval_MAGIC_singly_even (m, n);
1380 matrix_eval_MAGIC_doubly_even (m, n);
1385 matrix_eval_MAKE (double r, double c, double value)
1387 gsl_matrix *m = gsl_matrix_alloc (r, c);
1388 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1394 matrix_eval_MDIAG (gsl_vector *v)
1396 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1397 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1398 gsl_vector_memcpy (&diagonal, v);
1403 matrix_eval_MMAX (gsl_matrix *m)
1405 return gsl_matrix_max (m);
1409 matrix_eval_MMIN (gsl_matrix *m)
1411 return gsl_matrix_min (m);
1415 matrix_eval_MOD (gsl_matrix *m, double divisor)
1417 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1418 *d = fmod (*d, divisor);
1423 matrix_eval_MSSQ (gsl_matrix *m)
1426 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1432 matrix_eval_MSUM (gsl_matrix *m)
1435 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1441 matrix_eval_NCOL (gsl_matrix *m)
1447 matrix_eval_NROW (gsl_matrix *m)
1453 matrix_eval_RANK (gsl_matrix *m)
1455 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1456 gsl_linalg_QR_decomp (m, tau);
1457 gsl_vector_free (tau);
1459 return gsl_linalg_QRPT_rank (m, -1);
1463 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1465 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1467 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1472 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1474 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1475 "product of matrix dimensions."));
1479 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1482 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1484 gsl_matrix_set (dst, y1, x1, *d);
1495 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1500 return gsl_matrix_alloc (0, 1);
1502 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1503 for (size_t y = 0; y < m->size1; y++)
1505 double ext = gsl_matrix_get (m, y, 0);
1506 for (size_t x = 1; x < m->size2; x++)
1508 double value = gsl_matrix_get (m, y, x);
1509 if (min ? value < ext : value > ext)
1512 gsl_matrix_set (rext, y, 0, ext);
1518 matrix_eval_RMAX (gsl_matrix *m)
1520 return matrix_eval_row_extremum (m, false);
1524 matrix_eval_RMIN (gsl_matrix *m)
1526 return matrix_eval_row_extremum (m, true);
1530 matrix_eval_RND (double d)
1542 rank_compare_3way (const void *a_, const void *b_)
1544 const struct rank *a = a_;
1545 const struct rank *b = b_;
1547 return a->value < b->value ? -1 : a->value > b->value;
1551 matrix_eval_RNKORDER (gsl_matrix *m)
1553 size_t n = m->size1 * m->size2;
1554 struct rank *ranks = xmalloc (n * sizeof *ranks);
1556 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1557 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1558 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1560 for (size_t i = 0; i < n; )
1563 for (j = i + 1; j < n; j++)
1564 if (ranks[i].value != ranks[j].value)
1567 double rank = (i + j + 1.0) / 2.0;
1568 for (size_t k = i; k < j; k++)
1569 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1580 matrix_eval_row_sum (gsl_matrix *m, bool square)
1585 return gsl_matrix_alloc (0, 1);
1587 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1588 for (size_t y = 0; y < m->size1; y++)
1591 for (size_t x = 0; x < m->size2; x++)
1593 double d = gsl_matrix_get (m, y, x);
1594 sum += square ? pow2 (d) : d;
1596 gsl_matrix_set (result, y, 0, sum);
1602 matrix_eval_RSSQ (gsl_matrix *m)
1604 return matrix_eval_row_sum (m, true);
1608 matrix_eval_RSUM (gsl_matrix *m)
1610 return matrix_eval_row_sum (m, false);
1614 matrix_eval_SIN (double d)
1620 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1622 if (m1->size1 != m2->size1)
1624 msg (SE, _("SOLVE requires its arguments to have the same number of "
1625 "rows, but the first argument has dimensions %zu×%zu and "
1626 "the second %zu×%zu."),
1627 m1->size1, m1->size2,
1628 m2->size1, m2->size2);
1632 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1633 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1635 gsl_linalg_LU_decomp (m1, p, &signum);
1636 for (size_t i = 0; i < m2->size2; i++)
1638 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1639 gsl_vector xi = gsl_matrix_column (x, i).vector;
1640 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1642 gsl_permutation_free (p);
1647 matrix_eval_SQRT (double d)
1653 matrix_eval_SSCP (gsl_matrix *m)
1655 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1656 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1661 matrix_eval_SVAL (gsl_matrix *m)
1663 gsl_matrix *tmp_mat = NULL;
1664 if (m->size2 > m->size1)
1666 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1667 gsl_matrix_transpose_memcpy (tmp_mat, m);
1672 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1673 gsl_vector *S = gsl_vector_alloc (m->size2);
1674 gsl_vector *work = gsl_vector_alloc (m->size2);
1675 gsl_linalg_SV_decomp (m, V, S, work);
1677 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1678 for (size_t i = 0; i < m->size2; i++)
1679 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1681 gsl_matrix_free (V);
1682 gsl_vector_free (S);
1683 gsl_vector_free (work);
1684 gsl_matrix_free (tmp_mat);
1690 matrix_eval_SWEEP (gsl_matrix *m, double d)
1692 if (d < 1 || d > SIZE_MAX)
1694 msg (SE, _("Scalar argument to SWEEP must be integer."));
1698 if (k >= MIN (m->size1, m->size2))
1700 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1701 "equal to the smaller of the matrix argument's rows and "
1706 double m_kk = gsl_matrix_get (m, k, k);
1707 if (fabs (m_kk) > 1e-19)
1709 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1710 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1712 double m_ij = gsl_matrix_get (m, i, j);
1713 double m_ik = gsl_matrix_get (m, i, k);
1714 double m_kj = gsl_matrix_get (m, k, j);
1715 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1724 for (size_t i = 0; i < m->size1; i++)
1726 gsl_matrix_set (m, i, k, 0);
1727 gsl_matrix_set (m, k, i, 0);
1734 matrix_eval_TRACE (gsl_matrix *m)
1737 size_t n = MIN (m->size1, m->size2);
1738 for (size_t i = 0; i < n; i++)
1739 sum += gsl_matrix_get (m, i, i);
1744 matrix_eval_T (gsl_matrix *m)
1746 return matrix_eval_TRANSPOS (m);
1750 matrix_eval_TRANSPOS (gsl_matrix *m)
1752 if (m->size1 == m->size2)
1754 gsl_matrix_transpose (m);
1759 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1760 gsl_matrix_transpose_memcpy (t, m);
1766 matrix_eval_TRUNC (double d)
1772 matrix_eval_UNIFORM (double r_, double c_)
1776 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1778 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1782 gsl_matrix *m = gsl_matrix_alloc (r, c);
1783 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1784 *d = gsl_ran_flat (get_rng (), 0, 1);
1789 matrix_eval_PDF_BETA (double x, double a, double b)
1791 return gsl_ran_beta_pdf (x, a, b);
1795 matrix_eval_CDF_BETA (double x, double a, double b)
1797 return gsl_cdf_beta_P (x, a, b);
1801 matrix_eval_IDF_BETA (double P, double a, double b)
1803 return gsl_cdf_beta_Pinv (P, a, b);
1807 matrix_eval_RV_BETA (double a, double b)
1809 return gsl_ran_beta (get_rng (), a, b);
1813 matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
1815 return ncdf_beta (x, a, b, lambda);
1819 matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
1821 return npdf_beta (x, a, b, lambda);
1825 matrix_eval_CDF_BVNOR (double x0, double x1, double r)
1827 return cdf_bvnor (x0, x1, r);
1831 matrix_eval_PDF_BVNOR (double x0, double x1, double r)
1833 return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
1837 matrix_eval_CDF_CAUCHY (double x, double a, double b)
1839 return gsl_cdf_cauchy_P ((x - a) / b, 1);
1843 matrix_eval_IDF_CAUCHY (double P, double a, double b)
1845 return a + b * gsl_cdf_cauchy_Pinv (P, 1);
1849 matrix_eval_PDF_CAUCHY (double x, double a, double b)
1851 return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
1855 matrix_eval_RV_CAUCHY (double a, double b)
1857 return a + b * gsl_ran_cauchy (get_rng (), 1);
1861 matrix_eval_CDF_CHISQ (double x, double df)
1863 return gsl_cdf_chisq_P (x, df);
1867 matrix_eval_CHICDF (double x, double df)
1869 return matrix_eval_CDF_CHISQ (x, df);
1873 matrix_eval_IDF_CHISQ (double P, double df)
1875 return gsl_cdf_chisq_Pinv (P, df);
1879 matrix_eval_PDF_CHISQ (double x, double df)
1881 return gsl_ran_chisq_pdf (x, df);
1885 matrix_eval_RV_CHISQ (double df)
1887 return gsl_ran_chisq (get_rng (), df);
1891 matrix_eval_SIG_CHISQ (double x, double df)
1893 return gsl_cdf_chisq_Q (x, df);
1897 matrix_eval_CDF_EXP (double x, double a)
1899 return gsl_cdf_exponential_P (x, 1. / a);
1903 matrix_eval_IDF_EXP (double P, double a)
1905 return gsl_cdf_exponential_Pinv (P, 1. / a);
1909 matrix_eval_PDF_EXP (double x, double a)
1911 return gsl_ran_exponential_pdf (x, 1. / a);
1915 matrix_eval_RV_EXP (double a)
1917 return gsl_ran_exponential (get_rng (), 1. / a);
1921 matrix_eval_PDF_XPOWER (double x, double a, double b)
1923 return gsl_ran_exppow_pdf (x, a, b);
1927 matrix_eval_RV_XPOWER (double a, double b)
1929 return gsl_ran_exppow (get_rng (), a, b);
1933 matrix_eval_CDF_F (double x, double df1, double df2)
1935 return gsl_cdf_fdist_P (x, df1, df2);
1939 matrix_eval_FCDF (double x, double df1, double df2)
1941 return matrix_eval_CDF_F (x, df1, df2);
1945 matrix_eval_IDF_F (double P, double df1, double df2)
1947 return idf_fdist (P, df1, df2);
1951 matrix_eval_RV_F (double df1, double df2)
1953 return gsl_ran_fdist (get_rng (), df1, df2);
1957 matrix_eval_PDF_F (double x, double df1, double df2)
1959 return gsl_ran_fdist_pdf (x, df1, df2);
1963 matrix_eval_SIG_F (double x, double df1, double df2)
1965 return gsl_cdf_fdist_Q (x, df1, df2);
1969 matrix_eval_CDF_GAMMA (double x, double a, double b)
1971 return gsl_cdf_gamma_P (x, a, 1. / b);
1975 matrix_eval_IDF_GAMMA (double P, double a, double b)
1977 return gsl_cdf_gamma_Pinv (P, a, 1. / b);
1981 matrix_eval_PDF_GAMMA (double x, double a, double b)
1983 return gsl_ran_gamma_pdf (x, a, 1. / b);
1987 matrix_eval_RV_GAMMA (double a, double b)
1989 return gsl_ran_gamma (get_rng (), a, 1. / b);
1993 matrix_eval_PDF_LANDAU (double x)
1995 return gsl_ran_landau_pdf (x);
1999 matrix_eval_RV_LANDAU (void)
2001 return gsl_ran_landau (get_rng ());
2005 matrix_eval_CDF_LAPLACE (double x, double a, double b)
2007 return gsl_cdf_laplace_P ((x - a) / b, 1);
2011 matrix_eval_IDF_LAPLACE (double P, double a, double b)
2013 return a + b * gsl_cdf_laplace_Pinv (P, 1);
2017 matrix_eval_PDF_LAPLACE (double x, double a, double b)
2019 return gsl_ran_laplace_pdf ((x - a) / b, 1);
2023 matrix_eval_RV_LAPLACE (double a, double b)
2025 return a + b * gsl_ran_laplace (get_rng (), 1);
2029 matrix_eval_RV_LEVY (double c, double alpha)
2031 return gsl_ran_levy (get_rng (), c, alpha);
2035 matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
2037 return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
2041 matrix_eval_CDF_LOGISTIC (double x, double a, double b)
2043 return gsl_cdf_logistic_P ((x - a) / b, 1);
2047 matrix_eval_IDF_LOGISTIC (double P, double a, double b)
2049 return a + b * gsl_cdf_logistic_Pinv (P, 1);
2053 matrix_eval_PDF_LOGISTIC (double x, double a, double b)
2055 return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
2059 matrix_eval_RV_LOGISTIC (double a, double b)
2061 return a + b * gsl_ran_logistic (get_rng (), 1);
2065 matrix_eval_CDF_LNORMAL (double x, double m, double s)
2067 return gsl_cdf_lognormal_P (x, log (m), s);
2071 matrix_eval_IDF_LNORMAL (double P, double m, double s)
2073 return gsl_cdf_lognormal_Pinv (P, log (m), s);;
2077 matrix_eval_PDF_LNORMAL (double x, double m, double s)
2079 return gsl_ran_lognormal_pdf (x, log (m), s);
2083 matrix_eval_RV_LNORMAL (double m, double s)
2085 return gsl_ran_lognormal (get_rng (), log (m), s);
2089 matrix_eval_CDF_NORMAL (double x, double u, double s)
2091 return gsl_cdf_gaussian_P (x - u, s);
2095 matrix_eval_IDF_NORMAL (double P, double u, double s)
2097 return u + gsl_cdf_gaussian_Pinv (P, s);
2101 matrix_eval_PDF_NORMAL (double x, double u, double s)
2103 return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
2107 matrix_eval_RV_NORMAL (double u, double s)
2109 return u + gsl_ran_gaussian (get_rng (), s);
2113 matrix_eval_CDFNORM (double x)
2115 return gsl_cdf_ugaussian_P (x);
2119 matrix_eval_PROBIT (double P)
2121 return gsl_cdf_ugaussian_Pinv (P);
2125 matrix_eval_NORMAL (double s)
2127 return gsl_ran_gaussian (get_rng (), s);
2131 matrix_eval_PDF_NTAIL (double x, double a, double sigma)
2133 return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
2137 matrix_eval_RV_NTAIL (double a, double sigma)
2139 return gsl_ran_gaussian_tail (get_rng (), a, sigma);
2143 matrix_eval_CDF_PARETO (double x, double a, double b)
2145 return gsl_cdf_pareto_P (x, b, a);
2149 matrix_eval_IDF_PARETO (double P, double a, double b)
2151 return gsl_cdf_pareto_Pinv (P, b, a);
2155 matrix_eval_PDF_PARETO (double x, double a, double b)
2157 return gsl_ran_pareto_pdf (x, b, a);
2161 matrix_eval_RV_PARETO (double a, double b)
2163 return gsl_ran_pareto (get_rng (), b, a);
2167 matrix_eval_CDF_RAYLEIGH (double x, double sigma)
2169 return gsl_cdf_rayleigh_P (x, sigma);
2173 matrix_eval_IDF_RAYLEIGH (double P, double sigma)
2175 return gsl_cdf_rayleigh_Pinv (P, sigma);
2179 matrix_eval_PDF_RAYLEIGH (double x, double sigma)
2181 return gsl_ran_rayleigh_pdf (x, sigma);
2185 matrix_eval_RV_RAYLEIGH (double sigma)
2187 return gsl_ran_rayleigh (get_rng (), sigma);
2191 matrix_eval_PDF_RTAIL (double x, double a, double sigma)
2193 return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
2197 matrix_eval_RV_RTAIL (double a, double sigma)
2199 return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
2203 matrix_eval_CDF_T (double x, double df)
2205 return gsl_cdf_tdist_P (x, df);
2209 matrix_eval_TCDF (double x, double df)
2211 return matrix_eval_CDF_T (x, df);
2215 matrix_eval_IDF_T (double P, double df)
2217 return gsl_cdf_tdist_Pinv (P, df);
2221 matrix_eval_PDF_T (double x, double df)
2223 return gsl_ran_tdist_pdf (x, df);
2227 matrix_eval_RV_T (double df)
2229 return gsl_ran_tdist (get_rng (), df);
2233 matrix_eval_CDF_T1G (double x, double a, double b)
2235 return gsl_cdf_gumbel1_P (x, a, b);
2239 matrix_eval_IDF_T1G (double P, double a, double b)
2241 return gsl_cdf_gumbel1_Pinv (P, a, b);
2245 matrix_eval_PDF_T1G (double x, double a, double b)
2247 return gsl_ran_gumbel1_pdf (x, a, b);
2251 matrix_eval_RV_T1G (double a, double b)
2253 return gsl_ran_gumbel1 (get_rng (), a, b);
2257 matrix_eval_CDF_T2G (double x, double a, double b)
2259 return gsl_cdf_gumbel1_P (x, a, b);
2263 matrix_eval_IDF_T2G (double P, double a, double b)
2265 return gsl_cdf_gumbel1_Pinv (P, a, b);
2269 matrix_eval_PDF_T2G (double x, double a, double b)
2271 return gsl_ran_gumbel1_pdf (x, a, b);
2275 matrix_eval_RV_T2G (double a, double b)
2277 return gsl_ran_gumbel1 (get_rng (), a, b);
2281 matrix_eval_CDF_UNIFORM (double x, double a, double b)
2283 return gsl_cdf_flat_P (x, a, b);
2287 matrix_eval_IDF_UNIFORM (double P, double a, double b)
2289 return gsl_cdf_flat_Pinv (P, a, b);
2293 matrix_eval_PDF_UNIFORM (double x, double a, double b)
2295 return gsl_ran_flat_pdf (x, a, b);
2299 matrix_eval_RV_UNIFORM (double a, double b)
2301 return gsl_ran_flat (get_rng (), a, b);
2305 matrix_eval_CDF_WEIBULL (double x, double a, double b)
2307 return gsl_cdf_weibull_P (x, a, b);
2311 matrix_eval_IDF_WEIBULL (double P, double a, double b)
2313 return gsl_cdf_weibull_Pinv (P, a, b);
2317 matrix_eval_PDF_WEIBULL (double x, double a, double b)
2319 return gsl_ran_weibull_pdf (x, a, b);
2323 matrix_eval_RV_WEIBULL (double a, double b)
2325 return gsl_ran_weibull (get_rng (), a, b);
2329 matrix_eval_CDF_BERNOULLI (double k, double p)
2331 return k ? 1 : 1 - p;
2335 matrix_eval_PDF_BERNOULLI (double k, double p)
2337 return gsl_ran_bernoulli_pdf (k, p);
2341 matrix_eval_RV_BERNOULLI (double p)
2343 return gsl_ran_bernoulli (get_rng (), p);
2347 matrix_eval_CDF_BINOM (double k, double n, double p)
2349 return gsl_cdf_binomial_P (k, p, n);
2353 matrix_eval_PDF_BINOM (double k, double n, double p)
2355 return gsl_ran_binomial_pdf (k, p, n);
2359 matrix_eval_RV_BINOM (double n, double p)
2361 return gsl_ran_binomial (get_rng (), p, n);
2365 matrix_eval_CDF_GEOM (double k, double p)
2367 return gsl_cdf_geometric_P (k, p);
2371 matrix_eval_PDF_GEOM (double k, double p)
2373 return gsl_ran_geometric_pdf (k, p);
2377 matrix_eval_RV_GEOM (double p)
2379 return gsl_ran_geometric (get_rng (), p);
2383 matrix_eval_CDF_HYPER (double k, double a, double b, double c)
2385 return gsl_cdf_hypergeometric_P (k, c, a - c, b);
2389 matrix_eval_PDF_HYPER (double k, double a, double b, double c)
2391 return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
2395 matrix_eval_RV_HYPER (double a, double b, double c)
2397 return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
2401 matrix_eval_PDF_LOG (double k, double p)
2403 return gsl_ran_logarithmic_pdf (k, p);
2407 matrix_eval_RV_LOG (double p)
2409 return gsl_ran_logarithmic (get_rng (), p);
2413 matrix_eval_CDF_NEGBIN (double k, double n, double p)
2415 return gsl_cdf_negative_binomial_P (k, p, n);
2419 matrix_eval_PDF_NEGBIN (double k, double n, double p)
2421 return gsl_ran_negative_binomial_pdf (k, p, n);
2425 matrix_eval_RV_NEGBIN (double n, double p)
2427 return gsl_ran_negative_binomial (get_rng (), p, n);
2431 matrix_eval_CDF_POISSON (double k, double mu)
2433 return gsl_cdf_poisson_P (k, mu);
2437 matrix_eval_PDF_POISSON (double k, double mu)
2439 return gsl_ran_poisson_pdf (k, mu);
2443 matrix_eval_RV_POISSON (double mu)
2445 return gsl_ran_poisson (get_rng (), mu);
2448 struct matrix_function
2452 size_t min_args, max_args;
2455 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
2458 word_matches (const char **test, const char **name)
2460 size_t test_len = strcspn (*test, ".");
2461 size_t name_len = strcspn (*name, ".");
2462 if (test_len == name_len)
2464 if (buf_compare_case (*test, *name, test_len))
2467 else if (test_len < 3 || test_len > name_len)
2471 if (buf_compare_case (*test, *name, test_len))
2477 if (**test != **name)
2488 /* Returns 0 if TOKEN and FUNC do not match,
2489 1 if TOKEN is an acceptable abbreviation for FUNC,
2490 2 if TOKEN equals FUNC. */
2492 compare_function_names (const char *token_, const char *func_)
2494 const char *token = token_;
2495 const char *func = func_;
2496 while (*token || *func)
2497 if (!word_matches (&token, &func))
2499 return !c_strcasecmp (token_, func_) ? 2 : 1;
2502 static const struct matrix_function *
2503 matrix_parse_function_name (const char *token)
2505 static const struct matrix_function functions[] =
2507 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
2508 { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
2512 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
2514 for (size_t i = 0; i < N_FUNCTIONS; i++)
2516 if (compare_function_names (token, functions[i].name) > 0)
2517 return &functions[i];
2522 static struct read_file *
2523 read_file_create (struct matrix_state *s, struct file_handle *fh)
2525 for (size_t i = 0; i < s->n_read_files; i++)
2527 struct read_file *rf = s->read_files[i];
2535 struct read_file *rf = xmalloc (sizeof *rf);
2536 *rf = (struct read_file) { .file = fh };
2538 s->read_files = xrealloc (s->read_files,
2539 (s->n_read_files + 1) * sizeof *s->read_files);
2540 s->read_files[s->n_read_files++] = rf;
2544 static struct dfm_reader *
2545 read_file_open (struct read_file *rf)
2548 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
2553 read_file_destroy (struct read_file *rf)
2557 fh_unref (rf->file);
2558 dfm_close_reader (rf->reader);
2559 free (rf->encoding);
2564 static struct write_file *
2565 write_file_create (struct matrix_state *s, struct file_handle *fh)
2567 for (size_t i = 0; i < s->n_write_files; i++)
2569 struct write_file *wf = s->write_files[i];
2577 struct write_file *wf = xmalloc (sizeof *wf);
2578 *wf = (struct write_file) { .file = fh };
2580 s->write_files = xrealloc (s->write_files,
2581 (s->n_write_files + 1) * sizeof *s->write_files);
2582 s->write_files[s->n_write_files++] = wf;
2586 static struct dfm_writer *
2587 write_file_open (struct write_file *wf)
2590 wf->writer = dfm_open_writer (wf->file, wf->encoding);
2595 write_file_destroy (struct write_file *wf)
2601 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
2602 wf->held->s.ss.length);
2603 u8_line_destroy (wf->held);
2607 fh_unref (wf->file);
2608 dfm_close_writer (wf->writer);
2609 free (wf->encoding);
2615 matrix_parse_function (struct matrix_state *s, const char *token,
2616 struct matrix_expr **exprp)
2619 if (lex_next_token (s->lexer, 1) != T_LPAREN)
2622 if (lex_match_id (s->lexer, "EOF"))
2625 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
2629 if (!lex_force_match (s->lexer, T_RPAREN))
2635 struct read_file *rf = read_file_create (s, fh);
2637 struct matrix_expr *e = xmalloc (sizeof *e);
2638 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
2643 const struct matrix_function *f = matrix_parse_function_name (token);
2647 struct matrix_expr *e = xmalloc (sizeof *e);
2648 *e = (struct matrix_expr) {
2650 .location = lex_get_location (s->lexer, 0, 0)
2653 lex_get_n (s->lexer, 2);
2654 if (lex_token (s->lexer) != T_RPAREN)
2656 size_t allocated_subs = 0;
2659 struct msg_location *arg_location = lex_get_location (s->lexer, 0, 0);
2660 struct matrix_expr *sub = matrix_parse_expr (s);
2663 msg_location_destroy (arg_location);
2668 //lex_extend_location (s->lexer, 0, arg_location);
2669 sub->location = arg_location;
2672 msg_location_destroy (arg_location);
2674 if (e->n_subs >= allocated_subs)
2675 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
2676 e->subs[e->n_subs++] = sub;
2678 while (lex_match (s->lexer, T_COMMA));
2680 if (!lex_force_match (s->lexer, T_RPAREN))
2683 if (e->n_subs < f->min_args || e->n_subs > f->max_args)
2685 if (f->min_args == f->max_args)
2686 msg (SE, ngettext ("Matrix function %s requires %zu argument.",
2687 "Matrix function %s requires %zu arguments.",
2689 f->name, f->min_args);
2690 else if (f->min_args == 1 && f->max_args == 2)
2691 msg (SE, ngettext ("Matrix function %s requires 1 or 2 arguments, "
2692 "but %zu was provided.",
2693 "Matrix function %s requires 1 or 2 arguments, "
2694 "but %zu were provided.",
2696 f->name, e->n_subs);
2697 else if (f->min_args == 1 && f->max_args == INT_MAX)
2698 msg (SE, _("Matrix function %s requires at least one argument."),
2710 matrix_expr_destroy (e);
2714 static struct matrix_expr *
2715 matrix_parse_primary (struct matrix_state *s)
2717 if (lex_is_number (s->lexer))
2719 double number = lex_number (s->lexer);
2720 return matrix_expr_create_number (s->lexer, number);
2722 else if (lex_is_string (s->lexer))
2724 char string[sizeof (double)];
2725 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
2728 memcpy (&number, string, sizeof number);
2729 return matrix_expr_create_number (s->lexer, number);
2731 else if (lex_match (s->lexer, T_LPAREN))
2733 struct matrix_expr *e = matrix_parse_or_xor (s);
2734 if (!e || !lex_force_match (s->lexer, T_RPAREN))
2736 matrix_expr_destroy (e);
2741 else if (lex_match (s->lexer, T_LCURLY))
2743 struct matrix_expr *e = matrix_parse_curly_semi (s);
2744 if (!e || !lex_force_match (s->lexer, T_RCURLY))
2746 matrix_expr_destroy (e);
2751 else if (lex_token (s->lexer) == T_ID)
2753 struct matrix_expr *retval;
2754 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
2757 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
2760 lex_error (s->lexer, _("Unknown variable %s."),
2761 lex_tokcstr (s->lexer));
2766 struct matrix_expr *e = xmalloc (sizeof *e);
2767 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
2770 else if (lex_token (s->lexer) == T_ALL)
2772 struct matrix_expr *retval;
2773 if (matrix_parse_function (s, "ALL", &retval))
2777 lex_error (s->lexer, NULL);
2781 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
2784 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
2786 if (lex_match (s->lexer, T_COLON))
2793 *indexp = matrix_parse_expr (s);
2794 return *indexp != NULL;
2798 static struct matrix_expr *
2799 matrix_parse_postfix (struct matrix_state *s)
2801 struct matrix_expr *lhs = matrix_parse_primary (s);
2802 if (!lhs || !lex_match (s->lexer, T_LPAREN))
2805 struct matrix_expr *i0;
2806 if (!matrix_parse_index_expr (s, &i0))
2808 matrix_expr_destroy (lhs);
2811 if (lex_match (s->lexer, T_RPAREN))
2813 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
2814 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
2815 else if (lex_match (s->lexer, T_COMMA))
2817 struct matrix_expr *i1;
2818 if (!matrix_parse_index_expr (s, &i1)
2819 || !lex_force_match (s->lexer, T_RPAREN))
2821 matrix_expr_destroy (lhs);
2822 matrix_expr_destroy (i0);
2823 matrix_expr_destroy (i1);
2826 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
2827 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
2828 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
2833 lex_error_expecting (s->lexer, "`)'", "`,'");
2838 static struct matrix_expr *
2839 matrix_parse_unary (struct matrix_state *s)
2841 if (lex_match (s->lexer, T_DASH))
2843 struct matrix_expr *lhs = matrix_parse_unary (s);
2844 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
2846 else if (lex_match (s->lexer, T_PLUS))
2847 return matrix_parse_unary (s);
2849 return matrix_parse_postfix (s);
2852 static struct matrix_expr *
2853 matrix_parse_seq (struct matrix_state *s)
2855 struct matrix_expr *start = matrix_parse_unary (s);
2856 if (!start || !lex_match (s->lexer, T_COLON))
2859 struct matrix_expr *end = matrix_parse_unary (s);
2862 matrix_expr_destroy (start);
2866 if (lex_match (s->lexer, T_COLON))
2868 struct matrix_expr *increment = matrix_parse_unary (s);
2871 matrix_expr_destroy (start);
2872 matrix_expr_destroy (end);
2875 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2878 return matrix_expr_create_2 (MOP_SEQ, start, end);
2881 static struct matrix_expr *
2882 matrix_parse_exp (struct matrix_state *s)
2884 struct matrix_expr *lhs = matrix_parse_seq (s);
2891 if (lex_match (s->lexer, T_EXP))
2893 else if (lex_match_phrase (s->lexer, "&**"))
2898 struct matrix_expr *rhs = matrix_parse_seq (s);
2901 matrix_expr_destroy (lhs);
2904 lhs = matrix_expr_create_2 (op, lhs, rhs);
2908 static struct matrix_expr *
2909 matrix_parse_mul_div (struct matrix_state *s)
2911 struct matrix_expr *lhs = matrix_parse_exp (s);
2918 if (lex_match (s->lexer, T_ASTERISK))
2920 else if (lex_match (s->lexer, T_SLASH))
2922 else if (lex_match_phrase (s->lexer, "&*"))
2924 else if (lex_match_phrase (s->lexer, "&/"))
2929 struct matrix_expr *rhs = matrix_parse_exp (s);
2932 matrix_expr_destroy (lhs);
2935 lhs = matrix_expr_create_2 (op, lhs, rhs);
2939 static struct matrix_expr *
2940 matrix_parse_add_sub (struct matrix_state *s)
2942 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2949 if (lex_match (s->lexer, T_PLUS))
2951 else if (lex_match (s->lexer, T_DASH))
2953 else if (lex_token (s->lexer) == T_NEG_NUM)
2958 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2961 matrix_expr_destroy (lhs);
2964 lhs = matrix_expr_create_2 (op, lhs, rhs);
2968 static struct matrix_expr *
2969 matrix_parse_relational (struct matrix_state *s)
2971 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2978 if (lex_match (s->lexer, T_GT))
2980 else if (lex_match (s->lexer, T_GE))
2982 else if (lex_match (s->lexer, T_LT))
2984 else if (lex_match (s->lexer, T_LE))
2986 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2988 else if (lex_match (s->lexer, T_NE))
2993 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2996 matrix_expr_destroy (lhs);
2999 lhs = matrix_expr_create_2 (op, lhs, rhs);
3003 static struct matrix_expr *
3004 matrix_parse_not (struct matrix_state *s)
3006 if (lex_match (s->lexer, T_NOT))
3008 struct matrix_expr *lhs = matrix_parse_not (s);
3009 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
3012 return matrix_parse_relational (s);
3015 static struct matrix_expr *
3016 matrix_parse_and (struct matrix_state *s)
3018 struct matrix_expr *lhs = matrix_parse_not (s);
3022 while (lex_match (s->lexer, T_AND))
3024 struct matrix_expr *rhs = matrix_parse_not (s);
3027 matrix_expr_destroy (lhs);
3030 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
3035 static struct matrix_expr *
3036 matrix_parse_or_xor (struct matrix_state *s)
3038 struct matrix_expr *lhs = matrix_parse_and (s);
3045 if (lex_match (s->lexer, T_OR))
3047 else if (lex_match_id (s->lexer, "XOR"))
3052 struct matrix_expr *rhs = matrix_parse_and (s);
3055 matrix_expr_destroy (lhs);
3058 lhs = matrix_expr_create_2 (op, lhs, rhs);
3062 static struct matrix_expr *
3063 matrix_parse_expr (struct matrix_state *s)
3065 return matrix_parse_or_xor (s);
3068 /* Expression evaluation. */
3071 matrix_op_eval (enum matrix_op op, double a, double b)
3075 case MOP_ADD_ELEMS: return a + b;
3076 case MOP_SUB_ELEMS: return a - b;
3077 case MOP_MUL_ELEMS: return a * b;
3078 case MOP_DIV_ELEMS: return a / b;
3079 case MOP_EXP_ELEMS: return pow (a, b);
3080 case MOP_GT: return a > b;
3081 case MOP_GE: return a >= b;
3082 case MOP_LT: return a < b;
3083 case MOP_LE: return a <= b;
3084 case MOP_EQ: return a == b;
3085 case MOP_NE: return a != b;
3086 case MOP_AND: return (a > 0) && (b > 0);
3087 case MOP_OR: return (a > 0) || (b > 0);
3088 case MOP_XOR: return (a > 0) != (b > 0);
3090 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3099 case MOP_PASTE_HORZ:
3100 case MOP_PASTE_VERT:
3116 matrix_op_name (enum matrix_op op)
3120 case MOP_ADD_ELEMS: return "+";
3121 case MOP_SUB_ELEMS: return "-";
3122 case MOP_MUL_ELEMS: return "&*";
3123 case MOP_DIV_ELEMS: return "&/";
3124 case MOP_EXP_ELEMS: return "&**";
3125 case MOP_GT: return ">";
3126 case MOP_GE: return ">=";
3127 case MOP_LT: return "<";
3128 case MOP_LE: return "<=";
3129 case MOP_EQ: return "=";
3130 case MOP_NE: return "<>";
3131 case MOP_AND: return "AND";
3132 case MOP_OR: return "OR";
3133 case MOP_XOR: return "XOR";
3135 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3144 case MOP_PASTE_HORZ:
3145 case MOP_PASTE_VERT:
3161 is_scalar (const gsl_matrix *m)
3163 return m->size1 == 1 && m->size2 == 1;
3167 to_scalar (const gsl_matrix *m)
3169 assert (is_scalar (m));
3170 return gsl_matrix_get (m, 0, 0);
3174 matrix_expr_evaluate_elementwise (enum matrix_op op,
3175 gsl_matrix *a, gsl_matrix *b)
3179 double be = to_scalar (b);
3180 for (size_t r = 0; r < a->size1; r++)
3181 for (size_t c = 0; c < a->size2; c++)
3183 double *ae = gsl_matrix_ptr (a, r, c);
3184 *ae = matrix_op_eval (op, *ae, be);
3188 else if (is_scalar (a))
3190 double ae = to_scalar (a);
3191 for (size_t r = 0; r < b->size1; r++)
3192 for (size_t c = 0; c < b->size2; c++)
3194 double *be = gsl_matrix_ptr (b, r, c);
3195 *be = matrix_op_eval (op, ae, *be);
3199 else if (a->size1 == b->size1 && a->size2 == b->size2)
3201 for (size_t r = 0; r < a->size1; r++)
3202 for (size_t c = 0; c < a->size2; c++)
3204 double *ae = gsl_matrix_ptr (a, r, c);
3205 double be = gsl_matrix_get (b, r, c);
3206 *ae = matrix_op_eval (op, *ae, be);
3212 msg (SE, _("Operands to %s must have the same dimensions or one "
3213 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
3214 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
3220 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
3222 if (is_scalar (a) || is_scalar (b))
3223 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
3225 if (a->size2 != b->size1)
3227 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
3228 "not conformable for multiplication."),
3229 a->size1, a->size2, b->size1, b->size2);
3233 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3234 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3239 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3241 gsl_matrix *tmp = *a;
3247 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3250 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3251 swap_matrix (z, tmp);
3255 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3257 mul_matrix (x, *x, *x, tmp);
3261 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
3264 if (x->size1 != x->size2)
3266 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
3267 "the left-hand size, not one with dimensions %zu×%zu."),
3268 x->size1, x->size2);
3273 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
3274 "right-hand side, not a matrix with dimensions %zu×%zu."),
3275 b->size1, b->size2);
3278 double bf = to_scalar (b);
3279 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3281 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
3282 "or outside the valid range."), bf);
3287 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3289 gsl_matrix_set_identity (y);
3293 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3295 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3298 mul_matrix (&y, x, y, &t);
3299 square_matrix (&x, &t);
3302 square_matrix (&x, &t);
3304 mul_matrix (&y, x, y, &t);
3308 /* Garbage collection.
3310 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3311 a permutation of them. We are returning one of them; that one must not be
3312 destroyed. We must not destroy 'x_' because the caller owns it. */
3314 gsl_matrix_free (y_);
3316 gsl_matrix_free (t_);
3322 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
3325 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3327 msg (SE, _("All operands of : operator must be scalars."));
3331 long int start = to_scalar (start_);
3332 long int end = to_scalar (end_);
3333 long int by = by_ ? to_scalar (by_) : 1;
3337 msg (SE, _("The increment operand to : must be nonzero."));
3341 long int n = (end >= start && by > 0 ? (end - start + by) / by
3342 : end <= start && by < 0 ? (start - end - by) / -by
3344 gsl_matrix *m = gsl_matrix_alloc (1, n);
3345 for (long int i = 0; i < n; i++)
3346 gsl_matrix_set (m, 0, i, start + i * by);
3351 matrix_expr_evaluate_not (gsl_matrix *a)
3353 for (size_t r = 0; r < a->size1; r++)
3354 for (size_t c = 0; c < a->size2; c++)
3356 double *ae = gsl_matrix_ptr (a, r, c);
3363 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
3365 if (a->size1 != b->size1)
3367 if (!a->size1 || !a->size2)
3369 else if (!b->size1 || !b->size2)
3372 msg (SE, _("All columns in a matrix must have the same number of rows, "
3373 "but this tries to paste matrices with %zu and %zu rows."),
3374 a->size1, b->size1);
3378 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3379 for (size_t y = 0; y < a->size1; y++)
3381 for (size_t x = 0; x < a->size2; x++)
3382 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3383 for (size_t x = 0; x < b->size2; x++)
3384 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3390 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
3392 if (a->size2 != b->size2)
3394 if (!a->size1 || !a->size2)
3396 else if (!b->size1 || !b->size2)
3399 msg (SE, _("All rows in a matrix must have the same number of columns, "
3400 "but this tries to stack matrices with %zu and %zu columns."),
3401 a->size2, b->size2);
3405 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3406 for (size_t x = 0; x < a->size2; x++)
3408 for (size_t y = 0; y < a->size1; y++)
3409 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3410 for (size_t y = 0; y < b->size1; y++)
3411 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3417 matrix_to_vector (gsl_matrix *m)
3420 gsl_vector v = to_vector (m);
3421 assert (v.block == m->block || !v.block);
3425 gsl_matrix_free (m);
3426 return xmemdup (&v, sizeof v);
3440 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3443 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
3444 enum index_type index_type, size_t other_size,
3445 struct index_vector *iv)
3454 msg (SE, _("Vector index must be scalar or vector, not a "
3456 m->size1, m->size2);
3460 msg (SE, _("Matrix row index must be scalar or vector, not a "
3462 m->size1, m->size2);
3466 msg (SE, _("Matrix column index must be scalar or vector, not a "
3468 m->size1, m->size2);
3474 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3475 *iv = (struct index_vector) {
3476 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3479 for (size_t i = 0; i < v.size; i++)
3481 double index = gsl_vector_get (&v, i);
3482 if (index < 1 || index >= size + 1)
3487 msg (SE, _("Index %g is out of range for vector "
3488 "with %zu elements."), index, size);
3492 msg (SE, _("%g is not a valid row index for "
3493 "a %zu×%zu matrix."),
3494 index, size, other_size);
3498 msg (SE, _("%g is not a valid column index for "
3499 "a %zu×%zu matrix."),
3500 index, other_size, size);
3507 iv->indexes[i] = index - 1;
3513 *iv = (struct index_vector) {
3514 .indexes = xnmalloc (size, sizeof *iv->indexes),
3517 for (size_t i = 0; i < size; i++)
3524 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
3526 if (!is_vector (sm))
3528 msg (SE, _("Vector index operator may not be applied to "
3529 "a %zu×%zu matrix."),
3530 sm->size1, sm->size2);
3538 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
3540 if (!matrix_expr_evaluate_vec_all (sm))
3543 gsl_vector sv = to_vector (sm);
3544 struct index_vector iv;
3545 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
3548 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3549 sm->size1 == 1 ? iv.n : 1);
3550 gsl_vector dv = to_vector (dm);
3551 for (size_t dx = 0; dx < iv.n; dx++)
3553 size_t sx = iv.indexes[dx];
3554 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3562 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
3565 struct index_vector iv0;
3566 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
3569 struct index_vector iv1;
3570 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
3577 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3578 for (size_t dy = 0; dy < iv0.n; dy++)
3580 size_t sy = iv0.indexes[dy];
3582 for (size_t dx = 0; dx < iv1.n; dx++)
3584 size_t sx = iv1.indexes[dx];
3585 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3593 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3594 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3595 const struct matrix_function_properties *, gsl_matrix *[], \
3596 const struct matrix_expr *, matrix_proto_##PROTO *);
3601 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
3603 if (!is_scalar (subs[index]))
3605 msg (SE, _("Function %s argument %zu must be a scalar, "
3606 "not a %zu×%zu matrix."),
3607 name, index + 1, subs[index]->size1, subs[index]->size2);
3614 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
3616 if (!is_vector (subs[index]))
3618 msg (SE, _("Function %s argument %zu must be a vector, "
3619 "not a %zu×%zu matrix."),
3620 name, index + 1, subs[index]->size1, subs[index]->size2);
3627 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
3629 for (size_t i = 0; i < n_subs; i++)
3631 if (!check_scalar_arg (name, subs, i))
3633 d[i] = to_scalar (subs[i]);
3639 parse_constraint_value (const char **constraintsp)
3642 long retval = strtol (*constraintsp, &tail, 10);
3643 assert (tail > *constraintsp);
3644 *constraintsp = tail;
3649 argument_invalid (const struct matrix_function_properties *props,
3650 size_t arg_index, gsl_matrix *a, size_t y, size_t x,
3654 ds_put_format (out, _("Argument %zu to matrix function %s "
3655 "has invalid value %g."),
3656 arg_index, props->name, gsl_matrix_get (a, y, x));
3658 ds_put_format (out, _("Row %zu, column %zu of argument %zu "
3659 "to matrix function %s has "
3660 "invalid value %g."),
3661 y + 1, x + 1, arg_index, props->name,
3662 gsl_matrix_get (a, y, x));
3665 enum matrix_argument_relop
3675 argument_inequality_error (const struct matrix_function_properties *props,
3676 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3677 size_t b_index, double b,
3678 enum matrix_argument_relop relop)
3680 struct string s = DS_EMPTY_INITIALIZER;
3682 argument_invalid (props, a_index, a, y, x, &s);
3683 ds_put_cstr (&s, " ");
3687 ds_put_format (&s, _("This argument must be greater than or "
3688 "equal to argument %zu, but that has value %g."),
3693 ds_put_format (&s, _("This argument must be greater than argument %zu, "
3694 "but that has value %g."),
3699 ds_put_format (&s, _("This argument must be less than or "
3700 "equal to argument %zu, but that has value %g."),
3705 ds_put_format (&s, _("This argument must be less than argument %zu, "
3706 "but that has value %g."),
3712 _("This argument must not be the same as argument %zu."),
3717 msg (ME, ds_cstr (&s));
3722 argument_value_error (const struct matrix_function_properties *props,
3723 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3725 enum matrix_argument_relop relop)
3727 struct string s = DS_EMPTY_INITIALIZER;
3728 argument_invalid (props, a_index, a, y, x, &s);
3729 ds_put_cstr (&s, " ");
3734 _("This argument must be greater than or equal to %g."),
3739 ds_put_format (&s, _("This argument must be greater than %g."), b);
3743 ds_put_format (&s, _("This argument must be less than or equal to %g."),
3748 ds_put_format (&s, _("This argument must be less than %g."), b);
3752 ds_put_format (&s, _("This argument may not be %g."), b);
3756 msg (ME, ds_cstr (&s));
3761 matrix_argument_relop_is_satisfied (double a, double b,
3762 enum matrix_argument_relop relop)
3766 case MRR_GE: return a >= b;
3767 case MRR_GT: return a > b;
3768 case MRR_LE: return a <= b;
3769 case MRR_LT: return a < b;
3770 case MRR_NE: return a != b;
3776 static enum matrix_argument_relop
3777 matrix_argument_relop_flip (enum matrix_argument_relop relop)
3781 case MRR_GE: return MRR_LE;
3782 case MRR_GT: return MRR_LT;
3783 case MRR_LE: return MRR_GE;
3784 case MRR_LT: return MRR_GT;
3785 case MRR_NE: return MRR_NE;
3792 check_constraints (const struct matrix_function_properties *props,
3793 gsl_matrix *args[], size_t n_args)
3795 const char *constraints = props->constraints;
3799 size_t arg_index = SIZE_MAX;
3800 while (*constraints)
3802 if (*constraints >= 'a' && *constraints <= 'd')
3804 arg_index = *constraints++ - 'a';
3805 assert (arg_index < n_args);
3807 else if (*constraints == '[' || *constraints == '(')
3809 assert (arg_index < n_args);
3810 bool open_lower = *constraints++ == '(';
3811 int minimum = parse_constraint_value (&constraints);
3812 assert (*constraints == ',');
3814 int maximum = parse_constraint_value (&constraints);
3815 assert (*constraints == ']' || *constraints == ')');
3816 bool open_upper = *constraints++ == ')';
3818 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3819 if ((open_lower ? *d <= minimum : *d < minimum)
3820 || (open_upper ? *d >= maximum : *d > maximum))
3822 if (!is_scalar (args[arg_index]))
3823 msg (ME, _("Row %zu, column %zu of argument %zu to matrix "
3824 "function %s has value %g, which is outside "
3825 "the valid range %c%d,%d%c."),
3826 y + 1, x + 1, arg_index + 1, props->name, *d,
3827 open_lower ? '(' : '[',
3829 open_upper ? ')' : ']');
3831 msg (ME, _("Argument %zu to matrix function %s has value %g, "
3832 "which is outside the valid range %c%d,%d%c."),
3833 arg_index + 1, props->name, *d,
3834 open_lower ? '(' : '[',
3836 open_upper ? ')' : ']');
3840 else if (*constraints == '>'
3841 || *constraints == '<'
3842 || *constraints == '!')
3844 enum matrix_argument_relop relop;
3845 switch (*constraints++)
3848 if (*constraints == '=')
3858 if (*constraints == '=')
3868 assert (*constraints == '=');
3873 if (*constraints >= 'a' && *constraints <= 'd')
3875 size_t a_index = arg_index;
3876 size_t b_index = *constraints - 'a';
3877 assert (a_index < n_args);
3878 assert (b_index < n_args);
3880 /* We only support one of the two arguments being non-scalar.
3881 It's easier to support only the first one being non-scalar, so
3882 flip things around if it's the other way. */
3883 if (!is_scalar (args[b_index]))
3885 assert (is_scalar (args[a_index]));
3886 size_t tmp_index = a_index;
3888 b_index = tmp_index;
3889 relop = matrix_argument_relop_flip (relop);
3892 double b = to_scalar (args[b_index]);
3893 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
3894 if (!matrix_argument_relop_is_satisfied (*a, b, relop))
3896 argument_inequality_error (props,
3897 a_index + 1, args[a_index], y, x,
3905 int comparand = parse_constraint_value (&constraints);
3907 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3908 if (!matrix_argument_relop_is_satisfied (*d, comparand, relop))
3910 argument_value_error (props,
3911 arg_index + 1, args[arg_index], y, x,
3920 assert (*constraints == ' ');
3922 arg_index = SIZE_MAX;
3929 matrix_expr_evaluate_d_none (
3930 const struct matrix_function_properties *props UNUSED,
3931 gsl_matrix *subs[] UNUSED, const struct matrix_expr *e,
3932 matrix_proto_d_none *f)
3934 assert (e->n_subs == 0);
3936 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3937 gsl_matrix_set (m, 0, 0, f ());
3942 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
3943 gsl_matrix *subs[], const struct matrix_expr *e,
3944 matrix_proto_d_d *f)
3946 assert (e->n_subs == 1);
3949 if (!to_scalar_args (props->name, subs, e->n_subs, &d))
3952 if (!check_constraints (props, subs, e->n_subs))
3955 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3956 gsl_matrix_set (m, 0, 0, f (d));
3961 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
3962 gsl_matrix *subs[], const struct matrix_expr *e,
3963 matrix_proto_d_dd *f)
3965 assert (e->n_subs == 2);
3968 if (!to_scalar_args (props->name, subs, e->n_subs, d))
3971 if (!check_constraints (props, subs, e->n_subs))
3974 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3975 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
3980 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
3981 gsl_matrix *subs[], const struct matrix_expr *e,
3982 matrix_proto_d_ddd *f)
3984 assert (e->n_subs == 3);
3987 if (!to_scalar_args (props->name, subs, e->n_subs, d))
3990 if (!check_constraints (props, subs, e->n_subs))
3993 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3994 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
3999 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
4000 gsl_matrix *subs[], const struct matrix_expr *e,
4001 matrix_proto_m_d *f)
4003 assert (e->n_subs == 1);
4006 return to_scalar_args (props->name, subs, e->n_subs, &d) ? f (d) : NULL;
4010 matrix_expr_evaluate_m_dd (const struct matrix_function_properties *props,
4011 gsl_matrix *subs[], const struct matrix_expr *e,
4012 matrix_proto_m_dd *f)
4014 assert (e->n_subs == 2);
4017 return to_scalar_args (props->name, subs, e->n_subs, d) ? f(d[0], d[1]) : NULL;
4021 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
4022 gsl_matrix *subs[], const struct matrix_expr *e,
4023 matrix_proto_m_ddd *f)
4025 assert (e->n_subs == 3);
4028 return to_scalar_args (props->name, subs, e->n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
4032 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED,
4033 gsl_matrix *subs[], const struct matrix_expr *e,
4034 matrix_proto_m_m *f)
4036 assert (e->n_subs == 1);
4041 matrix_expr_evaluate_m_me (const struct matrix_function_properties *props UNUSED,
4042 gsl_matrix *subs[], const struct matrix_expr *e,
4043 matrix_proto_m_me *f)
4045 assert (e->n_subs == 1);
4046 return f (subs[0], e);
4050 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
4051 gsl_matrix *subs[], const struct matrix_expr *e,
4052 matrix_proto_m_e *f)
4054 assert (e->n_subs == 1);
4056 if (!check_constraints (props, subs, e->n_subs))
4059 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4065 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED,
4066 gsl_matrix *subs[], const struct matrix_expr *e,
4067 matrix_proto_m_md *f)
4069 assert (e->n_subs == 2);
4070 if (!check_scalar_arg (props->name, subs, 1))
4072 return f (subs[0], to_scalar (subs[1]));
4076 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
4077 gsl_matrix *subs[], const struct matrix_expr *e,
4078 matrix_proto_m_ed *f)
4080 assert (e->n_subs == 2);
4081 if (!check_scalar_arg (props->name, subs, 1))
4084 if (!check_constraints (props, subs, e->n_subs))
4087 double b = to_scalar (subs[1]);
4088 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4094 matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props UNUSED,
4095 gsl_matrix *subs[], const struct matrix_expr *e,
4096 matrix_proto_m_mdd *f)
4098 assert (e->n_subs == 3);
4099 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
4101 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
4105 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
4106 gsl_matrix *subs[], const struct matrix_expr *e,
4107 matrix_proto_m_edd *f)
4109 assert (e->n_subs == 3);
4110 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
4113 if (!check_constraints (props, subs, e->n_subs))
4116 double b = to_scalar (subs[1]);
4117 double c = to_scalar (subs[2]);
4118 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4124 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4125 gsl_matrix *subs[], const struct matrix_expr *e,
4126 matrix_proto_m_eddd *f)
4128 assert (e->n_subs == 4);
4129 for (size_t i = 1; i < 4; i++)
4130 if (!check_scalar_arg (props->name, subs, i))
4133 if (!check_constraints (props, subs, e->n_subs))
4136 double b = to_scalar (subs[1]);
4137 double c = to_scalar (subs[2]);
4138 double d = to_scalar (subs[3]);
4139 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4140 *a = f (*a, b, c, d);
4145 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4146 gsl_matrix *subs[], const struct matrix_expr *e,
4147 matrix_proto_m_eed *f)
4149 assert (e->n_subs == 3);
4150 if (!check_scalar_arg (props->name, subs, 2))
4153 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4154 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4156 msg (ME, _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4157 "%zu×%zu, but %s requires these arguments either to have "
4158 "the same dimensions or for one of them to be a scalar."),
4160 subs[0]->size1, subs[0]->size2,
4161 subs[1]->size1, subs[1]->size2,
4166 if (!check_constraints (props, subs, e->n_subs))
4169 double c = to_scalar (subs[2]);
4171 if (is_scalar (subs[0]))
4173 double a = to_scalar (subs[0]);
4174 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4180 double b = to_scalar (subs[1]);
4181 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4188 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED,
4189 gsl_matrix *subs[], const struct matrix_expr *e,
4190 matrix_proto_m_mm *f)
4192 assert (e->n_subs == 2);
4193 return f (subs[0], subs[1]);
4197 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4198 gsl_matrix *subs[], const struct matrix_expr *e,
4199 matrix_proto_m_v *f)
4201 assert (e->n_subs == 1);
4202 if (!check_vector_arg (props->name, subs, 0))
4204 gsl_vector v = to_vector (subs[0]);
4209 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED,
4210 gsl_matrix *subs[], const struct matrix_expr *e,
4211 matrix_proto_d_m *f)
4213 assert (e->n_subs == 1);
4215 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4216 gsl_matrix_set (m, 0, 0, f (subs[0]));
4221 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED,
4222 gsl_matrix *subs[], const struct matrix_expr *e,
4223 matrix_proto_m_any *f)
4225 return f (subs, e->n_subs);
4229 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props,
4230 gsl_matrix *subs[], const struct matrix_expr *e,
4231 matrix_proto_IDENT *f)
4233 assert (e->n_subs <= 2);
4236 if (!to_scalar_args (props->name, subs, e->n_subs, d))
4239 return f (d[0], d[e->n_subs - 1]);
4243 matrix_expr_evaluate (const struct matrix_expr *e)
4245 if (e->op == MOP_NUMBER)
4247 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4248 gsl_matrix_set (m, 0, 0, e->number);
4251 else if (e->op == MOP_VARIABLE)
4253 const gsl_matrix *src = e->variable->value;
4256 msg (SE, _("Uninitialized variable %s used in expression."),
4261 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4262 gsl_matrix_memcpy (dst, src);
4265 else if (e->op == MOP_EOF)
4267 struct dfm_reader *reader = read_file_open (e->eof);
4268 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4269 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4273 enum { N_LOCAL = 3 };
4274 gsl_matrix *local_subs[N_LOCAL];
4275 gsl_matrix **subs = (e->n_subs < N_LOCAL
4277 : xmalloc (e->n_subs * sizeof *subs));
4279 for (size_t i = 0; i < e->n_subs; i++)
4281 subs[i] = matrix_expr_evaluate (e->subs[i]);
4284 for (size_t j = 0; j < i; j++)
4285 gsl_matrix_free (subs[j]);
4286 if (subs != local_subs)
4292 gsl_matrix *result = NULL;
4295 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4296 case MOP_F_##ENUM: \
4298 static const struct matrix_function_properties props = { \
4300 .constraints = CONSTRAINTS, \
4302 result = matrix_expr_evaluate_##PROTO (&props, subs, e, \
4303 matrix_eval_##ENUM); \
4310 gsl_matrix_scale (subs[0], -1.0);
4328 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
4332 result = matrix_expr_evaluate_not (subs[0]);
4336 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
4340 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
4344 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
4348 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
4351 case MOP_PASTE_HORZ:
4352 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
4355 case MOP_PASTE_VERT:
4356 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
4360 result = gsl_matrix_alloc (0, 0);
4364 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
4368 result = matrix_expr_evaluate_vec_all (subs[0]);
4372 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
4376 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
4380 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
4389 for (size_t i = 0; i < e->n_subs; i++)
4390 if (subs[i] != result)
4391 gsl_matrix_free (subs[i]);
4392 if (subs != local_subs)
4398 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4401 gsl_matrix *m = matrix_expr_evaluate (e);
4407 msg (SE, _("Expression for %s must evaluate to scalar, "
4408 "not a %zu×%zu matrix."),
4409 context, m->size1, m->size2);
4410 gsl_matrix_free (m);
4415 gsl_matrix_free (m);
4420 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4424 if (!matrix_expr_evaluate_scalar (e, context, &d))
4428 if (d < LONG_MIN || d > LONG_MAX)
4430 msg (SE, _("Expression for %s is outside the integer range."), context);
4437 struct matrix_lvalue
4439 struct matrix_var *var;
4440 struct matrix_expr *indexes[2];
4443 struct msg_location *var_location; /* Variable name. */
4444 struct msg_location *index_location; /* Variable name plus indexing. */
4448 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4452 msg_location_destroy (lvalue->var_location);
4453 msg_location_destroy (lvalue->index_location);
4454 for (size_t i = 0; i < lvalue->n_indexes; i++)
4455 matrix_expr_destroy (lvalue->indexes[i]);
4460 static struct matrix_lvalue *
4461 matrix_lvalue_parse (struct matrix_state *s)
4463 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4464 if (!lex_force_id (s->lexer))
4466 int start_ofs = lex_ofs (s->lexer);
4467 lvalue->var_location = lex_get_location (s->lexer, 0, 0);
4468 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4469 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4473 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
4477 lex_get_n (s->lexer, 2);
4479 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
4481 lvalue->n_indexes++;
4483 if (lex_match (s->lexer, T_COMMA))
4485 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
4487 lvalue->n_indexes++;
4489 if (!lex_force_match (s->lexer, T_RPAREN))
4492 lvalue->index_location = lex_ofs_location (s->lexer, start_ofs,
4493 lex_ofs (s->lexer) - 1);
4498 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4504 matrix_lvalue_destroy (lvalue);
4509 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4510 enum index_type index_type, size_t other_size,
4511 struct index_vector *iv)
4516 m = matrix_expr_evaluate (e);
4523 bool ok = matrix_normalize_index_vector (m, size, index_type,
4525 gsl_matrix_free (m);
4530 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4531 struct index_vector *iv, gsl_matrix *sm)
4533 gsl_vector dv = to_vector (lvalue->var->value);
4535 /* Convert source matrix 'sm' to source vector 'sv'. */
4536 if (!is_vector (sm))
4538 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
4539 sm->size1, sm->size2);
4542 gsl_vector sv = to_vector (sm);
4543 if (iv->n != sv.size)
4545 msg (SE, _("Can't assign %zu-element vector "
4546 "to %zu-element subvector."), sv.size, iv->n);
4550 /* Assign elements. */
4551 for (size_t x = 0; x < iv->n; x++)
4552 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4557 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4558 struct index_vector *iv0,
4559 struct index_vector *iv1, gsl_matrix *sm)
4561 gsl_matrix *dm = lvalue->var->value;
4563 /* Convert source matrix 'sm' to source vector 'sv'. */
4564 if (iv0->n != sm->size1)
4566 msg (SE, _("Row index vector for assignment to %s has %zu elements "
4567 "but source matrix has %zu rows."),
4568 lvalue->var->name, iv0->n, sm->size1);
4571 if (iv1->n != sm->size2)
4573 msg (SE, _("Column index vector for assignment to %s has %zu elements "
4574 "but source matrix has %zu columns."),
4575 lvalue->var->name, iv1->n, sm->size2);
4579 /* Assign elements. */
4580 for (size_t y = 0; y < iv0->n; y++)
4582 size_t f0 = iv0->indexes[y];
4583 for (size_t x = 0; x < iv1->n; x++)
4585 size_t f1 = iv1->indexes[x];
4586 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4592 /* Takes ownership of M. */
4594 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4595 struct index_vector *iv0, struct index_vector *iv1,
4598 if (!lvalue->n_indexes)
4600 gsl_matrix_free (lvalue->var->value);
4601 lvalue->var->value = sm;
4606 bool ok = (lvalue->n_indexes == 1
4607 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
4608 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
4609 free (iv0->indexes);
4610 free (iv1->indexes);
4611 gsl_matrix_free (sm);
4617 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4618 struct index_vector *iv0,
4619 struct index_vector *iv1)
4621 *iv0 = INDEX_VECTOR_INIT;
4622 *iv1 = INDEX_VECTOR_INIT;
4623 if (!lvalue->n_indexes)
4626 /* Validate destination matrix exists and has the right shape. */
4627 gsl_matrix *dm = lvalue->var->value;
4630 msg_at (SE, lvalue->var_location,
4631 _("Undefined variable %s."), lvalue->var->name);
4634 else if (dm->size1 == 0 || dm->size2 == 0)
4636 msg_at (SE, lvalue->index_location,
4637 _("Cannot index %zu×%zu matrix."), dm->size1, dm->size2);
4640 else if (lvalue->n_indexes == 1)
4642 if (!is_vector (dm))
4644 msg_at (SE, lvalue->index_location,
4645 _("Can't use vector indexing on %zu×%zu matrix %s."),
4646 dm->size1, dm->size2, lvalue->var->name);
4649 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4650 MAX (dm->size1, dm->size2),
4655 assert (lvalue->n_indexes == 2);
4656 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4657 IV_ROW, dm->size2, iv0))
4660 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4661 IV_COLUMN, dm->size1, iv1))
4663 free (iv0->indexes);
4670 /* Takes ownership of M. */
4672 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
4674 struct index_vector iv0, iv1;
4675 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
4677 gsl_matrix_free (sm);
4681 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
4685 /* Matrix command. */
4689 struct matrix_cmd **commands;
4693 static void matrix_cmds_uninit (struct matrix_cmds *);
4697 enum matrix_cmd_type
4720 struct compute_command
4722 struct matrix_lvalue *lvalue;
4723 struct matrix_expr *rvalue;
4727 struct print_command
4729 struct matrix_expr *expression;
4730 bool use_default_format;
4731 struct fmt_spec format;
4733 int space; /* -1 means NEWPAGE. */
4737 struct string_array literals; /* CLABELS/RLABELS. */
4738 struct matrix_expr *expr; /* CNAMES/RNAMES. */
4744 struct do_if_command
4748 struct matrix_expr *condition;
4749 struct matrix_cmds commands;
4759 struct matrix_var *var;
4760 struct matrix_expr *start;
4761 struct matrix_expr *end;
4762 struct matrix_expr *increment;
4764 /* Loop conditions. */
4765 struct matrix_expr *top_condition;
4766 struct matrix_expr *bottom_condition;
4769 struct matrix_cmds commands;
4773 struct display_command
4775 struct matrix_state *state;
4779 struct release_command
4781 struct matrix_var **vars;
4788 struct matrix_expr *expression;
4789 struct save_file *sf;
4795 struct read_file *rf;
4796 struct matrix_lvalue *dst;
4797 struct matrix_expr *size;
4799 enum fmt_type format;
4806 struct write_command
4808 struct write_file *wf;
4809 struct matrix_expr *expression;
4812 /* If this is nonnull, WRITE uses this format.
4814 If this is NULL, WRITE uses free-field format with as many
4815 digits of precision as needed. */
4816 struct fmt_spec *format;
4825 struct matrix_lvalue *dst;
4826 struct dataset *dataset;
4827 struct file_handle *file;
4829 struct var_syntax *vars;
4831 struct matrix_var *names;
4833 /* Treatment of missing values. */
4838 MGET_ERROR, /* Flag error on command. */
4839 MGET_ACCEPT, /* Accept without change, user-missing only. */
4840 MGET_OMIT, /* Drop this case. */
4841 MGET_RECODE /* Recode to 'substitute'. */
4844 double substitute; /* MGET_RECODE only. */
4850 struct msave_command
4852 struct msave_common *common;
4853 struct matrix_expr *expr;
4854 const char *rowtype;
4855 const struct matrix_expr *factors;
4856 const struct matrix_expr *splits;
4862 struct matrix_state *state;
4863 struct file_handle *file;
4865 struct stringi_set rowtypes;
4869 struct eigen_command
4871 struct matrix_expr *expr;
4872 struct matrix_var *evec;
4873 struct matrix_var *eval;
4877 struct setdiag_command
4879 struct matrix_var *dst;
4880 struct matrix_expr *expr;
4886 struct matrix_expr *expr;
4887 struct matrix_var *u;
4888 struct matrix_var *s;
4889 struct matrix_var *v;
4895 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
4896 static bool matrix_cmd_execute (struct matrix_cmd *);
4897 static void matrix_cmd_destroy (struct matrix_cmd *);
4900 static struct matrix_cmd *
4901 matrix_parse_compute (struct matrix_state *s)
4903 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4904 *cmd = (struct matrix_cmd) {
4905 .type = MCMD_COMPUTE,
4906 .compute = { .lvalue = NULL },
4909 struct compute_command *compute = &cmd->compute;
4910 compute->lvalue = matrix_lvalue_parse (s);
4911 if (!compute->lvalue)
4914 if (!lex_force_match (s->lexer, T_EQUALS))
4917 compute->rvalue = matrix_parse_expr (s);
4918 if (!compute->rvalue)
4924 matrix_cmd_destroy (cmd);
4929 matrix_cmd_execute_compute (struct compute_command *compute)
4931 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
4935 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
4939 print_labels_destroy (struct print_labels *labels)
4943 string_array_destroy (&labels->literals);
4944 matrix_expr_destroy (labels->expr);
4950 parse_literal_print_labels (struct matrix_state *s,
4951 struct print_labels **labelsp)
4953 lex_match (s->lexer, T_EQUALS);
4954 print_labels_destroy (*labelsp);
4955 *labelsp = xzalloc (sizeof **labelsp);
4956 while (lex_token (s->lexer) != T_SLASH
4957 && lex_token (s->lexer) != T_ENDCMD
4958 && lex_token (s->lexer) != T_STOP)
4960 struct string label = DS_EMPTY_INITIALIZER;
4961 while (lex_token (s->lexer) != T_COMMA
4962 && lex_token (s->lexer) != T_SLASH
4963 && lex_token (s->lexer) != T_ENDCMD
4964 && lex_token (s->lexer) != T_STOP)
4966 if (!ds_is_empty (&label))
4967 ds_put_byte (&label, ' ');
4969 if (lex_token (s->lexer) == T_STRING)
4970 ds_put_cstr (&label, lex_tokcstr (s->lexer));
4973 char *rep = lex_next_representation (s->lexer, 0, 0);
4974 ds_put_cstr (&label, rep);
4979 string_array_append_nocopy (&(*labelsp)->literals,
4980 ds_steal_cstr (&label));
4982 if (!lex_match (s->lexer, T_COMMA))
4988 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
4990 lex_match (s->lexer, T_EQUALS);
4991 struct matrix_expr *e = matrix_parse_exp (s);
4995 print_labels_destroy (*labelsp);
4996 *labelsp = xzalloc (sizeof **labelsp);
4997 (*labelsp)->expr = e;
5001 static struct matrix_cmd *
5002 matrix_parse_print (struct matrix_state *s)
5004 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5005 *cmd = (struct matrix_cmd) {
5008 .use_default_format = true,
5012 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
5015 for (size_t i = 0; ; i++)
5017 enum token_type t = lex_next_token (s->lexer, i);
5018 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
5020 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
5022 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
5025 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
5030 cmd->print.expression = matrix_parse_exp (s);
5031 if (!cmd->print.expression)
5035 while (lex_match (s->lexer, T_SLASH))
5037 if (lex_match_id (s->lexer, "FORMAT"))
5039 lex_match (s->lexer, T_EQUALS);
5040 if (!parse_format_specifier (s->lexer, &cmd->print.format))
5042 cmd->print.use_default_format = false;
5044 else if (lex_match_id (s->lexer, "TITLE"))
5046 lex_match (s->lexer, T_EQUALS);
5047 if (!lex_force_string (s->lexer))
5049 free (cmd->print.title);
5050 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
5053 else if (lex_match_id (s->lexer, "SPACE"))
5055 lex_match (s->lexer, T_EQUALS);
5056 if (lex_match_id (s->lexer, "NEWPAGE"))
5057 cmd->print.space = -1;
5058 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
5060 cmd->print.space = lex_integer (s->lexer);
5066 else if (lex_match_id (s->lexer, "RLABELS"))
5067 parse_literal_print_labels (s, &cmd->print.rlabels);
5068 else if (lex_match_id (s->lexer, "CLABELS"))
5069 parse_literal_print_labels (s, &cmd->print.clabels);
5070 else if (lex_match_id (s->lexer, "RNAMES"))
5072 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
5075 else if (lex_match_id (s->lexer, "CNAMES"))
5077 if (!parse_expr_print_labels (s, &cmd->print.clabels))
5082 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
5083 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
5091 matrix_cmd_destroy (cmd);
5096 matrix_is_integer (const gsl_matrix *m)
5098 for (size_t y = 0; y < m->size1; y++)
5099 for (size_t x = 0; x < m->size2; x++)
5101 double d = gsl_matrix_get (m, y, x);
5109 matrix_max_magnitude (const gsl_matrix *m)
5112 for (size_t y = 0; y < m->size1; y++)
5113 for (size_t x = 0; x < m->size2; x++)
5115 double d = fabs (gsl_matrix_get (m, y, x));
5123 format_fits (struct fmt_spec format, double x)
5125 char *s = data_out (&(union value) { .f = x }, NULL,
5126 &format, settings_get_fmt_settings ());
5127 bool fits = *s != '*' && !strchr (s, 'E');
5132 static struct fmt_spec
5133 default_format (const gsl_matrix *m, int *log_scale)
5137 double max = matrix_max_magnitude (m);
5139 if (matrix_is_integer (m))
5140 for (int w = 1; w <= 12; w++)
5142 struct fmt_spec format = { .type = FMT_F, .w = w };
5143 if (format_fits (format, -max))
5147 if (max >= 1e9 || max <= 1e-4)
5150 snprintf (s, sizeof s, "%.1e", max);
5152 const char *e = strchr (s, 'e');
5154 *log_scale = atoi (e + 1);
5157 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5161 trimmed_string (double d)
5163 struct substring s = ss_buffer ((char *) &d, sizeof d);
5164 ss_rtrim (&s, ss_cstr (" "));
5165 return ss_xstrdup (s);
5168 static struct string_array *
5169 print_labels_get (const struct print_labels *labels, size_t n,
5170 const char *prefix, bool truncate)
5175 struct string_array *out = xzalloc (sizeof *out);
5176 if (labels->literals.n)
5177 string_array_clone (out, &labels->literals);
5178 else if (labels->expr)
5180 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5181 if (m && is_vector (m))
5183 gsl_vector v = to_vector (m);
5184 for (size_t i = 0; i < v.size; i++)
5185 string_array_append_nocopy (out, trimmed_string (
5186 gsl_vector_get (&v, i)));
5188 gsl_matrix_free (m);
5194 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5197 string_array_append (out, "");
5200 string_array_delete (out, out->n - 1);
5203 for (size_t i = 0; i < out->n; i++)
5205 char *s = out->strings[i];
5206 s[strnlen (s, 8)] = '\0';
5213 matrix_cmd_print_space (int space)
5216 output_item_submit (page_break_item_create ());
5217 for (int i = 0; i < space; i++)
5218 output_log ("%s", "");
5222 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
5223 struct fmt_spec format, int log_scale)
5225 matrix_cmd_print_space (print->space);
5227 output_log ("%s", print->title);
5229 output_log (" 10 ** %d X", log_scale);
5231 struct string_array *clabels = print_labels_get (print->clabels,
5232 m->size2, "col", true);
5233 if (clabels && format.w < 8)
5235 struct string_array *rlabels = print_labels_get (print->rlabels,
5236 m->size1, "row", true);
5240 struct string line = DS_EMPTY_INITIALIZER;
5242 ds_put_byte_multiple (&line, ' ', 8);
5243 for (size_t x = 0; x < m->size2; x++)
5244 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5245 output_log_nocopy (ds_steal_cstr (&line));
5248 double scale = pow (10.0, log_scale);
5249 bool numeric = fmt_is_numeric (format.type);
5250 for (size_t y = 0; y < m->size1; y++)
5252 struct string line = DS_EMPTY_INITIALIZER;
5254 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5256 for (size_t x = 0; x < m->size2; x++)
5258 double f = gsl_matrix_get (m, y, x);
5260 ? data_out (&(union value) { .f = f / scale}, NULL,
5261 &format, settings_get_fmt_settings ())
5262 : trimmed_string (f));
5263 ds_put_format (&line, " %s", s);
5266 output_log_nocopy (ds_steal_cstr (&line));
5269 string_array_destroy (rlabels);
5271 string_array_destroy (clabels);
5276 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5277 const struct print_labels *print_labels, size_t n,
5278 const char *name, const char *prefix)
5280 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5282 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5283 for (size_t i = 0; i < n; i++)
5284 pivot_category_create_leaf (
5286 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5287 : pivot_value_new_integer (i + 1)));
5289 d->hide_all_labels = true;
5290 string_array_destroy (labels);
5295 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
5296 struct fmt_spec format, int log_scale)
5298 struct pivot_table *table = pivot_table_create__ (
5299 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5302 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5304 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5305 N_("Columns"), "col");
5307 struct pivot_footnote *footnote = NULL;
5310 char *s = xasprintf ("× 10**%d", log_scale);
5311 footnote = pivot_table_create_footnote (
5312 table, pivot_value_new_user_text_nocopy (s));
5315 double scale = pow (10.0, log_scale);
5316 bool numeric = fmt_is_numeric (format.type);
5317 for (size_t y = 0; y < m->size1; y++)
5318 for (size_t x = 0; x < m->size2; x++)
5320 double f = gsl_matrix_get (m, y, x);
5321 struct pivot_value *v;
5324 v = pivot_value_new_number (f / scale);
5325 v->numeric.format = format;
5328 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5330 pivot_value_add_footnote (v, footnote);
5331 pivot_table_put2 (table, y, x, v);
5334 pivot_table_submit (table);
5338 matrix_cmd_execute_print (const struct print_command *print)
5340 if (print->expression)
5342 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5347 struct fmt_spec format = (print->use_default_format
5348 ? default_format (m, &log_scale)
5351 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5352 matrix_cmd_print_text (print, m, format, log_scale);
5354 matrix_cmd_print_tables (print, m, format, log_scale);
5356 gsl_matrix_free (m);
5360 matrix_cmd_print_space (print->space);
5362 output_log ("%s", print->title);
5369 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
5370 const char *command_name,
5371 const char *stop1, const char *stop2)
5373 lex_end_of_command (s->lexer);
5374 lex_discard_rest_of_command (s->lexer);
5376 size_t allocated = 0;
5379 while (lex_token (s->lexer) == T_ENDCMD)
5382 if (lex_at_phrase (s->lexer, stop1)
5383 || (stop2 && lex_at_phrase (s->lexer, stop2)))
5386 if (lex_at_phrase (s->lexer, "END MATRIX"))
5388 msg (SE, _("Premature END MATRIX within %s."), command_name);
5392 struct matrix_cmd *cmd = matrix_parse_command (s);
5396 if (c->n >= allocated)
5397 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
5398 c->commands[c->n++] = cmd;
5403 matrix_parse_do_if_clause (struct matrix_state *s,
5404 struct do_if_command *ifc,
5405 bool parse_condition,
5406 size_t *allocated_clauses)
5408 if (ifc->n_clauses >= *allocated_clauses)
5409 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5410 sizeof *ifc->clauses);
5411 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5412 *c = (struct do_if_clause) { .condition = NULL };
5414 if (parse_condition)
5416 c->condition = matrix_parse_expr (s);
5421 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
5424 static struct matrix_cmd *
5425 matrix_parse_do_if (struct matrix_state *s)
5427 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5428 *cmd = (struct matrix_cmd) {
5430 .do_if = { .n_clauses = 0 }
5433 size_t allocated_clauses = 0;
5436 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
5439 while (lex_match_phrase (s->lexer, "ELSE IF"));
5441 if (lex_match_id (s->lexer, "ELSE")
5442 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
5445 if (!lex_match_phrase (s->lexer, "END IF"))
5450 matrix_cmd_destroy (cmd);
5455 matrix_cmd_execute_do_if (struct do_if_command *cmd)
5457 for (size_t i = 0; i < cmd->n_clauses; i++)
5459 struct do_if_clause *c = &cmd->clauses[i];
5463 if (!matrix_expr_evaluate_scalar (c->condition,
5464 i ? "ELSE IF" : "DO IF",
5469 for (size_t j = 0; j < c->commands.n; j++)
5470 if (!matrix_cmd_execute (c->commands.commands[j]))
5477 static struct matrix_cmd *
5478 matrix_parse_loop (struct matrix_state *s)
5480 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5481 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5483 struct loop_command *loop = &cmd->loop;
5484 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5486 struct substring name = lex_tokss (s->lexer);
5487 loop->var = matrix_var_lookup (s, name);
5489 loop->var = matrix_var_create (s, name);
5494 loop->start = matrix_parse_expr (s);
5495 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5498 loop->end = matrix_parse_expr (s);
5502 if (lex_match (s->lexer, T_BY))
5504 loop->increment = matrix_parse_expr (s);
5505 if (!loop->increment)
5510 if (lex_match_id (s->lexer, "IF"))
5512 loop->top_condition = matrix_parse_expr (s);
5513 if (!loop->top_condition)
5517 bool was_in_loop = s->in_loop;
5519 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
5521 s->in_loop = was_in_loop;
5525 if (!lex_match_phrase (s->lexer, "END LOOP"))
5528 if (lex_match_id (s->lexer, "IF"))
5530 loop->bottom_condition = matrix_parse_expr (s);
5531 if (!loop->bottom_condition)
5538 matrix_cmd_destroy (cmd);
5543 matrix_cmd_execute_loop (struct loop_command *cmd)
5547 long int increment = 1;
5550 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5551 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5553 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5557 if (increment > 0 ? value > end
5558 : increment < 0 ? value < end
5563 int mxloops = settings_get_mxloops ();
5564 for (int i = 0; i < mxloops; i++)
5569 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5571 gsl_matrix_free (cmd->var->value);
5572 cmd->var->value = NULL;
5574 if (!cmd->var->value)
5575 cmd->var->value = gsl_matrix_alloc (1, 1);
5577 gsl_matrix_set (cmd->var->value, 0, 0, value);
5580 if (cmd->top_condition)
5583 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5588 for (size_t j = 0; j < cmd->commands.n; j++)
5589 if (!matrix_cmd_execute (cmd->commands.commands[j]))
5592 if (cmd->bottom_condition)
5595 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5596 "END LOOP IF", &d) || d > 0)
5603 if (increment > 0 ? value > end : value < end)
5609 static struct matrix_cmd *
5610 matrix_parse_break (struct matrix_state *s)
5614 msg (SE, _("BREAK not inside LOOP."));
5618 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5619 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
5623 static struct matrix_cmd *
5624 matrix_parse_display (struct matrix_state *s)
5626 while (lex_token (s->lexer) != T_ENDCMD)
5628 if (!lex_match_id (s->lexer, "DICTIONARY")
5629 && !lex_match_id (s->lexer, "STATUS"))
5631 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5636 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5637 *cmd = (struct matrix_cmd) { .type = MCMD_DISPLAY, .display = { s } };
5642 compare_matrix_var_pointers (const void *a_, const void *b_)
5644 const struct matrix_var *const *ap = a_;
5645 const struct matrix_var *const *bp = b_;
5646 const struct matrix_var *a = *ap;
5647 const struct matrix_var *b = *bp;
5648 return strcmp (a->name, b->name);
5652 matrix_cmd_execute_display (struct display_command *cmd)
5654 const struct matrix_state *s = cmd->state;
5656 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5657 struct pivot_dimension *attr_dimension
5658 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5659 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5660 N_("Rows"), N_("Columns"));
5661 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5663 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5666 const struct matrix_var *var;
5667 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5669 vars[n_vars++] = var;
5670 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5672 struct pivot_dimension *rows = pivot_dimension_create (
5673 table, PIVOT_AXIS_ROW, N_("Variable"));
5674 for (size_t i = 0; i < n_vars; i++)
5676 const struct matrix_var *var = vars[i];
5677 pivot_category_create_leaf (
5678 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
5680 size_t r = var->value->size1;
5681 size_t c = var->value->size2;
5682 double values[] = { r, c, r * c * sizeof (double) / 1024 };
5683 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
5684 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
5687 pivot_table_submit (table);
5690 static struct matrix_cmd *
5691 matrix_parse_release (struct matrix_state *s)
5693 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5694 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
5696 size_t allocated_vars = 0;
5697 while (lex_token (s->lexer) == T_ID)
5699 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
5702 if (cmd->release.n_vars >= allocated_vars)
5703 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
5704 sizeof *cmd->release.vars);
5705 cmd->release.vars[cmd->release.n_vars++] = var;
5708 lex_error (s->lexer, _("Variable name expected."));
5711 if (!lex_match (s->lexer, T_COMMA))
5719 matrix_cmd_execute_release (struct release_command *cmd)
5721 for (size_t i = 0; i < cmd->n_vars; i++)
5723 struct matrix_var *var = cmd->vars[i];
5724 gsl_matrix_free (var->value);
5729 static struct save_file *
5730 save_file_create (struct matrix_state *s, struct file_handle *fh,
5731 struct string_array *variables,
5732 struct matrix_expr *names,
5733 struct stringi_set *strings)
5735 for (size_t i = 0; i < s->n_save_files; i++)
5737 struct save_file *sf = s->save_files[i];
5738 if (fh_equal (sf->file, fh))
5742 string_array_destroy (variables);
5743 matrix_expr_destroy (names);
5744 stringi_set_destroy (strings);
5750 struct save_file *sf = xmalloc (sizeof *sf);
5751 *sf = (struct save_file) {
5753 .dataset = s->dataset,
5754 .variables = *variables,
5756 .strings = STRINGI_SET_INITIALIZER (sf->strings),
5759 stringi_set_swap (&sf->strings, strings);
5760 stringi_set_destroy (strings);
5762 s->save_files = xrealloc (s->save_files,
5763 (s->n_save_files + 1) * sizeof *s->save_files);
5764 s->save_files[s->n_save_files++] = sf;
5768 static struct casewriter *
5769 save_file_open (struct save_file *sf, gsl_matrix *m)
5771 if (sf->writer || sf->error)
5775 size_t n_variables = caseproto_get_n_widths (
5776 casewriter_get_proto (sf->writer));
5777 if (m->size2 != n_variables)
5779 msg (ME, _("The first SAVE to %s within this matrix program "
5780 "had %zu columns, so a %zu×%zu matrix cannot be "
5782 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
5783 n_variables, m->size1, m->size2);
5791 struct dictionary *dict = dict_create (get_default_encoding ());
5793 /* Fill 'names' with user-specified names if there were any, either from
5794 sf->variables or sf->names. */
5795 struct string_array names = { .n = 0 };
5796 if (sf->variables.n)
5797 string_array_clone (&names, &sf->variables);
5800 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
5801 if (nm && is_vector (nm))
5803 gsl_vector nv = to_vector (nm);
5804 for (size_t i = 0; i < nv.size; i++)
5806 char *name = trimmed_string (gsl_vector_get (&nv, i));
5807 if (dict_id_is_valid (dict, name, true))
5808 string_array_append_nocopy (&names, name);
5813 gsl_matrix_free (nm);
5816 struct stringi_set strings;
5817 stringi_set_clone (&strings, &sf->strings);
5819 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
5824 name = names.strings[i];
5827 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
5831 int width = stringi_set_delete (&strings, name) ? 8 : 0;
5832 struct variable *var = dict_create_var (dict, name, width);
5835 msg (ME, _("Duplicate variable name %s in SAVE statement."),
5841 if (!stringi_set_is_empty (&strings))
5843 size_t count = stringi_set_count (&strings);
5844 const char *example = stringi_set_node_get_string (
5845 stringi_set_first (&strings));
5848 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
5849 "variable %s."), example);
5851 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
5852 "unknown variable: %s.",
5853 "The SAVE command STRINGS subcommand specifies %zu "
5854 "unknown variables, including %s.",
5859 stringi_set_destroy (&strings);
5860 string_array_destroy (&names);
5869 if (sf->file == fh_inline_file ())
5870 sf->writer = autopaging_writer_create (dict_get_proto (dict));
5872 sf->writer = any_writer_open (sf->file, dict);
5884 save_file_destroy (struct save_file *sf)
5888 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
5890 dataset_set_dict (sf->dataset, sf->dict);
5891 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
5895 casewriter_destroy (sf->writer);
5896 dict_unref (sf->dict);
5898 fh_unref (sf->file);
5899 string_array_destroy (&sf->variables);
5900 matrix_expr_destroy (sf->names);
5901 stringi_set_destroy (&sf->strings);
5906 static struct matrix_cmd *
5907 matrix_parse_save (struct matrix_state *s)
5909 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5910 *cmd = (struct matrix_cmd) {
5912 .save = { .expression = NULL },
5915 struct file_handle *fh = NULL;
5916 struct save_command *save = &cmd->save;
5918 struct string_array variables = STRING_ARRAY_INITIALIZER;
5919 struct matrix_expr *names = NULL;
5920 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
5922 save->expression = matrix_parse_exp (s);
5923 if (!save->expression)
5926 while (lex_match (s->lexer, T_SLASH))
5928 if (lex_match_id (s->lexer, "OUTFILE"))
5930 lex_match (s->lexer, T_EQUALS);
5933 fh = (lex_match (s->lexer, T_ASTERISK)
5935 : fh_parse (s->lexer, FH_REF_FILE, s->session));
5939 else if (lex_match_id (s->lexer, "VARIABLES"))
5941 lex_match (s->lexer, T_EQUALS);
5945 struct dictionary *d = dict_create (get_default_encoding ());
5946 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
5947 PV_NO_SCRATCH | PV_NO_DUPLICATE);
5952 string_array_clear (&variables);
5953 variables = (struct string_array) {
5959 else if (lex_match_id (s->lexer, "NAMES"))
5961 lex_match (s->lexer, T_EQUALS);
5962 matrix_expr_destroy (names);
5963 names = matrix_parse_exp (s);
5967 else if (lex_match_id (s->lexer, "STRINGS"))
5969 lex_match (s->lexer, T_EQUALS);
5970 while (lex_token (s->lexer) == T_ID)
5972 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
5974 if (!lex_match (s->lexer, T_COMMA))
5980 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
5988 if (s->prev_save_file)
5989 fh = fh_ref (s->prev_save_file);
5992 lex_sbc_missing ("OUTFILE");
5996 fh_unref (s->prev_save_file);
5997 s->prev_save_file = fh_ref (fh);
5999 if (variables.n && names)
6001 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
6002 matrix_expr_destroy (names);
6006 save->sf = save_file_create (s, fh, &variables, names, &strings);
6010 string_array_destroy (&variables);
6011 matrix_expr_destroy (names);
6012 stringi_set_destroy (&strings);
6014 matrix_cmd_destroy (cmd);
6019 matrix_cmd_execute_save (const struct save_command *save)
6021 gsl_matrix *m = matrix_expr_evaluate (save->expression);
6025 struct casewriter *writer = save_file_open (save->sf, m);
6028 gsl_matrix_free (m);
6032 const struct caseproto *proto = casewriter_get_proto (writer);
6033 for (size_t y = 0; y < m->size1; y++)
6035 struct ccase *c = case_create (proto);
6036 for (size_t x = 0; x < m->size2; x++)
6038 double d = gsl_matrix_get (m, y, x);
6039 int width = caseproto_get_width (proto, x);
6040 union value *value = case_data_rw_idx (c, x);
6044 memcpy (value->s, &d, width);
6046 casewriter_write (writer, c);
6048 gsl_matrix_free (m);
6051 static struct matrix_cmd *
6052 matrix_parse_read (struct matrix_state *s)
6054 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6055 *cmd = (struct matrix_cmd) {
6057 .read = { .format = FMT_F },
6060 struct file_handle *fh = NULL;
6061 char *encoding = NULL;
6062 struct read_command *read = &cmd->read;
6063 read->dst = matrix_lvalue_parse (s);
6068 int repetitions = 0;
6069 int record_width = 0;
6070 bool seen_format = false;
6071 while (lex_match (s->lexer, T_SLASH))
6073 if (lex_match_id (s->lexer, "FILE"))
6075 lex_match (s->lexer, T_EQUALS);
6078 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6082 else if (lex_match_id (s->lexer, "ENCODING"))
6084 lex_match (s->lexer, T_EQUALS);
6085 if (!lex_force_string (s->lexer))
6089 encoding = ss_xstrdup (lex_tokss (s->lexer));
6093 else if (lex_match_id (s->lexer, "FIELD"))
6095 lex_match (s->lexer, T_EQUALS);
6097 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6099 read->c1 = lex_integer (s->lexer);
6101 if (!lex_force_match (s->lexer, T_TO)
6102 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
6104 read->c2 = lex_integer (s->lexer) + 1;
6107 record_width = read->c2 - read->c1;
6108 if (lex_match (s->lexer, T_BY))
6110 if (!lex_force_int_range (s->lexer, "BY", 1,
6111 read->c2 - read->c1))
6113 by = lex_integer (s->lexer);
6116 if (record_width % by)
6118 msg (SE, _("BY %d does not evenly divide record width %d."),
6126 else if (lex_match_id (s->lexer, "SIZE"))
6128 lex_match (s->lexer, T_EQUALS);
6129 matrix_expr_destroy (read->size);
6130 read->size = matrix_parse_exp (s);
6134 else if (lex_match_id (s->lexer, "MODE"))
6136 lex_match (s->lexer, T_EQUALS);
6137 if (lex_match_id (s->lexer, "RECTANGULAR"))
6138 read->symmetric = false;
6139 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6140 read->symmetric = true;
6143 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6147 else if (lex_match_id (s->lexer, "REREAD"))
6148 read->reread = true;
6149 else if (lex_match_id (s->lexer, "FORMAT"))
6153 lex_sbc_only_once ("FORMAT");
6158 lex_match (s->lexer, T_EQUALS);
6160 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6163 const char *p = lex_tokcstr (s->lexer);
6164 if (c_isdigit (p[0]))
6166 repetitions = atoi (p);
6167 p += strspn (p, "0123456789");
6168 if (!fmt_from_name (p, &read->format))
6170 lex_error (s->lexer, _("Unknown format %s."), p);
6175 else if (fmt_from_name (p, &read->format))
6179 struct fmt_spec format;
6180 if (!parse_format_specifier (s->lexer, &format))
6182 read->format = format.type;
6188 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6189 "REREAD", "FORMAT");
6196 lex_sbc_missing ("FIELD");
6200 if (!read->dst->n_indexes && !read->size)
6202 msg (SE, _("SIZE is required for reading data into a full matrix "
6203 "(as opposed to a submatrix)."));
6209 if (s->prev_read_file)
6210 fh = fh_ref (s->prev_read_file);
6213 lex_sbc_missing ("FILE");
6217 fh_unref (s->prev_read_file);
6218 s->prev_read_file = fh_ref (fh);
6220 read->rf = read_file_create (s, fh);
6224 free (read->rf->encoding);
6225 read->rf->encoding = encoding;
6229 /* Field width may be specified in multiple ways:
6232 2. The format on FORMAT.
6233 3. The repetition factor on FORMAT.
6235 (2) and (3) are mutually exclusive.
6237 If more than one of these is present, they must agree. If none of them is
6238 present, then free-field format is used.
6240 if (repetitions > record_width)
6242 msg (SE, _("%d repetitions cannot fit in record width %d."),
6243 repetitions, record_width);
6246 int w = (repetitions ? record_width / repetitions
6252 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6253 "which implies field width %d, "
6254 "but BY specifies field width %d."),
6255 repetitions, record_width, w, by);
6257 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6266 matrix_cmd_destroy (cmd);
6272 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6273 struct substring data, size_t y, size_t x,
6274 int first_column, int last_column, char *error)
6276 int line_number = dfm_get_line_number (reader);
6277 struct msg_location *location = xmalloc (sizeof *location);
6278 *location = (struct msg_location) {
6279 .file_name = intern_new (dfm_get_file_name (reader)),
6280 .p[0] = { .line = line_number, .column = first_column },
6281 .p[1] = { .line = line_number, .column = last_column },
6283 struct msg *m = xmalloc (sizeof *m);
6285 .category = MSG_C_DATA,
6286 .severity = MSG_S_WARNING,
6287 .location = location,
6288 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
6289 "for matrix row %zu, column %zu: %s"),
6290 (int) data.length, data.string, fmt_name (format),
6291 y + 1, x + 1, error),
6299 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
6300 gsl_matrix *m, struct substring p, size_t y, size_t x,
6301 const char *line_start)
6303 const char *input_encoding = dfm_reader_get_encoding (reader);
6306 if (fmt_is_numeric (read->format))
6309 error = data_in (p, input_encoding, read->format,
6310 settings_get_fmt_settings (), &v, 0, NULL);
6311 if (!error && v.f == SYSMIS)
6312 error = xstrdup (_("Matrix data may not contain missing value."));
6317 uint8_t s[sizeof (double)];
6318 union value v = { .s = s };
6319 error = data_in (p, input_encoding, read->format,
6320 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6321 memcpy (&f, s, sizeof f);
6326 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6327 int c2 = c1 + ss_utf8_count_columns (p);
6328 parse_error (reader, read->format, p, y, x, c1, c2, error);
6332 gsl_matrix_set (m, y, x, f);
6333 if (read->symmetric && x != y)
6334 gsl_matrix_set (m, x, y, f);
6339 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
6340 struct substring *line, const char **startp)
6342 if (dfm_eof (reader))
6344 msg (SE, _("Unexpected end of file reading matrix data."));
6347 dfm_expand_tabs (reader);
6348 struct substring record = dfm_get_record (reader);
6349 /* XXX need to recode record into UTF-8 */
6350 *startp = record.string;
6351 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6356 matrix_read (struct read_command *read, struct dfm_reader *reader,
6359 for (size_t y = 0; y < m->size1; y++)
6361 size_t nx = read->symmetric ? y + 1 : m->size2;
6363 struct substring line = ss_empty ();
6364 const char *line_start = line.string;
6365 for (size_t x = 0; x < nx; x++)
6372 ss_ltrim (&line, ss_cstr (" ,"));
6373 if (!ss_is_empty (line))
6375 if (!matrix_read_line (read, reader, &line, &line_start))
6377 dfm_forward_record (reader);
6380 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6384 if (!matrix_read_line (read, reader, &line, &line_start))
6386 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6387 int f = x % fields_per_line;
6388 if (f == fields_per_line - 1)
6389 dfm_forward_record (reader);
6391 p = ss_substr (line, read->w * f, read->w);
6394 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6398 dfm_forward_record (reader);
6401 ss_ltrim (&line, ss_cstr (" ,"));
6402 if (!ss_is_empty (line))
6405 msg (SW, _("Trailing garbage on line \"%.*s\""),
6406 (int) line.length, line.string);
6413 matrix_cmd_execute_read (struct read_command *read)
6415 struct index_vector iv0, iv1;
6416 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6419 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6422 gsl_matrix *m = matrix_expr_evaluate (read->size);
6428 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6429 "not a %zu×%zu matrix."), m->size1, m->size2);
6430 gsl_matrix_free (m);
6436 gsl_vector v = to_vector (m);
6440 d[0] = gsl_vector_get (&v, 0);
6443 else if (v.size == 2)
6445 d[0] = gsl_vector_get (&v, 0);
6446 d[1] = gsl_vector_get (&v, 1);
6450 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6451 "not a %zu×%zu matrix."),
6452 m->size1, m->size2),
6453 gsl_matrix_free (m);
6458 gsl_matrix_free (m);
6460 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6462 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
6463 "are outside valid range."),
6474 if (read->dst->n_indexes)
6476 size_t submatrix_size[2];
6477 if (read->dst->n_indexes == 2)
6479 submatrix_size[0] = iv0.n;
6480 submatrix_size[1] = iv1.n;
6482 else if (read->dst->var->value->size1 == 1)
6484 submatrix_size[0] = 1;
6485 submatrix_size[1] = iv0.n;
6489 submatrix_size[0] = iv0.n;
6490 submatrix_size[1] = 1;
6495 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6497 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
6498 "differ from submatrix dimensions %zu×%zu."),
6500 submatrix_size[0], submatrix_size[1]);
6508 size[0] = submatrix_size[0];
6509 size[1] = submatrix_size[1];
6513 struct dfm_reader *reader = read_file_open (read->rf);
6515 dfm_reread_record (reader, 1);
6517 if (read->symmetric && size[0] != size[1])
6519 msg (SE, _("Cannot read non-square %zu×%zu matrix "
6520 "using READ with MODE=SYMMETRIC."),
6526 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6527 matrix_read (read, reader, tmp);
6528 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
6531 static struct matrix_cmd *
6532 matrix_parse_write (struct matrix_state *s)
6534 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6535 *cmd = (struct matrix_cmd) {
6539 struct file_handle *fh = NULL;
6540 char *encoding = NULL;
6541 struct write_command *write = &cmd->write;
6542 write->expression = matrix_parse_exp (s);
6543 if (!write->expression)
6547 int repetitions = 0;
6548 int record_width = 0;
6549 enum fmt_type format = FMT_F;
6550 bool has_format = false;
6551 while (lex_match (s->lexer, T_SLASH))
6553 if (lex_match_id (s->lexer, "OUTFILE"))
6555 lex_match (s->lexer, T_EQUALS);
6558 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6562 else if (lex_match_id (s->lexer, "ENCODING"))
6564 lex_match (s->lexer, T_EQUALS);
6565 if (!lex_force_string (s->lexer))
6569 encoding = ss_xstrdup (lex_tokss (s->lexer));
6573 else if (lex_match_id (s->lexer, "FIELD"))
6575 lex_match (s->lexer, T_EQUALS);
6577 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6579 write->c1 = lex_integer (s->lexer);
6581 if (!lex_force_match (s->lexer, T_TO)
6582 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
6584 write->c2 = lex_integer (s->lexer) + 1;
6587 record_width = write->c2 - write->c1;
6588 if (lex_match (s->lexer, T_BY))
6590 if (!lex_force_int_range (s->lexer, "BY", 1,
6591 write->c2 - write->c1))
6593 by = lex_integer (s->lexer);
6596 if (record_width % by)
6598 msg (SE, _("BY %d does not evenly divide record width %d."),
6606 else if (lex_match_id (s->lexer, "MODE"))
6608 lex_match (s->lexer, T_EQUALS);
6609 if (lex_match_id (s->lexer, "RECTANGULAR"))
6610 write->triangular = false;
6611 else if (lex_match_id (s->lexer, "TRIANGULAR"))
6612 write->triangular = true;
6615 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
6619 else if (lex_match_id (s->lexer, "HOLD"))
6621 else if (lex_match_id (s->lexer, "FORMAT"))
6623 if (has_format || write->format)
6625 lex_sbc_only_once ("FORMAT");
6629 lex_match (s->lexer, T_EQUALS);
6631 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6634 const char *p = lex_tokcstr (s->lexer);
6635 if (c_isdigit (p[0]))
6637 repetitions = atoi (p);
6638 p += strspn (p, "0123456789");
6639 if (!fmt_from_name (p, &format))
6641 lex_error (s->lexer, _("Unknown format %s."), p);
6647 else if (fmt_from_name (p, &format))
6654 struct fmt_spec spec;
6655 if (!parse_format_specifier (s->lexer, &spec))
6657 write->format = xmemdup (&spec, sizeof spec);
6662 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
6670 lex_sbc_missing ("FIELD");
6676 if (s->prev_write_file)
6677 fh = fh_ref (s->prev_write_file);
6680 lex_sbc_missing ("OUTFILE");
6684 fh_unref (s->prev_write_file);
6685 s->prev_write_file = fh_ref (fh);
6687 write->wf = write_file_create (s, fh);
6691 free (write->wf->encoding);
6692 write->wf->encoding = encoding;
6696 /* Field width may be specified in multiple ways:
6699 2. The format on FORMAT.
6700 3. The repetition factor on FORMAT.
6702 (2) and (3) are mutually exclusive.
6704 If more than one of these is present, they must agree. If none of them is
6705 present, then free-field format is used.
6707 if (repetitions > record_width)
6709 msg (SE, _("%d repetitions cannot fit in record width %d."),
6710 repetitions, record_width);
6713 int w = (repetitions ? record_width / repetitions
6714 : write->format ? write->format->w
6719 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6720 "which implies field width %d, "
6721 "but BY specifies field width %d."),
6722 repetitions, record_width, w, by);
6724 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6728 if (w && !write->format)
6730 write->format = xmalloc (sizeof *write->format);
6731 *write->format = (struct fmt_spec) { .type = format, .w = w };
6733 if (!fmt_check_output (write->format))
6737 if (write->format && fmt_var_width (write->format) > sizeof (double))
6739 char s[FMT_STRING_LEN_MAX + 1];
6740 fmt_to_string (write->format, s);
6741 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
6742 s, sizeof (double));
6750 matrix_cmd_destroy (cmd);
6755 matrix_cmd_execute_write (struct write_command *write)
6757 gsl_matrix *m = matrix_expr_evaluate (write->expression);
6761 if (write->triangular && m->size1 != m->size2)
6763 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
6764 "the matrix to be written has dimensions %zu×%zu."),
6765 m->size1, m->size2);
6766 gsl_matrix_free (m);
6770 struct dfm_writer *writer = write_file_open (write->wf);
6771 if (!writer || !m->size1)
6773 gsl_matrix_free (m);
6777 const struct fmt_settings *settings = settings_get_fmt_settings ();
6778 struct u8_line *line = write->wf->held;
6779 for (size_t y = 0; y < m->size1; y++)
6783 line = xmalloc (sizeof *line);
6784 u8_line_init (line);
6786 size_t nx = write->triangular ? y + 1 : m->size2;
6788 for (size_t x = 0; x < nx; x++)
6791 double f = gsl_matrix_get (m, y, x);
6795 if (fmt_is_numeric (write->format->type))
6798 v.s = (uint8_t *) &f;
6799 s = data_out (&v, NULL, write->format, settings);
6803 s = xmalloc (DBL_BUFSIZE_BOUND);
6804 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
6805 >= DBL_BUFSIZE_BOUND)
6808 size_t len = strlen (s);
6809 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
6810 if (width + x0 > write->c2)
6812 dfm_put_record_utf8 (writer, line->s.ss.string,
6814 u8_line_clear (line);
6817 u8_line_put (line, x0, x0 + width, s, len);
6820 x0 += write->format ? write->format->w : width + 1;
6823 if (y + 1 >= m->size1 && write->hold)
6825 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
6826 u8_line_clear (line);
6830 u8_line_destroy (line);
6834 write->wf->held = line;
6836 gsl_matrix_free (m);
6839 static struct matrix_cmd *
6840 matrix_parse_get (struct matrix_state *s)
6842 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6843 *cmd = (struct matrix_cmd) {
6846 .dataset = s->dataset,
6847 .user = { .treatment = MGET_ERROR },
6848 .system = { .treatment = MGET_ERROR },
6852 struct get_command *get = &cmd->get;
6853 get->dst = matrix_lvalue_parse (s);
6857 while (lex_match (s->lexer, T_SLASH))
6859 if (lex_match_id (s->lexer, "FILE"))
6861 lex_match (s->lexer, T_EQUALS);
6863 fh_unref (get->file);
6864 if (lex_match (s->lexer, T_ASTERISK))
6868 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6873 else if (lex_match_id (s->lexer, "ENCODING"))
6875 lex_match (s->lexer, T_EQUALS);
6876 if (!lex_force_string (s->lexer))
6879 free (get->encoding);
6880 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
6884 else if (lex_match_id (s->lexer, "VARIABLES"))
6886 lex_match (s->lexer, T_EQUALS);
6890 lex_sbc_only_once ("VARIABLES");
6894 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
6897 else if (lex_match_id (s->lexer, "NAMES"))
6899 lex_match (s->lexer, T_EQUALS);
6900 if (!lex_force_id (s->lexer))
6903 struct substring name = lex_tokss (s->lexer);
6904 get->names = matrix_var_lookup (s, name);
6906 get->names = matrix_var_create (s, name);
6909 else if (lex_match_id (s->lexer, "MISSING"))
6911 lex_match (s->lexer, T_EQUALS);
6912 if (lex_match_id (s->lexer, "ACCEPT"))
6913 get->user.treatment = MGET_ACCEPT;
6914 else if (lex_match_id (s->lexer, "OMIT"))
6915 get->user.treatment = MGET_OMIT;
6916 else if (lex_is_number (s->lexer))
6918 get->user.treatment = MGET_RECODE;
6919 get->user.substitute = lex_number (s->lexer);
6924 lex_error (s->lexer, NULL);
6928 else if (lex_match_id (s->lexer, "SYSMIS"))
6930 lex_match (s->lexer, T_EQUALS);
6931 if (lex_match_id (s->lexer, "OMIT"))
6932 get->system.treatment = MGET_OMIT;
6933 else if (lex_is_number (s->lexer))
6935 get->system.treatment = MGET_RECODE;
6936 get->system.substitute = lex_number (s->lexer);
6941 lex_error (s->lexer, NULL);
6947 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
6948 "MISSING", "SYSMIS");
6953 if (get->user.treatment != MGET_ACCEPT)
6954 get->system.treatment = MGET_ERROR;
6959 matrix_cmd_destroy (cmd);
6964 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
6965 const struct dictionary *dict)
6967 struct variable **vars;
6972 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
6973 &vars, &n_vars, PV_NUMERIC))
6978 n_vars = dict_get_var_cnt (dict);
6979 vars = xnmalloc (n_vars, sizeof *vars);
6980 for (size_t i = 0; i < n_vars; i++)
6982 struct variable *var = dict_get_var (dict, i);
6983 if (!var_is_numeric (var))
6985 msg (SE, _("GET: Variable %s is not numeric."),
6986 var_get_name (var));
6996 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
6997 for (size_t i = 0; i < n_vars; i++)
6999 char s[sizeof (double)];
7001 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
7002 memcpy (&f, s, sizeof f);
7003 gsl_matrix_set (names, i, 0, f);
7006 gsl_matrix_free (get->names->value);
7007 get->names->value = names;
7011 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
7012 long long int casenum = 1;
7014 for (struct ccase *c = casereader_read (reader); c;
7015 c = casereader_read (reader), casenum++)
7017 if (n_rows >= m->size1)
7019 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
7020 for (size_t y = 0; y < n_rows; y++)
7021 for (size_t x = 0; x < n_vars; x++)
7022 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
7023 gsl_matrix_free (m);
7028 for (size_t x = 0; x < n_vars; x++)
7030 const struct variable *var = vars[x];
7031 double d = case_num (c, var);
7034 if (get->system.treatment == MGET_RECODE)
7035 d = get->system.substitute;
7036 else if (get->system.treatment == MGET_OMIT)
7040 msg (SE, _("GET: Variable %s in case %lld "
7041 "is system-missing."),
7042 var_get_name (var), casenum);
7046 else if (var_is_num_missing (var, d, MV_USER))
7048 if (get->user.treatment == MGET_RECODE)
7049 d = get->user.substitute;
7050 else if (get->user.treatment == MGET_OMIT)
7052 else if (get->user.treatment != MGET_ACCEPT)
7054 msg (SE, _("GET: Variable %s in case %lld has user-missing "
7056 var_get_name (var), casenum, d);
7060 gsl_matrix_set (m, n_rows, x, d);
7071 matrix_lvalue_evaluate_and_assign (get->dst, m);
7074 gsl_matrix_free (m);
7079 matrix_open_casereader (const char *command_name,
7080 struct file_handle *file, const char *encoding,
7081 struct dataset *dataset,
7082 struct casereader **readerp, struct dictionary **dictp)
7086 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
7087 return *readerp != NULL;
7091 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
7093 msg (ME, _("The %s command cannot read an empty active file."),
7097 *readerp = proc_open (dataset);
7098 *dictp = dict_ref (dataset_dict (dataset));
7104 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
7105 struct casereader *reader, struct dictionary *dict)
7108 casereader_destroy (reader);
7110 proc_commit (dataset);
7114 matrix_cmd_execute_get (struct get_command *get)
7116 struct casereader *r;
7117 struct dictionary *d;
7118 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
7121 matrix_cmd_execute_get__ (get, r, d);
7122 matrix_close_casereader (get->file, get->dataset, r, d);
7127 match_rowtype (struct lexer *lexer)
7129 static const char *rowtypes[] = {
7130 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7132 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7134 for (size_t i = 0; i < n_rowtypes; i++)
7135 if (lex_match_id (lexer, rowtypes[i]))
7138 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7143 parse_var_names (struct lexer *lexer, struct string_array *sa)
7145 lex_match (lexer, T_EQUALS);
7147 string_array_clear (sa);
7149 struct dictionary *dict = dict_create (get_default_encoding ());
7152 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7153 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7158 for (size_t i = 0; i < n_names; i++)
7159 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7160 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7162 msg (SE, _("Variable name %s is reserved."), names[i]);
7163 for (size_t j = 0; j < n_names; j++)
7169 string_array_clear (sa);
7170 sa->strings = names;
7171 sa->n = sa->allocated = n_names;
7177 msave_common_destroy (struct msave_common *common)
7181 fh_unref (common->outfile);
7182 string_array_destroy (&common->variables);
7183 string_array_destroy (&common->fnames);
7184 string_array_destroy (&common->snames);
7186 for (size_t i = 0; i < common->n_factors; i++)
7187 matrix_expr_destroy (common->factors[i]);
7188 free (common->factors);
7190 for (size_t i = 0; i < common->n_splits; i++)
7191 matrix_expr_destroy (common->splits[i]);
7192 free (common->splits);
7194 dict_unref (common->dict);
7195 casewriter_destroy (common->writer);
7202 compare_variables (const char *keyword,
7203 const struct string_array *new,
7204 const struct string_array *old)
7210 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7211 "on the first MSAVE within MATRIX."), keyword);
7214 else if (!string_array_equal_case (old, new))
7216 msg (SE, _("%s must specify the same variables each time within "
7217 "a given MATRIX."), keyword);
7224 static struct matrix_cmd *
7225 matrix_parse_msave (struct matrix_state *s)
7227 struct msave_common *common = xmalloc (sizeof *common);
7228 *common = (struct msave_common) { .outfile = NULL };
7230 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7231 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7233 struct matrix_expr *splits = NULL;
7234 struct matrix_expr *factors = NULL;
7236 struct msave_command *msave = &cmd->msave;
7237 msave->expr = matrix_parse_exp (s);
7241 while (lex_match (s->lexer, T_SLASH))
7243 if (lex_match_id (s->lexer, "TYPE"))
7245 lex_match (s->lexer, T_EQUALS);
7247 msave->rowtype = match_rowtype (s->lexer);
7248 if (!msave->rowtype)
7251 else if (lex_match_id (s->lexer, "OUTFILE"))
7253 lex_match (s->lexer, T_EQUALS);
7255 fh_unref (common->outfile);
7256 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7257 if (!common->outfile)
7260 else if (lex_match_id (s->lexer, "VARIABLES"))
7262 if (!parse_var_names (s->lexer, &common->variables))
7265 else if (lex_match_id (s->lexer, "FNAMES"))
7267 if (!parse_var_names (s->lexer, &common->fnames))
7270 else if (lex_match_id (s->lexer, "SNAMES"))
7272 if (!parse_var_names (s->lexer, &common->snames))
7275 else if (lex_match_id (s->lexer, "SPLIT"))
7277 lex_match (s->lexer, T_EQUALS);
7279 matrix_expr_destroy (splits);
7280 splits = matrix_parse_exp (s);
7284 else if (lex_match_id (s->lexer, "FACTOR"))
7286 lex_match (s->lexer, T_EQUALS);
7288 matrix_expr_destroy (factors);
7289 factors = matrix_parse_exp (s);
7295 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7296 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7300 if (!msave->rowtype)
7302 lex_sbc_missing ("TYPE");
7308 if (common->fnames.n && !factors)
7310 msg (SE, _("FNAMES requires FACTOR."));
7313 if (common->snames.n && !splits)
7315 msg (SE, _("SNAMES requires SPLIT."));
7318 if (!common->outfile)
7320 lex_sbc_missing ("OUTFILE");
7327 if (common->outfile)
7329 if (!fh_equal (common->outfile, s->common->outfile))
7331 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7332 "within a single MATRIX command."));
7336 if (!compare_variables ("VARIABLES",
7337 &common->variables, &s->common->variables)
7338 || !compare_variables ("FNAMES", &common->fnames, &s->common->fnames)
7339 || !compare_variables ("SNAMES", &common->snames, &s->common->snames))
7341 msave_common_destroy (common);
7343 msave->common = s->common;
7345 struct msave_common *c = s->common;
7348 if (c->n_factors >= c->allocated_factors)
7349 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7350 sizeof *c->factors);
7351 c->factors[c->n_factors++] = factors;
7353 if (c->n_factors > 0)
7354 msave->factors = c->factors[c->n_factors - 1];
7358 if (c->n_splits >= c->allocated_splits)
7359 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7361 c->splits[c->n_splits++] = splits;
7363 if (c->n_splits > 0)
7364 msave->splits = c->splits[c->n_splits - 1];
7369 matrix_expr_destroy (splits);
7370 matrix_expr_destroy (factors);
7371 msave_common_destroy (common);
7372 matrix_cmd_destroy (cmd);
7377 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7379 gsl_matrix *m = matrix_expr_evaluate (e);
7385 msg (SE, _("%s expression must evaluate to vector, "
7386 "not a %zu×%zu matrix."),
7387 name, m->size1, m->size2);
7388 gsl_matrix_free (m);
7392 return matrix_to_vector (m);
7396 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7398 for (size_t i = 0; i < vars->n; i++)
7399 if (!dict_create_var (d, vars->strings[i], 0))
7400 return vars->strings[i];
7404 static struct dictionary *
7405 msave_create_dict (const struct msave_common *common)
7407 struct dictionary *dict = dict_create (get_default_encoding ());
7409 const char *dup_split = msave_add_vars (dict, &common->snames);
7412 /* Should not be possible because the parser ensures that the names are
7417 dict_create_var_assert (dict, "ROWTYPE_", 8);
7419 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7422 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
7426 dict_create_var_assert (dict, "VARNAME_", 8);
7428 const char *dup_var = msave_add_vars (dict, &common->variables);
7431 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
7443 matrix_cmd_execute_msave (struct msave_command *msave)
7445 struct msave_common *common = msave->common;
7446 gsl_matrix *m = NULL;
7447 gsl_vector *factors = NULL;
7448 gsl_vector *splits = NULL;
7450 m = matrix_expr_evaluate (msave->expr);
7454 if (!common->variables.n)
7455 for (size_t i = 0; i < m->size2; i++)
7456 string_array_append_nocopy (&common->variables,
7457 xasprintf ("COL%zu", i + 1));
7458 else if (m->size2 != common->variables.n)
7461 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7462 m->size2, common->variables.n);
7468 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7472 if (!common->fnames.n)
7473 for (size_t i = 0; i < factors->size; i++)
7474 string_array_append_nocopy (&common->fnames,
7475 xasprintf ("FAC%zu", i + 1));
7476 else if (factors->size != common->fnames.n)
7478 msg (SE, _("There are %zu factor variables, "
7479 "but %zu factor values were supplied."),
7480 common->fnames.n, factors->size);
7487 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7491 if (!common->snames.n)
7492 for (size_t i = 0; i < splits->size; i++)
7493 string_array_append_nocopy (&common->snames,
7494 xasprintf ("SPL%zu", i + 1));
7495 else if (splits->size != common->snames.n)
7497 msg (SE, _("There are %zu split variables, "
7498 "but %zu split values were supplied."),
7499 common->snames.n, splits->size);
7504 if (!common->writer)
7506 struct dictionary *dict = msave_create_dict (common);
7510 common->writer = any_writer_open (common->outfile, dict);
7511 if (!common->writer)
7517 common->dict = dict;
7520 bool matrix = (!strcmp (msave->rowtype, "COV")
7521 || !strcmp (msave->rowtype, "CORR"));
7522 for (size_t y = 0; y < m->size1; y++)
7524 struct ccase *c = case_create (dict_get_proto (common->dict));
7527 /* Split variables */
7529 for (size_t i = 0; i < splits->size; i++)
7530 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
7533 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7534 msave->rowtype, ' ');
7538 for (size_t i = 0; i < factors->size; i++)
7539 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
7542 const char *varname_ = (matrix && y < common->variables.n
7543 ? common->variables.strings[y]
7545 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7548 /* Continuous variables. */
7549 for (size_t x = 0; x < m->size2; x++)
7550 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
7551 casewriter_write (common->writer, c);
7555 gsl_matrix_free (m);
7556 gsl_vector_free (factors);
7557 gsl_vector_free (splits);
7560 static struct matrix_cmd *
7561 matrix_parse_mget (struct matrix_state *s)
7563 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7564 *cmd = (struct matrix_cmd) {
7568 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
7572 struct mget_command *mget = &cmd->mget;
7574 lex_match (s->lexer, T_SLASH);
7575 while (lex_token (s->lexer) != T_ENDCMD)
7577 if (lex_match_id (s->lexer, "FILE"))
7579 lex_match (s->lexer, T_EQUALS);
7581 fh_unref (mget->file);
7582 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7586 else if (lex_match_id (s->lexer, "ENCODING"))
7588 lex_match (s->lexer, T_EQUALS);
7589 if (!lex_force_string (s->lexer))
7592 free (mget->encoding);
7593 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
7597 else if (lex_match_id (s->lexer, "TYPE"))
7599 lex_match (s->lexer, T_EQUALS);
7600 while (lex_token (s->lexer) != T_SLASH
7601 && lex_token (s->lexer) != T_ENDCMD)
7603 const char *rowtype = match_rowtype (s->lexer);
7607 stringi_set_insert (&mget->rowtypes, rowtype);
7612 lex_error_expecting (s->lexer, "FILE", "TYPE");
7615 lex_match (s->lexer, T_SLASH);
7620 matrix_cmd_destroy (cmd);
7624 static const struct variable *
7625 get_a8_var (const struct dictionary *d, const char *name)
7627 const struct variable *v = dict_lookup_var (d, name);
7630 msg (SE, _("Matrix data file lacks %s variable."), name);
7633 if (var_get_width (v) != 8)
7635 msg (SE, _("%s variable in matrix data file must be "
7636 "8-byte string, but it has width %d."),
7637 name, var_get_width (v));
7644 var_changed (const struct ccase *ca, const struct ccase *cb,
7645 const struct variable *var)
7648 ? !value_equal (case_data (ca, var), case_data (cb, var),
7649 var_get_width (var))
7654 vars_changed (const struct ccase *ca, const struct ccase *cb,
7655 const struct dictionary *d,
7656 size_t first_var, size_t n_vars)
7658 for (size_t i = 0; i < n_vars; i++)
7660 const struct variable *v = dict_get_var (d, first_var + i);
7661 if (var_changed (ca, cb, v))
7668 vars_all_missing (const struct ccase *c, const struct dictionary *d,
7669 size_t first_var, size_t n_vars)
7671 for (size_t i = 0; i < n_vars; i++)
7672 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
7678 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
7679 const struct dictionary *d,
7680 const struct variable *rowtype_var,
7681 const struct stringi_set *accepted_rowtypes,
7682 struct matrix_state *s,
7683 size_t ss, size_t sn, size_t si,
7684 size_t fs, size_t fn, size_t fi,
7685 size_t cs, size_t cn,
7686 struct pivot_table *pt,
7687 struct pivot_dimension *var_dimension)
7692 /* Is this a matrix for pooled data, either where there are no factor
7693 variables or the factor variables are missing? */
7694 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
7696 struct substring rowtype = case_ss (rows[0], rowtype_var);
7697 ss_rtrim (&rowtype, ss_cstr (" "));
7698 if (!stringi_set_is_empty (accepted_rowtypes)
7699 && !stringi_set_contains_len (accepted_rowtypes,
7700 rowtype.string, rowtype.length))
7703 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
7704 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
7705 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
7706 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
7707 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
7708 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
7712 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
7713 (int) rowtype.length, rowtype.string);
7717 struct string name = DS_EMPTY_INITIALIZER;
7718 ds_put_cstr (&name, prefix);
7720 ds_put_format (&name, "F%zu", fi);
7722 ds_put_format (&name, "S%zu", si);
7724 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
7726 mv = matrix_var_create (s, ds_ss (&name));
7729 msg (SW, _("Matrix data file contains variable with existing name %s."),
7731 goto exit_free_name;
7734 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
7735 size_t n_missing = 0;
7736 for (size_t y = 0; y < n_rows; y++)
7738 for (size_t x = 0; x < cn; x++)
7740 struct variable *var = dict_get_var (d, cs + x);
7741 double value = case_num (rows[y], var);
7742 if (var_is_num_missing (var, value, MV_ANY))
7747 gsl_matrix_set (m, y, x, value);
7751 int var_index = pivot_category_create_leaf (
7752 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
7753 double values[] = { n_rows, cn };
7754 for (size_t j = 0; j < sn; j++)
7756 struct variable *var = dict_get_var (d, ss + j);
7757 const union value *value = case_data (rows[0], var);
7758 pivot_table_put2 (pt, j, var_index,
7759 pivot_value_new_var_value (var, value));
7761 for (size_t j = 0; j < fn; j++)
7763 struct variable *var = dict_get_var (d, fs + j);
7764 const union value sysmis = { .f = SYSMIS };
7765 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
7766 pivot_table_put2 (pt, j + sn, var_index,
7767 pivot_value_new_var_value (var, value));
7769 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
7770 pivot_table_put2 (pt, j + sn + fn, var_index,
7771 pivot_value_new_integer (values[j]));
7774 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
7775 "value, which was treated as zero.",
7776 "Matrix data file variable %s contains %zu missing "
7777 "values, which were treated as zero.", n_missing),
7778 ds_cstr (&name), n_missing);
7785 for (size_t y = 0; y < n_rows; y++)
7786 case_unref (rows[y]);
7790 matrix_cmd_execute_mget__ (struct mget_command *mget,
7791 struct casereader *r, const struct dictionary *d)
7793 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
7794 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
7795 if (!rowtype_ || !varname_)
7798 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
7800 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
7803 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
7805 msg (SE, _("Matrix data file contains no continuous variables."));
7809 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
7811 const struct variable *v = dict_get_var (d, i);
7812 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
7815 _("Matrix data file contains unexpected string variable %s."),
7821 /* SPLIT variables. */
7823 size_t sn = var_get_dict_index (rowtype_);
7824 struct ccase *sc = NULL;
7827 /* FACTOR variables. */
7828 size_t fs = var_get_dict_index (rowtype_) + 1;
7829 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
7830 struct ccase *fc = NULL;
7833 /* Continuous variables. */
7834 size_t cs = var_get_dict_index (varname_) + 1;
7835 size_t cn = dict_get_var_cnt (d) - cs;
7836 struct ccase *cc = NULL;
7839 struct pivot_table *pt = pivot_table_create (
7840 N_("Matrix Variables Created by MGET"));
7841 struct pivot_dimension *attr_dimension = pivot_dimension_create (
7842 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
7843 struct pivot_dimension *var_dimension = pivot_dimension_create (
7844 pt, PIVOT_AXIS_ROW, N_("Variable"));
7847 struct pivot_category *splits = pivot_category_create_group (
7848 attr_dimension->root, N_("Split Values"));
7849 for (size_t i = 0; i < sn; i++)
7850 pivot_category_create_leaf (splits, pivot_value_new_variable (
7851 dict_get_var (d, ss + i)));
7855 struct pivot_category *factors = pivot_category_create_group (
7856 attr_dimension->root, N_("Factors"));
7857 for (size_t i = 0; i < fn; i++)
7858 pivot_category_create_leaf (factors, pivot_value_new_variable (
7859 dict_get_var (d, fs + i)));
7861 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
7862 N_("Rows"), N_("Columns"));
7865 struct ccase **rows = NULL;
7866 size_t allocated_rows = 0;
7870 while ((c = casereader_read (r)) != NULL)
7872 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
7882 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
7883 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
7884 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
7887 if (change != NOTHING_CHANGED)
7889 matrix_mget_commit_var (rows, n_rows, d,
7890 rowtype_, &mget->rowtypes,
7901 if (n_rows >= allocated_rows)
7902 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
7905 if (change == SPLITS_CHANGED)
7911 /* Reset the factor number, if there are factors. */
7915 if (row_has_factors)
7921 else if (change == FACTORS_CHANGED)
7923 if (row_has_factors)
7929 matrix_mget_commit_var (rows, n_rows, d,
7930 rowtype_, &mget->rowtypes,
7942 if (var_dimension->n_leaves)
7943 pivot_table_submit (pt);
7945 pivot_table_unref (pt);
7949 matrix_cmd_execute_mget (struct mget_command *mget)
7951 struct casereader *r;
7952 struct dictionary *d;
7953 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
7954 mget->state->dataset, &r, &d))
7956 matrix_cmd_execute_mget__ (mget, r, d);
7957 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
7962 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
7964 if (!lex_force_id (s->lexer))
7967 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
7969 *varp = matrix_var_create (s, lex_tokss (s->lexer));
7974 static struct matrix_cmd *
7975 matrix_parse_eigen (struct matrix_state *s)
7977 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7978 *cmd = (struct matrix_cmd) {
7980 .eigen = { .expr = NULL }
7983 struct eigen_command *eigen = &cmd->eigen;
7984 if (!lex_force_match (s->lexer, T_LPAREN))
7986 eigen->expr = matrix_parse_expr (s);
7988 || !lex_force_match (s->lexer, T_COMMA)
7989 || !matrix_parse_dst_var (s, &eigen->evec)
7990 || !lex_force_match (s->lexer, T_COMMA)
7991 || !matrix_parse_dst_var (s, &eigen->eval)
7992 || !lex_force_match (s->lexer, T_RPAREN))
7998 matrix_cmd_destroy (cmd);
8003 matrix_cmd_execute_eigen (struct eigen_command *eigen)
8005 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
8008 if (!is_symmetric (A))
8010 msg (SE, _("Argument of EIGEN must be symmetric."));
8011 gsl_matrix_free (A);
8015 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
8016 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
8017 gsl_vector v_eval = to_vector (eval);
8018 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
8019 gsl_eigen_symmv (A, &v_eval, evec, w);
8020 gsl_eigen_symmv_free (w);
8022 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
8024 gsl_matrix_free (A);
8026 gsl_matrix_free (eigen->eval->value);
8027 eigen->eval->value = eval;
8029 gsl_matrix_free (eigen->evec->value);
8030 eigen->evec->value = evec;
8033 static struct matrix_cmd *
8034 matrix_parse_setdiag (struct matrix_state *s)
8036 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
8037 *cmd = (struct matrix_cmd) {
8038 .type = MCMD_SETDIAG,
8039 .setdiag = { .dst = NULL }
8042 struct setdiag_command *setdiag = &cmd->setdiag;
8043 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
8046 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
8049 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
8054 if (!lex_force_match (s->lexer, T_COMMA))
8057 setdiag->expr = matrix_parse_expr (s);
8061 if (!lex_force_match (s->lexer, T_RPAREN))
8067 matrix_cmd_destroy (cmd);
8072 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
8074 gsl_matrix *dst = setdiag->dst->value;
8077 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
8078 setdiag->dst->name);
8082 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
8086 size_t n = MIN (dst->size1, dst->size2);
8087 if (is_scalar (src))
8089 double d = to_scalar (src);
8090 for (size_t i = 0; i < n; i++)
8091 gsl_matrix_set (dst, i, i, d);
8093 else if (is_vector (src))
8095 gsl_vector v = to_vector (src);
8096 for (size_t i = 0; i < n && i < v.size; i++)
8097 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
8100 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
8101 "not a %zu×%zu matrix."),
8102 src->size1, src->size2);
8103 gsl_matrix_free (src);
8106 static struct matrix_cmd *
8107 matrix_parse_svd (struct matrix_state *s)
8109 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
8110 *cmd = (struct matrix_cmd) {
8112 .svd = { .expr = NULL }
8115 struct svd_command *svd = &cmd->svd;
8116 if (!lex_force_match (s->lexer, T_LPAREN))
8118 svd->expr = matrix_parse_expr (s);
8120 || !lex_force_match (s->lexer, T_COMMA)
8121 || !matrix_parse_dst_var (s, &svd->u)
8122 || !lex_force_match (s->lexer, T_COMMA)
8123 || !matrix_parse_dst_var (s, &svd->s)
8124 || !lex_force_match (s->lexer, T_COMMA)
8125 || !matrix_parse_dst_var (s, &svd->v)
8126 || !lex_force_match (s->lexer, T_RPAREN))
8132 matrix_cmd_destroy (cmd);
8137 matrix_cmd_execute_svd (struct svd_command *svd)
8139 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8143 if (m->size1 >= m->size2)
8146 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8147 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8148 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8149 gsl_vector *work = gsl_vector_alloc (A->size2);
8150 gsl_linalg_SV_decomp (A, V, &Sv, work);
8151 gsl_vector_free (work);
8153 matrix_var_set (svd->u, A);
8154 matrix_var_set (svd->s, S);
8155 matrix_var_set (svd->v, V);
8159 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8160 gsl_matrix_transpose_memcpy (At, m);
8161 gsl_matrix_free (m);
8163 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8164 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8165 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8166 gsl_vector *work = gsl_vector_alloc (At->size2);
8167 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8168 gsl_vector_free (work);
8170 matrix_var_set (svd->v, At);
8171 matrix_var_set (svd->s, St);
8172 matrix_var_set (svd->u, Vt);
8177 matrix_cmd_execute (struct matrix_cmd *cmd)
8182 matrix_cmd_execute_compute (&cmd->compute);
8186 matrix_cmd_execute_print (&cmd->print);
8190 return matrix_cmd_execute_do_if (&cmd->do_if);
8193 matrix_cmd_execute_loop (&cmd->loop);
8200 matrix_cmd_execute_display (&cmd->display);
8204 matrix_cmd_execute_release (&cmd->release);
8208 matrix_cmd_execute_save (&cmd->save);
8212 matrix_cmd_execute_read (&cmd->read);
8216 matrix_cmd_execute_write (&cmd->write);
8220 matrix_cmd_execute_get (&cmd->get);
8224 matrix_cmd_execute_msave (&cmd->msave);
8228 matrix_cmd_execute_mget (&cmd->mget);
8232 matrix_cmd_execute_eigen (&cmd->eigen);
8236 matrix_cmd_execute_setdiag (&cmd->setdiag);
8240 matrix_cmd_execute_svd (&cmd->svd);
8248 matrix_cmds_uninit (struct matrix_cmds *cmds)
8250 for (size_t i = 0; i < cmds->n; i++)
8251 matrix_cmd_destroy (cmds->commands[i]);
8252 free (cmds->commands);
8256 matrix_cmd_destroy (struct matrix_cmd *cmd)
8264 matrix_lvalue_destroy (cmd->compute.lvalue);
8265 matrix_expr_destroy (cmd->compute.rvalue);
8269 matrix_expr_destroy (cmd->print.expression);
8270 free (cmd->print.title);
8271 print_labels_destroy (cmd->print.rlabels);
8272 print_labels_destroy (cmd->print.clabels);
8276 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8278 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8279 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
8281 free (cmd->do_if.clauses);
8285 matrix_expr_destroy (cmd->loop.start);
8286 matrix_expr_destroy (cmd->loop.end);
8287 matrix_expr_destroy (cmd->loop.increment);
8288 matrix_expr_destroy (cmd->loop.top_condition);
8289 matrix_expr_destroy (cmd->loop.bottom_condition);
8290 matrix_cmds_uninit (&cmd->loop.commands);
8300 free (cmd->release.vars);
8304 matrix_expr_destroy (cmd->save.expression);
8308 matrix_lvalue_destroy (cmd->read.dst);
8309 matrix_expr_destroy (cmd->read.size);
8313 matrix_expr_destroy (cmd->write.expression);
8314 free (cmd->write.format);
8318 matrix_lvalue_destroy (cmd->get.dst);
8319 fh_unref (cmd->get.file);
8320 free (cmd->get.encoding);
8321 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8325 matrix_expr_destroy (cmd->msave.expr);
8329 fh_unref (cmd->mget.file);
8330 stringi_set_destroy (&cmd->mget.rowtypes);
8334 matrix_expr_destroy (cmd->eigen.expr);
8338 matrix_expr_destroy (cmd->setdiag.expr);
8342 matrix_expr_destroy (cmd->svd.expr);
8348 struct matrix_command_name
8351 struct matrix_cmd *(*parse) (struct matrix_state *);
8354 static const struct matrix_command_name *
8355 matrix_parse_command_name (struct lexer *lexer)
8357 static const struct matrix_command_name commands[] = {
8358 { "COMPUTE", matrix_parse_compute },
8359 { "CALL EIGEN", matrix_parse_eigen },
8360 { "CALL SETDIAG", matrix_parse_setdiag },
8361 { "CALL SVD", matrix_parse_svd },
8362 { "PRINT", matrix_parse_print },
8363 { "DO IF", matrix_parse_do_if },
8364 { "LOOP", matrix_parse_loop },
8365 { "BREAK", matrix_parse_break },
8366 { "READ", matrix_parse_read },
8367 { "WRITE", matrix_parse_write },
8368 { "GET", matrix_parse_get },
8369 { "SAVE", matrix_parse_save },
8370 { "MGET", matrix_parse_mget },
8371 { "MSAVE", matrix_parse_msave },
8372 { "DISPLAY", matrix_parse_display },
8373 { "RELEASE", matrix_parse_release },
8375 static size_t n = sizeof commands / sizeof *commands;
8377 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8378 if (lex_match_phrase (lexer, c->name))
8383 static struct matrix_cmd *
8384 matrix_parse_command (struct matrix_state *s)
8386 size_t nesting_level = SIZE_MAX;
8388 struct matrix_cmd *c = NULL;
8389 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
8391 lex_error (s->lexer, _("Unknown matrix command."));
8392 else if (!cmd->parse)
8393 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8397 nesting_level = output_open_group (
8398 group_item_create_nocopy (utf8_to_title (cmd->name),
8399 utf8_to_title (cmd->name)));
8404 lex_end_of_command (s->lexer);
8405 lex_discard_rest_of_command (s->lexer);
8406 if (nesting_level != SIZE_MAX)
8407 output_close_groups (nesting_level);
8413 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8415 if (!lex_force_match (lexer, T_ENDCMD))
8418 struct matrix_state state = {
8420 .session = dataset_session (ds),
8422 .vars = HMAP_INITIALIZER (state.vars),
8427 while (lex_match (lexer, T_ENDCMD))
8429 if (lex_token (lexer) == T_STOP)
8431 msg (SE, _("Unexpected end of input expecting matrix command."));
8435 if (lex_match_phrase (lexer, "END MATRIX"))
8438 struct matrix_cmd *c = matrix_parse_command (&state);
8441 matrix_cmd_execute (c);
8442 matrix_cmd_destroy (c);
8446 struct matrix_var *var, *next;
8447 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8450 gsl_matrix_free (var->value);
8451 hmap_delete (&state.vars, &var->hmap_node);
8454 hmap_destroy (&state.vars);
8455 msave_common_destroy (state.common);
8456 fh_unref (state.prev_read_file);
8457 for (size_t i = 0; i < state.n_read_files; i++)
8458 read_file_destroy (state.read_files[i]);
8459 free (state.read_files);
8460 fh_unref (state.prev_write_file);
8461 for (size_t i = 0; i < state.n_write_files; i++)
8462 write_file_destroy (state.write_files[i]);
8463 free (state.write_files);
8464 fh_unref (state.prev_save_file);
8465 for (size_t i = 0; i < state.n_save_files; i++)
8466 save_file_destroy (state.save_files[i]);
8467 free (state.save_files);