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/misc.h"
53 #include "libpspp/str.h"
54 #include "libpspp/string-array.h"
55 #include "libpspp/stringi-set.h"
56 #include "libpspp/u8-line.h"
57 #include "math/distributions.h"
58 #include "math/random.h"
59 #include "output/driver.h"
60 #include "output/output-item.h"
61 #include "output/pivot-table.h"
63 #include "gl/c-ctype.h"
64 #include "gl/c-strcase.h"
65 #include "gl/ftoastr.h"
66 #include "gl/minmax.h"
70 #define _(msgid) gettext (msgid)
71 #define N_(msgid) (msgid)
75 struct hmap_node hmap_node;
83 struct file_handle *outfile;
84 struct string_array variables;
85 struct string_array fnames;
86 struct string_array snames;
89 /* Collects and owns factors and splits. The individual msave_command
90 structs point to these but do not own them. */
91 struct matrix_expr **factors;
92 size_t n_factors, allocated_factors;
93 struct matrix_expr **splits;
94 size_t n_splits, allocated_splits;
96 /* Execution state. */
97 struct dictionary *dict;
98 struct casewriter *writer;
103 struct file_handle *file;
104 struct dfm_reader *reader;
110 struct file_handle *file;
111 struct dfm_writer *writer;
113 struct u8_line *held;
118 struct file_handle *file;
119 struct dataset *dataset;
121 /* Parameters from parsing the first SAVE command for 'file'. */
122 struct string_array variables;
123 struct matrix_expr *names;
124 struct stringi_set strings;
126 /* Results from the first (only) attempt to open this save_file. */
128 struct casewriter *writer;
129 struct dictionary *dict;
134 struct dataset *dataset;
135 struct session *session;
139 struct msave_common *common;
141 struct file_handle *prev_read_file;
142 struct read_file **read_files;
145 struct file_handle *prev_write_file;
146 struct write_file **write_files;
147 size_t n_write_files;
149 struct file_handle *prev_save_file;
150 struct save_file **save_files;
154 static struct matrix_var *
155 matrix_var_lookup (struct matrix_state *s, struct substring name)
157 struct matrix_var *var;
159 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
160 utf8_hash_case_substring (name, 0), &s->vars)
161 if (!utf8_sscasecmp (ss_cstr (var->name), name))
167 static struct matrix_var *
168 matrix_var_create (struct matrix_state *s, struct substring name)
170 struct matrix_var *var = xmalloc (sizeof *var);
171 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
172 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
177 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
179 gsl_matrix_free (var->value);
183 /* The third argument to F() is a "prototype". For most prototypes, the first
184 letter (before the _) represents the return type and each other letter
185 (after the _) is an argument type. The types are:
187 - "m": A matrix of unrestricted dimensions.
191 - "v": A row or column vector.
193 - "e": Primarily for the first argument, this is a matrix with
194 unrestricted dimensions treated elementwise. Each element in the matrix
195 is passed to the implementation function separately.
197 The fourth argument is an optional constraints string. For this purpose the
198 first argument is named "a", the second "b", and so on. The following kinds
199 of constraints are supported. For matrix arguments, the constraints are
200 applied to each value in the matrix separately:
202 - "a(0,1)" or "a[0,1]": 0 < a < 1 or 0 <= a <= 1, respectively. Any
203 integer may substitute for 0 and 1. Half-open constraints (] and [) are
206 - "ai": Restrict a to integer values.
208 - "a>0", "a<0", "a>=0", "a<=0".
210 - "a<b", "a>b", "a<=b", "a>=b".
212 #define MATRIX_FUNCTIONS \
213 F(ABS, "ABS", m_e, NULL) \
214 F(ALL, "ALL", d_m, NULL) \
215 F(ANY, "ANY", d_m, NULL) \
216 F(ARSIN, "ARSIN", m_e, "a[-1,1]") \
217 F(ARTAN, "ARTAN", m_e, NULL) \
218 F(BLOCK, "BLOCK", m_any, NULL) \
219 F(CHOL, "CHOL", m_m, NULL) \
220 F(CMIN, "CMIN", m_m, NULL) \
221 F(CMAX, "CMAX", m_m, NULL) \
222 F(COS, "COS", m_e, NULL) \
223 F(CSSQ, "CSSQ", m_m, NULL) \
224 F(CSUM, "CSUM", m_m, NULL) \
225 F(DESIGN, "DESIGN", m_m, NULL) \
226 F(DET, "DET", d_m, NULL) \
227 F(DIAG, "DIAG", m_m, NULL) \
228 F(EVAL, "EVAL", m_m, NULL) \
229 F(EXP, "EXP", m_e, NULL) \
230 F(GINV, "GINV", m_m, NULL) \
231 F(GRADE, "GRADE", m_m, NULL) \
232 F(GSCH, "GSCH", m_m, NULL) \
233 F(IDENT, "IDENT", IDENT, NULL) \
234 F(INV, "INV", m_m, NULL) \
235 F(KRONEKER, "KRONEKER", m_mm, NULL) \
236 F(LG10, "LG10", m_e, "a>0") \
237 F(LN, "LN", m_e, "a>0") \
238 F(MAGIC, "MAGIC", m_d, "ai>=3") \
239 F(MAKE, "MAKE", m_ddd, NULL) \
240 F(MDIAG, "MDIAG", m_v, NULL) \
241 F(MMAX, "MMAX", d_m, NULL) \
242 F(MMIN, "MMIN", d_m, NULL) \
243 F(MOD, "MOD", m_md, NULL) \
244 F(MSSQ, "MSSQ", d_m, NULL) \
245 F(MSUM, "MSUM", d_m, NULL) \
246 F(NCOL, "NCOL", d_m, NULL) \
247 F(NROW, "NROW", d_m, NULL) \
248 F(RANK, "RANK", d_m, NULL) \
249 F(RESHAPE, "RESHAPE", m_mdd, NULL) \
250 F(RMAX, "RMAX", m_m, NULL) \
251 F(RMIN, "RMIN", m_m, NULL) \
252 F(RND, "RND", m_e, NULL) \
253 F(RNKORDER, "RNKORDER", m_m, NULL) \
254 F(RSSQ, "RSSQ", m_m, NULL) \
255 F(RSUM, "RSUM", m_m, NULL) \
256 F(SIN, "SIN", m_e, NULL) \
257 F(SOLVE, "SOLVE", m_mm, NULL) \
258 F(SQRT, "SQRT", m_e, "a>=0") \
259 F(SSCP, "SSCP", m_m, NULL) \
260 F(SVAL, "SVAL", m_m, NULL) \
261 F(SWEEP, "SWEEP", m_md, NULL) \
262 F(T, "T", m_m, NULL) \
263 F(TRACE, "TRACE", d_m, NULL) \
264 F(TRANSPOS, "TRANSPOS", m_m, NULL) \
265 F(TRUNC, "TRUNC", m_e, NULL) \
266 F(UNIFORM, "UNIFORM", m_dd, "ai>=0 bi>=0") \
267 F(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0") \
268 F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0") \
269 F(IDF_BETA, "IDF.BETA", m_edd, "a[0,1] b>0 c>0") \
270 F(RV_BETA, "RV.BETA", d_dd, "a>0 b>0") \
271 F(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
272 F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
273 F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]") \
274 F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]") \
275 F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0") \
276 F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0") \
277 F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "c>0") \
278 F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \
279 F(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0") \
280 F(CHICDF, "CHICDF", m_ed, "a>=0 b>0") \
281 F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0") \
282 F(PDF_CHISQ, "PDF.CHISQ", m_ed, "a>=0 b>0") \
283 F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \
284 F(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0") \
285 F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0") \
286 F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0") \
287 F(PDF_EXP, "PDF.EXP", m_ed, "a>=0 b>0") \
288 F(RV_EXP, "RV.EXP", d_d, "a>0") \
289 F(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0") \
290 F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \
291 F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0") \
292 F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0") \
293 F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0") \
294 F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0") \
295 F(RV_F, "RV.F", d_dd, "a>0 b>0") \
296 F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0") \
297 F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
298 F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0") \
299 F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
300 F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0") \
301 F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL) \
302 F(RV_LANDAU, "RV.LANDAU", d_none, NULL) \
303 F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0") \
304 F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0") \
305 F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0") \
306 F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0") \
307 F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]") \
308 F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]") \
309 F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0") \
310 F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0") \
311 F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0") \
312 F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0") \
313 F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
314 F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0") \
315 F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
316 F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0") \
317 F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0") \
318 F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0") \
319 F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0") \
320 F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0") \
321 F(CDFNORM, "CDFNORM", m_e, NULL) \
322 F(PROBIT, "PROBIT", m_e, "a(0,1)") \
323 F(NORMAL, "NORMAL", m_e, "a>0") \
324 F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0") \
325 F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0") \
326 F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0") \
327 F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0") \
328 F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0") \
329 F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0") \
330 F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0") \
331 F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0") \
332 F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0") \
333 F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0") \
334 F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL) \
335 F(RV_RTAIL, "RV.RTAIL", d_dd, NULL) \
336 F(CDF_T, "CDF.T", m_ed, "b>0") \
337 F(TCDF, "TCDF", m_ed, "b>0") \
338 F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0") \
339 F(PDF_T, "PDF.T", m_ed, "b>0") \
340 F(RV_T, "RV.T", d_d, "a>0") \
341 F(CDF_T1G, "CDF.T1G", m_edd, NULL) \
342 F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)") \
343 F(PDF_T1G, "PDF.T1G", m_edd, NULL) \
344 F(RV_T1G, "RV.T1G", d_dd, NULL) \
345 F(CDF_T2G, "CDF.T2G", m_edd, NULL) \
346 F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)") \
347 F(PDF_T2G, "PDF.T2G", m_edd, NULL) \
348 F(RV_T2G, "RV.T2G", d_dd, NULL) \
349 F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c") \
350 F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c") \
351 F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c") \
352 F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b") \
353 F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
354 F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0") \
355 F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
356 F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0") \
357 F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
358 F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
359 F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]") \
360 F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]") \
361 F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]") \
362 F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]") \
363 F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]") \
364 F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]") \
365 F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]") \
366 F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
367 F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
368 F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a") \
369 F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]") \
370 F(RV_LOG, "RV.LOG", d_d, "a(0,1]") \
371 F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
372 F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
373 F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]") \
374 F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0") \
375 F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0") \
376 F(RV_POISSON, "RV.POISSON", d_d, "a>0")
378 struct matrix_function_properties
381 const char *constraints;
384 enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 };
385 enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
386 enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
387 enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 };
388 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
389 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
390 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
391 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
392 enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 };
393 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
394 enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 };
395 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
396 enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
397 enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
398 enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
399 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
400 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
401 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
402 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
403 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
405 typedef double matrix_proto_d_none (void);
406 typedef double matrix_proto_d_d (double);
407 typedef double matrix_proto_d_dd (double, double);
408 typedef double matrix_proto_d_dd (double, double);
409 typedef double matrix_proto_d_ddd (double, double, double);
410 typedef gsl_matrix *matrix_proto_m_d (double);
411 typedef gsl_matrix *matrix_proto_m_dd (double, double);
412 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
413 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
414 typedef double matrix_proto_m_e (double);
415 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
416 typedef double matrix_proto_m_ed (double, double);
417 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
418 typedef double matrix_proto_m_edd (double, double, double);
419 typedef double matrix_proto_m_eddd (double, double, double, double);
420 typedef double matrix_proto_m_eed (double, double, double);
421 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
422 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
423 typedef double matrix_proto_d_m (gsl_matrix *);
424 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
425 typedef gsl_matrix *matrix_proto_IDENT (double, double);
427 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
428 static matrix_proto_##PROTO matrix_eval_##ENUM;
432 /* Matrix expressions. */
439 #define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
443 /* Elementwise and scalar arithmetic. */
444 MOP_NEGATE, /* unary - */
445 MOP_ADD_ELEMS, /* + */
446 MOP_SUB_ELEMS, /* - */
447 MOP_MUL_ELEMS, /* &* */
448 MOP_DIV_ELEMS, /* / and &/ */
449 MOP_EXP_ELEMS, /* &** */
451 MOP_SEQ_BY, /* a:b:c */
453 /* Matrix arithmetic. */
455 MOP_EXP_MAT, /* ** */
472 MOP_PASTE_HORZ, /* a, b, c, ... */
473 MOP_PASTE_VERT, /* a; b; c; ... */
477 MOP_VEC_INDEX, /* x(y) */
478 MOP_VEC_ALL, /* x(:) */
479 MOP_MAT_INDEX, /* x(y,z) */
480 MOP_ROW_INDEX, /* x(y,:) */
481 MOP_COL_INDEX, /* x(:,z) */
488 MOP_EOF, /* EOF('file') */
496 struct matrix_expr **subs;
501 struct matrix_var *variable;
502 struct read_file *eof;
505 struct msg_location *location;
509 matrix_expr_destroy (struct matrix_expr *e)
516 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
547 for (size_t i = 0; i < e->n_subs; i++)
548 matrix_expr_destroy (e->subs[i]);
557 msg_location_destroy (e->location);
561 static struct matrix_expr *
562 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
565 struct matrix_expr *e = xmalloc (sizeof *e);
566 *e = (struct matrix_expr) {
568 .subs = xmemdup (subs, n_subs * sizeof *subs),
572 for (size_t i = 0; i < n_subs; i++)
573 msg_location_merge (&e->location, subs[i]->location);
577 static struct matrix_expr *
578 matrix_expr_create_0 (enum matrix_op op)
580 struct matrix_expr *sub;
581 return matrix_expr_create_subs (op, &sub, 0);
584 static struct matrix_expr *
585 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
587 return matrix_expr_create_subs (op, &sub, 1);
590 static struct matrix_expr *
591 matrix_expr_create_2 (enum matrix_op op,
592 struct matrix_expr *sub0, struct matrix_expr *sub1)
594 struct matrix_expr *subs[] = { sub0, sub1 };
595 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
598 static struct matrix_expr *
599 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
600 struct matrix_expr *sub1, struct matrix_expr *sub2)
602 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
603 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
606 static struct matrix_expr *
607 matrix_expr_create_number (double number)
609 struct matrix_expr *e = xmalloc (sizeof *e);
610 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
614 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
616 static struct matrix_expr *
617 matrix_parse_curly_comma (struct matrix_state *s)
619 struct matrix_expr *lhs = matrix_parse_or_xor (s);
623 while (lex_match (s->lexer, T_COMMA))
625 struct matrix_expr *rhs = matrix_parse_or_xor (s);
628 matrix_expr_destroy (lhs);
631 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
636 static struct matrix_expr *
637 matrix_parse_curly_semi (struct matrix_state *s)
639 if (lex_token (s->lexer) == T_RCURLY)
640 return matrix_expr_create_0 (MOP_EMPTY);
642 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
646 while (lex_match (s->lexer, T_SEMICOLON))
648 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
651 matrix_expr_destroy (lhs);
654 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
659 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
660 for (size_t Y = 0; Y < (M)->size1; Y++) \
661 for (size_t X = 0; X < (M)->size2; X++) \
662 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
665 is_vector (const gsl_matrix *m)
667 return m->size1 <= 1 || m->size2 <= 1;
671 to_vector (gsl_matrix *m)
673 return (m->size1 == 1
674 ? gsl_matrix_row (m, 0).vector
675 : gsl_matrix_column (m, 0).vector);
680 matrix_eval_ABS (double d)
686 matrix_eval_ALL (gsl_matrix *m)
688 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
695 matrix_eval_ANY (gsl_matrix *m)
697 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
704 matrix_eval_ARSIN (double d)
710 matrix_eval_ARTAN (double d)
716 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
720 for (size_t i = 0; i < n; i++)
725 gsl_matrix *block = gsl_matrix_calloc (r, c);
727 for (size_t i = 0; i < n; i++)
729 for (size_t y = 0; y < m[i]->size1; y++)
730 for (size_t x = 0; x < m[i]->size2; x++)
731 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
739 matrix_eval_CHOL (gsl_matrix *m)
741 if (!gsl_linalg_cholesky_decomp1 (m))
743 for (size_t y = 0; y < m->size1; y++)
744 for (size_t x = y + 1; x < m->size2; x++)
745 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
747 for (size_t y = 0; y < m->size1; y++)
748 for (size_t x = 0; x < y; x++)
749 gsl_matrix_set (m, y, x, 0);
754 msg (SE, _("Input to CHOL function is not positive-definite."));
760 matrix_eval_col_extremum (gsl_matrix *m, bool min)
765 return gsl_matrix_alloc (1, 0);
767 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
768 for (size_t x = 0; x < m->size2; x++)
770 double ext = gsl_matrix_get (m, 0, x);
771 for (size_t y = 1; y < m->size1; y++)
773 double value = gsl_matrix_get (m, y, x);
774 if (min ? value < ext : value > ext)
777 gsl_matrix_set (cext, 0, x, ext);
783 matrix_eval_CMAX (gsl_matrix *m)
785 return matrix_eval_col_extremum (m, false);
789 matrix_eval_CMIN (gsl_matrix *m)
791 return matrix_eval_col_extremum (m, true);
795 matrix_eval_COS (double d)
801 matrix_eval_col_sum (gsl_matrix *m, bool square)
806 return gsl_matrix_alloc (1, 0);
808 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
809 for (size_t x = 0; x < m->size2; x++)
812 for (size_t y = 0; y < m->size1; y++)
814 double d = gsl_matrix_get (m, y, x);
815 sum += square ? pow2 (d) : d;
817 gsl_matrix_set (result, 0, x, sum);
823 matrix_eval_CSSQ (gsl_matrix *m)
825 return matrix_eval_col_sum (m, true);
829 matrix_eval_CSUM (gsl_matrix *m)
831 return matrix_eval_col_sum (m, false);
835 compare_double_3way (const void *a_, const void *b_)
837 const double *a = a_;
838 const double *b = b_;
839 return *a < *b ? -1 : *a > *b;
843 matrix_eval_DESIGN (gsl_matrix *m)
845 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
846 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
847 gsl_matrix_transpose_memcpy (&m2, m);
849 for (size_t y = 0; y < m2.size1; y++)
850 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
852 size_t *n = xcalloc (m2.size1, sizeof *n);
854 for (size_t i = 0; i < m2.size1; i++)
856 double *row = tmp + m2.size2 * i;
857 for (size_t j = 0; j < m2.size2; )
860 for (k = j + 1; k < m2.size2; k++)
861 if (row[j] != row[k])
863 row[n[i]++] = row[j];
868 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
873 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
875 for (size_t i = 0; i < m->size2; i++)
880 const double *unique = tmp + m2.size2 * i;
881 for (size_t j = 0; j < n[i]; j++, x++)
883 double value = unique[j];
884 for (size_t y = 0; y < m->size1; y++)
885 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
896 matrix_eval_DET (gsl_matrix *m)
898 gsl_permutation *p = gsl_permutation_alloc (m->size1);
900 gsl_linalg_LU_decomp (m, p, &signum);
901 gsl_permutation_free (p);
902 return gsl_linalg_LU_det (m, signum);
906 matrix_eval_DIAG (gsl_matrix *m)
908 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
909 for (size_t i = 0; i < diag->size1; i++)
910 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
915 is_symmetric (const gsl_matrix *m)
917 if (m->size1 != m->size2)
920 for (size_t y = 0; y < m->size1; y++)
921 for (size_t x = 0; x < y; x++)
922 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
929 compare_double_desc (const void *a_, const void *b_)
931 const double *a = a_;
932 const double *b = b_;
933 return *a > *b ? -1 : *a < *b;
937 matrix_eval_EVAL (gsl_matrix *m)
939 if (!is_symmetric (m))
941 msg (SE, _("Argument of EVAL must be symmetric."));
945 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
946 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
947 gsl_vector v_eval = to_vector (eval);
948 gsl_eigen_symm (m, &v_eval, w);
949 gsl_eigen_symm_free (w);
951 assert (v_eval.stride == 1);
952 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
958 matrix_eval_EXP (double d)
963 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
966 Charl Linssen <charl@itfromb.it>
970 matrix_eval_GINV (gsl_matrix *A)
975 gsl_matrix *tmp_mat = NULL;
978 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
979 tmp_mat = gsl_matrix_alloc (m, n);
980 gsl_matrix_transpose_memcpy (tmp_mat, A);
988 gsl_matrix *V = gsl_matrix_alloc (m, m);
989 gsl_vector *u = gsl_vector_alloc (m);
991 gsl_vector *tmp_vec = gsl_vector_alloc (m);
992 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
993 gsl_vector_free (tmp_vec);
996 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
997 gsl_matrix_set_zero (Sigma_pinv);
998 double cutoff = 1e-15 * gsl_vector_max (u);
1000 for (size_t i = 0; i < m; ++i)
1002 double x = gsl_vector_get (u, i);
1003 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
1006 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
1007 gsl_matrix *U = gsl_matrix_calloc (n, n);
1008 for (size_t i = 0; i < n; i++)
1009 for (size_t j = 0; j < m; j++)
1010 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
1012 /* two dot products to obtain pseudoinverse */
1013 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
1014 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
1019 A_pinv = gsl_matrix_alloc (n, m);
1020 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
1024 A_pinv = gsl_matrix_alloc (m, n);
1025 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
1028 gsl_matrix_free (tmp_mat);
1029 gsl_matrix_free (tmp_mat2);
1030 gsl_matrix_free (U);
1031 gsl_matrix_free (Sigma_pinv);
1032 gsl_vector_free (u);
1033 gsl_matrix_free (V);
1045 grade_compare_3way (const void *a_, const void *b_)
1047 const struct grade *a = a_;
1048 const struct grade *b = b_;
1050 return (a->value < b->value ? -1
1051 : a->value > b->value ? 1
1059 matrix_eval_GRADE (gsl_matrix *m)
1061 size_t n = m->size1 * m->size2;
1062 struct grade *grades = xmalloc (n * sizeof *grades);
1065 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1066 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
1067 qsort (grades, n, sizeof *grades, grade_compare_3way);
1069 for (size_t i = 0; i < n; i++)
1070 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
1078 dot (gsl_vector *a, gsl_vector *b)
1080 double result = 0.0;
1081 for (size_t i = 0; i < a->size; i++)
1082 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
1087 norm2 (gsl_vector *v)
1089 double result = 0.0;
1090 for (size_t i = 0; i < v->size; i++)
1091 result += pow2 (gsl_vector_get (v, i));
1096 norm (gsl_vector *v)
1098 return sqrt (norm2 (v));
1102 matrix_eval_GSCH (gsl_matrix *v)
1104 if (v->size2 < v->size1)
1106 msg (SE, _("GSCH requires its argument to have at least as many columns "
1107 "as rows, but it has dimensions %zu×%zu."),
1108 v->size1, v->size2);
1111 if (!v->size1 || !v->size2)
1114 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
1116 for (size_t vx = 0; vx < v->size2; vx++)
1118 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
1119 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
1121 gsl_vector_memcpy (&u_i, &v_i);
1122 for (size_t j = 0; j < ux; j++)
1124 gsl_vector u_j = gsl_matrix_column (u, j).vector;
1125 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
1126 for (size_t k = 0; k < u_i.size; k++)
1127 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
1128 - scale * gsl_vector_get (&u_j, k)));
1131 double len = norm (&u_i);
1134 gsl_vector_scale (&u_i, 1.0 / len);
1135 if (++ux >= v->size1)
1142 msg (SE, _("%zu×%zu argument to GSCH contains only "
1143 "%zu linearly independent columns."),
1144 v->size1, v->size2, ux);
1145 gsl_matrix_free (u);
1149 u->size2 = v->size1;
1154 matrix_eval_IDENT (double s1, double s2)
1156 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
1158 msg (SE, _("Arguments to IDENT must be non-negative integers."));
1162 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
1163 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1169 invert_matrix (gsl_matrix *x)
1171 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1173 gsl_linalg_LU_decomp (x, p, &signum);
1174 gsl_linalg_LU_invx (x, p);
1175 gsl_permutation_free (p);
1179 matrix_eval_INV (gsl_matrix *m)
1186 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1188 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1189 a->size2 * b->size2);
1191 for (size_t ar = 0; ar < a->size1; ar++)
1192 for (size_t br = 0; br < b->size1; br++, y++)
1195 for (size_t ac = 0; ac < a->size2; ac++)
1196 for (size_t bc = 0; bc < b->size2; bc++, x++)
1198 double av = gsl_matrix_get (a, ar, ac);
1199 double bv = gsl_matrix_get (b, br, bc);
1200 gsl_matrix_set (k, y, x, av * bv);
1207 matrix_eval_LG10 (double d)
1213 matrix_eval_LN (double d)
1219 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1221 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1224 for (size_t i = 1; i <= n * n; i++)
1226 gsl_matrix_set (m, y, x, i);
1228 size_t y1 = !y ? n - 1 : y - 1;
1229 size_t x1 = x + 1 >= n ? 0 : x + 1;
1230 if (gsl_matrix_get (m, y1, x1) == 0)
1236 y = y + 1 >= n ? 0 : y + 1;
1241 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1243 double a = gsl_matrix_get (m, y1, x1);
1244 double b = gsl_matrix_get (m, y2, x2);
1245 gsl_matrix_set (m, y1, x1, b);
1246 gsl_matrix_set (m, y2, x2, a);
1250 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1254 /* A. Umar, "On the Construction of Even Order Magic Squares",
1255 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1257 for (size_t i = 1; i <= n * n / 2; i++)
1259 gsl_matrix_set (m, y, x, i);
1269 for (size_t i = n * n; i > n * n / 2; i--)
1271 gsl_matrix_set (m, y, x, i);
1279 for (size_t y = 0; y < n; y++)
1280 for (size_t x = 0; x < n / 2; x++)
1282 unsigned int d = gsl_matrix_get (m, y, x);
1283 if (d % 2 != (y < n / 2))
1284 magic_exchange (m, y, x, y, n - x - 1);
1289 size_t x1 = n / 2 - 1;
1291 magic_exchange (m, y1, x1, y2, x1);
1292 magic_exchange (m, y1, x2, y2, x2);
1296 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1298 /* A. Umar, "On the Construction of Even Order Magic Squares",
1299 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1303 for (size_t i = 1; ; i++)
1305 gsl_matrix_set (m, y, x, i);
1306 if (++y == n / 2 - 1)
1318 for (size_t i = n * n; ; i--)
1320 gsl_matrix_set (m, y, x, i);
1321 if (++y == n / 2 - 1)
1330 for (size_t y = 0; y < n; y++)
1331 if (y != n / 2 - 1 && y != n / 2)
1332 for (size_t x = 0; x < n / 2; x++)
1334 unsigned int d = gsl_matrix_get (m, y, x);
1335 if (d % 2 != (y < n / 2))
1336 magic_exchange (m, y, x, y, n - x - 1);
1339 size_t a0 = (n * n - 2 * n) / 2 + 1;
1340 for (size_t i = 0; i < n / 2; i++)
1343 gsl_matrix_set (m, n / 2, i, a);
1344 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1346 for (size_t i = 0; i < n / 2; i++)
1348 size_t a = a0 + i + n / 2;
1349 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1350 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1352 for (size_t x = 1; x < n / 2; x += 2)
1353 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1354 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1355 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1356 size_t x1 = n / 2 - 2;
1357 size_t x2 = n / 2 + 1;
1358 size_t y1 = n / 2 - 2;
1359 size_t y2 = n / 2 + 1;
1360 magic_exchange (m, y1, x1, y2, x1);
1361 magic_exchange (m, y1, x2, y2, x2);
1365 matrix_eval_MAGIC (double n_)
1369 gsl_matrix *m = gsl_matrix_calloc (n, n);
1371 matrix_eval_MAGIC_odd (m, n);
1373 matrix_eval_MAGIC_singly_even (m, n);
1375 matrix_eval_MAGIC_doubly_even (m, n);
1380 matrix_eval_MAKE (double r, double c, double value)
1382 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1384 msg (SE, _("Size arguments to MAKE must be integers."));
1388 gsl_matrix *m = gsl_matrix_alloc (r, c);
1389 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1395 matrix_eval_MDIAG (gsl_vector *v)
1397 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1398 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1399 gsl_vector_memcpy (&diagonal, v);
1404 matrix_eval_MMAX (gsl_matrix *m)
1406 return gsl_matrix_max (m);
1410 matrix_eval_MMIN (gsl_matrix *m)
1412 return gsl_matrix_min (m);
1416 matrix_eval_MOD (gsl_matrix *m, double divisor)
1420 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1424 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1425 *d = fmod (*d, divisor);
1430 matrix_eval_MSSQ (gsl_matrix *m)
1433 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1439 matrix_eval_MSUM (gsl_matrix *m)
1442 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1448 matrix_eval_NCOL (gsl_matrix *m)
1454 matrix_eval_NROW (gsl_matrix *m)
1460 matrix_eval_RANK (gsl_matrix *m)
1462 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1463 gsl_linalg_QR_decomp (m, tau);
1464 gsl_vector_free (tau);
1466 return gsl_linalg_QRPT_rank (m, -1);
1470 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1472 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1474 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1479 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1481 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1482 "product of matrix dimensions."));
1486 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1489 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1491 gsl_matrix_set (dst, y1, x1, *d);
1502 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1507 return gsl_matrix_alloc (0, 1);
1509 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1510 for (size_t y = 0; y < m->size1; y++)
1512 double ext = gsl_matrix_get (m, y, 0);
1513 for (size_t x = 1; x < m->size2; x++)
1515 double value = gsl_matrix_get (m, y, x);
1516 if (min ? value < ext : value > ext)
1519 gsl_matrix_set (rext, y, 0, ext);
1525 matrix_eval_RMAX (gsl_matrix *m)
1527 return matrix_eval_row_extremum (m, false);
1531 matrix_eval_RMIN (gsl_matrix *m)
1533 return matrix_eval_row_extremum (m, true);
1537 matrix_eval_RND (double d)
1549 rank_compare_3way (const void *a_, const void *b_)
1551 const struct rank *a = a_;
1552 const struct rank *b = b_;
1554 return a->value < b->value ? -1 : a->value > b->value;
1558 matrix_eval_RNKORDER (gsl_matrix *m)
1560 size_t n = m->size1 * m->size2;
1561 struct rank *ranks = xmalloc (n * sizeof *ranks);
1563 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1564 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1565 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1567 for (size_t i = 0; i < n; )
1570 for (j = i + 1; j < n; j++)
1571 if (ranks[i].value != ranks[j].value)
1574 double rank = (i + j + 1.0) / 2.0;
1575 for (size_t k = i; k < j; k++)
1576 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1587 matrix_eval_row_sum (gsl_matrix *m, bool square)
1592 return gsl_matrix_alloc (0, 1);
1594 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1595 for (size_t y = 0; y < m->size1; y++)
1598 for (size_t x = 0; x < m->size2; x++)
1600 double d = gsl_matrix_get (m, y, x);
1601 sum += square ? pow2 (d) : d;
1603 gsl_matrix_set (result, y, 0, sum);
1609 matrix_eval_RSSQ (gsl_matrix *m)
1611 return matrix_eval_row_sum (m, true);
1615 matrix_eval_RSUM (gsl_matrix *m)
1617 return matrix_eval_row_sum (m, false);
1621 matrix_eval_SIN (double d)
1627 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1629 if (m1->size1 != m2->size1)
1631 msg (SE, _("SOLVE requires its arguments to have the same number of "
1632 "rows, but the first argument has dimensions %zu×%zu and "
1633 "the second %zu×%zu."),
1634 m1->size1, m1->size2,
1635 m2->size1, m2->size2);
1639 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1640 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1642 gsl_linalg_LU_decomp (m1, p, &signum);
1643 for (size_t i = 0; i < m2->size2; i++)
1645 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1646 gsl_vector xi = gsl_matrix_column (x, i).vector;
1647 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1649 gsl_permutation_free (p);
1654 matrix_eval_SQRT (double d)
1660 matrix_eval_SSCP (gsl_matrix *m)
1662 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1663 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1668 matrix_eval_SVAL (gsl_matrix *m)
1670 gsl_matrix *tmp_mat = NULL;
1671 if (m->size2 > m->size1)
1673 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1674 gsl_matrix_transpose_memcpy (tmp_mat, m);
1679 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1680 gsl_vector *S = gsl_vector_alloc (m->size2);
1681 gsl_vector *work = gsl_vector_alloc (m->size2);
1682 gsl_linalg_SV_decomp (m, V, S, work);
1684 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1685 for (size_t i = 0; i < m->size2; i++)
1686 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1688 gsl_matrix_free (V);
1689 gsl_vector_free (S);
1690 gsl_vector_free (work);
1691 gsl_matrix_free (tmp_mat);
1697 matrix_eval_SWEEP (gsl_matrix *m, double d)
1699 if (d < 1 || d > SIZE_MAX)
1701 msg (SE, _("Scalar argument to SWEEP must be integer."));
1705 if (k >= MIN (m->size1, m->size2))
1707 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1708 "equal to the smaller of the matrix argument's rows and "
1713 double m_kk = gsl_matrix_get (m, k, k);
1714 if (fabs (m_kk) > 1e-19)
1716 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1717 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1719 double m_ij = gsl_matrix_get (m, i, j);
1720 double m_ik = gsl_matrix_get (m, i, k);
1721 double m_kj = gsl_matrix_get (m, k, j);
1722 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1731 for (size_t i = 0; i < m->size1; i++)
1733 gsl_matrix_set (m, i, k, 0);
1734 gsl_matrix_set (m, k, i, 0);
1741 matrix_eval_TRACE (gsl_matrix *m)
1744 size_t n = MIN (m->size1, m->size2);
1745 for (size_t i = 0; i < n; i++)
1746 sum += gsl_matrix_get (m, i, i);
1751 matrix_eval_T (gsl_matrix *m)
1753 return matrix_eval_TRANSPOS (m);
1757 matrix_eval_TRANSPOS (gsl_matrix *m)
1759 if (m->size1 == m->size2)
1761 gsl_matrix_transpose (m);
1766 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1767 gsl_matrix_transpose_memcpy (t, m);
1773 matrix_eval_TRUNC (double d)
1779 matrix_eval_UNIFORM (double r_, double c_)
1783 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1785 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1789 gsl_matrix *m = gsl_matrix_alloc (r, c);
1790 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1791 *d = gsl_ran_flat (get_rng (), 0, 1);
1796 matrix_eval_PDF_BETA (double x, double a, double b)
1798 return gsl_ran_beta_pdf (x, a, b);
1802 matrix_eval_CDF_BETA (double x, double a, double b)
1804 return gsl_cdf_beta_P (x, a, b);
1808 matrix_eval_IDF_BETA (double P, double a, double b)
1810 return gsl_cdf_beta_Pinv (P, a, b);
1814 matrix_eval_RV_BETA (double a, double b)
1816 return gsl_ran_beta (get_rng (), a, b);
1820 matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
1822 return ncdf_beta (x, a, b, lambda);
1826 matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
1828 return npdf_beta (x, a, b, lambda);
1832 matrix_eval_CDF_BVNOR (double x0, double x1, double r)
1834 return cdf_bvnor (x0, x1, r);
1838 matrix_eval_PDF_BVNOR (double x0, double x1, double r)
1840 return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
1844 matrix_eval_CDF_CAUCHY (double x, double a, double b)
1846 return gsl_cdf_cauchy_P ((x - a) / b, 1);
1850 matrix_eval_IDF_CAUCHY (double P, double a, double b)
1852 return a + b * gsl_cdf_cauchy_Pinv (P, 1);
1856 matrix_eval_PDF_CAUCHY (double x, double a, double b)
1858 return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
1862 matrix_eval_RV_CAUCHY (double a, double b)
1864 return a + b * gsl_ran_cauchy (get_rng (), 1);
1868 matrix_eval_CDF_CHISQ (double x, double df)
1870 return gsl_cdf_chisq_P (x, df);
1874 matrix_eval_CHICDF (double x, double df)
1876 return matrix_eval_CDF_CHISQ (x, df);
1880 matrix_eval_IDF_CHISQ (double P, double df)
1882 return gsl_cdf_chisq_Pinv (P, df);
1886 matrix_eval_PDF_CHISQ (double x, double df)
1888 return gsl_ran_chisq_pdf (x, df);
1892 matrix_eval_RV_CHISQ (double df)
1894 return gsl_ran_chisq (get_rng (), df);
1898 matrix_eval_SIG_CHISQ (double x, double df)
1900 return gsl_cdf_chisq_Q (x, df);
1904 matrix_eval_CDF_EXP (double x, double a)
1906 return gsl_cdf_exponential_P (x, 1. / a);
1910 matrix_eval_IDF_EXP (double P, double a)
1912 return gsl_cdf_exponential_Pinv (P, 1. / a);
1916 matrix_eval_PDF_EXP (double x, double a)
1918 return gsl_ran_exponential_pdf (x, 1. / a);
1922 matrix_eval_RV_EXP (double a)
1924 return gsl_ran_exponential (get_rng (), 1. / a);
1928 matrix_eval_PDF_XPOWER (double x, double a, double b)
1930 return gsl_ran_exppow_pdf (x, a, b);
1934 matrix_eval_RV_XPOWER (double a, double b)
1936 return gsl_ran_exppow (get_rng (), a, b);
1940 matrix_eval_CDF_F (double x, double df1, double df2)
1942 return gsl_cdf_fdist_P (x, df1, df2);
1946 matrix_eval_FCDF (double x, double df1, double df2)
1948 return matrix_eval_CDF_F (x, df1, df2);
1952 matrix_eval_IDF_F (double P, double df1, double df2)
1954 return idf_fdist (P, df1, df2);
1958 matrix_eval_RV_F (double df1, double df2)
1960 return gsl_ran_fdist (get_rng (), df1, df2);
1964 matrix_eval_PDF_F (double x, double df1, double df2)
1966 return gsl_ran_fdist_pdf (x, df1, df2);
1970 matrix_eval_SIG_F (double x, double df1, double df2)
1972 return gsl_cdf_fdist_Q (x, df1, df2);
1976 matrix_eval_CDF_GAMMA (double x, double a, double b)
1978 return gsl_cdf_gamma_P (x, a, 1. / b);
1982 matrix_eval_IDF_GAMMA (double P, double a, double b)
1984 return gsl_cdf_gamma_Pinv (P, a, 1. / b);
1988 matrix_eval_PDF_GAMMA (double x, double a, double b)
1990 return gsl_ran_gamma_pdf (x, a, 1. / b);
1994 matrix_eval_RV_GAMMA (double a, double b)
1996 return gsl_ran_gamma (get_rng (), a, 1. / b);
2000 matrix_eval_PDF_LANDAU (double x)
2002 return gsl_ran_landau_pdf (x);
2006 matrix_eval_RV_LANDAU (void)
2008 return gsl_ran_landau (get_rng ());
2012 matrix_eval_CDF_LAPLACE (double x, double a, double b)
2014 return gsl_cdf_laplace_P ((x - a) / b, 1);
2018 matrix_eval_IDF_LAPLACE (double P, double a, double b)
2020 return a + b * gsl_cdf_laplace_Pinv (P, 1);
2024 matrix_eval_PDF_LAPLACE (double x, double a, double b)
2026 return gsl_ran_laplace_pdf ((x - a) / b, 1);
2030 matrix_eval_RV_LAPLACE (double a, double b)
2032 return a + b * gsl_ran_laplace (get_rng (), 1);
2036 matrix_eval_RV_LEVY (double c, double alpha)
2038 return gsl_ran_levy (get_rng (), c, alpha);
2042 matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
2044 return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
2048 matrix_eval_CDF_LOGISTIC (double x, double a, double b)
2050 return gsl_cdf_logistic_P ((x - a) / b, 1);
2054 matrix_eval_IDF_LOGISTIC (double P, double a, double b)
2056 return a + b * gsl_cdf_logistic_Pinv (P, 1);
2060 matrix_eval_PDF_LOGISTIC (double x, double a, double b)
2062 return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
2066 matrix_eval_RV_LOGISTIC (double a, double b)
2068 return a + b * gsl_ran_logistic (get_rng (), 1);
2072 matrix_eval_CDF_LNORMAL (double x, double m, double s)
2074 return gsl_cdf_lognormal_P (x, log (m), s);
2078 matrix_eval_IDF_LNORMAL (double P, double m, double s)
2080 return gsl_cdf_lognormal_Pinv (P, log (m), s);;
2084 matrix_eval_PDF_LNORMAL (double x, double m, double s)
2086 return gsl_ran_lognormal_pdf (x, log (m), s);
2090 matrix_eval_RV_LNORMAL (double m, double s)
2092 return gsl_ran_lognormal (get_rng (), log (m), s);
2096 matrix_eval_CDF_NORMAL (double x, double u, double s)
2098 return gsl_cdf_gaussian_P (x - u, s);
2102 matrix_eval_IDF_NORMAL (double P, double u, double s)
2104 return u + gsl_cdf_gaussian_Pinv (P, s);
2108 matrix_eval_PDF_NORMAL (double x, double u, double s)
2110 return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
2114 matrix_eval_RV_NORMAL (double u, double s)
2116 return u + gsl_ran_gaussian (get_rng (), s);
2120 matrix_eval_CDFNORM (double x)
2122 return gsl_cdf_ugaussian_P (x);
2126 matrix_eval_PROBIT (double P)
2128 return gsl_cdf_ugaussian_Pinv (P);
2132 matrix_eval_NORMAL (double s)
2134 return gsl_ran_gaussian (get_rng (), s);
2138 matrix_eval_PDF_NTAIL (double x, double a, double sigma)
2140 return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
2144 matrix_eval_RV_NTAIL (double a, double sigma)
2146 return gsl_ran_gaussian_tail (get_rng (), a, sigma);
2150 matrix_eval_CDF_PARETO (double x, double a, double b)
2152 return gsl_cdf_pareto_P (x, b, a);
2156 matrix_eval_IDF_PARETO (double P, double a, double b)
2158 return gsl_cdf_pareto_Pinv (P, b, a);
2162 matrix_eval_PDF_PARETO (double x, double a, double b)
2164 return gsl_ran_pareto_pdf (x, b, a);
2168 matrix_eval_RV_PARETO (double a, double b)
2170 return gsl_ran_pareto (get_rng (), b, a);
2174 matrix_eval_CDF_RAYLEIGH (double x, double sigma)
2176 return gsl_cdf_rayleigh_P (x, sigma);
2180 matrix_eval_IDF_RAYLEIGH (double P, double sigma)
2182 return gsl_cdf_rayleigh_Pinv (P, sigma);
2186 matrix_eval_PDF_RAYLEIGH (double x, double sigma)
2188 return gsl_ran_rayleigh_pdf (x, sigma);
2192 matrix_eval_RV_RAYLEIGH (double sigma)
2194 return gsl_ran_rayleigh (get_rng (), sigma);
2198 matrix_eval_PDF_RTAIL (double x, double a, double sigma)
2200 return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
2204 matrix_eval_RV_RTAIL (double a, double sigma)
2206 return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
2210 matrix_eval_CDF_T (double x, double df)
2212 return gsl_cdf_tdist_P (x, df);
2216 matrix_eval_TCDF (double x, double df)
2218 return matrix_eval_CDF_T (x, df);
2222 matrix_eval_IDF_T (double P, double df)
2224 return gsl_cdf_tdist_Pinv (P, df);
2228 matrix_eval_PDF_T (double x, double df)
2230 return gsl_ran_tdist_pdf (x, df);
2234 matrix_eval_RV_T (double df)
2236 return gsl_ran_tdist (get_rng (), df);
2240 matrix_eval_CDF_T1G (double x, double a, double b)
2242 return gsl_cdf_gumbel1_P (x, a, b);
2246 matrix_eval_IDF_T1G (double P, double a, double b)
2248 return gsl_cdf_gumbel1_Pinv (P, a, b);
2252 matrix_eval_PDF_T1G (double x, double a, double b)
2254 return gsl_ran_gumbel1_pdf (x, a, b);
2258 matrix_eval_RV_T1G (double a, double b)
2260 return gsl_ran_gumbel1 (get_rng (), a, b);
2264 matrix_eval_CDF_T2G (double x, double a, double b)
2266 return gsl_cdf_gumbel1_P (x, a, b);
2270 matrix_eval_IDF_T2G (double P, double a, double b)
2272 return gsl_cdf_gumbel1_Pinv (P, a, b);
2276 matrix_eval_PDF_T2G (double x, double a, double b)
2278 return gsl_ran_gumbel1_pdf (x, a, b);
2282 matrix_eval_RV_T2G (double a, double b)
2284 return gsl_ran_gumbel1 (get_rng (), a, b);
2288 matrix_eval_CDF_UNIFORM (double x, double a, double b)
2290 return gsl_cdf_flat_P (x, a, b);
2294 matrix_eval_IDF_UNIFORM (double P, double a, double b)
2296 return gsl_cdf_flat_Pinv (P, a, b);
2300 matrix_eval_PDF_UNIFORM (double x, double a, double b)
2302 return gsl_ran_flat_pdf (x, a, b);
2306 matrix_eval_RV_UNIFORM (double a, double b)
2308 return gsl_ran_flat (get_rng (), a, b);
2312 matrix_eval_CDF_WEIBULL (double x, double a, double b)
2314 return gsl_cdf_weibull_P (x, a, b);
2318 matrix_eval_IDF_WEIBULL (double P, double a, double b)
2320 return gsl_cdf_weibull_Pinv (P, a, b);
2324 matrix_eval_PDF_WEIBULL (double x, double a, double b)
2326 return gsl_ran_weibull_pdf (x, a, b);
2330 matrix_eval_RV_WEIBULL (double a, double b)
2332 return gsl_ran_weibull (get_rng (), a, b);
2336 matrix_eval_CDF_BERNOULLI (double k, double p)
2338 return k ? 1 : 1 - p;
2342 matrix_eval_PDF_BERNOULLI (double k, double p)
2344 return gsl_ran_bernoulli_pdf (k, p);
2348 matrix_eval_RV_BERNOULLI (double p)
2350 return gsl_ran_bernoulli (get_rng (), p);
2354 matrix_eval_CDF_BINOM (double k, double n, double p)
2356 return gsl_cdf_binomial_P (k, p, n);
2360 matrix_eval_PDF_BINOM (double k, double n, double p)
2362 return gsl_ran_binomial_pdf (k, p, n);
2366 matrix_eval_RV_BINOM (double n, double p)
2368 return gsl_ran_binomial (get_rng (), p, n);
2372 matrix_eval_CDF_GEOM (double k, double p)
2374 return gsl_cdf_geometric_P (k, p);
2378 matrix_eval_PDF_GEOM (double k, double p)
2380 return gsl_ran_geometric_pdf (k, p);
2384 matrix_eval_RV_GEOM (double p)
2386 return gsl_ran_geometric (get_rng (), p);
2390 matrix_eval_CDF_HYPER (double k, double a, double b, double c)
2392 return gsl_cdf_hypergeometric_P (k, c, a - c, b);
2396 matrix_eval_PDF_HYPER (double k, double a, double b, double c)
2398 return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
2402 matrix_eval_RV_HYPER (double a, double b, double c)
2404 return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
2408 matrix_eval_PDF_LOG (double k, double p)
2410 return gsl_ran_logarithmic_pdf (k, p);
2414 matrix_eval_RV_LOG (double p)
2416 return gsl_ran_logarithmic (get_rng (), p);
2420 matrix_eval_CDF_NEGBIN (double k, double n, double p)
2422 return gsl_cdf_negative_binomial_P (k, p, n);
2426 matrix_eval_PDF_NEGBIN (double k, double n, double p)
2428 return gsl_ran_negative_binomial_pdf (k, p, n);
2432 matrix_eval_RV_NEGBIN (double n, double p)
2434 return gsl_ran_negative_binomial (get_rng (), p, n);
2438 matrix_eval_CDF_POISSON (double k, double mu)
2440 return gsl_cdf_poisson_P (k, mu);
2444 matrix_eval_PDF_POISSON (double k, double mu)
2446 return gsl_ran_poisson_pdf (k, mu);
2450 matrix_eval_RV_POISSON (double mu)
2452 return gsl_ran_poisson (get_rng (), mu);
2455 struct matrix_function
2459 size_t min_args, max_args;
2462 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
2465 word_matches (const char **test, const char **name)
2467 size_t test_len = strcspn (*test, ".");
2468 size_t name_len = strcspn (*name, ".");
2469 if (test_len == name_len)
2471 if (buf_compare_case (*test, *name, test_len))
2474 else if (test_len < 3 || test_len > name_len)
2478 if (buf_compare_case (*test, *name, test_len))
2484 if (**test != **name)
2495 /* Returns 0 if TOKEN and FUNC do not match,
2496 1 if TOKEN is an acceptable abbreviation for FUNC,
2497 2 if TOKEN equals FUNC. */
2499 compare_function_names (const char *token_, const char *func_)
2501 const char *token = token_;
2502 const char *func = func_;
2503 while (*token || *func)
2504 if (!word_matches (&token, &func))
2506 return !c_strcasecmp (token_, func_) ? 2 : 1;
2509 static const struct matrix_function *
2510 matrix_parse_function_name (const char *token)
2512 static const struct matrix_function functions[] =
2514 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
2515 { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
2519 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
2521 for (size_t i = 0; i < N_FUNCTIONS; i++)
2523 if (compare_function_names (token, functions[i].name) > 0)
2524 return &functions[i];
2529 static struct read_file *
2530 read_file_create (struct matrix_state *s, struct file_handle *fh)
2532 for (size_t i = 0; i < s->n_read_files; i++)
2534 struct read_file *rf = s->read_files[i];
2542 struct read_file *rf = xmalloc (sizeof *rf);
2543 *rf = (struct read_file) { .file = fh };
2545 s->read_files = xrealloc (s->read_files,
2546 (s->n_read_files + 1) * sizeof *s->read_files);
2547 s->read_files[s->n_read_files++] = rf;
2551 static struct dfm_reader *
2552 read_file_open (struct read_file *rf)
2555 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
2560 read_file_destroy (struct read_file *rf)
2564 fh_unref (rf->file);
2565 dfm_close_reader (rf->reader);
2566 free (rf->encoding);
2571 static struct write_file *
2572 write_file_create (struct matrix_state *s, struct file_handle *fh)
2574 for (size_t i = 0; i < s->n_write_files; i++)
2576 struct write_file *wf = s->write_files[i];
2584 struct write_file *wf = xmalloc (sizeof *wf);
2585 *wf = (struct write_file) { .file = fh };
2587 s->write_files = xrealloc (s->write_files,
2588 (s->n_write_files + 1) * sizeof *s->write_files);
2589 s->write_files[s->n_write_files++] = wf;
2593 static struct dfm_writer *
2594 write_file_open (struct write_file *wf)
2597 wf->writer = dfm_open_writer (wf->file, wf->encoding);
2602 write_file_destroy (struct write_file *wf)
2608 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
2609 wf->held->s.ss.length);
2610 u8_line_destroy (wf->held);
2614 fh_unref (wf->file);
2615 dfm_close_writer (wf->writer);
2616 free (wf->encoding);
2622 matrix_parse_function (struct matrix_state *s, const char *token,
2623 struct matrix_expr **exprp)
2626 if (lex_next_token (s->lexer, 1) != T_LPAREN)
2629 if (lex_match_id (s->lexer, "EOF"))
2632 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
2636 if (!lex_force_match (s->lexer, T_RPAREN))
2642 struct read_file *rf = read_file_create (s, fh);
2644 struct matrix_expr *e = xmalloc (sizeof *e);
2645 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
2650 const struct matrix_function *f = matrix_parse_function_name (token);
2654 struct matrix_expr *e = xmalloc (sizeof *e);
2655 *e = (struct matrix_expr) {
2657 .location = lex_get_location (s->lexer, 0, 0)
2660 lex_get_n (s->lexer, 2);
2661 if (lex_token (s->lexer) != T_RPAREN)
2663 size_t allocated_subs = 0;
2666 struct msg_location *arg_location = lex_get_location (s->lexer, 0, 0);
2667 struct matrix_expr *sub = matrix_parse_expr (s);
2670 msg_location_destroy (arg_location);
2675 lex_extend_location (s->lexer, 0, arg_location);
2676 sub->location = arg_location;
2679 msg_location_destroy (arg_location);
2681 if (e->n_subs >= allocated_subs)
2682 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
2683 e->subs[e->n_subs++] = sub;
2685 while (lex_match (s->lexer, T_COMMA));
2687 if (!lex_force_match (s->lexer, T_RPAREN))
2690 if (e->n_subs < f->min_args || e->n_subs > f->max_args)
2692 if (f->min_args == f->max_args)
2693 msg (SE, ngettext ("Matrix function %s requires %zu argument.",
2694 "Matrix function %s requires %zu arguments.",
2696 f->name, f->min_args);
2697 else if (f->min_args == 1 && f->max_args == 2)
2698 msg (SE, ngettext ("Matrix function %s requires 1 or 2 arguments, "
2699 "but %zu was provided.",
2700 "Matrix function %s requires 1 or 2 arguments, "
2701 "but %zu were provided.",
2703 f->name, e->n_subs);
2704 else if (f->min_args == 1 && f->max_args == INT_MAX)
2705 msg (SE, _("Matrix function %s requires at least one argument."),
2717 matrix_expr_destroy (e);
2721 static struct matrix_expr *
2722 matrix_parse_primary (struct matrix_state *s)
2724 if (lex_is_number (s->lexer))
2726 double number = lex_number (s->lexer);
2728 return matrix_expr_create_number (number);
2730 else if (lex_is_string (s->lexer))
2732 char string[sizeof (double)];
2733 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
2737 memcpy (&number, string, sizeof number);
2738 return matrix_expr_create_number (number);
2740 else if (lex_match (s->lexer, T_LPAREN))
2742 struct matrix_expr *e = matrix_parse_or_xor (s);
2743 if (!e || !lex_force_match (s->lexer, T_RPAREN))
2745 matrix_expr_destroy (e);
2750 else if (lex_match (s->lexer, T_LCURLY))
2752 struct matrix_expr *e = matrix_parse_curly_semi (s);
2753 if (!e || !lex_force_match (s->lexer, T_RCURLY))
2755 matrix_expr_destroy (e);
2760 else if (lex_token (s->lexer) == T_ID)
2762 struct matrix_expr *retval;
2763 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
2766 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
2769 lex_error (s->lexer, _("Unknown variable %s."),
2770 lex_tokcstr (s->lexer));
2775 struct matrix_expr *e = xmalloc (sizeof *e);
2776 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
2779 else if (lex_token (s->lexer) == T_ALL)
2781 struct matrix_expr *retval;
2782 if (matrix_parse_function (s, "ALL", &retval))
2786 lex_error (s->lexer, NULL);
2790 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
2793 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
2795 if (lex_match (s->lexer, T_COLON))
2802 *indexp = matrix_parse_expr (s);
2803 return *indexp != NULL;
2807 static struct matrix_expr *
2808 matrix_parse_postfix (struct matrix_state *s)
2810 struct matrix_expr *lhs = matrix_parse_primary (s);
2811 if (!lhs || !lex_match (s->lexer, T_LPAREN))
2814 struct matrix_expr *i0;
2815 if (!matrix_parse_index_expr (s, &i0))
2817 matrix_expr_destroy (lhs);
2820 if (lex_match (s->lexer, T_RPAREN))
2822 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
2823 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
2824 else if (lex_match (s->lexer, T_COMMA))
2826 struct matrix_expr *i1;
2827 if (!matrix_parse_index_expr (s, &i1)
2828 || !lex_force_match (s->lexer, T_RPAREN))
2830 matrix_expr_destroy (lhs);
2831 matrix_expr_destroy (i0);
2832 matrix_expr_destroy (i1);
2835 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
2836 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
2837 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
2842 lex_error_expecting (s->lexer, "`)'", "`,'");
2847 static struct matrix_expr *
2848 matrix_parse_unary (struct matrix_state *s)
2850 if (lex_match (s->lexer, T_DASH))
2852 struct matrix_expr *lhs = matrix_parse_unary (s);
2853 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
2855 else if (lex_match (s->lexer, T_PLUS))
2856 return matrix_parse_unary (s);
2858 return matrix_parse_postfix (s);
2861 static struct matrix_expr *
2862 matrix_parse_seq (struct matrix_state *s)
2864 struct matrix_expr *start = matrix_parse_unary (s);
2865 if (!start || !lex_match (s->lexer, T_COLON))
2868 struct matrix_expr *end = matrix_parse_unary (s);
2871 matrix_expr_destroy (start);
2875 if (lex_match (s->lexer, T_COLON))
2877 struct matrix_expr *increment = matrix_parse_unary (s);
2880 matrix_expr_destroy (start);
2881 matrix_expr_destroy (end);
2884 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2887 return matrix_expr_create_2 (MOP_SEQ, start, end);
2890 static struct matrix_expr *
2891 matrix_parse_exp (struct matrix_state *s)
2893 struct matrix_expr *lhs = matrix_parse_seq (s);
2900 if (lex_match (s->lexer, T_EXP))
2902 else if (lex_match_phrase (s->lexer, "&**"))
2907 struct matrix_expr *rhs = matrix_parse_seq (s);
2910 matrix_expr_destroy (lhs);
2913 lhs = matrix_expr_create_2 (op, lhs, rhs);
2917 static struct matrix_expr *
2918 matrix_parse_mul_div (struct matrix_state *s)
2920 struct matrix_expr *lhs = matrix_parse_exp (s);
2927 if (lex_match (s->lexer, T_ASTERISK))
2929 else if (lex_match (s->lexer, T_SLASH))
2931 else if (lex_match_phrase (s->lexer, "&*"))
2933 else if (lex_match_phrase (s->lexer, "&/"))
2938 struct matrix_expr *rhs = matrix_parse_exp (s);
2941 matrix_expr_destroy (lhs);
2944 lhs = matrix_expr_create_2 (op, lhs, rhs);
2948 static struct matrix_expr *
2949 matrix_parse_add_sub (struct matrix_state *s)
2951 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2958 if (lex_match (s->lexer, T_PLUS))
2960 else if (lex_match (s->lexer, T_DASH))
2962 else if (lex_token (s->lexer) == T_NEG_NUM)
2967 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2970 matrix_expr_destroy (lhs);
2973 lhs = matrix_expr_create_2 (op, lhs, rhs);
2977 static struct matrix_expr *
2978 matrix_parse_relational (struct matrix_state *s)
2980 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2987 if (lex_match (s->lexer, T_GT))
2989 else if (lex_match (s->lexer, T_GE))
2991 else if (lex_match (s->lexer, T_LT))
2993 else if (lex_match (s->lexer, T_LE))
2995 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2997 else if (lex_match (s->lexer, T_NE))
3002 struct matrix_expr *rhs = matrix_parse_add_sub (s);
3005 matrix_expr_destroy (lhs);
3008 lhs = matrix_expr_create_2 (op, lhs, rhs);
3012 static struct matrix_expr *
3013 matrix_parse_not (struct matrix_state *s)
3015 if (lex_match (s->lexer, T_NOT))
3017 struct matrix_expr *lhs = matrix_parse_not (s);
3018 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
3021 return matrix_parse_relational (s);
3024 static struct matrix_expr *
3025 matrix_parse_and (struct matrix_state *s)
3027 struct matrix_expr *lhs = matrix_parse_not (s);
3031 while (lex_match (s->lexer, T_AND))
3033 struct matrix_expr *rhs = matrix_parse_not (s);
3036 matrix_expr_destroy (lhs);
3039 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
3044 static struct matrix_expr *
3045 matrix_parse_or_xor (struct matrix_state *s)
3047 struct matrix_expr *lhs = matrix_parse_and (s);
3054 if (lex_match (s->lexer, T_OR))
3056 else if (lex_match_id (s->lexer, "XOR"))
3061 struct matrix_expr *rhs = matrix_parse_and (s);
3064 matrix_expr_destroy (lhs);
3067 lhs = matrix_expr_create_2 (op, lhs, rhs);
3071 static struct matrix_expr *
3072 matrix_parse_expr (struct matrix_state *s)
3074 return matrix_parse_or_xor (s);
3077 /* Expression evaluation. */
3080 matrix_op_eval (enum matrix_op op, double a, double b)
3084 case MOP_ADD_ELEMS: return a + b;
3085 case MOP_SUB_ELEMS: return a - b;
3086 case MOP_MUL_ELEMS: return a * b;
3087 case MOP_DIV_ELEMS: return a / b;
3088 case MOP_EXP_ELEMS: return pow (a, b);
3089 case MOP_GT: return a > b;
3090 case MOP_GE: return a >= b;
3091 case MOP_LT: return a < b;
3092 case MOP_LE: return a <= b;
3093 case MOP_EQ: return a == b;
3094 case MOP_NE: return a != b;
3095 case MOP_AND: return (a > 0) && (b > 0);
3096 case MOP_OR: return (a > 0) || (b > 0);
3097 case MOP_XOR: return (a > 0) != (b > 0);
3099 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3108 case MOP_PASTE_HORZ:
3109 case MOP_PASTE_VERT:
3125 matrix_op_name (enum matrix_op op)
3129 case MOP_ADD_ELEMS: return "+";
3130 case MOP_SUB_ELEMS: return "-";
3131 case MOP_MUL_ELEMS: return "&*";
3132 case MOP_DIV_ELEMS: return "&/";
3133 case MOP_EXP_ELEMS: return "&**";
3134 case MOP_GT: return ">";
3135 case MOP_GE: return ">=";
3136 case MOP_LT: return "<";
3137 case MOP_LE: return "<=";
3138 case MOP_EQ: return "=";
3139 case MOP_NE: return "<>";
3140 case MOP_AND: return "AND";
3141 case MOP_OR: return "OR";
3142 case MOP_XOR: return "XOR";
3144 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3153 case MOP_PASTE_HORZ:
3154 case MOP_PASTE_VERT:
3170 is_scalar (const gsl_matrix *m)
3172 return m->size1 == 1 && m->size2 == 1;
3176 to_scalar (const gsl_matrix *m)
3178 assert (is_scalar (m));
3179 return gsl_matrix_get (m, 0, 0);
3183 matrix_expr_evaluate_elementwise (enum matrix_op op,
3184 gsl_matrix *a, gsl_matrix *b)
3188 double be = to_scalar (b);
3189 for (size_t r = 0; r < a->size1; r++)
3190 for (size_t c = 0; c < a->size2; c++)
3192 double *ae = gsl_matrix_ptr (a, r, c);
3193 *ae = matrix_op_eval (op, *ae, be);
3197 else if (is_scalar (a))
3199 double ae = to_scalar (a);
3200 for (size_t r = 0; r < b->size1; r++)
3201 for (size_t c = 0; c < b->size2; c++)
3203 double *be = gsl_matrix_ptr (b, r, c);
3204 *be = matrix_op_eval (op, ae, *be);
3208 else if (a->size1 == b->size1 && a->size2 == b->size2)
3210 for (size_t r = 0; r < a->size1; r++)
3211 for (size_t c = 0; c < a->size2; c++)
3213 double *ae = gsl_matrix_ptr (a, r, c);
3214 double be = gsl_matrix_get (b, r, c);
3215 *ae = matrix_op_eval (op, *ae, be);
3221 msg (SE, _("Operands to %s must have the same dimensions or one "
3222 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
3223 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
3229 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
3231 if (is_scalar (a) || is_scalar (b))
3232 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
3234 if (a->size2 != b->size1)
3236 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
3237 "not conformable for multiplication."),
3238 a->size1, a->size2, b->size1, b->size2);
3242 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3243 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3248 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3250 gsl_matrix *tmp = *a;
3256 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3259 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3260 swap_matrix (z, tmp);
3264 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3266 mul_matrix (x, *x, *x, tmp);
3270 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
3273 if (x->size1 != x->size2)
3275 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
3276 "the left-hand size, not one with dimensions %zu×%zu."),
3277 x->size1, x->size2);
3282 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
3283 "right-hand side, not a matrix with dimensions %zu×%zu."),
3284 b->size1, b->size2);
3287 double bf = to_scalar (b);
3288 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3290 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
3291 "or outside the valid range."), bf);
3296 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3298 gsl_matrix_set_identity (y);
3302 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3304 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3307 mul_matrix (&y, x, y, &t);
3308 square_matrix (&x, &t);
3311 square_matrix (&x, &t);
3313 mul_matrix (&y, x, y, &t);
3317 /* Garbage collection.
3319 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3320 a permutation of them. We are returning one of them; that one must not be
3321 destroyed. We must not destroy 'x_' because the caller owns it. */
3323 gsl_matrix_free (y_);
3325 gsl_matrix_free (t_);
3331 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
3334 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3336 msg (SE, _("All operands of : operator must be scalars."));
3340 long int start = to_scalar (start_);
3341 long int end = to_scalar (end_);
3342 long int by = by_ ? to_scalar (by_) : 1;
3346 msg (SE, _("The increment operand to : must be nonzero."));
3350 long int n = (end >= start && by > 0 ? (end - start + by) / by
3351 : end <= start && by < 0 ? (start - end - by) / -by
3353 gsl_matrix *m = gsl_matrix_alloc (1, n);
3354 for (long int i = 0; i < n; i++)
3355 gsl_matrix_set (m, 0, i, start + i * by);
3360 matrix_expr_evaluate_not (gsl_matrix *a)
3362 for (size_t r = 0; r < a->size1; r++)
3363 for (size_t c = 0; c < a->size2; c++)
3365 double *ae = gsl_matrix_ptr (a, r, c);
3372 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
3374 if (a->size1 != b->size1)
3376 if (!a->size1 || !a->size2)
3378 else if (!b->size1 || !b->size2)
3381 msg (SE, _("All columns in a matrix must have the same number of rows, "
3382 "but this tries to paste matrices with %zu and %zu rows."),
3383 a->size1, b->size1);
3387 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3388 for (size_t y = 0; y < a->size1; y++)
3390 for (size_t x = 0; x < a->size2; x++)
3391 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3392 for (size_t x = 0; x < b->size2; x++)
3393 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3399 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
3401 if (a->size2 != b->size2)
3403 if (!a->size1 || !a->size2)
3405 else if (!b->size1 || !b->size2)
3408 msg (SE, _("All rows in a matrix must have the same number of columns, "
3409 "but this tries to stack matrices with %zu and %zu columns."),
3410 a->size2, b->size2);
3414 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3415 for (size_t x = 0; x < a->size2; x++)
3417 for (size_t y = 0; y < a->size1; y++)
3418 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3419 for (size_t y = 0; y < b->size1; y++)
3420 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3426 matrix_to_vector (gsl_matrix *m)
3429 gsl_vector v = to_vector (m);
3430 assert (v.block == m->block || !v.block);
3434 gsl_matrix_free (m);
3435 return xmemdup (&v, sizeof v);
3449 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3452 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
3453 enum index_type index_type, size_t other_size,
3454 struct index_vector *iv)
3463 msg (SE, _("Vector index must be scalar or vector, not a "
3465 m->size1, m->size2);
3469 msg (SE, _("Matrix row index must be scalar or vector, not a "
3471 m->size1, m->size2);
3475 msg (SE, _("Matrix column index must be scalar or vector, not a "
3477 m->size1, m->size2);
3483 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3484 *iv = (struct index_vector) {
3485 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3488 for (size_t i = 0; i < v.size; i++)
3490 double index = gsl_vector_get (&v, i);
3491 if (index < 1 || index >= size + 1)
3496 msg (SE, _("Index %g is out of range for vector "
3497 "with %zu elements."), index, size);
3501 msg (SE, _("%g is not a valid row index for "
3502 "a %zu×%zu matrix."),
3503 index, size, other_size);
3507 msg (SE, _("%g is not a valid column index for "
3508 "a %zu×%zu matrix."),
3509 index, other_size, size);
3516 iv->indexes[i] = index - 1;
3522 *iv = (struct index_vector) {
3523 .indexes = xnmalloc (size, sizeof *iv->indexes),
3526 for (size_t i = 0; i < size; i++)
3533 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
3535 if (!is_vector (sm))
3537 msg (SE, _("Vector index operator may not be applied to "
3538 "a %zu×%zu matrix."),
3539 sm->size1, sm->size2);
3547 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
3549 if (!matrix_expr_evaluate_vec_all (sm))
3552 gsl_vector sv = to_vector (sm);
3553 struct index_vector iv;
3554 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
3557 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3558 sm->size1 == 1 ? iv.n : 1);
3559 gsl_vector dv = to_vector (dm);
3560 for (size_t dx = 0; dx < iv.n; dx++)
3562 size_t sx = iv.indexes[dx];
3563 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3571 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
3574 struct index_vector iv0;
3575 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
3578 struct index_vector iv1;
3579 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
3586 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3587 for (size_t dy = 0; dy < iv0.n; dy++)
3589 size_t sy = iv0.indexes[dy];
3591 for (size_t dx = 0; dx < iv1.n; dx++)
3593 size_t sx = iv1.indexes[dx];
3594 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3602 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3603 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3604 const struct matrix_function_properties *, gsl_matrix *[], size_t, \
3605 matrix_proto_##PROTO *);
3610 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
3612 if (!is_scalar (subs[index]))
3614 msg (SE, _("Function %s argument %zu must be a scalar, "
3615 "not a %zu×%zu matrix."),
3616 name, index + 1, subs[index]->size1, subs[index]->size2);
3623 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
3625 if (!is_vector (subs[index]))
3627 msg (SE, _("Function %s argument %zu must be a vector, "
3628 "not a %zu×%zu matrix."),
3629 name, index + 1, subs[index]->size1, subs[index]->size2);
3636 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
3638 for (size_t i = 0; i < n_subs; i++)
3640 if (!check_scalar_arg (name, subs, i))
3642 d[i] = to_scalar (subs[i]);
3648 parse_constraint_value (const char **constraintsp)
3651 long retval = strtol (*constraintsp, &tail, 10);
3652 assert (tail > *constraintsp);
3653 *constraintsp = tail;
3658 argument_invalid (const struct matrix_function_properties *props,
3659 size_t arg_index, gsl_matrix *a, size_t y, size_t x,
3663 ds_put_format (out, _("Argument %zu to matrix function %s "
3664 "has invalid value %g."),
3665 arg_index, props->name, gsl_matrix_get (a, y, x));
3667 ds_put_format (out, _("Row %zu, column %zu of argument %zu "
3668 "to matrix function %s has "
3669 "invalid value %g."),
3670 y + 1, x + 1, arg_index, props->name,
3671 gsl_matrix_get (a, y, x));
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 bool greater, bool equal)
3680 struct string s = DS_EMPTY_INITIALIZER;
3682 argument_invalid (props, a_index, a, y, x, &s);
3683 ds_put_cstr (&s, " ");
3684 if (greater && equal)
3685 ds_put_format (&s, _("This argument must be greater than or "
3686 "equal to argument %zu, but that has value %g."),
3688 else if (greater && !equal)
3689 ds_put_format (&s, _("This argument must be greater than argument %zu, "
3690 "but that has value %g."),
3693 ds_put_format (&s, _("This argument must be less than or "
3694 "equal to argument %zu, but that has value %g."),
3697 ds_put_format (&s, _("This argument must be less than argument %zu, "
3698 "but that has value %g."),
3700 msg (ME, ds_cstr (&s));
3705 argument_value_error (const struct matrix_function_properties *props,
3706 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3708 bool greater, bool equal)
3710 struct string s = DS_EMPTY_INITIALIZER;
3711 argument_invalid (props, a_index, a, y, x, &s);
3712 ds_put_cstr (&s, " ");
3713 if (greater && equal)
3714 ds_put_format (&s, _("This argument must be greater than or equal to %g."),
3716 else if (greater && !equal)
3717 ds_put_format (&s, _("This argument must be greater than %g."), b);
3719 ds_put_format (&s, _("This argument must be less than or equal to %g."), b);
3721 ds_put_format (&s, _("This argument must be less than %g."), b);
3722 msg (ME, ds_cstr (&s));
3727 check_constraints (const struct matrix_function_properties *props,
3728 gsl_matrix *args[], size_t n_args)
3730 const char *constraints = props->constraints;
3734 size_t arg_index = SIZE_MAX;
3735 while (*constraints)
3737 if (*constraints >= 'a' && *constraints <= 'd')
3739 arg_index = *constraints++ - 'a';
3740 assert (arg_index < n_args);
3742 else if (*constraints == '[' || *constraints == '(')
3744 assert (arg_index < n_args);
3745 bool open_lower = *constraints++ == '(';
3746 int minimum = parse_constraint_value (&constraints);
3747 assert (*constraints == ',');
3749 int maximum = parse_constraint_value (&constraints);
3750 assert (*constraints == ']' || *constraints == ')');
3751 bool open_upper = *constraints++ == ')';
3753 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3754 if ((open_lower ? *d <= minimum : *d < minimum)
3755 || (open_upper ? *d >= maximum : *d > maximum))
3757 if (!is_scalar (args[arg_index]))
3758 msg (ME, _("Row %zu, column %zu of argument %zu to matrix "
3759 "function %s has value %g, which is outside "
3760 "the valid range %c%d,%d%c."),
3761 y + 1, x + 1, arg_index + 1, props->name, *d,
3762 open_lower ? '(' : '[',
3764 open_upper ? ')' : ']');
3766 msg (ME, _("Argument %zu to matrix function %s has value %g, "
3767 "which is outside the valid range %c%d,%d%c."),
3768 arg_index + 1, props->name, *d,
3769 open_lower ? '(' : '[',
3771 open_upper ? ')' : ']');
3775 else if (*constraints == '>' || *constraints == '<')
3777 bool greater = *constraints++ == '>';
3779 if (*constraints == '=')
3787 if (*constraints >= 'a' && *constraints <= 'd')
3789 size_t a_index = arg_index;
3790 size_t b_index = *constraints - 'a';
3791 assert (a_index < n_args);
3792 assert (b_index < n_args);
3794 /* We only support one of the two arguments being non-scalar.
3795 It's easier to support only the first one being non-scalar, so
3796 flip things around if it's the other way. */
3797 if (!is_scalar (args[b_index]))
3799 assert (is_scalar (args[a_index]));
3800 size_t tmp_index = a_index;
3802 b_index = tmp_index;
3807 double b = to_scalar (args[b_index]);
3808 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
3810 ? (equal ? !(*a >= b) : !(*a > b))
3811 : (equal ? !(*a <= b) : !(*a < b)))
3813 argument_inequality_error (props,
3814 a_index + 1, args[a_index], y, x,
3822 int comparand = parse_constraint_value (&constraints);
3824 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3826 ? (equal ? !(*d >= comparand) : !(*d > comparand))
3827 : (equal ? !(*d <= comparand) : !(*d < comparand)))
3829 argument_value_error (props,
3830 arg_index + 1, args[arg_index], y, x,
3839 assert (*constraints == ' ');
3841 arg_index = SIZE_MAX;
3848 matrix_expr_evaluate_d_none (
3849 const struct matrix_function_properties *props UNUSED,
3850 gsl_matrix *subs[] UNUSED, size_t n_subs,
3851 matrix_proto_d_none *f)
3853 assert (n_subs == 0);
3855 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3856 gsl_matrix_set (m, 0, 0, f ());
3861 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
3862 gsl_matrix *subs[], size_t n_subs,
3863 matrix_proto_d_d *f)
3865 assert (n_subs == 1);
3868 if (!to_scalar_args (props->name, subs, n_subs, &d))
3871 if (!check_constraints (props, subs, n_subs))
3874 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3875 gsl_matrix_set (m, 0, 0, f (d));
3880 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
3881 gsl_matrix *subs[], size_t n_subs,
3882 matrix_proto_d_dd *f)
3884 assert (n_subs == 2);
3887 if (!to_scalar_args (props->name, subs, n_subs, d))
3890 if (!check_constraints (props, subs, n_subs))
3893 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3894 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
3899 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
3900 gsl_matrix *subs[], size_t n_subs,
3901 matrix_proto_d_ddd *f)
3903 assert (n_subs == 3);
3906 if (!to_scalar_args (props->name, subs, n_subs, d))
3909 if (!check_constraints (props, subs, n_subs))
3912 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3913 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
3918 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
3919 gsl_matrix *subs[], size_t n_subs,
3920 matrix_proto_m_d *f)
3922 assert (n_subs == 1);
3925 return to_scalar_args (props->name, subs, n_subs, &d) ? f (d) : NULL;
3929 matrix_expr_evaluate_m_dd (const struct matrix_function_properties *props,
3930 gsl_matrix *subs[], size_t n_subs,
3931 matrix_proto_m_dd *f)
3933 assert (n_subs == 2);
3936 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
3940 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
3941 gsl_matrix *subs[], size_t n_subs,
3942 matrix_proto_m_ddd *f)
3944 assert (n_subs == 3);
3947 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
3951 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED,
3952 gsl_matrix *subs[], size_t n_subs,
3953 matrix_proto_m_m *f)
3955 assert (n_subs == 1);
3960 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
3961 gsl_matrix *subs[], size_t n_subs,
3962 matrix_proto_m_e *f)
3964 assert (n_subs == 1);
3966 if (!check_constraints (props, subs, n_subs))
3969 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3975 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED,
3976 gsl_matrix *subs[], size_t n_subs,
3977 matrix_proto_m_md *f)
3979 assert (n_subs == 2);
3980 if (!check_scalar_arg (props->name, subs, 1))
3982 return f (subs[0], to_scalar (subs[1]));
3986 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
3987 gsl_matrix *subs[], size_t n_subs,
3988 matrix_proto_m_ed *f)
3990 assert (n_subs == 2);
3991 if (!check_scalar_arg (props->name, subs, 1))
3994 if (!check_constraints (props, subs, n_subs))
3997 double b = to_scalar (subs[1]);
3998 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4004 matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props UNUSED,
4005 gsl_matrix *subs[], size_t n_subs,
4006 matrix_proto_m_mdd *f)
4008 assert (n_subs == 3);
4009 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
4011 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
4015 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
4016 gsl_matrix *subs[], size_t n_subs,
4017 matrix_proto_m_edd *f)
4019 assert (n_subs == 3);
4020 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
4023 if (!check_constraints (props, subs, n_subs))
4026 double b = to_scalar (subs[1]);
4027 double c = to_scalar (subs[2]);
4028 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4034 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4035 gsl_matrix *subs[], size_t n_subs,
4036 matrix_proto_m_eddd *f)
4038 assert (n_subs == 4);
4039 for (size_t i = 1; i < 4; i++)
4040 if (!check_scalar_arg (props->name, subs, i))
4043 if (!check_constraints (props, subs, n_subs))
4046 double b = to_scalar (subs[1]);
4047 double c = to_scalar (subs[2]);
4048 double d = to_scalar (subs[3]);
4049 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4050 *a = f (*a, b, c, d);
4055 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4056 gsl_matrix *subs[], size_t n_subs,
4057 matrix_proto_m_eed *f)
4059 assert (n_subs == 3);
4060 if (!check_scalar_arg (props->name, subs, 2))
4063 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4064 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4066 msg (ME, _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4067 "%zu×%zu, but %s requires these arguments either to have "
4068 "the same dimensions or for one of them to be a scalar."),
4070 subs[0]->size1, subs[0]->size2,
4071 subs[1]->size1, subs[1]->size2,
4076 if (!check_constraints (props, subs, n_subs))
4079 double c = to_scalar (subs[2]);
4081 if (is_scalar (subs[0]))
4083 double a = to_scalar (subs[0]);
4084 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4090 double b = to_scalar (subs[1]);
4091 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4098 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED,
4099 gsl_matrix *subs[], size_t n_subs,
4100 matrix_proto_m_mm *f)
4102 assert (n_subs == 2);
4103 return f (subs[0], subs[1]);
4107 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4108 gsl_matrix *subs[], size_t n_subs,
4109 matrix_proto_m_v *f)
4111 assert (n_subs == 1);
4112 if (!check_vector_arg (props->name, subs, 0))
4114 gsl_vector v = to_vector (subs[0]);
4119 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED,
4120 gsl_matrix *subs[], size_t n_subs,
4121 matrix_proto_d_m *f)
4123 assert (n_subs == 1);
4125 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4126 gsl_matrix_set (m, 0, 0, f (subs[0]));
4131 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED,
4132 gsl_matrix *subs[], size_t n_subs,
4133 matrix_proto_m_any *f)
4135 return f (subs, n_subs);
4139 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props,
4140 gsl_matrix *subs[], size_t n_subs,
4141 matrix_proto_IDENT *f)
4143 assert (n_subs <= 2);
4146 if (!to_scalar_args (props->name, subs, n_subs, d))
4149 return f (d[0], d[n_subs - 1]);
4153 matrix_expr_evaluate (const struct matrix_expr *e)
4155 if (e->op == MOP_NUMBER)
4157 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4158 gsl_matrix_set (m, 0, 0, e->number);
4161 else if (e->op == MOP_VARIABLE)
4163 const gsl_matrix *src = e->variable->value;
4166 msg (SE, _("Uninitialized variable %s used in expression."),
4171 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4172 gsl_matrix_memcpy (dst, src);
4175 else if (e->op == MOP_EOF)
4177 struct dfm_reader *reader = read_file_open (e->eof);
4178 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4179 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4183 enum { N_LOCAL = 3 };
4184 gsl_matrix *local_subs[N_LOCAL];
4185 gsl_matrix **subs = (e->n_subs < N_LOCAL
4187 : xmalloc (e->n_subs * sizeof *subs));
4189 for (size_t i = 0; i < e->n_subs; i++)
4191 subs[i] = matrix_expr_evaluate (e->subs[i]);
4194 for (size_t j = 0; j < i; j++)
4195 gsl_matrix_free (subs[j]);
4196 if (subs != local_subs)
4202 gsl_matrix *result = NULL;
4205 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4206 case MOP_F_##ENUM: \
4208 static const struct matrix_function_properties props = { \
4210 .constraints = CONSTRAINTS, \
4212 result = matrix_expr_evaluate_##PROTO (&props, subs, e->n_subs, \
4213 matrix_eval_##ENUM); \
4220 gsl_matrix_scale (subs[0], -1.0);
4238 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
4242 result = matrix_expr_evaluate_not (subs[0]);
4246 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
4250 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
4254 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
4258 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
4261 case MOP_PASTE_HORZ:
4262 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
4265 case MOP_PASTE_VERT:
4266 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
4270 result = gsl_matrix_alloc (0, 0);
4274 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
4278 result = matrix_expr_evaluate_vec_all (subs[0]);
4282 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
4286 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
4290 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
4299 for (size_t i = 0; i < e->n_subs; i++)
4300 if (subs[i] != result)
4301 gsl_matrix_free (subs[i]);
4302 if (subs != local_subs)
4308 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4311 gsl_matrix *m = matrix_expr_evaluate (e);
4317 msg (SE, _("Expression for %s must evaluate to scalar, "
4318 "not a %zu×%zu matrix."),
4319 context, m->size1, m->size2);
4320 gsl_matrix_free (m);
4325 gsl_matrix_free (m);
4330 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4334 if (!matrix_expr_evaluate_scalar (e, context, &d))
4338 if (d < LONG_MIN || d > LONG_MAX)
4340 msg (SE, _("Expression for %s is outside the integer range."), context);
4347 struct matrix_lvalue
4349 struct matrix_var *var;
4350 struct msg_location *var_location;
4351 struct matrix_expr *indexes[2];
4356 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4360 for (size_t i = 0; i < lvalue->n_indexes; i++)
4361 matrix_expr_destroy (lvalue->indexes[i]);
4366 static struct matrix_lvalue *
4367 matrix_lvalue_parse (struct matrix_state *s)
4369 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4370 if (!lex_force_id (s->lexer))
4372 lvalue->var_location = lex_get_location (s->lexer, 0, 0);
4373 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4374 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4378 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
4383 lex_get_n (s->lexer, 2);
4385 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
4387 lvalue->n_indexes++;
4389 if (lex_match (s->lexer, T_COMMA))
4391 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
4393 lvalue->n_indexes++;
4395 if (!lex_force_match (s->lexer, T_RPAREN))
4401 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4407 matrix_lvalue_destroy (lvalue);
4412 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4413 enum index_type index_type, size_t other_size,
4414 struct index_vector *iv)
4419 m = matrix_expr_evaluate (e);
4426 bool ok = matrix_normalize_index_vector (m, size, index_type,
4428 gsl_matrix_free (m);
4433 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4434 struct index_vector *iv, gsl_matrix *sm)
4436 gsl_vector dv = to_vector (lvalue->var->value);
4438 /* Convert source matrix 'sm' to source vector 'sv'. */
4439 if (!is_vector (sm))
4441 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
4442 sm->size1, sm->size2);
4445 gsl_vector sv = to_vector (sm);
4446 if (iv->n != sv.size)
4448 msg (SE, _("Can't assign %zu-element vector "
4449 "to %zu-element subvector."), sv.size, iv->n);
4453 /* Assign elements. */
4454 for (size_t x = 0; x < iv->n; x++)
4455 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4460 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4461 struct index_vector *iv0,
4462 struct index_vector *iv1, gsl_matrix *sm)
4464 gsl_matrix *dm = lvalue->var->value;
4466 /* Convert source matrix 'sm' to source vector 'sv'. */
4467 if (iv0->n != sm->size1)
4469 msg (SE, _("Row index vector for assignment to %s has %zu elements "
4470 "but source matrix has %zu rows."),
4471 lvalue->var->name, iv0->n, sm->size1);
4474 if (iv1->n != sm->size2)
4476 msg (SE, _("Column index vector for assignment to %s has %zu elements "
4477 "but source matrix has %zu columns."),
4478 lvalue->var->name, iv1->n, sm->size2);
4482 /* Assign elements. */
4483 for (size_t y = 0; y < iv0->n; y++)
4485 size_t f0 = iv0->indexes[y];
4486 for (size_t x = 0; x < iv1->n; x++)
4488 size_t f1 = iv1->indexes[x];
4489 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4495 /* Takes ownership of M. */
4497 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4498 struct index_vector *iv0, struct index_vector *iv1,
4501 if (!lvalue->n_indexes)
4503 gsl_matrix_free (lvalue->var->value);
4504 lvalue->var->value = sm;
4509 bool ok = (lvalue->n_indexes == 1
4510 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
4511 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
4512 free (iv0->indexes);
4513 free (iv1->indexes);
4514 gsl_matrix_free (sm);
4520 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4521 struct index_vector *iv0,
4522 struct index_vector *iv1)
4524 *iv0 = INDEX_VECTOR_INIT;
4525 *iv1 = INDEX_VECTOR_INIT;
4526 if (!lvalue->n_indexes)
4529 /* Validate destination matrix exists and has the right shape. */
4530 gsl_matrix *dm = lvalue->var->value;
4533 msg_at (SE, lvalue->var_location,
4534 _("Undefined variable %s."), lvalue->var->name);
4537 else if (dm->size1 == 0 || dm->size2 == 0)
4539 msg_at (SE, lvalue->var_location,
4540 _("Cannot index %zu×%zu matrix."), dm->size1, dm->size2);
4543 else if (lvalue->n_indexes == 1)
4545 if (!is_vector (dm))
4547 msg_at (SE, lvalue->var_location,
4548 _("Can't use vector indexing on %zu×%zu matrix %s."),
4549 dm->size1, dm->size2, lvalue->var->name);
4552 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4553 MAX (dm->size1, dm->size2),
4558 assert (lvalue->n_indexes == 2);
4559 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4560 IV_ROW, dm->size2, iv0))
4563 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4564 IV_COLUMN, dm->size1, iv1))
4566 free (iv0->indexes);
4573 /* Takes ownership of M. */
4575 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
4577 struct index_vector iv0, iv1;
4578 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
4580 gsl_matrix_free (sm);
4584 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
4588 /* Matrix command. */
4592 struct matrix_cmd **commands;
4596 static void matrix_cmds_uninit (struct matrix_cmds *);
4600 enum matrix_cmd_type
4623 struct compute_command
4625 struct matrix_lvalue *lvalue;
4626 struct matrix_expr *rvalue;
4630 struct print_command
4632 struct matrix_expr *expression;
4633 bool use_default_format;
4634 struct fmt_spec format;
4636 int space; /* -1 means NEWPAGE. */
4640 struct string_array literals; /* CLABELS/RLABELS. */
4641 struct matrix_expr *expr; /* CNAMES/RNAMES. */
4647 struct do_if_command
4651 struct matrix_expr *condition;
4652 struct matrix_cmds commands;
4662 struct matrix_var *var;
4663 struct matrix_expr *start;
4664 struct matrix_expr *end;
4665 struct matrix_expr *increment;
4667 /* Loop conditions. */
4668 struct matrix_expr *top_condition;
4669 struct matrix_expr *bottom_condition;
4672 struct matrix_cmds commands;
4676 struct display_command
4678 struct matrix_state *state;
4682 struct release_command
4684 struct matrix_var **vars;
4691 struct matrix_expr *expression;
4692 struct save_file *sf;
4698 struct read_file *rf;
4699 struct matrix_lvalue *dst;
4700 struct matrix_expr *size;
4702 enum fmt_type format;
4709 struct write_command
4711 struct write_file *wf;
4712 struct matrix_expr *expression;
4715 /* If this is nonnull, WRITE uses this format.
4717 If this is NULL, WRITE uses free-field format with as many
4718 digits of precision as needed. */
4719 struct fmt_spec *format;
4728 struct matrix_lvalue *dst;
4729 struct dataset *dataset;
4730 struct file_handle *file;
4732 struct var_syntax *vars;
4734 struct matrix_var *names;
4736 /* Treatment of missing values. */
4741 MGET_ERROR, /* Flag error on command. */
4742 MGET_ACCEPT, /* Accept without change, user-missing only. */
4743 MGET_OMIT, /* Drop this case. */
4744 MGET_RECODE /* Recode to 'substitute'. */
4747 double substitute; /* MGET_RECODE only. */
4753 struct msave_command
4755 struct msave_common *common;
4756 struct matrix_expr *expr;
4757 const char *rowtype;
4758 const struct matrix_expr *factors;
4759 const struct matrix_expr *splits;
4765 struct matrix_state *state;
4766 struct file_handle *file;
4768 struct stringi_set rowtypes;
4772 struct eigen_command
4774 struct matrix_expr *expr;
4775 struct matrix_var *evec;
4776 struct matrix_var *eval;
4780 struct setdiag_command
4782 struct matrix_var *dst;
4783 struct matrix_expr *expr;
4789 struct matrix_expr *expr;
4790 struct matrix_var *u;
4791 struct matrix_var *s;
4792 struct matrix_var *v;
4798 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
4799 static bool matrix_cmd_execute (struct matrix_cmd *);
4800 static void matrix_cmd_destroy (struct matrix_cmd *);
4803 static struct matrix_cmd *
4804 matrix_parse_compute (struct matrix_state *s)
4806 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4807 *cmd = (struct matrix_cmd) {
4808 .type = MCMD_COMPUTE,
4809 .compute = { .lvalue = NULL },
4812 struct compute_command *compute = &cmd->compute;
4813 compute->lvalue = matrix_lvalue_parse (s);
4814 if (!compute->lvalue)
4817 if (!lex_force_match (s->lexer, T_EQUALS))
4820 compute->rvalue = matrix_parse_expr (s);
4821 if (!compute->rvalue)
4827 matrix_cmd_destroy (cmd);
4832 matrix_cmd_execute_compute (struct compute_command *compute)
4834 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
4838 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
4842 print_labels_destroy (struct print_labels *labels)
4846 string_array_destroy (&labels->literals);
4847 matrix_expr_destroy (labels->expr);
4853 parse_literal_print_labels (struct matrix_state *s,
4854 struct print_labels **labelsp)
4856 lex_match (s->lexer, T_EQUALS);
4857 print_labels_destroy (*labelsp);
4858 *labelsp = xzalloc (sizeof **labelsp);
4859 while (lex_token (s->lexer) != T_SLASH
4860 && lex_token (s->lexer) != T_ENDCMD
4861 && lex_token (s->lexer) != T_STOP)
4863 struct string label = DS_EMPTY_INITIALIZER;
4864 while (lex_token (s->lexer) != T_COMMA
4865 && lex_token (s->lexer) != T_SLASH
4866 && lex_token (s->lexer) != T_ENDCMD
4867 && lex_token (s->lexer) != T_STOP)
4869 if (!ds_is_empty (&label))
4870 ds_put_byte (&label, ' ');
4872 if (lex_token (s->lexer) == T_STRING)
4873 ds_put_cstr (&label, lex_tokcstr (s->lexer));
4876 char *rep = lex_next_representation (s->lexer, 0, 0);
4877 ds_put_cstr (&label, rep);
4882 string_array_append_nocopy (&(*labelsp)->literals,
4883 ds_steal_cstr (&label));
4885 if (!lex_match (s->lexer, T_COMMA))
4891 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
4893 lex_match (s->lexer, T_EQUALS);
4894 struct matrix_expr *e = matrix_parse_exp (s);
4898 print_labels_destroy (*labelsp);
4899 *labelsp = xzalloc (sizeof **labelsp);
4900 (*labelsp)->expr = e;
4904 static struct matrix_cmd *
4905 matrix_parse_print (struct matrix_state *s)
4907 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4908 *cmd = (struct matrix_cmd) {
4911 .use_default_format = true,
4915 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
4918 for (size_t i = 0; ; i++)
4920 enum token_type t = lex_next_token (s->lexer, i);
4921 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
4923 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
4925 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
4928 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
4933 cmd->print.expression = matrix_parse_exp (s);
4934 if (!cmd->print.expression)
4938 while (lex_match (s->lexer, T_SLASH))
4940 if (lex_match_id (s->lexer, "FORMAT"))
4942 lex_match (s->lexer, T_EQUALS);
4943 if (!parse_format_specifier (s->lexer, &cmd->print.format))
4945 cmd->print.use_default_format = false;
4947 else if (lex_match_id (s->lexer, "TITLE"))
4949 lex_match (s->lexer, T_EQUALS);
4950 if (!lex_force_string (s->lexer))
4952 free (cmd->print.title);
4953 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
4956 else if (lex_match_id (s->lexer, "SPACE"))
4958 lex_match (s->lexer, T_EQUALS);
4959 if (lex_match_id (s->lexer, "NEWPAGE"))
4960 cmd->print.space = -1;
4961 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
4963 cmd->print.space = lex_integer (s->lexer);
4969 else if (lex_match_id (s->lexer, "RLABELS"))
4970 parse_literal_print_labels (s, &cmd->print.rlabels);
4971 else if (lex_match_id (s->lexer, "CLABELS"))
4972 parse_literal_print_labels (s, &cmd->print.clabels);
4973 else if (lex_match_id (s->lexer, "RNAMES"))
4975 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
4978 else if (lex_match_id (s->lexer, "CNAMES"))
4980 if (!parse_expr_print_labels (s, &cmd->print.clabels))
4985 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
4986 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
4994 matrix_cmd_destroy (cmd);
4999 matrix_is_integer (const gsl_matrix *m)
5001 for (size_t y = 0; y < m->size1; y++)
5002 for (size_t x = 0; x < m->size2; x++)
5004 double d = gsl_matrix_get (m, y, x);
5012 matrix_max_magnitude (const gsl_matrix *m)
5015 for (size_t y = 0; y < m->size1; y++)
5016 for (size_t x = 0; x < m->size2; x++)
5018 double d = fabs (gsl_matrix_get (m, y, x));
5026 format_fits (struct fmt_spec format, double x)
5028 char *s = data_out (&(union value) { .f = x }, NULL,
5029 &format, settings_get_fmt_settings ());
5030 bool fits = *s != '*' && !strchr (s, 'E');
5035 static struct fmt_spec
5036 default_format (const gsl_matrix *m, int *log_scale)
5040 double max = matrix_max_magnitude (m);
5042 if (matrix_is_integer (m))
5043 for (int w = 1; w <= 12; w++)
5045 struct fmt_spec format = { .type = FMT_F, .w = w };
5046 if (format_fits (format, -max))
5050 if (max >= 1e9 || max <= 1e-4)
5053 snprintf (s, sizeof s, "%.1e", max);
5055 const char *e = strchr (s, 'e');
5057 *log_scale = atoi (e + 1);
5060 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5064 trimmed_string (double d)
5066 struct substring s = ss_buffer ((char *) &d, sizeof d);
5067 ss_rtrim (&s, ss_cstr (" "));
5068 return ss_xstrdup (s);
5071 static struct string_array *
5072 print_labels_get (const struct print_labels *labels, size_t n,
5073 const char *prefix, bool truncate)
5078 struct string_array *out = xzalloc (sizeof *out);
5079 if (labels->literals.n)
5080 string_array_clone (out, &labels->literals);
5081 else if (labels->expr)
5083 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5084 if (m && is_vector (m))
5086 gsl_vector v = to_vector (m);
5087 for (size_t i = 0; i < v.size; i++)
5088 string_array_append_nocopy (out, trimmed_string (
5089 gsl_vector_get (&v, i)));
5091 gsl_matrix_free (m);
5097 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5100 string_array_append (out, "");
5103 string_array_delete (out, out->n - 1);
5106 for (size_t i = 0; i < out->n; i++)
5108 char *s = out->strings[i];
5109 s[strnlen (s, 8)] = '\0';
5116 matrix_cmd_print_space (int space)
5119 output_item_submit (page_break_item_create ());
5120 for (int i = 0; i < space; i++)
5121 output_log ("%s", "");
5125 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
5126 struct fmt_spec format, int log_scale)
5128 matrix_cmd_print_space (print->space);
5130 output_log ("%s", print->title);
5132 output_log (" 10 ** %d X", log_scale);
5134 struct string_array *clabels = print_labels_get (print->clabels,
5135 m->size2, "col", true);
5136 if (clabels && format.w < 8)
5138 struct string_array *rlabels = print_labels_get (print->rlabels,
5139 m->size1, "row", true);
5143 struct string line = DS_EMPTY_INITIALIZER;
5145 ds_put_byte_multiple (&line, ' ', 8);
5146 for (size_t x = 0; x < m->size2; x++)
5147 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5148 output_log_nocopy (ds_steal_cstr (&line));
5151 double scale = pow (10.0, log_scale);
5152 bool numeric = fmt_is_numeric (format.type);
5153 for (size_t y = 0; y < m->size1; y++)
5155 struct string line = DS_EMPTY_INITIALIZER;
5157 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5159 for (size_t x = 0; x < m->size2; x++)
5161 double f = gsl_matrix_get (m, y, x);
5163 ? data_out (&(union value) { .f = f / scale}, NULL,
5164 &format, settings_get_fmt_settings ())
5165 : trimmed_string (f));
5166 ds_put_format (&line, " %s", s);
5169 output_log_nocopy (ds_steal_cstr (&line));
5172 string_array_destroy (rlabels);
5174 string_array_destroy (clabels);
5179 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5180 const struct print_labels *print_labels, size_t n,
5181 const char *name, const char *prefix)
5183 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5185 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5186 for (size_t i = 0; i < n; i++)
5187 pivot_category_create_leaf (
5189 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5190 : pivot_value_new_integer (i + 1)));
5192 d->hide_all_labels = true;
5193 string_array_destroy (labels);
5198 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
5199 struct fmt_spec format, int log_scale)
5201 struct pivot_table *table = pivot_table_create__ (
5202 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5205 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5207 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5208 N_("Columns"), "col");
5210 struct pivot_footnote *footnote = NULL;
5213 char *s = xasprintf ("× 10**%d", log_scale);
5214 footnote = pivot_table_create_footnote (
5215 table, pivot_value_new_user_text_nocopy (s));
5218 double scale = pow (10.0, log_scale);
5219 bool numeric = fmt_is_numeric (format.type);
5220 for (size_t y = 0; y < m->size1; y++)
5221 for (size_t x = 0; x < m->size2; x++)
5223 double f = gsl_matrix_get (m, y, x);
5224 struct pivot_value *v;
5227 v = pivot_value_new_number (f / scale);
5228 v->numeric.format = format;
5231 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5233 pivot_value_add_footnote (v, footnote);
5234 pivot_table_put2 (table, y, x, v);
5237 pivot_table_submit (table);
5241 matrix_cmd_execute_print (const struct print_command *print)
5243 if (print->expression)
5245 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5250 struct fmt_spec format = (print->use_default_format
5251 ? default_format (m, &log_scale)
5254 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5255 matrix_cmd_print_text (print, m, format, log_scale);
5257 matrix_cmd_print_tables (print, m, format, log_scale);
5259 gsl_matrix_free (m);
5263 matrix_cmd_print_space (print->space);
5265 output_log ("%s", print->title);
5272 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
5273 const char *command_name,
5274 const char *stop1, const char *stop2)
5276 lex_end_of_command (s->lexer);
5277 lex_discard_rest_of_command (s->lexer);
5279 size_t allocated = 0;
5282 while (lex_token (s->lexer) == T_ENDCMD)
5285 if (lex_at_phrase (s->lexer, stop1)
5286 || (stop2 && lex_at_phrase (s->lexer, stop2)))
5289 if (lex_at_phrase (s->lexer, "END MATRIX"))
5291 msg (SE, _("Premature END MATRIX within %s."), command_name);
5295 struct matrix_cmd *cmd = matrix_parse_command (s);
5299 if (c->n >= allocated)
5300 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
5301 c->commands[c->n++] = cmd;
5306 matrix_parse_do_if_clause (struct matrix_state *s,
5307 struct do_if_command *ifc,
5308 bool parse_condition,
5309 size_t *allocated_clauses)
5311 if (ifc->n_clauses >= *allocated_clauses)
5312 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5313 sizeof *ifc->clauses);
5314 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5315 *c = (struct do_if_clause) { .condition = NULL };
5317 if (parse_condition)
5319 c->condition = matrix_parse_expr (s);
5324 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
5327 static struct matrix_cmd *
5328 matrix_parse_do_if (struct matrix_state *s)
5330 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5331 *cmd = (struct matrix_cmd) {
5333 .do_if = { .n_clauses = 0 }
5336 size_t allocated_clauses = 0;
5339 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
5342 while (lex_match_phrase (s->lexer, "ELSE IF"));
5344 if (lex_match_id (s->lexer, "ELSE")
5345 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
5348 if (!lex_match_phrase (s->lexer, "END IF"))
5353 matrix_cmd_destroy (cmd);
5358 matrix_cmd_execute_do_if (struct do_if_command *cmd)
5360 for (size_t i = 0; i < cmd->n_clauses; i++)
5362 struct do_if_clause *c = &cmd->clauses[i];
5366 if (!matrix_expr_evaluate_scalar (c->condition,
5367 i ? "ELSE IF" : "DO IF",
5372 for (size_t j = 0; j < c->commands.n; j++)
5373 if (!matrix_cmd_execute (c->commands.commands[j]))
5380 static struct matrix_cmd *
5381 matrix_parse_loop (struct matrix_state *s)
5383 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5384 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5386 struct loop_command *loop = &cmd->loop;
5387 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5389 struct substring name = lex_tokss (s->lexer);
5390 loop->var = matrix_var_lookup (s, name);
5392 loop->var = matrix_var_create (s, name);
5397 loop->start = matrix_parse_expr (s);
5398 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5401 loop->end = matrix_parse_expr (s);
5405 if (lex_match (s->lexer, T_BY))
5407 loop->increment = matrix_parse_expr (s);
5408 if (!loop->increment)
5413 if (lex_match_id (s->lexer, "IF"))
5415 loop->top_condition = matrix_parse_expr (s);
5416 if (!loop->top_condition)
5420 bool was_in_loop = s->in_loop;
5422 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
5424 s->in_loop = was_in_loop;
5428 if (!lex_match_phrase (s->lexer, "END LOOP"))
5431 if (lex_match_id (s->lexer, "IF"))
5433 loop->bottom_condition = matrix_parse_expr (s);
5434 if (!loop->bottom_condition)
5441 matrix_cmd_destroy (cmd);
5446 matrix_cmd_execute_loop (struct loop_command *cmd)
5450 long int increment = 1;
5453 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5454 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5456 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5460 if (increment > 0 ? value > end
5461 : increment < 0 ? value < end
5466 int mxloops = settings_get_mxloops ();
5467 for (int i = 0; i < mxloops; i++)
5472 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5474 gsl_matrix_free (cmd->var->value);
5475 cmd->var->value = NULL;
5477 if (!cmd->var->value)
5478 cmd->var->value = gsl_matrix_alloc (1, 1);
5480 gsl_matrix_set (cmd->var->value, 0, 0, value);
5483 if (cmd->top_condition)
5486 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5491 for (size_t j = 0; j < cmd->commands.n; j++)
5492 if (!matrix_cmd_execute (cmd->commands.commands[j]))
5495 if (cmd->bottom_condition)
5498 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5499 "END LOOP IF", &d) || d > 0)
5506 if (increment > 0 ? value > end : value < end)
5512 static struct matrix_cmd *
5513 matrix_parse_break (struct matrix_state *s)
5517 msg (SE, _("BREAK not inside LOOP."));
5521 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5522 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
5526 static struct matrix_cmd *
5527 matrix_parse_display (struct matrix_state *s)
5529 while (lex_token (s->lexer) != T_ENDCMD)
5531 if (!lex_match_id (s->lexer, "DICTIONARY")
5532 && !lex_match_id (s->lexer, "STATUS"))
5534 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5539 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5540 *cmd = (struct matrix_cmd) { .type = MCMD_DISPLAY, .display = { s } };
5545 compare_matrix_var_pointers (const void *a_, const void *b_)
5547 const struct matrix_var *const *ap = a_;
5548 const struct matrix_var *const *bp = b_;
5549 const struct matrix_var *a = *ap;
5550 const struct matrix_var *b = *bp;
5551 return strcmp (a->name, b->name);
5555 matrix_cmd_execute_display (struct display_command *cmd)
5557 const struct matrix_state *s = cmd->state;
5559 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5560 struct pivot_dimension *attr_dimension
5561 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5562 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5563 N_("Rows"), N_("Columns"));
5564 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5566 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5569 const struct matrix_var *var;
5570 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5572 vars[n_vars++] = var;
5573 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5575 struct pivot_dimension *rows = pivot_dimension_create (
5576 table, PIVOT_AXIS_ROW, N_("Variable"));
5577 for (size_t i = 0; i < n_vars; i++)
5579 const struct matrix_var *var = vars[i];
5580 pivot_category_create_leaf (
5581 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
5583 size_t r = var->value->size1;
5584 size_t c = var->value->size2;
5585 double values[] = { r, c, r * c * sizeof (double) / 1024 };
5586 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
5587 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
5590 pivot_table_submit (table);
5593 static struct matrix_cmd *
5594 matrix_parse_release (struct matrix_state *s)
5596 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5597 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
5599 size_t allocated_vars = 0;
5600 while (lex_token (s->lexer) == T_ID)
5602 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
5605 if (cmd->release.n_vars >= allocated_vars)
5606 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
5607 sizeof *cmd->release.vars);
5608 cmd->release.vars[cmd->release.n_vars++] = var;
5611 lex_error (s->lexer, _("Variable name expected."));
5614 if (!lex_match (s->lexer, T_COMMA))
5622 matrix_cmd_execute_release (struct release_command *cmd)
5624 for (size_t i = 0; i < cmd->n_vars; i++)
5626 struct matrix_var *var = cmd->vars[i];
5627 gsl_matrix_free (var->value);
5632 static struct save_file *
5633 save_file_create (struct matrix_state *s, struct file_handle *fh,
5634 struct string_array *variables,
5635 struct matrix_expr *names,
5636 struct stringi_set *strings)
5638 for (size_t i = 0; i < s->n_save_files; i++)
5640 struct save_file *sf = s->save_files[i];
5641 if (fh_equal (sf->file, fh))
5645 string_array_destroy (variables);
5646 matrix_expr_destroy (names);
5647 stringi_set_destroy (strings);
5653 struct save_file *sf = xmalloc (sizeof *sf);
5654 *sf = (struct save_file) {
5656 .dataset = s->dataset,
5657 .variables = *variables,
5659 .strings = STRINGI_SET_INITIALIZER (sf->strings),
5662 stringi_set_swap (&sf->strings, strings);
5663 stringi_set_destroy (strings);
5665 s->save_files = xrealloc (s->save_files,
5666 (s->n_save_files + 1) * sizeof *s->save_files);
5667 s->save_files[s->n_save_files++] = sf;
5671 static struct casewriter *
5672 save_file_open (struct save_file *sf, gsl_matrix *m)
5674 if (sf->writer || sf->error)
5678 size_t n_variables = caseproto_get_n_widths (
5679 casewriter_get_proto (sf->writer));
5680 if (m->size2 != n_variables)
5682 msg (ME, _("The first SAVE to %s within this matrix program "
5683 "had %zu columns, so a %zu×%zu matrix cannot be "
5685 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
5686 n_variables, m->size1, m->size2);
5694 struct dictionary *dict = dict_create (get_default_encoding ());
5696 /* Fill 'names' with user-specified names if there were any, either from
5697 sf->variables or sf->names. */
5698 struct string_array names = { .n = 0 };
5699 if (sf->variables.n)
5700 string_array_clone (&names, &sf->variables);
5703 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
5704 if (nm && is_vector (nm))
5706 gsl_vector nv = to_vector (nm);
5707 for (size_t i = 0; i < nv.size; i++)
5709 char *name = trimmed_string (gsl_vector_get (&nv, i));
5710 if (dict_id_is_valid (dict, name, true))
5711 string_array_append_nocopy (&names, name);
5716 gsl_matrix_free (nm);
5719 struct stringi_set strings;
5720 stringi_set_clone (&strings, &sf->strings);
5722 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
5727 name = names.strings[i];
5730 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
5734 int width = stringi_set_delete (&strings, name) ? 8 : 0;
5735 struct variable *var = dict_create_var (dict, name, width);
5738 msg (ME, _("Duplicate variable name %s in SAVE statement."),
5744 if (!stringi_set_is_empty (&strings))
5746 size_t count = stringi_set_count (&strings);
5747 const char *example = stringi_set_node_get_string (
5748 stringi_set_first (&strings));
5751 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
5752 "variable %s."), example);
5754 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
5755 "unknown variable: %s.",
5756 "The SAVE command STRINGS subcommand specifies %zu "
5757 "unknown variables, including %s.",
5762 stringi_set_destroy (&strings);
5763 string_array_destroy (&names);
5772 if (sf->file == fh_inline_file ())
5773 sf->writer = autopaging_writer_create (dict_get_proto (dict));
5775 sf->writer = any_writer_open (sf->file, dict);
5787 save_file_destroy (struct save_file *sf)
5791 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
5793 dataset_set_dict (sf->dataset, sf->dict);
5794 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
5798 casewriter_destroy (sf->writer);
5799 dict_unref (sf->dict);
5801 fh_unref (sf->file);
5802 string_array_destroy (&sf->variables);
5803 matrix_expr_destroy (sf->names);
5804 stringi_set_destroy (&sf->strings);
5809 static struct matrix_cmd *
5810 matrix_parse_save (struct matrix_state *s)
5812 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5813 *cmd = (struct matrix_cmd) {
5815 .save = { .expression = NULL },
5818 struct file_handle *fh = NULL;
5819 struct save_command *save = &cmd->save;
5821 struct string_array variables = STRING_ARRAY_INITIALIZER;
5822 struct matrix_expr *names = NULL;
5823 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
5825 save->expression = matrix_parse_exp (s);
5826 if (!save->expression)
5829 while (lex_match (s->lexer, T_SLASH))
5831 if (lex_match_id (s->lexer, "OUTFILE"))
5833 lex_match (s->lexer, T_EQUALS);
5836 fh = (lex_match (s->lexer, T_ASTERISK)
5838 : fh_parse (s->lexer, FH_REF_FILE, s->session));
5842 else if (lex_match_id (s->lexer, "VARIABLES"))
5844 lex_match (s->lexer, T_EQUALS);
5848 struct dictionary *d = dict_create (get_default_encoding ());
5849 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
5850 PV_NO_SCRATCH | PV_NO_DUPLICATE);
5855 string_array_clear (&variables);
5856 variables = (struct string_array) {
5862 else if (lex_match_id (s->lexer, "NAMES"))
5864 lex_match (s->lexer, T_EQUALS);
5865 matrix_expr_destroy (names);
5866 names = matrix_parse_exp (s);
5870 else if (lex_match_id (s->lexer, "STRINGS"))
5872 lex_match (s->lexer, T_EQUALS);
5873 while (lex_token (s->lexer) == T_ID)
5875 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
5877 if (!lex_match (s->lexer, T_COMMA))
5883 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
5891 if (s->prev_save_file)
5892 fh = fh_ref (s->prev_save_file);
5895 lex_sbc_missing ("OUTFILE");
5899 fh_unref (s->prev_save_file);
5900 s->prev_save_file = fh_ref (fh);
5902 if (variables.n && names)
5904 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
5905 matrix_expr_destroy (names);
5909 save->sf = save_file_create (s, fh, &variables, names, &strings);
5913 string_array_destroy (&variables);
5914 matrix_expr_destroy (names);
5915 stringi_set_destroy (&strings);
5917 matrix_cmd_destroy (cmd);
5922 matrix_cmd_execute_save (const struct save_command *save)
5924 gsl_matrix *m = matrix_expr_evaluate (save->expression);
5928 struct casewriter *writer = save_file_open (save->sf, m);
5931 gsl_matrix_free (m);
5935 const struct caseproto *proto = casewriter_get_proto (writer);
5936 for (size_t y = 0; y < m->size1; y++)
5938 struct ccase *c = case_create (proto);
5939 for (size_t x = 0; x < m->size2; x++)
5941 double d = gsl_matrix_get (m, y, x);
5942 int width = caseproto_get_width (proto, x);
5943 union value *value = case_data_rw_idx (c, x);
5947 memcpy (value->s, &d, width);
5949 casewriter_write (writer, c);
5951 gsl_matrix_free (m);
5954 static struct matrix_cmd *
5955 matrix_parse_read (struct matrix_state *s)
5957 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5958 *cmd = (struct matrix_cmd) {
5960 .read = { .format = FMT_F },
5963 struct file_handle *fh = NULL;
5964 char *encoding = NULL;
5965 struct read_command *read = &cmd->read;
5966 read->dst = matrix_lvalue_parse (s);
5971 int repetitions = 0;
5972 int record_width = 0;
5973 bool seen_format = false;
5974 while (lex_match (s->lexer, T_SLASH))
5976 if (lex_match_id (s->lexer, "FILE"))
5978 lex_match (s->lexer, T_EQUALS);
5981 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5985 else if (lex_match_id (s->lexer, "ENCODING"))
5987 lex_match (s->lexer, T_EQUALS);
5988 if (!lex_force_string (s->lexer))
5992 encoding = ss_xstrdup (lex_tokss (s->lexer));
5996 else if (lex_match_id (s->lexer, "FIELD"))
5998 lex_match (s->lexer, T_EQUALS);
6000 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6002 read->c1 = lex_integer (s->lexer);
6004 if (!lex_force_match (s->lexer, T_TO)
6005 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
6007 read->c2 = lex_integer (s->lexer) + 1;
6010 record_width = read->c2 - read->c1;
6011 if (lex_match (s->lexer, T_BY))
6013 if (!lex_force_int_range (s->lexer, "BY", 1,
6014 read->c2 - read->c1))
6016 by = lex_integer (s->lexer);
6019 if (record_width % by)
6021 msg (SE, _("BY %d does not evenly divide record width %d."),
6029 else if (lex_match_id (s->lexer, "SIZE"))
6031 lex_match (s->lexer, T_EQUALS);
6032 matrix_expr_destroy (read->size);
6033 read->size = matrix_parse_exp (s);
6037 else if (lex_match_id (s->lexer, "MODE"))
6039 lex_match (s->lexer, T_EQUALS);
6040 if (lex_match_id (s->lexer, "RECTANGULAR"))
6041 read->symmetric = false;
6042 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6043 read->symmetric = true;
6046 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6050 else if (lex_match_id (s->lexer, "REREAD"))
6051 read->reread = true;
6052 else if (lex_match_id (s->lexer, "FORMAT"))
6056 lex_sbc_only_once ("FORMAT");
6061 lex_match (s->lexer, T_EQUALS);
6063 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6066 const char *p = lex_tokcstr (s->lexer);
6067 if (c_isdigit (p[0]))
6069 repetitions = atoi (p);
6070 p += strspn (p, "0123456789");
6071 if (!fmt_from_name (p, &read->format))
6073 lex_error (s->lexer, _("Unknown format %s."), p);
6078 else if (fmt_from_name (p, &read->format))
6082 struct fmt_spec format;
6083 if (!parse_format_specifier (s->lexer, &format))
6085 read->format = format.type;
6091 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6092 "REREAD", "FORMAT");
6099 lex_sbc_missing ("FIELD");
6103 if (!read->dst->n_indexes && !read->size)
6105 msg (SE, _("SIZE is required for reading data into a full matrix "
6106 "(as opposed to a submatrix)."));
6112 if (s->prev_read_file)
6113 fh = fh_ref (s->prev_read_file);
6116 lex_sbc_missing ("FILE");
6120 fh_unref (s->prev_read_file);
6121 s->prev_read_file = fh_ref (fh);
6123 read->rf = read_file_create (s, fh);
6127 free (read->rf->encoding);
6128 read->rf->encoding = encoding;
6132 /* Field width may be specified in multiple ways:
6135 2. The format on FORMAT.
6136 3. The repetition factor on FORMAT.
6138 (2) and (3) are mutually exclusive.
6140 If more than one of these is present, they must agree. If none of them is
6141 present, then free-field format is used.
6143 if (repetitions > record_width)
6145 msg (SE, _("%d repetitions cannot fit in record width %d."),
6146 repetitions, record_width);
6149 int w = (repetitions ? record_width / repetitions
6155 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6156 "which implies field width %d, "
6157 "but BY specifies field width %d."),
6158 repetitions, record_width, w, by);
6160 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6169 matrix_cmd_destroy (cmd);
6175 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6176 struct substring data, size_t y, size_t x,
6177 int first_column, int last_column, char *error)
6179 int line_number = dfm_get_line_number (reader);
6180 struct msg_location *location = xmalloc (sizeof *location);
6181 *location = (struct msg_location) {
6182 .file_name = xstrdup (dfm_get_file_name (reader)),
6183 .first_line = line_number,
6184 .last_line = line_number + 1,
6185 .first_column = first_column,
6186 .last_column = last_column,
6188 struct msg *m = xmalloc (sizeof *m);
6190 .category = MSG_C_DATA,
6191 .severity = MSG_S_WARNING,
6192 .location = location,
6193 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
6194 "for matrix row %zu, column %zu: %s"),
6195 (int) data.length, data.string, fmt_name (format),
6196 y + 1, x + 1, error),
6204 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
6205 gsl_matrix *m, struct substring p, size_t y, size_t x,
6206 const char *line_start)
6208 const char *input_encoding = dfm_reader_get_encoding (reader);
6211 if (fmt_is_numeric (read->format))
6214 error = data_in (p, input_encoding, read->format,
6215 settings_get_fmt_settings (), &v, 0, NULL);
6216 if (!error && v.f == SYSMIS)
6217 error = xstrdup (_("Matrix data may not contain missing value."));
6222 uint8_t s[sizeof (double)];
6223 union value v = { .s = s };
6224 error = data_in (p, input_encoding, read->format,
6225 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6226 memcpy (&f, s, sizeof f);
6231 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6232 int c2 = c1 + ss_utf8_count_columns (p) - 1;
6233 parse_error (reader, read->format, p, y, x, c1, c2, error);
6237 gsl_matrix_set (m, y, x, f);
6238 if (read->symmetric && x != y)
6239 gsl_matrix_set (m, x, y, f);
6244 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
6245 struct substring *line, const char **startp)
6247 if (dfm_eof (reader))
6249 msg (SE, _("Unexpected end of file reading matrix data."));
6252 dfm_expand_tabs (reader);
6253 struct substring record = dfm_get_record (reader);
6254 /* XXX need to recode record into UTF-8 */
6255 *startp = record.string;
6256 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6261 matrix_read (struct read_command *read, struct dfm_reader *reader,
6264 for (size_t y = 0; y < m->size1; y++)
6266 size_t nx = read->symmetric ? y + 1 : m->size2;
6268 struct substring line = ss_empty ();
6269 const char *line_start = line.string;
6270 for (size_t x = 0; x < nx; x++)
6277 ss_ltrim (&line, ss_cstr (" ,"));
6278 if (!ss_is_empty (line))
6280 if (!matrix_read_line (read, reader, &line, &line_start))
6282 dfm_forward_record (reader);
6285 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6289 if (!matrix_read_line (read, reader, &line, &line_start))
6291 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6292 int f = x % fields_per_line;
6293 if (f == fields_per_line - 1)
6294 dfm_forward_record (reader);
6296 p = ss_substr (line, read->w * f, read->w);
6299 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6303 dfm_forward_record (reader);
6306 ss_ltrim (&line, ss_cstr (" ,"));
6307 if (!ss_is_empty (line))
6310 msg (SW, _("Trailing garbage on line \"%.*s\""),
6311 (int) line.length, line.string);
6318 matrix_cmd_execute_read (struct read_command *read)
6320 struct index_vector iv0, iv1;
6321 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6324 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6327 gsl_matrix *m = matrix_expr_evaluate (read->size);
6333 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6334 "not a %zu×%zu matrix."), m->size1, m->size2);
6335 gsl_matrix_free (m);
6341 gsl_vector v = to_vector (m);
6345 d[0] = gsl_vector_get (&v, 0);
6348 else if (v.size == 2)
6350 d[0] = gsl_vector_get (&v, 0);
6351 d[1] = gsl_vector_get (&v, 1);
6355 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6356 "not a %zu×%zu matrix."),
6357 m->size1, m->size2),
6358 gsl_matrix_free (m);
6363 gsl_matrix_free (m);
6365 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6367 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
6368 "are outside valid range."),
6379 if (read->dst->n_indexes)
6381 size_t submatrix_size[2];
6382 if (read->dst->n_indexes == 2)
6384 submatrix_size[0] = iv0.n;
6385 submatrix_size[1] = iv1.n;
6387 else if (read->dst->var->value->size1 == 1)
6389 submatrix_size[0] = 1;
6390 submatrix_size[1] = iv0.n;
6394 submatrix_size[0] = iv0.n;
6395 submatrix_size[1] = 1;
6400 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6402 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
6403 "differ from submatrix dimensions %zu×%zu."),
6405 submatrix_size[0], submatrix_size[1]);
6413 size[0] = submatrix_size[0];
6414 size[1] = submatrix_size[1];
6418 struct dfm_reader *reader = read_file_open (read->rf);
6420 dfm_reread_record (reader, 1);
6422 if (read->symmetric && size[0] != size[1])
6424 msg (SE, _("Cannot read non-square %zu×%zu matrix "
6425 "using READ with MODE=SYMMETRIC."),
6431 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6432 matrix_read (read, reader, tmp);
6433 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
6436 static struct matrix_cmd *
6437 matrix_parse_write (struct matrix_state *s)
6439 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6440 *cmd = (struct matrix_cmd) {
6444 struct file_handle *fh = NULL;
6445 char *encoding = NULL;
6446 struct write_command *write = &cmd->write;
6447 write->expression = matrix_parse_exp (s);
6448 if (!write->expression)
6452 int repetitions = 0;
6453 int record_width = 0;
6454 enum fmt_type format = FMT_F;
6455 bool has_format = false;
6456 while (lex_match (s->lexer, T_SLASH))
6458 if (lex_match_id (s->lexer, "OUTFILE"))
6460 lex_match (s->lexer, T_EQUALS);
6463 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6467 else if (lex_match_id (s->lexer, "ENCODING"))
6469 lex_match (s->lexer, T_EQUALS);
6470 if (!lex_force_string (s->lexer))
6474 encoding = ss_xstrdup (lex_tokss (s->lexer));
6478 else if (lex_match_id (s->lexer, "FIELD"))
6480 lex_match (s->lexer, T_EQUALS);
6482 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6484 write->c1 = lex_integer (s->lexer);
6486 if (!lex_force_match (s->lexer, T_TO)
6487 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
6489 write->c2 = lex_integer (s->lexer) + 1;
6492 record_width = write->c2 - write->c1;
6493 if (lex_match (s->lexer, T_BY))
6495 if (!lex_force_int_range (s->lexer, "BY", 1,
6496 write->c2 - write->c1))
6498 by = lex_integer (s->lexer);
6501 if (record_width % by)
6503 msg (SE, _("BY %d does not evenly divide record width %d."),
6511 else if (lex_match_id (s->lexer, "MODE"))
6513 lex_match (s->lexer, T_EQUALS);
6514 if (lex_match_id (s->lexer, "RECTANGULAR"))
6515 write->triangular = false;
6516 else if (lex_match_id (s->lexer, "TRIANGULAR"))
6517 write->triangular = true;
6520 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
6524 else if (lex_match_id (s->lexer, "HOLD"))
6526 else if (lex_match_id (s->lexer, "FORMAT"))
6528 if (has_format || write->format)
6530 lex_sbc_only_once ("FORMAT");
6534 lex_match (s->lexer, T_EQUALS);
6536 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6539 const char *p = lex_tokcstr (s->lexer);
6540 if (c_isdigit (p[0]))
6542 repetitions = atoi (p);
6543 p += strspn (p, "0123456789");
6544 if (!fmt_from_name (p, &format))
6546 lex_error (s->lexer, _("Unknown format %s."), p);
6552 else if (fmt_from_name (p, &format))
6559 struct fmt_spec spec;
6560 if (!parse_format_specifier (s->lexer, &spec))
6562 write->format = xmemdup (&spec, sizeof spec);
6567 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
6575 lex_sbc_missing ("FIELD");
6581 if (s->prev_write_file)
6582 fh = fh_ref (s->prev_write_file);
6585 lex_sbc_missing ("OUTFILE");
6589 fh_unref (s->prev_write_file);
6590 s->prev_write_file = fh_ref (fh);
6592 write->wf = write_file_create (s, fh);
6596 free (write->wf->encoding);
6597 write->wf->encoding = encoding;
6601 /* Field width may be specified in multiple ways:
6604 2. The format on FORMAT.
6605 3. The repetition factor on FORMAT.
6607 (2) and (3) are mutually exclusive.
6609 If more than one of these is present, they must agree. If none of them is
6610 present, then free-field format is used.
6612 if (repetitions > record_width)
6614 msg (SE, _("%d repetitions cannot fit in record width %d."),
6615 repetitions, record_width);
6618 int w = (repetitions ? record_width / repetitions
6619 : write->format ? write->format->w
6624 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6625 "which implies field width %d, "
6626 "but BY specifies field width %d."),
6627 repetitions, record_width, w, by);
6629 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6633 if (w && !write->format)
6635 write->format = xmalloc (sizeof *write->format);
6636 *write->format = (struct fmt_spec) { .type = format, .w = w };
6638 if (!fmt_check_output (write->format))
6642 if (write->format && fmt_var_width (write->format) > sizeof (double))
6644 char s[FMT_STRING_LEN_MAX + 1];
6645 fmt_to_string (write->format, s);
6646 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
6647 s, sizeof (double));
6655 matrix_cmd_destroy (cmd);
6660 matrix_cmd_execute_write (struct write_command *write)
6662 gsl_matrix *m = matrix_expr_evaluate (write->expression);
6666 if (write->triangular && m->size1 != m->size2)
6668 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
6669 "the matrix to be written has dimensions %zu×%zu."),
6670 m->size1, m->size2);
6671 gsl_matrix_free (m);
6675 struct dfm_writer *writer = write_file_open (write->wf);
6676 if (!writer || !m->size1)
6678 gsl_matrix_free (m);
6682 const struct fmt_settings *settings = settings_get_fmt_settings ();
6683 struct u8_line *line = write->wf->held;
6684 for (size_t y = 0; y < m->size1; y++)
6688 line = xmalloc (sizeof *line);
6689 u8_line_init (line);
6691 size_t nx = write->triangular ? y + 1 : m->size2;
6693 for (size_t x = 0; x < nx; x++)
6696 double f = gsl_matrix_get (m, y, x);
6700 if (fmt_is_numeric (write->format->type))
6703 v.s = (uint8_t *) &f;
6704 s = data_out (&v, NULL, write->format, settings);
6708 s = xmalloc (DBL_BUFSIZE_BOUND);
6709 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
6710 >= DBL_BUFSIZE_BOUND)
6713 size_t len = strlen (s);
6714 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
6715 if (width + x0 > write->c2)
6717 dfm_put_record_utf8 (writer, line->s.ss.string,
6719 u8_line_clear (line);
6722 u8_line_put (line, x0, x0 + width, s, len);
6725 x0 += write->format ? write->format->w : width + 1;
6728 if (y + 1 >= m->size1 && write->hold)
6730 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
6731 u8_line_clear (line);
6735 u8_line_destroy (line);
6739 write->wf->held = line;
6741 gsl_matrix_free (m);
6744 static struct matrix_cmd *
6745 matrix_parse_get (struct matrix_state *s)
6747 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6748 *cmd = (struct matrix_cmd) {
6751 .dataset = s->dataset,
6752 .user = { .treatment = MGET_ERROR },
6753 .system = { .treatment = MGET_ERROR },
6757 struct get_command *get = &cmd->get;
6758 get->dst = matrix_lvalue_parse (s);
6762 while (lex_match (s->lexer, T_SLASH))
6764 if (lex_match_id (s->lexer, "FILE"))
6766 lex_match (s->lexer, T_EQUALS);
6768 fh_unref (get->file);
6769 if (lex_match (s->lexer, T_ASTERISK))
6773 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6778 else if (lex_match_id (s->lexer, "ENCODING"))
6780 lex_match (s->lexer, T_EQUALS);
6781 if (!lex_force_string (s->lexer))
6784 free (get->encoding);
6785 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
6789 else if (lex_match_id (s->lexer, "VARIABLES"))
6791 lex_match (s->lexer, T_EQUALS);
6795 lex_sbc_only_once ("VARIABLES");
6799 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
6802 else if (lex_match_id (s->lexer, "NAMES"))
6804 lex_match (s->lexer, T_EQUALS);
6805 if (!lex_force_id (s->lexer))
6808 struct substring name = lex_tokss (s->lexer);
6809 get->names = matrix_var_lookup (s, name);
6811 get->names = matrix_var_create (s, name);
6814 else if (lex_match_id (s->lexer, "MISSING"))
6816 lex_match (s->lexer, T_EQUALS);
6817 if (lex_match_id (s->lexer, "ACCEPT"))
6818 get->user.treatment = MGET_ACCEPT;
6819 else if (lex_match_id (s->lexer, "OMIT"))
6820 get->user.treatment = MGET_OMIT;
6821 else if (lex_is_number (s->lexer))
6823 get->user.treatment = MGET_RECODE;
6824 get->user.substitute = lex_number (s->lexer);
6829 lex_error (s->lexer, NULL);
6833 else if (lex_match_id (s->lexer, "SYSMIS"))
6835 lex_match (s->lexer, T_EQUALS);
6836 if (lex_match_id (s->lexer, "OMIT"))
6837 get->system.treatment = MGET_OMIT;
6838 else if (lex_is_number (s->lexer))
6840 get->system.treatment = MGET_RECODE;
6841 get->system.substitute = lex_number (s->lexer);
6846 lex_error (s->lexer, NULL);
6852 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
6853 "MISSING", "SYSMIS");
6858 if (get->user.treatment != MGET_ACCEPT)
6859 get->system.treatment = MGET_ERROR;
6864 matrix_cmd_destroy (cmd);
6869 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
6870 const struct dictionary *dict)
6872 struct variable **vars;
6877 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
6878 &vars, &n_vars, PV_NUMERIC))
6883 n_vars = dict_get_var_cnt (dict);
6884 vars = xnmalloc (n_vars, sizeof *vars);
6885 for (size_t i = 0; i < n_vars; i++)
6887 struct variable *var = dict_get_var (dict, i);
6888 if (!var_is_numeric (var))
6890 msg (SE, _("GET: Variable %s is not numeric."),
6891 var_get_name (var));
6901 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
6902 for (size_t i = 0; i < n_vars; i++)
6904 char s[sizeof (double)];
6906 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
6907 memcpy (&f, s, sizeof f);
6908 gsl_matrix_set (names, i, 0, f);
6911 gsl_matrix_free (get->names->value);
6912 get->names->value = names;
6916 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
6917 long long int casenum = 1;
6919 for (struct ccase *c = casereader_read (reader); c;
6920 c = casereader_read (reader), casenum++)
6922 if (n_rows >= m->size1)
6924 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
6925 for (size_t y = 0; y < n_rows; y++)
6926 for (size_t x = 0; x < n_vars; x++)
6927 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
6928 gsl_matrix_free (m);
6933 for (size_t x = 0; x < n_vars; x++)
6935 const struct variable *var = vars[x];
6936 double d = case_num (c, var);
6939 if (get->system.treatment == MGET_RECODE)
6940 d = get->system.substitute;
6941 else if (get->system.treatment == MGET_OMIT)
6945 msg (SE, _("GET: Variable %s in case %lld "
6946 "is system-missing."),
6947 var_get_name (var), casenum);
6951 else if (var_is_num_missing (var, d, MV_USER))
6953 if (get->user.treatment == MGET_RECODE)
6954 d = get->user.substitute;
6955 else if (get->user.treatment == MGET_OMIT)
6957 else if (get->user.treatment != MGET_ACCEPT)
6959 msg (SE, _("GET: Variable %s in case %lld has user-missing "
6961 var_get_name (var), casenum, d);
6965 gsl_matrix_set (m, n_rows, x, d);
6976 matrix_lvalue_evaluate_and_assign (get->dst, m);
6979 gsl_matrix_free (m);
6984 matrix_open_casereader (const char *command_name,
6985 struct file_handle *file, const char *encoding,
6986 struct dataset *dataset,
6987 struct casereader **readerp, struct dictionary **dictp)
6991 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
6992 return *readerp != NULL;
6996 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
6998 msg (ME, _("The %s command cannot read an empty active file."),
7002 *readerp = proc_open (dataset);
7003 *dictp = dict_ref (dataset_dict (dataset));
7009 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
7010 struct casereader *reader, struct dictionary *dict)
7013 casereader_destroy (reader);
7015 proc_commit (dataset);
7019 matrix_cmd_execute_get (struct get_command *get)
7021 struct casereader *r;
7022 struct dictionary *d;
7023 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
7026 matrix_cmd_execute_get__ (get, r, d);
7027 matrix_close_casereader (get->file, get->dataset, r, d);
7032 match_rowtype (struct lexer *lexer)
7034 static const char *rowtypes[] = {
7035 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7037 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7039 for (size_t i = 0; i < n_rowtypes; i++)
7040 if (lex_match_id (lexer, rowtypes[i]))
7043 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7048 parse_var_names (struct lexer *lexer, struct string_array *sa)
7050 lex_match (lexer, T_EQUALS);
7052 string_array_clear (sa);
7054 struct dictionary *dict = dict_create (get_default_encoding ());
7057 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7058 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7063 for (size_t i = 0; i < n_names; i++)
7064 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7065 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7067 msg (SE, _("Variable name %s is reserved."), names[i]);
7068 for (size_t j = 0; j < n_names; j++)
7074 string_array_clear (sa);
7075 sa->strings = names;
7076 sa->n = sa->allocated = n_names;
7082 msave_common_destroy (struct msave_common *common)
7086 fh_unref (common->outfile);
7087 string_array_destroy (&common->variables);
7088 string_array_destroy (&common->fnames);
7089 string_array_destroy (&common->snames);
7091 for (size_t i = 0; i < common->n_factors; i++)
7092 matrix_expr_destroy (common->factors[i]);
7093 free (common->factors);
7095 for (size_t i = 0; i < common->n_splits; i++)
7096 matrix_expr_destroy (common->splits[i]);
7097 free (common->splits);
7099 dict_unref (common->dict);
7100 casewriter_destroy (common->writer);
7107 compare_variables (const char *keyword,
7108 const struct string_array *new,
7109 const struct string_array *old)
7115 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7116 "on the first MSAVE within MATRIX."), keyword);
7119 else if (!string_array_equal_case (old, new))
7121 msg (SE, _("%s must specify the same variables each time within "
7122 "a given MATRIX."), keyword);
7129 static struct matrix_cmd *
7130 matrix_parse_msave (struct matrix_state *s)
7132 struct msave_common *common = xmalloc (sizeof *common);
7133 *common = (struct msave_common) { .outfile = NULL };
7135 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7136 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7138 struct matrix_expr *splits = NULL;
7139 struct matrix_expr *factors = NULL;
7141 struct msave_command *msave = &cmd->msave;
7142 msave->expr = matrix_parse_exp (s);
7146 while (lex_match (s->lexer, T_SLASH))
7148 if (lex_match_id (s->lexer, "TYPE"))
7150 lex_match (s->lexer, T_EQUALS);
7152 msave->rowtype = match_rowtype (s->lexer);
7153 if (!msave->rowtype)
7156 else if (lex_match_id (s->lexer, "OUTFILE"))
7158 lex_match (s->lexer, T_EQUALS);
7160 fh_unref (common->outfile);
7161 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7162 if (!common->outfile)
7165 else if (lex_match_id (s->lexer, "VARIABLES"))
7167 if (!parse_var_names (s->lexer, &common->variables))
7170 else if (lex_match_id (s->lexer, "FNAMES"))
7172 if (!parse_var_names (s->lexer, &common->fnames))
7175 else if (lex_match_id (s->lexer, "SNAMES"))
7177 if (!parse_var_names (s->lexer, &common->snames))
7180 else if (lex_match_id (s->lexer, "SPLIT"))
7182 lex_match (s->lexer, T_EQUALS);
7184 matrix_expr_destroy (splits);
7185 splits = matrix_parse_exp (s);
7189 else if (lex_match_id (s->lexer, "FACTOR"))
7191 lex_match (s->lexer, T_EQUALS);
7193 matrix_expr_destroy (factors);
7194 factors = matrix_parse_exp (s);
7200 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7201 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7205 if (!msave->rowtype)
7207 lex_sbc_missing ("TYPE");
7213 if (common->fnames.n && !factors)
7215 msg (SE, _("FNAMES requires FACTOR."));
7218 if (common->snames.n && !splits)
7220 msg (SE, _("SNAMES requires SPLIT."));
7223 if (!common->outfile)
7225 lex_sbc_missing ("OUTFILE");
7232 if (common->outfile)
7234 if (!fh_equal (common->outfile, s->common->outfile))
7236 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7237 "within a single MATRIX command."));
7241 if (!compare_variables ("VARIABLES",
7242 &common->variables, &s->common->variables)
7243 || !compare_variables ("FNAMES", &common->fnames, &s->common->fnames)
7244 || !compare_variables ("SNAMES", &common->snames, &s->common->snames))
7246 msave_common_destroy (common);
7248 msave->common = s->common;
7250 struct msave_common *c = s->common;
7253 if (c->n_factors >= c->allocated_factors)
7254 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7255 sizeof *c->factors);
7256 c->factors[c->n_factors++] = factors;
7258 if (c->n_factors > 0)
7259 msave->factors = c->factors[c->n_factors - 1];
7263 if (c->n_splits >= c->allocated_splits)
7264 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7266 c->splits[c->n_splits++] = splits;
7268 if (c->n_splits > 0)
7269 msave->splits = c->splits[c->n_splits - 1];
7274 matrix_expr_destroy (splits);
7275 matrix_expr_destroy (factors);
7276 msave_common_destroy (common);
7277 matrix_cmd_destroy (cmd);
7282 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7284 gsl_matrix *m = matrix_expr_evaluate (e);
7290 msg (SE, _("%s expression must evaluate to vector, "
7291 "not a %zu×%zu matrix."),
7292 name, m->size1, m->size2);
7293 gsl_matrix_free (m);
7297 return matrix_to_vector (m);
7301 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7303 for (size_t i = 0; i < vars->n; i++)
7304 if (!dict_create_var (d, vars->strings[i], 0))
7305 return vars->strings[i];
7309 static struct dictionary *
7310 msave_create_dict (const struct msave_common *common)
7312 struct dictionary *dict = dict_create (get_default_encoding ());
7314 const char *dup_split = msave_add_vars (dict, &common->snames);
7317 /* Should not be possible because the parser ensures that the names are
7322 dict_create_var_assert (dict, "ROWTYPE_", 8);
7324 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7327 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
7331 dict_create_var_assert (dict, "VARNAME_", 8);
7333 const char *dup_var = msave_add_vars (dict, &common->variables);
7336 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
7348 matrix_cmd_execute_msave (struct msave_command *msave)
7350 struct msave_common *common = msave->common;
7351 gsl_matrix *m = NULL;
7352 gsl_vector *factors = NULL;
7353 gsl_vector *splits = NULL;
7355 m = matrix_expr_evaluate (msave->expr);
7359 if (!common->variables.n)
7360 for (size_t i = 0; i < m->size2; i++)
7361 string_array_append_nocopy (&common->variables,
7362 xasprintf ("COL%zu", i + 1));
7363 else if (m->size2 != common->variables.n)
7366 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7367 m->size2, common->variables.n);
7373 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7377 if (!common->fnames.n)
7378 for (size_t i = 0; i < factors->size; i++)
7379 string_array_append_nocopy (&common->fnames,
7380 xasprintf ("FAC%zu", i + 1));
7381 else if (factors->size != common->fnames.n)
7383 msg (SE, _("There are %zu factor variables, "
7384 "but %zu factor values were supplied."),
7385 common->fnames.n, factors->size);
7392 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7396 if (!common->snames.n)
7397 for (size_t i = 0; i < splits->size; i++)
7398 string_array_append_nocopy (&common->snames,
7399 xasprintf ("SPL%zu", i + 1));
7400 else if (splits->size != common->snames.n)
7402 msg (SE, _("There are %zu split variables, "
7403 "but %zu split values were supplied."),
7404 common->snames.n, splits->size);
7409 if (!common->writer)
7411 struct dictionary *dict = msave_create_dict (common);
7415 common->writer = any_writer_open (common->outfile, dict);
7416 if (!common->writer)
7422 common->dict = dict;
7425 bool matrix = (!strcmp (msave->rowtype, "COV")
7426 || !strcmp (msave->rowtype, "CORR"));
7427 for (size_t y = 0; y < m->size1; y++)
7429 struct ccase *c = case_create (dict_get_proto (common->dict));
7432 /* Split variables */
7434 for (size_t i = 0; i < splits->size; i++)
7435 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
7438 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7439 msave->rowtype, ' ');
7443 for (size_t i = 0; i < factors->size; i++)
7444 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
7447 const char *varname_ = (matrix && y < common->variables.n
7448 ? common->variables.strings[y]
7450 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7453 /* Continuous variables. */
7454 for (size_t x = 0; x < m->size2; x++)
7455 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
7456 casewriter_write (common->writer, c);
7460 gsl_matrix_free (m);
7461 gsl_vector_free (factors);
7462 gsl_vector_free (splits);
7465 static struct matrix_cmd *
7466 matrix_parse_mget (struct matrix_state *s)
7468 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7469 *cmd = (struct matrix_cmd) {
7473 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
7477 struct mget_command *mget = &cmd->mget;
7479 lex_match (s->lexer, T_SLASH);
7480 while (lex_token (s->lexer) != T_ENDCMD)
7482 if (lex_match_id (s->lexer, "FILE"))
7484 lex_match (s->lexer, T_EQUALS);
7486 fh_unref (mget->file);
7487 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7491 else if (lex_match_id (s->lexer, "ENCODING"))
7493 lex_match (s->lexer, T_EQUALS);
7494 if (!lex_force_string (s->lexer))
7497 free (mget->encoding);
7498 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
7502 else if (lex_match_id (s->lexer, "TYPE"))
7504 lex_match (s->lexer, T_EQUALS);
7505 while (lex_token (s->lexer) != T_SLASH
7506 && lex_token (s->lexer) != T_ENDCMD)
7508 const char *rowtype = match_rowtype (s->lexer);
7512 stringi_set_insert (&mget->rowtypes, rowtype);
7517 lex_error_expecting (s->lexer, "FILE", "TYPE");
7520 lex_match (s->lexer, T_SLASH);
7525 matrix_cmd_destroy (cmd);
7529 static const struct variable *
7530 get_a8_var (const struct dictionary *d, const char *name)
7532 const struct variable *v = dict_lookup_var (d, name);
7535 msg (SE, _("Matrix data file lacks %s variable."), name);
7538 if (var_get_width (v) != 8)
7540 msg (SE, _("%s variable in matrix data file must be "
7541 "8-byte string, but it has width %d."),
7542 name, var_get_width (v));
7549 var_changed (const struct ccase *ca, const struct ccase *cb,
7550 const struct variable *var)
7553 ? !value_equal (case_data (ca, var), case_data (cb, var),
7554 var_get_width (var))
7559 vars_changed (const struct ccase *ca, const struct ccase *cb,
7560 const struct dictionary *d,
7561 size_t first_var, size_t n_vars)
7563 for (size_t i = 0; i < n_vars; i++)
7565 const struct variable *v = dict_get_var (d, first_var + i);
7566 if (var_changed (ca, cb, v))
7573 vars_all_missing (const struct ccase *c, const struct dictionary *d,
7574 size_t first_var, size_t n_vars)
7576 for (size_t i = 0; i < n_vars; i++)
7577 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
7583 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
7584 const struct dictionary *d,
7585 const struct variable *rowtype_var,
7586 const struct stringi_set *accepted_rowtypes,
7587 struct matrix_state *s,
7588 size_t ss, size_t sn, size_t si,
7589 size_t fs, size_t fn, size_t fi,
7590 size_t cs, size_t cn,
7591 struct pivot_table *pt,
7592 struct pivot_dimension *var_dimension)
7597 /* Is this a matrix for pooled data, either where there are no factor
7598 variables or the factor variables are missing? */
7599 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
7601 struct substring rowtype = case_ss (rows[0], rowtype_var);
7602 ss_rtrim (&rowtype, ss_cstr (" "));
7603 if (!stringi_set_is_empty (accepted_rowtypes)
7604 && !stringi_set_contains_len (accepted_rowtypes,
7605 rowtype.string, rowtype.length))
7608 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
7609 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
7610 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
7611 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
7612 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
7613 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
7617 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
7618 (int) rowtype.length, rowtype.string);
7622 struct string name = DS_EMPTY_INITIALIZER;
7623 ds_put_cstr (&name, prefix);
7625 ds_put_format (&name, "F%zu", fi);
7627 ds_put_format (&name, "S%zu", si);
7629 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
7631 mv = matrix_var_create (s, ds_ss (&name));
7634 msg (SW, _("Matrix data file contains variable with existing name %s."),
7636 goto exit_free_name;
7639 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
7640 size_t n_missing = 0;
7641 for (size_t y = 0; y < n_rows; y++)
7643 for (size_t x = 0; x < cn; x++)
7645 struct variable *var = dict_get_var (d, cs + x);
7646 double value = case_num (rows[y], var);
7647 if (var_is_num_missing (var, value, MV_ANY))
7652 gsl_matrix_set (m, y, x, value);
7656 int var_index = pivot_category_create_leaf (
7657 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
7658 double values[] = { n_rows, cn };
7659 for (size_t j = 0; j < sn; j++)
7661 struct variable *var = dict_get_var (d, ss + j);
7662 const union value *value = case_data (rows[0], var);
7663 pivot_table_put2 (pt, j, var_index,
7664 pivot_value_new_var_value (var, value));
7666 for (size_t j = 0; j < fn; j++)
7668 struct variable *var = dict_get_var (d, fs + j);
7669 const union value sysmis = { .f = SYSMIS };
7670 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
7671 pivot_table_put2 (pt, j + sn, var_index,
7672 pivot_value_new_var_value (var, value));
7674 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
7675 pivot_table_put2 (pt, j + sn + fn, var_index,
7676 pivot_value_new_integer (values[j]));
7679 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
7680 "value, which was treated as zero.",
7681 "Matrix data file variable %s contains %zu missing "
7682 "values, which were treated as zero.", n_missing),
7683 ds_cstr (&name), n_missing);
7690 for (size_t y = 0; y < n_rows; y++)
7691 case_unref (rows[y]);
7695 matrix_cmd_execute_mget__ (struct mget_command *mget,
7696 struct casereader *r, const struct dictionary *d)
7698 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
7699 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
7700 if (!rowtype_ || !varname_)
7703 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
7705 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
7708 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
7710 msg (SE, _("Matrix data file contains no continuous variables."));
7714 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
7716 const struct variable *v = dict_get_var (d, i);
7717 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
7720 _("Matrix data file contains unexpected string variable %s."),
7726 /* SPLIT variables. */
7728 size_t sn = var_get_dict_index (rowtype_);
7729 struct ccase *sc = NULL;
7732 /* FACTOR variables. */
7733 size_t fs = var_get_dict_index (rowtype_) + 1;
7734 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
7735 struct ccase *fc = NULL;
7738 /* Continuous variables. */
7739 size_t cs = var_get_dict_index (varname_) + 1;
7740 size_t cn = dict_get_var_cnt (d) - cs;
7741 struct ccase *cc = NULL;
7744 struct pivot_table *pt = pivot_table_create (
7745 N_("Matrix Variables Created by MGET"));
7746 struct pivot_dimension *attr_dimension = pivot_dimension_create (
7747 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
7748 struct pivot_dimension *var_dimension = pivot_dimension_create (
7749 pt, PIVOT_AXIS_ROW, N_("Variable"));
7752 struct pivot_category *splits = pivot_category_create_group (
7753 attr_dimension->root, N_("Split Values"));
7754 for (size_t i = 0; i < sn; i++)
7755 pivot_category_create_leaf (splits, pivot_value_new_variable (
7756 dict_get_var (d, ss + i)));
7760 struct pivot_category *factors = pivot_category_create_group (
7761 attr_dimension->root, N_("Factors"));
7762 for (size_t i = 0; i < fn; i++)
7763 pivot_category_create_leaf (factors, pivot_value_new_variable (
7764 dict_get_var (d, fs + i)));
7766 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
7767 N_("Rows"), N_("Columns"));
7770 struct ccase **rows = NULL;
7771 size_t allocated_rows = 0;
7775 while ((c = casereader_read (r)) != NULL)
7777 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
7787 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
7788 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
7789 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
7792 if (change != NOTHING_CHANGED)
7794 matrix_mget_commit_var (rows, n_rows, d,
7795 rowtype_, &mget->rowtypes,
7806 if (n_rows >= allocated_rows)
7807 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
7810 if (change == SPLITS_CHANGED)
7816 /* Reset the factor number, if there are factors. */
7820 if (row_has_factors)
7826 else if (change == FACTORS_CHANGED)
7828 if (row_has_factors)
7834 matrix_mget_commit_var (rows, n_rows, d,
7835 rowtype_, &mget->rowtypes,
7847 if (var_dimension->n_leaves)
7848 pivot_table_submit (pt);
7850 pivot_table_unref (pt);
7854 matrix_cmd_execute_mget (struct mget_command *mget)
7856 struct casereader *r;
7857 struct dictionary *d;
7858 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
7859 mget->state->dataset, &r, &d))
7861 matrix_cmd_execute_mget__ (mget, r, d);
7862 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
7867 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
7869 if (!lex_force_id (s->lexer))
7872 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
7874 *varp = matrix_var_create (s, lex_tokss (s->lexer));
7879 static struct matrix_cmd *
7880 matrix_parse_eigen (struct matrix_state *s)
7882 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7883 *cmd = (struct matrix_cmd) {
7885 .eigen = { .expr = NULL }
7888 struct eigen_command *eigen = &cmd->eigen;
7889 if (!lex_force_match (s->lexer, T_LPAREN))
7891 eigen->expr = matrix_parse_expr (s);
7893 || !lex_force_match (s->lexer, T_COMMA)
7894 || !matrix_parse_dst_var (s, &eigen->evec)
7895 || !lex_force_match (s->lexer, T_COMMA)
7896 || !matrix_parse_dst_var (s, &eigen->eval)
7897 || !lex_force_match (s->lexer, T_RPAREN))
7903 matrix_cmd_destroy (cmd);
7908 matrix_cmd_execute_eigen (struct eigen_command *eigen)
7910 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
7913 if (!is_symmetric (A))
7915 msg (SE, _("Argument of EIGEN must be symmetric."));
7916 gsl_matrix_free (A);
7920 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
7921 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
7922 gsl_vector v_eval = to_vector (eval);
7923 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
7924 gsl_eigen_symmv (A, &v_eval, evec, w);
7925 gsl_eigen_symmv_free (w);
7927 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
7929 gsl_matrix_free (A);
7931 gsl_matrix_free (eigen->eval->value);
7932 eigen->eval->value = eval;
7934 gsl_matrix_free (eigen->evec->value);
7935 eigen->evec->value = evec;
7938 static struct matrix_cmd *
7939 matrix_parse_setdiag (struct matrix_state *s)
7941 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7942 *cmd = (struct matrix_cmd) {
7943 .type = MCMD_SETDIAG,
7944 .setdiag = { .dst = NULL }
7947 struct setdiag_command *setdiag = &cmd->setdiag;
7948 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
7951 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
7954 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
7959 if (!lex_force_match (s->lexer, T_COMMA))
7962 setdiag->expr = matrix_parse_expr (s);
7966 if (!lex_force_match (s->lexer, T_RPAREN))
7972 matrix_cmd_destroy (cmd);
7977 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
7979 gsl_matrix *dst = setdiag->dst->value;
7982 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
7983 setdiag->dst->name);
7987 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
7991 size_t n = MIN (dst->size1, dst->size2);
7992 if (is_scalar (src))
7994 double d = to_scalar (src);
7995 for (size_t i = 0; i < n; i++)
7996 gsl_matrix_set (dst, i, i, d);
7998 else if (is_vector (src))
8000 gsl_vector v = to_vector (src);
8001 for (size_t i = 0; i < n && i < v.size; i++)
8002 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
8005 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
8006 "not a %zu×%zu matrix."),
8007 src->size1, src->size2);
8008 gsl_matrix_free (src);
8011 static struct matrix_cmd *
8012 matrix_parse_svd (struct matrix_state *s)
8014 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
8015 *cmd = (struct matrix_cmd) {
8017 .svd = { .expr = NULL }
8020 struct svd_command *svd = &cmd->svd;
8021 if (!lex_force_match (s->lexer, T_LPAREN))
8023 svd->expr = matrix_parse_expr (s);
8025 || !lex_force_match (s->lexer, T_COMMA)
8026 || !matrix_parse_dst_var (s, &svd->u)
8027 || !lex_force_match (s->lexer, T_COMMA)
8028 || !matrix_parse_dst_var (s, &svd->s)
8029 || !lex_force_match (s->lexer, T_COMMA)
8030 || !matrix_parse_dst_var (s, &svd->v)
8031 || !lex_force_match (s->lexer, T_RPAREN))
8037 matrix_cmd_destroy (cmd);
8042 matrix_cmd_execute_svd (struct svd_command *svd)
8044 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8048 if (m->size1 >= m->size2)
8051 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8052 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8053 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8054 gsl_vector *work = gsl_vector_alloc (A->size2);
8055 gsl_linalg_SV_decomp (A, V, &Sv, work);
8056 gsl_vector_free (work);
8058 matrix_var_set (svd->u, A);
8059 matrix_var_set (svd->s, S);
8060 matrix_var_set (svd->v, V);
8064 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8065 gsl_matrix_transpose_memcpy (At, m);
8066 gsl_matrix_free (m);
8068 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8069 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8070 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8071 gsl_vector *work = gsl_vector_alloc (At->size2);
8072 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8073 gsl_vector_free (work);
8075 matrix_var_set (svd->v, At);
8076 matrix_var_set (svd->s, St);
8077 matrix_var_set (svd->u, Vt);
8082 matrix_cmd_execute (struct matrix_cmd *cmd)
8087 matrix_cmd_execute_compute (&cmd->compute);
8091 matrix_cmd_execute_print (&cmd->print);
8095 return matrix_cmd_execute_do_if (&cmd->do_if);
8098 matrix_cmd_execute_loop (&cmd->loop);
8105 matrix_cmd_execute_display (&cmd->display);
8109 matrix_cmd_execute_release (&cmd->release);
8113 matrix_cmd_execute_save (&cmd->save);
8117 matrix_cmd_execute_read (&cmd->read);
8121 matrix_cmd_execute_write (&cmd->write);
8125 matrix_cmd_execute_get (&cmd->get);
8129 matrix_cmd_execute_msave (&cmd->msave);
8133 matrix_cmd_execute_mget (&cmd->mget);
8137 matrix_cmd_execute_eigen (&cmd->eigen);
8141 matrix_cmd_execute_setdiag (&cmd->setdiag);
8145 matrix_cmd_execute_svd (&cmd->svd);
8153 matrix_cmds_uninit (struct matrix_cmds *cmds)
8155 for (size_t i = 0; i < cmds->n; i++)
8156 matrix_cmd_destroy (cmds->commands[i]);
8157 free (cmds->commands);
8161 matrix_cmd_destroy (struct matrix_cmd *cmd)
8169 matrix_lvalue_destroy (cmd->compute.lvalue);
8170 matrix_expr_destroy (cmd->compute.rvalue);
8174 matrix_expr_destroy (cmd->print.expression);
8175 free (cmd->print.title);
8176 print_labels_destroy (cmd->print.rlabels);
8177 print_labels_destroy (cmd->print.clabels);
8181 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8183 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8184 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
8186 free (cmd->do_if.clauses);
8190 matrix_expr_destroy (cmd->loop.start);
8191 matrix_expr_destroy (cmd->loop.end);
8192 matrix_expr_destroy (cmd->loop.increment);
8193 matrix_expr_destroy (cmd->loop.top_condition);
8194 matrix_expr_destroy (cmd->loop.bottom_condition);
8195 matrix_cmds_uninit (&cmd->loop.commands);
8205 free (cmd->release.vars);
8209 matrix_expr_destroy (cmd->save.expression);
8213 matrix_lvalue_destroy (cmd->read.dst);
8214 matrix_expr_destroy (cmd->read.size);
8218 matrix_expr_destroy (cmd->write.expression);
8219 free (cmd->write.format);
8223 matrix_lvalue_destroy (cmd->get.dst);
8224 fh_unref (cmd->get.file);
8225 free (cmd->get.encoding);
8226 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8230 matrix_expr_destroy (cmd->msave.expr);
8234 fh_unref (cmd->mget.file);
8235 stringi_set_destroy (&cmd->mget.rowtypes);
8239 matrix_expr_destroy (cmd->eigen.expr);
8243 matrix_expr_destroy (cmd->setdiag.expr);
8247 matrix_expr_destroy (cmd->svd.expr);
8253 struct matrix_command_name
8256 struct matrix_cmd *(*parse) (struct matrix_state *);
8259 static const struct matrix_command_name *
8260 matrix_parse_command_name (struct lexer *lexer)
8262 static const struct matrix_command_name commands[] = {
8263 { "COMPUTE", matrix_parse_compute },
8264 { "CALL EIGEN", matrix_parse_eigen },
8265 { "CALL SETDIAG", matrix_parse_setdiag },
8266 { "CALL SVD", matrix_parse_svd },
8267 { "PRINT", matrix_parse_print },
8268 { "DO IF", matrix_parse_do_if },
8269 { "LOOP", matrix_parse_loop },
8270 { "BREAK", matrix_parse_break },
8271 { "READ", matrix_parse_read },
8272 { "WRITE", matrix_parse_write },
8273 { "GET", matrix_parse_get },
8274 { "SAVE", matrix_parse_save },
8275 { "MGET", matrix_parse_mget },
8276 { "MSAVE", matrix_parse_msave },
8277 { "DISPLAY", matrix_parse_display },
8278 { "RELEASE", matrix_parse_release },
8280 static size_t n = sizeof commands / sizeof *commands;
8282 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8283 if (lex_match_phrase (lexer, c->name))
8288 static struct matrix_cmd *
8289 matrix_parse_command (struct matrix_state *s)
8291 size_t nesting_level = SIZE_MAX;
8293 struct matrix_cmd *c = NULL;
8294 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
8296 lex_error (s->lexer, _("Unknown matrix command."));
8297 else if (!cmd->parse)
8298 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8302 nesting_level = output_open_group (
8303 group_item_create_nocopy (utf8_to_title (cmd->name),
8304 utf8_to_title (cmd->name)));
8309 lex_end_of_command (s->lexer);
8310 lex_discard_rest_of_command (s->lexer);
8311 if (nesting_level != SIZE_MAX)
8312 output_close_groups (nesting_level);
8318 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8320 if (!lex_force_match (lexer, T_ENDCMD))
8323 struct matrix_state state = {
8325 .session = dataset_session (ds),
8327 .vars = HMAP_INITIALIZER (state.vars),
8332 while (lex_match (lexer, T_ENDCMD))
8334 if (lex_token (lexer) == T_STOP)
8336 msg (SE, _("Unexpected end of input expecting matrix command."));
8340 if (lex_match_phrase (lexer, "END MATRIX"))
8343 struct matrix_cmd *c = matrix_parse_command (&state);
8346 matrix_cmd_execute (c);
8347 matrix_cmd_destroy (c);
8351 struct matrix_var *var, *next;
8352 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8355 gsl_matrix_free (var->value);
8356 hmap_delete (&state.vars, &var->hmap_node);
8359 hmap_destroy (&state.vars);
8360 msave_common_destroy (state.common);
8361 fh_unref (state.prev_read_file);
8362 for (size_t i = 0; i < state.n_read_files; i++)
8363 read_file_destroy (state.read_files[i]);
8364 free (state.read_files);
8365 fh_unref (state.prev_write_file);
8366 for (size_t i = 0; i < state.n_write_files; i++)
8367 write_file_destroy (state.write_files[i]);
8368 free (state.write_files);
8369 fh_unref (state.prev_save_file);
8370 for (size_t i = 0; i < state.n_save_files; i++)
8371 save_file_destroy (state.save_files[i]);
8372 free (state.save_files);