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;
507 matrix_expr_destroy (struct matrix_expr *e)
514 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
545 for (size_t i = 0; i < e->n_subs; i++)
546 matrix_expr_destroy (e->subs[i]);
558 static struct matrix_expr *
559 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
562 struct matrix_expr *e = xmalloc (sizeof *e);
563 *e = (struct matrix_expr) {
565 .subs = xmemdup (subs, n_subs * sizeof *subs),
571 static struct matrix_expr *
572 matrix_expr_create_0 (enum matrix_op op)
574 struct matrix_expr *sub;
575 return matrix_expr_create_subs (op, &sub, 0);
578 static struct matrix_expr *
579 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
581 return matrix_expr_create_subs (op, &sub, 1);
584 static struct matrix_expr *
585 matrix_expr_create_2 (enum matrix_op op,
586 struct matrix_expr *sub0, struct matrix_expr *sub1)
588 struct matrix_expr *subs[] = { sub0, sub1 };
589 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
592 static struct matrix_expr *
593 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
594 struct matrix_expr *sub1, struct matrix_expr *sub2)
596 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
597 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
600 static struct matrix_expr *
601 matrix_expr_create_number (double number)
603 struct matrix_expr *e = xmalloc (sizeof *e);
604 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
608 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
610 static struct matrix_expr *
611 matrix_parse_curly_comma (struct matrix_state *s)
613 struct matrix_expr *lhs = matrix_parse_or_xor (s);
617 while (lex_match (s->lexer, T_COMMA))
619 struct matrix_expr *rhs = matrix_parse_or_xor (s);
622 matrix_expr_destroy (lhs);
625 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
630 static struct matrix_expr *
631 matrix_parse_curly_semi (struct matrix_state *s)
633 if (lex_token (s->lexer) == T_RCURLY)
634 return matrix_expr_create_0 (MOP_EMPTY);
636 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
640 while (lex_match (s->lexer, T_SEMICOLON))
642 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
645 matrix_expr_destroy (lhs);
648 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
653 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
654 for (size_t Y = 0; Y < (M)->size1; Y++) \
655 for (size_t X = 0; X < (M)->size2; X++) \
656 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
659 is_vector (const gsl_matrix *m)
661 return m->size1 <= 1 || m->size2 <= 1;
665 to_vector (gsl_matrix *m)
667 return (m->size1 == 1
668 ? gsl_matrix_row (m, 0).vector
669 : gsl_matrix_column (m, 0).vector);
674 matrix_eval_ABS (double d)
680 matrix_eval_ALL (gsl_matrix *m)
682 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
689 matrix_eval_ANY (gsl_matrix *m)
691 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
698 matrix_eval_ARSIN (double d)
704 matrix_eval_ARTAN (double d)
710 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
714 for (size_t i = 0; i < n; i++)
719 gsl_matrix *block = gsl_matrix_calloc (r, c);
721 for (size_t i = 0; i < n; i++)
723 for (size_t y = 0; y < m[i]->size1; y++)
724 for (size_t x = 0; x < m[i]->size2; x++)
725 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
733 matrix_eval_CHOL (gsl_matrix *m)
735 if (!gsl_linalg_cholesky_decomp1 (m))
737 for (size_t y = 0; y < m->size1; y++)
738 for (size_t x = y + 1; x < m->size2; x++)
739 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
741 for (size_t y = 0; y < m->size1; y++)
742 for (size_t x = 0; x < y; x++)
743 gsl_matrix_set (m, y, x, 0);
748 msg (SE, _("Input to CHOL function is not positive-definite."));
754 matrix_eval_col_extremum (gsl_matrix *m, bool min)
759 return gsl_matrix_alloc (1, 0);
761 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
762 for (size_t x = 0; x < m->size2; x++)
764 double ext = gsl_matrix_get (m, 0, x);
765 for (size_t y = 1; y < m->size1; y++)
767 double value = gsl_matrix_get (m, y, x);
768 if (min ? value < ext : value > ext)
771 gsl_matrix_set (cext, 0, x, ext);
777 matrix_eval_CMAX (gsl_matrix *m)
779 return matrix_eval_col_extremum (m, false);
783 matrix_eval_CMIN (gsl_matrix *m)
785 return matrix_eval_col_extremum (m, true);
789 matrix_eval_COS (double d)
795 matrix_eval_col_sum (gsl_matrix *m, bool square)
800 return gsl_matrix_alloc (1, 0);
802 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
803 for (size_t x = 0; x < m->size2; x++)
806 for (size_t y = 0; y < m->size1; y++)
808 double d = gsl_matrix_get (m, y, x);
809 sum += square ? pow2 (d) : d;
811 gsl_matrix_set (result, 0, x, sum);
817 matrix_eval_CSSQ (gsl_matrix *m)
819 return matrix_eval_col_sum (m, true);
823 matrix_eval_CSUM (gsl_matrix *m)
825 return matrix_eval_col_sum (m, false);
829 compare_double_3way (const void *a_, const void *b_)
831 const double *a = a_;
832 const double *b = b_;
833 return *a < *b ? -1 : *a > *b;
837 matrix_eval_DESIGN (gsl_matrix *m)
839 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
840 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
841 gsl_matrix_transpose_memcpy (&m2, m);
843 for (size_t y = 0; y < m2.size1; y++)
844 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
846 size_t *n = xcalloc (m2.size1, sizeof *n);
848 for (size_t i = 0; i < m2.size1; i++)
850 double *row = tmp + m2.size2 * i;
851 for (size_t j = 0; j < m2.size2; )
854 for (k = j + 1; k < m2.size2; k++)
855 if (row[j] != row[k])
857 row[n[i]++] = row[j];
862 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
867 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
869 for (size_t i = 0; i < m->size2; i++)
874 const double *unique = tmp + m2.size2 * i;
875 for (size_t j = 0; j < n[i]; j++, x++)
877 double value = unique[j];
878 for (size_t y = 0; y < m->size1; y++)
879 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
890 matrix_eval_DET (gsl_matrix *m)
892 gsl_permutation *p = gsl_permutation_alloc (m->size1);
894 gsl_linalg_LU_decomp (m, p, &signum);
895 gsl_permutation_free (p);
896 return gsl_linalg_LU_det (m, signum);
900 matrix_eval_DIAG (gsl_matrix *m)
902 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
903 for (size_t i = 0; i < diag->size1; i++)
904 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
909 is_symmetric (const gsl_matrix *m)
911 if (m->size1 != m->size2)
914 for (size_t y = 0; y < m->size1; y++)
915 for (size_t x = 0; x < y; x++)
916 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
923 compare_double_desc (const void *a_, const void *b_)
925 const double *a = a_;
926 const double *b = b_;
927 return *a > *b ? -1 : *a < *b;
931 matrix_eval_EVAL (gsl_matrix *m)
933 if (!is_symmetric (m))
935 msg (SE, _("Argument of EVAL must be symmetric."));
939 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
940 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
941 gsl_vector v_eval = to_vector (eval);
942 gsl_eigen_symm (m, &v_eval, w);
943 gsl_eigen_symm_free (w);
945 assert (v_eval.stride == 1);
946 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
952 matrix_eval_EXP (double d)
957 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
960 Charl Linssen <charl@itfromb.it>
964 matrix_eval_GINV (gsl_matrix *A)
969 gsl_matrix *tmp_mat = NULL;
972 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
973 tmp_mat = gsl_matrix_alloc (m, n);
974 gsl_matrix_transpose_memcpy (tmp_mat, A);
982 gsl_matrix *V = gsl_matrix_alloc (m, m);
983 gsl_vector *u = gsl_vector_alloc (m);
985 gsl_vector *tmp_vec = gsl_vector_alloc (m);
986 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
987 gsl_vector_free (tmp_vec);
990 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
991 gsl_matrix_set_zero (Sigma_pinv);
992 double cutoff = 1e-15 * gsl_vector_max (u);
994 for (size_t i = 0; i < m; ++i)
996 double x = gsl_vector_get (u, i);
997 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
1000 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
1001 gsl_matrix *U = gsl_matrix_calloc (n, n);
1002 for (size_t i = 0; i < n; i++)
1003 for (size_t j = 0; j < m; j++)
1004 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
1006 /* two dot products to obtain pseudoinverse */
1007 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
1008 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
1013 A_pinv = gsl_matrix_alloc (n, m);
1014 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
1018 A_pinv = gsl_matrix_alloc (m, n);
1019 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
1022 gsl_matrix_free (tmp_mat);
1023 gsl_matrix_free (tmp_mat2);
1024 gsl_matrix_free (U);
1025 gsl_matrix_free (Sigma_pinv);
1026 gsl_vector_free (u);
1027 gsl_matrix_free (V);
1039 grade_compare_3way (const void *a_, const void *b_)
1041 const struct grade *a = a_;
1042 const struct grade *b = b_;
1044 return (a->value < b->value ? -1
1045 : a->value > b->value ? 1
1053 matrix_eval_GRADE (gsl_matrix *m)
1055 size_t n = m->size1 * m->size2;
1056 struct grade *grades = xmalloc (n * sizeof *grades);
1059 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1060 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
1061 qsort (grades, n, sizeof *grades, grade_compare_3way);
1063 for (size_t i = 0; i < n; i++)
1064 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
1072 dot (gsl_vector *a, gsl_vector *b)
1074 double result = 0.0;
1075 for (size_t i = 0; i < a->size; i++)
1076 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
1081 norm2 (gsl_vector *v)
1083 double result = 0.0;
1084 for (size_t i = 0; i < v->size; i++)
1085 result += pow2 (gsl_vector_get (v, i));
1090 norm (gsl_vector *v)
1092 return sqrt (norm2 (v));
1096 matrix_eval_GSCH (gsl_matrix *v)
1098 if (v->size2 < v->size1)
1100 msg (SE, _("GSCH requires its argument to have at least as many columns "
1101 "as rows, but it has dimensions %zu×%zu."),
1102 v->size1, v->size2);
1105 if (!v->size1 || !v->size2)
1108 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
1110 for (size_t vx = 0; vx < v->size2; vx++)
1112 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
1113 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
1115 gsl_vector_memcpy (&u_i, &v_i);
1116 for (size_t j = 0; j < ux; j++)
1118 gsl_vector u_j = gsl_matrix_column (u, j).vector;
1119 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
1120 for (size_t k = 0; k < u_i.size; k++)
1121 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
1122 - scale * gsl_vector_get (&u_j, k)));
1125 double len = norm (&u_i);
1128 gsl_vector_scale (&u_i, 1.0 / len);
1129 if (++ux >= v->size1)
1136 msg (SE, _("%zu×%zu argument to GSCH contains only "
1137 "%zu linearly independent columns."),
1138 v->size1, v->size2, ux);
1139 gsl_matrix_free (u);
1143 u->size2 = v->size1;
1148 matrix_eval_IDENT (double s1, double s2)
1150 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
1152 msg (SE, _("Arguments to IDENT must be non-negative integers."));
1156 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
1157 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1163 invert_matrix (gsl_matrix *x)
1165 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1167 gsl_linalg_LU_decomp (x, p, &signum);
1168 gsl_linalg_LU_invx (x, p);
1169 gsl_permutation_free (p);
1173 matrix_eval_INV (gsl_matrix *m)
1180 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1182 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1183 a->size2 * b->size2);
1185 for (size_t ar = 0; ar < a->size1; ar++)
1186 for (size_t br = 0; br < b->size1; br++, y++)
1189 for (size_t ac = 0; ac < a->size2; ac++)
1190 for (size_t bc = 0; bc < b->size2; bc++, x++)
1192 double av = gsl_matrix_get (a, ar, ac);
1193 double bv = gsl_matrix_get (b, br, bc);
1194 gsl_matrix_set (k, y, x, av * bv);
1201 matrix_eval_LG10 (double d)
1207 matrix_eval_LN (double d)
1213 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1215 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1218 for (size_t i = 1; i <= n * n; i++)
1220 gsl_matrix_set (m, y, x, i);
1222 size_t y1 = !y ? n - 1 : y - 1;
1223 size_t x1 = x + 1 >= n ? 0 : x + 1;
1224 if (gsl_matrix_get (m, y1, x1) == 0)
1230 y = y + 1 >= n ? 0 : y + 1;
1235 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1237 double a = gsl_matrix_get (m, y1, x1);
1238 double b = gsl_matrix_get (m, y2, x2);
1239 gsl_matrix_set (m, y1, x1, b);
1240 gsl_matrix_set (m, y2, x2, a);
1244 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1248 /* A. Umar, "On the Construction of Even Order Magic Squares",
1249 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1251 for (size_t i = 1; i <= n * n / 2; i++)
1253 gsl_matrix_set (m, y, x, i);
1263 for (size_t i = n * n; i > n * n / 2; i--)
1265 gsl_matrix_set (m, y, x, i);
1273 for (size_t y = 0; y < n; y++)
1274 for (size_t x = 0; x < n / 2; x++)
1276 unsigned int d = gsl_matrix_get (m, y, x);
1277 if (d % 2 != (y < n / 2))
1278 magic_exchange (m, y, x, y, n - x - 1);
1283 size_t x1 = n / 2 - 1;
1285 magic_exchange (m, y1, x1, y2, x1);
1286 magic_exchange (m, y1, x2, y2, x2);
1290 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1292 /* A. Umar, "On the Construction of Even Order Magic Squares",
1293 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1297 for (size_t i = 1; ; i++)
1299 gsl_matrix_set (m, y, x, i);
1300 if (++y == n / 2 - 1)
1312 for (size_t i = n * n; ; i--)
1314 gsl_matrix_set (m, y, x, i);
1315 if (++y == n / 2 - 1)
1324 for (size_t y = 0; y < n; y++)
1325 if (y != n / 2 - 1 && y != n / 2)
1326 for (size_t x = 0; x < n / 2; x++)
1328 unsigned int d = gsl_matrix_get (m, y, x);
1329 if (d % 2 != (y < n / 2))
1330 magic_exchange (m, y, x, y, n - x - 1);
1333 size_t a0 = (n * n - 2 * n) / 2 + 1;
1334 for (size_t i = 0; i < n / 2; i++)
1337 gsl_matrix_set (m, n / 2, i, a);
1338 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1340 for (size_t i = 0; i < n / 2; i++)
1342 size_t a = a0 + i + n / 2;
1343 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1344 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1346 for (size_t x = 1; x < n / 2; x += 2)
1347 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1348 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1349 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1350 size_t x1 = n / 2 - 2;
1351 size_t x2 = n / 2 + 1;
1352 size_t y1 = n / 2 - 2;
1353 size_t y2 = n / 2 + 1;
1354 magic_exchange (m, y1, x1, y2, x1);
1355 magic_exchange (m, y1, x2, y2, x2);
1359 matrix_eval_MAGIC (double n_)
1363 gsl_matrix *m = gsl_matrix_calloc (n, n);
1365 matrix_eval_MAGIC_odd (m, n);
1367 matrix_eval_MAGIC_singly_even (m, n);
1369 matrix_eval_MAGIC_doubly_even (m, n);
1374 matrix_eval_MAKE (double r, double c, double value)
1376 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1378 msg (SE, _("Size arguments to MAKE must be integers."));
1382 gsl_matrix *m = gsl_matrix_alloc (r, c);
1383 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1389 matrix_eval_MDIAG (gsl_vector *v)
1391 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1392 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1393 gsl_vector_memcpy (&diagonal, v);
1398 matrix_eval_MMAX (gsl_matrix *m)
1400 return gsl_matrix_max (m);
1404 matrix_eval_MMIN (gsl_matrix *m)
1406 return gsl_matrix_min (m);
1410 matrix_eval_MOD (gsl_matrix *m, double divisor)
1414 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1418 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1419 *d = fmod (*d, divisor);
1424 matrix_eval_MSSQ (gsl_matrix *m)
1427 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1433 matrix_eval_MSUM (gsl_matrix *m)
1436 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1442 matrix_eval_NCOL (gsl_matrix *m)
1448 matrix_eval_NROW (gsl_matrix *m)
1454 matrix_eval_RANK (gsl_matrix *m)
1456 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1457 gsl_linalg_QR_decomp (m, tau);
1458 gsl_vector_free (tau);
1460 return gsl_linalg_QRPT_rank (m, -1);
1464 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1466 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1468 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1473 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1475 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1476 "product of matrix dimensions."));
1480 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1483 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1485 gsl_matrix_set (dst, y1, x1, *d);
1496 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1501 return gsl_matrix_alloc (0, 1);
1503 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1504 for (size_t y = 0; y < m->size1; y++)
1506 double ext = gsl_matrix_get (m, y, 0);
1507 for (size_t x = 1; x < m->size2; x++)
1509 double value = gsl_matrix_get (m, y, x);
1510 if (min ? value < ext : value > ext)
1513 gsl_matrix_set (rext, y, 0, ext);
1519 matrix_eval_RMAX (gsl_matrix *m)
1521 return matrix_eval_row_extremum (m, false);
1525 matrix_eval_RMIN (gsl_matrix *m)
1527 return matrix_eval_row_extremum (m, true);
1531 matrix_eval_RND (double d)
1543 rank_compare_3way (const void *a_, const void *b_)
1545 const struct rank *a = a_;
1546 const struct rank *b = b_;
1548 return a->value < b->value ? -1 : a->value > b->value;
1552 matrix_eval_RNKORDER (gsl_matrix *m)
1554 size_t n = m->size1 * m->size2;
1555 struct rank *ranks = xmalloc (n * sizeof *ranks);
1557 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1558 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1559 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1561 for (size_t i = 0; i < n; )
1564 for (j = i + 1; j < n; j++)
1565 if (ranks[i].value != ranks[j].value)
1568 double rank = (i + j + 1.0) / 2.0;
1569 for (size_t k = i; k < j; k++)
1570 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1581 matrix_eval_row_sum (gsl_matrix *m, bool square)
1586 return gsl_matrix_alloc (0, 1);
1588 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1589 for (size_t y = 0; y < m->size1; y++)
1592 for (size_t x = 0; x < m->size2; x++)
1594 double d = gsl_matrix_get (m, y, x);
1595 sum += square ? pow2 (d) : d;
1597 gsl_matrix_set (result, y, 0, sum);
1603 matrix_eval_RSSQ (gsl_matrix *m)
1605 return matrix_eval_row_sum (m, true);
1609 matrix_eval_RSUM (gsl_matrix *m)
1611 return matrix_eval_row_sum (m, false);
1615 matrix_eval_SIN (double d)
1621 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1623 if (m1->size1 != m2->size1)
1625 msg (SE, _("SOLVE requires its arguments to have the same number of "
1626 "rows, but the first argument has dimensions %zu×%zu and "
1627 "the second %zu×%zu."),
1628 m1->size1, m1->size2,
1629 m2->size1, m2->size2);
1633 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1634 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1636 gsl_linalg_LU_decomp (m1, p, &signum);
1637 for (size_t i = 0; i < m2->size2; i++)
1639 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1640 gsl_vector xi = gsl_matrix_column (x, i).vector;
1641 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1643 gsl_permutation_free (p);
1648 matrix_eval_SQRT (double d)
1654 matrix_eval_SSCP (gsl_matrix *m)
1656 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1657 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1662 matrix_eval_SVAL (gsl_matrix *m)
1664 gsl_matrix *tmp_mat = NULL;
1665 if (m->size2 > m->size1)
1667 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1668 gsl_matrix_transpose_memcpy (tmp_mat, m);
1673 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1674 gsl_vector *S = gsl_vector_alloc (m->size2);
1675 gsl_vector *work = gsl_vector_alloc (m->size2);
1676 gsl_linalg_SV_decomp (m, V, S, work);
1678 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1679 for (size_t i = 0; i < m->size2; i++)
1680 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1682 gsl_matrix_free (V);
1683 gsl_vector_free (S);
1684 gsl_vector_free (work);
1685 gsl_matrix_free (tmp_mat);
1691 matrix_eval_SWEEP (gsl_matrix *m, double d)
1693 if (d < 1 || d > SIZE_MAX)
1695 msg (SE, _("Scalar argument to SWEEP must be integer."));
1699 if (k >= MIN (m->size1, m->size2))
1701 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1702 "equal to the smaller of the matrix argument's rows and "
1707 double m_kk = gsl_matrix_get (m, k, k);
1708 if (fabs (m_kk) > 1e-19)
1710 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1711 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1713 double m_ij = gsl_matrix_get (m, i, j);
1714 double m_ik = gsl_matrix_get (m, i, k);
1715 double m_kj = gsl_matrix_get (m, k, j);
1716 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1725 for (size_t i = 0; i < m->size1; i++)
1727 gsl_matrix_set (m, i, k, 0);
1728 gsl_matrix_set (m, k, i, 0);
1735 matrix_eval_TRACE (gsl_matrix *m)
1738 size_t n = MIN (m->size1, m->size2);
1739 for (size_t i = 0; i < n; i++)
1740 sum += gsl_matrix_get (m, i, i);
1745 matrix_eval_T (gsl_matrix *m)
1747 return matrix_eval_TRANSPOS (m);
1751 matrix_eval_TRANSPOS (gsl_matrix *m)
1753 if (m->size1 == m->size2)
1755 gsl_matrix_transpose (m);
1760 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1761 gsl_matrix_transpose_memcpy (t, m);
1767 matrix_eval_TRUNC (double d)
1773 matrix_eval_UNIFORM (double r_, double c_)
1777 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1779 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1783 gsl_matrix *m = gsl_matrix_alloc (r, c);
1784 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1785 *d = gsl_ran_flat (get_rng (), 0, 1);
1790 matrix_eval_PDF_BETA (double x, double a, double b)
1792 return gsl_ran_beta_pdf (x, a, b);
1796 matrix_eval_CDF_BETA (double x, double a, double b)
1798 return gsl_cdf_beta_P (x, a, b);
1802 matrix_eval_IDF_BETA (double P, double a, double b)
1804 return gsl_cdf_beta_Pinv (P, a, b);
1808 matrix_eval_RV_BETA (double a, double b)
1810 return gsl_ran_beta (get_rng (), a, b);
1814 matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
1816 return ncdf_beta (x, a, b, lambda);
1820 matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
1822 return npdf_beta (x, a, b, lambda);
1826 matrix_eval_CDF_BVNOR (double x0, double x1, double r)
1828 return cdf_bvnor (x0, x1, r);
1832 matrix_eval_PDF_BVNOR (double x0, double x1, double r)
1834 return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
1838 matrix_eval_CDF_CAUCHY (double x, double a, double b)
1840 return gsl_cdf_cauchy_P ((x - a) / b, 1);
1844 matrix_eval_IDF_CAUCHY (double P, double a, double b)
1846 return a + b * gsl_cdf_cauchy_Pinv (P, 1);
1850 matrix_eval_PDF_CAUCHY (double x, double a, double b)
1852 return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
1856 matrix_eval_RV_CAUCHY (double a, double b)
1858 return a + b * gsl_ran_cauchy (get_rng (), 1);
1862 matrix_eval_CDF_CHISQ (double x, double df)
1864 return gsl_cdf_chisq_P (x, df);
1868 matrix_eval_CHICDF (double x, double df)
1870 return matrix_eval_CDF_CHISQ (x, df);
1874 matrix_eval_IDF_CHISQ (double P, double df)
1876 return gsl_cdf_chisq_Pinv (P, df);
1880 matrix_eval_PDF_CHISQ (double x, double df)
1882 return gsl_ran_chisq_pdf (x, df);
1886 matrix_eval_RV_CHISQ (double df)
1888 return gsl_ran_chisq (get_rng (), df);
1892 matrix_eval_SIG_CHISQ (double x, double df)
1894 return gsl_cdf_chisq_Q (x, df);
1898 matrix_eval_CDF_EXP (double x, double a)
1900 return gsl_cdf_exponential_P (x, 1. / a);
1904 matrix_eval_IDF_EXP (double P, double a)
1906 return gsl_cdf_exponential_Pinv (P, 1. / a);
1910 matrix_eval_PDF_EXP (double x, double a)
1912 return gsl_ran_exponential_pdf (x, 1. / a);
1916 matrix_eval_RV_EXP (double a)
1918 return gsl_ran_exponential (get_rng (), 1. / a);
1922 matrix_eval_PDF_XPOWER (double x, double a, double b)
1924 return gsl_ran_exppow_pdf (x, a, b);
1928 matrix_eval_RV_XPOWER (double a, double b)
1930 return gsl_ran_exppow (get_rng (), a, b);
1934 matrix_eval_CDF_F (double x, double df1, double df2)
1936 return gsl_cdf_fdist_P (x, df1, df2);
1940 matrix_eval_FCDF (double x, double df1, double df2)
1942 return matrix_eval_CDF_F (x, df1, df2);
1946 matrix_eval_IDF_F (double P, double df1, double df2)
1948 return idf_fdist (P, df1, df2);
1952 matrix_eval_RV_F (double df1, double df2)
1954 return gsl_ran_fdist (get_rng (), df1, df2);
1958 matrix_eval_PDF_F (double x, double df1, double df2)
1960 return gsl_ran_fdist_pdf (x, df1, df2);
1964 matrix_eval_SIG_F (double x, double df1, double df2)
1966 return gsl_cdf_fdist_Q (x, df1, df2);
1970 matrix_eval_CDF_GAMMA (double x, double a, double b)
1972 return gsl_cdf_gamma_P (x, a, 1. / b);
1976 matrix_eval_IDF_GAMMA (double P, double a, double b)
1978 return gsl_cdf_gamma_Pinv (P, a, 1. / b);
1982 matrix_eval_PDF_GAMMA (double x, double a, double b)
1984 return gsl_ran_gamma_pdf (x, a, 1. / b);
1988 matrix_eval_RV_GAMMA (double a, double b)
1990 return gsl_ran_gamma (get_rng (), a, 1. / b);
1994 matrix_eval_PDF_LANDAU (double x)
1996 return gsl_ran_landau_pdf (x);
2000 matrix_eval_RV_LANDAU (void)
2002 return gsl_ran_landau (get_rng ());
2006 matrix_eval_CDF_LAPLACE (double x, double a, double b)
2008 return gsl_cdf_laplace_P ((x - a) / b, 1);
2012 matrix_eval_IDF_LAPLACE (double P, double a, double b)
2014 return a + b * gsl_cdf_laplace_Pinv (P, 1);
2018 matrix_eval_PDF_LAPLACE (double x, double a, double b)
2020 return gsl_ran_laplace_pdf ((x - a) / b, 1);
2024 matrix_eval_RV_LAPLACE (double a, double b)
2026 return a + b * gsl_ran_laplace (get_rng (), 1);
2030 matrix_eval_RV_LEVY (double c, double alpha)
2032 return gsl_ran_levy (get_rng (), c, alpha);
2036 matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
2038 return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
2042 matrix_eval_CDF_LOGISTIC (double x, double a, double b)
2044 return gsl_cdf_logistic_P ((x - a) / b, 1);
2048 matrix_eval_IDF_LOGISTIC (double P, double a, double b)
2050 return a + b * gsl_cdf_logistic_Pinv (P, 1);
2054 matrix_eval_PDF_LOGISTIC (double x, double a, double b)
2056 return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
2060 matrix_eval_RV_LOGISTIC (double a, double b)
2062 return a + b * gsl_ran_logistic (get_rng (), 1);
2066 matrix_eval_CDF_LNORMAL (double x, double m, double s)
2068 return gsl_cdf_lognormal_P (x, log (m), s);
2072 matrix_eval_IDF_LNORMAL (double P, double m, double s)
2074 return gsl_cdf_lognormal_Pinv (P, log (m), s);;
2078 matrix_eval_PDF_LNORMAL (double x, double m, double s)
2080 return gsl_ran_lognormal_pdf (x, log (m), s);
2084 matrix_eval_RV_LNORMAL (double m, double s)
2086 return gsl_ran_lognormal (get_rng (), log (m), s);
2090 matrix_eval_CDF_NORMAL (double x, double u, double s)
2092 return gsl_cdf_gaussian_P (x - u, s);
2096 matrix_eval_IDF_NORMAL (double P, double u, double s)
2098 return u + gsl_cdf_gaussian_Pinv (P, s);
2102 matrix_eval_PDF_NORMAL (double x, double u, double s)
2104 return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
2108 matrix_eval_RV_NORMAL (double u, double s)
2110 return u + gsl_ran_gaussian (get_rng (), s);
2114 matrix_eval_CDFNORM (double x)
2116 return gsl_cdf_ugaussian_P (x);
2120 matrix_eval_PROBIT (double P)
2122 return gsl_cdf_ugaussian_Pinv (P);
2126 matrix_eval_NORMAL (double s)
2128 return gsl_ran_gaussian (get_rng (), s);
2132 matrix_eval_PDF_NTAIL (double x, double a, double sigma)
2134 return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
2138 matrix_eval_RV_NTAIL (double a, double sigma)
2140 return gsl_ran_gaussian_tail (get_rng (), a, sigma);
2144 matrix_eval_CDF_PARETO (double x, double a, double b)
2146 return gsl_cdf_pareto_P (x, b, a);
2150 matrix_eval_IDF_PARETO (double P, double a, double b)
2152 return gsl_cdf_pareto_Pinv (P, b, a);
2156 matrix_eval_PDF_PARETO (double x, double a, double b)
2158 return gsl_ran_pareto_pdf (x, b, a);
2162 matrix_eval_RV_PARETO (double a, double b)
2164 return gsl_ran_pareto (get_rng (), b, a);
2168 matrix_eval_CDF_RAYLEIGH (double x, double sigma)
2170 return gsl_cdf_rayleigh_P (x, sigma);
2174 matrix_eval_IDF_RAYLEIGH (double P, double sigma)
2176 return gsl_cdf_rayleigh_Pinv (P, sigma);
2180 matrix_eval_PDF_RAYLEIGH (double x, double sigma)
2182 return gsl_ran_rayleigh_pdf (x, sigma);
2186 matrix_eval_RV_RAYLEIGH (double sigma)
2188 return gsl_ran_rayleigh (get_rng (), sigma);
2192 matrix_eval_PDF_RTAIL (double x, double a, double sigma)
2194 return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
2198 matrix_eval_RV_RTAIL (double a, double sigma)
2200 return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
2204 matrix_eval_CDF_T (double x, double df)
2206 return gsl_cdf_tdist_P (x, df);
2210 matrix_eval_TCDF (double x, double df)
2212 return matrix_eval_CDF_T (x, df);
2216 matrix_eval_IDF_T (double P, double df)
2218 return gsl_cdf_tdist_Pinv (P, df);
2222 matrix_eval_PDF_T (double x, double df)
2224 return gsl_ran_tdist_pdf (x, df);
2228 matrix_eval_RV_T (double df)
2230 return gsl_ran_tdist (get_rng (), df);
2234 matrix_eval_CDF_T1G (double x, double a, double b)
2236 return gsl_cdf_gumbel1_P (x, a, b);
2240 matrix_eval_IDF_T1G (double P, double a, double b)
2242 return gsl_cdf_gumbel1_Pinv (P, a, b);
2246 matrix_eval_PDF_T1G (double x, double a, double b)
2248 return gsl_ran_gumbel1_pdf (x, a, b);
2252 matrix_eval_RV_T1G (double a, double b)
2254 return gsl_ran_gumbel1 (get_rng (), a, b);
2258 matrix_eval_CDF_T2G (double x, double a, double b)
2260 return gsl_cdf_gumbel1_P (x, a, b);
2264 matrix_eval_IDF_T2G (double P, double a, double b)
2266 return gsl_cdf_gumbel1_Pinv (P, a, b);
2270 matrix_eval_PDF_T2G (double x, double a, double b)
2272 return gsl_ran_gumbel1_pdf (x, a, b);
2276 matrix_eval_RV_T2G (double a, double b)
2278 return gsl_ran_gumbel1 (get_rng (), a, b);
2282 matrix_eval_CDF_UNIFORM (double x, double a, double b)
2284 return gsl_cdf_flat_P (x, a, b);
2288 matrix_eval_IDF_UNIFORM (double P, double a, double b)
2290 return gsl_cdf_flat_Pinv (P, a, b);
2294 matrix_eval_PDF_UNIFORM (double x, double a, double b)
2296 return gsl_ran_flat_pdf (x, a, b);
2300 matrix_eval_RV_UNIFORM (double a, double b)
2302 return gsl_ran_flat (get_rng (), a, b);
2306 matrix_eval_CDF_WEIBULL (double x, double a, double b)
2308 return gsl_cdf_weibull_P (x, a, b);
2312 matrix_eval_IDF_WEIBULL (double P, double a, double b)
2314 return gsl_cdf_weibull_Pinv (P, a, b);
2318 matrix_eval_PDF_WEIBULL (double x, double a, double b)
2320 return gsl_ran_weibull_pdf (x, a, b);
2324 matrix_eval_RV_WEIBULL (double a, double b)
2326 return gsl_ran_weibull (get_rng (), a, b);
2330 matrix_eval_CDF_BERNOULLI (double k, double p)
2332 return k ? 1 : 1 - p;
2336 matrix_eval_PDF_BERNOULLI (double k, double p)
2338 return gsl_ran_bernoulli_pdf (k, p);
2342 matrix_eval_RV_BERNOULLI (double p)
2344 return gsl_ran_bernoulli (get_rng (), p);
2348 matrix_eval_CDF_BINOM (double k, double n, double p)
2350 return gsl_cdf_binomial_P (k, p, n);
2354 matrix_eval_PDF_BINOM (double k, double n, double p)
2356 return gsl_ran_binomial_pdf (k, p, n);
2360 matrix_eval_RV_BINOM (double n, double p)
2362 return gsl_ran_binomial (get_rng (), p, n);
2366 matrix_eval_CDF_GEOM (double k, double p)
2368 return gsl_cdf_geometric_P (k, p);
2372 matrix_eval_PDF_GEOM (double k, double p)
2374 return gsl_ran_geometric_pdf (k, p);
2378 matrix_eval_RV_GEOM (double p)
2380 return gsl_ran_geometric (get_rng (), p);
2384 matrix_eval_CDF_HYPER (double k, double a, double b, double c)
2386 return gsl_cdf_hypergeometric_P (k, c, a - c, b);
2390 matrix_eval_PDF_HYPER (double k, double a, double b, double c)
2392 return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
2396 matrix_eval_RV_HYPER (double a, double b, double c)
2398 return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
2402 matrix_eval_PDF_LOG (double k, double p)
2404 return gsl_ran_logarithmic_pdf (k, p);
2408 matrix_eval_RV_LOG (double p)
2410 return gsl_ran_logarithmic (get_rng (), p);
2414 matrix_eval_CDF_NEGBIN (double k, double n, double p)
2416 return gsl_cdf_negative_binomial_P (k, p, n);
2420 matrix_eval_PDF_NEGBIN (double k, double n, double p)
2422 return gsl_ran_negative_binomial_pdf (k, p, n);
2426 matrix_eval_RV_NEGBIN (double n, double p)
2428 return gsl_ran_negative_binomial (get_rng (), p, n);
2432 matrix_eval_CDF_POISSON (double k, double mu)
2434 return gsl_cdf_poisson_P (k, mu);
2438 matrix_eval_PDF_POISSON (double k, double mu)
2440 return gsl_ran_poisson_pdf (k, mu);
2444 matrix_eval_RV_POISSON (double mu)
2446 return gsl_ran_poisson (get_rng (), mu);
2449 struct matrix_function
2453 size_t min_args, max_args;
2456 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
2459 word_matches (const char **test, const char **name)
2461 size_t test_len = strcspn (*test, ".");
2462 size_t name_len = strcspn (*name, ".");
2463 if (test_len == name_len)
2465 if (buf_compare_case (*test, *name, test_len))
2468 else if (test_len < 3 || test_len > name_len)
2472 if (buf_compare_case (*test, *name, test_len))
2478 if (**test != **name)
2489 /* Returns 0 if TOKEN and FUNC do not match,
2490 1 if TOKEN is an acceptable abbreviation for FUNC,
2491 2 if TOKEN equals FUNC. */
2493 compare_function_names (const char *token_, const char *func_)
2495 const char *token = token_;
2496 const char *func = func_;
2497 while (*token || *func)
2498 if (!word_matches (&token, &func))
2500 return !c_strcasecmp (token_, func_) ? 2 : 1;
2503 static const struct matrix_function *
2504 matrix_parse_function_name (const char *token)
2506 static const struct matrix_function functions[] =
2508 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
2509 { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
2513 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
2515 for (size_t i = 0; i < N_FUNCTIONS; i++)
2517 if (compare_function_names (token, functions[i].name) > 0)
2518 return &functions[i];
2523 static struct read_file *
2524 read_file_create (struct matrix_state *s, struct file_handle *fh)
2526 for (size_t i = 0; i < s->n_read_files; i++)
2528 struct read_file *rf = s->read_files[i];
2536 struct read_file *rf = xmalloc (sizeof *rf);
2537 *rf = (struct read_file) { .file = fh };
2539 s->read_files = xrealloc (s->read_files,
2540 (s->n_read_files + 1) * sizeof *s->read_files);
2541 s->read_files[s->n_read_files++] = rf;
2545 static struct dfm_reader *
2546 read_file_open (struct read_file *rf)
2549 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
2554 read_file_destroy (struct read_file *rf)
2558 fh_unref (rf->file);
2559 dfm_close_reader (rf->reader);
2560 free (rf->encoding);
2565 static struct write_file *
2566 write_file_create (struct matrix_state *s, struct file_handle *fh)
2568 for (size_t i = 0; i < s->n_write_files; i++)
2570 struct write_file *wf = s->write_files[i];
2578 struct write_file *wf = xmalloc (sizeof *wf);
2579 *wf = (struct write_file) { .file = fh };
2581 s->write_files = xrealloc (s->write_files,
2582 (s->n_write_files + 1) * sizeof *s->write_files);
2583 s->write_files[s->n_write_files++] = wf;
2587 static struct dfm_writer *
2588 write_file_open (struct write_file *wf)
2591 wf->writer = dfm_open_writer (wf->file, wf->encoding);
2596 write_file_destroy (struct write_file *wf)
2602 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
2603 wf->held->s.ss.length);
2604 u8_line_destroy (wf->held);
2608 fh_unref (wf->file);
2609 dfm_close_writer (wf->writer);
2610 free (wf->encoding);
2616 matrix_parse_function (struct matrix_state *s, const char *token,
2617 struct matrix_expr **exprp)
2620 if (lex_next_token (s->lexer, 1) != T_LPAREN)
2623 if (lex_match_id (s->lexer, "EOF"))
2626 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
2630 if (!lex_force_match (s->lexer, T_RPAREN))
2636 struct read_file *rf = read_file_create (s, fh);
2638 struct matrix_expr *e = xmalloc (sizeof *e);
2639 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
2644 const struct matrix_function *f = matrix_parse_function_name (token);
2648 lex_get_n (s->lexer, 2);
2650 struct matrix_expr *e = xmalloc (sizeof *e);
2651 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
2653 if (lex_token (s->lexer) != T_RPAREN)
2655 size_t allocated_subs = 0;
2658 struct matrix_expr *sub = matrix_parse_expr (s);
2662 if (e->n_subs >= allocated_subs)
2663 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
2664 e->subs[e->n_subs++] = sub;
2666 while (lex_match (s->lexer, T_COMMA));
2668 if (!lex_force_match (s->lexer, T_RPAREN))
2671 if (e->n_subs < f->min_args || e->n_subs > f->max_args)
2673 if (f->min_args == f->max_args)
2674 msg (SE, ngettext ("Matrix function %s requires %zu argument.",
2675 "Matrix function %s requires %zu arguments.",
2677 f->name, f->min_args);
2678 else if (f->min_args == 1 && f->max_args == 2)
2679 msg (SE, ngettext ("Matrix function %s requires 1 or 2 arguments, "
2680 "but %zu was provided.",
2681 "Matrix function %s requires 1 or 2 arguments, "
2682 "but %zu were provided.",
2684 f->name, e->n_subs);
2685 else if (f->min_args == 1 && f->max_args == INT_MAX)
2686 msg (SE, _("Matrix function %s requires at least one argument."),
2698 matrix_expr_destroy (e);
2702 static struct matrix_expr *
2703 matrix_parse_primary (struct matrix_state *s)
2705 if (lex_is_number (s->lexer))
2707 double number = lex_number (s->lexer);
2709 return matrix_expr_create_number (number);
2711 else if (lex_is_string (s->lexer))
2713 char string[sizeof (double)];
2714 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
2718 memcpy (&number, string, sizeof number);
2719 return matrix_expr_create_number (number);
2721 else if (lex_match (s->lexer, T_LPAREN))
2723 struct matrix_expr *e = matrix_parse_or_xor (s);
2724 if (!e || !lex_force_match (s->lexer, T_RPAREN))
2726 matrix_expr_destroy (e);
2731 else if (lex_match (s->lexer, T_LCURLY))
2733 struct matrix_expr *e = matrix_parse_curly_semi (s);
2734 if (!e || !lex_force_match (s->lexer, T_RCURLY))
2736 matrix_expr_destroy (e);
2741 else if (lex_token (s->lexer) == T_ID)
2743 struct matrix_expr *retval;
2744 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
2747 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
2750 lex_error (s->lexer, _("Unknown variable %s."),
2751 lex_tokcstr (s->lexer));
2756 struct matrix_expr *e = xmalloc (sizeof *e);
2757 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
2760 else if (lex_token (s->lexer) == T_ALL)
2762 struct matrix_expr *retval;
2763 if (matrix_parse_function (s, "ALL", &retval))
2767 lex_error (s->lexer, NULL);
2771 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
2774 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
2776 if (lex_match (s->lexer, T_COLON))
2783 *indexp = matrix_parse_expr (s);
2784 return *indexp != NULL;
2788 static struct matrix_expr *
2789 matrix_parse_postfix (struct matrix_state *s)
2791 struct matrix_expr *lhs = matrix_parse_primary (s);
2792 if (!lhs || !lex_match (s->lexer, T_LPAREN))
2795 struct matrix_expr *i0;
2796 if (!matrix_parse_index_expr (s, &i0))
2798 matrix_expr_destroy (lhs);
2801 if (lex_match (s->lexer, T_RPAREN))
2803 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
2804 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
2805 else if (lex_match (s->lexer, T_COMMA))
2807 struct matrix_expr *i1;
2808 if (!matrix_parse_index_expr (s, &i1)
2809 || !lex_force_match (s->lexer, T_RPAREN))
2811 matrix_expr_destroy (lhs);
2812 matrix_expr_destroy (i0);
2813 matrix_expr_destroy (i1);
2816 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
2817 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
2818 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
2823 lex_error_expecting (s->lexer, "`)'", "`,'");
2828 static struct matrix_expr *
2829 matrix_parse_unary (struct matrix_state *s)
2831 if (lex_match (s->lexer, T_DASH))
2833 struct matrix_expr *lhs = matrix_parse_unary (s);
2834 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
2836 else if (lex_match (s->lexer, T_PLUS))
2837 return matrix_parse_unary (s);
2839 return matrix_parse_postfix (s);
2842 static struct matrix_expr *
2843 matrix_parse_seq (struct matrix_state *s)
2845 struct matrix_expr *start = matrix_parse_unary (s);
2846 if (!start || !lex_match (s->lexer, T_COLON))
2849 struct matrix_expr *end = matrix_parse_unary (s);
2852 matrix_expr_destroy (start);
2856 if (lex_match (s->lexer, T_COLON))
2858 struct matrix_expr *increment = matrix_parse_unary (s);
2861 matrix_expr_destroy (start);
2862 matrix_expr_destroy (end);
2865 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2868 return matrix_expr_create_2 (MOP_SEQ, start, end);
2871 static struct matrix_expr *
2872 matrix_parse_exp (struct matrix_state *s)
2874 struct matrix_expr *lhs = matrix_parse_seq (s);
2881 if (lex_match (s->lexer, T_EXP))
2883 else if (lex_match_phrase (s->lexer, "&**"))
2888 struct matrix_expr *rhs = matrix_parse_seq (s);
2891 matrix_expr_destroy (lhs);
2894 lhs = matrix_expr_create_2 (op, lhs, rhs);
2898 static struct matrix_expr *
2899 matrix_parse_mul_div (struct matrix_state *s)
2901 struct matrix_expr *lhs = matrix_parse_exp (s);
2908 if (lex_match (s->lexer, T_ASTERISK))
2910 else if (lex_match (s->lexer, T_SLASH))
2912 else if (lex_match_phrase (s->lexer, "&*"))
2914 else if (lex_match_phrase (s->lexer, "&/"))
2919 struct matrix_expr *rhs = matrix_parse_exp (s);
2922 matrix_expr_destroy (lhs);
2925 lhs = matrix_expr_create_2 (op, lhs, rhs);
2929 static struct matrix_expr *
2930 matrix_parse_add_sub (struct matrix_state *s)
2932 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2939 if (lex_match (s->lexer, T_PLUS))
2941 else if (lex_match (s->lexer, T_DASH))
2943 else if (lex_token (s->lexer) == T_NEG_NUM)
2948 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2951 matrix_expr_destroy (lhs);
2954 lhs = matrix_expr_create_2 (op, lhs, rhs);
2958 static struct matrix_expr *
2959 matrix_parse_relational (struct matrix_state *s)
2961 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2968 if (lex_match (s->lexer, T_GT))
2970 else if (lex_match (s->lexer, T_GE))
2972 else if (lex_match (s->lexer, T_LT))
2974 else if (lex_match (s->lexer, T_LE))
2976 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2978 else if (lex_match (s->lexer, T_NE))
2983 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2986 matrix_expr_destroy (lhs);
2989 lhs = matrix_expr_create_2 (op, lhs, rhs);
2993 static struct matrix_expr *
2994 matrix_parse_not (struct matrix_state *s)
2996 if (lex_match (s->lexer, T_NOT))
2998 struct matrix_expr *lhs = matrix_parse_not (s);
2999 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
3002 return matrix_parse_relational (s);
3005 static struct matrix_expr *
3006 matrix_parse_and (struct matrix_state *s)
3008 struct matrix_expr *lhs = matrix_parse_not (s);
3012 while (lex_match (s->lexer, T_AND))
3014 struct matrix_expr *rhs = matrix_parse_not (s);
3017 matrix_expr_destroy (lhs);
3020 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
3025 static struct matrix_expr *
3026 matrix_parse_or_xor (struct matrix_state *s)
3028 struct matrix_expr *lhs = matrix_parse_and (s);
3035 if (lex_match (s->lexer, T_OR))
3037 else if (lex_match_id (s->lexer, "XOR"))
3042 struct matrix_expr *rhs = matrix_parse_and (s);
3045 matrix_expr_destroy (lhs);
3048 lhs = matrix_expr_create_2 (op, lhs, rhs);
3052 static struct matrix_expr *
3053 matrix_parse_expr (struct matrix_state *s)
3055 return matrix_parse_or_xor (s);
3058 /* Expression evaluation. */
3061 matrix_op_eval (enum matrix_op op, double a, double b)
3065 case MOP_ADD_ELEMS: return a + b;
3066 case MOP_SUB_ELEMS: return a - b;
3067 case MOP_MUL_ELEMS: return a * b;
3068 case MOP_DIV_ELEMS: return a / b;
3069 case MOP_EXP_ELEMS: return pow (a, b);
3070 case MOP_GT: return a > b;
3071 case MOP_GE: return a >= b;
3072 case MOP_LT: return a < b;
3073 case MOP_LE: return a <= b;
3074 case MOP_EQ: return a == b;
3075 case MOP_NE: return a != b;
3076 case MOP_AND: return (a > 0) && (b > 0);
3077 case MOP_OR: return (a > 0) || (b > 0);
3078 case MOP_XOR: return (a > 0) != (b > 0);
3080 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3089 case MOP_PASTE_HORZ:
3090 case MOP_PASTE_VERT:
3106 matrix_op_name (enum matrix_op op)
3110 case MOP_ADD_ELEMS: return "+";
3111 case MOP_SUB_ELEMS: return "-";
3112 case MOP_MUL_ELEMS: return "&*";
3113 case MOP_DIV_ELEMS: return "&/";
3114 case MOP_EXP_ELEMS: return "&**";
3115 case MOP_GT: return ">";
3116 case MOP_GE: return ">=";
3117 case MOP_LT: return "<";
3118 case MOP_LE: return "<=";
3119 case MOP_EQ: return "=";
3120 case MOP_NE: return "<>";
3121 case MOP_AND: return "AND";
3122 case MOP_OR: return "OR";
3123 case MOP_XOR: return "XOR";
3125 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3134 case MOP_PASTE_HORZ:
3135 case MOP_PASTE_VERT:
3151 is_scalar (const gsl_matrix *m)
3153 return m->size1 == 1 && m->size2 == 1;
3157 to_scalar (const gsl_matrix *m)
3159 assert (is_scalar (m));
3160 return gsl_matrix_get (m, 0, 0);
3164 matrix_expr_evaluate_elementwise (enum matrix_op op,
3165 gsl_matrix *a, gsl_matrix *b)
3169 double be = to_scalar (b);
3170 for (size_t r = 0; r < a->size1; r++)
3171 for (size_t c = 0; c < a->size2; c++)
3173 double *ae = gsl_matrix_ptr (a, r, c);
3174 *ae = matrix_op_eval (op, *ae, be);
3178 else if (is_scalar (a))
3180 double ae = to_scalar (a);
3181 for (size_t r = 0; r < b->size1; r++)
3182 for (size_t c = 0; c < b->size2; c++)
3184 double *be = gsl_matrix_ptr (b, r, c);
3185 *be = matrix_op_eval (op, ae, *be);
3189 else if (a->size1 == b->size1 && a->size2 == b->size2)
3191 for (size_t r = 0; r < a->size1; r++)
3192 for (size_t c = 0; c < a->size2; c++)
3194 double *ae = gsl_matrix_ptr (a, r, c);
3195 double be = gsl_matrix_get (b, r, c);
3196 *ae = matrix_op_eval (op, *ae, be);
3202 msg (SE, _("Operands to %s must have the same dimensions or one "
3203 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
3204 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
3210 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
3212 if (is_scalar (a) || is_scalar (b))
3213 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
3215 if (a->size2 != b->size1)
3217 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
3218 "not conformable for multiplication."),
3219 a->size1, a->size2, b->size1, b->size2);
3223 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3224 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3229 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3231 gsl_matrix *tmp = *a;
3237 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3240 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3241 swap_matrix (z, tmp);
3245 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3247 mul_matrix (x, *x, *x, tmp);
3251 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
3254 if (x->size1 != x->size2)
3256 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
3257 "the left-hand size, not one with dimensions %zu×%zu."),
3258 x->size1, x->size2);
3263 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
3264 "right-hand side, not a matrix with dimensions %zu×%zu."),
3265 b->size1, b->size2);
3268 double bf = to_scalar (b);
3269 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3271 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
3272 "or outside the valid range."), bf);
3277 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3279 gsl_matrix_set_identity (y);
3283 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3285 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3288 mul_matrix (&y, x, y, &t);
3289 square_matrix (&x, &t);
3292 square_matrix (&x, &t);
3294 mul_matrix (&y, x, y, &t);
3298 /* Garbage collection.
3300 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3301 a permutation of them. We are returning one of them; that one must not be
3302 destroyed. We must not destroy 'x_' because the caller owns it. */
3304 gsl_matrix_free (y_);
3306 gsl_matrix_free (t_);
3312 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
3315 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3317 msg (SE, _("All operands of : operator must be scalars."));
3321 long int start = to_scalar (start_);
3322 long int end = to_scalar (end_);
3323 long int by = by_ ? to_scalar (by_) : 1;
3327 msg (SE, _("The increment operand to : must be nonzero."));
3331 long int n = (end >= start && by > 0 ? (end - start + by) / by
3332 : end <= start && by < 0 ? (start - end - by) / -by
3334 gsl_matrix *m = gsl_matrix_alloc (1, n);
3335 for (long int i = 0; i < n; i++)
3336 gsl_matrix_set (m, 0, i, start + i * by);
3341 matrix_expr_evaluate_not (gsl_matrix *a)
3343 for (size_t r = 0; r < a->size1; r++)
3344 for (size_t c = 0; c < a->size2; c++)
3346 double *ae = gsl_matrix_ptr (a, r, c);
3353 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
3355 if (a->size1 != b->size1)
3357 if (!a->size1 || !a->size2)
3359 else if (!b->size1 || !b->size2)
3362 msg (SE, _("All columns in a matrix must have the same number of rows, "
3363 "but this tries to paste matrices with %zu and %zu rows."),
3364 a->size1, b->size1);
3368 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3369 for (size_t y = 0; y < a->size1; y++)
3371 for (size_t x = 0; x < a->size2; x++)
3372 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3373 for (size_t x = 0; x < b->size2; x++)
3374 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3380 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
3382 if (a->size2 != b->size2)
3384 if (!a->size1 || !a->size2)
3386 else if (!b->size1 || !b->size2)
3389 msg (SE, _("All rows in a matrix must have the same number of columns, "
3390 "but this tries to stack matrices with %zu and %zu columns."),
3391 a->size2, b->size2);
3395 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3396 for (size_t x = 0; x < a->size2; x++)
3398 for (size_t y = 0; y < a->size1; y++)
3399 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3400 for (size_t y = 0; y < b->size1; y++)
3401 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3407 matrix_to_vector (gsl_matrix *m)
3410 gsl_vector v = to_vector (m);
3411 assert (v.block == m->block || !v.block);
3415 gsl_matrix_free (m);
3416 return xmemdup (&v, sizeof v);
3430 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3433 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
3434 enum index_type index_type, size_t other_size,
3435 struct index_vector *iv)
3444 msg (SE, _("Vector index must be scalar or vector, not a "
3446 m->size1, m->size2);
3450 msg (SE, _("Matrix row index must be scalar or vector, not a "
3452 m->size1, m->size2);
3456 msg (SE, _("Matrix column index must be scalar or vector, not a "
3458 m->size1, m->size2);
3464 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3465 *iv = (struct index_vector) {
3466 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3469 for (size_t i = 0; i < v.size; i++)
3471 double index = gsl_vector_get (&v, i);
3472 if (index < 1 || index >= size + 1)
3477 msg (SE, _("Index %g is out of range for vector "
3478 "with %zu elements."), index, size);
3482 msg (SE, _("%g is not a valid row index for "
3483 "a %zu×%zu matrix."),
3484 index, size, other_size);
3488 msg (SE, _("%g is not a valid column index for "
3489 "a %zu×%zu matrix."),
3490 index, other_size, size);
3497 iv->indexes[i] = index - 1;
3503 *iv = (struct index_vector) {
3504 .indexes = xnmalloc (size, sizeof *iv->indexes),
3507 for (size_t i = 0; i < size; i++)
3514 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
3516 if (!is_vector (sm))
3518 msg (SE, _("Vector index operator may not be applied to "
3519 "a %zu×%zu matrix."),
3520 sm->size1, sm->size2);
3528 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
3530 if (!matrix_expr_evaluate_vec_all (sm))
3533 gsl_vector sv = to_vector (sm);
3534 struct index_vector iv;
3535 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
3538 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3539 sm->size1 == 1 ? iv.n : 1);
3540 gsl_vector dv = to_vector (dm);
3541 for (size_t dx = 0; dx < iv.n; dx++)
3543 size_t sx = iv.indexes[dx];
3544 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3552 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
3555 struct index_vector iv0;
3556 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
3559 struct index_vector iv1;
3560 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
3567 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3568 for (size_t dy = 0; dy < iv0.n; dy++)
3570 size_t sy = iv0.indexes[dy];
3572 for (size_t dx = 0; dx < iv1.n; dx++)
3574 size_t sx = iv1.indexes[dx];
3575 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3583 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3584 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3585 const struct matrix_function_properties *, gsl_matrix *[], size_t, \
3586 matrix_proto_##PROTO *);
3591 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
3593 if (!is_scalar (subs[index]))
3595 msg (SE, _("Function %s argument %zu must be a scalar, "
3596 "not a %zu×%zu matrix."),
3597 name, index + 1, subs[index]->size1, subs[index]->size2);
3604 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
3606 if (!is_vector (subs[index]))
3608 msg (SE, _("Function %s argument %zu must be a vector, "
3609 "not a %zu×%zu matrix."),
3610 name, index + 1, subs[index]->size1, subs[index]->size2);
3617 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
3619 for (size_t i = 0; i < n_subs; i++)
3621 if (!check_scalar_arg (name, subs, i))
3623 d[i] = to_scalar (subs[i]);
3629 parse_constraint_value (const char **constraintsp)
3632 long retval = strtol (*constraintsp, &tail, 10);
3633 assert (tail > *constraintsp);
3634 *constraintsp = tail;
3639 argument_invalid (const struct matrix_function_properties *props,
3640 size_t arg_index, gsl_matrix *a, size_t y, size_t x,
3644 ds_put_format (out, _("Argument %zu to matrix function %s "
3645 "has invalid value %g."),
3646 arg_index, props->name, gsl_matrix_get (a, y, x));
3648 ds_put_format (out, _("Row %zu, column %zu of argument %zu "
3649 "to matrix function %s has "
3650 "invalid value %g."),
3651 y + 1, x + 1, arg_index, props->name,
3652 gsl_matrix_get (a, y, x));
3656 argument_inequality_error (const struct matrix_function_properties *props,
3657 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3658 size_t b_index, double b,
3659 bool greater, bool equal)
3661 struct string s = DS_EMPTY_INITIALIZER;
3663 argument_invalid (props, a_index, a, y, x, &s);
3664 ds_put_cstr (&s, " ");
3665 if (greater && equal)
3666 ds_put_format (&s, _("This argument must be greater than or "
3667 "equal to argument %zu, but that has value %g."),
3669 else if (greater && !equal)
3670 ds_put_format (&s, _("This argument must be greater than argument %zu, "
3671 "but that has value %g."),
3674 ds_put_format (&s, _("This argument must be less than or "
3675 "equal to argument %zu, but that has value %g."),
3678 ds_put_format (&s, _("This argument must be less than argument %zu, "
3679 "but that has value %g."),
3681 msg (ME, ds_cstr (&s));
3686 argument_value_error (const struct matrix_function_properties *props,
3687 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3689 bool greater, bool equal)
3691 struct string s = DS_EMPTY_INITIALIZER;
3692 argument_invalid (props, a_index, a, y, x, &s);
3693 ds_put_cstr (&s, " ");
3694 if (greater && equal)
3695 ds_put_format (&s, _("This argument must be greater than or equal to %g."),
3697 else if (greater && !equal)
3698 ds_put_format (&s, _("This argument must be greater than %g."), b);
3700 ds_put_format (&s, _("This argument must be less than or equal to %g."), b);
3702 ds_put_format (&s, _("This argument must be less than %g."), b);
3703 msg (ME, ds_cstr (&s));
3708 check_constraints (const struct matrix_function_properties *props,
3709 gsl_matrix *args[], size_t n_args)
3711 const char *constraints = props->constraints;
3715 size_t arg_index = SIZE_MAX;
3716 while (*constraints)
3718 if (*constraints >= 'a' && *constraints <= 'd')
3720 arg_index = *constraints++ - 'a';
3721 assert (arg_index < n_args);
3723 else if (*constraints == '[' || *constraints == '(')
3725 assert (arg_index < n_args);
3726 bool open_lower = *constraints++ == '(';
3727 int minimum = parse_constraint_value (&constraints);
3728 assert (*constraints == ',');
3730 int maximum = parse_constraint_value (&constraints);
3731 assert (*constraints == ']' || *constraints == ')');
3732 bool open_upper = *constraints++ == ')';
3734 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3735 if ((open_lower ? *d <= minimum : *d < minimum)
3736 || (open_upper ? *d >= maximum : *d > maximum))
3738 if (!is_scalar (args[arg_index]))
3739 msg (ME, _("Row %zu, column %zu of argument %zu to matrix "
3740 "function %s has value %g, which is outside "
3741 "the valid range %c%d,%d%c."),
3742 y + 1, x + 1, arg_index + 1, props->name, *d,
3743 open_lower ? '(' : '[',
3745 open_upper ? ')' : ']');
3747 msg (ME, _("Argument %zu to matrix function %s has value %g, "
3748 "which is outside the valid range %c%d,%d%c."),
3749 arg_index + 1, props->name, *d,
3750 open_lower ? '(' : '[',
3752 open_upper ? ')' : ']');
3756 else if (*constraints == '>' || *constraints == '<')
3758 bool greater = *constraints++ == '>';
3760 if (*constraints == '=')
3768 if (*constraints >= 'a' && *constraints <= 'd')
3770 size_t a_index = arg_index;
3771 size_t b_index = *constraints - 'a';
3772 assert (a_index < n_args);
3773 assert (b_index < n_args);
3775 /* We only support one of the two arguments being non-scalar.
3776 It's easier to support only the first one being non-scalar, so
3777 flip things around if it's the other way. */
3778 if (!is_scalar (args[b_index]))
3780 assert (is_scalar (args[a_index]));
3781 size_t tmp_index = a_index;
3783 b_index = tmp_index;
3788 double b = to_scalar (args[b_index]);
3789 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
3791 ? (equal ? !(*a >= b) : !(*a > b))
3792 : (equal ? !(*a <= b) : !(*a < b)))
3794 argument_inequality_error (props,
3795 a_index + 1, args[a_index], y, x,
3803 int comparand = parse_constraint_value (&constraints);
3805 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3807 ? (equal ? !(*d >= comparand) : !(*d > comparand))
3808 : (equal ? !(*d <= comparand) : !(*d < comparand)))
3810 argument_value_error (props,
3811 arg_index + 1, args[arg_index], y, x,
3820 assert (*constraints == ' ');
3822 arg_index = SIZE_MAX;
3829 matrix_expr_evaluate_d_none (
3830 const struct matrix_function_properties *props UNUSED,
3831 gsl_matrix *subs[] UNUSED, size_t n_subs,
3832 matrix_proto_d_none *f)
3834 assert (n_subs == 0);
3836 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3837 gsl_matrix_set (m, 0, 0, f ());
3842 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
3843 gsl_matrix *subs[], size_t n_subs,
3844 matrix_proto_d_d *f)
3846 assert (n_subs == 1);
3849 if (!to_scalar_args (props->name, subs, n_subs, &d))
3852 if (!check_constraints (props, subs, n_subs))
3855 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3856 gsl_matrix_set (m, 0, 0, f (d));
3861 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
3862 gsl_matrix *subs[], size_t n_subs,
3863 matrix_proto_d_dd *f)
3865 assert (n_subs == 2);
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[0], d[1]));
3880 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
3881 gsl_matrix *subs[], size_t n_subs,
3882 matrix_proto_d_ddd *f)
3884 assert (n_subs == 3);
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], d[2]));
3899 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
3900 gsl_matrix *subs[], size_t n_subs,
3901 matrix_proto_m_d *f)
3903 assert (n_subs == 1);
3906 return to_scalar_args (props->name, subs, n_subs, &d) ? f (d) : NULL;
3910 matrix_expr_evaluate_m_dd (const struct matrix_function_properties *props,
3911 gsl_matrix *subs[], size_t n_subs,
3912 matrix_proto_m_dd *f)
3914 assert (n_subs == 2);
3917 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
3921 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
3922 gsl_matrix *subs[], size_t n_subs,
3923 matrix_proto_m_ddd *f)
3925 assert (n_subs == 3);
3928 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
3932 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED,
3933 gsl_matrix *subs[], size_t n_subs,
3934 matrix_proto_m_m *f)
3936 assert (n_subs == 1);
3941 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
3942 gsl_matrix *subs[], size_t n_subs,
3943 matrix_proto_m_e *f)
3945 assert (n_subs == 1);
3947 if (!check_constraints (props, subs, n_subs))
3950 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3956 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED,
3957 gsl_matrix *subs[], size_t n_subs,
3958 matrix_proto_m_md *f)
3960 assert (n_subs == 2);
3961 if (!check_scalar_arg (props->name, subs, 1))
3963 return f (subs[0], to_scalar (subs[1]));
3967 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
3968 gsl_matrix *subs[], size_t n_subs,
3969 matrix_proto_m_ed *f)
3971 assert (n_subs == 2);
3972 if (!check_scalar_arg (props->name, subs, 1))
3975 if (!check_constraints (props, subs, n_subs))
3978 double b = to_scalar (subs[1]);
3979 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3985 matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props UNUSED,
3986 gsl_matrix *subs[], size_t n_subs,
3987 matrix_proto_m_mdd *f)
3989 assert (n_subs == 3);
3990 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
3992 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
3996 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
3997 gsl_matrix *subs[], size_t n_subs,
3998 matrix_proto_m_edd *f)
4000 assert (n_subs == 3);
4001 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
4004 if (!check_constraints (props, subs, n_subs))
4007 double b = to_scalar (subs[1]);
4008 double c = to_scalar (subs[2]);
4009 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4015 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4016 gsl_matrix *subs[], size_t n_subs,
4017 matrix_proto_m_eddd *f)
4019 assert (n_subs == 4);
4020 for (size_t i = 1; i < 4; i++)
4021 if (!check_scalar_arg (props->name, subs, i))
4024 if (!check_constraints (props, subs, n_subs))
4027 double b = to_scalar (subs[1]);
4028 double c = to_scalar (subs[2]);
4029 double d = to_scalar (subs[3]);
4030 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4031 *a = f (*a, b, c, d);
4036 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4037 gsl_matrix *subs[], size_t n_subs,
4038 matrix_proto_m_eed *f)
4040 assert (n_subs == 3);
4041 if (!check_scalar_arg (props->name, subs, 2))
4044 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4045 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4047 msg (ME, _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4048 "%zu×%zu, but %s requires these arguments either to have "
4049 "the same dimensions or for one of them to be a scalar."),
4051 subs[0]->size1, subs[0]->size2,
4052 subs[1]->size1, subs[1]->size2,
4057 if (!check_constraints (props, subs, n_subs))
4060 double c = to_scalar (subs[2]);
4062 if (is_scalar (subs[0]))
4064 double a = to_scalar (subs[0]);
4065 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4071 double b = to_scalar (subs[1]);
4072 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4079 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED,
4080 gsl_matrix *subs[], size_t n_subs,
4081 matrix_proto_m_mm *f)
4083 assert (n_subs == 2);
4084 return f (subs[0], subs[1]);
4088 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4089 gsl_matrix *subs[], size_t n_subs,
4090 matrix_proto_m_v *f)
4092 assert (n_subs == 1);
4093 if (!check_vector_arg (props->name, subs, 0))
4095 gsl_vector v = to_vector (subs[0]);
4100 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED,
4101 gsl_matrix *subs[], size_t n_subs,
4102 matrix_proto_d_m *f)
4104 assert (n_subs == 1);
4106 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4107 gsl_matrix_set (m, 0, 0, f (subs[0]));
4112 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED,
4113 gsl_matrix *subs[], size_t n_subs,
4114 matrix_proto_m_any *f)
4116 return f (subs, n_subs);
4120 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props,
4121 gsl_matrix *subs[], size_t n_subs,
4122 matrix_proto_IDENT *f)
4124 assert (n_subs <= 2);
4127 if (!to_scalar_args (props->name, subs, n_subs, d))
4130 return f (d[0], d[n_subs - 1]);
4134 matrix_expr_evaluate (const struct matrix_expr *e)
4136 if (e->op == MOP_NUMBER)
4138 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4139 gsl_matrix_set (m, 0, 0, e->number);
4142 else if (e->op == MOP_VARIABLE)
4144 const gsl_matrix *src = e->variable->value;
4147 msg (SE, _("Uninitialized variable %s used in expression."),
4152 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4153 gsl_matrix_memcpy (dst, src);
4156 else if (e->op == MOP_EOF)
4158 struct dfm_reader *reader = read_file_open (e->eof);
4159 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4160 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4164 enum { N_LOCAL = 3 };
4165 gsl_matrix *local_subs[N_LOCAL];
4166 gsl_matrix **subs = (e->n_subs < N_LOCAL
4168 : xmalloc (e->n_subs * sizeof *subs));
4170 for (size_t i = 0; i < e->n_subs; i++)
4172 subs[i] = matrix_expr_evaluate (e->subs[i]);
4175 for (size_t j = 0; j < i; j++)
4176 gsl_matrix_free (subs[j]);
4177 if (subs != local_subs)
4183 gsl_matrix *result = NULL;
4186 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4187 case MOP_F_##ENUM: \
4189 static const struct matrix_function_properties props = { \
4191 .constraints = CONSTRAINTS, \
4193 result = matrix_expr_evaluate_##PROTO (&props, subs, e->n_subs, \
4194 matrix_eval_##ENUM); \
4201 gsl_matrix_scale (subs[0], -1.0);
4219 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
4223 result = matrix_expr_evaluate_not (subs[0]);
4227 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
4231 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
4235 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
4239 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
4242 case MOP_PASTE_HORZ:
4243 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
4246 case MOP_PASTE_VERT:
4247 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
4251 result = gsl_matrix_alloc (0, 0);
4255 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
4259 result = matrix_expr_evaluate_vec_all (subs[0]);
4263 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
4267 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
4271 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
4280 for (size_t i = 0; i < e->n_subs; i++)
4281 if (subs[i] != result)
4282 gsl_matrix_free (subs[i]);
4283 if (subs != local_subs)
4289 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4292 gsl_matrix *m = matrix_expr_evaluate (e);
4298 msg (SE, _("Expression for %s must evaluate to scalar, "
4299 "not a %zu×%zu matrix."),
4300 context, m->size1, m->size2);
4301 gsl_matrix_free (m);
4306 gsl_matrix_free (m);
4311 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4315 if (!matrix_expr_evaluate_scalar (e, context, &d))
4319 if (d < LONG_MIN || d > LONG_MAX)
4321 msg (SE, _("Expression for %s is outside the integer range."), context);
4328 struct matrix_lvalue
4330 struct matrix_var *var;
4331 struct matrix_expr *indexes[2];
4336 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4340 for (size_t i = 0; i < lvalue->n_indexes; i++)
4341 matrix_expr_destroy (lvalue->indexes[i]);
4346 static struct matrix_lvalue *
4347 matrix_lvalue_parse (struct matrix_state *s)
4349 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4350 if (!lex_force_id (s->lexer))
4352 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4353 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4357 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
4362 lex_get_n (s->lexer, 2);
4364 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
4366 lvalue->n_indexes++;
4368 if (lex_match (s->lexer, T_COMMA))
4370 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
4372 lvalue->n_indexes++;
4374 if (!lex_force_match (s->lexer, T_RPAREN))
4380 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4386 matrix_lvalue_destroy (lvalue);
4391 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4392 enum index_type index_type, size_t other_size,
4393 struct index_vector *iv)
4398 m = matrix_expr_evaluate (e);
4405 bool ok = matrix_normalize_index_vector (m, size, index_type,
4407 gsl_matrix_free (m);
4412 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4413 struct index_vector *iv, gsl_matrix *sm)
4415 gsl_vector dv = to_vector (lvalue->var->value);
4417 /* Convert source matrix 'sm' to source vector 'sv'. */
4418 if (!is_vector (sm))
4420 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
4421 sm->size1, sm->size2);
4424 gsl_vector sv = to_vector (sm);
4425 if (iv->n != sv.size)
4427 msg (SE, _("Can't assign %zu-element vector "
4428 "to %zu-element subvector."), sv.size, iv->n);
4432 /* Assign elements. */
4433 for (size_t x = 0; x < iv->n; x++)
4434 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4439 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4440 struct index_vector *iv0,
4441 struct index_vector *iv1, gsl_matrix *sm)
4443 gsl_matrix *dm = lvalue->var->value;
4445 /* Convert source matrix 'sm' to source vector 'sv'. */
4446 if (iv0->n != sm->size1)
4448 msg (SE, _("Row index vector for assignment to %s has %zu elements "
4449 "but source matrix has %zu rows."),
4450 lvalue->var->name, iv0->n, sm->size1);
4453 if (iv1->n != sm->size2)
4455 msg (SE, _("Column index vector for assignment to %s has %zu elements "
4456 "but source matrix has %zu columns."),
4457 lvalue->var->name, iv1->n, sm->size2);
4461 /* Assign elements. */
4462 for (size_t y = 0; y < iv0->n; y++)
4464 size_t f0 = iv0->indexes[y];
4465 for (size_t x = 0; x < iv1->n; x++)
4467 size_t f1 = iv1->indexes[x];
4468 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4474 /* Takes ownership of M. */
4476 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4477 struct index_vector *iv0, struct index_vector *iv1,
4480 if (!lvalue->n_indexes)
4482 gsl_matrix_free (lvalue->var->value);
4483 lvalue->var->value = sm;
4488 bool ok = (lvalue->n_indexes == 1
4489 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
4490 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
4491 free (iv0->indexes);
4492 free (iv1->indexes);
4493 gsl_matrix_free (sm);
4499 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4500 struct index_vector *iv0,
4501 struct index_vector *iv1)
4503 *iv0 = INDEX_VECTOR_INIT;
4504 *iv1 = INDEX_VECTOR_INIT;
4505 if (!lvalue->n_indexes)
4508 /* Validate destination matrix exists and has the right shape. */
4509 gsl_matrix *dm = lvalue->var->value;
4512 msg (SE, _("Undefined variable %s."), lvalue->var->name);
4515 else if (dm->size1 == 0 || dm->size2 == 0)
4517 msg (SE, _("Cannot index %zu×%zu matrix."),
4518 dm->size1, dm->size2);
4521 else if (lvalue->n_indexes == 1)
4523 if (!is_vector (dm))
4525 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
4526 dm->size1, dm->size2, lvalue->var->name);
4529 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4530 MAX (dm->size1, dm->size2),
4535 assert (lvalue->n_indexes == 2);
4536 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4537 IV_ROW, dm->size2, iv0))
4540 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4541 IV_COLUMN, dm->size1, iv1))
4543 free (iv0->indexes);
4550 /* Takes ownership of M. */
4552 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
4554 struct index_vector iv0, iv1;
4555 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
4557 gsl_matrix_free (sm);
4561 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
4565 /* Matrix command. */
4569 struct matrix_cmd **commands;
4573 static void matrix_cmds_uninit (struct matrix_cmds *);
4577 enum matrix_cmd_type
4600 struct compute_command
4602 struct matrix_lvalue *lvalue;
4603 struct matrix_expr *rvalue;
4607 struct print_command
4609 struct matrix_expr *expression;
4610 bool use_default_format;
4611 struct fmt_spec format;
4613 int space; /* -1 means NEWPAGE. */
4617 struct string_array literals; /* CLABELS/RLABELS. */
4618 struct matrix_expr *expr; /* CNAMES/RNAMES. */
4624 struct do_if_command
4628 struct matrix_expr *condition;
4629 struct matrix_cmds commands;
4639 struct matrix_var *var;
4640 struct matrix_expr *start;
4641 struct matrix_expr *end;
4642 struct matrix_expr *increment;
4644 /* Loop conditions. */
4645 struct matrix_expr *top_condition;
4646 struct matrix_expr *bottom_condition;
4649 struct matrix_cmds commands;
4653 struct display_command
4655 struct matrix_state *state;
4659 struct release_command
4661 struct matrix_var **vars;
4668 struct matrix_expr *expression;
4669 struct save_file *sf;
4675 struct read_file *rf;
4676 struct matrix_lvalue *dst;
4677 struct matrix_expr *size;
4679 enum fmt_type format;
4686 struct write_command
4688 struct write_file *wf;
4689 struct matrix_expr *expression;
4692 /* If this is nonnull, WRITE uses this format.
4694 If this is NULL, WRITE uses free-field format with as many
4695 digits of precision as needed. */
4696 struct fmt_spec *format;
4705 struct matrix_lvalue *dst;
4706 struct dataset *dataset;
4707 struct file_handle *file;
4709 struct var_syntax *vars;
4711 struct matrix_var *names;
4713 /* Treatment of missing values. */
4718 MGET_ERROR, /* Flag error on command. */
4719 MGET_ACCEPT, /* Accept without change, user-missing only. */
4720 MGET_OMIT, /* Drop this case. */
4721 MGET_RECODE /* Recode to 'substitute'. */
4724 double substitute; /* MGET_RECODE only. */
4730 struct msave_command
4732 struct msave_common *common;
4733 struct matrix_expr *expr;
4734 const char *rowtype;
4735 const struct matrix_expr *factors;
4736 const struct matrix_expr *splits;
4742 struct matrix_state *state;
4743 struct file_handle *file;
4745 struct stringi_set rowtypes;
4749 struct eigen_command
4751 struct matrix_expr *expr;
4752 struct matrix_var *evec;
4753 struct matrix_var *eval;
4757 struct setdiag_command
4759 struct matrix_var *dst;
4760 struct matrix_expr *expr;
4766 struct matrix_expr *expr;
4767 struct matrix_var *u;
4768 struct matrix_var *s;
4769 struct matrix_var *v;
4775 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
4776 static bool matrix_cmd_execute (struct matrix_cmd *);
4777 static void matrix_cmd_destroy (struct matrix_cmd *);
4780 static struct matrix_cmd *
4781 matrix_parse_compute (struct matrix_state *s)
4783 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4784 *cmd = (struct matrix_cmd) {
4785 .type = MCMD_COMPUTE,
4786 .compute = { .lvalue = NULL },
4789 struct compute_command *compute = &cmd->compute;
4790 compute->lvalue = matrix_lvalue_parse (s);
4791 if (!compute->lvalue)
4794 if (!lex_force_match (s->lexer, T_EQUALS))
4797 compute->rvalue = matrix_parse_expr (s);
4798 if (!compute->rvalue)
4804 matrix_cmd_destroy (cmd);
4809 matrix_cmd_execute_compute (struct compute_command *compute)
4811 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
4815 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
4819 print_labels_destroy (struct print_labels *labels)
4823 string_array_destroy (&labels->literals);
4824 matrix_expr_destroy (labels->expr);
4830 parse_literal_print_labels (struct matrix_state *s,
4831 struct print_labels **labelsp)
4833 lex_match (s->lexer, T_EQUALS);
4834 print_labels_destroy (*labelsp);
4835 *labelsp = xzalloc (sizeof **labelsp);
4836 while (lex_token (s->lexer) != T_SLASH
4837 && lex_token (s->lexer) != T_ENDCMD
4838 && lex_token (s->lexer) != T_STOP)
4840 struct string label = DS_EMPTY_INITIALIZER;
4841 while (lex_token (s->lexer) != T_COMMA
4842 && lex_token (s->lexer) != T_SLASH
4843 && lex_token (s->lexer) != T_ENDCMD
4844 && lex_token (s->lexer) != T_STOP)
4846 if (!ds_is_empty (&label))
4847 ds_put_byte (&label, ' ');
4849 if (lex_token (s->lexer) == T_STRING)
4850 ds_put_cstr (&label, lex_tokcstr (s->lexer));
4853 char *rep = lex_next_representation (s->lexer, 0, 0);
4854 ds_put_cstr (&label, rep);
4859 string_array_append_nocopy (&(*labelsp)->literals,
4860 ds_steal_cstr (&label));
4862 if (!lex_match (s->lexer, T_COMMA))
4868 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
4870 lex_match (s->lexer, T_EQUALS);
4871 struct matrix_expr *e = matrix_parse_exp (s);
4875 print_labels_destroy (*labelsp);
4876 *labelsp = xzalloc (sizeof **labelsp);
4877 (*labelsp)->expr = e;
4881 static struct matrix_cmd *
4882 matrix_parse_print (struct matrix_state *s)
4884 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4885 *cmd = (struct matrix_cmd) {
4888 .use_default_format = true,
4892 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
4895 for (size_t i = 0; ; i++)
4897 enum token_type t = lex_next_token (s->lexer, i);
4898 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
4900 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
4902 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
4905 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
4910 cmd->print.expression = matrix_parse_exp (s);
4911 if (!cmd->print.expression)
4915 while (lex_match (s->lexer, T_SLASH))
4917 if (lex_match_id (s->lexer, "FORMAT"))
4919 lex_match (s->lexer, T_EQUALS);
4920 if (!parse_format_specifier (s->lexer, &cmd->print.format))
4922 cmd->print.use_default_format = false;
4924 else if (lex_match_id (s->lexer, "TITLE"))
4926 lex_match (s->lexer, T_EQUALS);
4927 if (!lex_force_string (s->lexer))
4929 free (cmd->print.title);
4930 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
4933 else if (lex_match_id (s->lexer, "SPACE"))
4935 lex_match (s->lexer, T_EQUALS);
4936 if (lex_match_id (s->lexer, "NEWPAGE"))
4937 cmd->print.space = -1;
4938 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
4940 cmd->print.space = lex_integer (s->lexer);
4946 else if (lex_match_id (s->lexer, "RLABELS"))
4947 parse_literal_print_labels (s, &cmd->print.rlabels);
4948 else if (lex_match_id (s->lexer, "CLABELS"))
4949 parse_literal_print_labels (s, &cmd->print.clabels);
4950 else if (lex_match_id (s->lexer, "RNAMES"))
4952 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
4955 else if (lex_match_id (s->lexer, "CNAMES"))
4957 if (!parse_expr_print_labels (s, &cmd->print.clabels))
4962 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
4963 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
4971 matrix_cmd_destroy (cmd);
4976 matrix_is_integer (const gsl_matrix *m)
4978 for (size_t y = 0; y < m->size1; y++)
4979 for (size_t x = 0; x < m->size2; x++)
4981 double d = gsl_matrix_get (m, y, x);
4989 matrix_max_magnitude (const gsl_matrix *m)
4992 for (size_t y = 0; y < m->size1; y++)
4993 for (size_t x = 0; x < m->size2; x++)
4995 double d = fabs (gsl_matrix_get (m, y, x));
5003 format_fits (struct fmt_spec format, double x)
5005 char *s = data_out (&(union value) { .f = x }, NULL,
5006 &format, settings_get_fmt_settings ());
5007 bool fits = *s != '*' && !strchr (s, 'E');
5012 static struct fmt_spec
5013 default_format (const gsl_matrix *m, int *log_scale)
5017 double max = matrix_max_magnitude (m);
5019 if (matrix_is_integer (m))
5020 for (int w = 1; w <= 12; w++)
5022 struct fmt_spec format = { .type = FMT_F, .w = w };
5023 if (format_fits (format, -max))
5027 if (max >= 1e9 || max <= 1e-4)
5030 snprintf (s, sizeof s, "%.1e", max);
5032 const char *e = strchr (s, 'e');
5034 *log_scale = atoi (e + 1);
5037 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5041 trimmed_string (double d)
5043 struct substring s = ss_buffer ((char *) &d, sizeof d);
5044 ss_rtrim (&s, ss_cstr (" "));
5045 return ss_xstrdup (s);
5048 static struct string_array *
5049 print_labels_get (const struct print_labels *labels, size_t n,
5050 const char *prefix, bool truncate)
5055 struct string_array *out = xzalloc (sizeof *out);
5056 if (labels->literals.n)
5057 string_array_clone (out, &labels->literals);
5058 else if (labels->expr)
5060 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5061 if (m && is_vector (m))
5063 gsl_vector v = to_vector (m);
5064 for (size_t i = 0; i < v.size; i++)
5065 string_array_append_nocopy (out, trimmed_string (
5066 gsl_vector_get (&v, i)));
5068 gsl_matrix_free (m);
5074 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5077 string_array_append (out, "");
5080 string_array_delete (out, out->n - 1);
5083 for (size_t i = 0; i < out->n; i++)
5085 char *s = out->strings[i];
5086 s[strnlen (s, 8)] = '\0';
5093 matrix_cmd_print_space (int space)
5096 output_item_submit (page_break_item_create ());
5097 for (int i = 0; i < space; i++)
5098 output_log ("%s", "");
5102 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
5103 struct fmt_spec format, int log_scale)
5105 matrix_cmd_print_space (print->space);
5107 output_log ("%s", print->title);
5109 output_log (" 10 ** %d X", log_scale);
5111 struct string_array *clabels = print_labels_get (print->clabels,
5112 m->size2, "col", true);
5113 if (clabels && format.w < 8)
5115 struct string_array *rlabels = print_labels_get (print->rlabels,
5116 m->size1, "row", true);
5120 struct string line = DS_EMPTY_INITIALIZER;
5122 ds_put_byte_multiple (&line, ' ', 8);
5123 for (size_t x = 0; x < m->size2; x++)
5124 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5125 output_log_nocopy (ds_steal_cstr (&line));
5128 double scale = pow (10.0, log_scale);
5129 bool numeric = fmt_is_numeric (format.type);
5130 for (size_t y = 0; y < m->size1; y++)
5132 struct string line = DS_EMPTY_INITIALIZER;
5134 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5136 for (size_t x = 0; x < m->size2; x++)
5138 double f = gsl_matrix_get (m, y, x);
5140 ? data_out (&(union value) { .f = f / scale}, NULL,
5141 &format, settings_get_fmt_settings ())
5142 : trimmed_string (f));
5143 ds_put_format (&line, " %s", s);
5146 output_log_nocopy (ds_steal_cstr (&line));
5149 string_array_destroy (rlabels);
5151 string_array_destroy (clabels);
5156 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5157 const struct print_labels *print_labels, size_t n,
5158 const char *name, const char *prefix)
5160 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5162 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5163 for (size_t i = 0; i < n; i++)
5164 pivot_category_create_leaf (
5166 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5167 : pivot_value_new_integer (i + 1)));
5169 d->hide_all_labels = true;
5170 string_array_destroy (labels);
5175 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
5176 struct fmt_spec format, int log_scale)
5178 struct pivot_table *table = pivot_table_create__ (
5179 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5182 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5184 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5185 N_("Columns"), "col");
5187 struct pivot_footnote *footnote = NULL;
5190 char *s = xasprintf ("× 10**%d", log_scale);
5191 footnote = pivot_table_create_footnote (
5192 table, pivot_value_new_user_text_nocopy (s));
5195 double scale = pow (10.0, log_scale);
5196 bool numeric = fmt_is_numeric (format.type);
5197 for (size_t y = 0; y < m->size1; y++)
5198 for (size_t x = 0; x < m->size2; x++)
5200 double f = gsl_matrix_get (m, y, x);
5201 struct pivot_value *v;
5204 v = pivot_value_new_number (f / scale);
5205 v->numeric.format = format;
5208 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5210 pivot_value_add_footnote (v, footnote);
5211 pivot_table_put2 (table, y, x, v);
5214 pivot_table_submit (table);
5218 matrix_cmd_execute_print (const struct print_command *print)
5220 if (print->expression)
5222 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5227 struct fmt_spec format = (print->use_default_format
5228 ? default_format (m, &log_scale)
5231 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5232 matrix_cmd_print_text (print, m, format, log_scale);
5234 matrix_cmd_print_tables (print, m, format, log_scale);
5236 gsl_matrix_free (m);
5240 matrix_cmd_print_space (print->space);
5242 output_log ("%s", print->title);
5249 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
5250 const char *command_name,
5251 const char *stop1, const char *stop2)
5253 lex_end_of_command (s->lexer);
5254 lex_discard_rest_of_command (s->lexer);
5256 size_t allocated = 0;
5259 while (lex_token (s->lexer) == T_ENDCMD)
5262 if (lex_at_phrase (s->lexer, stop1)
5263 || (stop2 && lex_at_phrase (s->lexer, stop2)))
5266 if (lex_at_phrase (s->lexer, "END MATRIX"))
5268 msg (SE, _("Premature END MATRIX within %s."), command_name);
5272 struct matrix_cmd *cmd = matrix_parse_command (s);
5276 if (c->n >= allocated)
5277 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
5278 c->commands[c->n++] = cmd;
5283 matrix_parse_do_if_clause (struct matrix_state *s,
5284 struct do_if_command *ifc,
5285 bool parse_condition,
5286 size_t *allocated_clauses)
5288 if (ifc->n_clauses >= *allocated_clauses)
5289 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5290 sizeof *ifc->clauses);
5291 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5292 *c = (struct do_if_clause) { .condition = NULL };
5294 if (parse_condition)
5296 c->condition = matrix_parse_expr (s);
5301 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
5304 static struct matrix_cmd *
5305 matrix_parse_do_if (struct matrix_state *s)
5307 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5308 *cmd = (struct matrix_cmd) {
5310 .do_if = { .n_clauses = 0 }
5313 size_t allocated_clauses = 0;
5316 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
5319 while (lex_match_phrase (s->lexer, "ELSE IF"));
5321 if (lex_match_id (s->lexer, "ELSE")
5322 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
5325 if (!lex_match_phrase (s->lexer, "END IF"))
5330 matrix_cmd_destroy (cmd);
5335 matrix_cmd_execute_do_if (struct do_if_command *cmd)
5337 for (size_t i = 0; i < cmd->n_clauses; i++)
5339 struct do_if_clause *c = &cmd->clauses[i];
5343 if (!matrix_expr_evaluate_scalar (c->condition,
5344 i ? "ELSE IF" : "DO IF",
5349 for (size_t j = 0; j < c->commands.n; j++)
5350 if (!matrix_cmd_execute (c->commands.commands[j]))
5357 static struct matrix_cmd *
5358 matrix_parse_loop (struct matrix_state *s)
5360 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5361 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5363 struct loop_command *loop = &cmd->loop;
5364 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5366 struct substring name = lex_tokss (s->lexer);
5367 loop->var = matrix_var_lookup (s, name);
5369 loop->var = matrix_var_create (s, name);
5374 loop->start = matrix_parse_expr (s);
5375 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5378 loop->end = matrix_parse_expr (s);
5382 if (lex_match (s->lexer, T_BY))
5384 loop->increment = matrix_parse_expr (s);
5385 if (!loop->increment)
5390 if (lex_match_id (s->lexer, "IF"))
5392 loop->top_condition = matrix_parse_expr (s);
5393 if (!loop->top_condition)
5397 bool was_in_loop = s->in_loop;
5399 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
5401 s->in_loop = was_in_loop;
5405 if (!lex_match_phrase (s->lexer, "END LOOP"))
5408 if (lex_match_id (s->lexer, "IF"))
5410 loop->bottom_condition = matrix_parse_expr (s);
5411 if (!loop->bottom_condition)
5418 matrix_cmd_destroy (cmd);
5423 matrix_cmd_execute_loop (struct loop_command *cmd)
5427 long int increment = 1;
5430 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5431 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5433 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5437 if (increment > 0 ? value > end
5438 : increment < 0 ? value < end
5443 int mxloops = settings_get_mxloops ();
5444 for (int i = 0; i < mxloops; i++)
5449 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5451 gsl_matrix_free (cmd->var->value);
5452 cmd->var->value = NULL;
5454 if (!cmd->var->value)
5455 cmd->var->value = gsl_matrix_alloc (1, 1);
5457 gsl_matrix_set (cmd->var->value, 0, 0, value);
5460 if (cmd->top_condition)
5463 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5468 for (size_t j = 0; j < cmd->commands.n; j++)
5469 if (!matrix_cmd_execute (cmd->commands.commands[j]))
5472 if (cmd->bottom_condition)
5475 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5476 "END LOOP IF", &d) || d > 0)
5483 if (increment > 0 ? value > end : value < end)
5489 static struct matrix_cmd *
5490 matrix_parse_break (struct matrix_state *s)
5494 msg (SE, _("BREAK not inside LOOP."));
5498 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5499 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
5503 static struct matrix_cmd *
5504 matrix_parse_display (struct matrix_state *s)
5506 while (lex_token (s->lexer) != T_ENDCMD)
5508 if (!lex_match_id (s->lexer, "DICTIONARY")
5509 && !lex_match_id (s->lexer, "STATUS"))
5511 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5516 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5517 *cmd = (struct matrix_cmd) { .type = MCMD_DISPLAY, .display = { s } };
5522 compare_matrix_var_pointers (const void *a_, const void *b_)
5524 const struct matrix_var *const *ap = a_;
5525 const struct matrix_var *const *bp = b_;
5526 const struct matrix_var *a = *ap;
5527 const struct matrix_var *b = *bp;
5528 return strcmp (a->name, b->name);
5532 matrix_cmd_execute_display (struct display_command *cmd)
5534 const struct matrix_state *s = cmd->state;
5536 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5537 struct pivot_dimension *attr_dimension
5538 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5539 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5540 N_("Rows"), N_("Columns"));
5541 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5543 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5546 const struct matrix_var *var;
5547 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5549 vars[n_vars++] = var;
5550 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5552 struct pivot_dimension *rows = pivot_dimension_create (
5553 table, PIVOT_AXIS_ROW, N_("Variable"));
5554 for (size_t i = 0; i < n_vars; i++)
5556 const struct matrix_var *var = vars[i];
5557 pivot_category_create_leaf (
5558 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
5560 size_t r = var->value->size1;
5561 size_t c = var->value->size2;
5562 double values[] = { r, c, r * c * sizeof (double) / 1024 };
5563 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
5564 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
5567 pivot_table_submit (table);
5570 static struct matrix_cmd *
5571 matrix_parse_release (struct matrix_state *s)
5573 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5574 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
5576 size_t allocated_vars = 0;
5577 while (lex_token (s->lexer) == T_ID)
5579 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
5582 if (cmd->release.n_vars >= allocated_vars)
5583 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
5584 sizeof *cmd->release.vars);
5585 cmd->release.vars[cmd->release.n_vars++] = var;
5588 lex_error (s->lexer, _("Variable name expected."));
5591 if (!lex_match (s->lexer, T_COMMA))
5599 matrix_cmd_execute_release (struct release_command *cmd)
5601 for (size_t i = 0; i < cmd->n_vars; i++)
5603 struct matrix_var *var = cmd->vars[i];
5604 gsl_matrix_free (var->value);
5609 static struct save_file *
5610 save_file_create (struct matrix_state *s, struct file_handle *fh,
5611 struct string_array *variables,
5612 struct matrix_expr *names,
5613 struct stringi_set *strings)
5615 for (size_t i = 0; i < s->n_save_files; i++)
5617 struct save_file *sf = s->save_files[i];
5618 if (fh_equal (sf->file, fh))
5622 string_array_destroy (variables);
5623 matrix_expr_destroy (names);
5624 stringi_set_destroy (strings);
5630 struct save_file *sf = xmalloc (sizeof *sf);
5631 *sf = (struct save_file) {
5633 .dataset = s->dataset,
5634 .variables = *variables,
5636 .strings = STRINGI_SET_INITIALIZER (sf->strings),
5639 stringi_set_swap (&sf->strings, strings);
5640 stringi_set_destroy (strings);
5642 s->save_files = xrealloc (s->save_files,
5643 (s->n_save_files + 1) * sizeof *s->save_files);
5644 s->save_files[s->n_save_files++] = sf;
5648 static struct casewriter *
5649 save_file_open (struct save_file *sf, gsl_matrix *m)
5651 if (sf->writer || sf->error)
5655 size_t n_variables = caseproto_get_n_widths (
5656 casewriter_get_proto (sf->writer));
5657 if (m->size2 != n_variables)
5659 msg (ME, _("The first SAVE to %s within this matrix program "
5660 "had %zu columns, so a %zu×%zu matrix cannot be "
5662 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
5663 n_variables, m->size1, m->size2);
5671 struct dictionary *dict = dict_create (get_default_encoding ());
5673 /* Fill 'names' with user-specified names if there were any, either from
5674 sf->variables or sf->names. */
5675 struct string_array names = { .n = 0 };
5676 if (sf->variables.n)
5677 string_array_clone (&names, &sf->variables);
5680 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
5681 if (nm && is_vector (nm))
5683 gsl_vector nv = to_vector (nm);
5684 for (size_t i = 0; i < nv.size; i++)
5686 char *name = trimmed_string (gsl_vector_get (&nv, i));
5687 if (dict_id_is_valid (dict, name, true))
5688 string_array_append_nocopy (&names, name);
5693 gsl_matrix_free (nm);
5696 struct stringi_set strings;
5697 stringi_set_clone (&strings, &sf->strings);
5699 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
5704 name = names.strings[i];
5707 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
5711 int width = stringi_set_delete (&strings, name) ? 8 : 0;
5712 struct variable *var = dict_create_var (dict, name, width);
5715 msg (ME, _("Duplicate variable name %s in SAVE statement."),
5721 if (!stringi_set_is_empty (&strings))
5723 size_t count = stringi_set_count (&strings);
5724 const char *example = stringi_set_node_get_string (
5725 stringi_set_first (&strings));
5728 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
5729 "variable %s."), example);
5731 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
5732 "unknown variable: %s.",
5733 "The SAVE command STRINGS subcommand specifies %zu "
5734 "unknown variables, including %s.",
5739 stringi_set_destroy (&strings);
5740 string_array_destroy (&names);
5749 if (sf->file == fh_inline_file ())
5750 sf->writer = autopaging_writer_create (dict_get_proto (dict));
5752 sf->writer = any_writer_open (sf->file, dict);
5764 save_file_destroy (struct save_file *sf)
5768 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
5770 dataset_set_dict (sf->dataset, sf->dict);
5771 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
5775 casewriter_destroy (sf->writer);
5776 dict_unref (sf->dict);
5778 fh_unref (sf->file);
5779 string_array_destroy (&sf->variables);
5780 matrix_expr_destroy (sf->names);
5781 stringi_set_destroy (&sf->strings);
5786 static struct matrix_cmd *
5787 matrix_parse_save (struct matrix_state *s)
5789 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5790 *cmd = (struct matrix_cmd) {
5792 .save = { .expression = NULL },
5795 struct file_handle *fh = NULL;
5796 struct save_command *save = &cmd->save;
5798 struct string_array variables = STRING_ARRAY_INITIALIZER;
5799 struct matrix_expr *names = NULL;
5800 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
5802 save->expression = matrix_parse_exp (s);
5803 if (!save->expression)
5806 while (lex_match (s->lexer, T_SLASH))
5808 if (lex_match_id (s->lexer, "OUTFILE"))
5810 lex_match (s->lexer, T_EQUALS);
5813 fh = (lex_match (s->lexer, T_ASTERISK)
5815 : fh_parse (s->lexer, FH_REF_FILE, s->session));
5819 else if (lex_match_id (s->lexer, "VARIABLES"))
5821 lex_match (s->lexer, T_EQUALS);
5825 struct dictionary *d = dict_create (get_default_encoding ());
5826 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
5827 PV_NO_SCRATCH | PV_NO_DUPLICATE);
5832 string_array_clear (&variables);
5833 variables = (struct string_array) {
5839 else if (lex_match_id (s->lexer, "NAMES"))
5841 lex_match (s->lexer, T_EQUALS);
5842 matrix_expr_destroy (names);
5843 names = matrix_parse_exp (s);
5847 else if (lex_match_id (s->lexer, "STRINGS"))
5849 lex_match (s->lexer, T_EQUALS);
5850 while (lex_token (s->lexer) == T_ID)
5852 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
5854 if (!lex_match (s->lexer, T_COMMA))
5860 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
5868 if (s->prev_save_file)
5869 fh = fh_ref (s->prev_save_file);
5872 lex_sbc_missing ("OUTFILE");
5876 fh_unref (s->prev_save_file);
5877 s->prev_save_file = fh_ref (fh);
5879 if (variables.n && names)
5881 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
5882 matrix_expr_destroy (names);
5886 save->sf = save_file_create (s, fh, &variables, names, &strings);
5890 string_array_destroy (&variables);
5891 matrix_expr_destroy (names);
5892 stringi_set_destroy (&strings);
5894 matrix_cmd_destroy (cmd);
5899 matrix_cmd_execute_save (const struct save_command *save)
5901 gsl_matrix *m = matrix_expr_evaluate (save->expression);
5905 struct casewriter *writer = save_file_open (save->sf, m);
5908 gsl_matrix_free (m);
5912 const struct caseproto *proto = casewriter_get_proto (writer);
5913 for (size_t y = 0; y < m->size1; y++)
5915 struct ccase *c = case_create (proto);
5916 for (size_t x = 0; x < m->size2; x++)
5918 double d = gsl_matrix_get (m, y, x);
5919 int width = caseproto_get_width (proto, x);
5920 union value *value = case_data_rw_idx (c, x);
5924 memcpy (value->s, &d, width);
5926 casewriter_write (writer, c);
5928 gsl_matrix_free (m);
5931 static struct matrix_cmd *
5932 matrix_parse_read (struct matrix_state *s)
5934 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5935 *cmd = (struct matrix_cmd) {
5937 .read = { .format = FMT_F },
5940 struct file_handle *fh = NULL;
5941 char *encoding = NULL;
5942 struct read_command *read = &cmd->read;
5943 read->dst = matrix_lvalue_parse (s);
5948 int repetitions = 0;
5949 int record_width = 0;
5950 bool seen_format = false;
5951 while (lex_match (s->lexer, T_SLASH))
5953 if (lex_match_id (s->lexer, "FILE"))
5955 lex_match (s->lexer, T_EQUALS);
5958 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5962 else if (lex_match_id (s->lexer, "ENCODING"))
5964 lex_match (s->lexer, T_EQUALS);
5965 if (!lex_force_string (s->lexer))
5969 encoding = ss_xstrdup (lex_tokss (s->lexer));
5973 else if (lex_match_id (s->lexer, "FIELD"))
5975 lex_match (s->lexer, T_EQUALS);
5977 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5979 read->c1 = lex_integer (s->lexer);
5981 if (!lex_force_match (s->lexer, T_TO)
5982 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
5984 read->c2 = lex_integer (s->lexer) + 1;
5987 record_width = read->c2 - read->c1;
5988 if (lex_match (s->lexer, T_BY))
5990 if (!lex_force_int_range (s->lexer, "BY", 1,
5991 read->c2 - read->c1))
5993 by = lex_integer (s->lexer);
5996 if (record_width % by)
5998 msg (SE, _("BY %d does not evenly divide record width %d."),
6006 else if (lex_match_id (s->lexer, "SIZE"))
6008 lex_match (s->lexer, T_EQUALS);
6009 matrix_expr_destroy (read->size);
6010 read->size = matrix_parse_exp (s);
6014 else if (lex_match_id (s->lexer, "MODE"))
6016 lex_match (s->lexer, T_EQUALS);
6017 if (lex_match_id (s->lexer, "RECTANGULAR"))
6018 read->symmetric = false;
6019 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6020 read->symmetric = true;
6023 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6027 else if (lex_match_id (s->lexer, "REREAD"))
6028 read->reread = true;
6029 else if (lex_match_id (s->lexer, "FORMAT"))
6033 lex_sbc_only_once ("FORMAT");
6038 lex_match (s->lexer, T_EQUALS);
6040 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6043 const char *p = lex_tokcstr (s->lexer);
6044 if (c_isdigit (p[0]))
6046 repetitions = atoi (p);
6047 p += strspn (p, "0123456789");
6048 if (!fmt_from_name (p, &read->format))
6050 lex_error (s->lexer, _("Unknown format %s."), p);
6055 else if (fmt_from_name (p, &read->format))
6059 struct fmt_spec format;
6060 if (!parse_format_specifier (s->lexer, &format))
6062 read->format = format.type;
6068 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6069 "REREAD", "FORMAT");
6076 lex_sbc_missing ("FIELD");
6080 if (!read->dst->n_indexes && !read->size)
6082 msg (SE, _("SIZE is required for reading data into a full matrix "
6083 "(as opposed to a submatrix)."));
6089 if (s->prev_read_file)
6090 fh = fh_ref (s->prev_read_file);
6093 lex_sbc_missing ("FILE");
6097 fh_unref (s->prev_read_file);
6098 s->prev_read_file = fh_ref (fh);
6100 read->rf = read_file_create (s, fh);
6104 free (read->rf->encoding);
6105 read->rf->encoding = encoding;
6109 /* Field width may be specified in multiple ways:
6112 2. The format on FORMAT.
6113 3. The repetition factor on FORMAT.
6115 (2) and (3) are mutually exclusive.
6117 If more than one of these is present, they must agree. If none of them is
6118 present, then free-field format is used.
6120 if (repetitions > record_width)
6122 msg (SE, _("%d repetitions cannot fit in record width %d."),
6123 repetitions, record_width);
6126 int w = (repetitions ? record_width / repetitions
6132 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6133 "which implies field width %d, "
6134 "but BY specifies field width %d."),
6135 repetitions, record_width, w, by);
6137 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6146 matrix_cmd_destroy (cmd);
6152 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6153 struct substring data, size_t y, size_t x,
6154 int first_column, int last_column, char *error)
6156 int line_number = dfm_get_line_number (reader);
6157 struct msg_location *location = xmalloc (sizeof *location);
6158 *location = (struct msg_location) {
6159 .file_name = xstrdup (dfm_get_file_name (reader)),
6160 .first_line = line_number,
6161 .last_line = line_number + 1,
6162 .first_column = first_column,
6163 .last_column = last_column,
6165 struct msg *m = xmalloc (sizeof *m);
6167 .category = MSG_C_DATA,
6168 .severity = MSG_S_WARNING,
6169 .location = location,
6170 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
6171 "for matrix row %zu, column %zu: %s"),
6172 (int) data.length, data.string, fmt_name (format),
6173 y + 1, x + 1, error),
6181 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
6182 gsl_matrix *m, struct substring p, size_t y, size_t x,
6183 const char *line_start)
6185 const char *input_encoding = dfm_reader_get_encoding (reader);
6188 if (fmt_is_numeric (read->format))
6191 error = data_in (p, input_encoding, read->format,
6192 settings_get_fmt_settings (), &v, 0, NULL);
6193 if (!error && v.f == SYSMIS)
6194 error = xstrdup (_("Matrix data may not contain missing value."));
6199 uint8_t s[sizeof (double)];
6200 union value v = { .s = s };
6201 error = data_in (p, input_encoding, read->format,
6202 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6203 memcpy (&f, s, sizeof f);
6208 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6209 int c2 = c1 + ss_utf8_count_columns (p) - 1;
6210 parse_error (reader, read->format, p, y, x, c1, c2, error);
6214 gsl_matrix_set (m, y, x, f);
6215 if (read->symmetric && x != y)
6216 gsl_matrix_set (m, x, y, f);
6221 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
6222 struct substring *line, const char **startp)
6224 if (dfm_eof (reader))
6226 msg (SE, _("Unexpected end of file reading matrix data."));
6229 dfm_expand_tabs (reader);
6230 struct substring record = dfm_get_record (reader);
6231 /* XXX need to recode record into UTF-8 */
6232 *startp = record.string;
6233 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6238 matrix_read (struct read_command *read, struct dfm_reader *reader,
6241 for (size_t y = 0; y < m->size1; y++)
6243 size_t nx = read->symmetric ? y + 1 : m->size2;
6245 struct substring line = ss_empty ();
6246 const char *line_start = line.string;
6247 for (size_t x = 0; x < nx; x++)
6254 ss_ltrim (&line, ss_cstr (" ,"));
6255 if (!ss_is_empty (line))
6257 if (!matrix_read_line (read, reader, &line, &line_start))
6259 dfm_forward_record (reader);
6262 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6266 if (!matrix_read_line (read, reader, &line, &line_start))
6268 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6269 int f = x % fields_per_line;
6270 if (f == fields_per_line - 1)
6271 dfm_forward_record (reader);
6273 p = ss_substr (line, read->w * f, read->w);
6276 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6280 dfm_forward_record (reader);
6283 ss_ltrim (&line, ss_cstr (" ,"));
6284 if (!ss_is_empty (line))
6287 msg (SW, _("Trailing garbage on line \"%.*s\""),
6288 (int) line.length, line.string);
6295 matrix_cmd_execute_read (struct read_command *read)
6297 struct index_vector iv0, iv1;
6298 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6301 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6304 gsl_matrix *m = matrix_expr_evaluate (read->size);
6310 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6311 "not a %zu×%zu matrix."), m->size1, m->size2);
6312 gsl_matrix_free (m);
6318 gsl_vector v = to_vector (m);
6322 d[0] = gsl_vector_get (&v, 0);
6325 else if (v.size == 2)
6327 d[0] = gsl_vector_get (&v, 0);
6328 d[1] = gsl_vector_get (&v, 1);
6332 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6333 "not a %zu×%zu matrix."),
6334 m->size1, m->size2),
6335 gsl_matrix_free (m);
6340 gsl_matrix_free (m);
6342 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6344 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
6345 "are outside valid range."),
6356 if (read->dst->n_indexes)
6358 size_t submatrix_size[2];
6359 if (read->dst->n_indexes == 2)
6361 submatrix_size[0] = iv0.n;
6362 submatrix_size[1] = iv1.n;
6364 else if (read->dst->var->value->size1 == 1)
6366 submatrix_size[0] = 1;
6367 submatrix_size[1] = iv0.n;
6371 submatrix_size[0] = iv0.n;
6372 submatrix_size[1] = 1;
6377 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6379 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
6380 "differ from submatrix dimensions %zu×%zu."),
6382 submatrix_size[0], submatrix_size[1]);
6390 size[0] = submatrix_size[0];
6391 size[1] = submatrix_size[1];
6395 struct dfm_reader *reader = read_file_open (read->rf);
6397 dfm_reread_record (reader, 1);
6399 if (read->symmetric && size[0] != size[1])
6401 msg (SE, _("Cannot read non-square %zu×%zu matrix "
6402 "using READ with MODE=SYMMETRIC."),
6408 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6409 matrix_read (read, reader, tmp);
6410 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
6413 static struct matrix_cmd *
6414 matrix_parse_write (struct matrix_state *s)
6416 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6417 *cmd = (struct matrix_cmd) {
6421 struct file_handle *fh = NULL;
6422 char *encoding = NULL;
6423 struct write_command *write = &cmd->write;
6424 write->expression = matrix_parse_exp (s);
6425 if (!write->expression)
6429 int repetitions = 0;
6430 int record_width = 0;
6431 enum fmt_type format = FMT_F;
6432 bool has_format = false;
6433 while (lex_match (s->lexer, T_SLASH))
6435 if (lex_match_id (s->lexer, "OUTFILE"))
6437 lex_match (s->lexer, T_EQUALS);
6440 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6444 else if (lex_match_id (s->lexer, "ENCODING"))
6446 lex_match (s->lexer, T_EQUALS);
6447 if (!lex_force_string (s->lexer))
6451 encoding = ss_xstrdup (lex_tokss (s->lexer));
6455 else if (lex_match_id (s->lexer, "FIELD"))
6457 lex_match (s->lexer, T_EQUALS);
6459 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6461 write->c1 = lex_integer (s->lexer);
6463 if (!lex_force_match (s->lexer, T_TO)
6464 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
6466 write->c2 = lex_integer (s->lexer) + 1;
6469 record_width = write->c2 - write->c1;
6470 if (lex_match (s->lexer, T_BY))
6472 if (!lex_force_int_range (s->lexer, "BY", 1,
6473 write->c2 - write->c1))
6475 by = lex_integer (s->lexer);
6478 if (record_width % by)
6480 msg (SE, _("BY %d does not evenly divide record width %d."),
6488 else if (lex_match_id (s->lexer, "MODE"))
6490 lex_match (s->lexer, T_EQUALS);
6491 if (lex_match_id (s->lexer, "RECTANGULAR"))
6492 write->triangular = false;
6493 else if (lex_match_id (s->lexer, "TRIANGULAR"))
6494 write->triangular = true;
6497 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
6501 else if (lex_match_id (s->lexer, "HOLD"))
6503 else if (lex_match_id (s->lexer, "FORMAT"))
6505 if (has_format || write->format)
6507 lex_sbc_only_once ("FORMAT");
6511 lex_match (s->lexer, T_EQUALS);
6513 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6516 const char *p = lex_tokcstr (s->lexer);
6517 if (c_isdigit (p[0]))
6519 repetitions = atoi (p);
6520 p += strspn (p, "0123456789");
6521 if (!fmt_from_name (p, &format))
6523 lex_error (s->lexer, _("Unknown format %s."), p);
6529 else if (fmt_from_name (p, &format))
6536 struct fmt_spec spec;
6537 if (!parse_format_specifier (s->lexer, &spec))
6539 write->format = xmemdup (&spec, sizeof spec);
6544 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
6552 lex_sbc_missing ("FIELD");
6558 if (s->prev_write_file)
6559 fh = fh_ref (s->prev_write_file);
6562 lex_sbc_missing ("OUTFILE");
6566 fh_unref (s->prev_write_file);
6567 s->prev_write_file = fh_ref (fh);
6569 write->wf = write_file_create (s, fh);
6573 free (write->wf->encoding);
6574 write->wf->encoding = encoding;
6578 /* Field width may be specified in multiple ways:
6581 2. The format on FORMAT.
6582 3. The repetition factor on FORMAT.
6584 (2) and (3) are mutually exclusive.
6586 If more than one of these is present, they must agree. If none of them is
6587 present, then free-field format is used.
6589 if (repetitions > record_width)
6591 msg (SE, _("%d repetitions cannot fit in record width %d."),
6592 repetitions, record_width);
6595 int w = (repetitions ? record_width / repetitions
6596 : write->format ? write->format->w
6601 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6602 "which implies field width %d, "
6603 "but BY specifies field width %d."),
6604 repetitions, record_width, w, by);
6606 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6610 if (w && !write->format)
6612 write->format = xmalloc (sizeof *write->format);
6613 *write->format = (struct fmt_spec) { .type = format, .w = w };
6615 if (!fmt_check_output (write->format))
6619 if (write->format && fmt_var_width (write->format) > sizeof (double))
6621 char s[FMT_STRING_LEN_MAX + 1];
6622 fmt_to_string (write->format, s);
6623 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
6624 s, sizeof (double));
6632 matrix_cmd_destroy (cmd);
6637 matrix_cmd_execute_write (struct write_command *write)
6639 gsl_matrix *m = matrix_expr_evaluate (write->expression);
6643 if (write->triangular && m->size1 != m->size2)
6645 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
6646 "the matrix to be written has dimensions %zu×%zu."),
6647 m->size1, m->size2);
6648 gsl_matrix_free (m);
6652 struct dfm_writer *writer = write_file_open (write->wf);
6653 if (!writer || !m->size1)
6655 gsl_matrix_free (m);
6659 const struct fmt_settings *settings = settings_get_fmt_settings ();
6660 struct u8_line *line = write->wf->held;
6661 for (size_t y = 0; y < m->size1; y++)
6665 line = xmalloc (sizeof *line);
6666 u8_line_init (line);
6668 size_t nx = write->triangular ? y + 1 : m->size2;
6670 for (size_t x = 0; x < nx; x++)
6673 double f = gsl_matrix_get (m, y, x);
6677 if (fmt_is_numeric (write->format->type))
6680 v.s = (uint8_t *) &f;
6681 s = data_out (&v, NULL, write->format, settings);
6685 s = xmalloc (DBL_BUFSIZE_BOUND);
6686 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
6687 >= DBL_BUFSIZE_BOUND)
6690 size_t len = strlen (s);
6691 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
6692 if (width + x0 > write->c2)
6694 dfm_put_record_utf8 (writer, line->s.ss.string,
6696 u8_line_clear (line);
6699 u8_line_put (line, x0, x0 + width, s, len);
6702 x0 += write->format ? write->format->w : width + 1;
6705 if (y + 1 >= m->size1 && write->hold)
6707 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
6708 u8_line_clear (line);
6712 u8_line_destroy (line);
6716 write->wf->held = line;
6718 gsl_matrix_free (m);
6721 static struct matrix_cmd *
6722 matrix_parse_get (struct matrix_state *s)
6724 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6725 *cmd = (struct matrix_cmd) {
6728 .dataset = s->dataset,
6729 .user = { .treatment = MGET_ERROR },
6730 .system = { .treatment = MGET_ERROR },
6734 struct get_command *get = &cmd->get;
6735 get->dst = matrix_lvalue_parse (s);
6739 while (lex_match (s->lexer, T_SLASH))
6741 if (lex_match_id (s->lexer, "FILE"))
6743 lex_match (s->lexer, T_EQUALS);
6745 fh_unref (get->file);
6746 if (lex_match (s->lexer, T_ASTERISK))
6750 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6755 else if (lex_match_id (s->lexer, "ENCODING"))
6757 lex_match (s->lexer, T_EQUALS);
6758 if (!lex_force_string (s->lexer))
6761 free (get->encoding);
6762 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
6766 else if (lex_match_id (s->lexer, "VARIABLES"))
6768 lex_match (s->lexer, T_EQUALS);
6772 lex_sbc_only_once ("VARIABLES");
6776 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
6779 else if (lex_match_id (s->lexer, "NAMES"))
6781 lex_match (s->lexer, T_EQUALS);
6782 if (!lex_force_id (s->lexer))
6785 struct substring name = lex_tokss (s->lexer);
6786 get->names = matrix_var_lookup (s, name);
6788 get->names = matrix_var_create (s, name);
6791 else if (lex_match_id (s->lexer, "MISSING"))
6793 lex_match (s->lexer, T_EQUALS);
6794 if (lex_match_id (s->lexer, "ACCEPT"))
6795 get->user.treatment = MGET_ACCEPT;
6796 else if (lex_match_id (s->lexer, "OMIT"))
6797 get->user.treatment = MGET_OMIT;
6798 else if (lex_is_number (s->lexer))
6800 get->user.treatment = MGET_RECODE;
6801 get->user.substitute = lex_number (s->lexer);
6806 lex_error (s->lexer, NULL);
6810 else if (lex_match_id (s->lexer, "SYSMIS"))
6812 lex_match (s->lexer, T_EQUALS);
6813 if (lex_match_id (s->lexer, "OMIT"))
6814 get->system.treatment = MGET_OMIT;
6815 else if (lex_is_number (s->lexer))
6817 get->system.treatment = MGET_RECODE;
6818 get->system.substitute = lex_number (s->lexer);
6823 lex_error (s->lexer, NULL);
6829 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
6830 "MISSING", "SYSMIS");
6835 if (get->user.treatment != MGET_ACCEPT)
6836 get->system.treatment = MGET_ERROR;
6841 matrix_cmd_destroy (cmd);
6846 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
6847 const struct dictionary *dict)
6849 struct variable **vars;
6854 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
6855 &vars, &n_vars, PV_NUMERIC))
6860 n_vars = dict_get_var_cnt (dict);
6861 vars = xnmalloc (n_vars, sizeof *vars);
6862 for (size_t i = 0; i < n_vars; i++)
6864 struct variable *var = dict_get_var (dict, i);
6865 if (!var_is_numeric (var))
6867 msg (SE, _("GET: Variable %s is not numeric."),
6868 var_get_name (var));
6878 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
6879 for (size_t i = 0; i < n_vars; i++)
6881 char s[sizeof (double)];
6883 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
6884 memcpy (&f, s, sizeof f);
6885 gsl_matrix_set (names, i, 0, f);
6888 gsl_matrix_free (get->names->value);
6889 get->names->value = names;
6893 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
6894 long long int casenum = 1;
6896 for (struct ccase *c = casereader_read (reader); c;
6897 c = casereader_read (reader), casenum++)
6899 if (n_rows >= m->size1)
6901 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
6902 for (size_t y = 0; y < n_rows; y++)
6903 for (size_t x = 0; x < n_vars; x++)
6904 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
6905 gsl_matrix_free (m);
6910 for (size_t x = 0; x < n_vars; x++)
6912 const struct variable *var = vars[x];
6913 double d = case_num (c, var);
6916 if (get->system.treatment == MGET_RECODE)
6917 d = get->system.substitute;
6918 else if (get->system.treatment == MGET_OMIT)
6922 msg (SE, _("GET: Variable %s in case %lld "
6923 "is system-missing."),
6924 var_get_name (var), casenum);
6928 else if (var_is_num_missing (var, d, MV_USER))
6930 if (get->user.treatment == MGET_RECODE)
6931 d = get->user.substitute;
6932 else if (get->user.treatment == MGET_OMIT)
6934 else if (get->user.treatment != MGET_ACCEPT)
6936 msg (SE, _("GET: Variable %s in case %lld has user-missing "
6938 var_get_name (var), casenum, d);
6942 gsl_matrix_set (m, n_rows, x, d);
6953 matrix_lvalue_evaluate_and_assign (get->dst, m);
6956 gsl_matrix_free (m);
6961 matrix_open_casereader (const char *command_name,
6962 struct file_handle *file, const char *encoding,
6963 struct dataset *dataset,
6964 struct casereader **readerp, struct dictionary **dictp)
6968 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
6969 return *readerp != NULL;
6973 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
6975 msg (ME, _("The %s command cannot read an empty active file."),
6979 *readerp = proc_open (dataset);
6980 *dictp = dict_ref (dataset_dict (dataset));
6986 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
6987 struct casereader *reader, struct dictionary *dict)
6990 casereader_destroy (reader);
6992 proc_commit (dataset);
6996 matrix_cmd_execute_get (struct get_command *get)
6998 struct casereader *r;
6999 struct dictionary *d;
7000 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
7003 matrix_cmd_execute_get__ (get, r, d);
7004 matrix_close_casereader (get->file, get->dataset, r, d);
7009 match_rowtype (struct lexer *lexer)
7011 static const char *rowtypes[] = {
7012 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7014 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7016 for (size_t i = 0; i < n_rowtypes; i++)
7017 if (lex_match_id (lexer, rowtypes[i]))
7020 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7025 parse_var_names (struct lexer *lexer, struct string_array *sa)
7027 lex_match (lexer, T_EQUALS);
7029 string_array_clear (sa);
7031 struct dictionary *dict = dict_create (get_default_encoding ());
7034 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7035 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7040 for (size_t i = 0; i < n_names; i++)
7041 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7042 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7044 msg (SE, _("Variable name %s is reserved."), names[i]);
7045 for (size_t j = 0; j < n_names; j++)
7051 string_array_clear (sa);
7052 sa->strings = names;
7053 sa->n = sa->allocated = n_names;
7059 msave_common_destroy (struct msave_common *common)
7063 fh_unref (common->outfile);
7064 string_array_destroy (&common->variables);
7065 string_array_destroy (&common->fnames);
7066 string_array_destroy (&common->snames);
7068 for (size_t i = 0; i < common->n_factors; i++)
7069 matrix_expr_destroy (common->factors[i]);
7070 free (common->factors);
7072 for (size_t i = 0; i < common->n_splits; i++)
7073 matrix_expr_destroy (common->splits[i]);
7074 free (common->splits);
7076 dict_unref (common->dict);
7077 casewriter_destroy (common->writer);
7084 compare_variables (const char *keyword,
7085 const struct string_array *new,
7086 const struct string_array *old)
7092 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7093 "on the first MSAVE within MATRIX."), keyword);
7096 else if (!string_array_equal_case (old, new))
7098 msg (SE, _("%s must specify the same variables each time within "
7099 "a given MATRIX."), keyword);
7106 static struct matrix_cmd *
7107 matrix_parse_msave (struct matrix_state *s)
7109 struct msave_common *common = xmalloc (sizeof *common);
7110 *common = (struct msave_common) { .outfile = NULL };
7112 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7113 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7115 struct matrix_expr *splits = NULL;
7116 struct matrix_expr *factors = NULL;
7118 struct msave_command *msave = &cmd->msave;
7119 msave->expr = matrix_parse_exp (s);
7123 while (lex_match (s->lexer, T_SLASH))
7125 if (lex_match_id (s->lexer, "TYPE"))
7127 lex_match (s->lexer, T_EQUALS);
7129 msave->rowtype = match_rowtype (s->lexer);
7130 if (!msave->rowtype)
7133 else if (lex_match_id (s->lexer, "OUTFILE"))
7135 lex_match (s->lexer, T_EQUALS);
7137 fh_unref (common->outfile);
7138 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7139 if (!common->outfile)
7142 else if (lex_match_id (s->lexer, "VARIABLES"))
7144 if (!parse_var_names (s->lexer, &common->variables))
7147 else if (lex_match_id (s->lexer, "FNAMES"))
7149 if (!parse_var_names (s->lexer, &common->fnames))
7152 else if (lex_match_id (s->lexer, "SNAMES"))
7154 if (!parse_var_names (s->lexer, &common->snames))
7157 else if (lex_match_id (s->lexer, "SPLIT"))
7159 lex_match (s->lexer, T_EQUALS);
7161 matrix_expr_destroy (splits);
7162 splits = matrix_parse_exp (s);
7166 else if (lex_match_id (s->lexer, "FACTOR"))
7168 lex_match (s->lexer, T_EQUALS);
7170 matrix_expr_destroy (factors);
7171 factors = matrix_parse_exp (s);
7177 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7178 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7182 if (!msave->rowtype)
7184 lex_sbc_missing ("TYPE");
7190 if (common->fnames.n && !factors)
7192 msg (SE, _("FNAMES requires FACTOR."));
7195 if (common->snames.n && !splits)
7197 msg (SE, _("SNAMES requires SPLIT."));
7200 if (!common->outfile)
7202 lex_sbc_missing ("OUTFILE");
7209 if (common->outfile)
7211 if (!fh_equal (common->outfile, s->common->outfile))
7213 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7214 "within a single MATRIX command."));
7218 if (!compare_variables ("VARIABLES",
7219 &common->variables, &s->common->variables)
7220 || !compare_variables ("FNAMES", &common->fnames, &s->common->fnames)
7221 || !compare_variables ("SNAMES", &common->snames, &s->common->snames))
7223 msave_common_destroy (common);
7225 msave->common = s->common;
7227 struct msave_common *c = s->common;
7230 if (c->n_factors >= c->allocated_factors)
7231 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7232 sizeof *c->factors);
7233 c->factors[c->n_factors++] = factors;
7235 if (c->n_factors > 0)
7236 msave->factors = c->factors[c->n_factors - 1];
7240 if (c->n_splits >= c->allocated_splits)
7241 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7243 c->splits[c->n_splits++] = splits;
7245 if (c->n_splits > 0)
7246 msave->splits = c->splits[c->n_splits - 1];
7251 matrix_expr_destroy (splits);
7252 matrix_expr_destroy (factors);
7253 msave_common_destroy (common);
7254 matrix_cmd_destroy (cmd);
7259 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7261 gsl_matrix *m = matrix_expr_evaluate (e);
7267 msg (SE, _("%s expression must evaluate to vector, "
7268 "not a %zu×%zu matrix."),
7269 name, m->size1, m->size2);
7270 gsl_matrix_free (m);
7274 return matrix_to_vector (m);
7278 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7280 for (size_t i = 0; i < vars->n; i++)
7281 if (!dict_create_var (d, vars->strings[i], 0))
7282 return vars->strings[i];
7286 static struct dictionary *
7287 msave_create_dict (const struct msave_common *common)
7289 struct dictionary *dict = dict_create (get_default_encoding ());
7291 const char *dup_split = msave_add_vars (dict, &common->snames);
7294 /* Should not be possible because the parser ensures that the names are
7299 dict_create_var_assert (dict, "ROWTYPE_", 8);
7301 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7304 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
7308 dict_create_var_assert (dict, "VARNAME_", 8);
7310 const char *dup_var = msave_add_vars (dict, &common->variables);
7313 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
7325 matrix_cmd_execute_msave (struct msave_command *msave)
7327 struct msave_common *common = msave->common;
7328 gsl_matrix *m = NULL;
7329 gsl_vector *factors = NULL;
7330 gsl_vector *splits = NULL;
7332 m = matrix_expr_evaluate (msave->expr);
7336 if (!common->variables.n)
7337 for (size_t i = 0; i < m->size2; i++)
7338 string_array_append_nocopy (&common->variables,
7339 xasprintf ("COL%zu", i + 1));
7340 else if (m->size2 != common->variables.n)
7343 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7344 m->size2, common->variables.n);
7350 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7354 if (!common->fnames.n)
7355 for (size_t i = 0; i < factors->size; i++)
7356 string_array_append_nocopy (&common->fnames,
7357 xasprintf ("FAC%zu", i + 1));
7358 else if (factors->size != common->fnames.n)
7360 msg (SE, _("There are %zu factor variables, "
7361 "but %zu factor values were supplied."),
7362 common->fnames.n, factors->size);
7369 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7373 if (!common->snames.n)
7374 for (size_t i = 0; i < splits->size; i++)
7375 string_array_append_nocopy (&common->snames,
7376 xasprintf ("SPL%zu", i + 1));
7377 else if (splits->size != common->snames.n)
7379 msg (SE, _("There are %zu split variables, "
7380 "but %zu split values were supplied."),
7381 common->snames.n, splits->size);
7386 if (!common->writer)
7388 struct dictionary *dict = msave_create_dict (common);
7392 common->writer = any_writer_open (common->outfile, dict);
7393 if (!common->writer)
7399 common->dict = dict;
7402 bool matrix = (!strcmp (msave->rowtype, "COV")
7403 || !strcmp (msave->rowtype, "CORR"));
7404 for (size_t y = 0; y < m->size1; y++)
7406 struct ccase *c = case_create (dict_get_proto (common->dict));
7409 /* Split variables */
7411 for (size_t i = 0; i < splits->size; i++)
7412 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
7415 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7416 msave->rowtype, ' ');
7420 for (size_t i = 0; i < factors->size; i++)
7421 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
7424 const char *varname_ = (matrix && y < common->variables.n
7425 ? common->variables.strings[y]
7427 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7430 /* Continuous variables. */
7431 for (size_t x = 0; x < m->size2; x++)
7432 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
7433 casewriter_write (common->writer, c);
7437 gsl_matrix_free (m);
7438 gsl_vector_free (factors);
7439 gsl_vector_free (splits);
7442 static struct matrix_cmd *
7443 matrix_parse_mget (struct matrix_state *s)
7445 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7446 *cmd = (struct matrix_cmd) {
7450 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
7454 struct mget_command *mget = &cmd->mget;
7456 lex_match (s->lexer, T_SLASH);
7457 while (lex_token (s->lexer) != T_ENDCMD)
7459 if (lex_match_id (s->lexer, "FILE"))
7461 lex_match (s->lexer, T_EQUALS);
7463 fh_unref (mget->file);
7464 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7468 else if (lex_match_id (s->lexer, "ENCODING"))
7470 lex_match (s->lexer, T_EQUALS);
7471 if (!lex_force_string (s->lexer))
7474 free (mget->encoding);
7475 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
7479 else if (lex_match_id (s->lexer, "TYPE"))
7481 lex_match (s->lexer, T_EQUALS);
7482 while (lex_token (s->lexer) != T_SLASH
7483 && lex_token (s->lexer) != T_ENDCMD)
7485 const char *rowtype = match_rowtype (s->lexer);
7489 stringi_set_insert (&mget->rowtypes, rowtype);
7494 lex_error_expecting (s->lexer, "FILE", "TYPE");
7497 lex_match (s->lexer, T_SLASH);
7502 matrix_cmd_destroy (cmd);
7506 static const struct variable *
7507 get_a8_var (const struct dictionary *d, const char *name)
7509 const struct variable *v = dict_lookup_var (d, name);
7512 msg (SE, _("Matrix data file lacks %s variable."), name);
7515 if (var_get_width (v) != 8)
7517 msg (SE, _("%s variable in matrix data file must be "
7518 "8-byte string, but it has width %d."),
7519 name, var_get_width (v));
7526 var_changed (const struct ccase *ca, const struct ccase *cb,
7527 const struct variable *var)
7530 ? !value_equal (case_data (ca, var), case_data (cb, var),
7531 var_get_width (var))
7536 vars_changed (const struct ccase *ca, const struct ccase *cb,
7537 const struct dictionary *d,
7538 size_t first_var, size_t n_vars)
7540 for (size_t i = 0; i < n_vars; i++)
7542 const struct variable *v = dict_get_var (d, first_var + i);
7543 if (var_changed (ca, cb, v))
7550 vars_all_missing (const struct ccase *c, const struct dictionary *d,
7551 size_t first_var, size_t n_vars)
7553 for (size_t i = 0; i < n_vars; i++)
7554 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
7560 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
7561 const struct dictionary *d,
7562 const struct variable *rowtype_var,
7563 const struct stringi_set *accepted_rowtypes,
7564 struct matrix_state *s,
7565 size_t ss, size_t sn, size_t si,
7566 size_t fs, size_t fn, size_t fi,
7567 size_t cs, size_t cn,
7568 struct pivot_table *pt,
7569 struct pivot_dimension *var_dimension)
7574 /* Is this a matrix for pooled data, either where there are no factor
7575 variables or the factor variables are missing? */
7576 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
7578 struct substring rowtype = case_ss (rows[0], rowtype_var);
7579 ss_rtrim (&rowtype, ss_cstr (" "));
7580 if (!stringi_set_is_empty (accepted_rowtypes)
7581 && !stringi_set_contains_len (accepted_rowtypes,
7582 rowtype.string, rowtype.length))
7585 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
7586 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
7587 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
7588 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
7589 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
7590 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
7594 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
7595 (int) rowtype.length, rowtype.string);
7599 struct string name = DS_EMPTY_INITIALIZER;
7600 ds_put_cstr (&name, prefix);
7602 ds_put_format (&name, "F%zu", fi);
7604 ds_put_format (&name, "S%zu", si);
7606 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
7608 mv = matrix_var_create (s, ds_ss (&name));
7611 msg (SW, _("Matrix data file contains variable with existing name %s."),
7613 goto exit_free_name;
7616 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
7617 size_t n_missing = 0;
7618 for (size_t y = 0; y < n_rows; y++)
7620 for (size_t x = 0; x < cn; x++)
7622 struct variable *var = dict_get_var (d, cs + x);
7623 double value = case_num (rows[y], var);
7624 if (var_is_num_missing (var, value, MV_ANY))
7629 gsl_matrix_set (m, y, x, value);
7633 int var_index = pivot_category_create_leaf (
7634 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
7635 double values[] = { n_rows, cn };
7636 for (size_t j = 0; j < sn; j++)
7638 struct variable *var = dict_get_var (d, ss + j);
7639 const union value *value = case_data (rows[0], var);
7640 pivot_table_put2 (pt, j, var_index,
7641 pivot_value_new_var_value (var, value));
7643 for (size_t j = 0; j < fn; j++)
7645 struct variable *var = dict_get_var (d, fs + j);
7646 const union value sysmis = { .f = SYSMIS };
7647 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
7648 pivot_table_put2 (pt, j + sn, var_index,
7649 pivot_value_new_var_value (var, value));
7651 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
7652 pivot_table_put2 (pt, j + sn + fn, var_index,
7653 pivot_value_new_integer (values[j]));
7656 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
7657 "value, which was treated as zero.",
7658 "Matrix data file variable %s contains %zu missing "
7659 "values, which were treated as zero.", n_missing),
7660 ds_cstr (&name), n_missing);
7667 for (size_t y = 0; y < n_rows; y++)
7668 case_unref (rows[y]);
7672 matrix_cmd_execute_mget__ (struct mget_command *mget,
7673 struct casereader *r, const struct dictionary *d)
7675 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
7676 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
7677 if (!rowtype_ || !varname_)
7680 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
7682 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
7685 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
7687 msg (SE, _("Matrix data file contains no continuous variables."));
7691 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
7693 const struct variable *v = dict_get_var (d, i);
7694 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
7697 _("Matrix data file contains unexpected string variable %s."),
7703 /* SPLIT variables. */
7705 size_t sn = var_get_dict_index (rowtype_);
7706 struct ccase *sc = NULL;
7709 /* FACTOR variables. */
7710 size_t fs = var_get_dict_index (rowtype_) + 1;
7711 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
7712 struct ccase *fc = NULL;
7715 /* Continuous variables. */
7716 size_t cs = var_get_dict_index (varname_) + 1;
7717 size_t cn = dict_get_var_cnt (d) - cs;
7718 struct ccase *cc = NULL;
7721 struct pivot_table *pt = pivot_table_create (
7722 N_("Matrix Variables Created by MGET"));
7723 struct pivot_dimension *attr_dimension = pivot_dimension_create (
7724 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
7725 struct pivot_dimension *var_dimension = pivot_dimension_create (
7726 pt, PIVOT_AXIS_ROW, N_("Variable"));
7729 struct pivot_category *splits = pivot_category_create_group (
7730 attr_dimension->root, N_("Split Values"));
7731 for (size_t i = 0; i < sn; i++)
7732 pivot_category_create_leaf (splits, pivot_value_new_variable (
7733 dict_get_var (d, ss + i)));
7737 struct pivot_category *factors = pivot_category_create_group (
7738 attr_dimension->root, N_("Factors"));
7739 for (size_t i = 0; i < fn; i++)
7740 pivot_category_create_leaf (factors, pivot_value_new_variable (
7741 dict_get_var (d, fs + i)));
7743 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
7744 N_("Rows"), N_("Columns"));
7747 struct ccase **rows = NULL;
7748 size_t allocated_rows = 0;
7752 while ((c = casereader_read (r)) != NULL)
7754 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
7764 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
7765 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
7766 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
7769 if (change != NOTHING_CHANGED)
7771 matrix_mget_commit_var (rows, n_rows, d,
7772 rowtype_, &mget->rowtypes,
7783 if (n_rows >= allocated_rows)
7784 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
7787 if (change == SPLITS_CHANGED)
7793 /* Reset the factor number, if there are factors. */
7797 if (row_has_factors)
7803 else if (change == FACTORS_CHANGED)
7805 if (row_has_factors)
7811 matrix_mget_commit_var (rows, n_rows, d,
7812 rowtype_, &mget->rowtypes,
7824 if (var_dimension->n_leaves)
7825 pivot_table_submit (pt);
7827 pivot_table_unref (pt);
7831 matrix_cmd_execute_mget (struct mget_command *mget)
7833 struct casereader *r;
7834 struct dictionary *d;
7835 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
7836 mget->state->dataset, &r, &d))
7838 matrix_cmd_execute_mget__ (mget, r, d);
7839 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
7844 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
7846 if (!lex_force_id (s->lexer))
7849 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
7851 *varp = matrix_var_create (s, lex_tokss (s->lexer));
7856 static struct matrix_cmd *
7857 matrix_parse_eigen (struct matrix_state *s)
7859 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7860 *cmd = (struct matrix_cmd) {
7862 .eigen = { .expr = NULL }
7865 struct eigen_command *eigen = &cmd->eigen;
7866 if (!lex_force_match (s->lexer, T_LPAREN))
7868 eigen->expr = matrix_parse_expr (s);
7870 || !lex_force_match (s->lexer, T_COMMA)
7871 || !matrix_parse_dst_var (s, &eigen->evec)
7872 || !lex_force_match (s->lexer, T_COMMA)
7873 || !matrix_parse_dst_var (s, &eigen->eval)
7874 || !lex_force_match (s->lexer, T_RPAREN))
7880 matrix_cmd_destroy (cmd);
7885 matrix_cmd_execute_eigen (struct eigen_command *eigen)
7887 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
7890 if (!is_symmetric (A))
7892 msg (SE, _("Argument of EIGEN must be symmetric."));
7893 gsl_matrix_free (A);
7897 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
7898 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
7899 gsl_vector v_eval = to_vector (eval);
7900 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
7901 gsl_eigen_symmv (A, &v_eval, evec, w);
7902 gsl_eigen_symmv_free (w);
7904 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
7906 gsl_matrix_free (A);
7908 gsl_matrix_free (eigen->eval->value);
7909 eigen->eval->value = eval;
7911 gsl_matrix_free (eigen->evec->value);
7912 eigen->evec->value = evec;
7915 static struct matrix_cmd *
7916 matrix_parse_setdiag (struct matrix_state *s)
7918 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7919 *cmd = (struct matrix_cmd) {
7920 .type = MCMD_SETDIAG,
7921 .setdiag = { .dst = NULL }
7924 struct setdiag_command *setdiag = &cmd->setdiag;
7925 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
7928 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
7931 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
7936 if (!lex_force_match (s->lexer, T_COMMA))
7939 setdiag->expr = matrix_parse_expr (s);
7943 if (!lex_force_match (s->lexer, T_RPAREN))
7949 matrix_cmd_destroy (cmd);
7954 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
7956 gsl_matrix *dst = setdiag->dst->value;
7959 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
7960 setdiag->dst->name);
7964 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
7968 size_t n = MIN (dst->size1, dst->size2);
7969 if (is_scalar (src))
7971 double d = to_scalar (src);
7972 for (size_t i = 0; i < n; i++)
7973 gsl_matrix_set (dst, i, i, d);
7975 else if (is_vector (src))
7977 gsl_vector v = to_vector (src);
7978 for (size_t i = 0; i < n && i < v.size; i++)
7979 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
7982 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
7983 "not a %zu×%zu matrix."),
7984 src->size1, src->size2);
7985 gsl_matrix_free (src);
7988 static struct matrix_cmd *
7989 matrix_parse_svd (struct matrix_state *s)
7991 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7992 *cmd = (struct matrix_cmd) {
7994 .svd = { .expr = NULL }
7997 struct svd_command *svd = &cmd->svd;
7998 if (!lex_force_match (s->lexer, T_LPAREN))
8000 svd->expr = matrix_parse_expr (s);
8002 || !lex_force_match (s->lexer, T_COMMA)
8003 || !matrix_parse_dst_var (s, &svd->u)
8004 || !lex_force_match (s->lexer, T_COMMA)
8005 || !matrix_parse_dst_var (s, &svd->s)
8006 || !lex_force_match (s->lexer, T_COMMA)
8007 || !matrix_parse_dst_var (s, &svd->v)
8008 || !lex_force_match (s->lexer, T_RPAREN))
8014 matrix_cmd_destroy (cmd);
8019 matrix_cmd_execute_svd (struct svd_command *svd)
8021 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8025 if (m->size1 >= m->size2)
8028 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8029 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8030 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8031 gsl_vector *work = gsl_vector_alloc (A->size2);
8032 gsl_linalg_SV_decomp (A, V, &Sv, work);
8033 gsl_vector_free (work);
8035 matrix_var_set (svd->u, A);
8036 matrix_var_set (svd->s, S);
8037 matrix_var_set (svd->v, V);
8041 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8042 gsl_matrix_transpose_memcpy (At, m);
8043 gsl_matrix_free (m);
8045 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8046 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8047 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8048 gsl_vector *work = gsl_vector_alloc (At->size2);
8049 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8050 gsl_vector_free (work);
8052 matrix_var_set (svd->v, At);
8053 matrix_var_set (svd->s, St);
8054 matrix_var_set (svd->u, Vt);
8059 matrix_cmd_execute (struct matrix_cmd *cmd)
8064 matrix_cmd_execute_compute (&cmd->compute);
8068 matrix_cmd_execute_print (&cmd->print);
8072 return matrix_cmd_execute_do_if (&cmd->do_if);
8075 matrix_cmd_execute_loop (&cmd->loop);
8082 matrix_cmd_execute_display (&cmd->display);
8086 matrix_cmd_execute_release (&cmd->release);
8090 matrix_cmd_execute_save (&cmd->save);
8094 matrix_cmd_execute_read (&cmd->read);
8098 matrix_cmd_execute_write (&cmd->write);
8102 matrix_cmd_execute_get (&cmd->get);
8106 matrix_cmd_execute_msave (&cmd->msave);
8110 matrix_cmd_execute_mget (&cmd->mget);
8114 matrix_cmd_execute_eigen (&cmd->eigen);
8118 matrix_cmd_execute_setdiag (&cmd->setdiag);
8122 matrix_cmd_execute_svd (&cmd->svd);
8130 matrix_cmds_uninit (struct matrix_cmds *cmds)
8132 for (size_t i = 0; i < cmds->n; i++)
8133 matrix_cmd_destroy (cmds->commands[i]);
8134 free (cmds->commands);
8138 matrix_cmd_destroy (struct matrix_cmd *cmd)
8146 matrix_lvalue_destroy (cmd->compute.lvalue);
8147 matrix_expr_destroy (cmd->compute.rvalue);
8151 matrix_expr_destroy (cmd->print.expression);
8152 free (cmd->print.title);
8153 print_labels_destroy (cmd->print.rlabels);
8154 print_labels_destroy (cmd->print.clabels);
8158 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8160 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8161 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
8163 free (cmd->do_if.clauses);
8167 matrix_expr_destroy (cmd->loop.start);
8168 matrix_expr_destroy (cmd->loop.end);
8169 matrix_expr_destroy (cmd->loop.increment);
8170 matrix_expr_destroy (cmd->loop.top_condition);
8171 matrix_expr_destroy (cmd->loop.bottom_condition);
8172 matrix_cmds_uninit (&cmd->loop.commands);
8182 free (cmd->release.vars);
8186 matrix_expr_destroy (cmd->save.expression);
8190 matrix_lvalue_destroy (cmd->read.dst);
8191 matrix_expr_destroy (cmd->read.size);
8195 matrix_expr_destroy (cmd->write.expression);
8196 free (cmd->write.format);
8200 matrix_lvalue_destroy (cmd->get.dst);
8201 fh_unref (cmd->get.file);
8202 free (cmd->get.encoding);
8203 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8207 matrix_expr_destroy (cmd->msave.expr);
8211 fh_unref (cmd->mget.file);
8212 stringi_set_destroy (&cmd->mget.rowtypes);
8216 matrix_expr_destroy (cmd->eigen.expr);
8220 matrix_expr_destroy (cmd->setdiag.expr);
8224 matrix_expr_destroy (cmd->svd.expr);
8230 struct matrix_command_name
8233 struct matrix_cmd *(*parse) (struct matrix_state *);
8236 static const struct matrix_command_name *
8237 matrix_parse_command_name (struct lexer *lexer)
8239 static const struct matrix_command_name commands[] = {
8240 { "COMPUTE", matrix_parse_compute },
8241 { "CALL EIGEN", matrix_parse_eigen },
8242 { "CALL SETDIAG", matrix_parse_setdiag },
8243 { "CALL SVD", matrix_parse_svd },
8244 { "PRINT", matrix_parse_print },
8245 { "DO IF", matrix_parse_do_if },
8246 { "LOOP", matrix_parse_loop },
8247 { "BREAK", matrix_parse_break },
8248 { "READ", matrix_parse_read },
8249 { "WRITE", matrix_parse_write },
8250 { "GET", matrix_parse_get },
8251 { "SAVE", matrix_parse_save },
8252 { "MGET", matrix_parse_mget },
8253 { "MSAVE", matrix_parse_msave },
8254 { "DISPLAY", matrix_parse_display },
8255 { "RELEASE", matrix_parse_release },
8257 static size_t n = sizeof commands / sizeof *commands;
8259 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8260 if (lex_match_phrase (lexer, c->name))
8265 static struct matrix_cmd *
8266 matrix_parse_command (struct matrix_state *s)
8268 size_t nesting_level = SIZE_MAX;
8270 struct matrix_cmd *c = NULL;
8271 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
8273 lex_error (s->lexer, _("Unknown matrix command."));
8274 else if (!cmd->parse)
8275 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8279 nesting_level = output_open_group (
8280 group_item_create_nocopy (utf8_to_title (cmd->name),
8281 utf8_to_title (cmd->name)));
8286 lex_end_of_command (s->lexer);
8287 lex_discard_rest_of_command (s->lexer);
8288 if (nesting_level != SIZE_MAX)
8289 output_close_groups (nesting_level);
8295 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8297 if (!lex_force_match (lexer, T_ENDCMD))
8300 struct matrix_state state = {
8302 .session = dataset_session (ds),
8304 .vars = HMAP_INITIALIZER (state.vars),
8309 while (lex_match (lexer, T_ENDCMD))
8311 if (lex_token (lexer) == T_STOP)
8313 msg (SE, _("Unexpected end of input expecting matrix command."));
8317 if (lex_match_phrase (lexer, "END MATRIX"))
8320 struct matrix_cmd *c = matrix_parse_command (&state);
8323 matrix_cmd_execute (c);
8324 matrix_cmd_destroy (c);
8328 struct matrix_var *var, *next;
8329 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8332 gsl_matrix_free (var->value);
8333 hmap_delete (&state.vars, &var->hmap_node);
8336 hmap_destroy (&state.vars);
8337 msave_common_destroy (state.common);
8338 fh_unref (state.prev_read_file);
8339 for (size_t i = 0; i < state.n_read_files; i++)
8340 read_file_destroy (state.read_files[i]);
8341 free (state.read_files);
8342 fh_unref (state.prev_write_file);
8343 for (size_t i = 0; i < state.n_write_files; i++)
8344 write_file_destroy (state.write_files[i]);
8345 free (state.write_files);
8346 fh_unref (state.prev_save_file);
8347 for (size_t i = 0; i < state.n_save_files; i++)
8348 save_file_destroy (state.save_files[i]);
8349 free (state.save_files);