1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2021 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include <gsl/gsl_blas.h>
20 #include <gsl/gsl_cdf.h>
21 #include <gsl/gsl_eigen.h>
22 #include <gsl/gsl_linalg.h>
23 #include <gsl/gsl_matrix.h>
24 #include <gsl/gsl_permutation.h>
25 #include <gsl/gsl_randist.h>
26 #include <gsl/gsl_vector.h>
31 #include "data/any-reader.h"
32 #include "data/any-writer.h"
33 #include "data/casereader.h"
34 #include "data/casewriter.h"
35 #include "data/data-in.h"
36 #include "data/data-out.h"
37 #include "data/dataset.h"
38 #include "data/dictionary.h"
39 #include "data/file-handle-def.h"
40 #include "language/command.h"
41 #include "language/data-io/data-reader.h"
42 #include "language/data-io/data-writer.h"
43 #include "language/data-io/file-handle.h"
44 #include "language/lexer/format-parser.h"
45 #include "language/lexer/lexer.h"
46 #include "language/lexer/variable-parser.h"
47 #include "libpspp/array.h"
48 #include "libpspp/assertion.h"
49 #include "libpspp/compiler.h"
50 #include "libpspp/hmap.h"
51 #include "libpspp/i18n.h"
52 #include "libpspp/intern.h"
53 #include "libpspp/misc.h"
54 #include "libpspp/str.h"
55 #include "libpspp/string-array.h"
56 #include "libpspp/stringi-set.h"
57 #include "libpspp/u8-line.h"
58 #include "math/distributions.h"
59 #include "math/random.h"
60 #include "output/driver.h"
61 #include "output/output-item.h"
62 #include "output/pivot-table.h"
64 #include "gl/c-ctype.h"
65 #include "gl/c-strcase.h"
66 #include "gl/ftoastr.h"
67 #include "gl/minmax.h"
71 #define _(msgid) gettext (msgid)
72 #define N_(msgid) (msgid)
76 struct hmap_node hmap_node;
84 struct file_handle *outfile;
85 struct string_array variables;
86 struct string_array fnames;
87 struct string_array snames;
90 /* Collects and owns factors and splits. The individual msave_command
91 structs point to these but do not own them. */
92 struct matrix_expr **factors;
93 size_t n_factors, allocated_factors;
94 struct matrix_expr **splits;
95 size_t n_splits, allocated_splits;
97 /* Execution state. */
98 struct dictionary *dict;
99 struct casewriter *writer;
104 struct file_handle *file;
105 struct dfm_reader *reader;
111 struct file_handle *file;
112 struct dfm_writer *writer;
114 struct u8_line *held;
119 struct file_handle *file;
120 struct dataset *dataset;
122 /* Parameters from parsing the first SAVE command for 'file'. */
123 struct string_array variables;
124 struct matrix_expr *names;
125 struct stringi_set strings;
127 /* Results from the first (only) attempt to open this save_file. */
129 struct casewriter *writer;
130 struct dictionary *dict;
135 struct dataset *dataset;
136 struct session *session;
140 struct msave_common *common;
142 struct file_handle *prev_read_file;
143 struct read_file **read_files;
146 struct file_handle *prev_write_file;
147 struct write_file **write_files;
148 size_t n_write_files;
150 struct file_handle *prev_save_file;
151 struct save_file **save_files;
155 static struct matrix_var *
156 matrix_var_lookup (struct matrix_state *s, struct substring name)
158 struct matrix_var *var;
160 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
161 utf8_hash_case_substring (name, 0), &s->vars)
162 if (!utf8_sscasecmp (ss_cstr (var->name), name))
168 static struct matrix_var *
169 matrix_var_create (struct matrix_state *s, struct substring name)
171 struct matrix_var *var = xmalloc (sizeof *var);
172 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
173 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
178 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
180 gsl_matrix_free (var->value);
184 /* The third argument to F() is a "prototype". For most prototypes, the first
185 letter (before the _) represents the return type and each other letter
186 (after the _) is an argument type. The types are:
188 - "m": A matrix of unrestricted dimensions.
192 - "v": A row or column vector.
194 - "e": Primarily for the first argument, this is a matrix with
195 unrestricted dimensions treated elementwise. Each element in the matrix
196 is passed to the implementation function separately.
198 - "n": This gets passed the "const struct matrix_expr *" that represents
199 the expression. This allows the evaluation function to grab the source
200 location of arguments so that it can report accurate error locations.
202 The fourth argument is an optional constraints string. For this purpose the
203 first argument is named "a", the second "b", and so on. The following kinds
204 of constraints are supported. For matrix arguments, the constraints are
205 applied to each value in the matrix separately:
207 - "a(0,1)" or "a[0,1]": 0 < a < 1 or 0 <= a <= 1, respectively. Any
208 integer may substitute for 0 and 1. Half-open constraints (] and [) are
211 - "ai": Restrict a to integer values.
213 - "a>0", "a<0", "a>=0", "a<=0", "a!=0".
215 - "a<b", "a>b", "a<=b", "a>=b", "b!=0".
217 #define MATRIX_FUNCTIONS \
218 F(ABS, "ABS", m_e, NULL) \
219 F(ALL, "ALL", d_m, NULL) \
220 F(ANY, "ANY", d_m, NULL) \
221 F(ARSIN, "ARSIN", m_e, "a[-1,1]") \
222 F(ARTAN, "ARTAN", m_e, NULL) \
223 F(BLOCK, "BLOCK", m_any, NULL) \
224 F(CHOL, "CHOL", m_mn, NULL) \
225 F(CMIN, "CMIN", m_m, NULL) \
226 F(CMAX, "CMAX", m_m, NULL) \
227 F(COS, "COS", m_e, NULL) \
228 F(CSSQ, "CSSQ", m_m, NULL) \
229 F(CSUM, "CSUM", m_m, NULL) \
230 F(DESIGN, "DESIGN", m_mn, NULL) \
231 F(DET, "DET", d_m, NULL) \
232 F(DIAG, "DIAG", m_m, NULL) \
233 F(EVAL, "EVAL", m_mn, NULL) \
234 F(EXP, "EXP", m_e, NULL) \
235 F(GINV, "GINV", m_m, NULL) \
236 F(GRADE, "GRADE", m_m, NULL) \
237 F(GSCH, "GSCH", m_mn, NULL) \
238 F(IDENT, "IDENT", IDENT, "ai>=0 bi>=0") \
239 F(INV, "INV", m_m, NULL) \
240 F(KRONEKER, "KRONEKER", m_mm, NULL) \
241 F(LG10, "LG10", m_e, "a>0") \
242 F(LN, "LN", m_e, "a>0") \
243 F(MAGIC, "MAGIC", m_d, "ai>=3") \
244 F(MAKE, "MAKE", m_ddd, "ai>=0 bi>=0") \
245 F(MDIAG, "MDIAG", m_v, NULL) \
246 F(MMAX, "MMAX", d_m, NULL) \
247 F(MMIN, "MMIN", d_m, NULL) \
248 F(MOD, "MOD", m_md, "b!=0") \
249 F(MSSQ, "MSSQ", d_m, NULL) \
250 F(MSUM, "MSUM", d_m, NULL) \
251 F(NCOL, "NCOL", d_m, NULL) \
252 F(NROW, "NROW", d_m, NULL) \
253 F(RANK, "RANK", d_m, NULL) \
254 F(RESHAPE, "RESHAPE", m_mddn, NULL) \
255 F(RMAX, "RMAX", m_m, NULL) \
256 F(RMIN, "RMIN", m_m, NULL) \
257 F(RND, "RND", m_e, NULL) \
258 F(RNKORDER, "RNKORDER", m_m, NULL) \
259 F(RSSQ, "RSSQ", m_m, NULL) \
260 F(RSUM, "RSUM", m_m, NULL) \
261 F(SIN, "SIN", m_e, NULL) \
262 F(SOLVE, "SOLVE", m_mmn, NULL) \
263 F(SQRT, "SQRT", m_e, "a>=0") \
264 F(SSCP, "SSCP", m_m, NULL) \
265 F(SVAL, "SVAL", m_m, NULL) \
266 F(SWEEP, "SWEEP", m_mdn, NULL) \
267 F(T, "T", m_m, NULL) \
268 F(TRACE, "TRACE", d_m, NULL) \
269 F(TRANSPOS, "TRANSPOS", m_m, NULL) \
270 F(TRUNC, "TRUNC", m_e, NULL) \
271 F(UNIFORM, "UNIFORM", m_ddn, "ai>=0 bi>=0") \
272 F(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0") \
273 F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0") \
274 F(IDF_BETA, "IDF.BETA", m_edd, "a[0,1] b>0 c>0") \
275 F(RV_BETA, "RV.BETA", d_dd, "a>0 b>0") \
276 F(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
277 F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
278 F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]") \
279 F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]") \
280 F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0") \
281 F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0") \
282 F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "c>0") \
283 F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \
284 F(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0") \
285 F(CHICDF, "CHICDF", m_ed, "a>=0 b>0") \
286 F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0") \
287 F(PDF_CHISQ, "PDF.CHISQ", m_ed, "a>=0 b>0") \
288 F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \
289 F(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0") \
290 F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0") \
291 F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0") \
292 F(PDF_EXP, "PDF.EXP", m_ed, "a>=0 b>0") \
293 F(RV_EXP, "RV.EXP", d_d, "a>0") \
294 F(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0") \
295 F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \
296 F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0") \
297 F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0") \
298 F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0") \
299 F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0") \
300 F(RV_F, "RV.F", d_dd, "a>0 b>0") \
301 F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0") \
302 F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
303 F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0") \
304 F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
305 F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0") \
306 F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL) \
307 F(RV_LANDAU, "RV.LANDAU", d_none, NULL) \
308 F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0") \
309 F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0") \
310 F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0") \
311 F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0") \
312 F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]") \
313 F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]") \
314 F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0") \
315 F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0") \
316 F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0") \
317 F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0") \
318 F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
319 F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0") \
320 F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
321 F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0") \
322 F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0") \
323 F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0") \
324 F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0") \
325 F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0") \
326 F(CDFNORM, "CDFNORM", m_e, NULL) \
327 F(PROBIT, "PROBIT", m_e, "a(0,1)") \
328 F(NORMAL, "NORMAL", m_e, "a>0") \
329 F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0") \
330 F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0") \
331 F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0") \
332 F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0") \
333 F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0") \
334 F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0") \
335 F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0") \
336 F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0") \
337 F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0") \
338 F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0") \
339 F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL) \
340 F(RV_RTAIL, "RV.RTAIL", d_dd, NULL) \
341 F(CDF_T, "CDF.T", m_ed, "b>0") \
342 F(TCDF, "TCDF", m_ed, "b>0") \
343 F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0") \
344 F(PDF_T, "PDF.T", m_ed, "b>0") \
345 F(RV_T, "RV.T", d_d, "a>0") \
346 F(CDF_T1G, "CDF.T1G", m_edd, NULL) \
347 F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)") \
348 F(PDF_T1G, "PDF.T1G", m_edd, NULL) \
349 F(RV_T1G, "RV.T1G", d_dd, NULL) \
350 F(CDF_T2G, "CDF.T2G", m_edd, NULL) \
351 F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)") \
352 F(PDF_T2G, "PDF.T2G", m_edd, NULL) \
353 F(RV_T2G, "RV.T2G", d_dd, NULL) \
354 F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c") \
355 F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c") \
356 F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c") \
357 F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b") \
358 F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
359 F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0") \
360 F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
361 F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0") \
362 F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
363 F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
364 F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]") \
365 F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]") \
366 F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]") \
367 F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]") \
368 F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]") \
369 F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]") \
370 F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]") \
371 F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
372 F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
373 F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a") \
374 F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]") \
375 F(RV_LOG, "RV.LOG", d_d, "a(0,1]") \
376 F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
377 F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
378 F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]") \
379 F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0") \
380 F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0") \
381 F(RV_POISSON, "RV.POISSON", d_d, "a>0")
383 struct matrix_function_properties
386 const char *constraints;
389 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
390 enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
391 enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
392 enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 };
393 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
394 enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 };
395 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
396 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
397 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
398 enum { m_ddn_MIN_ARGS = 2, m_ddn_MAX_ARGS = 2 };
399 enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 };
400 enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 };
401 enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
402 enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
403 enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
404 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
405 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
406 enum { m_mddn_MIN_ARGS = 3, m_mddn_MAX_ARGS = 3 };
407 enum { m_mdn_MIN_ARGS = 2, m_mdn_MAX_ARGS = 2 };
408 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
409 enum { m_mmn_MIN_ARGS = 2, m_mmn_MAX_ARGS = 2 };
410 enum { m_mn_MIN_ARGS = 1, m_mn_MAX_ARGS = 1 };
411 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
413 typedef double matrix_proto_d_none (void);
414 typedef double matrix_proto_d_d (double);
415 typedef double matrix_proto_d_dd (double, double);
416 typedef double matrix_proto_d_dd (double, double);
417 typedef double matrix_proto_d_ddd (double, double, double);
418 typedef gsl_matrix *matrix_proto_m_d (double);
419 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
420 typedef gsl_matrix *matrix_proto_m_ddn (double, double,
421 const struct matrix_expr *);
422 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
423 typedef gsl_matrix *matrix_proto_m_mn (gsl_matrix *,
424 const struct matrix_expr *);
425 typedef double matrix_proto_m_e (double);
426 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
427 typedef gsl_matrix *matrix_proto_m_mdn (gsl_matrix *, double,
428 const struct matrix_expr *);
429 typedef double matrix_proto_m_ed (double, double);
430 typedef gsl_matrix *matrix_proto_m_mddn (gsl_matrix *, double, double,
431 const struct matrix_expr *);
432 typedef double matrix_proto_m_edd (double, double, double);
433 typedef double matrix_proto_m_eddd (double, double, double, double);
434 typedef double matrix_proto_m_eed (double, double, double);
435 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
436 typedef gsl_matrix *matrix_proto_m_mmn (gsl_matrix *, gsl_matrix *,
437 const struct matrix_expr *);
438 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
439 typedef double matrix_proto_d_m (gsl_matrix *);
440 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
441 typedef gsl_matrix *matrix_proto_IDENT (double, double);
443 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
444 static matrix_proto_##PROTO matrix_eval_##ENUM;
448 /* Matrix expressions. */
455 #define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
459 /* Elementwise and scalar arithmetic. */
460 MOP_NEGATE, /* unary - */
461 MOP_ADD_ELEMS, /* + */
462 MOP_SUB_ELEMS, /* - */
463 MOP_MUL_ELEMS, /* &* */
464 MOP_DIV_ELEMS, /* / and &/ */
465 MOP_EXP_ELEMS, /* &** */
467 MOP_SEQ_BY, /* a:b:c */
469 /* Matrix arithmetic. */
471 MOP_EXP_MAT, /* ** */
488 MOP_PASTE_HORZ, /* a, b, c, ... */
489 MOP_PASTE_VERT, /* a; b; c; ... */
493 MOP_VEC_INDEX, /* x(y) */
494 MOP_VEC_ALL, /* x(:) */
495 MOP_MAT_INDEX, /* x(y,z) */
496 MOP_ROW_INDEX, /* x(y,:) */
497 MOP_COL_INDEX, /* x(:,z) */
504 MOP_EOF, /* EOF('file') */
512 struct matrix_expr **subs;
517 struct matrix_var *variable;
518 struct read_file *eof;
521 struct msg_location *location;
525 matrix_expr_location__ (const struct matrix_expr *e,
526 const struct msg_location **minp,
527 const struct msg_location **maxp)
529 struct msg_location *loc = e->location;
534 || loc->p[0].line < (*minp)->p[0].line
535 || (loc->p[0].line == (*minp)->p[0].line
536 && loc->p[0].column < (*minp)->p[0].column)))
541 || loc->p[1].line > (*maxp)->p[1].line
542 || (loc->p[1].line == (*maxp)->p[1].line
543 && loc->p[1].column > (*maxp)->p[1].column)))
548 assert (e->op != MOP_NUMBER && e->op != MOP_VARIABLE && e->op != MOP_EOF);
549 for (size_t i = 0; i < e->n_subs; i++)
550 matrix_expr_location__ (e->subs[i], minp, maxp);
553 static struct msg_location *
554 matrix_expr_location (const struct matrix_expr *e_)
556 struct matrix_expr *e = CONST_CAST (struct matrix_expr *, e_);
560 const struct msg_location *min = NULL;
561 const struct msg_location *max = NULL;
562 matrix_expr_location__ (e, &min, &max);
565 e->location = msg_location_dup (min);
566 e->location->p[1] = max->p[1];
573 static struct matrix_expr *
574 matrix_expr_wrap_location (struct matrix_state *s, int start_ofs,
575 struct matrix_expr *e)
577 if (e && !e->location)
578 e->location = lex_ofs_location (s->lexer, start_ofs,
579 lex_ofs (s->lexer) - 1);
584 matrix_expr_destroy (struct matrix_expr *e)
591 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
622 for (size_t i = 0; i < e->n_subs; i++)
623 matrix_expr_destroy (e->subs[i]);
632 msg_location_destroy (e->location);
636 static struct matrix_expr *
637 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
640 struct matrix_expr *e = xmalloc (sizeof *e);
641 *e = (struct matrix_expr) {
643 .subs = xmemdup (subs, n_subs * sizeof *subs),
650 static struct matrix_expr *
651 matrix_expr_create_0 (enum matrix_op op)
653 struct matrix_expr *sub;
654 return matrix_expr_create_subs (op, &sub, 0);
657 static struct matrix_expr *
658 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
660 return matrix_expr_create_subs (op, &sub, 1);
663 static struct matrix_expr *
664 matrix_expr_create_2 (enum matrix_op op,
665 struct matrix_expr *sub0, struct matrix_expr *sub1)
667 struct matrix_expr *subs[] = { sub0, sub1 };
668 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
671 static struct matrix_expr *
672 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
673 struct matrix_expr *sub1, struct matrix_expr *sub2)
675 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
676 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
679 static struct matrix_expr *
680 matrix_expr_create_number (struct lexer *lexer, double number)
682 struct matrix_expr *e = xmalloc (sizeof *e);
683 *e = (struct matrix_expr) {
691 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
693 static struct matrix_expr *
694 matrix_parse_curly_comma (struct matrix_state *s)
696 struct matrix_expr *lhs = matrix_parse_expr (s);
700 while (lex_match (s->lexer, T_COMMA))
702 struct matrix_expr *rhs = matrix_parse_expr (s);
705 matrix_expr_destroy (lhs);
708 lhs = matrix_expr_create_2 (MOP_PASTE_HORZ, lhs, rhs);
713 static struct matrix_expr *
714 matrix_parse_curly_semi (struct matrix_state *s)
716 if (lex_token (s->lexer) == T_RCURLY)
717 return matrix_expr_create_0 (MOP_EMPTY);
719 struct matrix_expr *lhs = matrix_parse_curly_comma (s);
723 while (lex_match (s->lexer, T_SEMICOLON))
725 struct matrix_expr *rhs = matrix_parse_curly_comma (s);
728 matrix_expr_destroy (lhs);
731 lhs = matrix_expr_create_2 (MOP_PASTE_VERT, lhs, rhs);
736 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
737 for (size_t Y = 0; Y < (M)->size1; Y++) \
738 for (size_t X = 0; X < (M)->size2; X++) \
739 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
742 is_vector (const gsl_matrix *m)
744 return m->size1 <= 1 || m->size2 <= 1;
748 to_vector (gsl_matrix *m)
750 return (m->size1 == 1
751 ? gsl_matrix_row (m, 0).vector
752 : gsl_matrix_column (m, 0).vector);
757 matrix_eval_ABS (double d)
763 matrix_eval_ALL (gsl_matrix *m)
765 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
772 matrix_eval_ANY (gsl_matrix *m)
774 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
781 matrix_eval_ARSIN (double d)
787 matrix_eval_ARTAN (double d)
793 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
797 for (size_t i = 0; i < n; i++)
802 gsl_matrix *block = gsl_matrix_calloc (r, c);
804 for (size_t i = 0; i < n; i++)
806 for (size_t y = 0; y < m[i]->size1; y++)
807 for (size_t x = 0; x < m[i]->size2; x++)
808 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
816 matrix_eval_CHOL (gsl_matrix *m, const struct matrix_expr *e)
818 if (!gsl_linalg_cholesky_decomp1 (m))
820 for (size_t y = 0; y < m->size1; y++)
821 for (size_t x = y + 1; x < m->size2; x++)
822 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
824 for (size_t y = 0; y < m->size1; y++)
825 for (size_t x = 0; x < y; x++)
826 gsl_matrix_set (m, y, x, 0);
831 msg_at (SE, e->subs[0]->location,
832 _("Input to CHOL function is not positive-definite."));
838 matrix_eval_col_extremum (gsl_matrix *m, bool min)
843 return gsl_matrix_alloc (1, 0);
845 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
846 for (size_t x = 0; x < m->size2; x++)
848 double ext = gsl_matrix_get (m, 0, x);
849 for (size_t y = 1; y < m->size1; y++)
851 double value = gsl_matrix_get (m, y, x);
852 if (min ? value < ext : value > ext)
855 gsl_matrix_set (cext, 0, x, ext);
861 matrix_eval_CMAX (gsl_matrix *m)
863 return matrix_eval_col_extremum (m, false);
867 matrix_eval_CMIN (gsl_matrix *m)
869 return matrix_eval_col_extremum (m, true);
873 matrix_eval_COS (double d)
879 matrix_eval_col_sum (gsl_matrix *m, bool square)
884 return gsl_matrix_alloc (1, 0);
886 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
887 for (size_t x = 0; x < m->size2; x++)
890 for (size_t y = 0; y < m->size1; y++)
892 double d = gsl_matrix_get (m, y, x);
893 sum += square ? pow2 (d) : d;
895 gsl_matrix_set (result, 0, x, sum);
901 matrix_eval_CSSQ (gsl_matrix *m)
903 return matrix_eval_col_sum (m, true);
907 matrix_eval_CSUM (gsl_matrix *m)
909 return matrix_eval_col_sum (m, false);
913 compare_double_3way (const void *a_, const void *b_)
915 const double *a = a_;
916 const double *b = b_;
917 return *a < *b ? -1 : *a > *b;
921 matrix_eval_DESIGN (gsl_matrix *m, const struct matrix_expr *e)
923 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
924 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
925 gsl_matrix_transpose_memcpy (&m2, m);
927 for (size_t y = 0; y < m2.size1; y++)
928 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
930 size_t *n = xcalloc (m2.size1, sizeof *n);
932 for (size_t i = 0; i < m2.size1; i++)
934 double *row = tmp + m2.size2 * i;
935 for (size_t j = 0; j < m2.size2; )
938 for (k = j + 1; k < m2.size2; k++)
939 if (row[j] != row[k])
941 row[n[i]++] = row[j];
946 msg_at (MW, e->subs[0]->location,
947 _("Column %zu in DESIGN argument has constant value."), i + 1);
952 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
954 for (size_t i = 0; i < m->size2; i++)
959 const double *unique = tmp + m2.size2 * i;
960 for (size_t j = 0; j < n[i]; j++, x++)
962 double value = unique[j];
963 for (size_t y = 0; y < m->size1; y++)
964 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
975 matrix_eval_DET (gsl_matrix *m)
977 gsl_permutation *p = gsl_permutation_alloc (m->size1);
979 gsl_linalg_LU_decomp (m, p, &signum);
980 gsl_permutation_free (p);
981 return gsl_linalg_LU_det (m, signum);
985 matrix_eval_DIAG (gsl_matrix *m)
987 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
988 for (size_t i = 0; i < diag->size1; i++)
989 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
994 is_symmetric (const gsl_matrix *m)
996 if (m->size1 != m->size2)
999 for (size_t y = 0; y < m->size1; y++)
1000 for (size_t x = 0; x < y; x++)
1001 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
1008 compare_double_desc (const void *a_, const void *b_)
1010 const double *a = a_;
1011 const double *b = b_;
1012 return *a > *b ? -1 : *a < *b;
1016 matrix_eval_EVAL (gsl_matrix *m, const struct matrix_expr *e)
1018 if (!is_symmetric (m))
1020 msg_at (SE, e->subs[0]->location,
1021 _("Argument of EVAL must be symmetric."));
1025 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
1026 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
1027 gsl_vector v_eval = to_vector (eval);
1028 gsl_eigen_symm (m, &v_eval, w);
1029 gsl_eigen_symm_free (w);
1031 assert (v_eval.stride == 1);
1032 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
1038 matrix_eval_EXP (double d)
1043 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
1046 Charl Linssen <charl@itfromb.it>
1050 matrix_eval_GINV (gsl_matrix *A)
1052 size_t n = A->size1;
1053 size_t m = A->size2;
1055 gsl_matrix *tmp_mat = NULL;
1058 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
1059 tmp_mat = gsl_matrix_alloc (m, n);
1060 gsl_matrix_transpose_memcpy (tmp_mat, A);
1068 gsl_matrix *V = gsl_matrix_alloc (m, m);
1069 gsl_vector *u = gsl_vector_alloc (m);
1071 gsl_vector *tmp_vec = gsl_vector_alloc (m);
1072 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
1073 gsl_vector_free (tmp_vec);
1076 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
1077 gsl_matrix_set_zero (Sigma_pinv);
1078 double cutoff = 1e-15 * gsl_vector_max (u);
1080 for (size_t i = 0; i < m; ++i)
1082 double x = gsl_vector_get (u, i);
1083 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
1086 /* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
1087 gsl_matrix *U = gsl_matrix_calloc (n, n);
1088 for (size_t i = 0; i < n; i++)
1089 for (size_t j = 0; j < m; j++)
1090 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
1092 /* two dot products to obtain pseudoinverse */
1093 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
1094 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
1099 A_pinv = gsl_matrix_alloc (n, m);
1100 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
1104 A_pinv = gsl_matrix_alloc (m, n);
1105 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
1108 gsl_matrix_free (tmp_mat);
1109 gsl_matrix_free (tmp_mat2);
1110 gsl_matrix_free (U);
1111 gsl_matrix_free (Sigma_pinv);
1112 gsl_vector_free (u);
1113 gsl_matrix_free (V);
1125 grade_compare_3way (const void *a_, const void *b_)
1127 const struct grade *a = a_;
1128 const struct grade *b = b_;
1130 return (a->value < b->value ? -1
1131 : a->value > b->value ? 1
1139 matrix_eval_GRADE (gsl_matrix *m)
1141 size_t n = m->size1 * m->size2;
1142 struct grade *grades = xmalloc (n * sizeof *grades);
1145 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1146 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
1147 qsort (grades, n, sizeof *grades, grade_compare_3way);
1149 for (size_t i = 0; i < n; i++)
1150 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
1158 dot (gsl_vector *a, gsl_vector *b)
1160 double result = 0.0;
1161 for (size_t i = 0; i < a->size; i++)
1162 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
1167 norm2 (gsl_vector *v)
1169 double result = 0.0;
1170 for (size_t i = 0; i < v->size; i++)
1171 result += pow2 (gsl_vector_get (v, i));
1176 norm (gsl_vector *v)
1178 return sqrt (norm2 (v));
1182 matrix_eval_GSCH (gsl_matrix *v, const struct matrix_expr *e)
1184 if (v->size2 < v->size1)
1186 msg_at (SE, e->subs[0]->location,
1187 _("GSCH requires its argument to have at least as many columns "
1188 "as rows, but it has dimensions %zu×%zu."),
1189 v->size1, v->size2);
1192 if (!v->size1 || !v->size2)
1195 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
1197 for (size_t vx = 0; vx < v->size2; vx++)
1199 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
1200 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
1202 gsl_vector_memcpy (&u_i, &v_i);
1203 for (size_t j = 0; j < ux; j++)
1205 gsl_vector u_j = gsl_matrix_column (u, j).vector;
1206 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
1207 for (size_t k = 0; k < u_i.size; k++)
1208 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
1209 - scale * gsl_vector_get (&u_j, k)));
1212 double len = norm (&u_i);
1215 gsl_vector_scale (&u_i, 1.0 / len);
1216 if (++ux >= v->size1)
1223 msg_at (SE, e->subs[0]->location,
1224 _("%zu×%zu argument to GSCH contains only "
1225 "%zu linearly independent columns."),
1226 v->size1, v->size2, ux);
1227 gsl_matrix_free (u);
1231 u->size2 = v->size1;
1236 matrix_eval_IDENT (double s1, double s2)
1238 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
1239 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1245 invert_matrix (gsl_matrix *x)
1247 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1249 gsl_linalg_LU_decomp (x, p, &signum);
1250 gsl_linalg_LU_invx (x, p);
1251 gsl_permutation_free (p);
1255 matrix_eval_INV (gsl_matrix *m)
1262 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1264 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1265 a->size2 * b->size2);
1267 for (size_t ar = 0; ar < a->size1; ar++)
1268 for (size_t br = 0; br < b->size1; br++, y++)
1271 for (size_t ac = 0; ac < a->size2; ac++)
1272 for (size_t bc = 0; bc < b->size2; bc++, x++)
1274 double av = gsl_matrix_get (a, ar, ac);
1275 double bv = gsl_matrix_get (b, br, bc);
1276 gsl_matrix_set (k, y, x, av * bv);
1283 matrix_eval_LG10 (double d)
1289 matrix_eval_LN (double d)
1295 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1297 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1300 for (size_t i = 1; i <= n * n; i++)
1302 gsl_matrix_set (m, y, x, i);
1304 size_t y1 = !y ? n - 1 : y - 1;
1305 size_t x1 = x + 1 >= n ? 0 : x + 1;
1306 if (gsl_matrix_get (m, y1, x1) == 0)
1312 y = y + 1 >= n ? 0 : y + 1;
1317 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1319 double a = gsl_matrix_get (m, y1, x1);
1320 double b = gsl_matrix_get (m, y2, x2);
1321 gsl_matrix_set (m, y1, x1, b);
1322 gsl_matrix_set (m, y2, x2, a);
1326 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1330 /* A. Umar, "On the Construction of Even Order Magic Squares",
1331 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1333 for (size_t i = 1; i <= n * n / 2; i++)
1335 gsl_matrix_set (m, y, x, i);
1345 for (size_t i = n * n; i > n * n / 2; i--)
1347 gsl_matrix_set (m, y, x, i);
1355 for (size_t y = 0; y < n; y++)
1356 for (size_t x = 0; x < n / 2; x++)
1358 unsigned int d = gsl_matrix_get (m, y, x);
1359 if (d % 2 != (y < n / 2))
1360 magic_exchange (m, y, x, y, n - x - 1);
1365 size_t x1 = n / 2 - 1;
1367 magic_exchange (m, y1, x1, y2, x1);
1368 magic_exchange (m, y1, x2, y2, x2);
1372 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1374 /* A. Umar, "On the Construction of Even Order Magic Squares",
1375 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1379 for (size_t i = 1; ; i++)
1381 gsl_matrix_set (m, y, x, i);
1382 if (++y == n / 2 - 1)
1394 for (size_t i = n * n; ; i--)
1396 gsl_matrix_set (m, y, x, i);
1397 if (++y == n / 2 - 1)
1406 for (size_t y = 0; y < n; y++)
1407 if (y != n / 2 - 1 && y != n / 2)
1408 for (size_t x = 0; x < n / 2; x++)
1410 unsigned int d = gsl_matrix_get (m, y, x);
1411 if (d % 2 != (y < n / 2))
1412 magic_exchange (m, y, x, y, n - x - 1);
1415 size_t a0 = (n * n - 2 * n) / 2 + 1;
1416 for (size_t i = 0; i < n / 2; i++)
1419 gsl_matrix_set (m, n / 2, i, a);
1420 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
1422 for (size_t i = 0; i < n / 2; i++)
1424 size_t a = a0 + i + n / 2;
1425 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
1426 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
1428 for (size_t x = 1; x < n / 2; x += 2)
1429 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1430 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
1431 magic_exchange (m, n / 2, x, n / 2 - 1, x);
1432 size_t x1 = n / 2 - 2;
1433 size_t x2 = n / 2 + 1;
1434 size_t y1 = n / 2 - 2;
1435 size_t y2 = n / 2 + 1;
1436 magic_exchange (m, y1, x1, y2, x1);
1437 magic_exchange (m, y1, x2, y2, x2);
1441 matrix_eval_MAGIC (double n_)
1445 gsl_matrix *m = gsl_matrix_calloc (n, n);
1447 matrix_eval_MAGIC_odd (m, n);
1449 matrix_eval_MAGIC_singly_even (m, n);
1451 matrix_eval_MAGIC_doubly_even (m, n);
1456 matrix_eval_MAKE (double r, double c, double value)
1458 gsl_matrix *m = gsl_matrix_alloc (r, c);
1459 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1465 matrix_eval_MDIAG (gsl_vector *v)
1467 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
1468 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
1469 gsl_vector_memcpy (&diagonal, v);
1474 matrix_eval_MMAX (gsl_matrix *m)
1476 return gsl_matrix_max (m);
1480 matrix_eval_MMIN (gsl_matrix *m)
1482 return gsl_matrix_min (m);
1486 matrix_eval_MOD (gsl_matrix *m, double divisor)
1488 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1489 *d = fmod (*d, divisor);
1494 matrix_eval_MSSQ (gsl_matrix *m)
1497 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1503 matrix_eval_MSUM (gsl_matrix *m)
1506 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1512 matrix_eval_NCOL (gsl_matrix *m)
1518 matrix_eval_NROW (gsl_matrix *m)
1524 matrix_eval_RANK (gsl_matrix *m)
1526 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
1527 gsl_linalg_QR_decomp (m, tau);
1528 gsl_vector_free (tau);
1530 return gsl_linalg_QRPT_rank (m, -1);
1534 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_,
1535 const struct matrix_expr *e)
1537 bool r_ok = r_ >= 0 && r_ < SIZE_MAX;
1538 bool c_ok = c_ >= 0 && c_ < SIZE_MAX;
1542 !r_ok ? e->subs[1]->location : e->subs[2]->location,
1543 _("Arguments 2 and 3 to RESHAPE must be integers."));
1548 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
1550 struct msg_location *loc = msg_location_dup (e->subs[1]->location);
1551 loc->p[1] = e->subs[2]->location->p[1];
1552 msg_at (SE, loc, _("Product of RESHAPE size arguments (%zu×%zu = %zu) "
1553 "differs from product of matrix dimensions "
1554 "(%zu×%zu = %zu)."),
1556 m->size1, m->size2, m->size1 * m->size2);
1557 msg_location_destroy (loc);
1561 gsl_matrix *dst = gsl_matrix_alloc (r, c);
1564 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
1566 gsl_matrix_set (dst, y1, x1, *d);
1577 matrix_eval_row_extremum (gsl_matrix *m, bool min)
1582 return gsl_matrix_alloc (0, 1);
1584 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
1585 for (size_t y = 0; y < m->size1; y++)
1587 double ext = gsl_matrix_get (m, y, 0);
1588 for (size_t x = 1; x < m->size2; x++)
1590 double value = gsl_matrix_get (m, y, x);
1591 if (min ? value < ext : value > ext)
1594 gsl_matrix_set (rext, y, 0, ext);
1600 matrix_eval_RMAX (gsl_matrix *m)
1602 return matrix_eval_row_extremum (m, false);
1606 matrix_eval_RMIN (gsl_matrix *m)
1608 return matrix_eval_row_extremum (m, true);
1612 matrix_eval_RND (double d)
1624 rank_compare_3way (const void *a_, const void *b_)
1626 const struct rank *a = a_;
1627 const struct rank *b = b_;
1629 return a->value < b->value ? -1 : a->value > b->value;
1633 matrix_eval_RNKORDER (gsl_matrix *m)
1635 size_t n = m->size1 * m->size2;
1636 struct rank *ranks = xmalloc (n * sizeof *ranks);
1638 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1639 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
1640 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
1642 for (size_t i = 0; i < n; )
1645 for (j = i + 1; j < n; j++)
1646 if (ranks[i].value != ranks[j].value)
1649 double rank = (i + j + 1.0) / 2.0;
1650 for (size_t k = i; k < j; k++)
1651 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
1662 matrix_eval_row_sum (gsl_matrix *m, bool square)
1667 return gsl_matrix_alloc (0, 1);
1669 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
1670 for (size_t y = 0; y < m->size1; y++)
1673 for (size_t x = 0; x < m->size2; x++)
1675 double d = gsl_matrix_get (m, y, x);
1676 sum += square ? pow2 (d) : d;
1678 gsl_matrix_set (result, y, 0, sum);
1684 matrix_eval_RSSQ (gsl_matrix *m)
1686 return matrix_eval_row_sum (m, true);
1690 matrix_eval_RSUM (gsl_matrix *m)
1692 return matrix_eval_row_sum (m, false);
1696 matrix_eval_SIN (double d)
1702 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2, const struct matrix_expr *e)
1704 if (m1->size1 != m2->size1)
1706 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
1707 loc->p[1] = e->subs[1]->location->p[1];
1711 _("SOLVE requires its arguments to have the same number of "
1712 "rows, but the first argument has dimensions %zu×%zu and "
1713 "the second %zu×%zu."),
1714 m1->size1, m1->size2,
1715 m2->size1, m2->size2);
1717 msg_at (SE, e->location, _("SOLVE requires its arguments to have the same "
1718 "number of rows."));
1719 msg_at (SN, e->subs[0]->location, _("Argument 1 has dimensions %zu×%zu."),
1720 m1->size1, m1->size2);
1721 msg_at (SN, e->subs[1]->location, _("Argument 2 has dimensions %zu×%zu."),
1722 m2->size1, m2->size2);
1725 msg_location_destroy (loc);
1729 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
1730 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
1732 gsl_linalg_LU_decomp (m1, p, &signum);
1733 for (size_t i = 0; i < m2->size2; i++)
1735 gsl_vector bi = gsl_matrix_column (m2, i).vector;
1736 gsl_vector xi = gsl_matrix_column (x, i).vector;
1737 gsl_linalg_LU_solve (m1, p, &bi, &xi);
1739 gsl_permutation_free (p);
1744 matrix_eval_SQRT (double d)
1750 matrix_eval_SSCP (gsl_matrix *m)
1752 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
1753 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
1758 matrix_eval_SVAL (gsl_matrix *m)
1760 gsl_matrix *tmp_mat = NULL;
1761 if (m->size2 > m->size1)
1763 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
1764 gsl_matrix_transpose_memcpy (tmp_mat, m);
1769 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
1770 gsl_vector *S = gsl_vector_alloc (m->size2);
1771 gsl_vector *work = gsl_vector_alloc (m->size2);
1772 gsl_linalg_SV_decomp (m, V, S, work);
1774 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
1775 for (size_t i = 0; i < m->size2; i++)
1776 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
1778 gsl_matrix_free (V);
1779 gsl_vector_free (S);
1780 gsl_vector_free (work);
1781 gsl_matrix_free (tmp_mat);
1787 matrix_eval_SWEEP (gsl_matrix *m, double d, const struct matrix_expr *e)
1789 if (d < 1 || d > SIZE_MAX)
1791 msg_at (SE, e->subs[1]->location,
1792 _("Scalar argument to SWEEP must be integer."));
1796 if (k >= MIN (m->size1, m->size2))
1798 msg_at (SE, e->subs[1]->location,
1799 _("Scalar argument to SWEEP must be integer less than or "
1800 "equal to the smaller of the matrix argument's rows and "
1805 double m_kk = gsl_matrix_get (m, k, k);
1806 if (fabs (m_kk) > 1e-19)
1808 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
1809 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
1811 double m_ij = gsl_matrix_get (m, i, j);
1812 double m_ik = gsl_matrix_get (m, i, k);
1813 double m_kj = gsl_matrix_get (m, k, j);
1814 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
1823 for (size_t i = 0; i < m->size1; i++)
1825 gsl_matrix_set (m, i, k, 0);
1826 gsl_matrix_set (m, k, i, 0);
1833 matrix_eval_TRACE (gsl_matrix *m)
1836 size_t n = MIN (m->size1, m->size2);
1837 for (size_t i = 0; i < n; i++)
1838 sum += gsl_matrix_get (m, i, i);
1843 matrix_eval_T (gsl_matrix *m)
1845 return matrix_eval_TRANSPOS (m);
1849 matrix_eval_TRANSPOS (gsl_matrix *m)
1851 if (m->size1 == m->size2)
1853 gsl_matrix_transpose (m);
1858 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
1859 gsl_matrix_transpose_memcpy (t, m);
1865 matrix_eval_TRUNC (double d)
1871 matrix_eval_UNIFORM (double r_, double c_, const struct matrix_expr *e)
1875 if (size_overflow_p (xtimes (r, xmax (c, 1))))
1877 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
1878 loc->p[1] = e->subs[1]->location->p[1];
1881 _("Product of arguments to UNIFORM exceeds memory size."));
1883 msg_location_destroy (loc);
1887 gsl_matrix *m = gsl_matrix_alloc (r, c);
1888 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1889 *d = gsl_ran_flat (get_rng (), 0, 1);
1894 matrix_eval_PDF_BETA (double x, double a, double b)
1896 return gsl_ran_beta_pdf (x, a, b);
1900 matrix_eval_CDF_BETA (double x, double a, double b)
1902 return gsl_cdf_beta_P (x, a, b);
1906 matrix_eval_IDF_BETA (double P, double a, double b)
1908 return gsl_cdf_beta_Pinv (P, a, b);
1912 matrix_eval_RV_BETA (double a, double b)
1914 return gsl_ran_beta (get_rng (), a, b);
1918 matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
1920 return ncdf_beta (x, a, b, lambda);
1924 matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
1926 return npdf_beta (x, a, b, lambda);
1930 matrix_eval_CDF_BVNOR (double x0, double x1, double r)
1932 return cdf_bvnor (x0, x1, r);
1936 matrix_eval_PDF_BVNOR (double x0, double x1, double r)
1938 return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
1942 matrix_eval_CDF_CAUCHY (double x, double a, double b)
1944 return gsl_cdf_cauchy_P ((x - a) / b, 1);
1948 matrix_eval_IDF_CAUCHY (double P, double a, double b)
1950 return a + b * gsl_cdf_cauchy_Pinv (P, 1);
1954 matrix_eval_PDF_CAUCHY (double x, double a, double b)
1956 return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
1960 matrix_eval_RV_CAUCHY (double a, double b)
1962 return a + b * gsl_ran_cauchy (get_rng (), 1);
1966 matrix_eval_CDF_CHISQ (double x, double df)
1968 return gsl_cdf_chisq_P (x, df);
1972 matrix_eval_CHICDF (double x, double df)
1974 return matrix_eval_CDF_CHISQ (x, df);
1978 matrix_eval_IDF_CHISQ (double P, double df)
1980 return gsl_cdf_chisq_Pinv (P, df);
1984 matrix_eval_PDF_CHISQ (double x, double df)
1986 return gsl_ran_chisq_pdf (x, df);
1990 matrix_eval_RV_CHISQ (double df)
1992 return gsl_ran_chisq (get_rng (), df);
1996 matrix_eval_SIG_CHISQ (double x, double df)
1998 return gsl_cdf_chisq_Q (x, df);
2002 matrix_eval_CDF_EXP (double x, double a)
2004 return gsl_cdf_exponential_P (x, 1. / a);
2008 matrix_eval_IDF_EXP (double P, double a)
2010 return gsl_cdf_exponential_Pinv (P, 1. / a);
2014 matrix_eval_PDF_EXP (double x, double a)
2016 return gsl_ran_exponential_pdf (x, 1. / a);
2020 matrix_eval_RV_EXP (double a)
2022 return gsl_ran_exponential (get_rng (), 1. / a);
2026 matrix_eval_PDF_XPOWER (double x, double a, double b)
2028 return gsl_ran_exppow_pdf (x, a, b);
2032 matrix_eval_RV_XPOWER (double a, double b)
2034 return gsl_ran_exppow (get_rng (), a, b);
2038 matrix_eval_CDF_F (double x, double df1, double df2)
2040 return gsl_cdf_fdist_P (x, df1, df2);
2044 matrix_eval_FCDF (double x, double df1, double df2)
2046 return matrix_eval_CDF_F (x, df1, df2);
2050 matrix_eval_IDF_F (double P, double df1, double df2)
2052 return idf_fdist (P, df1, df2);
2056 matrix_eval_RV_F (double df1, double df2)
2058 return gsl_ran_fdist (get_rng (), df1, df2);
2062 matrix_eval_PDF_F (double x, double df1, double df2)
2064 return gsl_ran_fdist_pdf (x, df1, df2);
2068 matrix_eval_SIG_F (double x, double df1, double df2)
2070 return gsl_cdf_fdist_Q (x, df1, df2);
2074 matrix_eval_CDF_GAMMA (double x, double a, double b)
2076 return gsl_cdf_gamma_P (x, a, 1. / b);
2080 matrix_eval_IDF_GAMMA (double P, double a, double b)
2082 return gsl_cdf_gamma_Pinv (P, a, 1. / b);
2086 matrix_eval_PDF_GAMMA (double x, double a, double b)
2088 return gsl_ran_gamma_pdf (x, a, 1. / b);
2092 matrix_eval_RV_GAMMA (double a, double b)
2094 return gsl_ran_gamma (get_rng (), a, 1. / b);
2098 matrix_eval_PDF_LANDAU (double x)
2100 return gsl_ran_landau_pdf (x);
2104 matrix_eval_RV_LANDAU (void)
2106 return gsl_ran_landau (get_rng ());
2110 matrix_eval_CDF_LAPLACE (double x, double a, double b)
2112 return gsl_cdf_laplace_P ((x - a) / b, 1);
2116 matrix_eval_IDF_LAPLACE (double P, double a, double b)
2118 return a + b * gsl_cdf_laplace_Pinv (P, 1);
2122 matrix_eval_PDF_LAPLACE (double x, double a, double b)
2124 return gsl_ran_laplace_pdf ((x - a) / b, 1);
2128 matrix_eval_RV_LAPLACE (double a, double b)
2130 return a + b * gsl_ran_laplace (get_rng (), 1);
2134 matrix_eval_RV_LEVY (double c, double alpha)
2136 return gsl_ran_levy (get_rng (), c, alpha);
2140 matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
2142 return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
2146 matrix_eval_CDF_LOGISTIC (double x, double a, double b)
2148 return gsl_cdf_logistic_P ((x - a) / b, 1);
2152 matrix_eval_IDF_LOGISTIC (double P, double a, double b)
2154 return a + b * gsl_cdf_logistic_Pinv (P, 1);
2158 matrix_eval_PDF_LOGISTIC (double x, double a, double b)
2160 return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
2164 matrix_eval_RV_LOGISTIC (double a, double b)
2166 return a + b * gsl_ran_logistic (get_rng (), 1);
2170 matrix_eval_CDF_LNORMAL (double x, double m, double s)
2172 return gsl_cdf_lognormal_P (x, log (m), s);
2176 matrix_eval_IDF_LNORMAL (double P, double m, double s)
2178 return gsl_cdf_lognormal_Pinv (P, log (m), s);;
2182 matrix_eval_PDF_LNORMAL (double x, double m, double s)
2184 return gsl_ran_lognormal_pdf (x, log (m), s);
2188 matrix_eval_RV_LNORMAL (double m, double s)
2190 return gsl_ran_lognormal (get_rng (), log (m), s);
2194 matrix_eval_CDF_NORMAL (double x, double u, double s)
2196 return gsl_cdf_gaussian_P (x - u, s);
2200 matrix_eval_IDF_NORMAL (double P, double u, double s)
2202 return u + gsl_cdf_gaussian_Pinv (P, s);
2206 matrix_eval_PDF_NORMAL (double x, double u, double s)
2208 return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
2212 matrix_eval_RV_NORMAL (double u, double s)
2214 return u + gsl_ran_gaussian (get_rng (), s);
2218 matrix_eval_CDFNORM (double x)
2220 return gsl_cdf_ugaussian_P (x);
2224 matrix_eval_PROBIT (double P)
2226 return gsl_cdf_ugaussian_Pinv (P);
2230 matrix_eval_NORMAL (double s)
2232 return gsl_ran_gaussian (get_rng (), s);
2236 matrix_eval_PDF_NTAIL (double x, double a, double sigma)
2238 return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
2242 matrix_eval_RV_NTAIL (double a, double sigma)
2244 return gsl_ran_gaussian_tail (get_rng (), a, sigma);
2248 matrix_eval_CDF_PARETO (double x, double a, double b)
2250 return gsl_cdf_pareto_P (x, b, a);
2254 matrix_eval_IDF_PARETO (double P, double a, double b)
2256 return gsl_cdf_pareto_Pinv (P, b, a);
2260 matrix_eval_PDF_PARETO (double x, double a, double b)
2262 return gsl_ran_pareto_pdf (x, b, a);
2266 matrix_eval_RV_PARETO (double a, double b)
2268 return gsl_ran_pareto (get_rng (), b, a);
2272 matrix_eval_CDF_RAYLEIGH (double x, double sigma)
2274 return gsl_cdf_rayleigh_P (x, sigma);
2278 matrix_eval_IDF_RAYLEIGH (double P, double sigma)
2280 return gsl_cdf_rayleigh_Pinv (P, sigma);
2284 matrix_eval_PDF_RAYLEIGH (double x, double sigma)
2286 return gsl_ran_rayleigh_pdf (x, sigma);
2290 matrix_eval_RV_RAYLEIGH (double sigma)
2292 return gsl_ran_rayleigh (get_rng (), sigma);
2296 matrix_eval_PDF_RTAIL (double x, double a, double sigma)
2298 return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
2302 matrix_eval_RV_RTAIL (double a, double sigma)
2304 return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
2308 matrix_eval_CDF_T (double x, double df)
2310 return gsl_cdf_tdist_P (x, df);
2314 matrix_eval_TCDF (double x, double df)
2316 return matrix_eval_CDF_T (x, df);
2320 matrix_eval_IDF_T (double P, double df)
2322 return gsl_cdf_tdist_Pinv (P, df);
2326 matrix_eval_PDF_T (double x, double df)
2328 return gsl_ran_tdist_pdf (x, df);
2332 matrix_eval_RV_T (double df)
2334 return gsl_ran_tdist (get_rng (), df);
2338 matrix_eval_CDF_T1G (double x, double a, double b)
2340 return gsl_cdf_gumbel1_P (x, a, b);
2344 matrix_eval_IDF_T1G (double P, double a, double b)
2346 return gsl_cdf_gumbel1_Pinv (P, a, b);
2350 matrix_eval_PDF_T1G (double x, double a, double b)
2352 return gsl_ran_gumbel1_pdf (x, a, b);
2356 matrix_eval_RV_T1G (double a, double b)
2358 return gsl_ran_gumbel1 (get_rng (), a, b);
2362 matrix_eval_CDF_T2G (double x, double a, double b)
2364 return gsl_cdf_gumbel1_P (x, a, b);
2368 matrix_eval_IDF_T2G (double P, double a, double b)
2370 return gsl_cdf_gumbel1_Pinv (P, a, b);
2374 matrix_eval_PDF_T2G (double x, double a, double b)
2376 return gsl_ran_gumbel1_pdf (x, a, b);
2380 matrix_eval_RV_T2G (double a, double b)
2382 return gsl_ran_gumbel1 (get_rng (), a, b);
2386 matrix_eval_CDF_UNIFORM (double x, double a, double b)
2388 return gsl_cdf_flat_P (x, a, b);
2392 matrix_eval_IDF_UNIFORM (double P, double a, double b)
2394 return gsl_cdf_flat_Pinv (P, a, b);
2398 matrix_eval_PDF_UNIFORM (double x, double a, double b)
2400 return gsl_ran_flat_pdf (x, a, b);
2404 matrix_eval_RV_UNIFORM (double a, double b)
2406 return gsl_ran_flat (get_rng (), a, b);
2410 matrix_eval_CDF_WEIBULL (double x, double a, double b)
2412 return gsl_cdf_weibull_P (x, a, b);
2416 matrix_eval_IDF_WEIBULL (double P, double a, double b)
2418 return gsl_cdf_weibull_Pinv (P, a, b);
2422 matrix_eval_PDF_WEIBULL (double x, double a, double b)
2424 return gsl_ran_weibull_pdf (x, a, b);
2428 matrix_eval_RV_WEIBULL (double a, double b)
2430 return gsl_ran_weibull (get_rng (), a, b);
2434 matrix_eval_CDF_BERNOULLI (double k, double p)
2436 return k ? 1 : 1 - p;
2440 matrix_eval_PDF_BERNOULLI (double k, double p)
2442 return gsl_ran_bernoulli_pdf (k, p);
2446 matrix_eval_RV_BERNOULLI (double p)
2448 return gsl_ran_bernoulli (get_rng (), p);
2452 matrix_eval_CDF_BINOM (double k, double n, double p)
2454 return gsl_cdf_binomial_P (k, p, n);
2458 matrix_eval_PDF_BINOM (double k, double n, double p)
2460 return gsl_ran_binomial_pdf (k, p, n);
2464 matrix_eval_RV_BINOM (double n, double p)
2466 return gsl_ran_binomial (get_rng (), p, n);
2470 matrix_eval_CDF_GEOM (double k, double p)
2472 return gsl_cdf_geometric_P (k, p);
2476 matrix_eval_PDF_GEOM (double k, double p)
2478 return gsl_ran_geometric_pdf (k, p);
2482 matrix_eval_RV_GEOM (double p)
2484 return gsl_ran_geometric (get_rng (), p);
2488 matrix_eval_CDF_HYPER (double k, double a, double b, double c)
2490 return gsl_cdf_hypergeometric_P (k, c, a - c, b);
2494 matrix_eval_PDF_HYPER (double k, double a, double b, double c)
2496 return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
2500 matrix_eval_RV_HYPER (double a, double b, double c)
2502 return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
2506 matrix_eval_PDF_LOG (double k, double p)
2508 return gsl_ran_logarithmic_pdf (k, p);
2512 matrix_eval_RV_LOG (double p)
2514 return gsl_ran_logarithmic (get_rng (), p);
2518 matrix_eval_CDF_NEGBIN (double k, double n, double p)
2520 return gsl_cdf_negative_binomial_P (k, p, n);
2524 matrix_eval_PDF_NEGBIN (double k, double n, double p)
2526 return gsl_ran_negative_binomial_pdf (k, p, n);
2530 matrix_eval_RV_NEGBIN (double n, double p)
2532 return gsl_ran_negative_binomial (get_rng (), p, n);
2536 matrix_eval_CDF_POISSON (double k, double mu)
2538 return gsl_cdf_poisson_P (k, mu);
2542 matrix_eval_PDF_POISSON (double k, double mu)
2544 return gsl_ran_poisson_pdf (k, mu);
2548 matrix_eval_RV_POISSON (double mu)
2550 return gsl_ran_poisson (get_rng (), mu);
2553 struct matrix_function
2557 size_t min_args, max_args;
2560 static struct matrix_expr *matrix_parse_expr (struct matrix_state *);
2563 word_matches (const char **test, const char **name)
2565 size_t test_len = strcspn (*test, ".");
2566 size_t name_len = strcspn (*name, ".");
2567 if (test_len == name_len)
2569 if (buf_compare_case (*test, *name, test_len))
2572 else if (test_len < 3 || test_len > name_len)
2576 if (buf_compare_case (*test, *name, test_len))
2582 if (**test != **name)
2593 /* Returns 0 if TOKEN and FUNC do not match,
2594 1 if TOKEN is an acceptable abbreviation for FUNC,
2595 2 if TOKEN equals FUNC. */
2597 compare_function_names (const char *token_, const char *func_)
2599 const char *token = token_;
2600 const char *func = func_;
2601 while (*token || *func)
2602 if (!word_matches (&token, &func))
2604 return !c_strcasecmp (token_, func_) ? 2 : 1;
2607 static const struct matrix_function *
2608 matrix_parse_function_name (const char *token)
2610 static const struct matrix_function functions[] =
2612 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
2613 { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
2617 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
2619 for (size_t i = 0; i < N_FUNCTIONS; i++)
2621 if (compare_function_names (token, functions[i].name) > 0)
2622 return &functions[i];
2627 static struct read_file *
2628 read_file_create (struct matrix_state *s, struct file_handle *fh)
2630 for (size_t i = 0; i < s->n_read_files; i++)
2632 struct read_file *rf = s->read_files[i];
2640 struct read_file *rf = xmalloc (sizeof *rf);
2641 *rf = (struct read_file) { .file = fh };
2643 s->read_files = xrealloc (s->read_files,
2644 (s->n_read_files + 1) * sizeof *s->read_files);
2645 s->read_files[s->n_read_files++] = rf;
2649 static struct dfm_reader *
2650 read_file_open (struct read_file *rf)
2653 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
2658 read_file_destroy (struct read_file *rf)
2662 fh_unref (rf->file);
2663 dfm_close_reader (rf->reader);
2664 free (rf->encoding);
2669 static struct write_file *
2670 write_file_create (struct matrix_state *s, struct file_handle *fh)
2672 for (size_t i = 0; i < s->n_write_files; i++)
2674 struct write_file *wf = s->write_files[i];
2682 struct write_file *wf = xmalloc (sizeof *wf);
2683 *wf = (struct write_file) { .file = fh };
2685 s->write_files = xrealloc (s->write_files,
2686 (s->n_write_files + 1) * sizeof *s->write_files);
2687 s->write_files[s->n_write_files++] = wf;
2691 static struct dfm_writer *
2692 write_file_open (struct write_file *wf)
2695 wf->writer = dfm_open_writer (wf->file, wf->encoding);
2700 write_file_destroy (struct write_file *wf)
2706 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
2707 wf->held->s.ss.length);
2708 u8_line_destroy (wf->held);
2712 fh_unref (wf->file);
2713 dfm_close_writer (wf->writer);
2714 free (wf->encoding);
2720 matrix_parse_function (struct matrix_state *s, const char *token,
2721 struct matrix_expr **exprp)
2724 if (lex_next_token (s->lexer, 1) != T_LPAREN)
2727 int start_ofs = lex_ofs (s->lexer);
2728 if (lex_match_id (s->lexer, "EOF"))
2731 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
2735 if (!lex_force_match (s->lexer, T_RPAREN))
2741 struct read_file *rf = read_file_create (s, fh);
2743 struct matrix_expr *e = xmalloc (sizeof *e);
2744 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
2745 matrix_expr_wrap_location (s, start_ofs, e);
2750 const struct matrix_function *f = matrix_parse_function_name (token);
2754 struct matrix_expr *e = xmalloc (sizeof *e);
2755 *e = (struct matrix_expr) { .op = f->op };
2757 lex_get_n (s->lexer, 2);
2758 if (lex_token (s->lexer) != T_RPAREN)
2760 size_t allocated_subs = 0;
2763 struct matrix_expr *sub = matrix_parse_expr (s);
2767 if (e->n_subs >= allocated_subs)
2768 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
2769 e->subs[e->n_subs++] = sub;
2771 while (lex_match (s->lexer, T_COMMA));
2773 if (!lex_force_match (s->lexer, T_RPAREN))
2776 if (e->n_subs < f->min_args || e->n_subs > f->max_args)
2778 if (f->min_args == f->max_args)
2779 msg_at (SE, e->location,
2780 ngettext ("Matrix function %s requires %zu argument.",
2781 "Matrix function %s requires %zu arguments.",
2783 f->name, f->min_args);
2784 else if (f->min_args == 1 && f->max_args == 2)
2785 msg_at (SE, e->location,
2786 ngettext ("Matrix function %s requires 1 or 2 arguments, "
2787 "but %zu was provided.",
2788 "Matrix function %s requires 1 or 2 arguments, "
2789 "but %zu were provided.",
2791 f->name, e->n_subs);
2792 else if (f->min_args == 1 && f->max_args == INT_MAX)
2793 msg_at (SE, e->location,
2794 _("Matrix function %s requires at least one argument."),
2802 matrix_expr_wrap_location (s, start_ofs, e);
2808 matrix_expr_destroy (e);
2812 static struct matrix_expr *
2813 matrix_parse_primary__ (struct matrix_state *s)
2815 if (lex_is_number (s->lexer))
2816 return matrix_expr_create_number (s->lexer, lex_number (s->lexer));
2817 else if (lex_is_string (s->lexer))
2819 char string[sizeof (double)];
2820 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
2823 memcpy (&number, string, sizeof number);
2824 return matrix_expr_create_number (s->lexer, number);
2826 else if (lex_match (s->lexer, T_LPAREN))
2828 struct matrix_expr *e = matrix_parse_expr (s);
2829 if (!e || !lex_force_match (s->lexer, T_RPAREN))
2831 matrix_expr_destroy (e);
2836 else if (lex_match (s->lexer, T_LCURLY))
2838 struct matrix_expr *e = matrix_parse_curly_semi (s);
2839 if (!e || !lex_force_match (s->lexer, T_RCURLY))
2841 matrix_expr_destroy (e);
2846 else if (lex_token (s->lexer) == T_ID)
2848 struct matrix_expr *retval;
2849 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
2852 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
2855 lex_error (s->lexer, _("Unknown variable %s."),
2856 lex_tokcstr (s->lexer));
2861 struct matrix_expr *e = xmalloc (sizeof *e);
2862 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
2865 else if (lex_token (s->lexer) == T_ALL)
2867 struct matrix_expr *retval;
2868 if (matrix_parse_function (s, "ALL", &retval))
2872 lex_error (s->lexer, NULL);
2876 static struct matrix_expr *
2877 matrix_parse_primary (struct matrix_state *s)
2879 int start_ofs = lex_ofs (s->lexer);
2880 return matrix_expr_wrap_location (s, start_ofs, matrix_parse_primary__ (s));
2883 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
2886 matrix_parse_index_expr (struct matrix_state *s, struct matrix_expr **indexp)
2888 if (lex_match (s->lexer, T_COLON))
2895 *indexp = matrix_parse_expr (s);
2896 return *indexp != NULL;
2900 static struct matrix_expr *
2901 matrix_parse_postfix (struct matrix_state *s)
2903 struct matrix_expr *lhs = matrix_parse_primary (s);
2904 if (!lhs || !lex_match (s->lexer, T_LPAREN))
2907 struct matrix_expr *i0;
2908 if (!matrix_parse_index_expr (s, &i0))
2910 matrix_expr_destroy (lhs);
2913 if (lex_match (s->lexer, T_RPAREN))
2915 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
2916 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
2917 else if (lex_match (s->lexer, T_COMMA))
2919 struct matrix_expr *i1;
2920 if (!matrix_parse_index_expr (s, &i1)
2921 || !lex_force_match (s->lexer, T_RPAREN))
2923 matrix_expr_destroy (lhs);
2924 matrix_expr_destroy (i0);
2925 matrix_expr_destroy (i1);
2928 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
2929 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
2930 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
2935 lex_error_expecting (s->lexer, "`)'", "`,'");
2940 static struct matrix_expr *
2941 matrix_parse_unary (struct matrix_state *s)
2943 int start_ofs = lex_ofs (s->lexer);
2945 struct matrix_expr *e;
2946 if (lex_match (s->lexer, T_DASH))
2948 struct matrix_expr *sub = matrix_parse_unary (s);
2951 e = matrix_expr_create_1 (MOP_NEGATE, sub);
2953 else if (lex_match (s->lexer, T_PLUS))
2955 e = matrix_parse_unary (s);
2960 return matrix_parse_postfix (s);
2962 matrix_expr_wrap_location (s, start_ofs, e);
2963 e->location->p[0] = lex_ofs_start_point (s->lexer, start_ofs);
2967 static struct matrix_expr *
2968 matrix_parse_seq (struct matrix_state *s)
2970 struct matrix_expr *start = matrix_parse_unary (s);
2971 if (!start || !lex_match (s->lexer, T_COLON))
2974 struct matrix_expr *end = matrix_parse_unary (s);
2977 matrix_expr_destroy (start);
2981 if (lex_match (s->lexer, T_COLON))
2983 struct matrix_expr *increment = matrix_parse_unary (s);
2986 matrix_expr_destroy (start);
2987 matrix_expr_destroy (end);
2990 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
2993 return matrix_expr_create_2 (MOP_SEQ, start, end);
2996 static struct matrix_expr *
2997 matrix_parse_exp (struct matrix_state *s)
2999 struct matrix_expr *lhs = matrix_parse_seq (s);
3006 if (lex_match (s->lexer, T_EXP))
3008 else if (lex_match_phrase (s->lexer, "&**"))
3013 struct matrix_expr *rhs = matrix_parse_seq (s);
3016 matrix_expr_destroy (lhs);
3019 lhs = matrix_expr_create_2 (op, lhs, rhs);
3023 static struct matrix_expr *
3024 matrix_parse_mul_div (struct matrix_state *s)
3026 struct matrix_expr *lhs = matrix_parse_exp (s);
3033 if (lex_match (s->lexer, T_ASTERISK))
3035 else if (lex_match (s->lexer, T_SLASH))
3037 else if (lex_match_phrase (s->lexer, "&*"))
3039 else if (lex_match_phrase (s->lexer, "&/"))
3044 struct matrix_expr *rhs = matrix_parse_exp (s);
3047 matrix_expr_destroy (lhs);
3050 lhs = matrix_expr_create_2 (op, lhs, rhs);
3054 static struct matrix_expr *
3055 matrix_parse_add_sub (struct matrix_state *s)
3057 struct matrix_expr *lhs = matrix_parse_mul_div (s);
3064 if (lex_match (s->lexer, T_PLUS))
3066 else if (lex_match (s->lexer, T_DASH))
3068 else if (lex_token (s->lexer) == T_NEG_NUM)
3073 struct matrix_expr *rhs = matrix_parse_mul_div (s);
3076 matrix_expr_destroy (lhs);
3079 lhs = matrix_expr_create_2 (op, lhs, rhs);
3083 static struct matrix_expr *
3084 matrix_parse_relational (struct matrix_state *s)
3086 struct matrix_expr *lhs = matrix_parse_add_sub (s);
3093 if (lex_match (s->lexer, T_GT))
3095 else if (lex_match (s->lexer, T_GE))
3097 else if (lex_match (s->lexer, T_LT))
3099 else if (lex_match (s->lexer, T_LE))
3101 else if (lex_match (s->lexer, T_EQUALS) || lex_match (s->lexer, T_EQ))
3103 else if (lex_match (s->lexer, T_NE))
3108 struct matrix_expr *rhs = matrix_parse_add_sub (s);
3111 matrix_expr_destroy (lhs);
3114 lhs = matrix_expr_create_2 (op, lhs, rhs);
3118 static struct matrix_expr *
3119 matrix_parse_not (struct matrix_state *s)
3121 int start_ofs = lex_ofs (s->lexer);
3122 if (lex_match (s->lexer, T_NOT))
3124 struct matrix_expr *sub = matrix_parse_not (s);
3128 matrix_expr_wrap_location (s, start_ofs, sub);
3129 sub->location->p[0] = lex_ofs_start_point (s->lexer, start_ofs);
3133 return matrix_parse_relational (s);
3136 static struct matrix_expr *
3137 matrix_parse_and (struct matrix_state *s)
3139 struct matrix_expr *lhs = matrix_parse_not (s);
3143 while (lex_match (s->lexer, T_AND))
3145 struct matrix_expr *rhs = matrix_parse_not (s);
3148 matrix_expr_destroy (lhs);
3151 lhs = matrix_expr_create_2 (MOP_AND, lhs, rhs);
3156 static struct matrix_expr *
3157 matrix_parse_expr__ (struct matrix_state *s)
3159 struct matrix_expr *lhs = matrix_parse_and (s);
3166 if (lex_match (s->lexer, T_OR))
3168 else if (lex_match_id (s->lexer, "XOR"))
3175 struct matrix_expr *rhs = matrix_parse_and (s);
3178 matrix_expr_destroy (lhs);
3181 lhs = matrix_expr_create_2 (op, lhs, rhs);
3185 static struct matrix_expr *
3186 matrix_parse_expr (struct matrix_state *s)
3188 int start_ofs = lex_ofs (s->lexer);
3189 return matrix_expr_wrap_location (s, start_ofs, matrix_parse_expr__ (s));;
3192 /* Expression evaluation. */
3195 matrix_op_eval (enum matrix_op op, double a, double b)
3199 case MOP_ADD_ELEMS: return a + b;
3200 case MOP_SUB_ELEMS: return a - b;
3201 case MOP_MUL_ELEMS: return a * b;
3202 case MOP_DIV_ELEMS: return a / b;
3203 case MOP_EXP_ELEMS: return pow (a, b);
3204 case MOP_GT: return a > b;
3205 case MOP_GE: return a >= b;
3206 case MOP_LT: return a < b;
3207 case MOP_LE: return a <= b;
3208 case MOP_EQ: return a == b;
3209 case MOP_NE: return a != b;
3210 case MOP_AND: return (a > 0) && (b > 0);
3211 case MOP_OR: return (a > 0) || (b > 0);
3212 case MOP_XOR: return (a > 0) != (b > 0);
3214 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3223 case MOP_PASTE_HORZ:
3224 case MOP_PASTE_VERT:
3240 matrix_op_name (enum matrix_op op)
3244 case MOP_ADD_ELEMS: return "+";
3245 case MOP_SUB_ELEMS: return "-";
3246 case MOP_MUL_ELEMS: return "&*";
3247 case MOP_DIV_ELEMS: return "&/";
3248 case MOP_EXP_ELEMS: return "&**";
3249 case MOP_GT: return ">";
3250 case MOP_GE: return ">=";
3251 case MOP_LT: return "<";
3252 case MOP_LE: return "<=";
3253 case MOP_EQ: return "=";
3254 case MOP_NE: return "<>";
3255 case MOP_AND: return "AND";
3256 case MOP_OR: return "OR";
3257 case MOP_XOR: return "XOR";
3259 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3268 case MOP_PASTE_HORZ:
3269 case MOP_PASTE_VERT:
3285 is_scalar (const gsl_matrix *m)
3287 return m->size1 == 1 && m->size2 == 1;
3291 to_scalar (const gsl_matrix *m)
3293 assert (is_scalar (m));
3294 return gsl_matrix_get (m, 0, 0);
3298 matrix_expr_evaluate_elementwise (const struct matrix_expr *e,
3300 gsl_matrix *a, gsl_matrix *b)
3304 double be = to_scalar (b);
3305 for (size_t r = 0; r < a->size1; r++)
3306 for (size_t c = 0; c < a->size2; c++)
3308 double *ae = gsl_matrix_ptr (a, r, c);
3309 *ae = matrix_op_eval (op, *ae, be);
3313 else if (is_scalar (a))
3315 double ae = to_scalar (a);
3316 for (size_t r = 0; r < b->size1; r++)
3317 for (size_t c = 0; c < b->size2; c++)
3319 double *be = gsl_matrix_ptr (b, r, c);
3320 *be = matrix_op_eval (op, ae, *be);
3324 else if (a->size1 == b->size1 && a->size2 == b->size2)
3326 for (size_t r = 0; r < a->size1; r++)
3327 for (size_t c = 0; c < a->size2; c++)
3329 double *ae = gsl_matrix_ptr (a, r, c);
3330 double be = gsl_matrix_get (b, r, c);
3331 *ae = matrix_op_eval (op, *ae, be);
3337 msg_at (SE, e->location,
3338 _("Operands to %s must have the same dimensions or one "
3339 "must be a scalar, not %zu×%zu and %zu×%zu matrices."),
3340 matrix_op_name (op), a->size1, a->size2, b->size1, b->size2);
3346 matrix_expr_evaluate_mul_mat (const struct matrix_expr *e,
3347 gsl_matrix *a, gsl_matrix *b)
3349 if (is_scalar (a) || is_scalar (b))
3350 return matrix_expr_evaluate_elementwise (e, MOP_MUL_ELEMS, a, b);
3352 if (a->size2 != b->size1)
3354 msg_at (SE, e->location,
3355 _("Matrices with dimensions %zu×%zu and %zu×%zu are "
3356 "not conformable for multiplication."),
3357 a->size1, a->size2, b->size1, b->size2);
3361 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3362 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3367 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3369 gsl_matrix *tmp = *a;
3375 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3378 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3379 swap_matrix (z, tmp);
3383 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3385 mul_matrix (x, *x, *x, tmp);
3389 matrix_expr_evaluate_exp_mat (const struct matrix_expr *e,
3390 gsl_matrix *x_, gsl_matrix *b)
3393 if (x->size1 != x->size2)
3395 msg_at (SE, matrix_expr_location (e->subs[0]),
3396 _("Matrix exponentation with ** requires a square matrix on "
3397 "the left-hand size, not one with dimensions %zu×%zu."),
3398 x->size1, x->size2);
3403 msg_at (SE, matrix_expr_location (e->subs[1]),
3404 _("Matrix exponentiation with ** requires a scalar on the "
3405 "right-hand side, not a matrix with dimensions %zu×%zu."),
3406 b->size1, b->size2);
3409 double bf = to_scalar (b);
3410 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3412 msg_at (SE, matrix_expr_location (e->subs[1]),
3413 _("Exponent %.1f in matrix multiplication is non-integer "
3414 "or outside the valid range."), bf);
3419 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3421 gsl_matrix_set_identity (y);
3425 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3427 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3430 mul_matrix (&y, x, y, &t);
3431 square_matrix (&x, &t);
3434 square_matrix (&x, &t);
3436 mul_matrix (&y, x, y, &t);
3440 /* Garbage collection.
3442 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3443 a permutation of them. We are returning one of them; that one must not be
3444 destroyed. We must not destroy 'x_' because the caller owns it. */
3446 gsl_matrix_free (y_);
3448 gsl_matrix_free (t_);
3454 note_nonscalar (gsl_matrix *m, const struct matrix_expr *e)
3457 msg_at (SN, matrix_expr_location (e),
3458 _("This operand is a %zu×%zu matrix."), m->size1, m->size2);
3462 matrix_expr_evaluate_seq (const struct matrix_expr *e,
3463 gsl_matrix *start_, gsl_matrix *end_,
3466 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3468 msg_at (SE, matrix_expr_location (e),
3469 _("All operands of : operator must be scalars."));
3471 note_nonscalar (start_, e->subs[0]);
3472 note_nonscalar (end_, e->subs[1]);
3474 note_nonscalar (by_, e->subs[2]);
3478 long int start = to_scalar (start_);
3479 long int end = to_scalar (end_);
3480 long int by = by_ ? to_scalar (by_) : 1;
3484 msg (SE, _("The increment operand to : must be nonzero."));
3488 long int n = (end >= start && by > 0 ? (end - start + by) / by
3489 : end <= start && by < 0 ? (start - end - by) / -by
3491 gsl_matrix *m = gsl_matrix_alloc (1, n);
3492 for (long int i = 0; i < n; i++)
3493 gsl_matrix_set (m, 0, i, start + i * by);
3498 matrix_expr_evaluate_not (gsl_matrix *a)
3500 for (size_t r = 0; r < a->size1; r++)
3501 for (size_t c = 0; c < a->size2; c++)
3503 double *ae = gsl_matrix_ptr (a, r, c);
3510 matrix_expr_evaluate_paste_horz (gsl_matrix *a, gsl_matrix *b)
3512 if (a->size1 != b->size1)
3514 if (!a->size1 || !a->size2)
3516 else if (!b->size1 || !b->size2)
3519 msg (SE, _("All columns in a matrix must have the same number of rows, "
3520 "but this tries to paste matrices with %zu and %zu rows."),
3521 a->size1, b->size1);
3525 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3526 for (size_t y = 0; y < a->size1; y++)
3528 for (size_t x = 0; x < a->size2; x++)
3529 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3530 for (size_t x = 0; x < b->size2; x++)
3531 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3537 matrix_expr_evaluate_paste_vert (gsl_matrix *a, gsl_matrix *b)
3539 if (a->size2 != b->size2)
3541 if (!a->size1 || !a->size2)
3543 else if (!b->size1 || !b->size2)
3546 msg (SE, _("All rows in a matrix must have the same number of columns, "
3547 "but this tries to stack matrices with %zu and %zu columns."),
3548 a->size2, b->size2);
3552 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3553 for (size_t x = 0; x < a->size2; x++)
3555 for (size_t y = 0; y < a->size1; y++)
3556 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3557 for (size_t y = 0; y < b->size1; y++)
3558 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3564 matrix_to_vector (gsl_matrix *m)
3567 gsl_vector v = to_vector (m);
3568 assert (v.block == m->block || !v.block);
3572 gsl_matrix_free (m);
3573 return xmemdup (&v, sizeof v);
3587 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3590 matrix_normalize_index_vector (const gsl_matrix *m, size_t size,
3591 enum index_type index_type, size_t other_size,
3592 struct index_vector *iv)
3601 msg (SE, _("Vector index must be scalar or vector, not a "
3603 m->size1, m->size2);
3607 msg (SE, _("Matrix row index must be scalar or vector, not a "
3609 m->size1, m->size2);
3613 msg (SE, _("Matrix column index must be scalar or vector, not a "
3615 m->size1, m->size2);
3621 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3622 *iv = (struct index_vector) {
3623 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3626 for (size_t i = 0; i < v.size; i++)
3628 double index = gsl_vector_get (&v, i);
3629 if (index < 1 || index >= size + 1)
3634 msg (SE, _("Index %g is out of range for vector "
3635 "with %zu elements."), index, size);
3639 msg (SE, _("%g is not a valid row index for "
3640 "a %zu×%zu matrix."),
3641 index, size, other_size);
3645 msg (SE, _("%g is not a valid column index for "
3646 "a %zu×%zu matrix."),
3647 index, other_size, size);
3654 iv->indexes[i] = index - 1;
3660 *iv = (struct index_vector) {
3661 .indexes = xnmalloc (size, sizeof *iv->indexes),
3664 for (size_t i = 0; i < size; i++)
3671 matrix_expr_evaluate_vec_all (gsl_matrix *sm)
3673 if (!is_vector (sm))
3675 msg (SE, _("Vector index operator may not be applied to "
3676 "a %zu×%zu matrix."),
3677 sm->size1, sm->size2);
3685 matrix_expr_evaluate_vec_index (gsl_matrix *sm, gsl_matrix *im)
3687 if (!matrix_expr_evaluate_vec_all (sm))
3690 gsl_vector sv = to_vector (sm);
3691 struct index_vector iv;
3692 if (!matrix_normalize_index_vector (im, sv.size, IV_VECTOR, 0, &iv))
3695 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3696 sm->size1 == 1 ? iv.n : 1);
3697 gsl_vector dv = to_vector (dm);
3698 for (size_t dx = 0; dx < iv.n; dx++)
3700 size_t sx = iv.indexes[dx];
3701 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3709 matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0,
3712 struct index_vector iv0;
3713 if (!matrix_normalize_index_vector (im0, sm->size1, IV_ROW, sm->size2, &iv0))
3716 struct index_vector iv1;
3717 if (!matrix_normalize_index_vector (im1, sm->size2, IV_COLUMN, sm->size1,
3724 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3725 for (size_t dy = 0; dy < iv0.n; dy++)
3727 size_t sy = iv0.indexes[dy];
3729 for (size_t dx = 0; dx < iv1.n; dx++)
3731 size_t sx = iv1.indexes[dx];
3732 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3740 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3741 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3742 const struct matrix_function_properties *, gsl_matrix *[], \
3743 const struct matrix_expr *, matrix_proto_##PROTO *);
3748 check_scalar_arg (const char *name, gsl_matrix *subs[], size_t index)
3750 if (!is_scalar (subs[index]))
3752 msg (SE, _("Function %s argument %zu must be a scalar, "
3753 "not a %zu×%zu matrix."),
3754 name, index + 1, subs[index]->size1, subs[index]->size2);
3761 check_vector_arg (const char *name, gsl_matrix *subs[], size_t index)
3763 if (!is_vector (subs[index]))
3765 msg (SE, _("Function %s argument %zu must be a vector, "
3766 "not a %zu×%zu matrix."),
3767 name, index + 1, subs[index]->size1, subs[index]->size2);
3774 to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[])
3776 for (size_t i = 0; i < n_subs; i++)
3778 if (!check_scalar_arg (name, subs, i))
3780 d[i] = to_scalar (subs[i]);
3786 parse_constraint_value (const char **constraintsp)
3789 long retval = strtol (*constraintsp, &tail, 10);
3790 assert (tail > *constraintsp);
3791 *constraintsp = tail;
3796 argument_invalid (const struct matrix_function_properties *props,
3797 size_t arg_index, gsl_matrix *a, size_t y, size_t x,
3801 ds_put_format (out, _("Argument %zu to matrix function %s "
3802 "has invalid value %g."),
3803 arg_index, props->name, gsl_matrix_get (a, y, x));
3805 ds_put_format (out, _("Row %zu, column %zu of argument %zu "
3806 "to matrix function %s has "
3807 "invalid value %g."),
3808 y + 1, x + 1, arg_index, props->name,
3809 gsl_matrix_get (a, y, x));
3812 enum matrix_argument_relop
3822 argument_inequality_error (const struct matrix_function_properties *props,
3823 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3824 size_t b_index, double b,
3825 enum matrix_argument_relop relop)
3827 struct string s = DS_EMPTY_INITIALIZER;
3829 argument_invalid (props, a_index, a, y, x, &s);
3830 ds_put_cstr (&s, " ");
3834 ds_put_format (&s, _("This argument must be greater than or "
3835 "equal to argument %zu, but that has value %g."),
3840 ds_put_format (&s, _("This argument must be greater than argument %zu, "
3841 "but that has value %g."),
3846 ds_put_format (&s, _("This argument must be less than or "
3847 "equal to argument %zu, but that has value %g."),
3852 ds_put_format (&s, _("This argument must be less than argument %zu, "
3853 "but that has value %g."),
3859 _("This argument must not be the same as argument %zu."),
3864 msg (ME, ds_cstr (&s));
3869 argument_value_error (const struct matrix_function_properties *props,
3870 size_t a_index, gsl_matrix *a, size_t y, size_t x,
3872 enum matrix_argument_relop relop)
3874 struct string s = DS_EMPTY_INITIALIZER;
3875 argument_invalid (props, a_index, a, y, x, &s);
3876 ds_put_cstr (&s, " ");
3881 _("This argument must be greater than or equal to %g."),
3886 ds_put_format (&s, _("This argument must be greater than %g."), b);
3890 ds_put_format (&s, _("This argument must be less than or equal to %g."),
3895 ds_put_format (&s, _("This argument must be less than %g."), b);
3899 ds_put_format (&s, _("This argument may not be %g."), b);
3903 msg (ME, ds_cstr (&s));
3908 matrix_argument_relop_is_satisfied (double a, double b,
3909 enum matrix_argument_relop relop)
3913 case MRR_GE: return a >= b;
3914 case MRR_GT: return a > b;
3915 case MRR_LE: return a <= b;
3916 case MRR_LT: return a < b;
3917 case MRR_NE: return a != b;
3923 static enum matrix_argument_relop
3924 matrix_argument_relop_flip (enum matrix_argument_relop relop)
3928 case MRR_GE: return MRR_LE;
3929 case MRR_GT: return MRR_LT;
3930 case MRR_LE: return MRR_GE;
3931 case MRR_LT: return MRR_GT;
3932 case MRR_NE: return MRR_NE;
3939 check_constraints (const struct matrix_function_properties *props,
3940 gsl_matrix *args[], size_t n_args)
3942 const char *constraints = props->constraints;
3946 size_t arg_index = SIZE_MAX;
3947 while (*constraints)
3949 if (*constraints >= 'a' && *constraints <= 'd')
3951 arg_index = *constraints++ - 'a';
3952 assert (arg_index < n_args);
3954 else if (*constraints == '[' || *constraints == '(')
3956 assert (arg_index < n_args);
3957 bool open_lower = *constraints++ == '(';
3958 int minimum = parse_constraint_value (&constraints);
3959 assert (*constraints == ',');
3961 int maximum = parse_constraint_value (&constraints);
3962 assert (*constraints == ']' || *constraints == ')');
3963 bool open_upper = *constraints++ == ')';
3965 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3966 if ((open_lower ? *d <= minimum : *d < minimum)
3967 || (open_upper ? *d >= maximum : *d > maximum))
3969 if (!is_scalar (args[arg_index]))
3970 msg (ME, _("Row %zu, column %zu of argument %zu to matrix "
3971 "function %s has value %g, which is outside "
3972 "the valid range %c%d,%d%c."),
3973 y + 1, x + 1, arg_index + 1, props->name, *d,
3974 open_lower ? '(' : '[',
3976 open_upper ? ')' : ']');
3978 msg (ME, _("Argument %zu to matrix function %s has value %g, "
3979 "which is outside the valid range %c%d,%d%c."),
3980 arg_index + 1, props->name, *d,
3981 open_lower ? '(' : '[',
3983 open_upper ? ')' : ']');
3987 else if (*constraints == '>'
3988 || *constraints == '<'
3989 || *constraints == '!')
3991 enum matrix_argument_relop relop;
3992 switch (*constraints++)
3995 if (*constraints == '=')
4005 if (*constraints == '=')
4015 assert (*constraints == '=');
4020 if (*constraints >= 'a' && *constraints <= 'd')
4022 size_t a_index = arg_index;
4023 size_t b_index = *constraints - 'a';
4024 assert (a_index < n_args);
4025 assert (b_index < n_args);
4027 /* We only support one of the two arguments being non-scalar.
4028 It's easier to support only the first one being non-scalar, so
4029 flip things around if it's the other way. */
4030 if (!is_scalar (args[b_index]))
4032 assert (is_scalar (args[a_index]));
4033 size_t tmp_index = a_index;
4035 b_index = tmp_index;
4036 relop = matrix_argument_relop_flip (relop);
4039 double b = to_scalar (args[b_index]);
4040 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
4041 if (!matrix_argument_relop_is_satisfied (*a, b, relop))
4043 argument_inequality_error (props,
4044 a_index + 1, args[a_index], y, x,
4052 int comparand = parse_constraint_value (&constraints);
4054 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
4055 if (!matrix_argument_relop_is_satisfied (*d, comparand, relop))
4057 argument_value_error (props,
4058 arg_index + 1, args[arg_index], y, x,
4067 assert (*constraints == ' ');
4069 arg_index = SIZE_MAX;
4076 matrix_expr_evaluate_d_none (
4077 const struct matrix_function_properties *props UNUSED,
4078 gsl_matrix *subs[] UNUSED, const struct matrix_expr *e,
4079 matrix_proto_d_none *f)
4081 assert (e->n_subs == 0);
4083 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4084 gsl_matrix_set (m, 0, 0, f ());
4089 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
4090 gsl_matrix *subs[], const struct matrix_expr *e,
4091 matrix_proto_d_d *f)
4093 assert (e->n_subs == 1);
4096 if (!to_scalar_args (props->name, subs, e->n_subs, &d))
4099 if (!check_constraints (props, subs, e->n_subs))
4102 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4103 gsl_matrix_set (m, 0, 0, f (d));
4108 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
4109 gsl_matrix *subs[], const struct matrix_expr *e,
4110 matrix_proto_d_dd *f)
4112 assert (e->n_subs == 2);
4115 if (!to_scalar_args (props->name, subs, e->n_subs, d))
4118 if (!check_constraints (props, subs, e->n_subs))
4121 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4122 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
4127 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
4128 gsl_matrix *subs[], const struct matrix_expr *e,
4129 matrix_proto_d_ddd *f)
4131 assert (e->n_subs == 3);
4134 if (!to_scalar_args (props->name, subs, e->n_subs, d))
4137 if (!check_constraints (props, subs, e->n_subs))
4140 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4141 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
4146 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
4147 gsl_matrix *subs[], const struct matrix_expr *e,
4148 matrix_proto_m_d *f)
4150 assert (e->n_subs == 1);
4153 return to_scalar_args (props->name, subs, e->n_subs, &d) ? f (d) : NULL;
4157 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
4158 gsl_matrix *subs[], const struct matrix_expr *e,
4159 matrix_proto_m_ddd *f)
4161 assert (e->n_subs == 3);
4164 return to_scalar_args (props->name, subs, e->n_subs, d) ? f(d[0], d[1], d[2]) : NULL;
4168 matrix_expr_evaluate_m_ddn (const struct matrix_function_properties *props,
4169 gsl_matrix *subs[], const struct matrix_expr *e,
4170 matrix_proto_m_ddn *f)
4172 assert (e->n_subs == 2);
4175 return (to_scalar_args (props->name, subs, e->n_subs, d)
4181 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED,
4182 gsl_matrix *subs[], const struct matrix_expr *e,
4183 matrix_proto_m_m *f)
4185 assert (e->n_subs == 1);
4190 matrix_expr_evaluate_m_mn (const struct matrix_function_properties *props UNUSED,
4191 gsl_matrix *subs[], const struct matrix_expr *e,
4192 matrix_proto_m_mn *f)
4194 assert (e->n_subs == 1);
4195 return f (subs[0], e);
4199 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
4200 gsl_matrix *subs[], const struct matrix_expr *e,
4201 matrix_proto_m_e *f)
4203 assert (e->n_subs == 1);
4205 if (!check_constraints (props, subs, e->n_subs))
4208 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4214 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED,
4215 gsl_matrix *subs[], const struct matrix_expr *e,
4216 matrix_proto_m_md *f)
4218 assert (e->n_subs == 2);
4219 if (!check_scalar_arg (props->name, subs, 1))
4221 return f (subs[0], to_scalar (subs[1]));
4225 matrix_expr_evaluate_m_mdn (const struct matrix_function_properties *props UNUSED,
4226 gsl_matrix *subs[], const struct matrix_expr *e,
4227 matrix_proto_m_mdn *f)
4229 assert (e->n_subs == 2);
4230 if (!check_scalar_arg (props->name, subs, 1))
4232 return f (subs[0], to_scalar (subs[1]), e);
4236 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
4237 gsl_matrix *subs[], const struct matrix_expr *e,
4238 matrix_proto_m_ed *f)
4240 assert (e->n_subs == 2);
4241 if (!check_scalar_arg (props->name, subs, 1))
4244 if (!check_constraints (props, subs, e->n_subs))
4247 double b = to_scalar (subs[1]);
4248 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4254 matrix_expr_evaluate_m_mddn (const struct matrix_function_properties *props UNUSED,
4255 gsl_matrix *subs[], const struct matrix_expr *e,
4256 matrix_proto_m_mddn *f)
4258 assert (e->n_subs == 3);
4259 if (!check_scalar_arg (props->name, subs, 1)
4260 || !check_scalar_arg (props->name, subs, 2))
4262 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]), e);
4266 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
4267 gsl_matrix *subs[], const struct matrix_expr *e,
4268 matrix_proto_m_edd *f)
4270 assert (e->n_subs == 3);
4271 if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
4274 if (!check_constraints (props, subs, e->n_subs))
4277 double b = to_scalar (subs[1]);
4278 double c = to_scalar (subs[2]);
4279 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4285 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4286 gsl_matrix *subs[], const struct matrix_expr *e,
4287 matrix_proto_m_eddd *f)
4289 assert (e->n_subs == 4);
4290 for (size_t i = 1; i < 4; i++)
4291 if (!check_scalar_arg (props->name, subs, i))
4294 if (!check_constraints (props, subs, e->n_subs))
4297 double b = to_scalar (subs[1]);
4298 double c = to_scalar (subs[2]);
4299 double d = to_scalar (subs[3]);
4300 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4301 *a = f (*a, b, c, d);
4306 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4307 gsl_matrix *subs[], const struct matrix_expr *e,
4308 matrix_proto_m_eed *f)
4310 assert (e->n_subs == 3);
4311 if (!check_scalar_arg (props->name, subs, 2))
4314 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4315 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4317 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
4318 loc->p[1] = e->subs[1]->location->p[1];
4321 _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4322 "%zu×%zu, but %s requires these arguments either to have "
4323 "the same dimensions or for one of them to be a scalar."),
4325 subs[0]->size1, subs[0]->size2,
4326 subs[1]->size1, subs[1]->size2,
4329 msg_location_destroy (loc);
4333 if (!check_constraints (props, subs, e->n_subs))
4336 double c = to_scalar (subs[2]);
4338 if (is_scalar (subs[0]))
4340 double a = to_scalar (subs[0]);
4341 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4347 double b = to_scalar (subs[1]);
4348 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4355 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED,
4356 gsl_matrix *subs[], const struct matrix_expr *e,
4357 matrix_proto_m_mm *f)
4359 assert (e->n_subs == 2);
4360 return f (subs[0], subs[1]);
4364 matrix_expr_evaluate_m_mmn (const struct matrix_function_properties *props UNUSED,
4365 gsl_matrix *subs[], const struct matrix_expr *e,
4366 matrix_proto_m_mmn *f)
4368 assert (e->n_subs == 2);
4369 return f (subs[0], subs[1], e);
4373 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4374 gsl_matrix *subs[], const struct matrix_expr *e,
4375 matrix_proto_m_v *f)
4377 assert (e->n_subs == 1);
4378 if (!check_vector_arg (props->name, subs, 0))
4380 gsl_vector v = to_vector (subs[0]);
4385 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED,
4386 gsl_matrix *subs[], const struct matrix_expr *e,
4387 matrix_proto_d_m *f)
4389 assert (e->n_subs == 1);
4391 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4392 gsl_matrix_set (m, 0, 0, f (subs[0]));
4397 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED,
4398 gsl_matrix *subs[], const struct matrix_expr *e,
4399 matrix_proto_m_any *f)
4401 return f (subs, e->n_subs);
4405 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props,
4406 gsl_matrix *subs[], const struct matrix_expr *e,
4407 matrix_proto_IDENT *f)
4409 assert (e->n_subs <= 2);
4412 if (!to_scalar_args (props->name, subs, e->n_subs, d))
4415 return f (d[0], d[e->n_subs - 1]);
4419 matrix_expr_evaluate (const struct matrix_expr *e)
4421 if (e->op == MOP_NUMBER)
4423 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4424 gsl_matrix_set (m, 0, 0, e->number);
4427 else if (e->op == MOP_VARIABLE)
4429 const gsl_matrix *src = e->variable->value;
4432 msg_at (SE, e->location,
4433 _("Uninitialized variable %s used in expression."),
4438 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4439 gsl_matrix_memcpy (dst, src);
4442 else if (e->op == MOP_EOF)
4444 struct dfm_reader *reader = read_file_open (e->eof);
4445 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4446 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4450 enum { N_LOCAL = 3 };
4451 gsl_matrix *local_subs[N_LOCAL];
4452 gsl_matrix **subs = (e->n_subs < N_LOCAL
4454 : xmalloc (e->n_subs * sizeof *subs));
4456 for (size_t i = 0; i < e->n_subs; i++)
4458 subs[i] = matrix_expr_evaluate (e->subs[i]);
4461 for (size_t j = 0; j < i; j++)
4462 gsl_matrix_free (subs[j]);
4463 if (subs != local_subs)
4469 gsl_matrix *result = NULL;
4472 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4473 case MOP_F_##ENUM: \
4475 static const struct matrix_function_properties props = { \
4477 .constraints = CONSTRAINTS, \
4479 result = matrix_expr_evaluate_##PROTO (&props, subs, e, \
4480 matrix_eval_##ENUM); \
4487 gsl_matrix_scale (subs[0], -1.0);
4505 result = matrix_expr_evaluate_elementwise (e, e->op, subs[0], subs[1]);
4509 result = matrix_expr_evaluate_not (subs[0]);
4513 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], NULL);
4517 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], subs[2]);
4521 result = matrix_expr_evaluate_mul_mat (e, subs[0], subs[1]);
4525 result = matrix_expr_evaluate_exp_mat (e, subs[0], subs[1]);
4528 case MOP_PASTE_HORZ:
4529 result = matrix_expr_evaluate_paste_horz (subs[0], subs[1]);
4532 case MOP_PASTE_VERT:
4533 result = matrix_expr_evaluate_paste_vert (subs[0], subs[1]);
4537 result = gsl_matrix_alloc (0, 0);
4541 result = matrix_expr_evaluate_vec_index (subs[0], subs[1]);
4545 result = matrix_expr_evaluate_vec_all (subs[0]);
4549 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], subs[2]);
4553 result = matrix_expr_evaluate_mat_index (subs[0], subs[1], NULL);
4557 result = matrix_expr_evaluate_mat_index (subs[0], NULL, subs[1]);
4566 for (size_t i = 0; i < e->n_subs; i++)
4567 if (subs[i] != result)
4568 gsl_matrix_free (subs[i]);
4569 if (subs != local_subs)
4575 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4578 gsl_matrix *m = matrix_expr_evaluate (e);
4584 msg (SE, _("Expression for %s must evaluate to scalar, "
4585 "not a %zu×%zu matrix."),
4586 context, m->size1, m->size2);
4587 gsl_matrix_free (m);
4592 gsl_matrix_free (m);
4597 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4601 if (!matrix_expr_evaluate_scalar (e, context, &d))
4605 if (d < LONG_MIN || d > LONG_MAX)
4607 msg (SE, _("Expression for %s is outside the integer range."), context);
4614 struct matrix_lvalue
4616 struct matrix_var *var;
4617 struct matrix_expr *indexes[2];
4620 struct msg_location *var_location; /* Variable name. */
4621 struct msg_location *index_location; /* Variable name plus indexing. */
4625 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4629 msg_location_destroy (lvalue->var_location);
4630 msg_location_destroy (lvalue->index_location);
4631 for (size_t i = 0; i < lvalue->n_indexes; i++)
4632 matrix_expr_destroy (lvalue->indexes[i]);
4637 static struct matrix_lvalue *
4638 matrix_lvalue_parse (struct matrix_state *s)
4640 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4641 if (!lex_force_id (s->lexer))
4643 int start_ofs = lex_ofs (s->lexer);
4644 lvalue->var_location = lex_get_location (s->lexer, 0, 0);
4645 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4646 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4650 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
4654 lex_get_n (s->lexer, 2);
4656 if (!matrix_parse_index_expr (s, &lvalue->indexes[0]))
4658 lvalue->n_indexes++;
4660 if (lex_match (s->lexer, T_COMMA))
4662 if (!matrix_parse_index_expr (s, &lvalue->indexes[1]))
4664 lvalue->n_indexes++;
4666 if (!lex_force_match (s->lexer, T_RPAREN))
4669 lvalue->index_location = lex_ofs_location (s->lexer, start_ofs,
4670 lex_ofs (s->lexer) - 1);
4675 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4681 matrix_lvalue_destroy (lvalue);
4686 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4687 enum index_type index_type, size_t other_size,
4688 struct index_vector *iv)
4693 m = matrix_expr_evaluate (e);
4700 bool ok = matrix_normalize_index_vector (m, size, index_type,
4702 gsl_matrix_free (m);
4707 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4708 struct index_vector *iv, gsl_matrix *sm)
4710 gsl_vector dv = to_vector (lvalue->var->value);
4712 /* Convert source matrix 'sm' to source vector 'sv'. */
4713 if (!is_vector (sm))
4715 msg (SE, _("Can't assign %zu×%zu matrix to subvector."),
4716 sm->size1, sm->size2);
4719 gsl_vector sv = to_vector (sm);
4720 if (iv->n != sv.size)
4722 msg (SE, _("Can't assign %zu-element vector "
4723 "to %zu-element subvector."), sv.size, iv->n);
4727 /* Assign elements. */
4728 for (size_t x = 0; x < iv->n; x++)
4729 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4734 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4735 struct index_vector *iv0,
4736 struct index_vector *iv1, gsl_matrix *sm)
4738 gsl_matrix *dm = lvalue->var->value;
4740 /* Convert source matrix 'sm' to source vector 'sv'. */
4741 if (iv0->n != sm->size1)
4743 msg (SE, _("Row index vector for assignment to %s has %zu elements "
4744 "but source matrix has %zu rows."),
4745 lvalue->var->name, iv0->n, sm->size1);
4748 if (iv1->n != sm->size2)
4750 msg (SE, _("Column index vector for assignment to %s has %zu elements "
4751 "but source matrix has %zu columns."),
4752 lvalue->var->name, iv1->n, sm->size2);
4756 /* Assign elements. */
4757 for (size_t y = 0; y < iv0->n; y++)
4759 size_t f0 = iv0->indexes[y];
4760 for (size_t x = 0; x < iv1->n; x++)
4762 size_t f1 = iv1->indexes[x];
4763 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4769 /* Takes ownership of M. */
4771 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4772 struct index_vector *iv0, struct index_vector *iv1,
4775 if (!lvalue->n_indexes)
4777 gsl_matrix_free (lvalue->var->value);
4778 lvalue->var->value = sm;
4783 bool ok = (lvalue->n_indexes == 1
4784 ? matrix_lvalue_assign_vector (lvalue, iv0, sm)
4785 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm));
4786 free (iv0->indexes);
4787 free (iv1->indexes);
4788 gsl_matrix_free (sm);
4794 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4795 struct index_vector *iv0,
4796 struct index_vector *iv1)
4798 *iv0 = INDEX_VECTOR_INIT;
4799 *iv1 = INDEX_VECTOR_INIT;
4800 if (!lvalue->n_indexes)
4803 /* Validate destination matrix exists and has the right shape. */
4804 gsl_matrix *dm = lvalue->var->value;
4807 msg_at (SE, lvalue->var_location,
4808 _("Undefined variable %s."), lvalue->var->name);
4811 else if (dm->size1 == 0 || dm->size2 == 0)
4813 msg_at (SE, lvalue->index_location,
4814 _("Cannot index %zu×%zu matrix %s."),
4815 dm->size1, dm->size2, lvalue->var->name);
4818 else if (lvalue->n_indexes == 1)
4820 if (!is_vector (dm))
4822 msg_at (SE, lvalue->index_location,
4823 _("Can't use vector indexing on %zu×%zu matrix %s."),
4824 dm->size1, dm->size2, lvalue->var->name);
4827 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4828 MAX (dm->size1, dm->size2),
4833 assert (lvalue->n_indexes == 2);
4834 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4835 IV_ROW, dm->size2, iv0))
4838 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4839 IV_COLUMN, dm->size1, iv1))
4841 free (iv0->indexes);
4848 /* Takes ownership of M. */
4850 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue, gsl_matrix *sm)
4852 struct index_vector iv0, iv1;
4853 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
4855 gsl_matrix_free (sm);
4859 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm);
4863 /* Matrix command. */
4867 struct matrix_cmd **commands;
4871 static void matrix_cmds_uninit (struct matrix_cmds *);
4875 enum matrix_cmd_type
4898 struct compute_command
4900 struct matrix_lvalue *lvalue;
4901 struct matrix_expr *rvalue;
4905 struct print_command
4907 struct matrix_expr *expression;
4908 bool use_default_format;
4909 struct fmt_spec format;
4911 int space; /* -1 means NEWPAGE. */
4915 struct string_array literals; /* CLABELS/RLABELS. */
4916 struct matrix_expr *expr; /* CNAMES/RNAMES. */
4922 struct do_if_command
4926 struct matrix_expr *condition;
4927 struct matrix_cmds commands;
4937 struct matrix_var *var;
4938 struct matrix_expr *start;
4939 struct matrix_expr *end;
4940 struct matrix_expr *increment;
4942 /* Loop conditions. */
4943 struct matrix_expr *top_condition;
4944 struct matrix_expr *bottom_condition;
4947 struct matrix_cmds commands;
4951 struct display_command
4953 struct matrix_state *state;
4957 struct release_command
4959 struct matrix_var **vars;
4966 struct matrix_expr *expression;
4967 struct save_file *sf;
4973 struct read_file *rf;
4974 struct matrix_lvalue *dst;
4975 struct matrix_expr *size;
4977 enum fmt_type format;
4984 struct write_command
4986 struct write_file *wf;
4987 struct matrix_expr *expression;
4990 /* If this is nonnull, WRITE uses this format.
4992 If this is NULL, WRITE uses free-field format with as many
4993 digits of precision as needed. */
4994 struct fmt_spec *format;
5003 struct matrix_lvalue *dst;
5004 struct dataset *dataset;
5005 struct file_handle *file;
5007 struct var_syntax *vars;
5009 struct matrix_var *names;
5011 /* Treatment of missing values. */
5016 MGET_ERROR, /* Flag error on command. */
5017 MGET_ACCEPT, /* Accept without change, user-missing only. */
5018 MGET_OMIT, /* Drop this case. */
5019 MGET_RECODE /* Recode to 'substitute'. */
5022 double substitute; /* MGET_RECODE only. */
5028 struct msave_command
5030 struct msave_common *common;
5031 struct matrix_expr *expr;
5032 const char *rowtype;
5033 const struct matrix_expr *factors;
5034 const struct matrix_expr *splits;
5040 struct matrix_state *state;
5041 struct file_handle *file;
5043 struct stringi_set rowtypes;
5047 struct eigen_command
5049 struct matrix_expr *expr;
5050 struct matrix_var *evec;
5051 struct matrix_var *eval;
5055 struct setdiag_command
5057 struct matrix_var *dst;
5058 struct matrix_expr *expr;
5064 struct matrix_expr *expr;
5065 struct matrix_var *u;
5066 struct matrix_var *s;
5067 struct matrix_var *v;
5073 static struct matrix_cmd *matrix_parse_command (struct matrix_state *);
5074 static bool matrix_cmd_execute (struct matrix_cmd *);
5075 static void matrix_cmd_destroy (struct matrix_cmd *);
5078 static struct matrix_cmd *
5079 matrix_parse_compute (struct matrix_state *s)
5081 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5082 *cmd = (struct matrix_cmd) {
5083 .type = MCMD_COMPUTE,
5084 .compute = { .lvalue = NULL },
5087 struct compute_command *compute = &cmd->compute;
5088 compute->lvalue = matrix_lvalue_parse (s);
5089 if (!compute->lvalue)
5092 if (!lex_force_match (s->lexer, T_EQUALS))
5095 compute->rvalue = matrix_parse_expr (s);
5096 if (!compute->rvalue)
5102 matrix_cmd_destroy (cmd);
5107 matrix_cmd_execute_compute (struct compute_command *compute)
5109 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
5113 matrix_lvalue_evaluate_and_assign (compute->lvalue, value);
5117 print_labels_destroy (struct print_labels *labels)
5121 string_array_destroy (&labels->literals);
5122 matrix_expr_destroy (labels->expr);
5128 parse_literal_print_labels (struct matrix_state *s,
5129 struct print_labels **labelsp)
5131 lex_match (s->lexer, T_EQUALS);
5132 print_labels_destroy (*labelsp);
5133 *labelsp = xzalloc (sizeof **labelsp);
5134 while (lex_token (s->lexer) != T_SLASH
5135 && lex_token (s->lexer) != T_ENDCMD
5136 && lex_token (s->lexer) != T_STOP)
5138 struct string label = DS_EMPTY_INITIALIZER;
5139 while (lex_token (s->lexer) != T_COMMA
5140 && lex_token (s->lexer) != T_SLASH
5141 && lex_token (s->lexer) != T_ENDCMD
5142 && lex_token (s->lexer) != T_STOP)
5144 if (!ds_is_empty (&label))
5145 ds_put_byte (&label, ' ');
5147 if (lex_token (s->lexer) == T_STRING)
5148 ds_put_cstr (&label, lex_tokcstr (s->lexer));
5151 char *rep = lex_next_representation (s->lexer, 0, 0);
5152 ds_put_cstr (&label, rep);
5157 string_array_append_nocopy (&(*labelsp)->literals,
5158 ds_steal_cstr (&label));
5160 if (!lex_match (s->lexer, T_COMMA))
5166 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
5168 lex_match (s->lexer, T_EQUALS);
5169 struct matrix_expr *e = matrix_parse_exp (s);
5173 print_labels_destroy (*labelsp);
5174 *labelsp = xzalloc (sizeof **labelsp);
5175 (*labelsp)->expr = e;
5179 static struct matrix_cmd *
5180 matrix_parse_print (struct matrix_state *s)
5182 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5183 *cmd = (struct matrix_cmd) {
5186 .use_default_format = true,
5190 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
5193 for (size_t i = 0; ; i++)
5195 enum token_type t = lex_next_token (s->lexer, i);
5196 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
5198 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
5200 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
5203 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
5208 cmd->print.expression = matrix_parse_exp (s);
5209 if (!cmd->print.expression)
5213 while (lex_match (s->lexer, T_SLASH))
5215 if (lex_match_id (s->lexer, "FORMAT"))
5217 lex_match (s->lexer, T_EQUALS);
5218 if (!parse_format_specifier (s->lexer, &cmd->print.format))
5220 cmd->print.use_default_format = false;
5222 else if (lex_match_id (s->lexer, "TITLE"))
5224 lex_match (s->lexer, T_EQUALS);
5225 if (!lex_force_string (s->lexer))
5227 free (cmd->print.title);
5228 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
5231 else if (lex_match_id (s->lexer, "SPACE"))
5233 lex_match (s->lexer, T_EQUALS);
5234 if (lex_match_id (s->lexer, "NEWPAGE"))
5235 cmd->print.space = -1;
5236 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
5238 cmd->print.space = lex_integer (s->lexer);
5244 else if (lex_match_id (s->lexer, "RLABELS"))
5245 parse_literal_print_labels (s, &cmd->print.rlabels);
5246 else if (lex_match_id (s->lexer, "CLABELS"))
5247 parse_literal_print_labels (s, &cmd->print.clabels);
5248 else if (lex_match_id (s->lexer, "RNAMES"))
5250 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
5253 else if (lex_match_id (s->lexer, "CNAMES"))
5255 if (!parse_expr_print_labels (s, &cmd->print.clabels))
5260 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
5261 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
5269 matrix_cmd_destroy (cmd);
5274 matrix_is_integer (const gsl_matrix *m)
5276 for (size_t y = 0; y < m->size1; y++)
5277 for (size_t x = 0; x < m->size2; x++)
5279 double d = gsl_matrix_get (m, y, x);
5287 matrix_max_magnitude (const gsl_matrix *m)
5290 for (size_t y = 0; y < m->size1; y++)
5291 for (size_t x = 0; x < m->size2; x++)
5293 double d = fabs (gsl_matrix_get (m, y, x));
5301 format_fits (struct fmt_spec format, double x)
5303 char *s = data_out (&(union value) { .f = x }, NULL,
5304 &format, settings_get_fmt_settings ());
5305 bool fits = *s != '*' && !strchr (s, 'E');
5310 static struct fmt_spec
5311 default_format (const gsl_matrix *m, int *log_scale)
5315 double max = matrix_max_magnitude (m);
5317 if (matrix_is_integer (m))
5318 for (int w = 1; w <= 12; w++)
5320 struct fmt_spec format = { .type = FMT_F, .w = w };
5321 if (format_fits (format, -max))
5325 if (max >= 1e9 || max <= 1e-4)
5328 snprintf (s, sizeof s, "%.1e", max);
5330 const char *e = strchr (s, 'e');
5332 *log_scale = atoi (e + 1);
5335 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5339 trimmed_string (double d)
5341 struct substring s = ss_buffer ((char *) &d, sizeof d);
5342 ss_rtrim (&s, ss_cstr (" "));
5343 return ss_xstrdup (s);
5346 static struct string_array *
5347 print_labels_get (const struct print_labels *labels, size_t n,
5348 const char *prefix, bool truncate)
5353 struct string_array *out = xzalloc (sizeof *out);
5354 if (labels->literals.n)
5355 string_array_clone (out, &labels->literals);
5356 else if (labels->expr)
5358 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5359 if (m && is_vector (m))
5361 gsl_vector v = to_vector (m);
5362 for (size_t i = 0; i < v.size; i++)
5363 string_array_append_nocopy (out, trimmed_string (
5364 gsl_vector_get (&v, i)));
5366 gsl_matrix_free (m);
5372 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5375 string_array_append (out, "");
5378 string_array_delete (out, out->n - 1);
5381 for (size_t i = 0; i < out->n; i++)
5383 char *s = out->strings[i];
5384 s[strnlen (s, 8)] = '\0';
5391 matrix_cmd_print_space (int space)
5394 output_item_submit (page_break_item_create ());
5395 for (int i = 0; i < space; i++)
5396 output_log ("%s", "");
5400 matrix_cmd_print_text (const struct print_command *print, const gsl_matrix *m,
5401 struct fmt_spec format, int log_scale)
5403 matrix_cmd_print_space (print->space);
5405 output_log ("%s", print->title);
5407 output_log (" 10 ** %d X", log_scale);
5409 struct string_array *clabels = print_labels_get (print->clabels,
5410 m->size2, "col", true);
5411 if (clabels && format.w < 8)
5413 struct string_array *rlabels = print_labels_get (print->rlabels,
5414 m->size1, "row", true);
5418 struct string line = DS_EMPTY_INITIALIZER;
5420 ds_put_byte_multiple (&line, ' ', 8);
5421 for (size_t x = 0; x < m->size2; x++)
5422 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5423 output_log_nocopy (ds_steal_cstr (&line));
5426 double scale = pow (10.0, log_scale);
5427 bool numeric = fmt_is_numeric (format.type);
5428 for (size_t y = 0; y < m->size1; y++)
5430 struct string line = DS_EMPTY_INITIALIZER;
5432 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5434 for (size_t x = 0; x < m->size2; x++)
5436 double f = gsl_matrix_get (m, y, x);
5438 ? data_out (&(union value) { .f = f / scale}, NULL,
5439 &format, settings_get_fmt_settings ())
5440 : trimmed_string (f));
5441 ds_put_format (&line, " %s", s);
5444 output_log_nocopy (ds_steal_cstr (&line));
5447 string_array_destroy (rlabels);
5449 string_array_destroy (clabels);
5454 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5455 const struct print_labels *print_labels, size_t n,
5456 const char *name, const char *prefix)
5458 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5460 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5461 for (size_t i = 0; i < n; i++)
5462 pivot_category_create_leaf (
5464 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5465 : pivot_value_new_integer (i + 1)));
5467 d->hide_all_labels = true;
5468 string_array_destroy (labels);
5473 matrix_cmd_print_tables (const struct print_command *print, const gsl_matrix *m,
5474 struct fmt_spec format, int log_scale)
5476 struct pivot_table *table = pivot_table_create__ (
5477 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5480 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5482 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5483 N_("Columns"), "col");
5485 struct pivot_footnote *footnote = NULL;
5488 char *s = xasprintf ("× 10**%d", log_scale);
5489 footnote = pivot_table_create_footnote (
5490 table, pivot_value_new_user_text_nocopy (s));
5493 double scale = pow (10.0, log_scale);
5494 bool numeric = fmt_is_numeric (format.type);
5495 for (size_t y = 0; y < m->size1; y++)
5496 for (size_t x = 0; x < m->size2; x++)
5498 double f = gsl_matrix_get (m, y, x);
5499 struct pivot_value *v;
5502 v = pivot_value_new_number (f / scale);
5503 v->numeric.format = format;
5506 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5508 pivot_value_add_footnote (v, footnote);
5509 pivot_table_put2 (table, y, x, v);
5512 pivot_table_submit (table);
5516 matrix_cmd_execute_print (const struct print_command *print)
5518 if (print->expression)
5520 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5525 struct fmt_spec format = (print->use_default_format
5526 ? default_format (m, &log_scale)
5529 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5530 matrix_cmd_print_text (print, m, format, log_scale);
5532 matrix_cmd_print_tables (print, m, format, log_scale);
5534 gsl_matrix_free (m);
5538 matrix_cmd_print_space (print->space);
5540 output_log ("%s", print->title);
5547 matrix_parse_commands (struct matrix_state *s, struct matrix_cmds *c,
5548 const char *command_name,
5549 const char *stop1, const char *stop2)
5551 lex_end_of_command (s->lexer);
5552 lex_discard_rest_of_command (s->lexer);
5554 size_t allocated = 0;
5557 while (lex_token (s->lexer) == T_ENDCMD)
5560 if (lex_at_phrase (s->lexer, stop1)
5561 || (stop2 && lex_at_phrase (s->lexer, stop2)))
5564 if (lex_at_phrase (s->lexer, "END MATRIX"))
5566 msg (SE, _("Premature END MATRIX within %s."), command_name);
5570 struct matrix_cmd *cmd = matrix_parse_command (s);
5574 if (c->n >= allocated)
5575 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
5576 c->commands[c->n++] = cmd;
5581 matrix_parse_do_if_clause (struct matrix_state *s,
5582 struct do_if_command *ifc,
5583 bool parse_condition,
5584 size_t *allocated_clauses)
5586 if (ifc->n_clauses >= *allocated_clauses)
5587 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5588 sizeof *ifc->clauses);
5589 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5590 *c = (struct do_if_clause) { .condition = NULL };
5592 if (parse_condition)
5594 c->condition = matrix_parse_expr (s);
5599 return matrix_parse_commands (s, &c->commands, "DO IF", "ELSE", "END IF");
5602 static struct matrix_cmd *
5603 matrix_parse_do_if (struct matrix_state *s)
5605 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5606 *cmd = (struct matrix_cmd) {
5608 .do_if = { .n_clauses = 0 }
5611 size_t allocated_clauses = 0;
5614 if (!matrix_parse_do_if_clause (s, &cmd->do_if, true, &allocated_clauses))
5617 while (lex_match_phrase (s->lexer, "ELSE IF"));
5619 if (lex_match_id (s->lexer, "ELSE")
5620 && !matrix_parse_do_if_clause (s, &cmd->do_if, false, &allocated_clauses))
5623 if (!lex_match_phrase (s->lexer, "END IF"))
5628 matrix_cmd_destroy (cmd);
5633 matrix_cmd_execute_do_if (struct do_if_command *cmd)
5635 for (size_t i = 0; i < cmd->n_clauses; i++)
5637 struct do_if_clause *c = &cmd->clauses[i];
5641 if (!matrix_expr_evaluate_scalar (c->condition,
5642 i ? "ELSE IF" : "DO IF",
5647 for (size_t j = 0; j < c->commands.n; j++)
5648 if (!matrix_cmd_execute (c->commands.commands[j]))
5655 static struct matrix_cmd *
5656 matrix_parse_loop (struct matrix_state *s)
5658 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5659 *cmd = (struct matrix_cmd) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5661 struct loop_command *loop = &cmd->loop;
5662 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5664 struct substring name = lex_tokss (s->lexer);
5665 loop->var = matrix_var_lookup (s, name);
5667 loop->var = matrix_var_create (s, name);
5672 loop->start = matrix_parse_expr (s);
5673 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5676 loop->end = matrix_parse_expr (s);
5680 if (lex_match (s->lexer, T_BY))
5682 loop->increment = matrix_parse_expr (s);
5683 if (!loop->increment)
5688 if (lex_match_id (s->lexer, "IF"))
5690 loop->top_condition = matrix_parse_expr (s);
5691 if (!loop->top_condition)
5695 bool was_in_loop = s->in_loop;
5697 bool ok = matrix_parse_commands (s, &loop->commands, "LOOP",
5699 s->in_loop = was_in_loop;
5703 if (!lex_match_phrase (s->lexer, "END LOOP"))
5706 if (lex_match_id (s->lexer, "IF"))
5708 loop->bottom_condition = matrix_parse_expr (s);
5709 if (!loop->bottom_condition)
5716 matrix_cmd_destroy (cmd);
5721 matrix_cmd_execute_loop (struct loop_command *cmd)
5725 long int increment = 1;
5728 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5729 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5731 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5735 if (increment > 0 ? value > end
5736 : increment < 0 ? value < end
5741 int mxloops = settings_get_mxloops ();
5742 for (int i = 0; i < mxloops; i++)
5747 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5749 gsl_matrix_free (cmd->var->value);
5750 cmd->var->value = NULL;
5752 if (!cmd->var->value)
5753 cmd->var->value = gsl_matrix_alloc (1, 1);
5755 gsl_matrix_set (cmd->var->value, 0, 0, value);
5758 if (cmd->top_condition)
5761 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5766 for (size_t j = 0; j < cmd->commands.n; j++)
5767 if (!matrix_cmd_execute (cmd->commands.commands[j]))
5770 if (cmd->bottom_condition)
5773 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5774 "END LOOP IF", &d) || d > 0)
5781 if (increment > 0 ? value > end : value < end)
5787 static struct matrix_cmd *
5788 matrix_parse_break (struct matrix_state *s)
5792 msg (SE, _("BREAK not inside LOOP."));
5796 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5797 *cmd = (struct matrix_cmd) { .type = MCMD_BREAK };
5801 static struct matrix_cmd *
5802 matrix_parse_display (struct matrix_state *s)
5804 while (lex_token (s->lexer) != T_ENDCMD)
5806 if (!lex_match_id (s->lexer, "DICTIONARY")
5807 && !lex_match_id (s->lexer, "STATUS"))
5809 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5814 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5815 *cmd = (struct matrix_cmd) { .type = MCMD_DISPLAY, .display = { s } };
5820 compare_matrix_var_pointers (const void *a_, const void *b_)
5822 const struct matrix_var *const *ap = a_;
5823 const struct matrix_var *const *bp = b_;
5824 const struct matrix_var *a = *ap;
5825 const struct matrix_var *b = *bp;
5826 return strcmp (a->name, b->name);
5830 matrix_cmd_execute_display (struct display_command *cmd)
5832 const struct matrix_state *s = cmd->state;
5834 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5835 struct pivot_dimension *attr_dimension
5836 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5837 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5838 N_("Rows"), N_("Columns"));
5839 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5841 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5844 const struct matrix_var *var;
5845 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5847 vars[n_vars++] = var;
5848 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5850 struct pivot_dimension *rows = pivot_dimension_create (
5851 table, PIVOT_AXIS_ROW, N_("Variable"));
5852 for (size_t i = 0; i < n_vars; i++)
5854 const struct matrix_var *var = vars[i];
5855 pivot_category_create_leaf (
5856 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
5858 size_t r = var->value->size1;
5859 size_t c = var->value->size2;
5860 double values[] = { r, c, r * c * sizeof (double) / 1024 };
5861 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
5862 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
5865 pivot_table_submit (table);
5868 static struct matrix_cmd *
5869 matrix_parse_release (struct matrix_state *s)
5871 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
5872 *cmd = (struct matrix_cmd) { .type = MCMD_RELEASE };
5874 size_t allocated_vars = 0;
5875 while (lex_token (s->lexer) == T_ID)
5877 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
5880 if (cmd->release.n_vars >= allocated_vars)
5881 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
5882 sizeof *cmd->release.vars);
5883 cmd->release.vars[cmd->release.n_vars++] = var;
5886 lex_error (s->lexer, _("Variable name expected."));
5889 if (!lex_match (s->lexer, T_COMMA))
5897 matrix_cmd_execute_release (struct release_command *cmd)
5899 for (size_t i = 0; i < cmd->n_vars; i++)
5901 struct matrix_var *var = cmd->vars[i];
5902 gsl_matrix_free (var->value);
5907 static struct save_file *
5908 save_file_create (struct matrix_state *s, struct file_handle *fh,
5909 struct string_array *variables,
5910 struct matrix_expr *names,
5911 struct stringi_set *strings)
5913 for (size_t i = 0; i < s->n_save_files; i++)
5915 struct save_file *sf = s->save_files[i];
5916 if (fh_equal (sf->file, fh))
5920 string_array_destroy (variables);
5921 matrix_expr_destroy (names);
5922 stringi_set_destroy (strings);
5928 struct save_file *sf = xmalloc (sizeof *sf);
5929 *sf = (struct save_file) {
5931 .dataset = s->dataset,
5932 .variables = *variables,
5934 .strings = STRINGI_SET_INITIALIZER (sf->strings),
5937 stringi_set_swap (&sf->strings, strings);
5938 stringi_set_destroy (strings);
5940 s->save_files = xrealloc (s->save_files,
5941 (s->n_save_files + 1) * sizeof *s->save_files);
5942 s->save_files[s->n_save_files++] = sf;
5946 static struct casewriter *
5947 save_file_open (struct save_file *sf, gsl_matrix *m)
5949 if (sf->writer || sf->error)
5953 size_t n_variables = caseproto_get_n_widths (
5954 casewriter_get_proto (sf->writer));
5955 if (m->size2 != n_variables)
5957 msg (ME, _("The first SAVE to %s within this matrix program "
5958 "had %zu columns, so a %zu×%zu matrix cannot be "
5960 sf->file == fh_inline_file () ? "*" : fh_get_name (sf->file),
5961 n_variables, m->size1, m->size2);
5969 struct dictionary *dict = dict_create (get_default_encoding ());
5971 /* Fill 'names' with user-specified names if there were any, either from
5972 sf->variables or sf->names. */
5973 struct string_array names = { .n = 0 };
5974 if (sf->variables.n)
5975 string_array_clone (&names, &sf->variables);
5978 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
5979 if (nm && is_vector (nm))
5981 gsl_vector nv = to_vector (nm);
5982 for (size_t i = 0; i < nv.size; i++)
5984 char *name = trimmed_string (gsl_vector_get (&nv, i));
5985 if (dict_id_is_valid (dict, name, true))
5986 string_array_append_nocopy (&names, name);
5991 gsl_matrix_free (nm);
5994 struct stringi_set strings;
5995 stringi_set_clone (&strings, &sf->strings);
5997 for (size_t i = 0; dict_get_var_cnt (dict) < m->size2; i++)
6002 name = names.strings[i];
6005 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
6009 int width = stringi_set_delete (&strings, name) ? 8 : 0;
6010 struct variable *var = dict_create_var (dict, name, width);
6013 msg (ME, _("Duplicate variable name %s in SAVE statement."),
6019 if (!stringi_set_is_empty (&strings))
6021 size_t count = stringi_set_count (&strings);
6022 const char *example = stringi_set_node_get_string (
6023 stringi_set_first (&strings));
6026 msg (ME, _("The SAVE command STRINGS subcommand specifies an unknown "
6027 "variable %s."), example);
6029 msg (ME, ngettext ("The SAVE command STRINGS subcommand specifies %zu "
6030 "unknown variable: %s.",
6031 "The SAVE command STRINGS subcommand specifies %zu "
6032 "unknown variables, including %s.",
6037 stringi_set_destroy (&strings);
6038 string_array_destroy (&names);
6047 if (sf->file == fh_inline_file ())
6048 sf->writer = autopaging_writer_create (dict_get_proto (dict));
6050 sf->writer = any_writer_open (sf->file, dict);
6062 save_file_destroy (struct save_file *sf)
6066 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
6068 dataset_set_dict (sf->dataset, sf->dict);
6069 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
6073 casewriter_destroy (sf->writer);
6074 dict_unref (sf->dict);
6076 fh_unref (sf->file);
6077 string_array_destroy (&sf->variables);
6078 matrix_expr_destroy (sf->names);
6079 stringi_set_destroy (&sf->strings);
6084 static struct matrix_cmd *
6085 matrix_parse_save (struct matrix_state *s)
6087 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6088 *cmd = (struct matrix_cmd) {
6090 .save = { .expression = NULL },
6093 struct file_handle *fh = NULL;
6094 struct save_command *save = &cmd->save;
6096 struct string_array variables = STRING_ARRAY_INITIALIZER;
6097 struct matrix_expr *names = NULL;
6098 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
6100 save->expression = matrix_parse_exp (s);
6101 if (!save->expression)
6104 while (lex_match (s->lexer, T_SLASH))
6106 if (lex_match_id (s->lexer, "OUTFILE"))
6108 lex_match (s->lexer, T_EQUALS);
6111 fh = (lex_match (s->lexer, T_ASTERISK)
6113 : fh_parse (s->lexer, FH_REF_FILE, s->session));
6117 else if (lex_match_id (s->lexer, "VARIABLES"))
6119 lex_match (s->lexer, T_EQUALS);
6123 struct dictionary *d = dict_create (get_default_encoding ());
6124 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
6125 PV_NO_SCRATCH | PV_NO_DUPLICATE);
6130 string_array_clear (&variables);
6131 variables = (struct string_array) {
6137 else if (lex_match_id (s->lexer, "NAMES"))
6139 lex_match (s->lexer, T_EQUALS);
6140 matrix_expr_destroy (names);
6141 names = matrix_parse_exp (s);
6145 else if (lex_match_id (s->lexer, "STRINGS"))
6147 lex_match (s->lexer, T_EQUALS);
6148 while (lex_token (s->lexer) == T_ID)
6150 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
6152 if (!lex_match (s->lexer, T_COMMA))
6158 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
6166 if (s->prev_save_file)
6167 fh = fh_ref (s->prev_save_file);
6170 lex_sbc_missing ("OUTFILE");
6174 fh_unref (s->prev_save_file);
6175 s->prev_save_file = fh_ref (fh);
6177 if (variables.n && names)
6179 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
6180 matrix_expr_destroy (names);
6184 save->sf = save_file_create (s, fh, &variables, names, &strings);
6188 string_array_destroy (&variables);
6189 matrix_expr_destroy (names);
6190 stringi_set_destroy (&strings);
6192 matrix_cmd_destroy (cmd);
6197 matrix_cmd_execute_save (const struct save_command *save)
6199 gsl_matrix *m = matrix_expr_evaluate (save->expression);
6203 struct casewriter *writer = save_file_open (save->sf, m);
6206 gsl_matrix_free (m);
6210 const struct caseproto *proto = casewriter_get_proto (writer);
6211 for (size_t y = 0; y < m->size1; y++)
6213 struct ccase *c = case_create (proto);
6214 for (size_t x = 0; x < m->size2; x++)
6216 double d = gsl_matrix_get (m, y, x);
6217 int width = caseproto_get_width (proto, x);
6218 union value *value = case_data_rw_idx (c, x);
6222 memcpy (value->s, &d, width);
6224 casewriter_write (writer, c);
6226 gsl_matrix_free (m);
6229 static struct matrix_cmd *
6230 matrix_parse_read (struct matrix_state *s)
6232 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6233 *cmd = (struct matrix_cmd) {
6235 .read = { .format = FMT_F },
6238 struct file_handle *fh = NULL;
6239 char *encoding = NULL;
6240 struct read_command *read = &cmd->read;
6241 read->dst = matrix_lvalue_parse (s);
6246 int repetitions = 0;
6247 int record_width = 0;
6248 bool seen_format = false;
6249 while (lex_match (s->lexer, T_SLASH))
6251 if (lex_match_id (s->lexer, "FILE"))
6253 lex_match (s->lexer, T_EQUALS);
6256 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6260 else if (lex_match_id (s->lexer, "ENCODING"))
6262 lex_match (s->lexer, T_EQUALS);
6263 if (!lex_force_string (s->lexer))
6267 encoding = ss_xstrdup (lex_tokss (s->lexer));
6271 else if (lex_match_id (s->lexer, "FIELD"))
6273 lex_match (s->lexer, T_EQUALS);
6275 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6277 read->c1 = lex_integer (s->lexer);
6279 if (!lex_force_match (s->lexer, T_TO)
6280 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
6282 read->c2 = lex_integer (s->lexer) + 1;
6285 record_width = read->c2 - read->c1;
6286 if (lex_match (s->lexer, T_BY))
6288 if (!lex_force_int_range (s->lexer, "BY", 1,
6289 read->c2 - read->c1))
6291 by = lex_integer (s->lexer);
6294 if (record_width % by)
6296 msg (SE, _("BY %d does not evenly divide record width %d."),
6304 else if (lex_match_id (s->lexer, "SIZE"))
6306 lex_match (s->lexer, T_EQUALS);
6307 matrix_expr_destroy (read->size);
6308 read->size = matrix_parse_exp (s);
6312 else if (lex_match_id (s->lexer, "MODE"))
6314 lex_match (s->lexer, T_EQUALS);
6315 if (lex_match_id (s->lexer, "RECTANGULAR"))
6316 read->symmetric = false;
6317 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6318 read->symmetric = true;
6321 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6325 else if (lex_match_id (s->lexer, "REREAD"))
6326 read->reread = true;
6327 else if (lex_match_id (s->lexer, "FORMAT"))
6331 lex_sbc_only_once ("FORMAT");
6336 lex_match (s->lexer, T_EQUALS);
6338 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6341 const char *p = lex_tokcstr (s->lexer);
6342 if (c_isdigit (p[0]))
6344 repetitions = atoi (p);
6345 p += strspn (p, "0123456789");
6346 if (!fmt_from_name (p, &read->format))
6348 lex_error (s->lexer, _("Unknown format %s."), p);
6353 else if (fmt_from_name (p, &read->format))
6357 struct fmt_spec format;
6358 if (!parse_format_specifier (s->lexer, &format))
6360 read->format = format.type;
6366 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6367 "REREAD", "FORMAT");
6374 lex_sbc_missing ("FIELD");
6378 if (!read->dst->n_indexes && !read->size)
6380 msg (SE, _("SIZE is required for reading data into a full matrix "
6381 "(as opposed to a submatrix)."));
6387 if (s->prev_read_file)
6388 fh = fh_ref (s->prev_read_file);
6391 lex_sbc_missing ("FILE");
6395 fh_unref (s->prev_read_file);
6396 s->prev_read_file = fh_ref (fh);
6398 read->rf = read_file_create (s, fh);
6402 free (read->rf->encoding);
6403 read->rf->encoding = encoding;
6407 /* Field width may be specified in multiple ways:
6410 2. The format on FORMAT.
6411 3. The repetition factor on FORMAT.
6413 (2) and (3) are mutually exclusive.
6415 If more than one of these is present, they must agree. If none of them is
6416 present, then free-field format is used.
6418 if (repetitions > record_width)
6420 msg (SE, _("%d repetitions cannot fit in record width %d."),
6421 repetitions, record_width);
6424 int w = (repetitions ? record_width / repetitions
6430 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6431 "which implies field width %d, "
6432 "but BY specifies field width %d."),
6433 repetitions, record_width, w, by);
6435 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6444 matrix_cmd_destroy (cmd);
6450 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6451 struct substring data, size_t y, size_t x,
6452 int first_column, int last_column, char *error)
6454 int line_number = dfm_get_line_number (reader);
6455 struct msg_location *location = xmalloc (sizeof *location);
6456 *location = (struct msg_location) {
6457 .file_name = intern_new (dfm_get_file_name (reader)),
6458 .p[0] = { .line = line_number, .column = first_column },
6459 .p[1] = { .line = line_number, .column = last_column },
6461 struct msg *m = xmalloc (sizeof *m);
6463 .category = MSG_C_DATA,
6464 .severity = MSG_S_WARNING,
6465 .location = location,
6466 .text = xasprintf (_("Error reading \"%.*s\" as format %s "
6467 "for matrix row %zu, column %zu: %s"),
6468 (int) data.length, data.string, fmt_name (format),
6469 y + 1, x + 1, error),
6477 matrix_read_set_field (struct read_command *read, struct dfm_reader *reader,
6478 gsl_matrix *m, struct substring p, size_t y, size_t x,
6479 const char *line_start)
6481 const char *input_encoding = dfm_reader_get_encoding (reader);
6484 if (fmt_is_numeric (read->format))
6487 error = data_in (p, input_encoding, read->format,
6488 settings_get_fmt_settings (), &v, 0, NULL);
6489 if (!error && v.f == SYSMIS)
6490 error = xstrdup (_("Matrix data may not contain missing value."));
6495 uint8_t s[sizeof (double)];
6496 union value v = { .s = s };
6497 error = data_in (p, input_encoding, read->format,
6498 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6499 memcpy (&f, s, sizeof f);
6504 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6505 int c2 = c1 + ss_utf8_count_columns (p);
6506 parse_error (reader, read->format, p, y, x, c1, c2, error);
6510 gsl_matrix_set (m, y, x, f);
6511 if (read->symmetric && x != y)
6512 gsl_matrix_set (m, x, y, f);
6517 matrix_read_line (struct read_command *read, struct dfm_reader *reader,
6518 struct substring *line, const char **startp)
6520 if (dfm_eof (reader))
6522 msg (SE, _("Unexpected end of file reading matrix data."));
6525 dfm_expand_tabs (reader);
6526 struct substring record = dfm_get_record (reader);
6527 /* XXX need to recode record into UTF-8 */
6528 *startp = record.string;
6529 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6534 matrix_read (struct read_command *read, struct dfm_reader *reader,
6537 for (size_t y = 0; y < m->size1; y++)
6539 size_t nx = read->symmetric ? y + 1 : m->size2;
6541 struct substring line = ss_empty ();
6542 const char *line_start = line.string;
6543 for (size_t x = 0; x < nx; x++)
6550 ss_ltrim (&line, ss_cstr (" ,"));
6551 if (!ss_is_empty (line))
6553 if (!matrix_read_line (read, reader, &line, &line_start))
6555 dfm_forward_record (reader);
6558 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6562 if (!matrix_read_line (read, reader, &line, &line_start))
6564 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6565 int f = x % fields_per_line;
6566 if (f == fields_per_line - 1)
6567 dfm_forward_record (reader);
6569 p = ss_substr (line, read->w * f, read->w);
6572 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6576 dfm_forward_record (reader);
6579 ss_ltrim (&line, ss_cstr (" ,"));
6580 if (!ss_is_empty (line))
6583 msg (SW, _("Trailing garbage on line \"%.*s\""),
6584 (int) line.length, line.string);
6591 matrix_cmd_execute_read (struct read_command *read)
6593 struct index_vector iv0, iv1;
6594 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6597 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6600 gsl_matrix *m = matrix_expr_evaluate (read->size);
6606 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6607 "not a %zu×%zu matrix."), m->size1, m->size2);
6608 gsl_matrix_free (m);
6614 gsl_vector v = to_vector (m);
6618 d[0] = gsl_vector_get (&v, 0);
6621 else if (v.size == 2)
6623 d[0] = gsl_vector_get (&v, 0);
6624 d[1] = gsl_vector_get (&v, 1);
6628 msg (SE, _("SIZE must evaluate to a scalar or a 2-element vector, "
6629 "not a %zu×%zu matrix."),
6630 m->size1, m->size2),
6631 gsl_matrix_free (m);
6636 gsl_matrix_free (m);
6638 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6640 msg (SE, _("Matrix dimensions %g×%g specified on SIZE "
6641 "are outside valid range."),
6652 if (read->dst->n_indexes)
6654 size_t submatrix_size[2];
6655 if (read->dst->n_indexes == 2)
6657 submatrix_size[0] = iv0.n;
6658 submatrix_size[1] = iv1.n;
6660 else if (read->dst->var->value->size1 == 1)
6662 submatrix_size[0] = 1;
6663 submatrix_size[1] = iv0.n;
6667 submatrix_size[0] = iv0.n;
6668 submatrix_size[1] = 1;
6673 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6675 msg (SE, _("Matrix dimensions %zu×%zu specified on SIZE "
6676 "differ from submatrix dimensions %zu×%zu."),
6678 submatrix_size[0], submatrix_size[1]);
6686 size[0] = submatrix_size[0];
6687 size[1] = submatrix_size[1];
6691 struct dfm_reader *reader = read_file_open (read->rf);
6693 dfm_reread_record (reader, 1);
6695 if (read->symmetric && size[0] != size[1])
6697 msg (SE, _("Cannot read non-square %zu×%zu matrix "
6698 "using READ with MODE=SYMMETRIC."),
6704 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6705 matrix_read (read, reader, tmp);
6706 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp);
6709 static struct matrix_cmd *
6710 matrix_parse_write (struct matrix_state *s)
6712 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
6713 *cmd = (struct matrix_cmd) {
6717 struct file_handle *fh = NULL;
6718 char *encoding = NULL;
6719 struct write_command *write = &cmd->write;
6720 write->expression = matrix_parse_exp (s);
6721 if (!write->expression)
6725 int repetitions = 0;
6726 int record_width = 0;
6727 enum fmt_type format = FMT_F;
6728 bool has_format = false;
6729 while (lex_match (s->lexer, T_SLASH))
6731 if (lex_match_id (s->lexer, "OUTFILE"))
6733 lex_match (s->lexer, T_EQUALS);
6736 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6740 else if (lex_match_id (s->lexer, "ENCODING"))
6742 lex_match (s->lexer, T_EQUALS);
6743 if (!lex_force_string (s->lexer))
6747 encoding = ss_xstrdup (lex_tokss (s->lexer));
6751 else if (lex_match_id (s->lexer, "FIELD"))
6753 lex_match (s->lexer, T_EQUALS);
6755 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6757 write->c1 = lex_integer (s->lexer);
6759 if (!lex_force_match (s->lexer, T_TO)
6760 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
6762 write->c2 = lex_integer (s->lexer) + 1;
6765 record_width = write->c2 - write->c1;
6766 if (lex_match (s->lexer, T_BY))
6768 if (!lex_force_int_range (s->lexer, "BY", 1,
6769 write->c2 - write->c1))
6771 by = lex_integer (s->lexer);
6774 if (record_width % by)
6776 msg (SE, _("BY %d does not evenly divide record width %d."),
6784 else if (lex_match_id (s->lexer, "MODE"))
6786 lex_match (s->lexer, T_EQUALS);
6787 if (lex_match_id (s->lexer, "RECTANGULAR"))
6788 write->triangular = false;
6789 else if (lex_match_id (s->lexer, "TRIANGULAR"))
6790 write->triangular = true;
6793 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
6797 else if (lex_match_id (s->lexer, "HOLD"))
6799 else if (lex_match_id (s->lexer, "FORMAT"))
6801 if (has_format || write->format)
6803 lex_sbc_only_once ("FORMAT");
6807 lex_match (s->lexer, T_EQUALS);
6809 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6812 const char *p = lex_tokcstr (s->lexer);
6813 if (c_isdigit (p[0]))
6815 repetitions = atoi (p);
6816 p += strspn (p, "0123456789");
6817 if (!fmt_from_name (p, &format))
6819 lex_error (s->lexer, _("Unknown format %s."), p);
6825 else if (fmt_from_name (p, &format))
6832 struct fmt_spec spec;
6833 if (!parse_format_specifier (s->lexer, &spec))
6835 write->format = xmemdup (&spec, sizeof spec);
6840 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
6848 lex_sbc_missing ("FIELD");
6854 if (s->prev_write_file)
6855 fh = fh_ref (s->prev_write_file);
6858 lex_sbc_missing ("OUTFILE");
6862 fh_unref (s->prev_write_file);
6863 s->prev_write_file = fh_ref (fh);
6865 write->wf = write_file_create (s, fh);
6869 free (write->wf->encoding);
6870 write->wf->encoding = encoding;
6874 /* Field width may be specified in multiple ways:
6877 2. The format on FORMAT.
6878 3. The repetition factor on FORMAT.
6880 (2) and (3) are mutually exclusive.
6882 If more than one of these is present, they must agree. If none of them is
6883 present, then free-field format is used.
6885 if (repetitions > record_width)
6887 msg (SE, _("%d repetitions cannot fit in record width %d."),
6888 repetitions, record_width);
6891 int w = (repetitions ? record_width / repetitions
6892 : write->format ? write->format->w
6897 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6898 "which implies field width %d, "
6899 "but BY specifies field width %d."),
6900 repetitions, record_width, w, by);
6902 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6906 if (w && !write->format)
6908 write->format = xmalloc (sizeof *write->format);
6909 *write->format = (struct fmt_spec) { .type = format, .w = w };
6911 if (!fmt_check_output (write->format))
6915 if (write->format && fmt_var_width (write->format) > sizeof (double))
6917 char s[FMT_STRING_LEN_MAX + 1];
6918 fmt_to_string (write->format, s);
6919 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
6920 s, sizeof (double));
6928 matrix_cmd_destroy (cmd);
6933 matrix_cmd_execute_write (struct write_command *write)
6935 gsl_matrix *m = matrix_expr_evaluate (write->expression);
6939 if (write->triangular && m->size1 != m->size2)
6941 msg (SE, _("WRITE with MODE=TRIANGULAR requires a square matrix but "
6942 "the matrix to be written has dimensions %zu×%zu."),
6943 m->size1, m->size2);
6944 gsl_matrix_free (m);
6948 struct dfm_writer *writer = write_file_open (write->wf);
6949 if (!writer || !m->size1)
6951 gsl_matrix_free (m);
6955 const struct fmt_settings *settings = settings_get_fmt_settings ();
6956 struct u8_line *line = write->wf->held;
6957 for (size_t y = 0; y < m->size1; y++)
6961 line = xmalloc (sizeof *line);
6962 u8_line_init (line);
6964 size_t nx = write->triangular ? y + 1 : m->size2;
6966 for (size_t x = 0; x < nx; x++)
6969 double f = gsl_matrix_get (m, y, x);
6973 if (fmt_is_numeric (write->format->type))
6976 v.s = (uint8_t *) &f;
6977 s = data_out (&v, NULL, write->format, settings);
6981 s = xmalloc (DBL_BUFSIZE_BOUND);
6982 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
6983 >= DBL_BUFSIZE_BOUND)
6986 size_t len = strlen (s);
6987 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
6988 if (width + x0 > write->c2)
6990 dfm_put_record_utf8 (writer, line->s.ss.string,
6992 u8_line_clear (line);
6995 u8_line_put (line, x0, x0 + width, s, len);
6998 x0 += write->format ? write->format->w : width + 1;
7001 if (y + 1 >= m->size1 && write->hold)
7003 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
7004 u8_line_clear (line);
7008 u8_line_destroy (line);
7012 write->wf->held = line;
7014 gsl_matrix_free (m);
7017 static struct matrix_cmd *
7018 matrix_parse_get (struct matrix_state *s)
7020 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7021 *cmd = (struct matrix_cmd) {
7024 .dataset = s->dataset,
7025 .user = { .treatment = MGET_ERROR },
7026 .system = { .treatment = MGET_ERROR },
7030 struct get_command *get = &cmd->get;
7031 get->dst = matrix_lvalue_parse (s);
7035 while (lex_match (s->lexer, T_SLASH))
7037 if (lex_match_id (s->lexer, "FILE"))
7039 lex_match (s->lexer, T_EQUALS);
7041 fh_unref (get->file);
7042 if (lex_match (s->lexer, T_ASTERISK))
7046 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7051 else if (lex_match_id (s->lexer, "ENCODING"))
7053 lex_match (s->lexer, T_EQUALS);
7054 if (!lex_force_string (s->lexer))
7057 free (get->encoding);
7058 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
7062 else if (lex_match_id (s->lexer, "VARIABLES"))
7064 lex_match (s->lexer, T_EQUALS);
7068 lex_sbc_only_once ("VARIABLES");
7072 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
7075 else if (lex_match_id (s->lexer, "NAMES"))
7077 lex_match (s->lexer, T_EQUALS);
7078 if (!lex_force_id (s->lexer))
7081 struct substring name = lex_tokss (s->lexer);
7082 get->names = matrix_var_lookup (s, name);
7084 get->names = matrix_var_create (s, name);
7087 else if (lex_match_id (s->lexer, "MISSING"))
7089 lex_match (s->lexer, T_EQUALS);
7090 if (lex_match_id (s->lexer, "ACCEPT"))
7091 get->user.treatment = MGET_ACCEPT;
7092 else if (lex_match_id (s->lexer, "OMIT"))
7093 get->user.treatment = MGET_OMIT;
7094 else if (lex_is_number (s->lexer))
7096 get->user.treatment = MGET_RECODE;
7097 get->user.substitute = lex_number (s->lexer);
7102 lex_error (s->lexer, NULL);
7106 else if (lex_match_id (s->lexer, "SYSMIS"))
7108 lex_match (s->lexer, T_EQUALS);
7109 if (lex_match_id (s->lexer, "OMIT"))
7110 get->system.treatment = MGET_OMIT;
7111 else if (lex_is_number (s->lexer))
7113 get->system.treatment = MGET_RECODE;
7114 get->system.substitute = lex_number (s->lexer);
7119 lex_error (s->lexer, NULL);
7125 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
7126 "MISSING", "SYSMIS");
7131 if (get->user.treatment != MGET_ACCEPT)
7132 get->system.treatment = MGET_ERROR;
7137 matrix_cmd_destroy (cmd);
7142 matrix_cmd_execute_get__ (struct get_command *get, struct casereader *reader,
7143 const struct dictionary *dict)
7145 struct variable **vars;
7150 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
7151 &vars, &n_vars, PV_NUMERIC))
7156 n_vars = dict_get_var_cnt (dict);
7157 vars = xnmalloc (n_vars, sizeof *vars);
7158 for (size_t i = 0; i < n_vars; i++)
7160 struct variable *var = dict_get_var (dict, i);
7161 if (!var_is_numeric (var))
7163 msg (SE, _("GET: Variable %s is not numeric."),
7164 var_get_name (var));
7174 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
7175 for (size_t i = 0; i < n_vars; i++)
7177 char s[sizeof (double)];
7179 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
7180 memcpy (&f, s, sizeof f);
7181 gsl_matrix_set (names, i, 0, f);
7184 gsl_matrix_free (get->names->value);
7185 get->names->value = names;
7189 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
7190 long long int casenum = 1;
7192 for (struct ccase *c = casereader_read (reader); c;
7193 c = casereader_read (reader), casenum++)
7195 if (n_rows >= m->size1)
7197 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
7198 for (size_t y = 0; y < n_rows; y++)
7199 for (size_t x = 0; x < n_vars; x++)
7200 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
7201 gsl_matrix_free (m);
7206 for (size_t x = 0; x < n_vars; x++)
7208 const struct variable *var = vars[x];
7209 double d = case_num (c, var);
7212 if (get->system.treatment == MGET_RECODE)
7213 d = get->system.substitute;
7214 else if (get->system.treatment == MGET_OMIT)
7218 msg (SE, _("GET: Variable %s in case %lld "
7219 "is system-missing."),
7220 var_get_name (var), casenum);
7224 else if (var_is_num_missing (var, d, MV_USER))
7226 if (get->user.treatment == MGET_RECODE)
7227 d = get->user.substitute;
7228 else if (get->user.treatment == MGET_OMIT)
7230 else if (get->user.treatment != MGET_ACCEPT)
7232 msg (SE, _("GET: Variable %s in case %lld has user-missing "
7234 var_get_name (var), casenum, d);
7238 gsl_matrix_set (m, n_rows, x, d);
7249 matrix_lvalue_evaluate_and_assign (get->dst, m);
7252 gsl_matrix_free (m);
7257 matrix_open_casereader (const char *command_name,
7258 struct file_handle *file, const char *encoding,
7259 struct dataset *dataset,
7260 struct casereader **readerp, struct dictionary **dictp)
7264 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
7265 return *readerp != NULL;
7269 if (dict_get_var_cnt (dataset_dict (dataset)) == 0)
7271 msg (ME, _("The %s command cannot read an empty active file."),
7275 *readerp = proc_open (dataset);
7276 *dictp = dict_ref (dataset_dict (dataset));
7282 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
7283 struct casereader *reader, struct dictionary *dict)
7286 casereader_destroy (reader);
7288 proc_commit (dataset);
7292 matrix_cmd_execute_get (struct get_command *get)
7294 struct casereader *r;
7295 struct dictionary *d;
7296 if (matrix_open_casereader ("GET", get->file, get->encoding, get->dataset,
7299 matrix_cmd_execute_get__ (get, r, d);
7300 matrix_close_casereader (get->file, get->dataset, r, d);
7305 match_rowtype (struct lexer *lexer)
7307 static const char *rowtypes[] = {
7308 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7310 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7312 for (size_t i = 0; i < n_rowtypes; i++)
7313 if (lex_match_id (lexer, rowtypes[i]))
7316 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7321 parse_var_names (struct lexer *lexer, struct string_array *sa)
7323 lex_match (lexer, T_EQUALS);
7325 string_array_clear (sa);
7327 struct dictionary *dict = dict_create (get_default_encoding ());
7330 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7331 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7336 for (size_t i = 0; i < n_names; i++)
7337 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7338 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7340 msg (SE, _("Variable name %s is reserved."), names[i]);
7341 for (size_t j = 0; j < n_names; j++)
7347 string_array_clear (sa);
7348 sa->strings = names;
7349 sa->n = sa->allocated = n_names;
7355 msave_common_destroy (struct msave_common *common)
7359 fh_unref (common->outfile);
7360 string_array_destroy (&common->variables);
7361 string_array_destroy (&common->fnames);
7362 string_array_destroy (&common->snames);
7364 for (size_t i = 0; i < common->n_factors; i++)
7365 matrix_expr_destroy (common->factors[i]);
7366 free (common->factors);
7368 for (size_t i = 0; i < common->n_splits; i++)
7369 matrix_expr_destroy (common->splits[i]);
7370 free (common->splits);
7372 dict_unref (common->dict);
7373 casewriter_destroy (common->writer);
7380 compare_variables (const char *keyword,
7381 const struct string_array *new,
7382 const struct string_array *old)
7388 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7389 "on the first MSAVE within MATRIX."), keyword);
7392 else if (!string_array_equal_case (old, new))
7394 msg (SE, _("%s must specify the same variables each time within "
7395 "a given MATRIX."), keyword);
7402 static struct matrix_cmd *
7403 matrix_parse_msave (struct matrix_state *s)
7405 struct msave_common *common = xmalloc (sizeof *common);
7406 *common = (struct msave_common) { .outfile = NULL };
7408 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7409 *cmd = (struct matrix_cmd) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7411 struct matrix_expr *splits = NULL;
7412 struct matrix_expr *factors = NULL;
7414 struct msave_command *msave = &cmd->msave;
7415 msave->expr = matrix_parse_exp (s);
7419 while (lex_match (s->lexer, T_SLASH))
7421 if (lex_match_id (s->lexer, "TYPE"))
7423 lex_match (s->lexer, T_EQUALS);
7425 msave->rowtype = match_rowtype (s->lexer);
7426 if (!msave->rowtype)
7429 else if (lex_match_id (s->lexer, "OUTFILE"))
7431 lex_match (s->lexer, T_EQUALS);
7433 fh_unref (common->outfile);
7434 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7435 if (!common->outfile)
7438 else if (lex_match_id (s->lexer, "VARIABLES"))
7440 if (!parse_var_names (s->lexer, &common->variables))
7443 else if (lex_match_id (s->lexer, "FNAMES"))
7445 if (!parse_var_names (s->lexer, &common->fnames))
7448 else if (lex_match_id (s->lexer, "SNAMES"))
7450 if (!parse_var_names (s->lexer, &common->snames))
7453 else if (lex_match_id (s->lexer, "SPLIT"))
7455 lex_match (s->lexer, T_EQUALS);
7457 matrix_expr_destroy (splits);
7458 splits = matrix_parse_exp (s);
7462 else if (lex_match_id (s->lexer, "FACTOR"))
7464 lex_match (s->lexer, T_EQUALS);
7466 matrix_expr_destroy (factors);
7467 factors = matrix_parse_exp (s);
7473 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7474 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7478 if (!msave->rowtype)
7480 lex_sbc_missing ("TYPE");
7486 if (common->fnames.n && !factors)
7488 msg (SE, _("FNAMES requires FACTOR."));
7491 if (common->snames.n && !splits)
7493 msg (SE, _("SNAMES requires SPLIT."));
7496 if (!common->outfile)
7498 lex_sbc_missing ("OUTFILE");
7505 if (common->outfile)
7507 if (!fh_equal (common->outfile, s->common->outfile))
7509 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7510 "within a single MATRIX command."));
7514 if (!compare_variables ("VARIABLES",
7515 &common->variables, &s->common->variables)
7516 || !compare_variables ("FNAMES", &common->fnames, &s->common->fnames)
7517 || !compare_variables ("SNAMES", &common->snames, &s->common->snames))
7519 msave_common_destroy (common);
7521 msave->common = s->common;
7523 struct msave_common *c = s->common;
7526 if (c->n_factors >= c->allocated_factors)
7527 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7528 sizeof *c->factors);
7529 c->factors[c->n_factors++] = factors;
7531 if (c->n_factors > 0)
7532 msave->factors = c->factors[c->n_factors - 1];
7536 if (c->n_splits >= c->allocated_splits)
7537 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7539 c->splits[c->n_splits++] = splits;
7541 if (c->n_splits > 0)
7542 msave->splits = c->splits[c->n_splits - 1];
7547 matrix_expr_destroy (splits);
7548 matrix_expr_destroy (factors);
7549 msave_common_destroy (common);
7550 matrix_cmd_destroy (cmd);
7555 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7557 gsl_matrix *m = matrix_expr_evaluate (e);
7563 msg (SE, _("%s expression must evaluate to vector, "
7564 "not a %zu×%zu matrix."),
7565 name, m->size1, m->size2);
7566 gsl_matrix_free (m);
7570 return matrix_to_vector (m);
7574 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7576 for (size_t i = 0; i < vars->n; i++)
7577 if (!dict_create_var (d, vars->strings[i], 0))
7578 return vars->strings[i];
7582 static struct dictionary *
7583 msave_create_dict (const struct msave_common *common)
7585 struct dictionary *dict = dict_create (get_default_encoding ());
7587 const char *dup_split = msave_add_vars (dict, &common->snames);
7590 /* Should not be possible because the parser ensures that the names are
7595 dict_create_var_assert (dict, "ROWTYPE_", 8);
7597 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7600 msg (SE, _("Duplicate or invalid FACTOR variable name %s."), dup_factor);
7604 dict_create_var_assert (dict, "VARNAME_", 8);
7606 const char *dup_var = msave_add_vars (dict, &common->variables);
7609 msg (SE, _("Duplicate or invalid variable name %s."), dup_var);
7621 matrix_cmd_execute_msave (struct msave_command *msave)
7623 struct msave_common *common = msave->common;
7624 gsl_matrix *m = NULL;
7625 gsl_vector *factors = NULL;
7626 gsl_vector *splits = NULL;
7628 m = matrix_expr_evaluate (msave->expr);
7632 if (!common->variables.n)
7633 for (size_t i = 0; i < m->size2; i++)
7634 string_array_append_nocopy (&common->variables,
7635 xasprintf ("COL%zu", i + 1));
7636 else if (m->size2 != common->variables.n)
7639 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7640 m->size2, common->variables.n);
7646 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7650 if (!common->fnames.n)
7651 for (size_t i = 0; i < factors->size; i++)
7652 string_array_append_nocopy (&common->fnames,
7653 xasprintf ("FAC%zu", i + 1));
7654 else if (factors->size != common->fnames.n)
7656 msg (SE, _("There are %zu factor variables, "
7657 "but %zu factor values were supplied."),
7658 common->fnames.n, factors->size);
7665 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7669 if (!common->snames.n)
7670 for (size_t i = 0; i < splits->size; i++)
7671 string_array_append_nocopy (&common->snames,
7672 xasprintf ("SPL%zu", i + 1));
7673 else if (splits->size != common->snames.n)
7675 msg (SE, _("There are %zu split variables, "
7676 "but %zu split values were supplied."),
7677 common->snames.n, splits->size);
7682 if (!common->writer)
7684 struct dictionary *dict = msave_create_dict (common);
7688 common->writer = any_writer_open (common->outfile, dict);
7689 if (!common->writer)
7695 common->dict = dict;
7698 bool matrix = (!strcmp (msave->rowtype, "COV")
7699 || !strcmp (msave->rowtype, "CORR"));
7700 for (size_t y = 0; y < m->size1; y++)
7702 struct ccase *c = case_create (dict_get_proto (common->dict));
7705 /* Split variables */
7707 for (size_t i = 0; i < splits->size; i++)
7708 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
7711 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7712 msave->rowtype, ' ');
7716 for (size_t i = 0; i < factors->size; i++)
7717 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
7720 const char *varname_ = (matrix && y < common->variables.n
7721 ? common->variables.strings[y]
7723 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
7726 /* Continuous variables. */
7727 for (size_t x = 0; x < m->size2; x++)
7728 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
7729 casewriter_write (common->writer, c);
7733 gsl_matrix_free (m);
7734 gsl_vector_free (factors);
7735 gsl_vector_free (splits);
7738 static struct matrix_cmd *
7739 matrix_parse_mget (struct matrix_state *s)
7741 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
7742 *cmd = (struct matrix_cmd) {
7746 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
7750 struct mget_command *mget = &cmd->mget;
7752 lex_match (s->lexer, T_SLASH);
7753 while (lex_token (s->lexer) != T_ENDCMD)
7755 if (lex_match_id (s->lexer, "FILE"))
7757 lex_match (s->lexer, T_EQUALS);
7759 fh_unref (mget->file);
7760 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7764 else if (lex_match_id (s->lexer, "ENCODING"))
7766 lex_match (s->lexer, T_EQUALS);
7767 if (!lex_force_string (s->lexer))
7770 free (mget->encoding);
7771 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
7775 else if (lex_match_id (s->lexer, "TYPE"))
7777 lex_match (s->lexer, T_EQUALS);
7778 while (lex_token (s->lexer) != T_SLASH
7779 && lex_token (s->lexer) != T_ENDCMD)
7781 const char *rowtype = match_rowtype (s->lexer);
7785 stringi_set_insert (&mget->rowtypes, rowtype);
7790 lex_error_expecting (s->lexer, "FILE", "TYPE");
7793 lex_match (s->lexer, T_SLASH);
7798 matrix_cmd_destroy (cmd);
7802 static const struct variable *
7803 get_a8_var (const struct dictionary *d, const char *name)
7805 const struct variable *v = dict_lookup_var (d, name);
7808 msg (SE, _("Matrix data file lacks %s variable."), name);
7811 if (var_get_width (v) != 8)
7813 msg (SE, _("%s variable in matrix data file must be "
7814 "8-byte string, but it has width %d."),
7815 name, var_get_width (v));
7822 var_changed (const struct ccase *ca, const struct ccase *cb,
7823 const struct variable *var)
7826 ? !value_equal (case_data (ca, var), case_data (cb, var),
7827 var_get_width (var))
7832 vars_changed (const struct ccase *ca, const struct ccase *cb,
7833 const struct dictionary *d,
7834 size_t first_var, size_t n_vars)
7836 for (size_t i = 0; i < n_vars; i++)
7838 const struct variable *v = dict_get_var (d, first_var + i);
7839 if (var_changed (ca, cb, v))
7846 vars_all_missing (const struct ccase *c, const struct dictionary *d,
7847 size_t first_var, size_t n_vars)
7849 for (size_t i = 0; i < n_vars; i++)
7850 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
7856 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
7857 const struct dictionary *d,
7858 const struct variable *rowtype_var,
7859 const struct stringi_set *accepted_rowtypes,
7860 struct matrix_state *s,
7861 size_t ss, size_t sn, size_t si,
7862 size_t fs, size_t fn, size_t fi,
7863 size_t cs, size_t cn,
7864 struct pivot_table *pt,
7865 struct pivot_dimension *var_dimension)
7870 /* Is this a matrix for pooled data, either where there are no factor
7871 variables or the factor variables are missing? */
7872 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
7874 struct substring rowtype = case_ss (rows[0], rowtype_var);
7875 ss_rtrim (&rowtype, ss_cstr (" "));
7876 if (!stringi_set_is_empty (accepted_rowtypes)
7877 && !stringi_set_contains_len (accepted_rowtypes,
7878 rowtype.string, rowtype.length))
7881 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
7882 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
7883 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
7884 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
7885 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
7886 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
7890 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
7891 (int) rowtype.length, rowtype.string);
7895 struct string name = DS_EMPTY_INITIALIZER;
7896 ds_put_cstr (&name, prefix);
7898 ds_put_format (&name, "F%zu", fi);
7900 ds_put_format (&name, "S%zu", si);
7902 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
7904 mv = matrix_var_create (s, ds_ss (&name));
7907 msg (SW, _("Matrix data file contains variable with existing name %s."),
7909 goto exit_free_name;
7912 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
7913 size_t n_missing = 0;
7914 for (size_t y = 0; y < n_rows; y++)
7916 for (size_t x = 0; x < cn; x++)
7918 struct variable *var = dict_get_var (d, cs + x);
7919 double value = case_num (rows[y], var);
7920 if (var_is_num_missing (var, value, MV_ANY))
7925 gsl_matrix_set (m, y, x, value);
7929 int var_index = pivot_category_create_leaf (
7930 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
7931 double values[] = { n_rows, cn };
7932 for (size_t j = 0; j < sn; j++)
7934 struct variable *var = dict_get_var (d, ss + j);
7935 const union value *value = case_data (rows[0], var);
7936 pivot_table_put2 (pt, j, var_index,
7937 pivot_value_new_var_value (var, value));
7939 for (size_t j = 0; j < fn; j++)
7941 struct variable *var = dict_get_var (d, fs + j);
7942 const union value sysmis = { .f = SYSMIS };
7943 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
7944 pivot_table_put2 (pt, j + sn, var_index,
7945 pivot_value_new_var_value (var, value));
7947 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
7948 pivot_table_put2 (pt, j + sn + fn, var_index,
7949 pivot_value_new_integer (values[j]));
7952 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
7953 "value, which was treated as zero.",
7954 "Matrix data file variable %s contains %zu missing "
7955 "values, which were treated as zero.", n_missing),
7956 ds_cstr (&name), n_missing);
7963 for (size_t y = 0; y < n_rows; y++)
7964 case_unref (rows[y]);
7968 matrix_cmd_execute_mget__ (struct mget_command *mget,
7969 struct casereader *r, const struct dictionary *d)
7971 const struct variable *rowtype_ = get_a8_var (d, "ROWTYPE_");
7972 const struct variable *varname_ = get_a8_var (d, "VARNAME_");
7973 if (!rowtype_ || !varname_)
7976 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
7978 msg (SE, _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
7981 if (var_get_dict_index (varname_) + 1 >= dict_get_var_cnt (d))
7983 msg (SE, _("Matrix data file contains no continuous variables."));
7987 for (size_t i = 0; i < dict_get_var_cnt (d); i++)
7989 const struct variable *v = dict_get_var (d, i);
7990 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
7993 _("Matrix data file contains unexpected string variable %s."),
7999 /* SPLIT variables. */
8001 size_t sn = var_get_dict_index (rowtype_);
8002 struct ccase *sc = NULL;
8005 /* FACTOR variables. */
8006 size_t fs = var_get_dict_index (rowtype_) + 1;
8007 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
8008 struct ccase *fc = NULL;
8011 /* Continuous variables. */
8012 size_t cs = var_get_dict_index (varname_) + 1;
8013 size_t cn = dict_get_var_cnt (d) - cs;
8014 struct ccase *cc = NULL;
8017 struct pivot_table *pt = pivot_table_create (
8018 N_("Matrix Variables Created by MGET"));
8019 struct pivot_dimension *attr_dimension = pivot_dimension_create (
8020 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
8021 struct pivot_dimension *var_dimension = pivot_dimension_create (
8022 pt, PIVOT_AXIS_ROW, N_("Variable"));
8025 struct pivot_category *splits = pivot_category_create_group (
8026 attr_dimension->root, N_("Split Values"));
8027 for (size_t i = 0; i < sn; i++)
8028 pivot_category_create_leaf (splits, pivot_value_new_variable (
8029 dict_get_var (d, ss + i)));
8033 struct pivot_category *factors = pivot_category_create_group (
8034 attr_dimension->root, N_("Factors"));
8035 for (size_t i = 0; i < fn; i++)
8036 pivot_category_create_leaf (factors, pivot_value_new_variable (
8037 dict_get_var (d, fs + i)));
8039 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
8040 N_("Rows"), N_("Columns"));
8043 struct ccase **rows = NULL;
8044 size_t allocated_rows = 0;
8048 while ((c = casereader_read (r)) != NULL)
8050 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
8060 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
8061 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
8062 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
8065 if (change != NOTHING_CHANGED)
8067 matrix_mget_commit_var (rows, n_rows, d,
8068 rowtype_, &mget->rowtypes,
8079 if (n_rows >= allocated_rows)
8080 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
8083 if (change == SPLITS_CHANGED)
8089 /* Reset the factor number, if there are factors. */
8093 if (row_has_factors)
8099 else if (change == FACTORS_CHANGED)
8101 if (row_has_factors)
8107 matrix_mget_commit_var (rows, n_rows, d,
8108 rowtype_, &mget->rowtypes,
8120 if (var_dimension->n_leaves)
8121 pivot_table_submit (pt);
8123 pivot_table_unref (pt);
8127 matrix_cmd_execute_mget (struct mget_command *mget)
8129 struct casereader *r;
8130 struct dictionary *d;
8131 if (matrix_open_casereader ("MGET", mget->file, mget->encoding,
8132 mget->state->dataset, &r, &d))
8134 matrix_cmd_execute_mget__ (mget, r, d);
8135 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
8140 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
8142 if (!lex_force_id (s->lexer))
8145 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
8147 *varp = matrix_var_create (s, lex_tokss (s->lexer));
8152 static struct matrix_cmd *
8153 matrix_parse_eigen (struct matrix_state *s)
8155 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
8156 *cmd = (struct matrix_cmd) {
8158 .eigen = { .expr = NULL }
8161 struct eigen_command *eigen = &cmd->eigen;
8162 if (!lex_force_match (s->lexer, T_LPAREN))
8164 eigen->expr = matrix_parse_expr (s);
8166 || !lex_force_match (s->lexer, T_COMMA)
8167 || !matrix_parse_dst_var (s, &eigen->evec)
8168 || !lex_force_match (s->lexer, T_COMMA)
8169 || !matrix_parse_dst_var (s, &eigen->eval)
8170 || !lex_force_match (s->lexer, T_RPAREN))
8176 matrix_cmd_destroy (cmd);
8181 matrix_cmd_execute_eigen (struct eigen_command *eigen)
8183 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
8186 if (!is_symmetric (A))
8188 msg (SE, _("Argument of EIGEN must be symmetric."));
8189 gsl_matrix_free (A);
8193 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
8194 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
8195 gsl_vector v_eval = to_vector (eval);
8196 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
8197 gsl_eigen_symmv (A, &v_eval, evec, w);
8198 gsl_eigen_symmv_free (w);
8200 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
8202 gsl_matrix_free (A);
8204 gsl_matrix_free (eigen->eval->value);
8205 eigen->eval->value = eval;
8207 gsl_matrix_free (eigen->evec->value);
8208 eigen->evec->value = evec;
8211 static struct matrix_cmd *
8212 matrix_parse_setdiag (struct matrix_state *s)
8214 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
8215 *cmd = (struct matrix_cmd) {
8216 .type = MCMD_SETDIAG,
8217 .setdiag = { .dst = NULL }
8220 struct setdiag_command *setdiag = &cmd->setdiag;
8221 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
8224 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
8227 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
8232 if (!lex_force_match (s->lexer, T_COMMA))
8235 setdiag->expr = matrix_parse_expr (s);
8239 if (!lex_force_match (s->lexer, T_RPAREN))
8245 matrix_cmd_destroy (cmd);
8250 matrix_cmd_execute_setdiag (struct setdiag_command *setdiag)
8252 gsl_matrix *dst = setdiag->dst->value;
8255 msg (SE, _("SETDIAG destination matrix %s is uninitialized."),
8256 setdiag->dst->name);
8260 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
8264 size_t n = MIN (dst->size1, dst->size2);
8265 if (is_scalar (src))
8267 double d = to_scalar (src);
8268 for (size_t i = 0; i < n; i++)
8269 gsl_matrix_set (dst, i, i, d);
8271 else if (is_vector (src))
8273 gsl_vector v = to_vector (src);
8274 for (size_t i = 0; i < n && i < v.size; i++)
8275 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
8278 msg (SE, _("SETDIAG argument 2 must be a scalar or a vector, "
8279 "not a %zu×%zu matrix."),
8280 src->size1, src->size2);
8281 gsl_matrix_free (src);
8284 static struct matrix_cmd *
8285 matrix_parse_svd (struct matrix_state *s)
8287 struct matrix_cmd *cmd = xmalloc (sizeof *cmd);
8288 *cmd = (struct matrix_cmd) {
8290 .svd = { .expr = NULL }
8293 struct svd_command *svd = &cmd->svd;
8294 if (!lex_force_match (s->lexer, T_LPAREN))
8296 svd->expr = matrix_parse_expr (s);
8298 || !lex_force_match (s->lexer, T_COMMA)
8299 || !matrix_parse_dst_var (s, &svd->u)
8300 || !lex_force_match (s->lexer, T_COMMA)
8301 || !matrix_parse_dst_var (s, &svd->s)
8302 || !lex_force_match (s->lexer, T_COMMA)
8303 || !matrix_parse_dst_var (s, &svd->v)
8304 || !lex_force_match (s->lexer, T_RPAREN))
8310 matrix_cmd_destroy (cmd);
8315 matrix_cmd_execute_svd (struct svd_command *svd)
8317 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8321 if (m->size1 >= m->size2)
8324 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8325 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8326 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8327 gsl_vector *work = gsl_vector_alloc (A->size2);
8328 gsl_linalg_SV_decomp (A, V, &Sv, work);
8329 gsl_vector_free (work);
8331 matrix_var_set (svd->u, A);
8332 matrix_var_set (svd->s, S);
8333 matrix_var_set (svd->v, V);
8337 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8338 gsl_matrix_transpose_memcpy (At, m);
8339 gsl_matrix_free (m);
8341 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8342 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8343 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8344 gsl_vector *work = gsl_vector_alloc (At->size2);
8345 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8346 gsl_vector_free (work);
8348 matrix_var_set (svd->v, At);
8349 matrix_var_set (svd->s, St);
8350 matrix_var_set (svd->u, Vt);
8355 matrix_cmd_execute (struct matrix_cmd *cmd)
8360 matrix_cmd_execute_compute (&cmd->compute);
8364 matrix_cmd_execute_print (&cmd->print);
8368 return matrix_cmd_execute_do_if (&cmd->do_if);
8371 matrix_cmd_execute_loop (&cmd->loop);
8378 matrix_cmd_execute_display (&cmd->display);
8382 matrix_cmd_execute_release (&cmd->release);
8386 matrix_cmd_execute_save (&cmd->save);
8390 matrix_cmd_execute_read (&cmd->read);
8394 matrix_cmd_execute_write (&cmd->write);
8398 matrix_cmd_execute_get (&cmd->get);
8402 matrix_cmd_execute_msave (&cmd->msave);
8406 matrix_cmd_execute_mget (&cmd->mget);
8410 matrix_cmd_execute_eigen (&cmd->eigen);
8414 matrix_cmd_execute_setdiag (&cmd->setdiag);
8418 matrix_cmd_execute_svd (&cmd->svd);
8426 matrix_cmds_uninit (struct matrix_cmds *cmds)
8428 for (size_t i = 0; i < cmds->n; i++)
8429 matrix_cmd_destroy (cmds->commands[i]);
8430 free (cmds->commands);
8434 matrix_cmd_destroy (struct matrix_cmd *cmd)
8442 matrix_lvalue_destroy (cmd->compute.lvalue);
8443 matrix_expr_destroy (cmd->compute.rvalue);
8447 matrix_expr_destroy (cmd->print.expression);
8448 free (cmd->print.title);
8449 print_labels_destroy (cmd->print.rlabels);
8450 print_labels_destroy (cmd->print.clabels);
8454 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8456 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8457 matrix_cmds_uninit (&cmd->do_if.clauses[i].commands);
8459 free (cmd->do_if.clauses);
8463 matrix_expr_destroy (cmd->loop.start);
8464 matrix_expr_destroy (cmd->loop.end);
8465 matrix_expr_destroy (cmd->loop.increment);
8466 matrix_expr_destroy (cmd->loop.top_condition);
8467 matrix_expr_destroy (cmd->loop.bottom_condition);
8468 matrix_cmds_uninit (&cmd->loop.commands);
8478 free (cmd->release.vars);
8482 matrix_expr_destroy (cmd->save.expression);
8486 matrix_lvalue_destroy (cmd->read.dst);
8487 matrix_expr_destroy (cmd->read.size);
8491 matrix_expr_destroy (cmd->write.expression);
8492 free (cmd->write.format);
8496 matrix_lvalue_destroy (cmd->get.dst);
8497 fh_unref (cmd->get.file);
8498 free (cmd->get.encoding);
8499 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8503 matrix_expr_destroy (cmd->msave.expr);
8507 fh_unref (cmd->mget.file);
8508 stringi_set_destroy (&cmd->mget.rowtypes);
8512 matrix_expr_destroy (cmd->eigen.expr);
8516 matrix_expr_destroy (cmd->setdiag.expr);
8520 matrix_expr_destroy (cmd->svd.expr);
8526 struct matrix_command_name
8529 struct matrix_cmd *(*parse) (struct matrix_state *);
8532 static const struct matrix_command_name *
8533 matrix_parse_command_name (struct lexer *lexer)
8535 static const struct matrix_command_name commands[] = {
8536 { "COMPUTE", matrix_parse_compute },
8537 { "CALL EIGEN", matrix_parse_eigen },
8538 { "CALL SETDIAG", matrix_parse_setdiag },
8539 { "CALL SVD", matrix_parse_svd },
8540 { "PRINT", matrix_parse_print },
8541 { "DO IF", matrix_parse_do_if },
8542 { "LOOP", matrix_parse_loop },
8543 { "BREAK", matrix_parse_break },
8544 { "READ", matrix_parse_read },
8545 { "WRITE", matrix_parse_write },
8546 { "GET", matrix_parse_get },
8547 { "SAVE", matrix_parse_save },
8548 { "MGET", matrix_parse_mget },
8549 { "MSAVE", matrix_parse_msave },
8550 { "DISPLAY", matrix_parse_display },
8551 { "RELEASE", matrix_parse_release },
8553 static size_t n = sizeof commands / sizeof *commands;
8555 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8556 if (lex_match_phrase (lexer, c->name))
8561 static struct matrix_cmd *
8562 matrix_parse_command (struct matrix_state *s)
8564 size_t nesting_level = SIZE_MAX;
8566 struct matrix_cmd *c = NULL;
8567 const struct matrix_command_name *cmd = matrix_parse_command_name (s->lexer);
8569 lex_error (s->lexer, _("Unknown matrix command."));
8570 else if (!cmd->parse)
8571 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8575 nesting_level = output_open_group (
8576 group_item_create_nocopy (utf8_to_title (cmd->name),
8577 utf8_to_title (cmd->name)));
8582 lex_end_of_command (s->lexer);
8583 lex_discard_rest_of_command (s->lexer);
8584 if (nesting_level != SIZE_MAX)
8585 output_close_groups (nesting_level);
8591 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8593 if (!lex_force_match (lexer, T_ENDCMD))
8596 struct matrix_state state = {
8598 .session = dataset_session (ds),
8600 .vars = HMAP_INITIALIZER (state.vars),
8605 while (lex_match (lexer, T_ENDCMD))
8607 if (lex_token (lexer) == T_STOP)
8609 msg (SE, _("Unexpected end of input expecting matrix command."));
8613 if (lex_match_phrase (lexer, "END MATRIX"))
8616 struct matrix_cmd *c = matrix_parse_command (&state);
8619 matrix_cmd_execute (c);
8620 matrix_cmd_destroy (c);
8624 struct matrix_var *var, *next;
8625 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8628 gsl_matrix_free (var->value);
8629 hmap_delete (&state.vars, &var->hmap_node);
8632 hmap_destroy (&state.vars);
8633 msave_common_destroy (state.common);
8634 fh_unref (state.prev_read_file);
8635 for (size_t i = 0; i < state.n_read_files; i++)
8636 read_file_destroy (state.read_files[i]);
8637 free (state.read_files);
8638 fh_unref (state.prev_write_file);
8639 for (size_t i = 0; i < state.n_write_files; i++)
8640 write_file_destroy (state.write_files[i]);
8641 free (state.write_files);
8642 fh_unref (state.prev_save_file);
8643 for (size_t i = 0; i < state.n_save_files; i++)
8644 save_file_destroy (state.save_files[i]);
8645 free (state.save_files);