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 #define MATRIX_FUNCTIONS \
184 F(ABS, "ABS", m_e, NULL) \
185 F(ALL, "ALL", d_m, NULL) \
186 F(ANY, "ANY", d_m, NULL) \
187 F(ARSIN, "ARSIN", m_e, "a[-1,1]") \
188 F(ARTAN, "ARTAN", m_e, NULL) \
189 F(BLOCK, "BLOCK", m_any, NULL) \
190 F(CHOL, "CHOL", m_m, NULL) \
191 F(CMIN, "CMIN", m_m, NULL) \
192 F(CMAX, "CMAX", m_m, NULL) \
193 F(COS, "COS", m_e, NULL) \
194 F(CSSQ, "CSSQ", m_m, NULL) \
195 F(CSUM, "CSUM", m_m, NULL) \
196 F(DESIGN, "DESIGN", m_m, NULL) \
197 F(DET, "DET", d_m, NULL) \
198 F(DIAG, "DIAG", m_m, NULL) \
199 F(EVAL, "EVAL", m_m, NULL) \
200 F(EXP, "EXP", m_e, NULL) \
201 F(GINV, "GINV", m_m, NULL) \
202 F(GRADE, "GRADE", m_m, NULL) \
203 F(GSCH, "GSCH", m_m, NULL) \
204 F(IDENT, "IDENT", IDENT, NULL) \
205 F(INV, "INV", m_m, NULL) \
206 F(KRONEKER, "KRONEKER", m_mm, NULL) \
207 F(LG10, "LG10", m_e, "a>0") \
208 F(LN, "LN", m_e, "a>0") \
209 F(MAGIC, "MAGIC", m_d, NULL) \
210 F(MAKE, "MAKE", m_ddd, NULL) \
211 F(MDIAG, "MDIAG", m_v, NULL) \
212 F(MMAX, "MMAX", d_m, NULL) \
213 F(MMIN, "MMIN", d_m, NULL) \
214 F(MOD, "MOD", m_md, NULL) \
215 F(MSSQ, "MSSQ", d_m, NULL) \
216 F(MSUM, "MSUM", d_m, NULL) \
217 F(NCOL, "NCOL", d_m, NULL) \
218 F(NROW, "NROW", d_m, NULL) \
219 F(RANK, "RANK", d_m, NULL) \
220 F(RESHAPE, "RESHAPE", m_mdd, NULL) \
221 F(RMAX, "RMAX", m_m, NULL) \
222 F(RMIN, "RMIN", m_m, NULL) \
223 F(RND, "RND", m_e, NULL) \
224 F(RNKORDER, "RNKORDER", m_m, NULL) \
225 F(RSSQ, "RSSQ", m_m, NULL) \
226 F(RSUM, "RSUM", m_m, NULL) \
227 F(SIN, "SIN", m_e, NULL) \
228 F(SOLVE, "SOLVE", m_mm, NULL) \
229 F(SQRT, "SQRT", m_m, NULL) \
230 F(SSCP, "SSCP", m_m, NULL) \
231 F(SVAL, "SVAL", m_m, NULL) \
232 F(SWEEP, "SWEEP", m_md, NULL) \
233 F(T, "T", m_m, NULL) \
234 F(TRACE, "TRACE", d_m, NULL) \
235 F(TRANSPOS, "TRANSPOS", m_m, NULL) \
236 F(TRUNC, "TRUNC", m_e, NULL) \
237 F(UNIFORM, "UNIFORM", m_dd, NULL) \
238 F(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0") \
239 F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0") \
240 F(IDF_BETA, "IDF.BETA", m_edd, "a[0,1] b>0 c>0") \
241 F(RV_BETA, "RV.BETA", d_dd, "a>0 b>0") \
242 F(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
243 F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
244 F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]") \
245 F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]") \
246 F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0") \
247 F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0") \
248 F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "c>0") \
249 F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \
250 F(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0") \
251 F(CHICDF, "CHICDF", m_ed, "a>=0 b>0") \
252 F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0") \
253 F(PDF_CHISQ, "PDF.CHISQ", m_ed, "a>=0 b>0") \
254 F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \
255 F(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0") \
256 F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0") \
257 F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0") \
258 F(PDF_EXP, "PDF.EXP", m_ed, "a>=0 b>0") \
259 F(RV_EXP, "RV.EXP", d_d, "a>0") \
260 F(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0") \
261 F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \
262 F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0") \
263 F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0") \
264 F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0") \
265 F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0") \
266 F(RV_F, "RV.F", d_dd, "a>0 b>0") \
267 F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0") \
268 F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
269 F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0") \
270 F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
271 F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0") \
272 F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL) \
273 F(RV_LANDAU, "RV.LANDAU", d_none, NULL) \
274 F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0") \
275 F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0") \
276 F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0") \
277 F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0") \
278 F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]") \
279 F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]") \
280 F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0") \
281 F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0") \
282 F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0") \
283 F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0") \
284 F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
285 F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0") \
286 F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
287 F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0") \
288 F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0") \
289 F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0") \
290 F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0") \
291 F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0") \
292 F(CDFNORM, "CDFNORM", m_e, NULL) \
293 F(PROBIT, "PROBIT", m_e, "a(0,1)") \
294 F(NORMAL, "NORMAL", m_e, "a>0") \
295 F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0") \
296 F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0") \
297 F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0") \
298 F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0") \
299 F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0") \
300 F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0") \
301 F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0") \
302 F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0") \
303 F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0") \
304 F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0") \
305 F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL) \
306 F(RV_RTAIL, "RV.RTAIL", d_dd, NULL) \
307 F(CDF_T, "CDF.T", m_ed, "b>0") \
308 F(TCDF, "TCDF", m_ed, "b>0") \
309 F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0") \
310 F(PDF_T, "PDF.T", m_ed, "b>0") \
311 F(RV_T, "RV.T", d_d, "a>0") \
312 F(CDF_T1G, "CDF.T1G", m_edd, NULL) \
313 F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)") \
314 F(PDF_T1G, "PDF.T1G", m_edd, NULL) \
315 F(RV_T1G, "RV.T1G", d_dd, NULL) \
316 F(CDF_T2G, "CDF.T2G", m_edd, NULL) \
317 F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)") \
318 F(PDF_T2G, "PDF.T2G", m_edd, NULL) \
319 F(RV_T2G, "RV.T2G", d_dd, NULL) \
320 F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c") \
321 F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c") \
322 F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c") \
323 F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b") \
324 F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
325 F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0") \
326 F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
327 F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0") \
328 F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
329 F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
330 F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]") \
331 F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]") \
332 F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]") \
333 F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]") \
334 F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]") \
335 F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]") \
336 F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]") \
337 F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
338 F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
339 F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a") \
340 F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]") \
341 F(RV_LOG, "RV.LOG", d_d, "a(0,1]") \
342 F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
343 F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
344 F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]") \
345 F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0") \
346 F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0") \
347 F(RV_POISSON, "RV.POISSON", d_d, "a>0")
349 struct matrix_function_properties
352 const char *constraints;
355 enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 };
356 enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
357 enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
358 enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 };
359 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
360 enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 };
361 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
362 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
363 enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 };
364 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
365 enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 };
366 enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 };
367 enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
368 enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
369 enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
370 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
371 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
372 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
373 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
374 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
376 typedef double matrix_proto_d_none (void);
377 typedef double matrix_proto_d_d (double);
378 typedef double matrix_proto_d_dd (double, double);
379 typedef double matrix_proto_d_dd (double, double);
380 typedef double matrix_proto_d_ddd (double, double, double);
381 typedef gsl_matrix *matrix_proto_m_d (double);
382 typedef gsl_matrix *matrix_proto_m_dd (double, double);
383 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
384 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
385 typedef double matrix_proto_m_e (double);
386 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
387 typedef double matrix_proto_m_ed (double, double);
388 typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double);
389 typedef double matrix_proto_m_edd (double, double, double);
390 typedef double matrix_proto_m_eddd (double, double, double, double);
391 typedef double matrix_proto_m_eed (double, double, double);
392 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
393 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
394 typedef double matrix_proto_d_m (gsl_matrix *);
395 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
396 typedef gsl_matrix *matrix_proto_IDENT (double, double);
398 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
399 static matrix_proto_##PROTO matrix_eval_##ENUM;
403 /* Matrix expressions. */
410 #define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
414 /* Elementwise and scalar arithmetic. */
415 MOP_NEGATE, /* unary - */
416 MOP_ADD_ELEMS, /* + */
417 MOP_SUB_ELEMS, /* - */
418 MOP_MUL_ELEMS, /* &* */
419 MOP_DIV_ELEMS, /* / and &/ */
420 MOP_EXP_ELEMS, /* &** */
422 MOP_SEQ_BY, /* a:b:c */
424 /* Matrix arithmetic. */
426 MOP_EXP_MAT, /* ** */
443 MOP_PASTE_HORZ, /* a, b, c, ... */
444 MOP_PASTE_VERT, /* a; b; c; ... */
448 MOP_VEC_INDEX, /* x(y) */
449 MOP_VEC_ALL, /* x(:) */
450 MOP_MAT_INDEX, /* x(y,z) */
451 MOP_ROW_INDEX, /* x(y,:) */
452 MOP_COL_INDEX, /* x(:,z) */
459 MOP_EOF, /* EOF('file') */
467 struct matrix_expr **subs;
472 struct matrix_var *variable;
473 struct read_file *eof;
478 matrix_expr_destroy (struct matrix_expr *e)
485 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
516 for (size_t i = 0; i < e->n_subs; i++)
517 matrix_expr_destroy (e->subs[i]);
529 static struct matrix_expr *
530 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
533 struct matrix_expr *e = xmalloc (sizeof *e);
534 *e = (struct matrix_expr) {
536 .subs = xmemdup (subs, n_subs * sizeof *subs),
542 static struct matrix_expr *
543 matrix_expr_create_0 (enum matrix_op op)
545 struct matrix_expr *sub;
546 return matrix_expr_create_subs (op, &sub, 0);
549 static struct matrix_expr *
550 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
552 return matrix_expr_create_subs (op, &sub, 1);
555 static struct matrix_expr *
556 matrix_expr_create_2 (enum matrix_op op,
557 struct matrix_expr *sub0, struct matrix_expr *sub1)
559 struct matrix_expr *subs[] = { sub0, sub1 };
560 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
563 static struct matrix_expr *
564 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
565 struct matrix_expr *sub1, struct matrix_expr *sub2)
567 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
568 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
571 static struct matrix_expr *
572 matrix_expr_create_number (double number)
574 struct matrix_expr *e = xmalloc (sizeof *e);
575 *e = (struct matrix_expr) { .op = MOP_NUMBER, .number = number };
579 static struct matrix_expr *matrix_parse_or_xor (struct matrix_state *);
581 static struct matrix_expr *
582 matrix_parse_curly_comma (struct matrix_state *s)
584 struct matrix_expr *lhs = matrix_parse_or_xor (s);
588 while (lex_match (s->lexer, T_COMMA))
590 struct matrix_expr *rhs = matrix_parse_or_xor (s);
593 matrix_expr_destroy (lhs);
596 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
601 static struct matrix_expr *
602 matrix_parse_curly_semi (struct matrix_state *s)
604 if (lex_token (s->lexer) == T_RCURLY)
605 return matrix_expr_create_0 (MOP_EMPTY);
607 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
611 while (lex_match (s->lexer, T_SEMICOLON))
613 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
616 matrix_expr_destroy (lhs);
619 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
624 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
625 for (size_t Y = 0; Y < (M)->size1; Y++) \
626 for (size_t X = 0; X < (M)->size2; X++) \
627 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
630 is_vector (const gsl_matrix *m)
632 return m->size1 <= 1 || m->size2 <= 1;
636 to_vector (gsl_matrix *m)
638 return (m->size1 == 1
639 ? gsl_matrix_row (m, 0).vector
640 : gsl_matrix_column (m, 0).vector);
645 matrix_eval_ABS (double d)
651 matrix_eval_ALL (gsl_matrix *m)
653 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
660 matrix_eval_ANY (gsl_matrix *m)
662 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
669 matrix_eval_ARSIN (double d)
675 matrix_eval_ARTAN (double d)
681 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
685 for (size_t i = 0; i < n; i++)
690 gsl_matrix *block = gsl_matrix_calloc (r, c);
692 for (size_t i = 0; i < n; i++)
694 for (size_t y = 0; y < m[i]->size1; y++)
695 for (size_t x = 0; x < m[i]->size2; x++)
696 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
704 matrix_eval_CHOL (gsl_matrix *m)
706 if (!gsl_linalg_cholesky_decomp1 (m))
708 for (size_t y = 0; y < m->size1; y++)
709 for (size_t x = y + 1; x < m->size2; x++)
710 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
712 for (size_t y = 0; y < m->size1; y++)
713 for (size_t x = 0; x < y; x++)
714 gsl_matrix_set (m, y, x, 0);
719 msg (SE, _("Input to CHOL function is not positive-definite."));
725 matrix_eval_col_extremum (gsl_matrix *m, bool min)
730 return gsl_matrix_alloc (1, 0);
732 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
733 for (size_t x = 0; x < m->size2; x++)
735 double ext = gsl_matrix_get (m, 0, x);
736 for (size_t y = 1; y < m->size1; y++)
738 double value = gsl_matrix_get (m, y, x);
739 if (min ? value < ext : value > ext)
742 gsl_matrix_set (cext, 0, x, ext);
748 matrix_eval_CMAX (gsl_matrix *m)
750 return matrix_eval_col_extremum (m, false);
754 matrix_eval_CMIN (gsl_matrix *m)
756 return matrix_eval_col_extremum (m, true);
760 matrix_eval_COS (double d)
766 matrix_eval_col_sum (gsl_matrix *m, bool square)
771 return gsl_matrix_alloc (1, 0);
773 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
774 for (size_t x = 0; x < m->size2; x++)
777 for (size_t y = 0; y < m->size1; y++)
779 double d = gsl_matrix_get (m, y, x);
780 sum += square ? pow2 (d) : d;
782 gsl_matrix_set (result, 0, x, sum);
788 matrix_eval_CSSQ (gsl_matrix *m)
790 return matrix_eval_col_sum (m, true);
794 matrix_eval_CSUM (gsl_matrix *m)
796 return matrix_eval_col_sum (m, false);
800 compare_double_3way (const void *a_, const void *b_)
802 const double *a = a_;
803 const double *b = b_;
804 return *a < *b ? -1 : *a > *b;
808 matrix_eval_DESIGN (gsl_matrix *m)
810 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
811 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
812 gsl_matrix_transpose_memcpy (&m2, m);
814 for (size_t y = 0; y < m2.size1; y++)
815 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
817 size_t *n = xcalloc (m2.size1, sizeof *n);
819 for (size_t i = 0; i < m2.size1; i++)
821 double *row = tmp + m2.size2 * i;
822 for (size_t j = 0; j < m2.size2; )
825 for (k = j + 1; k < m2.size2; k++)
826 if (row[j] != row[k])
828 row[n[i]++] = row[j];
833 msg (MW, _("Column %zu in DESIGN argument has constant value."), i + 1);
838 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
840 for (size_t i = 0; i < m->size2; i++)
845 const double *unique = tmp + m2.size2 * i;
846 for (size_t j = 0; j < n[i]; j++, x++)
848 double value = unique[j];
849 for (size_t y = 0; y < m->size1; y++)
850 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
861 matrix_eval_DET (gsl_matrix *m)
863 gsl_permutation *p = gsl_permutation_alloc (m->size1);
865 gsl_linalg_LU_decomp (m, p, &signum);
866 gsl_permutation_free (p);
867 return gsl_linalg_LU_det (m, signum);
871 matrix_eval_DIAG (gsl_matrix *m)
873 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
874 for (size_t i = 0; i < diag->size1; i++)
875 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
880 is_symmetric (const gsl_matrix *m)
882 if (m->size1 != m->size2)
885 for (size_t y = 0; y < m->size1; y++)
886 for (size_t x = 0; x < y; x++)
887 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
894 compare_double_desc (const void *a_, const void *b_)
896 const double *a = a_;
897 const double *b = b_;
898 return *a > *b ? -1 : *a < *b;
902 matrix_eval_EVAL (gsl_matrix *m)
904 if (!is_symmetric (m))
906 msg (SE, _("Argument of EVAL must be symmetric."));
910 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
911 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
912 gsl_vector v_eval = to_vector (eval);
913 gsl_eigen_symm (m, &v_eval, w);
914 gsl_eigen_symm_free (w);
916 assert (v_eval.stride == 1);
917 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
923 matrix_eval_EXP (double d)
928 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
931 Charl Linssen <charl@itfromb.it>
935 matrix_eval_GINV (gsl_matrix *A)
940 gsl_matrix *tmp_mat = NULL;
943 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
944 tmp_mat = gsl_matrix_alloc (m, n);
945 gsl_matrix_transpose_memcpy (tmp_mat, A);
953 gsl_matrix *V = gsl_matrix_alloc (m, m);
954 gsl_vector *u = gsl_vector_alloc (m);
956 gsl_vector *tmp_vec = gsl_vector_alloc (m);
957 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
958 gsl_vector_free (tmp_vec);
961 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
962 gsl_matrix_set_zero (Sigma_pinv);
963 double cutoff = 1e-15 * gsl_vector_max (u);
965 for (size_t i = 0; i < m; ++i)
967 double x = gsl_vector_get (u, i);
968 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
971 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
972 gsl_matrix *U = gsl_matrix_calloc (n, n);
973 for (size_t i = 0; i < n; i++)
974 for (size_t j = 0; j < m; j++)
975 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
977 /* two dot products to obtain pseudoinverse */
978 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
979 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
984 A_pinv = gsl_matrix_alloc (n, m);
985 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
989 A_pinv = gsl_matrix_alloc (m, n);
990 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
993 gsl_matrix_free (tmp_mat);
994 gsl_matrix_free (tmp_mat2);
996 gsl_matrix_free (Sigma_pinv);
1010 grade_compare_3way (const void *a_, const void *b_)
1012 const struct grade *a = a_;
1013 const struct grade *b = b_;
1015 return (a->value < b->value ? -1
1016 : a->value > b->value ? 1
1024 matrix_eval_GRADE (gsl_matrix *m)
1026 size_t n = m->size1 * m->size2;
1027 struct grade *grades = xmalloc (n * sizeof *grades);
1030 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1031 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
1032 qsort (grades, n, sizeof *grades, grade_compare_3way);
1034 for (size_t i = 0; i < n; i++)
1035 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
1043 dot (gsl_vector *a, gsl_vector *b)
1045 double result = 0.0;
1046 for (size_t i = 0; i < a->size; i++)
1047 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
1052 norm2 (gsl_vector *v)
1054 double result = 0.0;
1055 for (size_t i = 0; i < v->size; i++)
1056 result += pow2 (gsl_vector_get (v, i));
1061 norm (gsl_vector *v)
1063 return sqrt (norm2 (v));
1067 matrix_eval_GSCH (gsl_matrix *v)
1069 if (v->size2 < v->size1)
1071 msg (SE, _("GSCH requires its argument to have at least as many columns "
1072 "as rows, but it has dimensions %zu×%zu."),
1073 v->size1, v->size2);
1076 if (!v->size1 || !v->size2)
1079 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
1081 for (size_t vx = 0; vx < v->size2; vx++)
1083 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
1084 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
1086 gsl_vector_memcpy (&u_i, &v_i);
1087 for (size_t j = 0; j < ux; j++)
1089 gsl_vector u_j = gsl_matrix_column (u, j).vector;
1090 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
1091 for (size_t k = 0; k < u_i.size; k++)
1092 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
1093 - scale * gsl_vector_get (&u_j, k)));
1096 double len = norm (&u_i);
1099 gsl_vector_scale (&u_i, 1.0 / len);
1100 if (++ux >= v->size1)
1107 msg (SE, _("%zu×%zu argument to GSCH contains only "
1108 "%zu linearly independent columns."),
1109 v->size1, v->size2, ux);
1110 gsl_matrix_free (u);
1114 u->size2 = v->size1;
1119 matrix_eval_IDENT (double s1, double s2)
1121 if (s1 < 0 || s1 > SIZE_MAX || s2 < 0 || s2 > SIZE_MAX)
1123 msg (SE, _("Arguments to IDENT must be non-negative integers."));
1127 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
1128 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1134 invert_matrix (gsl_matrix *x)
1136 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1138 gsl_linalg_LU_decomp (x, p, &signum);
1139 gsl_linalg_LU_invx (x, p);
1140 gsl_permutation_free (p);
1144 matrix_eval_INV (gsl_matrix *m)
1151 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1153 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1154 a->size2 * b->size2);
1156 for (size_t ar = 0; ar < a->size1; ar++)
1157 for (size_t br = 0; br < b->size1; br++, y++)
1160 for (size_t ac = 0; ac < a->size2; ac++)
1161 for (size_t bc = 0; bc < b->size2; bc++, x++)
1163 double av = gsl_matrix_get (a, ar, ac);
1164 double bv = gsl_matrix_get (b, br, bc);
1165 gsl_matrix_set (k, y, x, av * bv);
1172 matrix_eval_LG10 (double d)
1178 matrix_eval_LN (double d)
1184 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1186 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1189 for (size_t i = 1; i <= n * n; i++)
1191 gsl_matrix_set (m, y, x, i);
1193 size_t y1 = !y ? n - 1 : y - 1;
1194 size_t x1 = x + 1 >= n ? 0 : x + 1;
1195 if (gsl_matrix_get (m, y1, x1) == 0)
1201 y = y + 1 >= n ? 0 : y + 1;
1206 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1208 double a = gsl_matrix_get (m, y1, x1);
1209 double b = gsl_matrix_get (m, y2, x2);
1210 gsl_matrix_set (m, y1, x1, b);
1211 gsl_matrix_set (m, y2, x2, a);
1215 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1219 /* A. Umar, "On the Construction of Even Order Magic Squares",
1220 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1222 for (size_t i = 1; i <= n * n / 2; i++)
1224 gsl_matrix_set (m, y, x, i);
1234 for (size_t i = n * n; i > n * n / 2; i--)
1236 gsl_matrix_set (m, y, x, i);
1244 for (size_t y = 0; y < n; y++)
1245 for (size_t x = 0; x < n / 2; x++)
1247 unsigned int d = gsl_matrix_get (m, y, x);
1248 if (d % 2 != (y < n / 2))
1249 magic_exchange (m, y, x, y, n - x - 1);
1254 size_t x1 = n / 2 - 1;
1256 magic_exchange (m, y1, x1, y2, x1);
1257 magic_exchange (m, y1, x2, y2, x2);
1261 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1263 /* A. Umar, "On the Construction of Even Order Magic Squares",
1264 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1268 for (size_t i = 1; ; i++)
1270 gsl_matrix_set (m, y, x, i);
1271 if (++y == n / 2 - 1)
1283 for (size_t i = n * n; ; i--)
1285 gsl_matrix_set (m, y, x, i);
1286 if (++y == n / 2 - 1)
1295 for (size_t y = 0; y < n; y++)
1296 if (y != n / 2 - 1 && y != n / 2)
1297 for (size_t x = 0; x < n / 2; x++)
1299 unsigned int d = gsl_matrix_get (m, y, x);
1300 if (d % 2 != (y < n / 2))
1301 magic_exchange (m, y, x, y, n - x - 1);
1304 size_t a0 = (n * n - 2 * n) / 2 + 1;
1305 for (size_t i = 0; i < n / 2; i++)
1308 gsl_matrix_set (m, n / 2, i, a);
1309 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1311 for (size_t i = 0; i < n / 2; i++)
1313 size_t a = a0 + i + n / 2;
1314 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1315 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1317 for (size_t x = 1; x < n / 2; x += 2)
1318 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1319 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1320 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1321 size_t x1 = n / 2 - 2;
1322 size_t x2 = n / 2 + 1;
1323 size_t y1 = n / 2 - 2;
1324 size_t y2 = n / 2 + 1;
1325 magic_exchange (m, y1, x1, y2, x1);
1326 magic_exchange (m, y1, x2, y2, x2);
1330 matrix_eval_MAGIC (double n_)
1332 if (n_ < 3 || n_ >= sqrt (SIZE_MAX))
1334 msg (SE, _("MAGIC argument must be an integer 3 or greater."));
1339 gsl_matrix *m = gsl_matrix_calloc (n, n);
1341 matrix_eval_MAGIC_odd (m, n);
1343 matrix_eval_MAGIC_singly_even (m, n);
1345 matrix_eval_MAGIC_doubly_even (m, n);
1350 matrix_eval_MAKE (double r, double c, double value)
1352 if (r < 0 || r >= SIZE_MAX || c < 0 || c >= SIZE_MAX)
1354 msg (SE, _("Size arguments to MAKE must be integers."));
1358 gsl_matrix *m = gsl_matrix_alloc (r, c);
1359 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1365 matrix_eval_MDIAG (gsl_vector *v)
1367 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1368 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1369 gsl_vector_memcpy (&diagonal, v);
1374 matrix_eval_MMAX (gsl_matrix *m)
1376 return gsl_matrix_max (m);
1380 matrix_eval_MMIN (gsl_matrix *m)
1382 return gsl_matrix_min (m);
1386 matrix_eval_MOD (gsl_matrix *m, double divisor)
1390 msg (SE, _("Divisor argument to MOD function must be nonzero."));
1394 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1395 *d = fmod (*d, divisor);
1400 matrix_eval_MSSQ (gsl_matrix *m)
1403 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1409 matrix_eval_MSUM (gsl_matrix *m)
1412 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1418 matrix_eval_NCOL (gsl_matrix *m)
1424 matrix_eval_NROW (gsl_matrix *m)
1430 matrix_eval_RANK (gsl_matrix *m)
1432 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1433 gsl_linalg_QR_decomp (m, tau);
1434 gsl_vector_free (tau);
1436 return gsl_linalg_QRPT_rank (m, -1);
1440 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_)
1442 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1444 msg (SE, _("Arguments 2 and 3 to RESHAPE must be integers."));
1449 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1451 msg (SE, _("Product of RESHAPE arguments 2 and 3 differ from "
1452 "product of matrix dimensions."));
1456 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1459 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1461 gsl_matrix_set (dst, y1, x1, *d);
1472 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1477 return gsl_matrix_alloc (0, 1);
1479 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1480 for (size_t y = 0; y < m->size1; y++)
1482 double ext = gsl_matrix_get (m, y, 0);
1483 for (size_t x = 1; x < m->size2; x++)
1485 double value = gsl_matrix_get (m, y, x);
1486 if (min ? value < ext : value > ext)
1489 gsl_matrix_set (rext, y, 0, ext);
1495 matrix_eval_RMAX (gsl_matrix *m)
1497 return matrix_eval_row_extremum (m, false);
1501 matrix_eval_RMIN (gsl_matrix *m)
1503 return matrix_eval_row_extremum (m, true);
1507 matrix_eval_RND (double d)
1519 rank_compare_3way (const void *a_, const void *b_)
1521 const struct rank *a = a_;
1522 const struct rank *b = b_;
1524 return a->value < b->value ? -1 : a->value > b->value;
1528 matrix_eval_RNKORDER (gsl_matrix *m)
1530 size_t n = m->size1 * m->size2;
1531 struct rank *ranks = xmalloc (n * sizeof *ranks);
1533 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1534 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1535 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1537 for (size_t i = 0; i < n; )
1540 for (j = i + 1; j < n; j++)
1541 if (ranks[i].value != ranks[j].value)
1544 double rank = (i + j + 1.0) / 2.0;
1545 for (size_t k = i; k < j; k++)
1546 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1557 matrix_eval_row_sum (gsl_matrix *m, bool square)
1562 return gsl_matrix_alloc (0, 1);
1564 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1565 for (size_t y = 0; y < m->size1; y++)
1568 for (size_t x = 0; x < m->size2; x++)
1570 double d = gsl_matrix_get (m, y, x);
1571 sum += square ? pow2 (d) : d;
1573 gsl_matrix_set (result, y, 0, sum);
1579 matrix_eval_RSSQ (gsl_matrix *m)
1581 return matrix_eval_row_sum (m, true);
1585 matrix_eval_RSUM (gsl_matrix *m)
1587 return matrix_eval_row_sum (m, false);
1591 matrix_eval_SIN (double d)
1597 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2)
1599 if (m1->size1 != m2->size1)
1601 msg (SE, _("SOLVE requires its arguments to have the same number of "
1602 "rows, but the first argument has dimensions %zu×%zu and "
1603 "the second %zu×%zu."),
1604 m1->size1, m1->size2,
1605 m2->size1, m2->size2);
1609 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1610 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1612 gsl_linalg_LU_decomp (m1, p, &signum);
1613 for (size_t i = 0; i < m2->size2; i++)
1615 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1616 gsl_vector xi = gsl_matrix_column (x, i).vector;
1617 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1619 gsl_permutation_free (p);
1624 matrix_eval_SQRT (gsl_matrix *m)
1626 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1630 msg (SE, _("Argument to SQRT must be nonnegative."));
1639 matrix_eval_SSCP (gsl_matrix *m)
1641 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1642 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1647 matrix_eval_SVAL (gsl_matrix *m)
1649 gsl_matrix *tmp_mat = NULL;
1650 if (m->size2 > m->size1)
1652 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1653 gsl_matrix_transpose_memcpy (tmp_mat, m);
1658 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1659 gsl_vector *S = gsl_vector_alloc (m->size2);
1660 gsl_vector *work = gsl_vector_alloc (m->size2);
1661 gsl_linalg_SV_decomp (m, V, S, work);
1663 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1664 for (size_t i = 0; i < m->size2; i++)
1665 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1667 gsl_matrix_free (V);
1668 gsl_vector_free (S);
1669 gsl_vector_free (work);
1670 gsl_matrix_free (tmp_mat);
1676 matrix_eval_SWEEP (gsl_matrix *m, double d)
1678 if (d < 1 || d > SIZE_MAX)
1680 msg (SE, _("Scalar argument to SWEEP must be integer."));
1684 if (k >= MIN (m->size1, m->size2))
1686 msg (SE, _("Scalar argument to SWEEP must be integer less than or "
1687 "equal to the smaller of the matrix argument's rows and "
1692 double m_kk = gsl_matrix_get (m, k, k);
1693 if (fabs (m_kk) > 1e-19)
1695 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1696 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1698 double m_ij = gsl_matrix_get (m, i, j);
1699 double m_ik = gsl_matrix_get (m, i, k);
1700 double m_kj = gsl_matrix_get (m, k, j);
1701 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1710 for (size_t i = 0; i < m->size1; i++)
1712 gsl_matrix_set (m, i, k, 0);
1713 gsl_matrix_set (m, k, i, 0);
1720 matrix_eval_TRACE (gsl_matrix *m)
1723 size_t n = MIN (m->size1, m->size2);
1724 for (size_t i = 0; i < n; i++)
1725 sum += gsl_matrix_get (m, i, i);
1730 matrix_eval_T (gsl_matrix *m)
1732 return matrix_eval_TRANSPOS (m);
1736 matrix_eval_TRANSPOS (gsl_matrix *m)
1738 if (m->size1 == m->size2)
1740 gsl_matrix_transpose (m);
1745 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1746 gsl_matrix_transpose_memcpy (t, m);
1752 matrix_eval_TRUNC (double d)
1758 matrix_eval_UNIFORM (double r_, double c_)
1760 if (r_ < 0 || r_ >= SIZE_MAX || c_ < 0 || c_ >= SIZE_MAX)
1762 msg (SE, _("Arguments to UNIFORM must be integers."));
1767 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1769 msg (SE, _("Product of arguments to UNIFORM exceeds memory size."));
1773 gsl_matrix *m = gsl_matrix_alloc (r, c);
1774 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1775 *d = gsl_ran_flat (get_rng (), 0, 1);
1780 matrix_eval_PDF_BETA (double x, double a, double b)
1782 return gsl_ran_beta_pdf (x, a, b);
1786 matrix_eval_CDF_BETA (double x, double a, double b)
1788 return gsl_cdf_beta_P (x, a, b);
1792 matrix_eval_IDF_BETA (double P, double a, double b)
1794 return gsl_cdf_beta_Pinv (P, a, b);
1798 matrix_eval_RV_BETA (double a, double b)
1800 return gsl_ran_beta (get_rng (), a, b);
1804 matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
1806 return ncdf_beta (x, a, b, lambda);
1810 matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
1812 return npdf_beta (x, a, b, lambda);
1816 matrix_eval_CDF_BVNOR (double x0, double x1, double r)
1818 return cdf_bvnor (x0, x1, r);
1822 matrix_eval_PDF_BVNOR (double x0, double x1, double r)
1824 return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
1828 matrix_eval_CDF_CAUCHY (double x, double a, double b)
1830 return gsl_cdf_cauchy_P ((x - a) / b, 1);
1834 matrix_eval_IDF_CAUCHY (double P, double a, double b)
1836 return a + b * gsl_cdf_cauchy_Pinv (P, 1);
1840 matrix_eval_PDF_CAUCHY (double x, double a, double b)
1842 return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
1846 matrix_eval_RV_CAUCHY (double a, double b)
1848 return a + b * gsl_ran_cauchy (get_rng (), 1);
1852 matrix_eval_CDF_CHISQ (double x, double df)
1854 return gsl_cdf_chisq_P (x, df);
1858 matrix_eval_CHICDF (double x, double df)
1860 return matrix_eval_CDF_CHISQ (x, df);
1864 matrix_eval_IDF_CHISQ (double P, double df)
1866 return gsl_cdf_chisq_Pinv (P, df);
1870 matrix_eval_PDF_CHISQ (double x, double df)
1872 return gsl_ran_chisq_pdf (x, df);
1876 matrix_eval_RV_CHISQ (double df)
1878 return gsl_ran_chisq (get_rng (), df);
1882 matrix_eval_SIG_CHISQ (double x, double df)
1884 return gsl_cdf_chisq_Q (x, df);
1888 matrix_eval_CDF_EXP (double x, double a)
1890 return gsl_cdf_exponential_P (x, 1. / a);
1894 matrix_eval_IDF_EXP (double P, double a)
1896 return gsl_cdf_exponential_Pinv (P, 1. / a);
1900 matrix_eval_PDF_EXP (double x, double a)
1902 return gsl_ran_exponential_pdf (x, 1. / a);
1906 matrix_eval_RV_EXP (double a)
1908 return gsl_ran_exponential (get_rng (), 1. / a);
1912 matrix_eval_PDF_XPOWER (double x, double a, double b)
1914 return gsl_ran_exppow_pdf (x, a, b);
1918 matrix_eval_RV_XPOWER (double a, double b)
1920 return gsl_ran_exppow (get_rng (), a, b);
1924 matrix_eval_CDF_F (double x, double df1, double df2)
1926 return gsl_cdf_fdist_P (x, df1, df2);
1930 matrix_eval_FCDF (double x, double df1, double df2)
1932 return matrix_eval_CDF_F (x, df1, df2);
1936 matrix_eval_IDF_F (double P, double df1, double df2)
1938 return idf_fdist (P, df1, df2);
1942 matrix_eval_RV_F (double df1, double df2)
1944 return gsl_ran_fdist (get_rng (), df1, df2);
1948 matrix_eval_PDF_F (double x, double df1, double df2)
1950 return gsl_ran_fdist_pdf (x, df1, df2);
1954 matrix_eval_SIG_F (double x, double df1, double df2)
1956 return gsl_cdf_fdist_Q (x, df1, df2);
1960 matrix_eval_CDF_GAMMA (double x, double a, double b)
1962 return gsl_cdf_gamma_P (x, a, 1. / b);
1966 matrix_eval_IDF_GAMMA (double P, double a, double b)
1968 return gsl_cdf_gamma_Pinv (P, a, 1. / b);
1972 matrix_eval_PDF_GAMMA (double x, double a, double b)
1974 return gsl_ran_gamma_pdf (x, a, 1. / b);
1978 matrix_eval_RV_GAMMA (double a, double b)
1980 return gsl_ran_gamma (get_rng (), a, 1. / b);
1984 matrix_eval_PDF_LANDAU (double x)
1986 return gsl_ran_landau_pdf (x);
1990 matrix_eval_RV_LANDAU (void)
1992 return gsl_ran_landau (get_rng ());
1996 matrix_eval_CDF_LAPLACE (double x, double a, double b)
1998 return gsl_cdf_laplace_P ((x - a) / b, 1);
2002 matrix_eval_IDF_LAPLACE (double P, double a, double b)
2004 return a + b * gsl_cdf_laplace_Pinv (P, 1);
2008 matrix_eval_PDF_LAPLACE (double x, double a, double b)
2010 return gsl_ran_laplace_pdf ((x - a) / b, 1);
2014 matrix_eval_RV_LAPLACE (double a, double b)
2016 return a + b * gsl_ran_laplace (get_rng (), 1);
2020 matrix_eval_RV_LEVY (double c, double alpha)
2022 return gsl_ran_levy (get_rng (), c, alpha);
2026 matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
2028 return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
2032 matrix_eval_CDF_LOGISTIC (double x, double a, double b)
2034 return gsl_cdf_logistic_P ((x - a) / b, 1);
2038 matrix_eval_IDF_LOGISTIC (double P, double a, double b)
2040 return a + b * gsl_cdf_logistic_Pinv (P, 1);
2044 matrix_eval_PDF_LOGISTIC (double x, double a, double b)
2046 return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
2050 matrix_eval_RV_LOGISTIC (double a, double b)
2052 return a + b * gsl_ran_logistic (get_rng (), 1);
2056 matrix_eval_CDF_LNORMAL (double x, double m, double s)
2058 return gsl_cdf_lognormal_P (x, log (m), s);
2062 matrix_eval_IDF_LNORMAL (double P, double m, double s)
2064 return gsl_cdf_lognormal_Pinv (P, log (m), s);;
2068 matrix_eval_PDF_LNORMAL (double x, double m, double s)
2070 return gsl_ran_lognormal_pdf (x, log (m), s);
2074 matrix_eval_RV_LNORMAL (double m, double s)
2076 return gsl_ran_lognormal (get_rng (), log (m), s);
2080 matrix_eval_CDF_NORMAL (double x, double u, double s)
2082 return gsl_cdf_gaussian_P (x - u, s);
2086 matrix_eval_IDF_NORMAL (double P, double u, double s)
2088 return u + gsl_cdf_gaussian_Pinv (P, s);
2092 matrix_eval_PDF_NORMAL (double x, double u, double s)
2094 return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
2098 matrix_eval_RV_NORMAL (double u, double s)
2100 return u + gsl_ran_gaussian (get_rng (), s);
2104 matrix_eval_CDFNORM (double x)
2106 return gsl_cdf_ugaussian_P (x);
2110 matrix_eval_PROBIT (double P)
2112 return gsl_cdf_ugaussian_Pinv (P);
2116 matrix_eval_NORMAL (double s)
2118 return gsl_ran_gaussian (get_rng (), s);
2122 matrix_eval_PDF_NTAIL (double x, double a, double sigma)
2124 return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
2128 matrix_eval_RV_NTAIL (double a, double sigma)
2130 return gsl_ran_gaussian_tail (get_rng (), a, sigma);
2134 matrix_eval_CDF_PARETO (double x, double a, double b)
2136 return gsl_cdf_pareto_P (x, b, a);
2140 matrix_eval_IDF_PARETO (double P, double a, double b)
2142 return gsl_cdf_pareto_Pinv (P, b, a);
2146 matrix_eval_PDF_PARETO (double x, double a, double b)
2148 return gsl_ran_pareto_pdf (x, b, a);
2152 matrix_eval_RV_PARETO (double a, double b)
2154 return gsl_ran_pareto (get_rng (), b, a);
2158 matrix_eval_CDF_RAYLEIGH (double x, double sigma)
2160 return gsl_cdf_rayleigh_P (x, sigma);
2164 matrix_eval_IDF_RAYLEIGH (double P, double sigma)
2166 return gsl_cdf_rayleigh_Pinv (P, sigma);
2170 matrix_eval_PDF_RAYLEIGH (double x, double sigma)
2172 return gsl_ran_rayleigh_pdf (x, sigma);
2176 matrix_eval_RV_RAYLEIGH (double sigma)
2178 return gsl_ran_rayleigh (get_rng (), sigma);
2182 matrix_eval_PDF_RTAIL (double x, double a, double sigma)
2184 return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
2188 matrix_eval_RV_RTAIL (double a, double sigma)
2190 return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
2194 matrix_eval_CDF_T (double x, double df)
2196 return gsl_cdf_tdist_P (x, df);
2200 matrix_eval_TCDF (double x, double df)
2202 return matrix_eval_CDF_T (x, df);
2206 matrix_eval_IDF_T (double P, double df)
2208 return gsl_cdf_tdist_Pinv (P, df);
2212 matrix_eval_PDF_T (double x, double df)
2214 return gsl_ran_tdist_pdf (x, df);
2218 matrix_eval_RV_T (double df)
2220 return gsl_ran_tdist (get_rng (), df);
2224 matrix_eval_CDF_T1G (double x, double a, double b)
2226 return gsl_cdf_gumbel1_P (x, a, b);
2230 matrix_eval_IDF_T1G (double P, double a, double b)
2232 return gsl_cdf_gumbel1_Pinv (P, a, b);
2236 matrix_eval_PDF_T1G (double x, double a, double b)
2238 return gsl_ran_gumbel1_pdf (x, a, b);
2242 matrix_eval_RV_T1G (double a, double b)
2244 return gsl_ran_gumbel1 (get_rng (), a, b);
2248 matrix_eval_CDF_T2G (double x, double a, double b)
2250 return gsl_cdf_gumbel1_P (x, a, b);
2254 matrix_eval_IDF_T2G (double P, double a, double b)
2256 return gsl_cdf_gumbel1_Pinv (P, a, b);
2260 matrix_eval_PDF_T2G (double x, double a, double b)
2262 return gsl_ran_gumbel1_pdf (x, a, b);
2266 matrix_eval_RV_T2G (double a, double b)
2268 return gsl_ran_gumbel1 (get_rng (), a, b);
2272 matrix_eval_CDF_UNIFORM (double x, double a, double b)
2274 return gsl_cdf_flat_P (x, a, b);
2278 matrix_eval_IDF_UNIFORM (double P, double a, double b)
2280 return gsl_cdf_flat_Pinv (P, a, b);
2284 matrix_eval_PDF_UNIFORM (double x, double a, double b)
2286 return gsl_ran_flat_pdf (x, a, b);
2290 matrix_eval_RV_UNIFORM (double a, double b)
2292 return gsl_ran_flat (get_rng (), a, b);
2296 matrix_eval_CDF_WEIBULL (double x, double a, double b)
2298 return gsl_cdf_weibull_P (x, a, b);
2302 matrix_eval_IDF_WEIBULL (double P, double a, double b)
2304 return gsl_cdf_weibull_Pinv (P, a, b);
2308 matrix_eval_PDF_WEIBULL (double x, double a, double b)
2310 return gsl_ran_weibull_pdf (x, a, b);
2314 matrix_eval_RV_WEIBULL (double a, double b)
2316 return gsl_ran_weibull (get_rng (), a, b);
2320 matrix_eval_CDF_BERNOULLI (double k, double p)
2322 return k ? 1 : 1 - p;
2326 matrix_eval_PDF_BERNOULLI (double k, double p)
2328 return gsl_ran_bernoulli_pdf (k, p);
2332 matrix_eval_RV_BERNOULLI (double p)
2334 return gsl_ran_bernoulli (get_rng (), p);
2338 matrix_eval_CDF_BINOM (double k, double n, double p)
2340 return gsl_cdf_binomial_P (k, p, n);
2344 matrix_eval_PDF_BINOM (double k, double n, double p)
2346 return gsl_ran_binomial_pdf (k, p, n);
2350 matrix_eval_RV_BINOM (double n, double p)
2352 return gsl_ran_binomial (get_rng (), p, n);
2356 matrix_eval_CDF_GEOM (double k, double p)
2358 return gsl_cdf_geometric_P (k, p);
2362 matrix_eval_PDF_GEOM (double k, double p)
2364 return gsl_ran_geometric_pdf (k, p);
2368 matrix_eval_RV_GEOM (double p)
2370 return gsl_ran_geometric (get_rng (), p);
2374 matrix_eval_CDF_HYPER (double k, double a, double b, double c)
2376 return gsl_cdf_hypergeometric_P (k, c, a - c, b);
2380 matrix_eval_PDF_HYPER (double k, double a, double b, double c)
2382 return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
2386 matrix_eval_RV_HYPER (double a, double b, double c)
2388 return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
2392 matrix_eval_PDF_LOG (double k, double p)
2394 return gsl_ran_logarithmic_pdf (k, p);
2398 matrix_eval_RV_LOG (double p)
2400 return gsl_ran_logarithmic (get_rng (), p);
2404 matrix_eval_CDF_NEGBIN (double k, double n, double p)
2406 return gsl_cdf_negative_binomial_P (k, p, n);
2410 matrix_eval_PDF_NEGBIN (double k, double n, double p)
2412 return gsl_ran_negative_binomial_pdf (k, p, n);
2416 matrix_eval_RV_NEGBIN (double n, double p)
2418 return gsl_ran_negative_binomial (get_rng (), p, n);
2422 matrix_eval_CDF_POISSON (double k, double mu)
2424 return gsl_cdf_poisson_P (k, mu);
2428 matrix_eval_PDF_POISSON (double k, double mu)
2430 return gsl_ran_poisson_pdf (k, mu);
2434 matrix_eval_RV_POISSON (double mu)
2436 return gsl_ran_poisson (get_rng (), mu);
2439 struct matrix_function
2443 size_t min_args, max_args;
2446 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
2449 word_matches (const char **test, const char **name)
2451 size_t test_len = strcspn (*test, ".");
2452 size_t name_len = strcspn (*name, ".");
2453 if (test_len == name_len)
2455 if (buf_compare_case (*test, *name, test_len))
2458 else if (test_len < 3 || test_len > name_len)
2462 if (buf_compare_case (*test, *name, test_len))
2468 if (**test != **name)
2479 /* Returns 0 if TOKEN and FUNC do not match,
2480 1 if TOKEN is an acceptable abbreviation for FUNC,
2481 2 if TOKEN equals FUNC. */
2483 compare_function_names (const char *token_, const char *func_)
2485 const char *token = token_;
2486 const char *func = func_;
2487 while (*token || *func)
2488 if (!word_matches (&token, &func))
2490 return !c_strcasecmp (token_, func_) ? 2 : 1;
2493 static const struct matrix_function *
2494 matrix_parse_function_name (const char *token)
2496 static const struct matrix_function functions[] =
2498 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
2499 { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
2503 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
2505 for (size_t i = 0; i < N_FUNCTIONS; i++)
2507 if (compare_function_names (token, functions[i].name) > 0)
2508 return &functions[i];
2513 static struct read_file *
2514 read_file_create (struct matrix_state *s, struct file_handle *fh)
2516 for (size_t i = 0; i < s->n_read_files; i++)
2518 struct read_file *rf = s->read_files[i];
2526 struct read_file *rf = xmalloc (sizeof *rf);
2527 *rf = (struct read_file) { .file = fh };
2529 s->read_files = xrealloc (s->read_files,
2530 (s->n_read_files + 1) * sizeof *s->read_files);
2531 s->read_files[s->n_read_files++] = rf;
2535 static struct dfm_reader *
2536 read_file_open (struct read_file *rf)
2539 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
2544 read_file_destroy (struct read_file *rf)
2548 fh_unref (rf->file);
2549 dfm_close_reader (rf->reader);
2550 free (rf->encoding);
2555 static struct write_file *
2556 write_file_create (struct matrix_state *s, struct file_handle *fh)
2558 for (size_t i = 0; i < s->n_write_files; i++)
2560 struct write_file *wf = s->write_files[i];
2568 struct write_file *wf = xmalloc (sizeof *wf);
2569 *wf = (struct write_file) { .file = fh };
2571 s->write_files = xrealloc (s->write_files,
2572 (s->n_write_files + 1) * sizeof *s->write_files);
2573 s->write_files[s->n_write_files++] = wf;
2577 static struct dfm_writer *
2578 write_file_open (struct write_file *wf)
2581 wf->writer = dfm_open_writer (wf->file, wf->encoding);
2586 write_file_destroy (struct write_file *wf)
2592 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
2593 wf->held->s.ss.length);
2594 u8_line_destroy (wf->held);
2598 fh_unref (wf->file);
2599 dfm_close_writer (wf->writer);
2600 free (wf->encoding);
2606 matrix_parse_function (struct matrix_state *s, const char *token,
2607 struct matrix_expr **exprp)
2610 if (lex_next_token (s->lexer, 1) != T_LPAREN)
2613 if (lex_match_id (s->lexer, "EOF"))
2616 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
2620 if (!lex_force_match (s->lexer, T_RPAREN))
2626 struct read_file *rf = read_file_create (s, fh);
2628 struct matrix_expr *e = xmalloc (sizeof *e);
2629 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
2634 const struct matrix_function *f = matrix_parse_function_name (token);
2638 lex_get_n (s->lexer, 2);
2640 struct matrix_expr *e = xmalloc (sizeof *e);
2641 *e = (struct matrix_expr) { .op = f->op, .subs = NULL };
2643 if (lex_token (s->lexer) != T_RPAREN)
2645 size_t allocated_subs = 0;
2648 struct matrix_expr *sub = matrix_parse_expr (s);
2652 if (e->n_subs >= allocated_subs)
2653 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
2654 e->subs[e->n_subs++] = sub;
2656 while (lex_match (s->lexer, T_COMMA));
2658 if (!lex_force_match (s->lexer, T_RPAREN))
2661 if (e->n_subs < f->min_args || e->n_subs > f->max_args)
2663 if (f->min_args == f->max_args)
2664 msg (SE, ngettext ("Matrix function %s requires %zu argument.",
2665 "Matrix function %s requires %zu arguments.",
2667 f->name, f->min_args);
2668 else if (f->min_args == 1 && f->max_args == 2)
2669 msg (SE, ngettext ("Matrix function %s requires 1 or 2 arguments, "
2670 "but %zu was provided.",
2671 "Matrix function %s requires 1 or 2 arguments, "
2672 "but %zu were provided.",
2674 f->name, e->n_subs);
2675 else if (f->min_args == 1 && f->max_args == INT_MAX)
2676 msg (SE, _("Matrix function %s requires at least one argument."),
2688 matrix_expr_destroy (e);
2692 static struct matrix_expr *
2693 matrix_parse_primary (struct matrix_state *s)
2695 if (lex_is_number (s->lexer))
2697 double number = lex_number (s->lexer);
2699 return matrix_expr_create_number (number);
2701 else if (lex_is_string (s->lexer))
2703 char string[sizeof (double)];
2704 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
2708 memcpy (&number, string, sizeof number);
2709 return matrix_expr_create_number (number);
2711 else if (lex_match (s->lexer, T_LPAREN))
2713 struct matrix_expr *e = matrix_parse_or_xor (s);
2714 if (!e || !lex_force_match (s->lexer, T_RPAREN))
2716 matrix_expr_destroy (e);
2721 else if (lex_match (s->lexer, T_LCURLY))
2723 struct matrix_expr *e = matrix_parse_curly_semi (s);
2724 if (!e || !lex_force_match (s->lexer, T_RCURLY))
2726 matrix_expr_destroy (e);
2731 else if (lex_token (s->lexer) == T_ID)
2733 struct matrix_expr *retval;
2734 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
2737 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
2740 lex_error (s->lexer, _("Unknown variable %s."),
2741 lex_tokcstr (s->lexer));
2746 struct matrix_expr *e = xmalloc (sizeof *e);
2747 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
2750 else if (lex_token (s->lexer) == T_ALL)
2752 struct matrix_expr *retval;
2753 if (matrix_parse_function (s, "ALL", &retval))
2757 lex_error (s->lexer, NULL);
2761 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
2764 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
2766 if (lex_match (s->lexer, T_COLON))
2773 *indexp = matrix_parse_expr (s);
2774 return *indexp != NULL;
2778 static struct matrix_expr *
2779 matrix_parse_postfix (struct matrix_state *s)
2781 struct matrix_expr *lhs = matrix_parse_primary (s);
2782 if (!lhs || !lex_match (s->lexer, T_LPAREN))
2785 struct matrix_expr *i0;
2786 if (!matrix_parse_index_expr (s, &i0))
2788 matrix_expr_destroy (lhs);
2791 if (lex_match (s->lexer, T_RPAREN))
2793 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
2794 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
2795 else if (lex_match (s->lexer, T_COMMA))
2797 struct matrix_expr *i1;
2798 if (!matrix_parse_index_expr (s, &i1)
2799 || !lex_force_match (s->lexer, T_RPAREN))
2801 matrix_expr_destroy (lhs);
2802 matrix_expr_destroy (i0);
2803 matrix_expr_destroy (i1);
2806 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
2807 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
2808 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
2813 lex_error_expecting (s->lexer, "`)'", "`,'");
2818 static struct matrix_expr *
2819 matrix_parse_unary (struct matrix_state *s)
2821 if (lex_match (s->lexer, T_DASH))
2823 struct matrix_expr *lhs = matrix_parse_unary (s);
2824 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
2826 else if (lex_match (s->lexer, T_PLUS))
2827 return matrix_parse_unary (s);
2829 return matrix_parse_postfix (s);
2832 static struct matrix_expr *
2833 matrix_parse_seq (struct matrix_state *s)
2835 struct matrix_expr *start = matrix_parse_unary (s);
2836 if (!start || !lex_match (s->lexer, T_COLON))
2839 struct matrix_expr *end = matrix_parse_unary (s);
2842 matrix_expr_destroy (start);
2846 if (lex_match (s->lexer, T_COLON))
2848 struct matrix_expr *increment = matrix_parse_unary (s);
2851 matrix_expr_destroy (start);
2852 matrix_expr_destroy (end);
2855 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2858 return matrix_expr_create_2 (MOP_SEQ, start, end);
2861 static struct matrix_expr *
2862 matrix_parse_exp (struct matrix_state *s)
2864 struct matrix_expr *lhs = matrix_parse_seq (s);
2871 if (lex_match (s->lexer, T_EXP))
2873 else if (lex_match_phrase (s->lexer, "&**"))
2878 struct matrix_expr *rhs = matrix_parse_seq (s);
2881 matrix_expr_destroy (lhs);
2884 lhs = matrix_expr_create_2 (op, lhs, rhs);
2888 static struct matrix_expr *
2889 matrix_parse_mul_div (struct matrix_state *s)
2891 struct matrix_expr *lhs = matrix_parse_exp (s);
2898 if (lex_match (s->lexer, T_ASTERISK))
2900 else if (lex_match (s->lexer, T_SLASH))
2902 else if (lex_match_phrase (s->lexer, "&*"))
2904 else if (lex_match_phrase (s->lexer, "&/"))
2909 struct matrix_expr *rhs = matrix_parse_exp (s);
2912 matrix_expr_destroy (lhs);
2915 lhs = matrix_expr_create_2 (op, lhs, rhs);
2919 static struct matrix_expr *
2920 matrix_parse_add_sub (struct matrix_state *s)
2922 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2929 if (lex_match (s->lexer, T_PLUS))
2931 else if (lex_match (s->lexer, T_DASH))
2933 else if (lex_token (s->lexer) == T_NEG_NUM)
2938 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2941 matrix_expr_destroy (lhs);
2944 lhs = matrix_expr_create_2 (op, lhs, rhs);
2948 static struct matrix_expr *
2949 matrix_parse_relational (struct matrix_state *s)
2951 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2958 if (lex_match (s->lexer, T_GT))
2960 else if (lex_match (s->lexer, T_GE))
2962 else if (lex_match (s->lexer, T_LT))
2964 else if (lex_match (s->lexer, T_LE))
2966 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2968 else if (lex_match (s->lexer, T_NE))
2973 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2976 matrix_expr_destroy (lhs);
2979 lhs = matrix_expr_create_2 (op, lhs, rhs);
2983 static struct matrix_expr *
2984 matrix_parse_not (struct matrix_state *s)
2986 if (lex_match (s->lexer, T_NOT))
2988 struct matrix_expr *lhs = matrix_parse_not (s);
2989 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2992 return matrix_parse_relational (s);
2995 static struct matrix_expr *
2996 matrix_parse_and (struct matrix_state *s)
2998 struct matrix_expr *lhs = matrix_parse_not (s);
3002 while (lex_match (s->lexer, T_AND))
3004 struct matrix_expr *rhs = matrix_parse_not (s);
3007 matrix_expr_destroy (lhs);
3010 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
3015 static struct matrix_expr *
3016 matrix_parse_or_xor (struct matrix_state *s)
3018 struct matrix_expr *lhs = matrix_parse_and (s);
3025 if (lex_match (s->lexer, T_OR))
3027 else if (lex_match_id (s->lexer, "XOR"))
3032 struct matrix_expr *rhs = matrix_parse_and (s);
3035 matrix_expr_destroy (lhs);
3038 lhs = matrix_expr_create_2 (op, lhs, rhs);
3042 static struct matrix_expr *
3043 matrix_parse_expr (struct matrix_state *s)
3045 return matrix_parse_or_xor (s);
3048 /* Expression evaluation. */
3051 matrix_op_eval (enum matrix_op op, double a, double b)
3055 case MOP_ADD_ELEMS: return a + b;
3056 case MOP_SUB_ELEMS: return a - b;
3057 case MOP_MUL_ELEMS: return a * b;
3058 case MOP_DIV_ELEMS: return a / b;
3059 case MOP_EXP_ELEMS: return pow (a, b);
3060 case MOP_GT: return a > b;
3061 case MOP_GE: return a >= b;
3062 case MOP_LT: return a < b;
3063 case MOP_LE: return a <= b;
3064 case MOP_EQ: return a == b;
3065 case MOP_NE: return a != b;
3066 case MOP_AND: return (a > 0) && (b > 0);
3067 case MOP_OR: return (a > 0) || (b > 0);
3068 case MOP_XOR: return (a > 0) != (b > 0);
3070 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3079 case MOP_PASTE_HORZ:
3080 case MOP_PASTE_VERT:
3096 matrix_op_name (enum matrix_op op)
3100 case MOP_ADD_ELEMS: return "+";
3101 case MOP_SUB_ELEMS: return "-";
3102 case MOP_MUL_ELEMS: return "&*";
3103 case MOP_DIV_ELEMS: return "&/";
3104 case MOP_EXP_ELEMS: return "&**";
3105 case MOP_GT: return ">";
3106 case MOP_GE: return ">=";
3107 case MOP_LT: return "<";
3108 case MOP_LE: return "<=";
3109 case MOP_EQ: return "=";
3110 case MOP_NE: return "<>";
3111 case MOP_AND: return "AND";
3112 case MOP_OR: return "OR";
3113 case MOP_XOR: return "XOR";
3115 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3124 case MOP_PASTE_HORZ:
3125 case MOP_PASTE_VERT:
3141 is_scalar (const gsl_matrix *m)
3143 return m->size1 == 1 && m->size2 == 1;
3147 to_scalar (const gsl_matrix *m)
3149 assert (is_scalar (m));
3150 return gsl_matrix_get (m, 0, 0);
3154 matrix_expr_evaluate_elementwise (enum matrix_op op,
3155 gsl_matrix *a, gsl_matrix *b)
3159 double be = to_scalar (b);
3160 for (size_t r = 0; r < a->size1; r++)
3161 for (size_t c = 0; c < a->size2; c++)
3163 double *ae = gsl_matrix_ptr (a, r, c);
3164 *ae = matrix_op_eval (op, *ae, be);
3168 else if (is_scalar (a))
3170 double ae = to_scalar (a);
3171 for (size_t r = 0; r < b->size1; r++)
3172 for (size_t c = 0; c < b->size2; c++)
3174 double *be = gsl_matrix_ptr (b, r, c);
3175 *be = matrix_op_eval (op, ae, *be);
3179 else if (a->size1 == b->size1 && a->size2 == b->size2)
3181 for (size_t r = 0; r < a->size1; r++)
3182 for (size_t c = 0; c < a->size2; c++)
3184 double *ae = gsl_matrix_ptr (a, r, c);
3185 double be = gsl_matrix_get (b, r, c);
3186 *ae = matrix_op_eval (op, *ae, be);
3192 msg (SE, _("Operands to %s must have the same dimensions or one "
3193 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
3194 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
3200 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
3202 if (is_scalar (a) || is_scalar (b))
3203 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
3205 if (a->size2 != b->size1)
3207 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
3208 "not conformable for multiplication."),
3209 a->size1, a->size2, b->size1, b->size2);
3213 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3214 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3219 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3221 gsl_matrix *tmp = *a;
3227 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3230 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3231 swap_matrix (z, tmp);
3235 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3237 mul_matrix (x, *x, *x, tmp);
3241 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
3244 if (x->size1 != x->size2)
3246 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
3247 "the left-hand size, not one with dimensions %zu×%zu."),
3248 x->size1, x->size2);
3253 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
3254 "right-hand side, not a matrix with dimensions %zu×%zu."),
3255 b->size1, b->size2);
3258 double bf = to_scalar (b);
3259 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3261 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
3262 "or outside the valid range."), bf);
3267 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3269 gsl_matrix_set_identity (y);
3273 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3275 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3278 mul_matrix (&y, x, y, &t);
3279 square_matrix (&x, &t);
3282 square_matrix (&x, &t);
3284 mul_matrix (&y, x, y, &t);
3288 /* Garbage collection.
3290 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3291 a permutation of them. We are returning one of them; that one must not be
3292 destroyed. We must not destroy 'x_' because the caller owns it. */
3294 gsl_matrix_free (y_);
3296 gsl_matrix_free (t_);
3302 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
3305 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3307 msg (SE, _("All operands of : operator must be scalars."));
3311 long int start = to_scalar (start_);
3312 long int end = to_scalar (end_);
3313 long int by = by_ ? to_scalar (by_) : 1;
3317 msg (SE, _("The increment operand to : must be nonzero."));
3321 long int n = (end >= start && by > 0 ? (end - start + by) / by
3322 : end <= start && by < 0 ? (start - end - by) / -by
3324 gsl_matrix *m = gsl_matrix_alloc (1, n);
3325 for (long int i = 0; i < n; i++)
3326 gsl_matrix_set (m, 0, i, start + i * by);
3331 matrix_expr_evaluate_not (gsl_matrix *a)
3333 for (size_t r = 0; r < a->size1; r++)
3334 for (size_t c = 0; c < a->size2; c++)
3336 double *ae = gsl_matrix_ptr (a, r, c);
3343 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
3345 if (a->size1 != b->size1)
3347 if (!a->size1 || !a->size2)
3349 else if (!b->size1 || !b->size2)
3352 msg (SE, _("All columns in a matrix must have the same number of rows, "
3353 "but this tries to paste matrices with %zu and %zu rows."),
3354 a->size1, b->size1);
3358 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3359 for (size_t y = 0; y < a->size1; y++)
3361 for (size_t x = 0; x < a->size2; x++)
3362 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3363 for (size_t x = 0; x < b->size2; x++)
3364 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3370 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
3372 if (a->size2 != b->size2)
3374 if (!a->size1 || !a->size2)
3376 else if (!b->size1 || !b->size2)
3379 msg (SE, _("All rows in a matrix must have the same number of columns, "
3380 "but this tries to stack matrices with %zu and %zu columns."),
3381 a->size2, b->size2);
3385 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3386 for (size_t x = 0; x < a->size2; x++)
3388 for (size_t y = 0; y < a->size1; y++)
3389 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3390 for (size_t y = 0; y < b->size1; y++)
3391 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3397 matrix_to_vector (gsl_matrix *m)
3400 gsl_vector v = to_vector (m);
3401 assert (v.block == m->block || !v.block);
3405 gsl_matrix_free (m);
3406 return xmemdup (&v, sizeof v);
3420 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3423 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
3424 enum index_type index_type, size_t other_size,
3425 struct index_vector *iv)
3434 msg (SE, _("Vector index must be scalar or vector, not a "
3436 m->size1, m->size2);
3440 msg (SE, _("Matrix row index must be scalar or vector, not a "
3442 m->size1, m->size2);
3446 msg (SE, _("Matrix column index must be scalar or vector, not a "
3448 m->size1, m->size2);
3454 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3455 *iv = (struct index_vector) {
3456 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3459 for (size_t i = 0; i < v.size; i++)
3461 double index = gsl_vector_get (&v, i);
3462 if (index < 1 || index >= size + 1)
3467 msg (SE, _("Index %g is out of range for vector "
3468 "with %zu elements."), index, size);
3472 msg (SE, _("%g is not a valid row index for "
3473 "a %zu×%zu matrix."),
3474 index, size, other_size);
3478 msg (SE, _("%g is not a valid column index for "
3479 "a %zu×%zu matrix."),
3480 index, other_size, size);
3487 iv->indexes[i] = index - 1;
3493 *iv = (struct index_vector) {
3494 .indexes = xnmalloc (size, sizeof *iv->indexes),
3497 for (size_t i = 0; i < size; i++)
3504 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
3506 if (!is_vector (sm))
3508 msg (SE, _("Vector index operator may not be applied to "
3509 "a %zu×%zu matrix."),
3510 sm->size1, sm->size2);
3518 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
3520 if (!matrix_expr_evaluate_vec_all (sm))
3523 gsl_vector sv = to_vector (sm);
3524 struct index_vector iv;
3525 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
3528 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3529 sm->size1 == 1 ? iv.n : 1);
3530 gsl_vector dv = to_vector (dm);
3531 for (size_t dx = 0; dx < iv.n; dx++)
3533 size_t sx = iv.indexes[dx];
3534 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3542 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
3545 struct index_vector iv0;
3546 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
3549 struct index_vector iv1;
3550 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
3557 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3558 for (size_t dy = 0; dy < iv0.n; dy++)
3560 size_t sy = iv0.indexes[dy];
3562 for (size_t dx = 0; dx < iv1.n; dx++)
3564 size_t sx = iv1.indexes[dx];
3565 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3573 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3574 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3575 const struct matrix_function_properties *, gsl_matrix *[], size_t, \
3576 matrix_proto_##PROTO *);
3581 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
3583 if (!is_scalar (subs[index]))
3585 msg (SE, _("Function %s argument %zu must be a scalar, "
3586 "not a %zu×%zu matrix."),
3587 name, index + 1, subs[index]->size1, subs[index]->size2);
3594 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
3596 if (!is_vector (subs[index]))
3598 msg (SE, _("Function %s argument %zu must be a vector, "
3599 "not a %zu×%zu matrix."),
3600 name, index + 1, subs[index]->size1, subs[index]->size2);
3607 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
3609 for (size_t i = 0; i < n_subs; i++)
3611 if (!check_scalar_arg (name, subs, i))
3613 d[i] = to_scalar (subs[i]);
3619 parse_constraint_value (const char **constraintsp)
3622 long retval = strtol (*constraintsp, &tail, 10);
3623 assert (tail > *constraintsp);
3624 *constraintsp = tail;
3629 argument_invalid (const struct matrix_function_properties *props,
3630 size_t arg_index, gsl_matrix *a, size_t y, size_t x,
3634 ds_put_format (out, _("Argument %zu to matrix function %s "
3635 "has invalid value %g."),
3636 arg_index, props->name, gsl_matrix_get (a, y, x));
3638 ds_put_format (out, _("Row %zu, column %zu of argument %zu "
3639 "to matrix function %s has "
3640 "invalid value %g."),
3641 y + 1, x + 1, arg_index, props->name,
3642 gsl_matrix_get (a, y, x));
3646 argument_inequality_error (const struct matrix_function_properties *props,
3647 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3648 size_t b_index, double b,
3649 bool greater, bool equal)
3651 struct string s = DS_EMPTY_INITIALIZER;
3653 argument_invalid (props, a_index, a, y, x, &s);
3654 ds_put_cstr (&s, " ");
3655 if (greater && equal)
3656 ds_put_format (&s, _("This argument must be greater than or "
3657 "equal to argument %zu, but that has value %g."),
3659 else if (greater && !equal)
3660 ds_put_format (&s, _("This argument must be greater than argument %zu, "
3661 "but that has value %g."),
3664 ds_put_format (&s, _("This argument must be less than or "
3665 "equal to argument %zu, but that has value %g."),
3668 ds_put_format (&s, _("This argument must be less than argument %zu, "
3669 "but that has value %g."),
3671 msg (ME, ds_cstr (&s));
3676 argument_value_error (const struct matrix_function_properties *props,
3677 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3679 bool greater, bool equal)
3681 struct string s = DS_EMPTY_INITIALIZER;
3682 argument_invalid (props, a_index, a, y, x, &s);
3683 ds_put_cstr (&s, " ");
3684 if (greater && equal)
3685 ds_put_format (&s, _("This argument must be greater than or equal to %g."),
3687 else if (greater && !equal)
3688 ds_put_format (&s, _("This argument must be greater than %g."), b);
3690 ds_put_format (&s, _("This argument must be less than or equal to %g."), b);
3692 ds_put_format (&s, _("This argument must be less than %g."), b);
3693 msg (ME, ds_cstr (&s));
3698 check_constraints (const struct matrix_function_properties *props,
3699 gsl_matrix *args[], size_t n_args)
3701 const char *constraints = props->constraints;
3705 size_t arg_index = SIZE_MAX;
3706 while (*constraints)
3708 if (*constraints >= 'a' && *constraints <= 'd')
3710 arg_index = *constraints++ - 'a';
3711 assert (arg_index < n_args);
3713 else if (*constraints == '[' || *constraints == '(')
3715 assert (arg_index < n_args);
3716 bool open_lower = *constraints++ == '(';
3717 int minimum = parse_constraint_value (&constraints);
3718 assert (*constraints == ',');
3720 int maximum = parse_constraint_value (&constraints);
3721 assert (*constraints == ']' || *constraints == ')');
3722 bool open_upper = *constraints++ == ')';
3724 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3725 if ((open_lower ? *d <= minimum : *d < minimum)
3726 || (open_upper ? *d >= maximum : *d > maximum))
3728 if (!is_scalar (args[arg_index]))
3729 msg (ME, _("Row %zu, column %zu of argument %zu to matrix "
3730 "function %s has value %g, which is outside "
3731 "the valid range %c%d,%d%c."),
3732 y + 1, x + 1, arg_index + 1, props->name, *d,
3733 open_lower ? '(' : '[',
3735 open_upper ? ')' : ']');
3737 msg (ME, _("Argument %zu to matrix function %s has value %g, "
3738 "which is outside the valid range %c%d,%d%c."),
3739 arg_index + 1, props->name, *d,
3740 open_lower ? '(' : '[',
3742 open_upper ? ')' : ']');
3746 else if (*constraints == '>' || *constraints == '<')
3748 bool greater = *constraints++ == '>';
3750 if (*constraints == '=')
3758 if (*constraints >= 'a' && *constraints <= 'd')
3760 size_t a_index = arg_index;
3761 size_t b_index = *constraints - 'a';
3762 assert (a_index < n_args);
3763 assert (b_index < n_args);
3765 /* We only support one of the two arguments being non-scalar.
3766 It's easier to support only the first one being non-scalar, so
3767 flip things around if it's the other way. */
3768 if (!is_scalar (args[b_index]))
3770 assert (is_scalar (args[a_index]));
3771 size_t tmp_index = a_index;
3773 b_index = tmp_index;
3778 double b = to_scalar (args[b_index]);
3779 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
3781 ? (equal ? !(*a >= b) : !(*a > b))
3782 : (equal ? !(*a <= b) : !(*a < b)))
3784 argument_inequality_error (props,
3785 a_index + 1, args[a_index], y, x,
3793 int comparand = parse_constraint_value (&constraints);
3795 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3797 ? (equal ? !(*d >= comparand) : !(*d > comparand))
3798 : (equal ? !(*d <= comparand) : !(*d < comparand)))
3800 argument_value_error (props,
3801 arg_index + 1, args[arg_index], y, x,
3810 assert (*constraints == ' ');
3812 arg_index = SIZE_MAX;
3819 matrix_expr_evaluate_d_none (
3820 const struct matrix_function_properties *props UNUSED,
3821 gsl_matrix *subs[] UNUSED, size_t n_subs,
3822 matrix_proto_d_none *f)
3824 assert (n_subs == 0);
3826 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3827 gsl_matrix_set (m, 0, 0, f ());
3832 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
3833 gsl_matrix *subs[], size_t n_subs,
3834 matrix_proto_d_d *f)
3836 assert (n_subs == 1);
3839 if (!to_scalar_args (props->name, subs, n_subs, &d))
3842 if (!check_constraints (props, subs, n_subs))
3845 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3846 gsl_matrix_set (m, 0, 0, f (d));
3851 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
3852 gsl_matrix *subs[], size_t n_subs,
3853 matrix_proto_d_dd *f)
3855 assert (n_subs == 2);
3858 if (!to_scalar_args (props->name, subs, n_subs, d))
3861 if (!check_constraints (props, subs, n_subs))
3864 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3865 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
3870 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
3871 gsl_matrix *subs[], size_t n_subs,
3872 matrix_proto_d_ddd *f)
3874 assert (n_subs == 3);
3877 if (!to_scalar_args (props->name, subs, n_subs, d))
3880 if (!check_constraints (props, subs, n_subs))
3883 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3884 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
3889 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
3890 gsl_matrix *subs[], size_t n_subs,
3891 matrix_proto_m_d *f)
3893 assert (n_subs == 1);
3896 return to_scalar_args (props->name, subs, n_subs, &d) ? f (d) : NULL;
3900 matrix_expr_evaluate_m_dd (const struct matrix_function_properties *props,
3901 gsl_matrix *subs[], size_t n_subs,
3902 matrix_proto_m_dd *f)
3904 assert (n_subs == 2);
3907 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
3911 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
3912 gsl_matrix *subs[], size_t n_subs,
3913 matrix_proto_m_ddd *f)
3915 assert (n_subs == 3);
3918 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
3922 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED,
3923 gsl_matrix *subs[], size_t n_subs,
3924 matrix_proto_m_m *f)
3926 assert (n_subs == 1);
3931 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
3932 gsl_matrix *subs[], size_t n_subs,
3933 matrix_proto_m_e *f)
3935 assert (n_subs == 1);
3937 if (!check_constraints (props, subs, n_subs))
3940 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3946 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED,
3947 gsl_matrix *subs[], size_t n_subs,
3948 matrix_proto_m_md *f)
3950 assert (n_subs == 2);
3951 if (!check_scalar_arg (props->name, subs, 1))
3953 return f (subs[0], to_scalar (subs[1]));
3957 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
3958 gsl_matrix *subs[], size_t n_subs,
3959 matrix_proto_m_ed *f)
3961 assert (n_subs == 2);
3962 if (!check_scalar_arg (props->name, subs, 1))
3965 if (!check_constraints (props, subs, n_subs))
3968 double b = to_scalar (subs[1]);
3969 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3975 matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props UNUSED,
3976 gsl_matrix *subs[], size_t n_subs,
3977 matrix_proto_m_mdd *f)
3979 assert (n_subs == 3);
3980 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
3982 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
3986 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
3987 gsl_matrix *subs[], size_t n_subs,
3988 matrix_proto_m_edd *f)
3990 assert (n_subs == 3);
3991 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
3994 if (!check_constraints (props, subs, n_subs))
3997 double b = to_scalar (subs[1]);
3998 double c = to_scalar (subs[2]);
3999 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4005 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4006 gsl_matrix *subs[], size_t n_subs,
4007 matrix_proto_m_eddd *f)
4009 assert (n_subs == 4);
4010 for (size_t i = 1; i < 4; i++)
4011 if (!check_scalar_arg (props->name, subs, i))
4014 if (!check_constraints (props, subs, n_subs))
4017 double b = to_scalar (subs[1]);
4018 double c = to_scalar (subs[2]);
4019 double d = to_scalar (subs[3]);
4020 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4021 *a = f (*a, b, c, d);
4026 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4027 gsl_matrix *subs[], size_t n_subs,
4028 matrix_proto_m_eed *f)
4030 assert (n_subs == 3);
4031 if (!check_scalar_arg (props->name, subs, 2))
4034 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4035 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4037 msg (ME, _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4038 "%zu×%zu, but %s requires these arguments either to have "
4039 "the same dimensions or for one of them to be a scalar."),
4041 subs[0]->size1, subs[0]->size2,
4042 subs[1]->size1, subs[1]->size2,
4047 if (!check_constraints (props, subs, n_subs))
4050 double c = to_scalar (subs[2]);
4052 if (is_scalar (subs[0]))
4054 double a = to_scalar (subs[0]);
4055 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4061 double b = to_scalar (subs[1]);
4062 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4069 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED,
4070 gsl_matrix *subs[], size_t n_subs,
4071 matrix_proto_m_mm *f)
4073 assert (n_subs == 2);
4074 return f (subs[0], subs[1]);
4078 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4079 gsl_matrix *subs[], size_t n_subs,
4080 matrix_proto_m_v *f)
4082 assert (n_subs == 1);
4083 if (!check_vector_arg (props->name, subs, 0))
4085 gsl_vector v = to_vector (subs[0]);
4090 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED,
4091 gsl_matrix *subs[], size_t n_subs,
4092 matrix_proto_d_m *f)
4094 assert (n_subs == 1);
4096 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4097 gsl_matrix_set (m, 0, 0, f (subs[0]));
4102 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED,
4103 gsl_matrix *subs[], size_t n_subs,
4104 matrix_proto_m_any *f)
4106 return f (subs, n_subs);
4110 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props,
4111 gsl_matrix *subs[], size_t n_subs,
4112 matrix_proto_IDENT *f)
4114 assert (n_subs <= 2);
4117 if (!to_scalar_args (props->name, subs, n_subs, d))
4120 return f (d[0], d[n_subs - 1]);
4124 matrix_expr_evaluate (const struct matrix_expr *e)
4126 if (e->op == MOP_NUMBER)
4128 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4129 gsl_matrix_set (m, 0, 0, e->number);
4132 else if (e->op == MOP_VARIABLE)
4134 const gsl_matrix *src = e->variable->value;
4137 msg (SE, _("Uninitialized variable %s used in expression."),
4142 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4143 gsl_matrix_memcpy (dst, src);
4146 else if (e->op == MOP_EOF)
4148 struct dfm_reader *reader = read_file_open (e->eof);
4149 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4150 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4154 enum { N_LOCAL = 3 };
4155 gsl_matrix *local_subs[N_LOCAL];
4156 gsl_matrix **subs = (e->n_subs < N_LOCAL
4158 : xmalloc (e->n_subs * sizeof *subs));
4160 for (size_t i = 0; i < e->n_subs; i++)
4162 subs[i] = matrix_expr_evaluate (e->subs[i]);
4165 for (size_t j = 0; j < i; j++)
4166 gsl_matrix_free (subs[j]);
4167 if (subs != local_subs)
4173 gsl_matrix *result = NULL;
4176 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4177 case MOP_F_##ENUM: \
4179 static const struct matrix_function_properties props = { \
4181 .constraints = CONSTRAINTS, \
4183 result = matrix_expr_evaluate_##PROTO (&props, subs, e->n_subs, \
4184 matrix_eval_##ENUM); \
4191 gsl_matrix_scale (subs[0], -1.0);
4209 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
4213 result = matrix_expr_evaluate_not (subs[0]);
4217 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
4221 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
4225 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
4229 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
4232 case MOP_PASTE_HORZ:
4233 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
4236 case MOP_PASTE_VERT:
4237 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
4241 result = gsl_matrix_alloc (0, 0);
4245 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
4249 result = matrix_expr_evaluate_vec_all (subs[0]);
4253 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
4257 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
4261 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
4270 for (size_t i = 0; i < e->n_subs; i++)
4271 if (subs[i] != result)
4272 gsl_matrix_free (subs[i]);
4273 if (subs != local_subs)
4279 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4282 gsl_matrix *m = matrix_expr_evaluate (e);
4288 msg (SE, _("Expression for %s must evaluate to scalar, "
4289 "not a %zu×%zu matrix."),
4290 context, m->size1, m->size2);
4291 gsl_matrix_free (m);
4296 gsl_matrix_free (m);
4301 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4305 if (!matrix_expr_evaluate_scalar (e, context, &d))
4309 if (d < LONG_MIN || d > LONG_MAX)
4311 msg (SE, _("Expression for %s is outside the integer range."), context);
4318 struct matrix_lvalue
4320 struct matrix_var *var;
4321 struct matrix_expr *indexes[2];
4326 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4330 for (size_t i = 0; i < lvalue->n_indexes; i++)
4331 matrix_expr_destroy (lvalue->indexes[i]);
4336 static struct matrix_lvalue *
4337 matrix_lvalue_parse (struct matrix_state *s)
4339 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4340 if (!lex_force_id (s->lexer))
4342 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4343 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4347 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
4352 lex_get_n (s->lexer, 2);
4354 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
4356 lvalue->n_indexes++;
4358 if (lex_match (s->lexer, T_COMMA))
4360 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
4362 lvalue->n_indexes++;
4364 if (!lex_force_match (s->lexer, T_RPAREN))
4370 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4376 matrix_lvalue_destroy (lvalue);
4381 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4382 enum index_type index_type, size_t other_size,
4383 struct index_vector *iv)
4388 m = matrix_expr_evaluate (e);
4395 bool ok = matrix_normalize_index_vector (m, size, index_type,
4397 gsl_matrix_free (m);
4402 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4403 struct index_vector *iv, gsl_matrix *sm)
4405 gsl_vector dv = to_vector (lvalue->var->value);
4407 /* Convert source matrix 'sm' to source vector 'sv'. */
4408 if (!is_vector (sm))
4410 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
4411 sm->size1, sm->size2);
4414 gsl_vector sv = to_vector (sm);
4415 if (iv->n != sv.size)
4417 msg (SE, _("Can't assign %zu-element vector "
4418 "to %zu-element subvector."), sv.size, iv->n);
4422 /* Assign elements. */
4423 for (size_t x = 0; x < iv->n; x++)
4424 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4429 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4430 struct index_vector *iv0,
4431 struct index_vector *iv1, gsl_matrix *sm)
4433 gsl_matrix *dm = lvalue->var->value;
4435 /* Convert source matrix 'sm' to source vector 'sv'. */
4436 if (iv0->n != sm->size1)
4438 msg (SE, _("Row index vector for assignment to %s has %zu elements "
4439 "but source matrix has %zu rows."),
4440 lvalue->var->name, iv0->n, sm->size1);
4443 if (iv1->n != sm->size2)
4445 msg (SE, _("Column index vector for assignment to %s has %zu elements "
4446 "but source matrix has %zu columns."),
4447 lvalue->var->name, iv1->n, sm->size2);
4451 /* Assign elements. */
4452 for (size_t y = 0; y < iv0->n; y++)
4454 size_t f0 = iv0->indexes[y];
4455 for (size_t x = 0; x < iv1->n; x++)
4457 size_t f1 = iv1->indexes[x];
4458 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4464 /* Takes ownership of M. */
4466 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4467 struct index_vector *iv0, struct index_vector *iv1,
4470 if (!lvalue->n_indexes)
4472 gsl_matrix_free (lvalue->var->value);
4473 lvalue->var->value = sm;
4478 bool ok = (lvalue->n_indexes == 1
4479 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
4480 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
4481 free (iv0->indexes);
4482 free (iv1->indexes);
4483 gsl_matrix_free (sm);
4489 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4490 struct index_vector *iv0,
4491 struct index_vector *iv1)
4493 *iv0 = INDEX_VECTOR_INIT;
4494 *iv1 = INDEX_VECTOR_INIT;
4495 if (!lvalue->n_indexes)
4498 /* Validate destination matrix exists and has the right shape. */
4499 gsl_matrix *dm = lvalue->var->value;
4502 msg (SE, _("Undefined variable %s."), lvalue->var->name);
4505 else if (dm->size1 == 0 || dm->size2 == 0)
4507 msg (SE, _("Cannot index %zu×%zu matrix."),
4508 dm->size1, dm->size2);
4511 else if (lvalue->n_indexes == 1)
4513 if (!is_vector (dm))
4515 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
4516 dm->size1, dm->size2, lvalue->var->name);
4519 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4520 MAX (dm->size1, dm->size2),
4525 assert (lvalue->n_indexes == 2);
4526 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4527 IV_ROW, dm->size2, iv0))
4530 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4531 IV_COLUMN, dm->size1, iv1))
4533 free (iv0->indexes);
4540 /* Takes ownership of M. */
4542 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
4544 struct index_vector iv0, iv1;
4545 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
4547 gsl_matrix_free (sm);
4551 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
4555 /* Matrix command. */
4559 struct matrix_cmd **commands;
4563 static void matrix_cmds_uninit (struct matrix_cmds *);
4567 enum matrix_cmd_type
4590 struct compute_command
4592 struct matrix_lvalue *lvalue;
4593 struct matrix_expr *rvalue;
4597 struct print_command
4599 struct matrix_expr *expression;
4600 bool use_default_format;
4601 struct fmt_spec format;
4603 int space; /* -1 means NEWPAGE. */
4607 struct string_array literals; /* CLABELS/RLABELS. */
4608 struct matrix_expr *expr; /* CNAMES/RNAMES. */
4614 struct do_if_command
4618 struct matrix_expr *condition;
4619 struct matrix_cmds commands;
4629 struct matrix_var *var;
4630 struct matrix_expr *start;
4631 struct matrix_expr *end;
4632 struct matrix_expr *increment;
4634 /* Loop conditions. */
4635 struct matrix_expr *top_condition;
4636 struct matrix_expr *bottom_condition;
4639 struct matrix_cmds commands;
4643 struct display_command
4645 struct matrix_state *state;
4649 struct release_command
4651 struct matrix_var **vars;
4658 struct matrix_expr *expression;
4659 struct save_file *sf;
4665 struct read_file *rf;
4666 struct matrix_lvalue *dst;
4667 struct matrix_expr *size;
4669 enum fmt_type format;
4676 struct write_command
4678 struct write_file *wf;
4679 struct matrix_expr *expression;
4682 /* If this is nonnull, WRITE uses this format.
4684 If this is NULL, WRITE uses free-field format with as many
4685 digits of precision as needed. */
4686 struct fmt_spec *format;
4695 struct matrix_lvalue *dst;
4696 struct dataset *dataset;
4697 struct file_handle *file;
4699 struct var_syntax *vars;
4701 struct matrix_var *names;
4703 /* Treatment of missing values. */
4708 MGET_ERROR, /* Flag error on command. */
4709 MGET_ACCEPT, /* Accept without change, user-missing only. */
4710 MGET_OMIT, /* Drop this case. */
4711 MGET_RECODE /* Recode to 'substitute'. */
4714 double substitute; /* MGET_RECODE only. */
4720 struct msave_command
4722 struct msave_common *common;
4723 struct matrix_expr *expr;
4724 const char *rowtype;
4725 const struct matrix_expr *factors;
4726 const struct matrix_expr *splits;
4732 struct matrix_state *state;
4733 struct file_handle *file;
4735 struct stringi_set rowtypes;
4739 struct eigen_command
4741 struct matrix_expr *expr;
4742 struct matrix_var *evec;
4743 struct matrix_var *eval;
4747 struct setdiag_command
4749 struct matrix_var *dst;
4750 struct matrix_expr *expr;
4756 struct matrix_expr *expr;
4757 struct matrix_var *u;
4758 struct matrix_var *s;
4759 struct matrix_var *v;
4765 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
4766 static bool matrix_cmd_execute (struct matrix_cmd *);
4767 static void matrix_cmd_destroy (struct matrix_cmd *);
4770 static struct matrix_cmd *
4771 matrix_parse_compute (struct matrix_state *s)
4773 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4774 *cmd = (struct matrix_cmd) {
4775 .type = MCMD_COMPUTE,
4776 .compute = { .lvalue = NULL },
4779 struct compute_command *compute = &cmd->compute;
4780 compute->lvalue = matrix_lvalue_parse (s);
4781 if (!compute->lvalue)
4784 if (!lex_force_match (s->lexer, T_EQUALS))
4787 compute->rvalue = matrix_parse_expr (s);
4788 if (!compute->rvalue)
4794 matrix_cmd_destroy (cmd);
4799 matrix_cmd_execute_compute (struct compute_command *compute)
4801 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
4805 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
4809 print_labels_destroy (struct print_labels *labels)
4813 string_array_destroy (&labels->literals);
4814 matrix_expr_destroy (labels->expr);
4820 parse_literal_print_labels (struct matrix_state *s,
4821 struct print_labels **labelsp)
4823 lex_match (s->lexer, T_EQUALS);
4824 print_labels_destroy (*labelsp);
4825 *labelsp = xzalloc (sizeof **labelsp);
4826 while (lex_token (s->lexer) != T_SLASH
4827 && lex_token (s->lexer) != T_ENDCMD
4828 && lex_token (s->lexer) != T_STOP)
4830 struct string label = DS_EMPTY_INITIALIZER;
4831 while (lex_token (s->lexer) != T_COMMA
4832 && lex_token (s->lexer) != T_SLASH
4833 && lex_token (s->lexer) != T_ENDCMD
4834 && lex_token (s->lexer) != T_STOP)
4836 if (!ds_is_empty (&label))
4837 ds_put_byte (&label, ' ');
4839 if (lex_token (s->lexer) == T_STRING)
4840 ds_put_cstr (&label, lex_tokcstr (s->lexer));
4843 char *rep = lex_next_representation (s->lexer, 0, 0);
4844 ds_put_cstr (&label, rep);
4849 string_array_append_nocopy (&(*labelsp)->literals,
4850 ds_steal_cstr (&label));
4852 if (!lex_match (s->lexer, T_COMMA))
4858 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
4860 lex_match (s->lexer, T_EQUALS);
4861 struct matrix_expr *e = matrix_parse_exp (s);
4865 print_labels_destroy (*labelsp);
4866 *labelsp = xzalloc (sizeof **labelsp);
4867 (*labelsp)->expr = e;
4871 static struct matrix_cmd *
4872 matrix_parse_print (struct matrix_state *s)
4874 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4875 *cmd = (struct matrix_cmd) {
4878 .use_default_format = true,
4882 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
4885 for (size_t i = 0; ; i++)
4887 enum token_type t = lex_next_token (s->lexer, i);
4888 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
4890 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
4892 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
4895 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
4900 cmd->print.expression = matrix_parse_exp (s);
4901 if (!cmd->print.expression)
4905 while (lex_match (s->lexer, T_SLASH))
4907 if (lex_match_id (s->lexer, "FORMAT"))
4909 lex_match (s->lexer, T_EQUALS);
4910 if (!parse_format_specifier (s->lexer, &cmd->print.format))
4912 cmd->print.use_default_format = false;
4914 else if (lex_match_id (s->lexer, "TITLE"))
4916 lex_match (s->lexer, T_EQUALS);
4917 if (!lex_force_string (s->lexer))
4919 free (cmd->print.title);
4920 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
4923 else if (lex_match_id (s->lexer, "SPACE"))
4925 lex_match (s->lexer, T_EQUALS);
4926 if (lex_match_id (s->lexer, "NEWPAGE"))
4927 cmd->print.space = -1;
4928 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
4930 cmd->print.space = lex_integer (s->lexer);
4936 else if (lex_match_id (s->lexer, "RLABELS"))
4937 parse_literal_print_labels (s, &cmd->print.rlabels);
4938 else if (lex_match_id (s->lexer, "CLABELS"))
4939 parse_literal_print_labels (s, &cmd->print.clabels);
4940 else if (lex_match_id (s->lexer, "RNAMES"))
4942 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
4945 else if (lex_match_id (s->lexer, "CNAMES"))
4947 if (!parse_expr_print_labels (s, &cmd->print.clabels))
4952 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
4953 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
4961 matrix_cmd_destroy (cmd);
4966 matrix_is_integer (const gsl_matrix *m)
4968 for (size_t y = 0; y < m->size1; y++)
4969 for (size_t x = 0; x < m->size2; x++)
4971 double d = gsl_matrix_get (m, y, x);
4979 matrix_max_magnitude (const gsl_matrix *m)
4982 for (size_t y = 0; y < m->size1; y++)
4983 for (size_t x = 0; x < m->size2; x++)
4985 double d = fabs (gsl_matrix_get (m, y, x));
4993 format_fits (struct fmt_spec format, double x)
4995 char *s = data_out (&(union value) { .f = x }, NULL,
4996 &format, settings_get_fmt_settings ());
4997 bool fits = *s != '*' && !strchr (s, 'E');
5002 static struct fmt_spec
5003 default_format (const gsl_matrix *m, int *log_scale)
5007 double max = matrix_max_magnitude (m);
5009 if (matrix_is_integer (m))
5010 for (int w = 1; w <= 12; w++)
5012 struct fmt_spec format = { .type = FMT_F, .w = w };
5013 if (format_fits (format, -max))
5017 if (max >= 1e9 || max <= 1e-4)
5020 snprintf (s, sizeof s, "%.1e", max);
5022 const char *e = strchr (s, 'e');
5024 *log_scale = atoi (e + 1);
5027 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5031 trimmed_string (double d)
5033 struct substring s = ss_buffer ((char *) &d, sizeof d);
5034 ss_rtrim (&s, ss_cstr (" "));
5035 return ss_xstrdup (s);
5038 static struct string_array *
5039 print_labels_get (const struct print_labels *labels, size_t n,
5040 const char *prefix, bool truncate)
5045 struct string_array *out = xzalloc (sizeof *out);
5046 if (labels->literals.n)
5047 string_array_clone (out, &labels->literals);
5048 else if (labels->expr)
5050 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5051 if (m && is_vector (m))
5053 gsl_vector v = to_vector (m);
5054 for (size_t i = 0; i < v.size; i++)
5055 string_array_append_nocopy (out, trimmed_string (
5056 gsl_vector_get (&v, i)));
5058 gsl_matrix_free (m);
5064 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5067 string_array_append (out, "");
5070 string_array_delete (out, out->n - 1);
5073 for (size_t i = 0; i < out->n; i++)
5075 char *s = out->strings[i];
5076 s[strnlen (s, 8)] = '\0';
5083 matrix_cmd_print_space (int space)
5086 output_item_submit (page_break_item_create ());
5087 for (int i = 0; i < space; i++)
5088 output_log ("%s", "");
5092 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
5093 struct fmt_spec format, int log_scale)
5095 matrix_cmd_print_space (print->space);
5097 output_log ("%s", print->title);
5099 output_log (" 10 ** %d X", log_scale);
5101 struct string_array *clabels = print_labels_get (print->clabels,
5102 m->size2, "col", true);
5103 if (clabels && format.w < 8)
5105 struct string_array *rlabels = print_labels_get (print->rlabels,
5106 m->size1, "row", true);
5110 struct string line = DS_EMPTY_INITIALIZER;
5112 ds_put_byte_multiple (&line, ' ', 8);
5113 for (size_t x = 0; x < m->size2; x++)
5114 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5115 output_log_nocopy (ds_steal_cstr (&line));
5118 double scale = pow (10.0, log_scale);
5119 bool numeric = fmt_is_numeric (format.type);
5120 for (size_t y = 0; y < m->size1; y++)
5122 struct string line = DS_EMPTY_INITIALIZER;
5124 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5126 for (size_t x = 0; x < m->size2; x++)
5128 double f = gsl_matrix_get (m, y, x);
5130 ? data_out (&(union value) { .f = f / scale}, NULL,
5131 &format, settings_get_fmt_settings ())
5132 : trimmed_string (f));
5133 ds_put_format (&line, " %s", s);
5136 output_log_nocopy (ds_steal_cstr (&line));
5139 string_array_destroy (rlabels);
5141 string_array_destroy (clabels);
5146 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5147 const struct print_labels *print_labels, size_t n,
5148 const char *name, const char *prefix)
5150 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5152 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5153 for (size_t i = 0; i < n; i++)
5154 pivot_category_create_leaf (
5156 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5157 : pivot_value_new_integer (i + 1)));
5159 d->hide_all_labels = true;
5160 string_array_destroy (labels);
5165 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
5166 struct fmt_spec format, int log_scale)
5168 struct pivot_table *table = pivot_table_create__ (
5169 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5172 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5174 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5175 N_("Columns"), "col");
5177 struct pivot_footnote *footnote = NULL;
5180 char *s = xasprintf ("× 10**%d", log_scale);
5181 footnote = pivot_table_create_footnote (
5182 table, pivot_value_new_user_text_nocopy (s));
5185 double scale = pow (10.0, log_scale);
5186 bool numeric = fmt_is_numeric (format.type);
5187 for (size_t y = 0; y < m->size1; y++)
5188 for (size_t x = 0; x < m->size2; x++)
5190 double f = gsl_matrix_get (m, y, x);
5191 struct pivot_value *v;
5194 v = pivot_value_new_number (f / scale);
5195 v->numeric.format = format;
5198 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5200 pivot_value_add_footnote (v, footnote);
5201 pivot_table_put2 (table, y, x, v);
5204 pivot_table_submit (table);
5208 matrix_cmd_execute_print (const struct print_command *print)
5210 if (print->expression)
5212 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5217 struct fmt_spec format = (print->use_default_format
5218 ? default_format (m, &log_scale)
5221 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5222 matrix_cmd_print_text (print, m, format, log_scale);
5224 matrix_cmd_print_tables (print, m, format, log_scale);
5226 gsl_matrix_free (m);
5230 matrix_cmd_print_space (print->space);
5232 output_log ("%s", print->title);
5239 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
5240 const char *command_name,
5241 const char *stop1, const char *stop2)
5243 lex_end_of_command (s->lexer);
5244 lex_discard_rest_of_command (s->lexer);
5246 size_t allocated = 0;
5249 while (lex_token (s->lexer) == T_ENDCMD)
5252 if (lex_at_phrase (s->lexer, stop1)
5253 || (stop2 && lex_at_phrase (s->lexer, stop2)))
5256 if (lex_at_phrase (s->lexer, "END MATRIX"))
5258 msg (SE, _("Premature END MATRIX within %s."), command_name);
5262 struct matrix_cmd *cmd = matrix_parse_command (s);
5266 if (c->n >= allocated)
5267 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
5268 c->commands[c->n++] = cmd;
5273 matrix_parse_do_if_clause (struct matrix_state *s,
5274 struct do_if_command *ifc,
5275 bool parse_condition,
5276 size_t *allocated_clauses)
5278 if (ifc->n_clauses >= *allocated_clauses)
5279 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5280 sizeof *ifc->clauses);
5281 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5282 *c = (struct do_if_clause) { .condition = NULL };
5284 if (parse_condition)
5286 c->condition = matrix_parse_expr (s);
5291 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
5294 static struct matrix_cmd *
5295 matrix_parse_do_if (struct matrix_state *s)
5297 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5298 *cmd = (struct matrix_cmd) {
5300 .do_if = { .n_clauses = 0 }
5303 size_t allocated_clauses = 0;
5306 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
5309 while (lex_match_phrase (s->lexer, "ELSE IF"));
5311 if (lex_match_id (s->lexer, "ELSE")
5312 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
5315 if (!lex_match_phrase (s->lexer, "END IF"))
5320 matrix_cmd_destroy (cmd);
5325 matrix_cmd_execute_do_if (struct do_if_command *cmd)
5327 for (size_t i = 0; i < cmd->n_clauses; i++)
5329 struct do_if_clause *c = &cmd->clauses[i];
5333 if (!matrix_expr_evaluate_scalar (c->condition,
5334 i ? "ELSE IF" : "DO IF",
5339 for (size_t j = 0; j < c->commands.n; j++)
5340 if (!matrix_cmd_execute (c->commands.commands[j]))
5347 static struct matrix_cmd *
5348 matrix_parse_loop (struct matrix_state *s)
5350 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5351 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5353 struct loop_command *loop = &cmd->loop;
5354 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5356 struct substring name = lex_tokss (s->lexer);
5357 loop->var = matrix_var_lookup (s, name);
5359 loop->var = matrix_var_create (s, name);
5364 loop->start = matrix_parse_expr (s);
5365 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5368 loop->end = matrix_parse_expr (s);
5372 if (lex_match (s->lexer, T_BY))
5374 loop->increment = matrix_parse_expr (s);
5375 if (!loop->increment)
5380 if (lex_match_id (s->lexer, "IF"))
5382 loop->top_condition = matrix_parse_expr (s);
5383 if (!loop->top_condition)
5387 bool was_in_loop = s->in_loop;
5389 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
5391 s->in_loop = was_in_loop;
5395 if (!lex_match_phrase (s->lexer, "END LOOP"))
5398 if (lex_match_id (s->lexer, "IF"))
5400 loop->bottom_condition = matrix_parse_expr (s);
5401 if (!loop->bottom_condition)
5408 matrix_cmd_destroy (cmd);
5413 matrix_cmd_execute_loop (struct loop_command *cmd)
5417 long int increment = 1;
5420 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5421 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5423 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5427 if (increment > 0 ? value > end
5428 : increment < 0 ? value < end
5433 int mxloops = settings_get_mxloops ();
5434 for (int i = 0; i < mxloops; i++)
5439 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5441 gsl_matrix_free (cmd->var->value);
5442 cmd->var->value = NULL;
5444 if (!cmd->var->value)
5445 cmd->var->value = gsl_matrix_alloc (1, 1);
5447 gsl_matrix_set (cmd->var->value, 0, 0, value);
5450 if (cmd->top_condition)
5453 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5458 for (size_t j = 0; j < cmd->commands.n; j++)
5459 if (!matrix_cmd_execute (cmd->commands.commands[j]))
5462 if (cmd->bottom_condition)
5465 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5466 "END LOOP IF", &d) || d > 0)
5473 if (increment > 0 ? value > end : value < end)
5479 static struct matrix_cmd *
5480 matrix_parse_break (struct matrix_state *s)
5484 msg (SE, _("BREAK not inside LOOP."));
5488 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5489 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
5493 static struct matrix_cmd *
5494 matrix_parse_display (struct matrix_state *s)
5496 while (lex_token (s->lexer) != T_ENDCMD)
5498 if (!lex_match_id (s->lexer, "DICTIONARY")
5499 && !lex_match_id (s->lexer, "STATUS"))
5501 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5506 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5507 *cmd = (struct matrix_cmd) { .type = MCMD_DISPLAY, .display = { s } };
5512 compare_matrix_var_pointers (const void *a_, const void *b_)
5514 const struct matrix_var *const *ap = a_;
5515 const struct matrix_var *const *bp = b_;
5516 const struct matrix_var *a = *ap;
5517 const struct matrix_var *b = *bp;
5518 return strcmp (a->name, b->name);
5522 matrix_cmd_execute_display (struct display_command *cmd)
5524 const struct matrix_state *s = cmd->state;
5526 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5527 struct pivot_dimension *attr_dimension
5528 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5529 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5530 N_("Rows"), N_("Columns"));
5531 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5533 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5536 const struct matrix_var *var;
5537 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5539 vars[n_vars++] = var;
5540 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5542 struct pivot_dimension *rows = pivot_dimension_create (
5543 table, PIVOT_AXIS_ROW, N_("Variable"));
5544 for (size_t i = 0; i < n_vars; i++)
5546 const struct matrix_var *var = vars[i];
5547 pivot_category_create_leaf (
5548 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
5550 size_t r = var->value->size1;
5551 size_t c = var->value->size2;
5552 double values[] = { r, c, r * c * sizeof (double) / 1024 };
5553 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
5554 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
5557 pivot_table_submit (table);
5560 static struct matrix_cmd *
5561 matrix_parse_release (struct matrix_state *s)
5563 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5564 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
5566 size_t allocated_vars = 0;
5567 while (lex_token (s->lexer) == T_ID)
5569 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
5572 if (cmd->release.n_vars >= allocated_vars)
5573 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
5574 sizeof *cmd->release.vars);
5575 cmd->release.vars[cmd->release.n_vars++] = var;
5578 lex_error (s->lexer, _("Variable name expected."));
5581 if (!lex_match (s->lexer, T_COMMA))
5589 matrix_cmd_execute_release (struct release_command *cmd)
5591 for (size_t i = 0; i < cmd->n_vars; i++)
5593 struct matrix_var *var = cmd->vars[i];
5594 gsl_matrix_free (var->value);
5599 static struct save_file *
5600 save_file_create (struct matrix_state *s, struct file_handle *fh,
5601 struct string_array *variables,
5602 struct matrix_expr *names,
5603 struct stringi_set *strings)
5605 for (size_t i = 0; i < s->n_save_files; i++)
5607 struct save_file *sf = s->save_files[i];
5608 if (fh_equal (sf->file, fh))
5612 string_array_destroy (variables);
5613 matrix_expr_destroy (names);
5614 stringi_set_destroy (strings);
5620 struct save_file *sf = xmalloc (sizeof *sf);
5621 *sf = (struct save_file) {
5623 .dataset = s->dataset,
5624 .variables = *variables,
5626 .strings = STRINGI_SET_INITIALIZER (sf->strings),
5629 stringi_set_swap (&sf->strings, strings);
5630 stringi_set_destroy (strings);
5632 s->save_files = xrealloc (s->save_files,
5633 (s->n_save_files + 1) * sizeof *s->save_files);
5634 s->save_files[s->n_save_files++] = sf;
5638 static struct casewriter *
5639 save_file_open (struct save_file *sf, gsl_matrix *m)
5641 if (sf->writer || sf->error)
5645 size_t n_variables = caseproto_get_n_widths (
5646 casewriter_get_proto (sf->writer));
5647 if (m->size2 != n_variables)
5649 msg (ME, _("The first SAVE to %s within this matrix program "
5650 "had %zu columns, so a %zu×%zu matrix cannot be "
5652 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
5653 n_variables, m->size1, m->size2);
5661 struct dictionary *dict = dict_create (get_default_encoding ());
5663 /* Fill 'names' with user-specified names if there were any, either from
5664 sf->variables or sf->names. */
5665 struct string_array names = { .n = 0 };
5666 if (sf->variables.n)
5667 string_array_clone (&names, &sf->variables);
5670 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
5671 if (nm && is_vector (nm))
5673 gsl_vector nv = to_vector (nm);
5674 for (size_t i = 0; i < nv.size; i++)
5676 char *name = trimmed_string (gsl_vector_get (&nv, i));
5677 if (dict_id_is_valid (dict, name, true))
5678 string_array_append_nocopy (&names, name);
5683 gsl_matrix_free (nm);
5686 struct stringi_set strings;
5687 stringi_set_clone (&strings, &sf->strings);
5689 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
5694 name = names.strings[i];
5697 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
5701 int width = stringi_set_delete (&strings, name) ? 8 : 0;
5702 struct variable *var = dict_create_var (dict, name, width);
5705 msg (ME, _("Duplicate variable name %s in SAVE statement."),
5711 if (!stringi_set_is_empty (&strings))
5713 size_t count = stringi_set_count (&strings);
5714 const char *example = stringi_set_node_get_string (
5715 stringi_set_first (&strings));
5718 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
5719 "variable %s."), example);
5721 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
5722 "unknown variable: %s.",
5723 "The SAVE command STRINGS subcommand specifies %zu "
5724 "unknown variables, including %s.",
5729 stringi_set_destroy (&strings);
5730 string_array_destroy (&names);
5739 if (sf->file == fh_inline_file ())
5740 sf->writer = autopaging_writer_create (dict_get_proto (dict));
5742 sf->writer = any_writer_open (sf->file, dict);
5754 save_file_destroy (struct save_file *sf)
5758 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
5760 dataset_set_dict (sf->dataset, sf->dict);
5761 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
5765 casewriter_destroy (sf->writer);
5766 dict_unref (sf->dict);
5768 fh_unref (sf->file);
5769 string_array_destroy (&sf->variables);
5770 matrix_expr_destroy (sf->names);
5771 stringi_set_destroy (&sf->strings);
5776 static struct matrix_cmd *
5777 matrix_parse_save (struct matrix_state *s)
5779 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5780 *cmd = (struct matrix_cmd) {
5782 .save = { .expression = NULL },
5785 struct file_handle *fh = NULL;
5786 struct save_command *save = &cmd->save;
5788 struct string_array variables = STRING_ARRAY_INITIALIZER;
5789 struct matrix_expr *names = NULL;
5790 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
5792 save->expression = matrix_parse_exp (s);
5793 if (!save->expression)
5796 while (lex_match (s->lexer, T_SLASH))
5798 if (lex_match_id (s->lexer, "OUTFILE"))
5800 lex_match (s->lexer, T_EQUALS);
5803 fh = (lex_match (s->lexer, T_ASTERISK)
5805 : fh_parse (s->lexer, FH_REF_FILE, s->session));
5809 else if (lex_match_id (s->lexer, "VARIABLES"))
5811 lex_match (s->lexer, T_EQUALS);
5815 struct dictionary *d = dict_create (get_default_encoding ());
5816 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
5817 PV_NO_SCRATCH | PV_NO_DUPLICATE);
5822 string_array_clear (&variables);
5823 variables = (struct string_array) {
5829 else if (lex_match_id (s->lexer, "NAMES"))
5831 lex_match (s->lexer, T_EQUALS);
5832 matrix_expr_destroy (names);
5833 names = matrix_parse_exp (s);
5837 else if (lex_match_id (s->lexer, "STRINGS"))
5839 lex_match (s->lexer, T_EQUALS);
5840 while (lex_token (s->lexer) == T_ID)
5842 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
5844 if (!lex_match (s->lexer, T_COMMA))
5850 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
5858 if (s->prev_save_file)
5859 fh = fh_ref (s->prev_save_file);
5862 lex_sbc_missing ("OUTFILE");
5866 fh_unref (s->prev_save_file);
5867 s->prev_save_file = fh_ref (fh);
5869 if (variables.n && names)
5871 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
5872 matrix_expr_destroy (names);
5876 save->sf = save_file_create (s, fh, &variables, names, &strings);
5880 string_array_destroy (&variables);
5881 matrix_expr_destroy (names);
5882 stringi_set_destroy (&strings);
5884 matrix_cmd_destroy (cmd);
5889 matrix_cmd_execute_save (const struct save_command *save)
5891 gsl_matrix *m = matrix_expr_evaluate (save->expression);
5895 struct casewriter *writer = save_file_open (save->sf, m);
5898 gsl_matrix_free (m);
5902 const struct caseproto *proto = casewriter_get_proto (writer);
5903 for (size_t y = 0; y < m->size1; y++)
5905 struct ccase *c = case_create (proto);
5906 for (size_t x = 0; x < m->size2; x++)
5908 double d = gsl_matrix_get (m, y, x);
5909 int width = caseproto_get_width (proto, x);
5910 union value *value = case_data_rw_idx (c, x);
5914 memcpy (value->s, &d, width);
5916 casewriter_write (writer, c);
5918 gsl_matrix_free (m);
5921 static struct matrix_cmd *
5922 matrix_parse_read (struct matrix_state *s)
5924 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5925 *cmd = (struct matrix_cmd) {
5927 .read = { .format = FMT_F },
5930 struct file_handle *fh = NULL;
5931 char *encoding = NULL;
5932 struct read_command *read = &cmd->read;
5933 read->dst = matrix_lvalue_parse (s);
5938 int repetitions = 0;
5939 int record_width = 0;
5940 bool seen_format = false;
5941 while (lex_match (s->lexer, T_SLASH))
5943 if (lex_match_id (s->lexer, "FILE"))
5945 lex_match (s->lexer, T_EQUALS);
5948 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5952 else if (lex_match_id (s->lexer, "ENCODING"))
5954 lex_match (s->lexer, T_EQUALS);
5955 if (!lex_force_string (s->lexer))
5959 encoding = ss_xstrdup (lex_tokss (s->lexer));
5963 else if (lex_match_id (s->lexer, "FIELD"))
5965 lex_match (s->lexer, T_EQUALS);
5967 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5969 read->c1 = lex_integer (s->lexer);
5971 if (!lex_force_match (s->lexer, T_TO)
5972 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
5974 read->c2 = lex_integer (s->lexer) + 1;
5977 record_width = read->c2 - read->c1;
5978 if (lex_match (s->lexer, T_BY))
5980 if (!lex_force_int_range (s->lexer, "BY", 1,
5981 read->c2 - read->c1))
5983 by = lex_integer (s->lexer);
5986 if (record_width % by)
5988 msg (SE, _("BY %d does not evenly divide record width %d."),
5996 else if (lex_match_id (s->lexer, "SIZE"))
5998 lex_match (s->lexer, T_EQUALS);
5999 matrix_expr_destroy (read->size);
6000 read->size = matrix_parse_exp (s);
6004 else if (lex_match_id (s->lexer, "MODE"))
6006 lex_match (s->lexer, T_EQUALS);
6007 if (lex_match_id (s->lexer, "RECTANGULAR"))
6008 read->symmetric = false;
6009 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6010 read->symmetric = true;
6013 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6017 else if (lex_match_id (s->lexer, "REREAD"))
6018 read->reread = true;
6019 else if (lex_match_id (s->lexer, "FORMAT"))
6023 lex_sbc_only_once ("FORMAT");
6028 lex_match (s->lexer, T_EQUALS);
6030 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6033 const char *p = lex_tokcstr (s->lexer);
6034 if (c_isdigit (p[0]))
6036 repetitions = atoi (p);
6037 p += strspn (p, "0123456789");
6038 if (!fmt_from_name (p, &read->format))
6040 lex_error (s->lexer, _("Unknown format %s."), p);
6045 else if (fmt_from_name (p, &read->format))
6049 struct fmt_spec format;
6050 if (!parse_format_specifier (s->lexer, &format))
6052 read->format = format.type;
6058 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6059 "REREAD", "FORMAT");
6066 lex_sbc_missing ("FIELD");
6070 if (!read->dst->n_indexes && !read->size)
6072 msg (SE, _("SIZE is required for reading data into a full matrix "
6073 "(as opposed to a submatrix)."));
6079 if (s->prev_read_file)
6080 fh = fh_ref (s->prev_read_file);
6083 lex_sbc_missing ("FILE");
6087 fh_unref (s->prev_read_file);
6088 s->prev_read_file = fh_ref (fh);
6090 read->rf = read_file_create (s, fh);
6094 free (read->rf->encoding);
6095 read->rf->encoding = encoding;
6099 /* Field width may be specified in multiple ways:
6102 2. The format on FORMAT.
6103 3. The repetition factor on FORMAT.
6105 (2) and (3) are mutually exclusive.
6107 If more than one of these is present, they must agree. If none of them is
6108 present, then free-field format is used.
6110 if (repetitions > record_width)
6112 msg (SE, _("%d repetitions cannot fit in record width %d."),
6113 repetitions, record_width);
6116 int w = (repetitions ? record_width / repetitions
6122 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6123 "which implies field width %d, "
6124 "but BY specifies field width %d."),
6125 repetitions, record_width, w, by);
6127 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6136 matrix_cmd_destroy (cmd);
6142 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6143 struct substring data, size_t y, size_t x,
6144 int first_column, int last_column, char *error)
6146 int line_number = dfm_get_line_number (reader);
6147 struct msg_location *location = xmalloc (sizeof *location);
6148 *location = (struct msg_location) {
6149 .file_name = xstrdup (dfm_get_file_name (reader)),
6150 .first_line = line_number,
6151 .last_line = line_number + 1,
6152 .first_column = first_column,
6153 .last_column = last_column,
6155 struct msg *m = xmalloc (sizeof *m);
6157 .category = MSG_C_DATA,
6158 .severity = MSG_S_WARNING,
6159 .location = location,
6160 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
6161 "for matrix row %zu, column %zu: %s"),
6162 (int) data.length, data.string, fmt_name (format),
6163 y + 1, x + 1, error),
6171 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
6172 gsl_matrix *m, struct substring p, size_t y, size_t x,
6173 const char *line_start)
6175 const char *input_encoding = dfm_reader_get_encoding (reader);
6178 if (fmt_is_numeric (read->format))
6181 error = data_in (p, input_encoding, read->format,
6182 settings_get_fmt_settings (), &v, 0, NULL);
6183 if (!error && v.f == SYSMIS)
6184 error = xstrdup (_("Matrix data may not contain missing value."));
6189 uint8_t s[sizeof (double)];
6190 union value v = { .s = s };
6191 error = data_in (p, input_encoding, read->format,
6192 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6193 memcpy (&f, s, sizeof f);
6198 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6199 int c2 = c1 + ss_utf8_count_columns (p) - 1;
6200 parse_error (reader, read->format, p, y, x, c1, c2, error);
6204 gsl_matrix_set (m, y, x, f);
6205 if (read->symmetric && x != y)
6206 gsl_matrix_set (m, x, y, f);
6211 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
6212 struct substring *line, const char **startp)
6214 if (dfm_eof (reader))
6216 msg (SE, _("Unexpected end of file reading matrix data."));
6219 dfm_expand_tabs (reader);
6220 struct substring record = dfm_get_record (reader);
6221 /* XXX need to recode record into UTF-8 */
6222 *startp = record.string;
6223 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6228 matrix_read (struct read_command *read, struct dfm_reader *reader,
6231 for (size_t y = 0; y < m->size1; y++)
6233 size_t nx = read->symmetric ? y + 1 : m->size2;
6235 struct substring line = ss_empty ();
6236 const char *line_start = line.string;
6237 for (size_t x = 0; x < nx; x++)
6244 ss_ltrim (&line, ss_cstr (" ,"));
6245 if (!ss_is_empty (line))
6247 if (!matrix_read_line (read, reader, &line, &line_start))
6249 dfm_forward_record (reader);
6252 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6256 if (!matrix_read_line (read, reader, &line, &line_start))
6258 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6259 int f = x % fields_per_line;
6260 if (f == fields_per_line - 1)
6261 dfm_forward_record (reader);
6263 p = ss_substr (line, read->w * f, read->w);
6266 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6270 dfm_forward_record (reader);
6273 ss_ltrim (&line, ss_cstr (" ,"));
6274 if (!ss_is_empty (line))
6277 msg (SW, _("Trailing garbage on line \"%.*s\""),
6278 (int) line.length, line.string);
6285 matrix_cmd_execute_read (struct read_command *read)
6287 struct index_vector iv0, iv1;
6288 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6291 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6294 gsl_matrix *m = matrix_expr_evaluate (read->size);
6300 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6301 "not a %zu×%zu matrix."), m->size1, m->size2);
6302 gsl_matrix_free (m);
6308 gsl_vector v = to_vector (m);
6312 d[0] = gsl_vector_get (&v, 0);
6315 else if (v.size == 2)
6317 d[0] = gsl_vector_get (&v, 0);
6318 d[1] = gsl_vector_get (&v, 1);
6322 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6323 "not a %zu×%zu matrix."),
6324 m->size1, m->size2),
6325 gsl_matrix_free (m);
6330 gsl_matrix_free (m);
6332 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6334 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
6335 "are outside valid range."),
6346 if (read->dst->n_indexes)
6348 size_t submatrix_size[2];
6349 if (read->dst->n_indexes == 2)
6351 submatrix_size[0] = iv0.n;
6352 submatrix_size[1] = iv1.n;
6354 else if (read->dst->var->value->size1 == 1)
6356 submatrix_size[0] = 1;
6357 submatrix_size[1] = iv0.n;
6361 submatrix_size[0] = iv0.n;
6362 submatrix_size[1] = 1;
6367 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6369 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
6370 "differ from submatrix dimensions %zu×%zu."),
6372 submatrix_size[0], submatrix_size[1]);
6380 size[0] = submatrix_size[0];
6381 size[1] = submatrix_size[1];
6385 struct dfm_reader *reader = read_file_open (read->rf);
6387 dfm_reread_record (reader, 1);
6389 if (read->symmetric && size[0] != size[1])
6391 msg (SE, _("Cannot read non-square %zu×%zu matrix "
6392 "using READ with MODE=SYMMETRIC."),
6398 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6399 matrix_read (read, reader, tmp);
6400 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
6403 static struct matrix_cmd *
6404 matrix_parse_write (struct matrix_state *s)
6406 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6407 *cmd = (struct matrix_cmd) {
6411 struct file_handle *fh = NULL;
6412 char *encoding = NULL;
6413 struct write_command *write = &cmd->write;
6414 write->expression = matrix_parse_exp (s);
6415 if (!write->expression)
6419 int repetitions = 0;
6420 int record_width = 0;
6421 enum fmt_type format = FMT_F;
6422 bool has_format = false;
6423 while (lex_match (s->lexer, T_SLASH))
6425 if (lex_match_id (s->lexer, "OUTFILE"))
6427 lex_match (s->lexer, T_EQUALS);
6430 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6434 else if (lex_match_id (s->lexer, "ENCODING"))
6436 lex_match (s->lexer, T_EQUALS);
6437 if (!lex_force_string (s->lexer))
6441 encoding = ss_xstrdup (lex_tokss (s->lexer));
6445 else if (lex_match_id (s->lexer, "FIELD"))
6447 lex_match (s->lexer, T_EQUALS);
6449 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6451 write->c1 = lex_integer (s->lexer);
6453 if (!lex_force_match (s->lexer, T_TO)
6454 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
6456 write->c2 = lex_integer (s->lexer) + 1;
6459 record_width = write->c2 - write->c1;
6460 if (lex_match (s->lexer, T_BY))
6462 if (!lex_force_int_range (s->lexer, "BY", 1,
6463 write->c2 - write->c1))
6465 by = lex_integer (s->lexer);
6468 if (record_width % by)
6470 msg (SE, _("BY %d does not evenly divide record width %d."),
6478 else if (lex_match_id (s->lexer, "MODE"))
6480 lex_match (s->lexer, T_EQUALS);
6481 if (lex_match_id (s->lexer, "RECTANGULAR"))
6482 write->triangular = false;
6483 else if (lex_match_id (s->lexer, "TRIANGULAR"))
6484 write->triangular = true;
6487 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
6491 else if (lex_match_id (s->lexer, "HOLD"))
6493 else if (lex_match_id (s->lexer, "FORMAT"))
6495 if (has_format || write->format)
6497 lex_sbc_only_once ("FORMAT");
6501 lex_match (s->lexer, T_EQUALS);
6503 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6506 const char *p = lex_tokcstr (s->lexer);
6507 if (c_isdigit (p[0]))
6509 repetitions = atoi (p);
6510 p += strspn (p, "0123456789");
6511 if (!fmt_from_name (p, &format))
6513 lex_error (s->lexer, _("Unknown format %s."), p);
6519 else if (fmt_from_name (p, &format))
6526 struct fmt_spec spec;
6527 if (!parse_format_specifier (s->lexer, &spec))
6529 write->format = xmemdup (&spec, sizeof spec);
6534 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
6542 lex_sbc_missing ("FIELD");
6548 if (s->prev_write_file)
6549 fh = fh_ref (s->prev_write_file);
6552 lex_sbc_missing ("OUTFILE");
6556 fh_unref (s->prev_write_file);
6557 s->prev_write_file = fh_ref (fh);
6559 write->wf = write_file_create (s, fh);
6563 free (write->wf->encoding);
6564 write->wf->encoding = encoding;
6568 /* Field width may be specified in multiple ways:
6571 2. The format on FORMAT.
6572 3. The repetition factor on FORMAT.
6574 (2) and (3) are mutually exclusive.
6576 If more than one of these is present, they must agree. If none of them is
6577 present, then free-field format is used.
6579 if (repetitions > record_width)
6581 msg (SE, _("%d repetitions cannot fit in record width %d."),
6582 repetitions, record_width);
6585 int w = (repetitions ? record_width / repetitions
6586 : write->format ? write->format->w
6591 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6592 "which implies field width %d, "
6593 "but BY specifies field width %d."),
6594 repetitions, record_width, w, by);
6596 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6600 if (w && !write->format)
6602 write->format = xmalloc (sizeof *write->format);
6603 *write->format = (struct fmt_spec) { .type = format, .w = w };
6605 if (!fmt_check_output (write->format))
6609 if (write->format && fmt_var_width (write->format) > sizeof (double))
6611 char s[FMT_STRING_LEN_MAX + 1];
6612 fmt_to_string (write->format, s);
6613 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
6614 s, sizeof (double));
6622 matrix_cmd_destroy (cmd);
6627 matrix_cmd_execute_write (struct write_command *write)
6629 gsl_matrix *m = matrix_expr_evaluate (write->expression);
6633 if (write->triangular && m->size1 != m->size2)
6635 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
6636 "the matrix to be written has dimensions %zu×%zu."),
6637 m->size1, m->size2);
6638 gsl_matrix_free (m);
6642 struct dfm_writer *writer = write_file_open (write->wf);
6643 if (!writer || !m->size1)
6645 gsl_matrix_free (m);
6649 const struct fmt_settings *settings = settings_get_fmt_settings ();
6650 struct u8_line *line = write->wf->held;
6651 for (size_t y = 0; y < m->size1; y++)
6655 line = xmalloc (sizeof *line);
6656 u8_line_init (line);
6658 size_t nx = write->triangular ? y + 1 : m->size2;
6660 for (size_t x = 0; x < nx; x++)
6663 double f = gsl_matrix_get (m, y, x);
6667 if (fmt_is_numeric (write->format->type))
6670 v.s = (uint8_t *) &f;
6671 s = data_out (&v, NULL, write->format, settings);
6675 s = xmalloc (DBL_BUFSIZE_BOUND);
6676 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
6677 >= DBL_BUFSIZE_BOUND)
6680 size_t len = strlen (s);
6681 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
6682 if (width + x0 > write->c2)
6684 dfm_put_record_utf8 (writer, line->s.ss.string,
6686 u8_line_clear (line);
6689 u8_line_put (line, x0, x0 + width, s, len);
6692 x0 += write->format ? write->format->w : width + 1;
6695 if (y + 1 >= m->size1 && write->hold)
6697 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
6698 u8_line_clear (line);
6702 u8_line_destroy (line);
6706 write->wf->held = line;
6708 gsl_matrix_free (m);
6711 static struct matrix_cmd *
6712 matrix_parse_get (struct matrix_state *s)
6714 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6715 *cmd = (struct matrix_cmd) {
6718 .dataset = s->dataset,
6719 .user = { .treatment = MGET_ERROR },
6720 .system = { .treatment = MGET_ERROR },
6724 struct get_command *get = &cmd->get;
6725 get->dst = matrix_lvalue_parse (s);
6729 while (lex_match (s->lexer, T_SLASH))
6731 if (lex_match_id (s->lexer, "FILE"))
6733 lex_match (s->lexer, T_EQUALS);
6735 fh_unref (get->file);
6736 if (lex_match (s->lexer, T_ASTERISK))
6740 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6745 else if (lex_match_id (s->lexer, "ENCODING"))
6747 lex_match (s->lexer, T_EQUALS);
6748 if (!lex_force_string (s->lexer))
6751 free (get->encoding);
6752 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
6756 else if (lex_match_id (s->lexer, "VARIABLES"))
6758 lex_match (s->lexer, T_EQUALS);
6762 lex_sbc_only_once ("VARIABLES");
6766 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
6769 else if (lex_match_id (s->lexer, "NAMES"))
6771 lex_match (s->lexer, T_EQUALS);
6772 if (!lex_force_id (s->lexer))
6775 struct substring name = lex_tokss (s->lexer);
6776 get->names = matrix_var_lookup (s, name);
6778 get->names = matrix_var_create (s, name);
6781 else if (lex_match_id (s->lexer, "MISSING"))
6783 lex_match (s->lexer, T_EQUALS);
6784 if (lex_match_id (s->lexer, "ACCEPT"))
6785 get->user.treatment = MGET_ACCEPT;
6786 else if (lex_match_id (s->lexer, "OMIT"))
6787 get->user.treatment = MGET_OMIT;
6788 else if (lex_is_number (s->lexer))
6790 get->user.treatment = MGET_RECODE;
6791 get->user.substitute = lex_number (s->lexer);
6796 lex_error (s->lexer, NULL);
6800 else if (lex_match_id (s->lexer, "SYSMIS"))
6802 lex_match (s->lexer, T_EQUALS);
6803 if (lex_match_id (s->lexer, "OMIT"))
6804 get->system.treatment = MGET_OMIT;
6805 else if (lex_is_number (s->lexer))
6807 get->system.treatment = MGET_RECODE;
6808 get->system.substitute = lex_number (s->lexer);
6813 lex_error (s->lexer, NULL);
6819 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
6820 "MISSING", "SYSMIS");
6825 if (get->user.treatment != MGET_ACCEPT)
6826 get->system.treatment = MGET_ERROR;
6831 matrix_cmd_destroy (cmd);
6836 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
6837 const struct dictionary *dict)
6839 struct variable **vars;
6844 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
6845 &vars, &n_vars, PV_NUMERIC))
6850 n_vars = dict_get_var_cnt (dict);
6851 vars = xnmalloc (n_vars, sizeof *vars);
6852 for (size_t i = 0; i < n_vars; i++)
6854 struct variable *var = dict_get_var (dict, i);
6855 if (!var_is_numeric (var))
6857 msg (SE, _("GET: Variable %s is not numeric."),
6858 var_get_name (var));
6868 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
6869 for (size_t i = 0; i < n_vars; i++)
6871 char s[sizeof (double)];
6873 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
6874 memcpy (&f, s, sizeof f);
6875 gsl_matrix_set (names, i, 0, f);
6878 gsl_matrix_free (get->names->value);
6879 get->names->value = names;
6883 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
6884 long long int casenum = 1;
6886 for (struct ccase *c = casereader_read (reader); c;
6887 c = casereader_read (reader), casenum++)
6889 if (n_rows >= m->size1)
6891 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
6892 for (size_t y = 0; y < n_rows; y++)
6893 for (size_t x = 0; x < n_vars; x++)
6894 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
6895 gsl_matrix_free (m);
6900 for (size_t x = 0; x < n_vars; x++)
6902 const struct variable *var = vars[x];
6903 double d = case_num (c, var);
6906 if (get->system.treatment == MGET_RECODE)
6907 d = get->system.substitute;
6908 else if (get->system.treatment == MGET_OMIT)
6912 msg (SE, _("GET: Variable %s in case %lld "
6913 "is system-missing."),
6914 var_get_name (var), casenum);
6918 else if (var_is_num_missing (var, d, MV_USER))
6920 if (get->user.treatment == MGET_RECODE)
6921 d = get->user.substitute;
6922 else if (get->user.treatment == MGET_OMIT)
6924 else if (get->user.treatment != MGET_ACCEPT)
6926 msg (SE, _("GET: Variable %s in case %lld has user-missing "
6928 var_get_name (var), casenum, d);
6932 gsl_matrix_set (m, n_rows, x, d);
6943 matrix_lvalue_evaluate_and_assign (get->dst, m);
6946 gsl_matrix_free (m);
6951 matrix_open_casereader (const char *command_name,
6952 struct file_handle *file, const char *encoding,
6953 struct dataset *dataset,
6954 struct casereader **readerp, struct dictionary **dictp)
6958 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
6959 return *readerp != NULL;
6963 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
6965 msg (ME, _("The %s command cannot read an empty active file."),
6969 *readerp = proc_open (dataset);
6970 *dictp = dict_ref (dataset_dict (dataset));
6976 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
6977 struct casereader *reader, struct dictionary *dict)
6980 casereader_destroy (reader);
6982 proc_commit (dataset);
6986 matrix_cmd_execute_get (struct get_command *get)
6988 struct casereader *r;
6989 struct dictionary *d;
6990 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
6993 matrix_cmd_execute_get__ (get, r, d);
6994 matrix_close_casereader (get->file, get->dataset, r, d);
6999 match_rowtype (struct lexer *lexer)
7001 static const char *rowtypes[] = {
7002 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7004 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7006 for (size_t i = 0; i < n_rowtypes; i++)
7007 if (lex_match_id (lexer, rowtypes[i]))
7010 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7015 parse_var_names (struct lexer *lexer, struct string_array *sa)
7017 lex_match (lexer, T_EQUALS);
7019 string_array_clear (sa);
7021 struct dictionary *dict = dict_create (get_default_encoding ());
7024 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7025 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7030 for (size_t i = 0; i < n_names; i++)
7031 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7032 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7034 msg (SE, _("Variable name %s is reserved."), names[i]);
7035 for (size_t j = 0; j < n_names; j++)
7041 string_array_clear (sa);
7042 sa->strings = names;
7043 sa->n = sa->allocated = n_names;
7049 msave_common_destroy (struct msave_common *common)
7053 fh_unref (common->outfile);
7054 string_array_destroy (&common->variables);
7055 string_array_destroy (&common->fnames);
7056 string_array_destroy (&common->snames);
7058 for (size_t i = 0; i < common->n_factors; i++)
7059 matrix_expr_destroy (common->factors[i]);
7060 free (common->factors);
7062 for (size_t i = 0; i < common->n_splits; i++)
7063 matrix_expr_destroy (common->splits[i]);
7064 free (common->splits);
7066 dict_unref (common->dict);
7067 casewriter_destroy (common->writer);
7074 compare_variables (const char *keyword,
7075 const struct string_array *new,
7076 const struct string_array *old)
7082 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7083 "on the first MSAVE within MATRIX."), keyword);
7086 else if (!string_array_equal_case (old, new))
7088 msg (SE, _("%s must specify the same variables each time within "
7089 "a given MATRIX."), keyword);
7096 static struct matrix_cmd *
7097 matrix_parse_msave (struct matrix_state *s)
7099 struct msave_common *common = xmalloc (sizeof *common);
7100 *common = (struct msave_common) { .outfile = NULL };
7102 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7103 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7105 struct matrix_expr *splits = NULL;
7106 struct matrix_expr *factors = NULL;
7108 struct msave_command *msave = &cmd->msave;
7109 msave->expr = matrix_parse_exp (s);
7113 while (lex_match (s->lexer, T_SLASH))
7115 if (lex_match_id (s->lexer, "TYPE"))
7117 lex_match (s->lexer, T_EQUALS);
7119 msave->rowtype = match_rowtype (s->lexer);
7120 if (!msave->rowtype)
7123 else if (lex_match_id (s->lexer, "OUTFILE"))
7125 lex_match (s->lexer, T_EQUALS);
7127 fh_unref (common->outfile);
7128 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7129 if (!common->outfile)
7132 else if (lex_match_id (s->lexer, "VARIABLES"))
7134 if (!parse_var_names (s->lexer, &common->variables))
7137 else if (lex_match_id (s->lexer, "FNAMES"))
7139 if (!parse_var_names (s->lexer, &common->fnames))
7142 else if (lex_match_id (s->lexer, "SNAMES"))
7144 if (!parse_var_names (s->lexer, &common->snames))
7147 else if (lex_match_id (s->lexer, "SPLIT"))
7149 lex_match (s->lexer, T_EQUALS);
7151 matrix_expr_destroy (splits);
7152 splits = matrix_parse_exp (s);
7156 else if (lex_match_id (s->lexer, "FACTOR"))
7158 lex_match (s->lexer, T_EQUALS);
7160 matrix_expr_destroy (factors);
7161 factors = matrix_parse_exp (s);
7167 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7168 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7172 if (!msave->rowtype)
7174 lex_sbc_missing ("TYPE");
7180 if (common->fnames.n && !factors)
7182 msg (SE, _("FNAMES requires FACTOR."));
7185 if (common->snames.n && !splits)
7187 msg (SE, _("SNAMES requires SPLIT."));
7190 if (!common->outfile)
7192 lex_sbc_missing ("OUTFILE");
7199 if (common->outfile)
7201 if (!fh_equal (common->outfile, s->common->outfile))
7203 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7204 "within a single MATRIX command."));
7208 if (!compare_variables ("VARIABLES",
7209 &common->variables, &s->common->variables)
7210 || !compare_variables ("FNAMES", &common->fnames, &s->common->fnames)
7211 || !compare_variables ("SNAMES", &common->snames, &s->common->snames))
7213 msave_common_destroy (common);
7215 msave->common = s->common;
7217 struct msave_common *c = s->common;
7220 if (c->n_factors >= c->allocated_factors)
7221 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7222 sizeof *c->factors);
7223 c->factors[c->n_factors++] = factors;
7225 if (c->n_factors > 0)
7226 msave->factors = c->factors[c->n_factors - 1];
7230 if (c->n_splits >= c->allocated_splits)
7231 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7233 c->splits[c->n_splits++] = splits;
7235 if (c->n_splits > 0)
7236 msave->splits = c->splits[c->n_splits - 1];
7241 matrix_expr_destroy (splits);
7242 matrix_expr_destroy (factors);
7243 msave_common_destroy (common);
7244 matrix_cmd_destroy (cmd);
7249 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7251 gsl_matrix *m = matrix_expr_evaluate (e);
7257 msg (SE, _("%s expression must evaluate to vector, "
7258 "not a %zu×%zu matrix."),
7259 name, m->size1, m->size2);
7260 gsl_matrix_free (m);
7264 return matrix_to_vector (m);
7268 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7270 for (size_t i = 0; i < vars->n; i++)
7271 if (!dict_create_var (d, vars->strings[i], 0))
7272 return vars->strings[i];
7276 static struct dictionary *
7277 msave_create_dict (const struct msave_common *common)
7279 struct dictionary *dict = dict_create (get_default_encoding ());
7281 const char *dup_split = msave_add_vars (dict, &common->snames);
7284 /* Should not be possible because the parser ensures that the names are
7289 dict_create_var_assert (dict, "ROWTYPE_", 8);
7291 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7294 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
7298 dict_create_var_assert (dict, "VARNAME_", 8);
7300 const char *dup_var = msave_add_vars (dict, &common->variables);
7303 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
7315 matrix_cmd_execute_msave (struct msave_command *msave)
7317 struct msave_common *common = msave->common;
7318 gsl_matrix *m = NULL;
7319 gsl_vector *factors = NULL;
7320 gsl_vector *splits = NULL;
7322 m = matrix_expr_evaluate (msave->expr);
7326 if (!common->variables.n)
7327 for (size_t i = 0; i < m->size2; i++)
7328 string_array_append_nocopy (&common->variables,
7329 xasprintf ("COL%zu", i + 1));
7330 else if (m->size2 != common->variables.n)
7333 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7334 m->size2, common->variables.n);
7340 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7344 if (!common->fnames.n)
7345 for (size_t i = 0; i < factors->size; i++)
7346 string_array_append_nocopy (&common->fnames,
7347 xasprintf ("FAC%zu", i + 1));
7348 else if (factors->size != common->fnames.n)
7350 msg (SE, _("There are %zu factor variables, "
7351 "but %zu factor values were supplied."),
7352 common->fnames.n, factors->size);
7359 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7363 if (!common->snames.n)
7364 for (size_t i = 0; i < splits->size; i++)
7365 string_array_append_nocopy (&common->snames,
7366 xasprintf ("SPL%zu", i + 1));
7367 else if (splits->size != common->snames.n)
7369 msg (SE, _("There are %zu split variables, "
7370 "but %zu split values were supplied."),
7371 common->snames.n, splits->size);
7376 if (!common->writer)
7378 struct dictionary *dict = msave_create_dict (common);
7382 common->writer = any_writer_open (common->outfile, dict);
7383 if (!common->writer)
7389 common->dict = dict;
7392 bool matrix = (!strcmp (msave->rowtype, "COV")
7393 || !strcmp (msave->rowtype, "CORR"));
7394 for (size_t y = 0; y < m->size1; y++)
7396 struct ccase *c = case_create (dict_get_proto (common->dict));
7399 /* Split variables */
7401 for (size_t i = 0; i < splits->size; i++)
7402 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
7405 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7406 msave->rowtype, ' ');
7410 for (size_t i = 0; i < factors->size; i++)
7411 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
7414 const char *varname_ = (matrix && y < common->variables.n
7415 ? common->variables.strings[y]
7417 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7420 /* Continuous variables. */
7421 for (size_t x = 0; x < m->size2; x++)
7422 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
7423 casewriter_write (common->writer, c);
7427 gsl_matrix_free (m);
7428 gsl_vector_free (factors);
7429 gsl_vector_free (splits);
7432 static struct matrix_cmd *
7433 matrix_parse_mget (struct matrix_state *s)
7435 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7436 *cmd = (struct matrix_cmd) {
7440 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
7444 struct mget_command *mget = &cmd->mget;
7446 lex_match (s->lexer, T_SLASH);
7447 while (lex_token (s->lexer) != T_ENDCMD)
7449 if (lex_match_id (s->lexer, "FILE"))
7451 lex_match (s->lexer, T_EQUALS);
7453 fh_unref (mget->file);
7454 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7458 else if (lex_match_id (s->lexer, "ENCODING"))
7460 lex_match (s->lexer, T_EQUALS);
7461 if (!lex_force_string (s->lexer))
7464 free (mget->encoding);
7465 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
7469 else if (lex_match_id (s->lexer, "TYPE"))
7471 lex_match (s->lexer, T_EQUALS);
7472 while (lex_token (s->lexer) != T_SLASH
7473 && lex_token (s->lexer) != T_ENDCMD)
7475 const char *rowtype = match_rowtype (s->lexer);
7479 stringi_set_insert (&mget->rowtypes, rowtype);
7484 lex_error_expecting (s->lexer, "FILE", "TYPE");
7487 lex_match (s->lexer, T_SLASH);
7492 matrix_cmd_destroy (cmd);
7496 static const struct variable *
7497 get_a8_var (const struct dictionary *d, const char *name)
7499 const struct variable *v = dict_lookup_var (d, name);
7502 msg (SE, _("Matrix data file lacks %s variable."), name);
7505 if (var_get_width (v) != 8)
7507 msg (SE, _("%s variable in matrix data file must be "
7508 "8-byte string, but it has width %d."),
7509 name, var_get_width (v));
7516 var_changed (const struct ccase *ca, const struct ccase *cb,
7517 const struct variable *var)
7520 ? !value_equal (case_data (ca, var), case_data (cb, var),
7521 var_get_width (var))
7526 vars_changed (const struct ccase *ca, const struct ccase *cb,
7527 const struct dictionary *d,
7528 size_t first_var, size_t n_vars)
7530 for (size_t i = 0; i < n_vars; i++)
7532 const struct variable *v = dict_get_var (d, first_var + i);
7533 if (var_changed (ca, cb, v))
7540 vars_all_missing (const struct ccase *c, const struct dictionary *d,
7541 size_t first_var, size_t n_vars)
7543 for (size_t i = 0; i < n_vars; i++)
7544 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
7550 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
7551 const struct dictionary *d,
7552 const struct variable *rowtype_var,
7553 const struct stringi_set *accepted_rowtypes,
7554 struct matrix_state *s,
7555 size_t ss, size_t sn, size_t si,
7556 size_t fs, size_t fn, size_t fi,
7557 size_t cs, size_t cn,
7558 struct pivot_table *pt,
7559 struct pivot_dimension *var_dimension)
7564 /* Is this a matrix for pooled data, either where there are no factor
7565 variables or the factor variables are missing? */
7566 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
7568 struct substring rowtype = case_ss (rows[0], rowtype_var);
7569 ss_rtrim (&rowtype, ss_cstr (" "));
7570 if (!stringi_set_is_empty (accepted_rowtypes)
7571 && !stringi_set_contains_len (accepted_rowtypes,
7572 rowtype.string, rowtype.length))
7575 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
7576 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
7577 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
7578 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
7579 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
7580 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
7584 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
7585 (int) rowtype.length, rowtype.string);
7589 struct string name = DS_EMPTY_INITIALIZER;
7590 ds_put_cstr (&name, prefix);
7592 ds_put_format (&name, "F%zu", fi);
7594 ds_put_format (&name, "S%zu", si);
7596 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
7598 mv = matrix_var_create (s, ds_ss (&name));
7601 msg (SW, _("Matrix data file contains variable with existing name %s."),
7603 goto exit_free_name;
7606 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
7607 size_t n_missing = 0;
7608 for (size_t y = 0; y < n_rows; y++)
7610 for (size_t x = 0; x < cn; x++)
7612 struct variable *var = dict_get_var (d, cs + x);
7613 double value = case_num (rows[y], var);
7614 if (var_is_num_missing (var, value, MV_ANY))
7619 gsl_matrix_set (m, y, x, value);
7623 int var_index = pivot_category_create_leaf (
7624 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
7625 double values[] = { n_rows, cn };
7626 for (size_t j = 0; j < sn; j++)
7628 struct variable *var = dict_get_var (d, ss + j);
7629 const union value *value = case_data (rows[0], var);
7630 pivot_table_put2 (pt, j, var_index,
7631 pivot_value_new_var_value (var, value));
7633 for (size_t j = 0; j < fn; j++)
7635 struct variable *var = dict_get_var (d, fs + j);
7636 const union value sysmis = { .f = SYSMIS };
7637 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
7638 pivot_table_put2 (pt, j + sn, var_index,
7639 pivot_value_new_var_value (var, value));
7641 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
7642 pivot_table_put2 (pt, j + sn + fn, var_index,
7643 pivot_value_new_integer (values[j]));
7646 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
7647 "value, which was treated as zero.",
7648 "Matrix data file variable %s contains %zu missing "
7649 "values, which were treated as zero.", n_missing),
7650 ds_cstr (&name), n_missing);
7657 for (size_t y = 0; y < n_rows; y++)
7658 case_unref (rows[y]);
7662 matrix_cmd_execute_mget__ (struct mget_command *mget,
7663 struct casereader *r, const struct dictionary *d)
7665 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
7666 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
7667 if (!rowtype_ || !varname_)
7670 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
7672 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
7675 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
7677 msg (SE, _("Matrix data file contains no continuous variables."));
7681 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
7683 const struct variable *v = dict_get_var (d, i);
7684 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
7687 _("Matrix data file contains unexpected string variable %s."),
7693 /* SPLIT variables. */
7695 size_t sn = var_get_dict_index (rowtype_);
7696 struct ccase *sc = NULL;
7699 /* FACTOR variables. */
7700 size_t fs = var_get_dict_index (rowtype_) + 1;
7701 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
7702 struct ccase *fc = NULL;
7705 /* Continuous variables. */
7706 size_t cs = var_get_dict_index (varname_) + 1;
7707 size_t cn = dict_get_var_cnt (d) - cs;
7708 struct ccase *cc = NULL;
7711 struct pivot_table *pt = pivot_table_create (
7712 N_("Matrix Variables Created by MGET"));
7713 struct pivot_dimension *attr_dimension = pivot_dimension_create (
7714 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
7715 struct pivot_dimension *var_dimension = pivot_dimension_create (
7716 pt, PIVOT_AXIS_ROW, N_("Variable"));
7719 struct pivot_category *splits = pivot_category_create_group (
7720 attr_dimension->root, N_("Split Values"));
7721 for (size_t i = 0; i < sn; i++)
7722 pivot_category_create_leaf (splits, pivot_value_new_variable (
7723 dict_get_var (d, ss + i)));
7727 struct pivot_category *factors = pivot_category_create_group (
7728 attr_dimension->root, N_("Factors"));
7729 for (size_t i = 0; i < fn; i++)
7730 pivot_category_create_leaf (factors, pivot_value_new_variable (
7731 dict_get_var (d, fs + i)));
7733 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
7734 N_("Rows"), N_("Columns"));
7737 struct ccase **rows = NULL;
7738 size_t allocated_rows = 0;
7742 while ((c = casereader_read (r)) != NULL)
7744 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
7754 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
7755 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
7756 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
7759 if (change != NOTHING_CHANGED)
7761 matrix_mget_commit_var (rows, n_rows, d,
7762 rowtype_, &mget->rowtypes,
7773 if (n_rows >= allocated_rows)
7774 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
7777 if (change == SPLITS_CHANGED)
7783 /* Reset the factor number, if there are factors. */
7787 if (row_has_factors)
7793 else if (change == FACTORS_CHANGED)
7795 if (row_has_factors)
7801 matrix_mget_commit_var (rows, n_rows, d,
7802 rowtype_, &mget->rowtypes,
7814 if (var_dimension->n_leaves)
7815 pivot_table_submit (pt);
7817 pivot_table_unref (pt);
7821 matrix_cmd_execute_mget (struct mget_command *mget)
7823 struct casereader *r;
7824 struct dictionary *d;
7825 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
7826 mget->state->dataset, &r, &d))
7828 matrix_cmd_execute_mget__ (mget, r, d);
7829 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
7834 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
7836 if (!lex_force_id (s->lexer))
7839 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
7841 *varp = matrix_var_create (s, lex_tokss (s->lexer));
7846 static struct matrix_cmd *
7847 matrix_parse_eigen (struct matrix_state *s)
7849 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7850 *cmd = (struct matrix_cmd) {
7852 .eigen = { .expr = NULL }
7855 struct eigen_command *eigen = &cmd->eigen;
7856 if (!lex_force_match (s->lexer, T_LPAREN))
7858 eigen->expr = matrix_parse_expr (s);
7860 || !lex_force_match (s->lexer, T_COMMA)
7861 || !matrix_parse_dst_var (s, &eigen->evec)
7862 || !lex_force_match (s->lexer, T_COMMA)
7863 || !matrix_parse_dst_var (s, &eigen->eval)
7864 || !lex_force_match (s->lexer, T_RPAREN))
7870 matrix_cmd_destroy (cmd);
7875 matrix_cmd_execute_eigen (struct eigen_command *eigen)
7877 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
7880 if (!is_symmetric (A))
7882 msg (SE, _("Argument of EIGEN must be symmetric."));
7883 gsl_matrix_free (A);
7887 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
7888 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
7889 gsl_vector v_eval = to_vector (eval);
7890 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
7891 gsl_eigen_symmv (A, &v_eval, evec, w);
7892 gsl_eigen_symmv_free (w);
7894 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
7896 gsl_matrix_free (A);
7898 gsl_matrix_free (eigen->eval->value);
7899 eigen->eval->value = eval;
7901 gsl_matrix_free (eigen->evec->value);
7902 eigen->evec->value = evec;
7905 static struct matrix_cmd *
7906 matrix_parse_setdiag (struct matrix_state *s)
7908 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7909 *cmd = (struct matrix_cmd) {
7910 .type = MCMD_SETDIAG,
7911 .setdiag = { .dst = NULL }
7914 struct setdiag_command *setdiag = &cmd->setdiag;
7915 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
7918 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
7921 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
7926 if (!lex_force_match (s->lexer, T_COMMA))
7929 setdiag->expr = matrix_parse_expr (s);
7933 if (!lex_force_match (s->lexer, T_RPAREN))
7939 matrix_cmd_destroy (cmd);
7944 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
7946 gsl_matrix *dst = setdiag->dst->value;
7949 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
7950 setdiag->dst->name);
7954 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
7958 size_t n = MIN (dst->size1, dst->size2);
7959 if (is_scalar (src))
7961 double d = to_scalar (src);
7962 for (size_t i = 0; i < n; i++)
7963 gsl_matrix_set (dst, i, i, d);
7965 else if (is_vector (src))
7967 gsl_vector v = to_vector (src);
7968 for (size_t i = 0; i < n && i < v.size; i++)
7969 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
7972 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
7973 "not a %zu×%zu matrix."),
7974 src->size1, src->size2);
7975 gsl_matrix_free (src);
7978 static struct matrix_cmd *
7979 matrix_parse_svd (struct matrix_state *s)
7981 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7982 *cmd = (struct matrix_cmd) {
7984 .svd = { .expr = NULL }
7987 struct svd_command *svd = &cmd->svd;
7988 if (!lex_force_match (s->lexer, T_LPAREN))
7990 svd->expr = matrix_parse_expr (s);
7992 || !lex_force_match (s->lexer, T_COMMA)
7993 || !matrix_parse_dst_var (s, &svd->u)
7994 || !lex_force_match (s->lexer, T_COMMA)
7995 || !matrix_parse_dst_var (s, &svd->s)
7996 || !lex_force_match (s->lexer, T_COMMA)
7997 || !matrix_parse_dst_var (s, &svd->v)
7998 || !lex_force_match (s->lexer, T_RPAREN))
8004 matrix_cmd_destroy (cmd);
8009 matrix_cmd_execute_svd (struct svd_command *svd)
8011 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8015 if (m->size1 >= m->size2)
8018 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8019 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8020 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8021 gsl_vector *work = gsl_vector_alloc (A->size2);
8022 gsl_linalg_SV_decomp (A, V, &Sv, work);
8023 gsl_vector_free (work);
8025 matrix_var_set (svd->u, A);
8026 matrix_var_set (svd->s, S);
8027 matrix_var_set (svd->v, V);
8031 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8032 gsl_matrix_transpose_memcpy (At, m);
8033 gsl_matrix_free (m);
8035 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8036 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8037 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8038 gsl_vector *work = gsl_vector_alloc (At->size2);
8039 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8040 gsl_vector_free (work);
8042 matrix_var_set (svd->v, At);
8043 matrix_var_set (svd->s, St);
8044 matrix_var_set (svd->u, Vt);
8049 matrix_cmd_execute (struct matrix_cmd *cmd)
8054 matrix_cmd_execute_compute (&cmd->compute);
8058 matrix_cmd_execute_print (&cmd->print);
8062 return matrix_cmd_execute_do_if (&cmd->do_if);
8065 matrix_cmd_execute_loop (&cmd->loop);
8072 matrix_cmd_execute_display (&cmd->display);
8076 matrix_cmd_execute_release (&cmd->release);
8080 matrix_cmd_execute_save (&cmd->save);
8084 matrix_cmd_execute_read (&cmd->read);
8088 matrix_cmd_execute_write (&cmd->write);
8092 matrix_cmd_execute_get (&cmd->get);
8096 matrix_cmd_execute_msave (&cmd->msave);
8100 matrix_cmd_execute_mget (&cmd->mget);
8104 matrix_cmd_execute_eigen (&cmd->eigen);
8108 matrix_cmd_execute_setdiag (&cmd->setdiag);
8112 matrix_cmd_execute_svd (&cmd->svd);
8120 matrix_cmds_uninit (struct matrix_cmds *cmds)
8122 for (size_t i = 0; i < cmds->n; i++)
8123 matrix_cmd_destroy (cmds->commands[i]);
8124 free (cmds->commands);
8128 matrix_cmd_destroy (struct matrix_cmd *cmd)
8136 matrix_lvalue_destroy (cmd->compute.lvalue);
8137 matrix_expr_destroy (cmd->compute.rvalue);
8141 matrix_expr_destroy (cmd->print.expression);
8142 free (cmd->print.title);
8143 print_labels_destroy (cmd->print.rlabels);
8144 print_labels_destroy (cmd->print.clabels);
8148 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8150 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8151 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
8153 free (cmd->do_if.clauses);
8157 matrix_expr_destroy (cmd->loop.start);
8158 matrix_expr_destroy (cmd->loop.end);
8159 matrix_expr_destroy (cmd->loop.increment);
8160 matrix_expr_destroy (cmd->loop.top_condition);
8161 matrix_expr_destroy (cmd->loop.bottom_condition);
8162 matrix_cmds_uninit (&cmd->loop.commands);
8172 free (cmd->release.vars);
8176 matrix_expr_destroy (cmd->save.expression);
8180 matrix_lvalue_destroy (cmd->read.dst);
8181 matrix_expr_destroy (cmd->read.size);
8185 matrix_expr_destroy (cmd->write.expression);
8186 free (cmd->write.format);
8190 matrix_lvalue_destroy (cmd->get.dst);
8191 fh_unref (cmd->get.file);
8192 free (cmd->get.encoding);
8193 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8197 matrix_expr_destroy (cmd->msave.expr);
8201 fh_unref (cmd->mget.file);
8202 stringi_set_destroy (&cmd->mget.rowtypes);
8206 matrix_expr_destroy (cmd->eigen.expr);
8210 matrix_expr_destroy (cmd->setdiag.expr);
8214 matrix_expr_destroy (cmd->svd.expr);
8220 struct matrix_command_name
8223 struct matrix_cmd *(*parse) (struct matrix_state *);
8226 static const struct matrix_command_name *
8227 matrix_parse_command_name (struct lexer *lexer)
8229 static const struct matrix_command_name commands[] = {
8230 { "COMPUTE", matrix_parse_compute },
8231 { "CALL EIGEN", matrix_parse_eigen },
8232 { "CALL SETDIAG", matrix_parse_setdiag },
8233 { "CALL SVD", matrix_parse_svd },
8234 { "PRINT", matrix_parse_print },
8235 { "DO IF", matrix_parse_do_if },
8236 { "LOOP", matrix_parse_loop },
8237 { "BREAK", matrix_parse_break },
8238 { "READ", matrix_parse_read },
8239 { "WRITE", matrix_parse_write },
8240 { "GET", matrix_parse_get },
8241 { "SAVE", matrix_parse_save },
8242 { "MGET", matrix_parse_mget },
8243 { "MSAVE", matrix_parse_msave },
8244 { "DISPLAY", matrix_parse_display },
8245 { "RELEASE", matrix_parse_release },
8247 static size_t n = sizeof commands / sizeof *commands;
8249 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8250 if (lex_match_phrase (lexer, c->name))
8255 static struct matrix_cmd *
8256 matrix_parse_command (struct matrix_state *s)
8258 size_t nesting_level = SIZE_MAX;
8260 struct matrix_cmd *c = NULL;
8261 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
8263 lex_error (s->lexer, _("Unknown matrix command."));
8264 else if (!cmd->parse)
8265 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8269 nesting_level = output_open_group (
8270 group_item_create_nocopy (utf8_to_title (cmd->name),
8271 utf8_to_title (cmd->name)));
8276 lex_end_of_command (s->lexer);
8277 lex_discard_rest_of_command (s->lexer);
8278 if (nesting_level != SIZE_MAX)
8279 output_close_groups (nesting_level);
8285 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8287 if (!lex_force_match (lexer, T_ENDCMD))
8290 struct matrix_state state = {
8292 .session = dataset_session (ds),
8294 .vars = HMAP_INITIALIZER (state.vars),
8299 while (lex_match (lexer, T_ENDCMD))
8301 if (lex_token (lexer) == T_STOP)
8303 msg (SE, _("Unexpected end of input expecting matrix command."));
8307 if (lex_match_phrase (lexer, "END MATRIX"))
8310 struct matrix_cmd *c = matrix_parse_command (&state);
8313 matrix_cmd_execute (c);
8314 matrix_cmd_destroy (c);
8318 struct matrix_var *var, *next;
8319 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8322 gsl_matrix_free (var->value);
8323 hmap_delete (&state.vars, &var->hmap_node);
8326 hmap_destroy (&state.vars);
8327 msave_common_destroy (state.common);
8328 fh_unref (state.prev_read_file);
8329 for (size_t i = 0; i < state.n_read_files; i++)
8330 read_file_destroy (state.read_files[i]);
8331 free (state.read_files);
8332 fh_unref (state.prev_write_file);
8333 for (size_t i = 0; i < state.n_write_files; i++)
8334 write_file_destroy (state.write_files[i]);
8335 free (state.write_files);
8336 fh_unref (state.prev_save_file);
8337 for (size_t i = 0; i < state.n_save_files; i++)
8338 save_file_destroy (state.save_files[i]);
8339 free (state.save_files);