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 size_t allocated_subs = 0;
2646 struct matrix_expr *sub = matrix_parse_expr (s);
2650 if (e->n_subs >= allocated_subs)
2651 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
2652 e->subs[e->n_subs++] = sub;
2654 while (lex_match (s->lexer, T_COMMA));
2655 if (!lex_force_match (s->lexer, T_RPAREN))
2662 matrix_expr_destroy (e);
2666 static struct matrix_expr *
2667 matrix_parse_primary (struct matrix_state *s)
2669 if (lex_is_number (s->lexer))
2671 double number = lex_number (s->lexer);
2673 return matrix_expr_create_number (number);
2675 else if (lex_is_string (s->lexer))
2677 char string[sizeof (double)];
2678 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
2682 memcpy (&number, string, sizeof number);
2683 return matrix_expr_create_number (number);
2685 else if (lex_match (s->lexer, T_LPAREN))
2687 struct matrix_expr *e = matrix_parse_or_xor (s);
2688 if (!e || !lex_force_match (s->lexer, T_RPAREN))
2690 matrix_expr_destroy (e);
2695 else if (lex_match (s->lexer, T_LCURLY))
2697 struct matrix_expr *e = matrix_parse_curly_semi (s);
2698 if (!e || !lex_force_match (s->lexer, T_RCURLY))
2700 matrix_expr_destroy (e);
2705 else if (lex_token (s->lexer) == T_ID)
2707 struct matrix_expr *retval;
2708 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
2711 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
2714 lex_error (s->lexer, _("Unknown variable %s."),
2715 lex_tokcstr (s->lexer));
2720 struct matrix_expr *e = xmalloc (sizeof *e);
2721 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
2724 else if (lex_token (s->lexer) == T_ALL)
2726 struct matrix_expr *retval;
2727 if (matrix_parse_function (s, "ALL", &retval))
2731 lex_error (s->lexer, NULL);
2735 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
2738 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
2740 if (lex_match (s->lexer, T_COLON))
2747 *indexp = matrix_parse_expr (s);
2748 return *indexp != NULL;
2752 static struct matrix_expr *
2753 matrix_parse_postfix (struct matrix_state *s)
2755 struct matrix_expr *lhs = matrix_parse_primary (s);
2756 if (!lhs || !lex_match (s->lexer, T_LPAREN))
2759 struct matrix_expr *i0;
2760 if (!matrix_parse_index_expr (s, &i0))
2762 matrix_expr_destroy (lhs);
2765 if (lex_match (s->lexer, T_RPAREN))
2767 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
2768 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
2769 else if (lex_match (s->lexer, T_COMMA))
2771 struct matrix_expr *i1;
2772 if (!matrix_parse_index_expr (s, &i1)
2773 || !lex_force_match (s->lexer, T_RPAREN))
2775 matrix_expr_destroy (lhs);
2776 matrix_expr_destroy (i0);
2777 matrix_expr_destroy (i1);
2780 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
2781 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
2782 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
2787 lex_error_expecting (s->lexer, "`)'", "`,'");
2792 static struct matrix_expr *
2793 matrix_parse_unary (struct matrix_state *s)
2795 if (lex_match (s->lexer, T_DASH))
2797 struct matrix_expr *lhs = matrix_parse_unary (s);
2798 return lhs ? matrix_expr_create_1 (MOP_NEGATE, lhs) : NULL;
2800 else if (lex_match (s->lexer, T_PLUS))
2801 return matrix_parse_unary (s);
2803 return matrix_parse_postfix (s);
2806 static struct matrix_expr *
2807 matrix_parse_seq (struct matrix_state *s)
2809 struct matrix_expr *start = matrix_parse_unary (s);
2810 if (!start || !lex_match (s->lexer, T_COLON))
2813 struct matrix_expr *end = matrix_parse_unary (s);
2816 matrix_expr_destroy (start);
2820 if (lex_match (s->lexer, T_COLON))
2822 struct matrix_expr *increment = matrix_parse_unary (s);
2825 matrix_expr_destroy (start);
2826 matrix_expr_destroy (end);
2829 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2832 return matrix_expr_create_2 (MOP_SEQ, start, end);
2835 static struct matrix_expr *
2836 matrix_parse_exp (struct matrix_state *s)
2838 struct matrix_expr *lhs = matrix_parse_seq (s);
2845 if (lex_match (s->lexer, T_EXP))
2847 else if (lex_match_phrase (s->lexer, "&**"))
2852 struct matrix_expr *rhs = matrix_parse_seq (s);
2855 matrix_expr_destroy (lhs);
2858 lhs = matrix_expr_create_2 (op, lhs, rhs);
2862 static struct matrix_expr *
2863 matrix_parse_mul_div (struct matrix_state *s)
2865 struct matrix_expr *lhs = matrix_parse_exp (s);
2872 if (lex_match (s->lexer, T_ASTERISK))
2874 else if (lex_match (s->lexer, T_SLASH))
2876 else if (lex_match_phrase (s->lexer, "&*"))
2878 else if (lex_match_phrase (s->lexer, "&/"))
2883 struct matrix_expr *rhs = matrix_parse_exp (s);
2886 matrix_expr_destroy (lhs);
2889 lhs = matrix_expr_create_2 (op, lhs, rhs);
2893 static struct matrix_expr *
2894 matrix_parse_add_sub (struct matrix_state *s)
2896 struct matrix_expr *lhs = matrix_parse_mul_div (s);
2903 if (lex_match (s->lexer, T_PLUS))
2905 else if (lex_match (s->lexer, T_DASH))
2907 else if (lex_token (s->lexer) == T_NEG_NUM)
2912 struct matrix_expr *rhs = matrix_parse_mul_div (s);
2915 matrix_expr_destroy (lhs);
2918 lhs = matrix_expr_create_2 (op, lhs, rhs);
2922 static struct matrix_expr *
2923 matrix_parse_relational (struct matrix_state *s)
2925 struct matrix_expr *lhs = matrix_parse_add_sub (s);
2932 if (lex_match (s->lexer, T_GT))
2934 else if (lex_match (s->lexer, T_GE))
2936 else if (lex_match (s->lexer, T_LT))
2938 else if (lex_match (s->lexer, T_LE))
2940 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
2942 else if (lex_match (s->lexer, T_NE))
2947 struct matrix_expr *rhs = matrix_parse_add_sub (s);
2950 matrix_expr_destroy (lhs);
2953 lhs = matrix_expr_create_2 (op, lhs, rhs);
2957 static struct matrix_expr *
2958 matrix_parse_not (struct matrix_state *s)
2960 if (lex_match (s->lexer, T_NOT))
2962 struct matrix_expr *lhs = matrix_parse_not (s);
2963 return lhs ? matrix_expr_create_1 (MOP_NOT, lhs) : NULL;
2966 return matrix_parse_relational (s);
2969 static struct matrix_expr *
2970 matrix_parse_and (struct matrix_state *s)
2972 struct matrix_expr *lhs = matrix_parse_not (s);
2976 while (lex_match (s->lexer, T_AND))
2978 struct matrix_expr *rhs = matrix_parse_not (s);
2981 matrix_expr_destroy (lhs);
2984 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
2989 static struct matrix_expr *
2990 matrix_parse_or_xor (struct matrix_state *s)
2992 struct matrix_expr *lhs = matrix_parse_and (s);
2999 if (lex_match (s->lexer, T_OR))
3001 else if (lex_match_id (s->lexer, "XOR"))
3006 struct matrix_expr *rhs = matrix_parse_and (s);
3009 matrix_expr_destroy (lhs);
3012 lhs = matrix_expr_create_2 (op, lhs, rhs);
3016 static struct matrix_expr *
3017 matrix_parse_expr (struct matrix_state *s)
3019 return matrix_parse_or_xor (s);
3022 /* Expression evaluation. */
3025 matrix_op_eval (enum matrix_op op, double a, double b)
3029 case MOP_ADD_ELEMS: return a + b;
3030 case MOP_SUB_ELEMS: return a - b;
3031 case MOP_MUL_ELEMS: return a * b;
3032 case MOP_DIV_ELEMS: return a / b;
3033 case MOP_EXP_ELEMS: return pow (a, b);
3034 case MOP_GT: return a > b;
3035 case MOP_GE: return a >= b;
3036 case MOP_LT: return a < b;
3037 case MOP_LE: return a <= b;
3038 case MOP_EQ: return a == b;
3039 case MOP_NE: return a != b;
3040 case MOP_AND: return (a > 0) && (b > 0);
3041 case MOP_OR: return (a > 0) || (b > 0);
3042 case MOP_XOR: return (a > 0) != (b > 0);
3044 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3053 case MOP_PASTE_HORZ:
3054 case MOP_PASTE_VERT:
3070 matrix_op_name (enum matrix_op op)
3074 case MOP_ADD_ELEMS: return "+";
3075 case MOP_SUB_ELEMS: return "-";
3076 case MOP_MUL_ELEMS: return "&*";
3077 case MOP_DIV_ELEMS: return "&/";
3078 case MOP_EXP_ELEMS: return "&**";
3079 case MOP_GT: return ">";
3080 case MOP_GE: return ">=";
3081 case MOP_LT: return "<";
3082 case MOP_LE: return "<=";
3083 case MOP_EQ: return "=";
3084 case MOP_NE: return "<>";
3085 case MOP_AND: return "AND";
3086 case MOP_OR: return "OR";
3087 case MOP_XOR: return "XOR";
3089 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3098 case MOP_PASTE_HORZ:
3099 case MOP_PASTE_VERT:
3115 is_scalar (const gsl_matrix *m)
3117 return m->size1 == 1 && m->size2 == 1;
3121 to_scalar (const gsl_matrix *m)
3123 assert (is_scalar (m));
3124 return gsl_matrix_get (m, 0, 0);
3128 matrix_expr_evaluate_elementwise (enum matrix_op op,
3129 gsl_matrix *a, gsl_matrix *b)
3133 double be = to_scalar (b);
3134 for (size_t r = 0; r < a->size1; r++)
3135 for (size_t c = 0; c < a->size2; c++)
3137 double *ae = gsl_matrix_ptr (a, r, c);
3138 *ae = matrix_op_eval (op, *ae, be);
3142 else if (is_scalar (a))
3144 double ae = to_scalar (a);
3145 for (size_t r = 0; r < b->size1; r++)
3146 for (size_t c = 0; c < b->size2; c++)
3148 double *be = gsl_matrix_ptr (b, r, c);
3149 *be = matrix_op_eval (op, ae, *be);
3153 else if (a->size1 == b->size1 && a->size2 == b->size2)
3155 for (size_t r = 0; r < a->size1; r++)
3156 for (size_t c = 0; c < a->size2; c++)
3158 double *ae = gsl_matrix_ptr (a, r, c);
3159 double be = gsl_matrix_get (b, r, c);
3160 *ae = matrix_op_eval (op, *ae, be);
3166 msg (SE, _("Operands to %s must have the same dimensions or one "
3167 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
3168 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
3174 matrix_expr_evaluate_mul_mat (gsl_matrix *a, gsl_matrix *b)
3176 if (is_scalar (a) || is_scalar (b))
3177 return matrix_expr_evaluate_elementwise (MOP_MUL_ELEMS, a, b);
3179 if (a->size2 != b->size1)
3181 msg (SE, _("Matrices with dimensions %zu×%zu and %zu×%zu are "
3182 "not conformable for multiplication."),
3183 a->size1, a->size2, b->size1, b->size2);
3187 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3188 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3193 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3195 gsl_matrix *tmp = *a;
3201 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3204 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3205 swap_matrix (z, tmp);
3209 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3211 mul_matrix (x, *x, *x, tmp);
3215 matrix_expr_evaluate_exp_mat (gsl_matrix *x_, gsl_matrix *b)
3218 if (x->size1 != x->size2)
3220 msg (SE, _("Matrix exponentation with ** requires a square matrix on "
3221 "the left-hand size, not one with dimensions %zu×%zu."),
3222 x->size1, x->size2);
3227 msg (SE, _("Matrix exponentiation with ** requires a scalar on the "
3228 "right-hand side, not a matrix with dimensions %zu×%zu."),
3229 b->size1, b->size2);
3232 double bf = to_scalar (b);
3233 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3235 msg (SE, _("Exponent %.1f in matrix multiplication is non-integer "
3236 "or outside the valid range."), bf);
3241 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3243 gsl_matrix_set_identity (y);
3247 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3249 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3252 mul_matrix (&y, x, y, &t);
3253 square_matrix (&x, &t);
3256 square_matrix (&x, &t);
3258 mul_matrix (&y, x, y, &t);
3262 /* Garbage collection.
3264 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3265 a permutation of them. We are returning one of them; that one must not be
3266 destroyed. We must not destroy 'x_' because the caller owns it. */
3268 gsl_matrix_free (y_);
3270 gsl_matrix_free (t_);
3276 matrix_expr_evaluate_seq (gsl_matrix *start_, gsl_matrix *end_,
3279 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3281 msg (SE, _("All operands of : operator must be scalars."));
3285 long int start = to_scalar (start_);
3286 long int end = to_scalar (end_);
3287 long int by = by_ ? to_scalar (by_) : 1;
3291 msg (SE, _("The increment operand to : must be nonzero."));
3295 long int n = (end >= start && by > 0 ? (end - start + by) / by
3296 : end <= start && by < 0 ? (start - end - by) / -by
3298 gsl_matrix *m = gsl_matrix_alloc (1, n);
3299 for (long int i = 0; i < n; i++)
3300 gsl_matrix_set (m, 0, i, start + i * by);
3305 matrix_expr_evaluate_not (gsl_matrix *a)
3307 for (size_t r = 0; r < a->size1; r++)
3308 for (size_t c = 0; c < a->size2; c++)
3310 double *ae = gsl_matrix_ptr (a, r, c);
3317 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
3319 if (a->size1 != b->size1)
3321 if (!a->size1 || !a->size2)
3323 else if (!b->size1 || !b->size2)
3326 msg (SE, _("All columns in a matrix must have the same number of rows, "
3327 "but this tries to paste matrices with %zu and %zu rows."),
3328 a->size1, b->size1);
3332 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3333 for (size_t y = 0; y < a->size1; y++)
3335 for (size_t x = 0; x < a->size2; x++)
3336 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3337 for (size_t x = 0; x < b->size2; x++)
3338 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3344 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
3346 if (a->size2 != b->size2)
3348 if (!a->size1 || !a->size2)
3350 else if (!b->size1 || !b->size2)
3353 msg (SE, _("All rows in a matrix must have the same number of columns, "
3354 "but this tries to stack matrices with %zu and %zu columns."),
3355 a->size2, b->size2);
3359 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3360 for (size_t x = 0; x < a->size2; x++)
3362 for (size_t y = 0; y < a->size1; y++)
3363 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3364 for (size_t y = 0; y < b->size1; y++)
3365 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3371 matrix_to_vector (gsl_matrix *m)
3374 gsl_vector v = to_vector (m);
3375 assert (v.block == m->block || !v.block);
3379 gsl_matrix_free (m);
3380 return xmemdup (&v, sizeof v);
3394 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3397 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
3398 enum index_type index_type, size_t other_size,
3399 struct index_vector *iv)
3408 msg (SE, _("Vector index must be scalar or vector, not a "
3410 m->size1, m->size2);
3414 msg (SE, _("Matrix row index must be scalar or vector, not a "
3416 m->size1, m->size2);
3420 msg (SE, _("Matrix column index must be scalar or vector, not a "
3422 m->size1, m->size2);
3428 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3429 *iv = (struct index_vector) {
3430 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3433 for (size_t i = 0; i < v.size; i++)
3435 double index = gsl_vector_get (&v, i);
3436 if (index < 1 || index >= size + 1)
3441 msg (SE, _("Index %g is out of range for vector "
3442 "with %zu elements."), index, size);
3446 msg (SE, _("%g is not a valid row index for "
3447 "a %zu×%zu matrix."),
3448 index, size, other_size);
3452 msg (SE, _("%g is not a valid column index for "
3453 "a %zu×%zu matrix."),
3454 index, other_size, size);
3461 iv->indexes[i] = index - 1;
3467 *iv = (struct index_vector) {
3468 .indexes = xnmalloc (size, sizeof *iv->indexes),
3471 for (size_t i = 0; i < size; i++)
3478 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
3480 if (!is_vector (sm))
3482 msg (SE, _("Vector index operator may not be applied to "
3483 "a %zu×%zu matrix."),
3484 sm->size1, sm->size2);
3492 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
3494 if (!matrix_expr_evaluate_vec_all (sm))
3497 gsl_vector sv = to_vector (sm);
3498 struct index_vector iv;
3499 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
3502 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3503 sm->size1 == 1 ? iv.n : 1);
3504 gsl_vector dv = to_vector (dm);
3505 for (size_t dx = 0; dx < iv.n; dx++)
3507 size_t sx = iv.indexes[dx];
3508 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3516 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
3519 struct index_vector iv0;
3520 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
3523 struct index_vector iv1;
3524 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
3531 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3532 for (size_t dy = 0; dy < iv0.n; dy++)
3534 size_t sy = iv0.indexes[dy];
3536 for (size_t dx = 0; dx < iv1.n; dx++)
3538 size_t sx = iv1.indexes[dx];
3539 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3547 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3548 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3549 const struct matrix_function_properties *, gsl_matrix *[], size_t, \
3550 matrix_proto_##PROTO *);
3555 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
3557 if (!is_scalar (subs[index]))
3559 msg (SE, _("Function %s argument %zu must be a scalar, "
3560 "not a %zu×%zu matrix."),
3561 name, index + 1, subs[index]->size1, subs[index]->size2);
3568 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
3570 if (!is_vector (subs[index]))
3572 msg (SE, _("Function %s argument %zu must be a vector, "
3573 "not a %zu×%zu matrix."),
3574 name, index + 1, subs[index]->size1, subs[index]->size2);
3581 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
3583 for (size_t i = 0; i < n_subs; i++)
3585 if (!check_scalar_arg (name, subs, i))
3587 d[i] = to_scalar (subs[i]);
3593 parse_constraint_value (const char **constraintsp)
3596 long retval = strtol (*constraintsp, &tail, 10);
3597 assert (tail > *constraintsp);
3598 *constraintsp = tail;
3603 argument_invalid (const struct matrix_function_properties *props,
3604 size_t arg_index, gsl_matrix *a, size_t y, size_t x,
3608 ds_put_format (out, _("Argument %zu to matrix function %s "
3609 "has invalid value %g."),
3610 arg_index, props->name, gsl_matrix_get (a, y, x));
3612 ds_put_format (out, _("Row %zu, column %zu of argument %zu "
3613 "to matrix function %s has "
3614 "invalid value %g."),
3615 y + 1, x + 1, arg_index, props->name,
3616 gsl_matrix_get (a, y, x));
3620 argument_inequality_error (const struct matrix_function_properties *props,
3621 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3622 size_t b_index, double b,
3623 bool greater, bool equal)
3625 struct string s = DS_EMPTY_INITIALIZER;
3627 argument_invalid (props, a_index, a, y, x, &s);
3628 ds_put_cstr (&s, " ");
3629 if (greater && equal)
3630 ds_put_format (&s, _("This argument must be greater than or "
3631 "equal to argument %zu, but that has value %g."),
3633 else if (greater && !equal)
3634 ds_put_format (&s, _("This argument must be greater than argument %zu, "
3635 "but that has value %g."),
3638 ds_put_format (&s, _("This argument must be less than or "
3639 "equal to argument %zu, but that has value %g."),
3642 ds_put_format (&s, _("This argument must be less than argument %zu, "
3643 "but that has value %g."),
3645 msg (ME, ds_cstr (&s));
3650 argument_value_error (const struct matrix_function_properties *props,
3651 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3653 bool greater, bool equal)
3655 struct string s = DS_EMPTY_INITIALIZER;
3656 argument_invalid (props, a_index, a, y, x, &s);
3657 ds_put_cstr (&s, " ");
3658 if (greater && equal)
3659 ds_put_format (&s, _("This argument must be greater than or equal to %g."),
3661 else if (greater && !equal)
3662 ds_put_format (&s, _("This argument must be greater than %g."), b);
3664 ds_put_format (&s, _("This argument must be less than or equal to %g."), b);
3666 ds_put_format (&s, _("This argument must be less than %g."), b);
3667 msg (ME, ds_cstr (&s));
3672 check_constraints (const struct matrix_function_properties *props,
3673 gsl_matrix *args[], size_t n_args)
3675 const char *constraints = props->constraints;
3679 size_t arg_index = SIZE_MAX;
3680 while (*constraints)
3682 if (*constraints >= 'a' && *constraints <= 'd')
3684 arg_index = *constraints++ - 'a';
3685 assert (arg_index < n_args);
3687 else if (*constraints == '[' || *constraints == '(')
3689 assert (arg_index < n_args);
3690 bool open_lower = *constraints++ == '(';
3691 int minimum = parse_constraint_value (&constraints);
3692 assert (*constraints == ',');
3694 int maximum = parse_constraint_value (&constraints);
3695 assert (*constraints == ']' || *constraints == ')');
3696 bool open_upper = *constraints++ == ')';
3698 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3699 if ((open_lower ? *d <= minimum : *d < minimum)
3700 || (open_upper ? *d >= maximum : *d > maximum))
3702 if (!is_scalar (args[arg_index]))
3703 msg (ME, _("Row %zu, column %zu of argument %zu to matrix "
3704 "function %s has value %g, which is outside "
3705 "the valid range %c%d,%d%c."),
3706 y + 1, x + 1, arg_index + 1, props->name, *d,
3707 open_lower ? '(' : '[',
3709 open_upper ? ')' : ']');
3711 msg (ME, _("Argument %zu to matrix function %s has value %g, "
3712 "which is outside the valid range %c%d,%d%c."),
3713 arg_index + 1, props->name, *d,
3714 open_lower ? '(' : '[',
3716 open_upper ? ')' : ']');
3720 else if (*constraints == '>' || *constraints == '<')
3722 bool greater = *constraints++ == '>';
3724 if (*constraints == '=')
3732 if (*constraints >= 'a' && *constraints <= 'd')
3734 size_t a_index = arg_index;
3735 size_t b_index = *constraints - 'a';
3736 assert (a_index < n_args);
3737 assert (b_index < n_args);
3739 /* We only support one of the two arguments being non-scalar.
3740 It's easier to support only the first one being non-scalar, so
3741 flip things around if it's the other way. */
3742 if (!is_scalar (args[b_index]))
3744 assert (is_scalar (args[a_index]));
3745 size_t tmp_index = a_index;
3747 b_index = tmp_index;
3752 double b = to_scalar (args[b_index]);
3753 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
3755 ? (equal ? !(*a >= b) : !(*a > b))
3756 : (equal ? !(*a <= b) : !(*a < b)))
3758 argument_inequality_error (props,
3759 a_index + 1, args[a_index], y, x,
3767 int comparand = parse_constraint_value (&constraints);
3769 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3771 ? (equal ? !(*d >= comparand) : !(*d > comparand))
3772 : (equal ? !(*d <= comparand) : !(*d < comparand)))
3774 argument_value_error (props,
3775 arg_index + 1, args[arg_index], y, x,
3784 assert (*constraints == ' ');
3786 arg_index = SIZE_MAX;
3793 matrix_expr_evaluate_d_none (
3794 const struct matrix_function_properties *props UNUSED,
3795 gsl_matrix *subs[] UNUSED, size_t n_subs,
3796 matrix_proto_d_none *f)
3798 assert (n_subs == 0);
3800 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3801 gsl_matrix_set (m, 0, 0, f ());
3806 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
3807 gsl_matrix *subs[], size_t n_subs,
3808 matrix_proto_d_d *f)
3810 assert (n_subs == 1);
3813 if (!to_scalar_args (props->name, subs, n_subs, &d))
3816 if (!check_constraints (props, subs, n_subs))
3819 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3820 gsl_matrix_set (m, 0, 0, f (d));
3825 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
3826 gsl_matrix *subs[], size_t n_subs,
3827 matrix_proto_d_dd *f)
3829 assert (n_subs == 2);
3832 if (!to_scalar_args (props->name, subs, n_subs, d))
3835 if (!check_constraints (props, subs, n_subs))
3838 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3839 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
3844 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
3845 gsl_matrix *subs[], size_t n_subs,
3846 matrix_proto_d_ddd *f)
3848 assert (n_subs == 3);
3851 if (!to_scalar_args (props->name, subs, n_subs, d))
3854 if (!check_constraints (props, subs, n_subs))
3857 gsl_matrix *m = gsl_matrix_alloc (1, 1);
3858 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
3863 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
3864 gsl_matrix *subs[], size_t n_subs,
3865 matrix_proto_m_d *f)
3867 assert (n_subs == 1);
3870 return to_scalar_args (props->name, subs, n_subs, &d) ? f (d) : NULL;
3874 matrix_expr_evaluate_m_dd (const struct matrix_function_properties *props,
3875 gsl_matrix *subs[], size_t n_subs,
3876 matrix_proto_m_dd *f)
3878 assert (n_subs == 2);
3881 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1]) : NULL;
3885 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
3886 gsl_matrix *subs[], size_t n_subs,
3887 matrix_proto_m_ddd *f)
3889 assert (n_subs == 3);
3892 return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
3896 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED,
3897 gsl_matrix *subs[], size_t n_subs,
3898 matrix_proto_m_m *f)
3900 assert (n_subs == 1);
3905 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
3906 gsl_matrix *subs[], size_t n_subs,
3907 matrix_proto_m_e *f)
3909 assert (n_subs == 1);
3911 if (!check_constraints (props, subs, n_subs))
3914 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3920 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED,
3921 gsl_matrix *subs[], size_t n_subs,
3922 matrix_proto_m_md *f)
3924 assert (n_subs == 2);
3925 if (!check_scalar_arg (props->name, subs, 1))
3927 return f (subs[0], to_scalar (subs[1]));
3931 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
3932 gsl_matrix *subs[], size_t n_subs,
3933 matrix_proto_m_ed *f)
3935 assert (n_subs == 2);
3936 if (!check_scalar_arg (props->name, subs, 1))
3939 if (!check_constraints (props, subs, n_subs))
3942 double b = to_scalar (subs[1]);
3943 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3949 matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props UNUSED,
3950 gsl_matrix *subs[], size_t n_subs,
3951 matrix_proto_m_mdd *f)
3953 assert (n_subs == 3);
3954 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
3956 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]));
3960 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
3961 gsl_matrix *subs[], size_t n_subs,
3962 matrix_proto_m_edd *f)
3964 assert (n_subs == 3);
3965 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
3968 if (!check_constraints (props, subs, n_subs))
3971 double b = to_scalar (subs[1]);
3972 double c = to_scalar (subs[2]);
3973 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3979 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
3980 gsl_matrix *subs[], size_t n_subs,
3981 matrix_proto_m_eddd *f)
3983 assert (n_subs == 4);
3984 for (size_t i = 1; i < 4; i++)
3985 if (!check_scalar_arg (props->name, subs, i))
3988 if (!check_constraints (props, subs, n_subs))
3991 double b = to_scalar (subs[1]);
3992 double c = to_scalar (subs[2]);
3993 double d = to_scalar (subs[3]);
3994 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
3995 *a = f (*a, b, c, d);
4000 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4001 gsl_matrix *subs[], size_t n_subs,
4002 matrix_proto_m_eed *f)
4004 assert (n_subs == 3);
4005 if (!check_scalar_arg (props->name, subs, 2))
4008 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4009 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4011 msg (ME, _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4012 "%zu×%zu, but %s requires these arguments either to have "
4013 "the same dimensions or for one of them to be a scalar."),
4015 subs[0]->size1, subs[0]->size2,
4016 subs[1]->size1, subs[1]->size2,
4021 if (!check_constraints (props, subs, n_subs))
4024 double c = to_scalar (subs[2]);
4026 if (is_scalar (subs[0]))
4028 double a = to_scalar (subs[0]);
4029 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4035 double b = to_scalar (subs[1]);
4036 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4043 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED,
4044 gsl_matrix *subs[], size_t n_subs,
4045 matrix_proto_m_mm *f)
4047 assert (n_subs == 2);
4048 return f (subs[0], subs[1]);
4052 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4053 gsl_matrix *subs[], size_t n_subs,
4054 matrix_proto_m_v *f)
4056 assert (n_subs == 1);
4057 if (!check_vector_arg (props->name, subs, 0))
4059 gsl_vector v = to_vector (subs[0]);
4064 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED,
4065 gsl_matrix *subs[], size_t n_subs,
4066 matrix_proto_d_m *f)
4068 assert (n_subs == 1);
4070 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4071 gsl_matrix_set (m, 0, 0, f (subs[0]));
4076 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED,
4077 gsl_matrix *subs[], size_t n_subs,
4078 matrix_proto_m_any *f)
4080 return f (subs, n_subs);
4084 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props,
4085 gsl_matrix *subs[], size_t n_subs,
4086 matrix_proto_IDENT *f)
4088 assert (n_subs <= 2);
4091 if (!to_scalar_args (props->name, subs, n_subs, d))
4094 return f (d[0], d[n_subs - 1]);
4098 matrix_expr_evaluate (const struct matrix_expr *e)
4100 if (e->op == MOP_NUMBER)
4102 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4103 gsl_matrix_set (m, 0, 0, e->number);
4106 else if (e->op == MOP_VARIABLE)
4108 const gsl_matrix *src = e->variable->value;
4111 msg (SE, _("Uninitialized variable %s used in expression."),
4116 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4117 gsl_matrix_memcpy (dst, src);
4120 else if (e->op == MOP_EOF)
4122 struct dfm_reader *reader = read_file_open (e->eof);
4123 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4124 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4128 enum { N_LOCAL = 3 };
4129 gsl_matrix *local_subs[N_LOCAL];
4130 gsl_matrix **subs = (e->n_subs < N_LOCAL
4132 : xmalloc (e->n_subs * sizeof *subs));
4134 for (size_t i = 0; i < e->n_subs; i++)
4136 subs[i] = matrix_expr_evaluate (e->subs[i]);
4139 for (size_t j = 0; j < i; j++)
4140 gsl_matrix_free (subs[j]);
4141 if (subs != local_subs)
4147 gsl_matrix *result = NULL;
4150 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4151 case MOP_F_##ENUM: \
4153 static const struct matrix_function_properties props = { \
4155 .constraints = CONSTRAINTS, \
4157 result = matrix_expr_evaluate_##PROTO (&props, subs, e->n_subs, \
4158 matrix_eval_##ENUM); \
4165 gsl_matrix_scale (subs[0], -1.0);
4183 result = matrix_expr_evaluate_elementwise (e->op, subs[0], subs[1]);
4187 result = matrix_expr_evaluate_not (subs[0]);
4191 result = matrix_expr_evaluate_seq (subs[0], subs[1], NULL);
4195 result = matrix_expr_evaluate_seq (subs[0], subs[1], subs[2]);
4199 result = matrix_expr_evaluate_mul_mat (subs[0], subs[1]);
4203 result = matrix_expr_evaluate_exp_mat (subs[0], subs[1]);
4206 case MOP_PASTE_HORZ:
4207 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
4210 case MOP_PASTE_VERT:
4211 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
4215 result = gsl_matrix_alloc (0, 0);
4219 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
4223 result = matrix_expr_evaluate_vec_all (subs[0]);
4227 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
4231 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
4235 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
4244 for (size_t i = 0; i < e->n_subs; i++)
4245 if (subs[i] != result)
4246 gsl_matrix_free (subs[i]);
4247 if (subs != local_subs)
4253 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4256 gsl_matrix *m = matrix_expr_evaluate (e);
4262 msg (SE, _("Expression for %s must evaluate to scalar, "
4263 "not a %zu×%zu matrix."),
4264 context, m->size1, m->size2);
4265 gsl_matrix_free (m);
4270 gsl_matrix_free (m);
4275 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4279 if (!matrix_expr_evaluate_scalar (e, context, &d))
4283 if (d < LONG_MIN || d > LONG_MAX)
4285 msg (SE, _("Expression for %s is outside the integer range."), context);
4292 struct matrix_lvalue
4294 struct matrix_var *var;
4295 struct matrix_expr *indexes[2];
4300 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4304 for (size_t i = 0; i < lvalue->n_indexes; i++)
4305 matrix_expr_destroy (lvalue->indexes[i]);
4310 static struct matrix_lvalue *
4311 matrix_lvalue_parse (struct matrix_state *s)
4313 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4314 if (!lex_force_id (s->lexer))
4316 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4317 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4321 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
4326 lex_get_n (s->lexer, 2);
4328 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
4330 lvalue->n_indexes++;
4332 if (lex_match (s->lexer, T_COMMA))
4334 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
4336 lvalue->n_indexes++;
4338 if (!lex_force_match (s->lexer, T_RPAREN))
4344 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4350 matrix_lvalue_destroy (lvalue);
4355 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4356 enum index_type index_type, size_t other_size,
4357 struct index_vector *iv)
4362 m = matrix_expr_evaluate (e);
4369 bool ok = matrix_normalize_index_vector (m, size, index_type,
4371 gsl_matrix_free (m);
4376 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4377 struct index_vector *iv, gsl_matrix *sm)
4379 gsl_vector dv = to_vector (lvalue->var->value);
4381 /* Convert source matrix 'sm' to source vector 'sv'. */
4382 if (!is_vector (sm))
4384 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
4385 sm->size1, sm->size2);
4388 gsl_vector sv = to_vector (sm);
4389 if (iv->n != sv.size)
4391 msg (SE, _("Can't assign %zu-element vector "
4392 "to %zu-element subvector."), sv.size, iv->n);
4396 /* Assign elements. */
4397 for (size_t x = 0; x < iv->n; x++)
4398 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4403 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4404 struct index_vector *iv0,
4405 struct index_vector *iv1, gsl_matrix *sm)
4407 gsl_matrix *dm = lvalue->var->value;
4409 /* Convert source matrix 'sm' to source vector 'sv'. */
4410 if (iv0->n != sm->size1)
4412 msg (SE, _("Row index vector for assignment to %s has %zu elements "
4413 "but source matrix has %zu rows."),
4414 lvalue->var->name, iv0->n, sm->size1);
4417 if (iv1->n != sm->size2)
4419 msg (SE, _("Column index vector for assignment to %s has %zu elements "
4420 "but source matrix has %zu columns."),
4421 lvalue->var->name, iv1->n, sm->size2);
4425 /* Assign elements. */
4426 for (size_t y = 0; y < iv0->n; y++)
4428 size_t f0 = iv0->indexes[y];
4429 for (size_t x = 0; x < iv1->n; x++)
4431 size_t f1 = iv1->indexes[x];
4432 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4438 /* Takes ownership of M. */
4440 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4441 struct index_vector *iv0, struct index_vector *iv1,
4444 if (!lvalue->n_indexes)
4446 gsl_matrix_free (lvalue->var->value);
4447 lvalue->var->value = sm;
4452 bool ok = (lvalue->n_indexes == 1
4453 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
4454 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
4455 free (iv0->indexes);
4456 free (iv1->indexes);
4457 gsl_matrix_free (sm);
4463 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4464 struct index_vector *iv0,
4465 struct index_vector *iv1)
4467 *iv0 = INDEX_VECTOR_INIT;
4468 *iv1 = INDEX_VECTOR_INIT;
4469 if (!lvalue->n_indexes)
4472 /* Validate destination matrix exists and has the right shape. */
4473 gsl_matrix *dm = lvalue->var->value;
4476 msg (SE, _("Undefined variable %s."), lvalue->var->name);
4479 else if (dm->size1 == 0 || dm->size2 == 0)
4481 msg (SE, _("Cannot index %zu×%zu matrix."),
4482 dm->size1, dm->size2);
4485 else if (lvalue->n_indexes == 1)
4487 if (!is_vector (dm))
4489 msg (SE, _("Can't use vector indexing on %zu×%zu matrix %s."),
4490 dm->size1, dm->size2, lvalue->var->name);
4493 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4494 MAX (dm->size1, dm->size2),
4499 assert (lvalue->n_indexes == 2);
4500 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4501 IV_ROW, dm->size2, iv0))
4504 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4505 IV_COLUMN, dm->size1, iv1))
4507 free (iv0->indexes);
4514 /* Takes ownership of M. */
4516 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
4518 struct index_vector iv0, iv1;
4519 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
4521 gsl_matrix_free (sm);
4525 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
4529 /* Matrix command. */
4533 struct matrix_cmd **commands;
4537 static void matrix_cmds_uninit (struct matrix_cmds *);
4541 enum matrix_cmd_type
4564 struct compute_command
4566 struct matrix_lvalue *lvalue;
4567 struct matrix_expr *rvalue;
4571 struct print_command
4573 struct matrix_expr *expression;
4574 bool use_default_format;
4575 struct fmt_spec format;
4577 int space; /* -1 means NEWPAGE. */
4581 struct string_array literals; /* CLABELS/RLABELS. */
4582 struct matrix_expr *expr; /* CNAMES/RNAMES. */
4588 struct do_if_command
4592 struct matrix_expr *condition;
4593 struct matrix_cmds commands;
4603 struct matrix_var *var;
4604 struct matrix_expr *start;
4605 struct matrix_expr *end;
4606 struct matrix_expr *increment;
4608 /* Loop conditions. */
4609 struct matrix_expr *top_condition;
4610 struct matrix_expr *bottom_condition;
4613 struct matrix_cmds commands;
4617 struct display_command
4619 struct matrix_state *state;
4623 struct release_command
4625 struct matrix_var **vars;
4632 struct matrix_expr *expression;
4633 struct save_file *sf;
4639 struct read_file *rf;
4640 struct matrix_lvalue *dst;
4641 struct matrix_expr *size;
4643 enum fmt_type format;
4650 struct write_command
4652 struct write_file *wf;
4653 struct matrix_expr *expression;
4656 /* If this is nonnull, WRITE uses this format.
4658 If this is NULL, WRITE uses free-field format with as many
4659 digits of precision as needed. */
4660 struct fmt_spec *format;
4669 struct matrix_lvalue *dst;
4670 struct dataset *dataset;
4671 struct file_handle *file;
4673 struct var_syntax *vars;
4675 struct matrix_var *names;
4677 /* Treatment of missing values. */
4682 MGET_ERROR, /* Flag error on command. */
4683 MGET_ACCEPT, /* Accept without change, user-missing only. */
4684 MGET_OMIT, /* Drop this case. */
4685 MGET_RECODE /* Recode to 'substitute'. */
4688 double substitute; /* MGET_RECODE only. */
4694 struct msave_command
4696 struct msave_common *common;
4697 struct matrix_expr *expr;
4698 const char *rowtype;
4699 const struct matrix_expr *factors;
4700 const struct matrix_expr *splits;
4706 struct matrix_state *state;
4707 struct file_handle *file;
4709 struct stringi_set rowtypes;
4713 struct eigen_command
4715 struct matrix_expr *expr;
4716 struct matrix_var *evec;
4717 struct matrix_var *eval;
4721 struct setdiag_command
4723 struct matrix_var *dst;
4724 struct matrix_expr *expr;
4730 struct matrix_expr *expr;
4731 struct matrix_var *u;
4732 struct matrix_var *s;
4733 struct matrix_var *v;
4739 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
4740 static bool matrix_cmd_execute (struct matrix_cmd *);
4741 static void matrix_cmd_destroy (struct matrix_cmd *);
4744 static struct matrix_cmd *
4745 matrix_parse_compute (struct matrix_state *s)
4747 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4748 *cmd = (struct matrix_cmd) {
4749 .type = MCMD_COMPUTE,
4750 .compute = { .lvalue = NULL },
4753 struct compute_command *compute = &cmd->compute;
4754 compute->lvalue = matrix_lvalue_parse (s);
4755 if (!compute->lvalue)
4758 if (!lex_force_match (s->lexer, T_EQUALS))
4761 compute->rvalue = matrix_parse_expr (s);
4762 if (!compute->rvalue)
4768 matrix_cmd_destroy (cmd);
4773 matrix_cmd_execute_compute (struct compute_command *compute)
4775 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
4779 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
4783 print_labels_destroy (struct print_labels *labels)
4787 string_array_destroy (&labels->literals);
4788 matrix_expr_destroy (labels->expr);
4794 parse_literal_print_labels (struct matrix_state *s,
4795 struct print_labels **labelsp)
4797 lex_match (s->lexer, T_EQUALS);
4798 print_labels_destroy (*labelsp);
4799 *labelsp = xzalloc (sizeof **labelsp);
4800 while (lex_token (s->lexer) != T_SLASH
4801 && lex_token (s->lexer) != T_ENDCMD
4802 && lex_token (s->lexer) != T_STOP)
4804 struct string label = DS_EMPTY_INITIALIZER;
4805 while (lex_token (s->lexer) != T_COMMA
4806 && lex_token (s->lexer) != T_SLASH
4807 && lex_token (s->lexer) != T_ENDCMD
4808 && lex_token (s->lexer) != T_STOP)
4810 if (!ds_is_empty (&label))
4811 ds_put_byte (&label, ' ');
4813 if (lex_token (s->lexer) == T_STRING)
4814 ds_put_cstr (&label, lex_tokcstr (s->lexer));
4817 char *rep = lex_next_representation (s->lexer, 0, 0);
4818 ds_put_cstr (&label, rep);
4823 string_array_append_nocopy (&(*labelsp)->literals,
4824 ds_steal_cstr (&label));
4826 if (!lex_match (s->lexer, T_COMMA))
4832 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
4834 lex_match (s->lexer, T_EQUALS);
4835 struct matrix_expr *e = matrix_parse_exp (s);
4839 print_labels_destroy (*labelsp);
4840 *labelsp = xzalloc (sizeof **labelsp);
4841 (*labelsp)->expr = e;
4845 static struct matrix_cmd *
4846 matrix_parse_print (struct matrix_state *s)
4848 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
4849 *cmd = (struct matrix_cmd) {
4852 .use_default_format = true,
4856 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
4859 for (size_t i = 0; ; i++)
4861 enum token_type t = lex_next_token (s->lexer, i);
4862 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
4864 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
4866 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
4869 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
4874 cmd->print.expression = matrix_parse_exp (s);
4875 if (!cmd->print.expression)
4879 while (lex_match (s->lexer, T_SLASH))
4881 if (lex_match_id (s->lexer, "FORMAT"))
4883 lex_match (s->lexer, T_EQUALS);
4884 if (!parse_format_specifier (s->lexer, &cmd->print.format))
4886 cmd->print.use_default_format = false;
4888 else if (lex_match_id (s->lexer, "TITLE"))
4890 lex_match (s->lexer, T_EQUALS);
4891 if (!lex_force_string (s->lexer))
4893 free (cmd->print.title);
4894 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
4897 else if (lex_match_id (s->lexer, "SPACE"))
4899 lex_match (s->lexer, T_EQUALS);
4900 if (lex_match_id (s->lexer, "NEWPAGE"))
4901 cmd->print.space = -1;
4902 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
4904 cmd->print.space = lex_integer (s->lexer);
4910 else if (lex_match_id (s->lexer, "RLABELS"))
4911 parse_literal_print_labels (s, &cmd->print.rlabels);
4912 else if (lex_match_id (s->lexer, "CLABELS"))
4913 parse_literal_print_labels (s, &cmd->print.clabels);
4914 else if (lex_match_id (s->lexer, "RNAMES"))
4916 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
4919 else if (lex_match_id (s->lexer, "CNAMES"))
4921 if (!parse_expr_print_labels (s, &cmd->print.clabels))
4926 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
4927 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
4935 matrix_cmd_destroy (cmd);
4940 matrix_is_integer (const gsl_matrix *m)
4942 for (size_t y = 0; y < m->size1; y++)
4943 for (size_t x = 0; x < m->size2; x++)
4945 double d = gsl_matrix_get (m, y, x);
4953 matrix_max_magnitude (const gsl_matrix *m)
4956 for (size_t y = 0; y < m->size1; y++)
4957 for (size_t x = 0; x < m->size2; x++)
4959 double d = fabs (gsl_matrix_get (m, y, x));
4967 format_fits (struct fmt_spec format, double x)
4969 char *s = data_out (&(union value) { .f = x }, NULL,
4970 &format, settings_get_fmt_settings ());
4971 bool fits = *s != '*' && !strchr (s, 'E');
4976 static struct fmt_spec
4977 default_format (const gsl_matrix *m, int *log_scale)
4981 double max = matrix_max_magnitude (m);
4983 if (matrix_is_integer (m))
4984 for (int w = 1; w <= 12; w++)
4986 struct fmt_spec format = { .type = FMT_F, .w = w };
4987 if (format_fits (format, -max))
4991 if (max >= 1e9 || max <= 1e-4)
4994 snprintf (s, sizeof s, "%.1e", max);
4996 const char *e = strchr (s, 'e');
4998 *log_scale = atoi (e + 1);
5001 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5005 trimmed_string (double d)
5007 struct substring s = ss_buffer ((char *) &d, sizeof d);
5008 ss_rtrim (&s, ss_cstr (" "));
5009 return ss_xstrdup (s);
5012 static struct string_array *
5013 print_labels_get (const struct print_labels *labels, size_t n,
5014 const char *prefix, bool truncate)
5019 struct string_array *out = xzalloc (sizeof *out);
5020 if (labels->literals.n)
5021 string_array_clone (out, &labels->literals);
5022 else if (labels->expr)
5024 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5025 if (m && is_vector (m))
5027 gsl_vector v = to_vector (m);
5028 for (size_t i = 0; i < v.size; i++)
5029 string_array_append_nocopy (out, trimmed_string (
5030 gsl_vector_get (&v, i)));
5032 gsl_matrix_free (m);
5038 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5041 string_array_append (out, "");
5044 string_array_delete (out, out->n - 1);
5047 for (size_t i = 0; i < out->n; i++)
5049 char *s = out->strings[i];
5050 s[strnlen (s, 8)] = '\0';
5057 matrix_cmd_print_space (int space)
5060 output_item_submit (page_break_item_create ());
5061 for (int i = 0; i < space; i++)
5062 output_log ("%s", "");
5066 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
5067 struct fmt_spec format, int log_scale)
5069 matrix_cmd_print_space (print->space);
5071 output_log ("%s", print->title);
5073 output_log (" 10 ** %d X", log_scale);
5075 struct string_array *clabels = print_labels_get (print->clabels,
5076 m->size2, "col", true);
5077 if (clabels && format.w < 8)
5079 struct string_array *rlabels = print_labels_get (print->rlabels,
5080 m->size1, "row", true);
5084 struct string line = DS_EMPTY_INITIALIZER;
5086 ds_put_byte_multiple (&line, ' ', 8);
5087 for (size_t x = 0; x < m->size2; x++)
5088 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5089 output_log_nocopy (ds_steal_cstr (&line));
5092 double scale = pow (10.0, log_scale);
5093 bool numeric = fmt_is_numeric (format.type);
5094 for (size_t y = 0; y < m->size1; y++)
5096 struct string line = DS_EMPTY_INITIALIZER;
5098 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5100 for (size_t x = 0; x < m->size2; x++)
5102 double f = gsl_matrix_get (m, y, x);
5104 ? data_out (&(union value) { .f = f / scale}, NULL,
5105 &format, settings_get_fmt_settings ())
5106 : trimmed_string (f));
5107 ds_put_format (&line, " %s", s);
5110 output_log_nocopy (ds_steal_cstr (&line));
5113 string_array_destroy (rlabels);
5115 string_array_destroy (clabels);
5120 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5121 const struct print_labels *print_labels, size_t n,
5122 const char *name, const char *prefix)
5124 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5126 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5127 for (size_t i = 0; i < n; i++)
5128 pivot_category_create_leaf (
5130 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5131 : pivot_value_new_integer (i + 1)));
5133 d->hide_all_labels = true;
5134 string_array_destroy (labels);
5139 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
5140 struct fmt_spec format, int log_scale)
5142 struct pivot_table *table = pivot_table_create__ (
5143 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5146 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5148 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5149 N_("Columns"), "col");
5151 struct pivot_footnote *footnote = NULL;
5154 char *s = xasprintf ("× 10**%d", log_scale);
5155 footnote = pivot_table_create_footnote (
5156 table, pivot_value_new_user_text_nocopy (s));
5159 double scale = pow (10.0, log_scale);
5160 bool numeric = fmt_is_numeric (format.type);
5161 for (size_t y = 0; y < m->size1; y++)
5162 for (size_t x = 0; x < m->size2; x++)
5164 double f = gsl_matrix_get (m, y, x);
5165 struct pivot_value *v;
5168 v = pivot_value_new_number (f / scale);
5169 v->numeric.format = format;
5172 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5174 pivot_value_add_footnote (v, footnote);
5175 pivot_table_put2 (table, y, x, v);
5178 pivot_table_submit (table);
5182 matrix_cmd_execute_print (const struct print_command *print)
5184 if (print->expression)
5186 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5191 struct fmt_spec format = (print->use_default_format
5192 ? default_format (m, &log_scale)
5195 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5196 matrix_cmd_print_text (print, m, format, log_scale);
5198 matrix_cmd_print_tables (print, m, format, log_scale);
5200 gsl_matrix_free (m);
5204 matrix_cmd_print_space (print->space);
5206 output_log ("%s", print->title);
5213 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
5214 const char *command_name,
5215 const char *stop1, const char *stop2)
5217 lex_end_of_command (s->lexer);
5218 lex_discard_rest_of_command (s->lexer);
5220 size_t allocated = 0;
5223 while (lex_token (s->lexer) == T_ENDCMD)
5226 if (lex_at_phrase (s->lexer, stop1)
5227 || (stop2 && lex_at_phrase (s->lexer, stop2)))
5230 if (lex_at_phrase (s->lexer, "END MATRIX"))
5232 msg (SE, _("Premature END MATRIX within %s."), command_name);
5236 struct matrix_cmd *cmd = matrix_parse_command (s);
5240 if (c->n >= allocated)
5241 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
5242 c->commands[c->n++] = cmd;
5247 matrix_parse_do_if_clause (struct matrix_state *s,
5248 struct do_if_command *ifc,
5249 bool parse_condition,
5250 size_t *allocated_clauses)
5252 if (ifc->n_clauses >= *allocated_clauses)
5253 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5254 sizeof *ifc->clauses);
5255 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5256 *c = (struct do_if_clause) { .condition = NULL };
5258 if (parse_condition)
5260 c->condition = matrix_parse_expr (s);
5265 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
5268 static struct matrix_cmd *
5269 matrix_parse_do_if (struct matrix_state *s)
5271 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5272 *cmd = (struct matrix_cmd) {
5274 .do_if = { .n_clauses = 0 }
5277 size_t allocated_clauses = 0;
5280 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
5283 while (lex_match_phrase (s->lexer, "ELSE IF"));
5285 if (lex_match_id (s->lexer, "ELSE")
5286 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
5289 if (!lex_match_phrase (s->lexer, "END IF"))
5294 matrix_cmd_destroy (cmd);
5299 matrix_cmd_execute_do_if (struct do_if_command *cmd)
5301 for (size_t i = 0; i < cmd->n_clauses; i++)
5303 struct do_if_clause *c = &cmd->clauses[i];
5307 if (!matrix_expr_evaluate_scalar (c->condition,
5308 i ? "ELSE IF" : "DO IF",
5313 for (size_t j = 0; j < c->commands.n; j++)
5314 if (!matrix_cmd_execute (c->commands.commands[j]))
5321 static struct matrix_cmd *
5322 matrix_parse_loop (struct matrix_state *s)
5324 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5325 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5327 struct loop_command *loop = &cmd->loop;
5328 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5330 struct substring name = lex_tokss (s->lexer);
5331 loop->var = matrix_var_lookup (s, name);
5333 loop->var = matrix_var_create (s, name);
5338 loop->start = matrix_parse_expr (s);
5339 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5342 loop->end = matrix_parse_expr (s);
5346 if (lex_match (s->lexer, T_BY))
5348 loop->increment = matrix_parse_expr (s);
5349 if (!loop->increment)
5354 if (lex_match_id (s->lexer, "IF"))
5356 loop->top_condition = matrix_parse_expr (s);
5357 if (!loop->top_condition)
5361 bool was_in_loop = s->in_loop;
5363 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
5365 s->in_loop = was_in_loop;
5369 if (!lex_match_phrase (s->lexer, "END LOOP"))
5372 if (lex_match_id (s->lexer, "IF"))
5374 loop->bottom_condition = matrix_parse_expr (s);
5375 if (!loop->bottom_condition)
5382 matrix_cmd_destroy (cmd);
5387 matrix_cmd_execute_loop (struct loop_command *cmd)
5391 long int increment = 1;
5394 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5395 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5397 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5401 if (increment > 0 ? value > end
5402 : increment < 0 ? value < end
5407 int mxloops = settings_get_mxloops ();
5408 for (int i = 0; i < mxloops; i++)
5413 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5415 gsl_matrix_free (cmd->var->value);
5416 cmd->var->value = NULL;
5418 if (!cmd->var->value)
5419 cmd->var->value = gsl_matrix_alloc (1, 1);
5421 gsl_matrix_set (cmd->var->value, 0, 0, value);
5424 if (cmd->top_condition)
5427 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5432 for (size_t j = 0; j < cmd->commands.n; j++)
5433 if (!matrix_cmd_execute (cmd->commands.commands[j]))
5436 if (cmd->bottom_condition)
5439 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5440 "END LOOP IF", &d) || d > 0)
5447 if (increment > 0 ? value > end : value < end)
5453 static struct matrix_cmd *
5454 matrix_parse_break (struct matrix_state *s)
5458 msg (SE, _("BREAK not inside LOOP."));
5462 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5463 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
5467 static struct matrix_cmd *
5468 matrix_parse_display (struct matrix_state *s)
5470 while (lex_token (s->lexer) != T_ENDCMD)
5472 if (!lex_match_id (s->lexer, "DICTIONARY")
5473 && !lex_match_id (s->lexer, "STATUS"))
5475 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5480 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5481 *cmd = (struct matrix_cmd) { .type = MCMD_DISPLAY, .display = { s } };
5486 compare_matrix_var_pointers (const void *a_, const void *b_)
5488 const struct matrix_var *const *ap = a_;
5489 const struct matrix_var *const *bp = b_;
5490 const struct matrix_var *a = *ap;
5491 const struct matrix_var *b = *bp;
5492 return strcmp (a->name, b->name);
5496 matrix_cmd_execute_display (struct display_command *cmd)
5498 const struct matrix_state *s = cmd->state;
5500 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5501 struct pivot_dimension *attr_dimension
5502 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5503 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5504 N_("Rows"), N_("Columns"));
5505 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5507 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5510 const struct matrix_var *var;
5511 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5513 vars[n_vars++] = var;
5514 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5516 struct pivot_dimension *rows = pivot_dimension_create (
5517 table, PIVOT_AXIS_ROW, N_("Variable"));
5518 for (size_t i = 0; i < n_vars; i++)
5520 const struct matrix_var *var = vars[i];
5521 pivot_category_create_leaf (
5522 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
5524 size_t r = var->value->size1;
5525 size_t c = var->value->size2;
5526 double values[] = { r, c, r * c * sizeof (double) / 1024 };
5527 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
5528 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
5531 pivot_table_submit (table);
5534 static struct matrix_cmd *
5535 matrix_parse_release (struct matrix_state *s)
5537 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5538 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
5540 size_t allocated_vars = 0;
5541 while (lex_token (s->lexer) == T_ID)
5543 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
5546 if (cmd->release.n_vars >= allocated_vars)
5547 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
5548 sizeof *cmd->release.vars);
5549 cmd->release.vars[cmd->release.n_vars++] = var;
5552 lex_error (s->lexer, _("Variable name expected."));
5555 if (!lex_match (s->lexer, T_COMMA))
5563 matrix_cmd_execute_release (struct release_command *cmd)
5565 for (size_t i = 0; i < cmd->n_vars; i++)
5567 struct matrix_var *var = cmd->vars[i];
5568 gsl_matrix_free (var->value);
5573 static struct save_file *
5574 save_file_create (struct matrix_state *s, struct file_handle *fh,
5575 struct string_array *variables,
5576 struct matrix_expr *names,
5577 struct stringi_set *strings)
5579 for (size_t i = 0; i < s->n_save_files; i++)
5581 struct save_file *sf = s->save_files[i];
5582 if (fh_equal (sf->file, fh))
5586 string_array_destroy (variables);
5587 matrix_expr_destroy (names);
5588 stringi_set_destroy (strings);
5594 struct save_file *sf = xmalloc (sizeof *sf);
5595 *sf = (struct save_file) {
5597 .dataset = s->dataset,
5598 .variables = *variables,
5600 .strings = STRINGI_SET_INITIALIZER (sf->strings),
5603 stringi_set_swap (&sf->strings, strings);
5604 stringi_set_destroy (strings);
5606 s->save_files = xrealloc (s->save_files,
5607 (s->n_save_files + 1) * sizeof *s->save_files);
5608 s->save_files[s->n_save_files++] = sf;
5612 static struct casewriter *
5613 save_file_open (struct save_file *sf, gsl_matrix *m)
5615 if (sf->writer || sf->error)
5619 size_t n_variables = caseproto_get_n_widths (
5620 casewriter_get_proto (sf->writer));
5621 if (m->size2 != n_variables)
5623 msg (ME, _("The first SAVE to %s within this matrix program "
5624 "had %zu columns, so a %zu×%zu matrix cannot be "
5626 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
5627 n_variables, m->size1, m->size2);
5635 struct dictionary *dict = dict_create (get_default_encoding ());
5637 /* Fill 'names' with user-specified names if there were any, either from
5638 sf->variables or sf->names. */
5639 struct string_array names = { .n = 0 };
5640 if (sf->variables.n)
5641 string_array_clone (&names, &sf->variables);
5644 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
5645 if (nm && is_vector (nm))
5647 gsl_vector nv = to_vector (nm);
5648 for (size_t i = 0; i < nv.size; i++)
5650 char *name = trimmed_string (gsl_vector_get (&nv, i));
5651 if (dict_id_is_valid (dict, name, true))
5652 string_array_append_nocopy (&names, name);
5657 gsl_matrix_free (nm);
5660 struct stringi_set strings;
5661 stringi_set_clone (&strings, &sf->strings);
5663 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
5668 name = names.strings[i];
5671 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
5675 int width = stringi_set_delete (&strings, name) ? 8 : 0;
5676 struct variable *var = dict_create_var (dict, name, width);
5679 msg (ME, _("Duplicate variable name %s in SAVE statement."),
5685 if (!stringi_set_is_empty (&strings))
5687 size_t count = stringi_set_count (&strings);
5688 const char *example = stringi_set_node_get_string (
5689 stringi_set_first (&strings));
5692 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
5693 "variable %s."), example);
5695 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
5696 "unknown variable: %s.",
5697 "The SAVE command STRINGS subcommand specifies %zu "
5698 "unknown variables, including %s.",
5703 stringi_set_destroy (&strings);
5704 string_array_destroy (&names);
5713 if (sf->file == fh_inline_file ())
5714 sf->writer = autopaging_writer_create (dict_get_proto (dict));
5716 sf->writer = any_writer_open (sf->file, dict);
5728 save_file_destroy (struct save_file *sf)
5732 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
5734 dataset_set_dict (sf->dataset, sf->dict);
5735 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
5739 casewriter_destroy (sf->writer);
5740 dict_unref (sf->dict);
5742 fh_unref (sf->file);
5743 string_array_destroy (&sf->variables);
5744 matrix_expr_destroy (sf->names);
5745 stringi_set_destroy (&sf->strings);
5750 static struct matrix_cmd *
5751 matrix_parse_save (struct matrix_state *s)
5753 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5754 *cmd = (struct matrix_cmd) {
5756 .save = { .expression = NULL },
5759 struct file_handle *fh = NULL;
5760 struct save_command *save = &cmd->save;
5762 struct string_array variables = STRING_ARRAY_INITIALIZER;
5763 struct matrix_expr *names = NULL;
5764 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
5766 save->expression = matrix_parse_exp (s);
5767 if (!save->expression)
5770 while (lex_match (s->lexer, T_SLASH))
5772 if (lex_match_id (s->lexer, "OUTFILE"))
5774 lex_match (s->lexer, T_EQUALS);
5777 fh = (lex_match (s->lexer, T_ASTERISK)
5779 : fh_parse (s->lexer, FH_REF_FILE, s->session));
5783 else if (lex_match_id (s->lexer, "VARIABLES"))
5785 lex_match (s->lexer, T_EQUALS);
5789 struct dictionary *d = dict_create (get_default_encoding ());
5790 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
5791 PV_NO_SCRATCH | PV_NO_DUPLICATE);
5796 string_array_clear (&variables);
5797 variables = (struct string_array) {
5803 else if (lex_match_id (s->lexer, "NAMES"))
5805 lex_match (s->lexer, T_EQUALS);
5806 matrix_expr_destroy (names);
5807 names = matrix_parse_exp (s);
5811 else if (lex_match_id (s->lexer, "STRINGS"))
5813 lex_match (s->lexer, T_EQUALS);
5814 while (lex_token (s->lexer) == T_ID)
5816 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
5818 if (!lex_match (s->lexer, T_COMMA))
5824 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
5832 if (s->prev_save_file)
5833 fh = fh_ref (s->prev_save_file);
5836 lex_sbc_missing ("OUTFILE");
5840 fh_unref (s->prev_save_file);
5841 s->prev_save_file = fh_ref (fh);
5843 if (variables.n && names)
5845 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
5846 matrix_expr_destroy (names);
5850 save->sf = save_file_create (s, fh, &variables, names, &strings);
5854 string_array_destroy (&variables);
5855 matrix_expr_destroy (names);
5856 stringi_set_destroy (&strings);
5858 matrix_cmd_destroy (cmd);
5863 matrix_cmd_execute_save (const struct save_command *save)
5865 gsl_matrix *m = matrix_expr_evaluate (save->expression);
5869 struct casewriter *writer = save_file_open (save->sf, m);
5872 gsl_matrix_free (m);
5876 const struct caseproto *proto = casewriter_get_proto (writer);
5877 for (size_t y = 0; y < m->size1; y++)
5879 struct ccase *c = case_create (proto);
5880 for (size_t x = 0; x < m->size2; x++)
5882 double d = gsl_matrix_get (m, y, x);
5883 int width = caseproto_get_width (proto, x);
5884 union value *value = case_data_rw_idx (c, x);
5888 memcpy (value->s, &d, width);
5890 casewriter_write (writer, c);
5892 gsl_matrix_free (m);
5895 static struct matrix_cmd *
5896 matrix_parse_read (struct matrix_state *s)
5898 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5899 *cmd = (struct matrix_cmd) {
5901 .read = { .format = FMT_F },
5904 struct file_handle *fh = NULL;
5905 char *encoding = NULL;
5906 struct read_command *read = &cmd->read;
5907 read->dst = matrix_lvalue_parse (s);
5912 int repetitions = 0;
5913 int record_width = 0;
5914 bool seen_format = false;
5915 while (lex_match (s->lexer, T_SLASH))
5917 if (lex_match_id (s->lexer, "FILE"))
5919 lex_match (s->lexer, T_EQUALS);
5922 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
5926 else if (lex_match_id (s->lexer, "ENCODING"))
5928 lex_match (s->lexer, T_EQUALS);
5929 if (!lex_force_string (s->lexer))
5933 encoding = ss_xstrdup (lex_tokss (s->lexer));
5937 else if (lex_match_id (s->lexer, "FIELD"))
5939 lex_match (s->lexer, T_EQUALS);
5941 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
5943 read->c1 = lex_integer (s->lexer);
5945 if (!lex_force_match (s->lexer, T_TO)
5946 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
5948 read->c2 = lex_integer (s->lexer) + 1;
5951 record_width = read->c2 - read->c1;
5952 if (lex_match (s->lexer, T_BY))
5954 if (!lex_force_int_range (s->lexer, "BY", 1,
5955 read->c2 - read->c1))
5957 by = lex_integer (s->lexer);
5960 if (record_width % by)
5962 msg (SE, _("BY %d does not evenly divide record width %d."),
5970 else if (lex_match_id (s->lexer, "SIZE"))
5972 lex_match (s->lexer, T_EQUALS);
5973 matrix_expr_destroy (read->size);
5974 read->size = matrix_parse_exp (s);
5978 else if (lex_match_id (s->lexer, "MODE"))
5980 lex_match (s->lexer, T_EQUALS);
5981 if (lex_match_id (s->lexer, "RECTANGULAR"))
5982 read->symmetric = false;
5983 else if (lex_match_id (s->lexer, "SYMMETRIC"))
5984 read->symmetric = true;
5987 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
5991 else if (lex_match_id (s->lexer, "REREAD"))
5992 read->reread = true;
5993 else if (lex_match_id (s->lexer, "FORMAT"))
5997 lex_sbc_only_once ("FORMAT");
6002 lex_match (s->lexer, T_EQUALS);
6004 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6007 const char *p = lex_tokcstr (s->lexer);
6008 if (c_isdigit (p[0]))
6010 repetitions = atoi (p);
6011 p += strspn (p, "0123456789");
6012 if (!fmt_from_name (p, &read->format))
6014 lex_error (s->lexer, _("Unknown format %s."), p);
6019 else if (fmt_from_name (p, &read->format))
6023 struct fmt_spec format;
6024 if (!parse_format_specifier (s->lexer, &format))
6026 read->format = format.type;
6032 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6033 "REREAD", "FORMAT");
6040 lex_sbc_missing ("FIELD");
6044 if (!read->dst->n_indexes && !read->size)
6046 msg (SE, _("SIZE is required for reading data into a full matrix "
6047 "(as opposed to a submatrix)."));
6053 if (s->prev_read_file)
6054 fh = fh_ref (s->prev_read_file);
6057 lex_sbc_missing ("FILE");
6061 fh_unref (s->prev_read_file);
6062 s->prev_read_file = fh_ref (fh);
6064 read->rf = read_file_create (s, fh);
6068 free (read->rf->encoding);
6069 read->rf->encoding = encoding;
6073 /* Field width may be specified in multiple ways:
6076 2. The format on FORMAT.
6077 3. The repetition factor on FORMAT.
6079 (2) and (3) are mutually exclusive.
6081 If more than one of these is present, they must agree. If none of them is
6082 present, then free-field format is used.
6084 if (repetitions > record_width)
6086 msg (SE, _("%d repetitions cannot fit in record width %d."),
6087 repetitions, record_width);
6090 int w = (repetitions ? record_width / repetitions
6096 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6097 "which implies field width %d, "
6098 "but BY specifies field width %d."),
6099 repetitions, record_width, w, by);
6101 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6110 matrix_cmd_destroy (cmd);
6116 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6117 struct substring data, size_t y, size_t x,
6118 int first_column, int last_column, char *error)
6120 int line_number = dfm_get_line_number (reader);
6121 struct msg_location *location = xmalloc (sizeof *location);
6122 *location = (struct msg_location) {
6123 .file_name = xstrdup (dfm_get_file_name (reader)),
6124 .first_line = line_number,
6125 .last_line = line_number + 1,
6126 .first_column = first_column,
6127 .last_column = last_column,
6129 struct msg *m = xmalloc (sizeof *m);
6131 .category = MSG_C_DATA,
6132 .severity = MSG_S_WARNING,
6133 .location = location,
6134 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
6135 "for matrix row %zu, column %zu: %s"),
6136 (int) data.length, data.string, fmt_name (format),
6137 y + 1, x + 1, error),
6145 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
6146 gsl_matrix *m, struct substring p, size_t y, size_t x,
6147 const char *line_start)
6149 const char *input_encoding = dfm_reader_get_encoding (reader);
6152 if (fmt_is_numeric (read->format))
6155 error = data_in (p, input_encoding, read->format,
6156 settings_get_fmt_settings (), &v, 0, NULL);
6157 if (!error && v.f == SYSMIS)
6158 error = xstrdup (_("Matrix data may not contain missing value."));
6163 uint8_t s[sizeof (double)];
6164 union value v = { .s = s };
6165 error = data_in (p, input_encoding, read->format,
6166 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6167 memcpy (&f, s, sizeof f);
6172 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6173 int c2 = c1 + ss_utf8_count_columns (p) - 1;
6174 parse_error (reader, read->format, p, y, x, c1, c2, error);
6178 gsl_matrix_set (m, y, x, f);
6179 if (read->symmetric && x != y)
6180 gsl_matrix_set (m, x, y, f);
6185 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
6186 struct substring *line, const char **startp)
6188 if (dfm_eof (reader))
6190 msg (SE, _("Unexpected end of file reading matrix data."));
6193 dfm_expand_tabs (reader);
6194 struct substring record = dfm_get_record (reader);
6195 /* XXX need to recode record into UTF-8 */
6196 *startp = record.string;
6197 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6202 matrix_read (struct read_command *read, struct dfm_reader *reader,
6205 for (size_t y = 0; y < m->size1; y++)
6207 size_t nx = read->symmetric ? y + 1 : m->size2;
6209 struct substring line = ss_empty ();
6210 const char *line_start = line.string;
6211 for (size_t x = 0; x < nx; x++)
6218 ss_ltrim (&line, ss_cstr (" ,"));
6219 if (!ss_is_empty (line))
6221 if (!matrix_read_line (read, reader, &line, &line_start))
6223 dfm_forward_record (reader);
6226 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6230 if (!matrix_read_line (read, reader, &line, &line_start))
6232 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6233 int f = x % fields_per_line;
6234 if (f == fields_per_line - 1)
6235 dfm_forward_record (reader);
6237 p = ss_substr (line, read->w * f, read->w);
6240 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6244 dfm_forward_record (reader);
6247 ss_ltrim (&line, ss_cstr (" ,"));
6248 if (!ss_is_empty (line))
6251 msg (SW, _("Trailing garbage on line \"%.*s\""),
6252 (int) line.length, line.string);
6259 matrix_cmd_execute_read (struct read_command *read)
6261 struct index_vector iv0, iv1;
6262 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6265 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6268 gsl_matrix *m = matrix_expr_evaluate (read->size);
6274 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6275 "not a %zu×%zu matrix."), m->size1, m->size2);
6276 gsl_matrix_free (m);
6282 gsl_vector v = to_vector (m);
6286 d[0] = gsl_vector_get (&v, 0);
6289 else if (v.size == 2)
6291 d[0] = gsl_vector_get (&v, 0);
6292 d[1] = gsl_vector_get (&v, 1);
6296 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6297 "not a %zu×%zu matrix."),
6298 m->size1, m->size2),
6299 gsl_matrix_free (m);
6304 gsl_matrix_free (m);
6306 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6308 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
6309 "are outside valid range."),
6320 if (read->dst->n_indexes)
6322 size_t submatrix_size[2];
6323 if (read->dst->n_indexes == 2)
6325 submatrix_size[0] = iv0.n;
6326 submatrix_size[1] = iv1.n;
6328 else if (read->dst->var->value->size1 == 1)
6330 submatrix_size[0] = 1;
6331 submatrix_size[1] = iv0.n;
6335 submatrix_size[0] = iv0.n;
6336 submatrix_size[1] = 1;
6341 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6343 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
6344 "differ from submatrix dimensions %zu×%zu."),
6346 submatrix_size[0], submatrix_size[1]);
6354 size[0] = submatrix_size[0];
6355 size[1] = submatrix_size[1];
6359 struct dfm_reader *reader = read_file_open (read->rf);
6361 dfm_reread_record (reader, 1);
6363 if (read->symmetric && size[0] != size[1])
6365 msg (SE, _("Cannot read non-square %zu×%zu matrix "
6366 "using READ with MODE=SYMMETRIC."),
6372 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6373 matrix_read (read, reader, tmp);
6374 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
6377 static struct matrix_cmd *
6378 matrix_parse_write (struct matrix_state *s)
6380 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6381 *cmd = (struct matrix_cmd) {
6385 struct file_handle *fh = NULL;
6386 char *encoding = NULL;
6387 struct write_command *write = &cmd->write;
6388 write->expression = matrix_parse_exp (s);
6389 if (!write->expression)
6393 int repetitions = 0;
6394 int record_width = 0;
6395 enum fmt_type format = FMT_F;
6396 bool has_format = false;
6397 while (lex_match (s->lexer, T_SLASH))
6399 if (lex_match_id (s->lexer, "OUTFILE"))
6401 lex_match (s->lexer, T_EQUALS);
6404 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6408 else if (lex_match_id (s->lexer, "ENCODING"))
6410 lex_match (s->lexer, T_EQUALS);
6411 if (!lex_force_string (s->lexer))
6415 encoding = ss_xstrdup (lex_tokss (s->lexer));
6419 else if (lex_match_id (s->lexer, "FIELD"))
6421 lex_match (s->lexer, T_EQUALS);
6423 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6425 write->c1 = lex_integer (s->lexer);
6427 if (!lex_force_match (s->lexer, T_TO)
6428 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
6430 write->c2 = lex_integer (s->lexer) + 1;
6433 record_width = write->c2 - write->c1;
6434 if (lex_match (s->lexer, T_BY))
6436 if (!lex_force_int_range (s->lexer, "BY", 1,
6437 write->c2 - write->c1))
6439 by = lex_integer (s->lexer);
6442 if (record_width % by)
6444 msg (SE, _("BY %d does not evenly divide record width %d."),
6452 else if (lex_match_id (s->lexer, "MODE"))
6454 lex_match (s->lexer, T_EQUALS);
6455 if (lex_match_id (s->lexer, "RECTANGULAR"))
6456 write->triangular = false;
6457 else if (lex_match_id (s->lexer, "TRIANGULAR"))
6458 write->triangular = true;
6461 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
6465 else if (lex_match_id (s->lexer, "HOLD"))
6467 else if (lex_match_id (s->lexer, "FORMAT"))
6469 if (has_format || write->format)
6471 lex_sbc_only_once ("FORMAT");
6475 lex_match (s->lexer, T_EQUALS);
6477 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6480 const char *p = lex_tokcstr (s->lexer);
6481 if (c_isdigit (p[0]))
6483 repetitions = atoi (p);
6484 p += strspn (p, "0123456789");
6485 if (!fmt_from_name (p, &format))
6487 lex_error (s->lexer, _("Unknown format %s."), p);
6493 else if (fmt_from_name (p, &format))
6500 struct fmt_spec spec;
6501 if (!parse_format_specifier (s->lexer, &spec))
6503 write->format = xmemdup (&spec, sizeof spec);
6508 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
6516 lex_sbc_missing ("FIELD");
6522 if (s->prev_write_file)
6523 fh = fh_ref (s->prev_write_file);
6526 lex_sbc_missing ("OUTFILE");
6530 fh_unref (s->prev_write_file);
6531 s->prev_write_file = fh_ref (fh);
6533 write->wf = write_file_create (s, fh);
6537 free (write->wf->encoding);
6538 write->wf->encoding = encoding;
6542 /* Field width may be specified in multiple ways:
6545 2. The format on FORMAT.
6546 3. The repetition factor on FORMAT.
6548 (2) and (3) are mutually exclusive.
6550 If more than one of these is present, they must agree. If none of them is
6551 present, then free-field format is used.
6553 if (repetitions > record_width)
6555 msg (SE, _("%d repetitions cannot fit in record width %d."),
6556 repetitions, record_width);
6559 int w = (repetitions ? record_width / repetitions
6560 : write->format ? write->format->w
6565 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6566 "which implies field width %d, "
6567 "but BY specifies field width %d."),
6568 repetitions, record_width, w, by);
6570 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6574 if (w && !write->format)
6576 write->format = xmalloc (sizeof *write->format);
6577 *write->format = (struct fmt_spec) { .type = format, .w = w };
6579 if (!fmt_check_output (write->format))
6583 if (write->format && fmt_var_width (write->format) > sizeof (double))
6585 char s[FMT_STRING_LEN_MAX + 1];
6586 fmt_to_string (write->format, s);
6587 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
6588 s, sizeof (double));
6596 matrix_cmd_destroy (cmd);
6601 matrix_cmd_execute_write (struct write_command *write)
6603 gsl_matrix *m = matrix_expr_evaluate (write->expression);
6607 if (write->triangular && m->size1 != m->size2)
6609 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
6610 "the matrix to be written has dimensions %zu×%zu."),
6611 m->size1, m->size2);
6612 gsl_matrix_free (m);
6616 struct dfm_writer *writer = write_file_open (write->wf);
6617 if (!writer || !m->size1)
6619 gsl_matrix_free (m);
6623 const struct fmt_settings *settings = settings_get_fmt_settings ();
6624 struct u8_line *line = write->wf->held;
6625 for (size_t y = 0; y < m->size1; y++)
6629 line = xmalloc (sizeof *line);
6630 u8_line_init (line);
6632 size_t nx = write->triangular ? y + 1 : m->size2;
6634 for (size_t x = 0; x < nx; x++)
6637 double f = gsl_matrix_get (m, y, x);
6641 if (fmt_is_numeric (write->format->type))
6644 v.s = (uint8_t *) &f;
6645 s = data_out (&v, NULL, write->format, settings);
6649 s = xmalloc (DBL_BUFSIZE_BOUND);
6650 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
6651 >= DBL_BUFSIZE_BOUND)
6654 size_t len = strlen (s);
6655 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
6656 if (width + x0 > write->c2)
6658 dfm_put_record_utf8 (writer, line->s.ss.string,
6660 u8_line_clear (line);
6663 u8_line_put (line, x0, x0 + width, s, len);
6666 x0 += write->format ? write->format->w : width + 1;
6669 if (y + 1 >= m->size1 && write->hold)
6671 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
6672 u8_line_clear (line);
6676 u8_line_destroy (line);
6680 write->wf->held = line;
6682 gsl_matrix_free (m);
6685 static struct matrix_cmd *
6686 matrix_parse_get (struct matrix_state *s)
6688 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6689 *cmd = (struct matrix_cmd) {
6692 .dataset = s->dataset,
6693 .user = { .treatment = MGET_ERROR },
6694 .system = { .treatment = MGET_ERROR },
6698 struct get_command *get = &cmd->get;
6699 get->dst = matrix_lvalue_parse (s);
6703 while (lex_match (s->lexer, T_SLASH))
6705 if (lex_match_id (s->lexer, "FILE"))
6707 lex_match (s->lexer, T_EQUALS);
6709 fh_unref (get->file);
6710 if (lex_match (s->lexer, T_ASTERISK))
6714 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
6719 else if (lex_match_id (s->lexer, "ENCODING"))
6721 lex_match (s->lexer, T_EQUALS);
6722 if (!lex_force_string (s->lexer))
6725 free (get->encoding);
6726 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
6730 else if (lex_match_id (s->lexer, "VARIABLES"))
6732 lex_match (s->lexer, T_EQUALS);
6736 lex_sbc_only_once ("VARIABLES");
6740 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
6743 else if (lex_match_id (s->lexer, "NAMES"))
6745 lex_match (s->lexer, T_EQUALS);
6746 if (!lex_force_id (s->lexer))
6749 struct substring name = lex_tokss (s->lexer);
6750 get->names = matrix_var_lookup (s, name);
6752 get->names = matrix_var_create (s, name);
6755 else if (lex_match_id (s->lexer, "MISSING"))
6757 lex_match (s->lexer, T_EQUALS);
6758 if (lex_match_id (s->lexer, "ACCEPT"))
6759 get->user.treatment = MGET_ACCEPT;
6760 else if (lex_match_id (s->lexer, "OMIT"))
6761 get->user.treatment = MGET_OMIT;
6762 else if (lex_is_number (s->lexer))
6764 get->user.treatment = MGET_RECODE;
6765 get->user.substitute = lex_number (s->lexer);
6770 lex_error (s->lexer, NULL);
6774 else if (lex_match_id (s->lexer, "SYSMIS"))
6776 lex_match (s->lexer, T_EQUALS);
6777 if (lex_match_id (s->lexer, "OMIT"))
6778 get->system.treatment = MGET_OMIT;
6779 else if (lex_is_number (s->lexer))
6781 get->system.treatment = MGET_RECODE;
6782 get->system.substitute = lex_number (s->lexer);
6787 lex_error (s->lexer, NULL);
6793 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
6794 "MISSING", "SYSMIS");
6799 if (get->user.treatment != MGET_ACCEPT)
6800 get->system.treatment = MGET_ERROR;
6805 matrix_cmd_destroy (cmd);
6810 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
6811 const struct dictionary *dict)
6813 struct variable **vars;
6818 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
6819 &vars, &n_vars, PV_NUMERIC))
6824 n_vars = dict_get_var_cnt (dict);
6825 vars = xnmalloc (n_vars, sizeof *vars);
6826 for (size_t i = 0; i < n_vars; i++)
6828 struct variable *var = dict_get_var (dict, i);
6829 if (!var_is_numeric (var))
6831 msg (SE, _("GET: Variable %s is not numeric."),
6832 var_get_name (var));
6842 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
6843 for (size_t i = 0; i < n_vars; i++)
6845 char s[sizeof (double)];
6847 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
6848 memcpy (&f, s, sizeof f);
6849 gsl_matrix_set (names, i, 0, f);
6852 gsl_matrix_free (get->names->value);
6853 get->names->value = names;
6857 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
6858 long long int casenum = 1;
6860 for (struct ccase *c = casereader_read (reader); c;
6861 c = casereader_read (reader), casenum++)
6863 if (n_rows >= m->size1)
6865 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
6866 for (size_t y = 0; y < n_rows; y++)
6867 for (size_t x = 0; x < n_vars; x++)
6868 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
6869 gsl_matrix_free (m);
6874 for (size_t x = 0; x < n_vars; x++)
6876 const struct variable *var = vars[x];
6877 double d = case_num (c, var);
6880 if (get->system.treatment == MGET_RECODE)
6881 d = get->system.substitute;
6882 else if (get->system.treatment == MGET_OMIT)
6886 msg (SE, _("GET: Variable %s in case %lld "
6887 "is system-missing."),
6888 var_get_name (var), casenum);
6892 else if (var_is_num_missing (var, d, MV_USER))
6894 if (get->user.treatment == MGET_RECODE)
6895 d = get->user.substitute;
6896 else if (get->user.treatment == MGET_OMIT)
6898 else if (get->user.treatment != MGET_ACCEPT)
6900 msg (SE, _("GET: Variable %s in case %lld has user-missing "
6902 var_get_name (var), casenum, d);
6906 gsl_matrix_set (m, n_rows, x, d);
6917 matrix_lvalue_evaluate_and_assign (get->dst, m);
6920 gsl_matrix_free (m);
6925 matrix_open_casereader (const char *command_name,
6926 struct file_handle *file, const char *encoding,
6927 struct dataset *dataset,
6928 struct casereader **readerp, struct dictionary **dictp)
6932 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
6933 return *readerp != NULL;
6937 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
6939 msg (ME, _("The %s command cannot read an empty active file."),
6943 *readerp = proc_open (dataset);
6944 *dictp = dict_ref (dataset_dict (dataset));
6950 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
6951 struct casereader *reader, struct dictionary *dict)
6954 casereader_destroy (reader);
6956 proc_commit (dataset);
6960 matrix_cmd_execute_get (struct get_command *get)
6962 struct casereader *r;
6963 struct dictionary *d;
6964 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
6967 matrix_cmd_execute_get__ (get, r, d);
6968 matrix_close_casereader (get->file, get->dataset, r, d);
6973 match_rowtype (struct lexer *lexer)
6975 static const char *rowtypes[] = {
6976 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
6978 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
6980 for (size_t i = 0; i < n_rowtypes; i++)
6981 if (lex_match_id (lexer, rowtypes[i]))
6984 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
6989 parse_var_names (struct lexer *lexer, struct string_array *sa)
6991 lex_match (lexer, T_EQUALS);
6993 string_array_clear (sa);
6995 struct dictionary *dict = dict_create (get_default_encoding ());
6998 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
6999 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7004 for (size_t i = 0; i < n_names; i++)
7005 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7006 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7008 msg (SE, _("Variable name %s is reserved."), names[i]);
7009 for (size_t j = 0; j < n_names; j++)
7015 string_array_clear (sa);
7016 sa->strings = names;
7017 sa->n = sa->allocated = n_names;
7023 msave_common_destroy (struct msave_common *common)
7027 fh_unref (common->outfile);
7028 string_array_destroy (&common->variables);
7029 string_array_destroy (&common->fnames);
7030 string_array_destroy (&common->snames);
7032 for (size_t i = 0; i < common->n_factors; i++)
7033 matrix_expr_destroy (common->factors[i]);
7034 free (common->factors);
7036 for (size_t i = 0; i < common->n_splits; i++)
7037 matrix_expr_destroy (common->splits[i]);
7038 free (common->splits);
7040 dict_unref (common->dict);
7041 casewriter_destroy (common->writer);
7048 compare_variables (const char *keyword,
7049 const struct string_array *new,
7050 const struct string_array *old)
7056 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7057 "on the first MSAVE within MATRIX."), keyword);
7060 else if (!string_array_equal_case (old, new))
7062 msg (SE, _("%s must specify the same variables each time within "
7063 "a given MATRIX."), keyword);
7070 static struct matrix_cmd *
7071 matrix_parse_msave (struct matrix_state *s)
7073 struct msave_common *common = xmalloc (sizeof *common);
7074 *common = (struct msave_common) { .outfile = NULL };
7076 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7077 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7079 struct matrix_expr *splits = NULL;
7080 struct matrix_expr *factors = NULL;
7082 struct msave_command *msave = &cmd->msave;
7083 msave->expr = matrix_parse_exp (s);
7087 while (lex_match (s->lexer, T_SLASH))
7089 if (lex_match_id (s->lexer, "TYPE"))
7091 lex_match (s->lexer, T_EQUALS);
7093 msave->rowtype = match_rowtype (s->lexer);
7094 if (!msave->rowtype)
7097 else if (lex_match_id (s->lexer, "OUTFILE"))
7099 lex_match (s->lexer, T_EQUALS);
7101 fh_unref (common->outfile);
7102 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7103 if (!common->outfile)
7106 else if (lex_match_id (s->lexer, "VARIABLES"))
7108 if (!parse_var_names (s->lexer, &common->variables))
7111 else if (lex_match_id (s->lexer, "FNAMES"))
7113 if (!parse_var_names (s->lexer, &common->fnames))
7116 else if (lex_match_id (s->lexer, "SNAMES"))
7118 if (!parse_var_names (s->lexer, &common->snames))
7121 else if (lex_match_id (s->lexer, "SPLIT"))
7123 lex_match (s->lexer, T_EQUALS);
7125 matrix_expr_destroy (splits);
7126 splits = matrix_parse_exp (s);
7130 else if (lex_match_id (s->lexer, "FACTOR"))
7132 lex_match (s->lexer, T_EQUALS);
7134 matrix_expr_destroy (factors);
7135 factors = matrix_parse_exp (s);
7141 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7142 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7146 if (!msave->rowtype)
7148 lex_sbc_missing ("TYPE");
7154 if (common->fnames.n && !factors)
7156 msg (SE, _("FNAMES requires FACTOR."));
7159 if (common->snames.n && !splits)
7161 msg (SE, _("SNAMES requires SPLIT."));
7164 if (!common->outfile)
7166 lex_sbc_missing ("OUTFILE");
7173 if (common->outfile)
7175 if (!fh_equal (common->outfile, s->common->outfile))
7177 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7178 "within a single MATRIX command."));
7182 if (!compare_variables ("VARIABLES",
7183 &common->variables, &s->common->variables)
7184 || !compare_variables ("FNAMES", &common->fnames, &s->common->fnames)
7185 || !compare_variables ("SNAMES", &common->snames, &s->common->snames))
7187 msave_common_destroy (common);
7189 msave->common = s->common;
7191 struct msave_common *c = s->common;
7194 if (c->n_factors >= c->allocated_factors)
7195 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7196 sizeof *c->factors);
7197 c->factors[c->n_factors++] = factors;
7199 if (c->n_factors > 0)
7200 msave->factors = c->factors[c->n_factors - 1];
7204 if (c->n_splits >= c->allocated_splits)
7205 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7207 c->splits[c->n_splits++] = splits;
7209 if (c->n_splits > 0)
7210 msave->splits = c->splits[c->n_splits - 1];
7215 matrix_expr_destroy (splits);
7216 matrix_expr_destroy (factors);
7217 msave_common_destroy (common);
7218 matrix_cmd_destroy (cmd);
7223 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7225 gsl_matrix *m = matrix_expr_evaluate (e);
7231 msg (SE, _("%s expression must evaluate to vector, "
7232 "not a %zu×%zu matrix."),
7233 name, m->size1, m->size2);
7234 gsl_matrix_free (m);
7238 return matrix_to_vector (m);
7242 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7244 for (size_t i = 0; i < vars->n; i++)
7245 if (!dict_create_var (d, vars->strings[i], 0))
7246 return vars->strings[i];
7250 static struct dictionary *
7251 msave_create_dict (const struct msave_common *common)
7253 struct dictionary *dict = dict_create (get_default_encoding ());
7255 const char *dup_split = msave_add_vars (dict, &common->snames);
7258 /* Should not be possible because the parser ensures that the names are
7263 dict_create_var_assert (dict, "ROWTYPE_", 8);
7265 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7268 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
7272 dict_create_var_assert (dict, "VARNAME_", 8);
7274 const char *dup_var = msave_add_vars (dict, &common->variables);
7277 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
7289 matrix_cmd_execute_msave (struct msave_command *msave)
7291 struct msave_common *common = msave->common;
7292 gsl_matrix *m = NULL;
7293 gsl_vector *factors = NULL;
7294 gsl_vector *splits = NULL;
7296 m = matrix_expr_evaluate (msave->expr);
7300 if (!common->variables.n)
7301 for (size_t i = 0; i < m->size2; i++)
7302 string_array_append_nocopy (&common->variables,
7303 xasprintf ("COL%zu", i + 1));
7304 else if (m->size2 != common->variables.n)
7307 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7308 m->size2, common->variables.n);
7314 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7318 if (!common->fnames.n)
7319 for (size_t i = 0; i < factors->size; i++)
7320 string_array_append_nocopy (&common->fnames,
7321 xasprintf ("FAC%zu", i + 1));
7322 else if (factors->size != common->fnames.n)
7324 msg (SE, _("There are %zu factor variables, "
7325 "but %zu factor values were supplied."),
7326 common->fnames.n, factors->size);
7333 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7337 if (!common->snames.n)
7338 for (size_t i = 0; i < splits->size; i++)
7339 string_array_append_nocopy (&common->snames,
7340 xasprintf ("SPL%zu", i + 1));
7341 else if (splits->size != common->snames.n)
7343 msg (SE, _("There are %zu split variables, "
7344 "but %zu split values were supplied."),
7345 common->snames.n, splits->size);
7350 if (!common->writer)
7352 struct dictionary *dict = msave_create_dict (common);
7356 common->writer = any_writer_open (common->outfile, dict);
7357 if (!common->writer)
7363 common->dict = dict;
7366 bool matrix = (!strcmp (msave->rowtype, "COV")
7367 || !strcmp (msave->rowtype, "CORR"));
7368 for (size_t y = 0; y < m->size1; y++)
7370 struct ccase *c = case_create (dict_get_proto (common->dict));
7373 /* Split variables */
7375 for (size_t i = 0; i < splits->size; i++)
7376 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
7379 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7380 msave->rowtype, ' ');
7384 for (size_t i = 0; i < factors->size; i++)
7385 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
7388 const char *varname_ = (matrix && y < common->variables.n
7389 ? common->variables.strings[y]
7391 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7394 /* Continuous variables. */
7395 for (size_t x = 0; x < m->size2; x++)
7396 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
7397 casewriter_write (common->writer, c);
7401 gsl_matrix_free (m);
7402 gsl_vector_free (factors);
7403 gsl_vector_free (splits);
7406 static struct matrix_cmd *
7407 matrix_parse_mget (struct matrix_state *s)
7409 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7410 *cmd = (struct matrix_cmd) {
7414 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
7418 struct mget_command *mget = &cmd->mget;
7420 lex_match (s->lexer, T_SLASH);
7421 while (lex_token (s->lexer) != T_ENDCMD)
7423 if (lex_match_id (s->lexer, "FILE"))
7425 lex_match (s->lexer, T_EQUALS);
7427 fh_unref (mget->file);
7428 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7432 else if (lex_match_id (s->lexer, "ENCODING"))
7434 lex_match (s->lexer, T_EQUALS);
7435 if (!lex_force_string (s->lexer))
7438 free (mget->encoding);
7439 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
7443 else if (lex_match_id (s->lexer, "TYPE"))
7445 lex_match (s->lexer, T_EQUALS);
7446 while (lex_token (s->lexer) != T_SLASH
7447 && lex_token (s->lexer) != T_ENDCMD)
7449 const char *rowtype = match_rowtype (s->lexer);
7453 stringi_set_insert (&mget->rowtypes, rowtype);
7458 lex_error_expecting (s->lexer, "FILE", "TYPE");
7461 lex_match (s->lexer, T_SLASH);
7466 matrix_cmd_destroy (cmd);
7470 static const struct variable *
7471 get_a8_var (const struct dictionary *d, const char *name)
7473 const struct variable *v = dict_lookup_var (d, name);
7476 msg (SE, _("Matrix data file lacks %s variable."), name);
7479 if (var_get_width (v) != 8)
7481 msg (SE, _("%s variable in matrix data file must be "
7482 "8-byte string, but it has width %d."),
7483 name, var_get_width (v));
7490 var_changed (const struct ccase *ca, const struct ccase *cb,
7491 const struct variable *var)
7494 ? !value_equal (case_data (ca, var), case_data (cb, var),
7495 var_get_width (var))
7500 vars_changed (const struct ccase *ca, const struct ccase *cb,
7501 const struct dictionary *d,
7502 size_t first_var, size_t n_vars)
7504 for (size_t i = 0; i < n_vars; i++)
7506 const struct variable *v = dict_get_var (d, first_var + i);
7507 if (var_changed (ca, cb, v))
7514 vars_all_missing (const struct ccase *c, const struct dictionary *d,
7515 size_t first_var, size_t n_vars)
7517 for (size_t i = 0; i < n_vars; i++)
7518 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
7524 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
7525 const struct dictionary *d,
7526 const struct variable *rowtype_var,
7527 const struct stringi_set *accepted_rowtypes,
7528 struct matrix_state *s,
7529 size_t ss, size_t sn, size_t si,
7530 size_t fs, size_t fn, size_t fi,
7531 size_t cs, size_t cn,
7532 struct pivot_table *pt,
7533 struct pivot_dimension *var_dimension)
7538 /* Is this a matrix for pooled data, either where there are no factor
7539 variables or the factor variables are missing? */
7540 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
7542 struct substring rowtype = case_ss (rows[0], rowtype_var);
7543 ss_rtrim (&rowtype, ss_cstr (" "));
7544 if (!stringi_set_is_empty (accepted_rowtypes)
7545 && !stringi_set_contains_len (accepted_rowtypes,
7546 rowtype.string, rowtype.length))
7549 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
7550 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
7551 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
7552 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
7553 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
7554 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
7558 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
7559 (int) rowtype.length, rowtype.string);
7563 struct string name = DS_EMPTY_INITIALIZER;
7564 ds_put_cstr (&name, prefix);
7566 ds_put_format (&name, "F%zu", fi);
7568 ds_put_format (&name, "S%zu", si);
7570 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
7572 mv = matrix_var_create (s, ds_ss (&name));
7575 msg (SW, _("Matrix data file contains variable with existing name %s."),
7577 goto exit_free_name;
7580 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
7581 size_t n_missing = 0;
7582 for (size_t y = 0; y < n_rows; y++)
7584 for (size_t x = 0; x < cn; x++)
7586 struct variable *var = dict_get_var (d, cs + x);
7587 double value = case_num (rows[y], var);
7588 if (var_is_num_missing (var, value, MV_ANY))
7593 gsl_matrix_set (m, y, x, value);
7597 int var_index = pivot_category_create_leaf (
7598 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
7599 double values[] = { n_rows, cn };
7600 for (size_t j = 0; j < sn; j++)
7602 struct variable *var = dict_get_var (d, ss + j);
7603 const union value *value = case_data (rows[0], var);
7604 pivot_table_put2 (pt, j, var_index,
7605 pivot_value_new_var_value (var, value));
7607 for (size_t j = 0; j < fn; j++)
7609 struct variable *var = dict_get_var (d, fs + j);
7610 const union value sysmis = { .f = SYSMIS };
7611 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
7612 pivot_table_put2 (pt, j + sn, var_index,
7613 pivot_value_new_var_value (var, value));
7615 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
7616 pivot_table_put2 (pt, j + sn + fn, var_index,
7617 pivot_value_new_integer (values[j]));
7620 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
7621 "value, which was treated as zero.",
7622 "Matrix data file variable %s contains %zu missing "
7623 "values, which were treated as zero.", n_missing),
7624 ds_cstr (&name), n_missing);
7631 for (size_t y = 0; y < n_rows; y++)
7632 case_unref (rows[y]);
7636 matrix_cmd_execute_mget__ (struct mget_command *mget,
7637 struct casereader *r, const struct dictionary *d)
7639 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
7640 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
7641 if (!rowtype_ || !varname_)
7644 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
7646 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
7649 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
7651 msg (SE, _("Matrix data file contains no continuous variables."));
7655 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
7657 const struct variable *v = dict_get_var (d, i);
7658 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
7661 _("Matrix data file contains unexpected string variable %s."),
7667 /* SPLIT variables. */
7669 size_t sn = var_get_dict_index (rowtype_);
7670 struct ccase *sc = NULL;
7673 /* FACTOR variables. */
7674 size_t fs = var_get_dict_index (rowtype_) + 1;
7675 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
7676 struct ccase *fc = NULL;
7679 /* Continuous variables. */
7680 size_t cs = var_get_dict_index (varname_) + 1;
7681 size_t cn = dict_get_var_cnt (d) - cs;
7682 struct ccase *cc = NULL;
7685 struct pivot_table *pt = pivot_table_create (
7686 N_("Matrix Variables Created by MGET"));
7687 struct pivot_dimension *attr_dimension = pivot_dimension_create (
7688 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
7689 struct pivot_dimension *var_dimension = pivot_dimension_create (
7690 pt, PIVOT_AXIS_ROW, N_("Variable"));
7693 struct pivot_category *splits = pivot_category_create_group (
7694 attr_dimension->root, N_("Split Values"));
7695 for (size_t i = 0; i < sn; i++)
7696 pivot_category_create_leaf (splits, pivot_value_new_variable (
7697 dict_get_var (d, ss + i)));
7701 struct pivot_category *factors = pivot_category_create_group (
7702 attr_dimension->root, N_("Factors"));
7703 for (size_t i = 0; i < fn; i++)
7704 pivot_category_create_leaf (factors, pivot_value_new_variable (
7705 dict_get_var (d, fs + i)));
7707 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
7708 N_("Rows"), N_("Columns"));
7711 struct ccase **rows = NULL;
7712 size_t allocated_rows = 0;
7716 while ((c = casereader_read (r)) != NULL)
7718 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
7728 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
7729 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
7730 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
7733 if (change != NOTHING_CHANGED)
7735 matrix_mget_commit_var (rows, n_rows, d,
7736 rowtype_, &mget->rowtypes,
7747 if (n_rows >= allocated_rows)
7748 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
7751 if (change == SPLITS_CHANGED)
7757 /* Reset the factor number, if there are factors. */
7761 if (row_has_factors)
7767 else if (change == FACTORS_CHANGED)
7769 if (row_has_factors)
7775 matrix_mget_commit_var (rows, n_rows, d,
7776 rowtype_, &mget->rowtypes,
7788 if (var_dimension->n_leaves)
7789 pivot_table_submit (pt);
7791 pivot_table_unref (pt);
7795 matrix_cmd_execute_mget (struct mget_command *mget)
7797 struct casereader *r;
7798 struct dictionary *d;
7799 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
7800 mget->state->dataset, &r, &d))
7802 matrix_cmd_execute_mget__ (mget, r, d);
7803 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
7808 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
7810 if (!lex_force_id (s->lexer))
7813 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
7815 *varp = matrix_var_create (s, lex_tokss (s->lexer));
7820 static struct matrix_cmd *
7821 matrix_parse_eigen (struct matrix_state *s)
7823 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7824 *cmd = (struct matrix_cmd) {
7826 .eigen = { .expr = NULL }
7829 struct eigen_command *eigen = &cmd->eigen;
7830 if (!lex_force_match (s->lexer, T_LPAREN))
7832 eigen->expr = matrix_parse_expr (s);
7834 || !lex_force_match (s->lexer, T_COMMA)
7835 || !matrix_parse_dst_var (s, &eigen->evec)
7836 || !lex_force_match (s->lexer, T_COMMA)
7837 || !matrix_parse_dst_var (s, &eigen->eval)
7838 || !lex_force_match (s->lexer, T_RPAREN))
7844 matrix_cmd_destroy (cmd);
7849 matrix_cmd_execute_eigen (struct eigen_command *eigen)
7851 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
7854 if (!is_symmetric (A))
7856 msg (SE, _("Argument of EIGEN must be symmetric."));
7857 gsl_matrix_free (A);
7861 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
7862 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
7863 gsl_vector v_eval = to_vector (eval);
7864 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
7865 gsl_eigen_symmv (A, &v_eval, evec, w);
7866 gsl_eigen_symmv_free (w);
7868 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
7870 gsl_matrix_free (A);
7872 gsl_matrix_free (eigen->eval->value);
7873 eigen->eval->value = eval;
7875 gsl_matrix_free (eigen->evec->value);
7876 eigen->evec->value = evec;
7879 static struct matrix_cmd *
7880 matrix_parse_setdiag (struct matrix_state *s)
7882 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7883 *cmd = (struct matrix_cmd) {
7884 .type = MCMD_SETDIAG,
7885 .setdiag = { .dst = NULL }
7888 struct setdiag_command *setdiag = &cmd->setdiag;
7889 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
7892 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
7895 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
7900 if (!lex_force_match (s->lexer, T_COMMA))
7903 setdiag->expr = matrix_parse_expr (s);
7907 if (!lex_force_match (s->lexer, T_RPAREN))
7913 matrix_cmd_destroy (cmd);
7918 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
7920 gsl_matrix *dst = setdiag->dst->value;
7923 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
7924 setdiag->dst->name);
7928 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
7932 size_t n = MIN (dst->size1, dst->size2);
7933 if (is_scalar (src))
7935 double d = to_scalar (src);
7936 for (size_t i = 0; i < n; i++)
7937 gsl_matrix_set (dst, i, i, d);
7939 else if (is_vector (src))
7941 gsl_vector v = to_vector (src);
7942 for (size_t i = 0; i < n && i < v.size; i++)
7943 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
7946 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
7947 "not a %zu×%zu matrix."),
7948 src->size1, src->size2);
7949 gsl_matrix_free (src);
7952 static struct matrix_cmd *
7953 matrix_parse_svd (struct matrix_state *s)
7955 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7956 *cmd = (struct matrix_cmd) {
7958 .svd = { .expr = NULL }
7961 struct svd_command *svd = &cmd->svd;
7962 if (!lex_force_match (s->lexer, T_LPAREN))
7964 svd->expr = matrix_parse_expr (s);
7966 || !lex_force_match (s->lexer, T_COMMA)
7967 || !matrix_parse_dst_var (s, &svd->u)
7968 || !lex_force_match (s->lexer, T_COMMA)
7969 || !matrix_parse_dst_var (s, &svd->s)
7970 || !lex_force_match (s->lexer, T_COMMA)
7971 || !matrix_parse_dst_var (s, &svd->v)
7972 || !lex_force_match (s->lexer, T_RPAREN))
7978 matrix_cmd_destroy (cmd);
7983 matrix_cmd_execute_svd (struct svd_command *svd)
7985 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
7989 if (m->size1 >= m->size2)
7992 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
7993 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
7994 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
7995 gsl_vector *work = gsl_vector_alloc (A->size2);
7996 gsl_linalg_SV_decomp (A, V, &Sv, work);
7997 gsl_vector_free (work);
7999 matrix_var_set (svd->u, A);
8000 matrix_var_set (svd->s, S);
8001 matrix_var_set (svd->v, V);
8005 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8006 gsl_matrix_transpose_memcpy (At, m);
8007 gsl_matrix_free (m);
8009 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8010 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8011 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8012 gsl_vector *work = gsl_vector_alloc (At->size2);
8013 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8014 gsl_vector_free (work);
8016 matrix_var_set (svd->v, At);
8017 matrix_var_set (svd->s, St);
8018 matrix_var_set (svd->u, Vt);
8023 matrix_cmd_execute (struct matrix_cmd *cmd)
8028 matrix_cmd_execute_compute (&cmd->compute);
8032 matrix_cmd_execute_print (&cmd->print);
8036 return matrix_cmd_execute_do_if (&cmd->do_if);
8039 matrix_cmd_execute_loop (&cmd->loop);
8046 matrix_cmd_execute_display (&cmd->display);
8050 matrix_cmd_execute_release (&cmd->release);
8054 matrix_cmd_execute_save (&cmd->save);
8058 matrix_cmd_execute_read (&cmd->read);
8062 matrix_cmd_execute_write (&cmd->write);
8066 matrix_cmd_execute_get (&cmd->get);
8070 matrix_cmd_execute_msave (&cmd->msave);
8074 matrix_cmd_execute_mget (&cmd->mget);
8078 matrix_cmd_execute_eigen (&cmd->eigen);
8082 matrix_cmd_execute_setdiag (&cmd->setdiag);
8086 matrix_cmd_execute_svd (&cmd->svd);
8094 matrix_cmds_uninit (struct matrix_cmds *cmds)
8096 for (size_t i = 0; i < cmds->n; i++)
8097 matrix_cmd_destroy (cmds->commands[i]);
8098 free (cmds->commands);
8102 matrix_cmd_destroy (struct matrix_cmd *cmd)
8110 matrix_lvalue_destroy (cmd->compute.lvalue);
8111 matrix_expr_destroy (cmd->compute.rvalue);
8115 matrix_expr_destroy (cmd->print.expression);
8116 free (cmd->print.title);
8117 print_labels_destroy (cmd->print.rlabels);
8118 print_labels_destroy (cmd->print.clabels);
8122 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8124 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8125 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
8127 free (cmd->do_if.clauses);
8131 matrix_expr_destroy (cmd->loop.start);
8132 matrix_expr_destroy (cmd->loop.end);
8133 matrix_expr_destroy (cmd->loop.increment);
8134 matrix_expr_destroy (cmd->loop.top_condition);
8135 matrix_expr_destroy (cmd->loop.bottom_condition);
8136 matrix_cmds_uninit (&cmd->loop.commands);
8146 free (cmd->release.vars);
8150 matrix_expr_destroy (cmd->save.expression);
8154 matrix_lvalue_destroy (cmd->read.dst);
8155 matrix_expr_destroy (cmd->read.size);
8159 matrix_expr_destroy (cmd->write.expression);
8160 free (cmd->write.format);
8164 matrix_lvalue_destroy (cmd->get.dst);
8165 fh_unref (cmd->get.file);
8166 free (cmd->get.encoding);
8167 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8171 matrix_expr_destroy (cmd->msave.expr);
8175 fh_unref (cmd->mget.file);
8176 stringi_set_destroy (&cmd->mget.rowtypes);
8180 matrix_expr_destroy (cmd->eigen.expr);
8184 matrix_expr_destroy (cmd->setdiag.expr);
8188 matrix_expr_destroy (cmd->svd.expr);
8194 struct matrix_command_name
8197 struct matrix_cmd *(*parse) (struct matrix_state *);
8200 static const struct matrix_command_name *
8201 matrix_parse_command_name (struct lexer *lexer)
8203 static const struct matrix_command_name commands[] = {
8204 { "COMPUTE", matrix_parse_compute },
8205 { "CALL EIGEN", matrix_parse_eigen },
8206 { "CALL SETDIAG", matrix_parse_setdiag },
8207 { "CALL SVD", matrix_parse_svd },
8208 { "PRINT", matrix_parse_print },
8209 { "DO IF", matrix_parse_do_if },
8210 { "LOOP", matrix_parse_loop },
8211 { "BREAK", matrix_parse_break },
8212 { "READ", matrix_parse_read },
8213 { "WRITE", matrix_parse_write },
8214 { "GET", matrix_parse_get },
8215 { "SAVE", matrix_parse_save },
8216 { "MGET", matrix_parse_mget },
8217 { "MSAVE", matrix_parse_msave },
8218 { "DISPLAY", matrix_parse_display },
8219 { "RELEASE", matrix_parse_release },
8221 static size_t n = sizeof commands / sizeof *commands;
8223 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8224 if (lex_match_phrase (lexer, c->name))
8229 static struct matrix_cmd *
8230 matrix_parse_command (struct matrix_state *s)
8232 size_t nesting_level = SIZE_MAX;
8234 struct matrix_cmd *c = NULL;
8235 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
8237 lex_error (s->lexer, _("Unknown matrix command."));
8238 else if (!cmd->parse)
8239 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8243 nesting_level = output_open_group (
8244 group_item_create_nocopy (utf8_to_title (cmd->name),
8245 utf8_to_title (cmd->name)));
8250 lex_end_of_command (s->lexer);
8251 lex_discard_rest_of_command (s->lexer);
8252 if (nesting_level != SIZE_MAX)
8253 output_close_groups (nesting_level);
8259 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8261 if (!lex_force_match (lexer, T_ENDCMD))
8264 struct matrix_state state = {
8266 .session = dataset_session (ds),
8268 .vars = HMAP_INITIALIZER (state.vars),
8273 while (lex_match (lexer, T_ENDCMD))
8275 if (lex_token (lexer) == T_STOP)
8277 msg (SE, _("Unexpected end of input expecting matrix command."));
8281 if (lex_match_phrase (lexer, "END MATRIX"))
8284 struct matrix_cmd *c = matrix_parse_command (&state);
8287 matrix_cmd_execute (c);
8288 matrix_cmd_destroy (c);
8292 struct matrix_var *var, *next;
8293 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8296 gsl_matrix_free (var->value);
8297 hmap_delete (&state.vars, &var->hmap_node);
8300 hmap_destroy (&state.vars);
8301 msave_common_destroy (state.common);
8302 fh_unref (state.prev_read_file);
8303 for (size_t i = 0; i < state.n_read_files; i++)
8304 read_file_destroy (state.read_files[i]);
8305 free (state.read_files);
8306 fh_unref (state.prev_write_file);
8307 for (size_t i = 0; i < state.n_write_files; i++)
8308 write_file_destroy (state.write_files[i]);
8309 free (state.write_files);
8310 fh_unref (state.prev_save_file);
8311 for (size_t i = 0; i < state.n_save_files; i++)
8312 save_file_destroy (state.save_files[i]);
8313 free (state.save_files);