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/commands/data-reader.h"
42 #include "language/commands/data-writer.h"
43 #include "language/commands/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 /* A variable in the matrix language. */
79 struct hmap_node hmap_node; /* In matrix_state's 'vars' hmap. */
80 char *name; /* UTF-8. */
81 gsl_matrix *value; /* NULL, if the variable is uninitialized. */
84 /* All the MSAVE commands within a matrix program share common configuration,
85 provided by the first MSAVE command within the program. This structure
86 encapsulates this configuration. */
89 /* Common configuration for all MSAVEs. */
90 struct msg_location *location; /* Range of lines for first MSAVE. */
91 struct file_handle *outfile; /* Output file for all the MSAVEs. */
92 struct msg_location *outfile_location;
93 struct string_array variables; /* VARIABLES subcommand. */
94 struct msg_location *variables_location;
95 struct string_array fnames; /* FNAMES subcommand. */
96 struct msg_location *fnames_location;
97 struct string_array snames; /* SNAMES subcommand. */
98 struct msg_location *snames_location;
100 /* Collects and owns factors and splits. The individual msave_command
101 structs point to these but do not own them. (This is because factors
102 and splits can be carried over from one MSAVE to the next, so it's
103 easiest to just take the most recent.) */
104 struct matrix_expr **factors;
105 size_t n_factors, allocated_factors;
106 struct matrix_expr **splits;
107 size_t n_splits, allocated_splits;
109 /* Execution state. */
110 struct dictionary *dict;
111 struct casewriter *writer;
114 /* A file used by one or more READ commands. */
118 struct file_handle *file;
120 /* Execution state. */
121 struct dfm_reader *reader;
125 static struct read_file *read_file_create (struct matrix_state *,
126 struct file_handle *);
127 static struct dfm_reader *read_file_open (struct read_file *);
129 /* A file used by one or more WRITE comamnds. */
133 struct file_handle *file;
135 /* Execution state. */
136 struct dfm_writer *writer;
138 struct u8_line *held; /* Output held by a previous WRITE with /HOLD. */
141 static struct write_file *write_file_create (struct matrix_state *,
142 struct file_handle *);
143 static struct dfm_writer *write_file_open (struct write_file *);
144 static void write_file_destroy (struct write_file *);
146 /* A file used by one or more SAVE commands. */
150 struct file_handle *file;
151 struct dataset *dataset;
152 struct string_array variables;
153 struct matrix_expr *names;
154 struct stringi_set strings;
156 /* Execution state. */
158 struct casewriter *writer;
159 struct dictionary *dict;
160 struct msg_location *location;
163 /* State of an entire matrix program. */
166 /* State passed into MATRIX from outside. */
167 struct dataset *dataset;
168 struct session *session;
171 /* Matrix program's own state. */
172 struct hmap vars; /* Dictionary of matrix variables. */
173 bool in_loop; /* True if parsing within a LOOP. */
176 struct msave_common *msave_common;
179 struct file_handle *prev_read_file;
180 struct read_file **read_files;
184 struct file_handle *prev_write_file;
185 struct write_file **write_files;
186 size_t n_write_files;
189 struct file_handle *prev_save_file;
190 struct save_file **save_files;
194 /* Finds and returns the variable with the given NAME (case-insensitive) within
195 S, if there is one, or a null pointer if there is not. */
196 static struct matrix_var *
197 matrix_var_lookup (struct matrix_state *s, struct substring name)
199 struct matrix_var *var;
201 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
202 utf8_hash_case_substring (name, 0), &s->vars)
203 if (!utf8_sscasecmp (ss_cstr (var->name), name))
209 /* Creates and returns a new variable named NAME within S. There must not
210 already be a variable with the same (case-insensitive) name. The variable
211 is created uninitialized. */
212 static struct matrix_var *
213 matrix_var_create (struct matrix_state *s, struct substring name)
215 struct matrix_var *var = xmalloc (sizeof *var);
216 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
217 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
221 /* Replaces VAR's value by VALUE. Takes ownership of VALUE. */
223 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
225 gsl_matrix_free (var->value);
229 /* Matrix function catalog. */
231 /* The third argument to F() is a "prototype". For most prototypes, the first
232 letter (before the _) represents the return type and each other letter
233 (after the _) is an argument type. The types are:
235 - "m": A matrix of unrestricted dimensions.
239 - "v": A row or column vector.
241 - "e": Primarily for the first argument, this is a matrix with
242 unrestricted dimensions treated elementwise. Each element in the matrix
243 is passed to the implementation function separately.
245 - "n": This gets passed the "const struct matrix_expr *" that represents
246 the expression. This allows the evaluation function to grab the source
247 location of arguments so that it can report accurate error locations.
248 This type doesn't correspond to an argument passed in by the user.
250 The fourth argument is an optional constraints string. For this purpose the
251 first argument is named "a", the second "b", and so on. The following kinds
252 of constraints are supported. For matrix arguments, the constraints are
253 applied to each value in the matrix separately:
255 - "a(0,1)" or "a[0,1]": 0 < a < 1 or 0 <= a <= 1, respectively. Any
256 integer may substitute for 0 and 1. Half-open constraints (] and [) are
259 - "ai": Restrict a to integer values.
261 - "a>0", "a<0", "a>=0", "a<=0", "a!=0".
263 - "a<b", "a>b", "a<=b", "a>=b", "b!=0".
265 #define MATRIX_FUNCTIONS \
266 F(ABS, "ABS", m_e, NULL) \
267 F(ALL, "ALL", d_m, NULL) \
268 F(ANY, "ANY", d_m, NULL) \
269 F(ARSIN, "ARSIN", m_e, "a[-1,1]") \
270 F(ARTAN, "ARTAN", m_e, NULL) \
271 F(BLOCK, "BLOCK", m_any, NULL) \
272 F(CHOL, "CHOL", m_mn, NULL) \
273 F(CMIN, "CMIN", m_m, NULL) \
274 F(CMAX, "CMAX", m_m, NULL) \
275 F(COS, "COS", m_e, NULL) \
276 F(CSSQ, "CSSQ", m_m, NULL) \
277 F(CSUM, "CSUM", m_m, NULL) \
278 F(DESIGN, "DESIGN", m_mn, NULL) \
279 F(DET, "DET", d_m, NULL) \
280 F(DIAG, "DIAG", m_m, NULL) \
281 F(EVAL, "EVAL", m_mn, NULL) \
282 F(EXP, "EXP", m_e, NULL) \
283 F(GINV, "GINV", m_m, NULL) \
284 F(GRADE, "GRADE", m_m, NULL) \
285 F(GSCH, "GSCH", m_mn, NULL) \
286 F(IDENT, "IDENT", IDENT, NULL) \
287 F(INV, "INV", m_m, NULL) \
288 F(KRONEKER, "KRONEKER", m_mm, NULL) \
289 F(LG10, "LG10", m_e, "a>0") \
290 F(LN, "LN", m_e, "a>0") \
291 F(MAGIC, "MAGIC", m_d, "ai>=3") \
292 F(MAKE, "MAKE", m_ddd, "ai>=0 bi>=0") \
293 F(MDIAG, "MDIAG", m_v, NULL) \
294 F(MMAX, "MMAX", d_m, NULL) \
295 F(MMIN, "MMIN", d_m, NULL) \
296 F(MOD, "MOD", m_md, "b!=0") \
297 F(MSSQ, "MSSQ", d_m, NULL) \
298 F(MSUM, "MSUM", d_m, NULL) \
299 F(NCOL, "NCOL", d_m, NULL) \
300 F(NROW, "NROW", d_m, NULL) \
301 F(RANK, "RANK", d_m, NULL) \
302 F(RESHAPE, "RESHAPE", m_mddn, NULL) \
303 F(RMAX, "RMAX", m_m, NULL) \
304 F(RMIN, "RMIN", m_m, NULL) \
305 F(RND, "RND", m_e, NULL) \
306 F(RNKORDER, "RNKORDER", m_m, NULL) \
307 F(RSSQ, "RSSQ", m_m, NULL) \
308 F(RSUM, "RSUM", m_m, NULL) \
309 F(SIN, "SIN", m_e, NULL) \
310 F(SOLVE, "SOLVE", m_mmn, NULL) \
311 F(SQRT, "SQRT", m_e, "a>=0") \
312 F(SSCP, "SSCP", m_m, NULL) \
313 F(SVAL, "SVAL", m_m, NULL) \
314 F(SWEEP, "SWEEP", m_mdn, NULL) \
315 F(T, "T", m_m, NULL) \
316 F(TRACE, "TRACE", d_m, NULL) \
317 F(TRANSPOS, "TRANSPOS", m_m, NULL) \
318 F(TRUNC, "TRUNC", m_e, NULL) \
319 F(UNIFORM, "UNIFORM", m_ddn, "ai>=0 bi>=0") \
320 F(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0") \
321 F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0") \
322 F(IDF_BETA, "IDF.BETA", m_edd, "a[0,1] b>0 c>0") \
323 F(RV_BETA, "RV.BETA", d_dd, "a>0 b>0") \
324 F(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
325 F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
326 F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]") \
327 F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]") \
328 F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0") \
329 F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0") \
330 F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "c>0") \
331 F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \
332 F(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0") \
333 F(CHICDF, "CHICDF", m_ed, "a>=0 b>0") \
334 F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0") \
335 F(PDF_CHISQ, "PDF.CHISQ", m_ed, "a>=0 b>0") \
336 F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \
337 F(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0") \
338 F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0") \
339 F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0") \
340 F(PDF_EXP, "PDF.EXP", m_ed, "a>=0 b>0") \
341 F(RV_EXP, "RV.EXP", d_d, "a>0") \
342 F(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0") \
343 F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \
344 F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0") \
345 F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0") \
346 F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0") \
347 F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0") \
348 F(RV_F, "RV.F", d_dd, "a>0 b>0") \
349 F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0") \
350 F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
351 F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0") \
352 F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
353 F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0") \
354 F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL) \
355 F(RV_LANDAU, "RV.LANDAU", d_none, NULL) \
356 F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0") \
357 F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0") \
358 F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0") \
359 F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0") \
360 F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]") \
361 F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]") \
362 F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0") \
363 F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0") \
364 F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0") \
365 F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0") \
366 F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
367 F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0") \
368 F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
369 F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0") \
370 F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0") \
371 F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0") \
372 F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0") \
373 F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0") \
374 F(CDFNORM, "CDFNORM", m_e, NULL) \
375 F(PROBIT, "PROBIT", m_e, "a(0,1)") \
376 F(NORMAL, "NORMAL", m_e, "a>0") \
377 F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0") \
378 F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0") \
379 F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0") \
380 F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0") \
381 F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0") \
382 F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0") \
383 F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0") \
384 F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0") \
385 F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0") \
386 F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0") \
387 F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL) \
388 F(RV_RTAIL, "RV.RTAIL", d_dd, NULL) \
389 F(CDF_T, "CDF.T", m_ed, "b>0") \
390 F(TCDF, "TCDF", m_ed, "b>0") \
391 F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0") \
392 F(PDF_T, "PDF.T", m_ed, "b>0") \
393 F(RV_T, "RV.T", d_d, "a>0") \
394 F(CDF_T1G, "CDF.T1G", m_edd, NULL) \
395 F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)") \
396 F(PDF_T1G, "PDF.T1G", m_edd, NULL) \
397 F(RV_T1G, "RV.T1G", d_dd, NULL) \
398 F(CDF_T2G, "CDF.T2G", m_edd, NULL) \
399 F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)") \
400 F(PDF_T2G, "PDF.T2G", m_edd, NULL) \
401 F(RV_T2G, "RV.T2G", d_dd, NULL) \
402 F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c") \
403 F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c") \
404 F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c") \
405 F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b") \
406 F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
407 F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0") \
408 F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
409 F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0") \
410 F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
411 F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
412 F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]") \
413 F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]") \
414 F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]") \
415 F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]") \
416 F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]") \
417 F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]") \
418 F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]") \
419 F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
420 F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
421 F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a") \
422 F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]") \
423 F(RV_LOG, "RV.LOG", d_d, "a(0,1]") \
424 F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
425 F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
426 F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]") \
427 F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0") \
428 F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0") \
429 F(RV_POISSON, "RV.POISSON", d_d, "a>0")
431 /* Properties of a matrix function.
433 These come straight from the macro invocations above. */
434 struct matrix_function_properties
437 const char *constraints;
440 /* Minimum and maximum argument counts for each matrix function prototype. */
441 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
442 enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
443 enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
444 enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 };
445 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
446 enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 };
447 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
448 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
449 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
450 enum { m_ddn_MIN_ARGS = 2, m_ddn_MAX_ARGS = 2 };
451 enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 };
452 enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 };
453 enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
454 enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
455 enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
456 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
457 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
458 enum { m_mddn_MIN_ARGS = 3, m_mddn_MAX_ARGS = 3 };
459 enum { m_mdn_MIN_ARGS = 2, m_mdn_MAX_ARGS = 2 };
460 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
461 enum { m_mmn_MIN_ARGS = 2, m_mmn_MAX_ARGS = 2 };
462 enum { m_mn_MIN_ARGS = 1, m_mn_MAX_ARGS = 1 };
463 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
465 /* C function prototype for each matrix function prototype. */
466 typedef double matrix_proto_d_none (void);
467 typedef double matrix_proto_d_d (double);
468 typedef double matrix_proto_d_dd (double, double);
469 typedef double matrix_proto_d_dd (double, double);
470 typedef double matrix_proto_d_ddd (double, double, double);
471 typedef gsl_matrix *matrix_proto_m_d (double);
472 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
473 typedef gsl_matrix *matrix_proto_m_ddn (double, double,
474 const struct matrix_expr *);
475 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
476 typedef gsl_matrix *matrix_proto_m_mn (gsl_matrix *,
477 const struct matrix_expr *);
478 typedef double matrix_proto_m_e (double);
479 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
480 typedef gsl_matrix *matrix_proto_m_mdn (gsl_matrix *, double,
481 const struct matrix_expr *);
482 typedef double matrix_proto_m_ed (double, double);
483 typedef gsl_matrix *matrix_proto_m_mddn (gsl_matrix *, double, double,
484 const struct matrix_expr *);
485 typedef double matrix_proto_m_edd (double, double, double);
486 typedef double matrix_proto_m_eddd (double, double, double, double);
487 typedef double matrix_proto_m_eed (double, double, double);
488 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
489 typedef gsl_matrix *matrix_proto_m_mmn (gsl_matrix *, gsl_matrix *,
490 const struct matrix_expr *);
491 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
492 typedef double matrix_proto_d_m (gsl_matrix *);
493 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
494 typedef gsl_matrix *matrix_proto_IDENT (double, double);
496 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
497 static matrix_proto_##PROTO matrix_eval_##ENUM;
501 /* Matrix expression data structure and parsing. */
503 /* A node in a matrix expression. */
509 #define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
513 /* Elementwise and scalar arithmetic. */
514 MOP_NEGATE, /* unary - */
515 MOP_ADD_ELEMS, /* + */
516 MOP_SUB_ELEMS, /* - */
517 MOP_MUL_ELEMS, /* &* */
518 MOP_DIV_ELEMS, /* / and &/ */
519 MOP_EXP_ELEMS, /* &** */
521 MOP_SEQ_BY, /* a:b:c */
523 /* Matrix arithmetic. */
525 MOP_EXP_MAT, /* ** */
542 MOP_PASTE_HORZ, /* a, b, c, ... */
543 MOP_PASTE_VERT, /* a; b; c; ... */
547 MOP_VEC_INDEX, /* x(y) */
548 MOP_VEC_ALL, /* x(:) */
549 MOP_MAT_INDEX, /* x(y,z) */
550 MOP_ROW_INDEX, /* x(y,:) */
551 MOP_COL_INDEX, /* x(:,z) */
558 MOP_EOF, /* EOF('file') */
564 /* Nonterminal expression nodes. */
567 struct matrix_expr **subs;
571 /* Terminal expression nodes. */
572 double number; /* MOP_NUMBER. */
573 struct matrix_var *variable; /* MOP_VARIABLE. */
574 struct read_file *eof; /* MOP_EOF. */
577 /* The syntax location corresponding to this expression node, for use in
578 error messages. This is always nonnull for terminal expression nodes.
579 For most others, it is null because it can be computed lazily if and
582 Use matrix_expr_location() instead of using this member directly, so
583 that it gets computed lazily if needed. */
584 struct msg_location *location;
588 matrix_expr_location__ (const struct matrix_expr *e,
589 const struct msg_location **minp,
590 const struct msg_location **maxp)
592 struct msg_location *loc = e->location;
595 const struct msg_location *min = *minp;
598 || loc->start.line < min->start.line
599 || (loc->start.line == min->start.line
600 && loc->start.column < min->start.column)))
603 const struct msg_location *max = *maxp;
606 || loc->end.line > max->end.line
607 || (loc->end.line == max->end.line
608 && loc->end.column > max->end.column)))
614 assert (e->op != MOP_NUMBER && e->op != MOP_VARIABLE && e->op != MOP_EOF);
615 for (size_t i = 0; i < e->n_subs; i++)
616 matrix_expr_location__ (e->subs[i], minp, maxp);
619 /* Returns the source code location corresponding to expression E, computing it
621 static const struct msg_location *
622 matrix_expr_location (const struct matrix_expr *e_)
624 struct matrix_expr *e = CONST_CAST (struct matrix_expr *, e_);
630 const struct msg_location *min = NULL;
631 const struct msg_location *max = NULL;
632 matrix_expr_location__ (e, &min, &max);
635 e->location = msg_location_dup (min);
636 e->location->end = max->end;
642 /* Sets e->location to the tokens in S's lexer from offset START_OFS to the
643 token before the current one. Has no effect if E already has a location or
646 matrix_expr_add_location (struct matrix_state *s, int start_ofs,
647 struct matrix_expr *e)
649 if (e && !e->location)
650 e->location = lex_ofs_location (s->lexer, start_ofs,
651 lex_ofs (s->lexer) - 1);
654 /* Frees E and all the data and sub-expressions that it references. */
656 matrix_expr_destroy (struct matrix_expr *e)
663 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
694 for (size_t i = 0; i < e->n_subs; i++)
695 matrix_expr_destroy (e->subs[i]);
704 msg_location_destroy (e->location);
708 /* Creates and returns a new matrix_expr with type OP, which must be a
709 nonterminal type. Initializes the new matrix_expr with the N_SUBS
710 expressions in SUBS as subexpressions. */
711 static struct matrix_expr *
712 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
715 struct matrix_expr *e = xmalloc (sizeof *e);
716 *e = (struct matrix_expr) {
718 .subs = xmemdup (subs, n_subs * sizeof *subs),
724 static struct matrix_expr *
725 matrix_expr_create_0 (enum matrix_op op)
727 struct matrix_expr *sub;
728 return matrix_expr_create_subs (op, &sub, 0);
731 static struct matrix_expr *
732 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
734 return matrix_expr_create_subs (op, &sub, 1);
737 static struct matrix_expr *
738 matrix_expr_create_2 (enum matrix_op op,
739 struct matrix_expr *sub0, struct matrix_expr *sub1)
741 struct matrix_expr *subs[] = { sub0, sub1 };
742 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
745 static struct matrix_expr *
746 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
747 struct matrix_expr *sub1, struct matrix_expr *sub2)
749 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
750 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
753 /* Creates and returns a new MOP_NUMBER expression node to contain NUMBER. */
754 static struct matrix_expr *
755 matrix_expr_create_number (double number)
757 struct matrix_expr *e = xmalloc (sizeof *e);
758 *e = (struct matrix_expr) {
765 static struct matrix_expr *matrix_expr_parse (struct matrix_state *);
767 /* A binary operator for matrix_parse_binary_operator(). */
768 struct matrix_operator_syntax
770 /* Exactly one of these specifies the operator syntax. */
771 enum token_type token; /* A token, e.g. T_ASTERISK. */
772 const char *id; /* An identifier, e.g. "XOR". */
773 const char *phrase; /* A token phrase, e.g. "&**". */
775 /* The matrix operator corresponding to the syntax. */
780 matrix_operator_syntax_match (struct lexer *lexer,
781 const struct matrix_operator_syntax *syntax,
782 size_t n_syntax, enum matrix_op *op)
784 const struct matrix_operator_syntax *end = &syntax[n_syntax];
785 for (const struct matrix_operator_syntax *syn = syntax; syn < end; syn++)
786 if (syn->id ? lex_match_id (lexer, syn->id)
787 : syn->phrase ? lex_match_phrase (lexer, syn->phrase)
788 : lex_match (lexer, syn->token))
796 /* Parses a binary operator level in the recursive descent parser, returning a
797 matrix expression if successful or a null pointer otherwise. PARSE_NEXT
798 must be the function to parse the next level of precedence. The N_SYNTAX
799 elements of SYNTAX must specify the syntax and matrix_expr node type to
800 parse at this level. */
801 static struct matrix_expr *
802 matrix_parse_binary_operator (
803 struct matrix_state *s,
804 struct matrix_expr *(*parse_next) (struct matrix_state *),
805 const struct matrix_operator_syntax *syntax, size_t n_syntax)
807 struct matrix_expr *lhs = parse_next (s);
814 if (!matrix_operator_syntax_match (s->lexer, syntax, n_syntax, &op))
817 struct matrix_expr *rhs = parse_next (s);
820 matrix_expr_destroy (lhs);
823 lhs = matrix_expr_create_2 (op, lhs, rhs);
827 /* Parses a comma-separated list of expressions within {}, transforming them
828 into MOP_PASTE_HORZ operators. Returns the new expression or NULL on
830 static struct matrix_expr *
831 matrix_parse_curly_comma (struct matrix_state *s)
833 static const struct matrix_operator_syntax op = {
834 .token = T_COMMA, .op = MOP_PASTE_HORZ
836 return matrix_parse_binary_operator (s, matrix_expr_parse, &op, 1);
839 /* Parses a semicolon-separated list of expressions within {}, transforming
840 them into MOP_PASTE_VERT operators. Returns the new expression or NULL on
842 static struct matrix_expr *
843 matrix_parse_curly_semi (struct matrix_state *s)
845 if (lex_token (s->lexer) == T_RCURLY)
847 /* {} is a special case for a 0×0 matrix. */
848 return matrix_expr_create_0 (MOP_EMPTY);
851 static const struct matrix_operator_syntax op = {
852 .token = T_SEMICOLON, .op = MOP_PASTE_VERT
854 return matrix_parse_binary_operator (s, matrix_parse_curly_comma, &op, 1);
857 struct matrix_function
861 size_t min_args, max_args;
864 static struct matrix_expr *matrix_expr_parse (struct matrix_state *);
867 word_matches (const char **test, const char **name)
869 size_t test_len = strcspn (*test, ".");
870 size_t name_len = strcspn (*name, ".");
871 if (test_len == name_len)
873 if (buf_compare_case (*test, *name, test_len))
876 else if (test_len < 3 || test_len > name_len)
880 if (buf_compare_case (*test, *name, test_len))
886 if (**test != **name)
897 /* Returns 0 if TOKEN and FUNC do not match,
898 1 if TOKEN is an acceptable abbreviation for FUNC,
899 2 if TOKEN equals FUNC. */
901 compare_function_names (const char *token_, const char *func_)
903 const char *token = token_;
904 const char *func = func_;
905 while (*token || *func)
906 if (!word_matches (&token, &func))
908 return !c_strcasecmp (token_, func_) ? 2 : 1;
911 static const struct matrix_function *
912 matrix_parse_function_name (const char *token)
914 static const struct matrix_function functions[] =
916 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
917 { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
921 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
923 for (size_t i = 0; i < N_FUNCTIONS; i++)
925 if (compare_function_names (token, functions[i].name) > 0)
926 return &functions[i];
932 matrix_parse_function (struct matrix_state *s, const char *token,
933 struct matrix_expr **exprp)
936 if (lex_next_token (s->lexer, 1) != T_LPAREN)
939 int start_ofs = lex_ofs (s->lexer);
940 if (lex_match_id (s->lexer, "EOF"))
943 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
947 if (!lex_force_match (s->lexer, T_RPAREN))
953 struct read_file *rf = read_file_create (s, fh);
955 struct matrix_expr *e = xmalloc (sizeof *e);
956 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
957 matrix_expr_add_location (s, start_ofs, e);
962 const struct matrix_function *f = matrix_parse_function_name (token);
966 struct matrix_expr *e = xmalloc (sizeof *e);
967 *e = (struct matrix_expr) { .op = f->op };
969 lex_get_n (s->lexer, 2);
970 if (lex_token (s->lexer) != T_RPAREN)
972 size_t allocated_subs = 0;
975 struct matrix_expr *sub = matrix_expr_parse (s);
979 if (e->n_subs >= allocated_subs)
980 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
981 e->subs[e->n_subs++] = sub;
983 while (lex_match (s->lexer, T_COMMA));
985 if (!lex_force_match (s->lexer, T_RPAREN))
988 if (e->n_subs < f->min_args || e->n_subs > f->max_args)
990 if (f->min_args == f->max_args)
991 msg_at (SE, e->location,
992 ngettext ("Matrix function %s requires %zu argument.",
993 "Matrix function %s requires %zu arguments.",
995 f->name, f->min_args);
996 else if (f->min_args == 1 && f->max_args == 2)
997 msg_at (SE, e->location,
998 ngettext ("Matrix function %s requires 1 or 2 arguments, "
999 "but %zu was provided.",
1000 "Matrix function %s requires 1 or 2 arguments, "
1001 "but %zu were provided.",
1003 f->name, e->n_subs);
1004 else if (f->min_args == 1 && f->max_args == INT_MAX)
1005 msg_at (SE, e->location,
1006 _("Matrix function %s requires at least one argument."),
1014 matrix_expr_add_location (s, start_ofs, e);
1020 matrix_expr_destroy (e);
1024 static struct matrix_expr *
1025 matrix_parse_primary__ (struct matrix_state *s)
1027 if (lex_is_number (s->lexer))
1029 double number = lex_number (s->lexer);
1032 return matrix_expr_create_number (number);
1034 else if (lex_is_string (s->lexer))
1036 char string[sizeof (double)];
1037 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1041 memcpy (&number, string, sizeof number);
1043 return matrix_expr_create_number (number);
1045 else if (lex_match (s->lexer, T_LPAREN))
1047 struct matrix_expr *e = matrix_expr_parse (s);
1048 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1050 matrix_expr_destroy (e);
1055 else if (lex_match (s->lexer, T_LCURLY))
1057 struct matrix_expr *e = matrix_parse_curly_semi (s);
1058 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1060 matrix_expr_destroy (e);
1065 else if (lex_token (s->lexer) == T_ID)
1067 struct matrix_expr *retval;
1068 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1071 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1074 lex_error (s->lexer, _("Unknown variable %s."),
1075 lex_tokcstr (s->lexer));
1080 struct matrix_expr *e = xmalloc (sizeof *e);
1081 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1084 else if (lex_token (s->lexer) == T_ALL)
1086 struct matrix_expr *retval;
1087 if (matrix_parse_function (s, "ALL", &retval))
1091 lex_error (s->lexer, _("Syntax error expecting matrix expression."));
1095 static struct matrix_expr *
1096 matrix_parse_primary (struct matrix_state *s)
1098 int start_ofs = lex_ofs (s->lexer);
1099 struct matrix_expr *e = matrix_parse_primary__ (s);
1100 matrix_expr_add_location (s, start_ofs, e);
1104 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1107 matrix_parse_index_expr (struct matrix_state *s,
1108 struct matrix_expr **indexp,
1109 struct msg_location **locationp)
1111 if (lex_match (s->lexer, T_COLON))
1114 *locationp = lex_get_location (s->lexer, -1, -1);
1120 *indexp = matrix_expr_parse (s);
1121 if (locationp && *indexp)
1122 *locationp = msg_location_dup (matrix_expr_location (*indexp));
1123 return *indexp != NULL;
1127 static struct matrix_expr *
1128 matrix_parse_postfix (struct matrix_state *s)
1130 struct matrix_expr *lhs = matrix_parse_primary (s);
1131 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1134 struct matrix_expr *i0;
1135 if (!matrix_parse_index_expr (s, &i0, NULL))
1137 matrix_expr_destroy (lhs);
1140 if (lex_match (s->lexer, T_RPAREN))
1142 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1143 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1144 else if (lex_match (s->lexer, T_COMMA))
1146 struct matrix_expr *i1;
1147 if (!matrix_parse_index_expr (s, &i1, NULL)
1148 || !lex_force_match (s->lexer, T_RPAREN))
1150 matrix_expr_destroy (lhs);
1151 matrix_expr_destroy (i0);
1152 matrix_expr_destroy (i1);
1155 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1156 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1157 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1162 lex_error_expecting (s->lexer, "`)'", "`,'");
1167 static struct matrix_expr *
1168 matrix_parse_unary (struct matrix_state *s)
1170 int start_ofs = lex_ofs (s->lexer);
1172 struct matrix_expr *e;
1173 if (lex_match (s->lexer, T_DASH))
1175 struct matrix_expr *sub = matrix_parse_unary (s);
1178 e = matrix_expr_create_1 (MOP_NEGATE, sub);
1180 else if (lex_match (s->lexer, T_PLUS))
1182 e = matrix_parse_unary (s);
1187 return matrix_parse_postfix (s);
1189 matrix_expr_add_location (s, start_ofs, e);
1190 e->location->start = lex_ofs_start_point (s->lexer, start_ofs);
1194 static struct matrix_expr *
1195 matrix_parse_seq (struct matrix_state *s)
1197 struct matrix_expr *start = matrix_parse_unary (s);
1198 if (!start || !lex_match (s->lexer, T_COLON))
1201 struct matrix_expr *end = matrix_parse_unary (s);
1204 matrix_expr_destroy (start);
1208 if (lex_match (s->lexer, T_COLON))
1210 struct matrix_expr *increment = matrix_parse_unary (s);
1213 matrix_expr_destroy (start);
1214 matrix_expr_destroy (end);
1217 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
1220 return matrix_expr_create_2 (MOP_SEQ, start, end);
1223 static struct matrix_expr *
1224 matrix_parse_exp (struct matrix_state *s)
1226 static const struct matrix_operator_syntax syntax[] = {
1227 { .token = T_EXP, .op = MOP_EXP_MAT },
1228 { .phrase = "&**", .op = MOP_EXP_ELEMS },
1230 size_t n_syntax = sizeof syntax / sizeof *syntax;
1232 return matrix_parse_binary_operator (s, matrix_parse_seq, syntax, n_syntax);
1235 static struct matrix_expr *
1236 matrix_parse_mul_div (struct matrix_state *s)
1238 static const struct matrix_operator_syntax syntax[] = {
1239 { .token = T_ASTERISK, .op = MOP_MUL_MAT },
1240 { .token = T_SLASH, .op = MOP_DIV_ELEMS },
1241 { .phrase = "&*", .op = MOP_MUL_ELEMS },
1242 { .phrase = "&/", .op = MOP_DIV_ELEMS },
1244 size_t n_syntax = sizeof syntax / sizeof *syntax;
1246 return matrix_parse_binary_operator (s, matrix_parse_exp, syntax, n_syntax);
1249 static struct matrix_expr *
1250 matrix_parse_add_sub (struct matrix_state *s)
1252 struct matrix_expr *lhs = matrix_parse_mul_div (s);
1259 if (lex_match (s->lexer, T_PLUS))
1261 else if (lex_match (s->lexer, T_DASH))
1263 else if (lex_token (s->lexer) == T_NEG_NUM)
1268 struct matrix_expr *rhs = matrix_parse_mul_div (s);
1271 matrix_expr_destroy (lhs);
1274 lhs = matrix_expr_create_2 (op, lhs, rhs);
1278 static struct matrix_expr *
1279 matrix_parse_relational (struct matrix_state *s)
1281 static const struct matrix_operator_syntax syntax[] = {
1282 { .token = T_GT, .op = MOP_GT },
1283 { .token = T_GE, .op = MOP_GE },
1284 { .token = T_LT, .op = MOP_LT },
1285 { .token = T_LE, .op = MOP_LE },
1286 { .token = T_EQUALS, .op = MOP_EQ },
1287 { .token = T_EQ, .op = MOP_EQ },
1288 { .token = T_NE, .op = MOP_NE },
1290 size_t n_syntax = sizeof syntax / sizeof *syntax;
1292 return matrix_parse_binary_operator (s, matrix_parse_add_sub,
1296 static struct matrix_expr *
1297 matrix_parse_not (struct matrix_state *s)
1299 int start_ofs = lex_ofs (s->lexer);
1300 if (lex_match (s->lexer, T_NOT))
1302 struct matrix_expr *sub = matrix_parse_not (s);
1306 struct matrix_expr *e = matrix_expr_create_1 (MOP_NOT, sub);
1307 matrix_expr_add_location (s, start_ofs, e);
1308 e->location->start = lex_ofs_start_point (s->lexer, start_ofs);
1312 return matrix_parse_relational (s);
1315 static struct matrix_expr *
1316 matrix_parse_and (struct matrix_state *s)
1318 static const struct matrix_operator_syntax op = {
1319 .token = T_AND, .op = MOP_AND
1322 return matrix_parse_binary_operator (s, matrix_parse_not, &op, 1);
1325 static struct matrix_expr *
1326 matrix_expr_parse__ (struct matrix_state *s)
1328 static const struct matrix_operator_syntax syntax[] = {
1329 { .token = T_OR, .op = MOP_OR },
1330 { .id = "XOR", .op = MOP_XOR },
1332 size_t n_syntax = sizeof syntax / sizeof *syntax;
1334 return matrix_parse_binary_operator (s, matrix_parse_and, syntax, n_syntax);
1337 static struct matrix_expr *
1338 matrix_expr_parse (struct matrix_state *s)
1340 int start_ofs = lex_ofs (s->lexer);
1341 struct matrix_expr *e = matrix_expr_parse__ (s);
1342 matrix_expr_add_location (s, start_ofs, e);
1346 /* Matrix expression evaluation. */
1348 /* Iterates over all the elements in matrix M, setting Y and X to the row and
1349 column indexes, respectively, and pointing D to the entry at each
1351 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
1352 for (size_t Y = 0; Y < (M)->size1; Y++) \
1353 for (size_t X = 0; X < (M)->size2; X++) \
1354 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
1357 is_vector (const gsl_matrix *m)
1359 return m->size1 <= 1 || m->size2 <= 1;
1363 to_vector (gsl_matrix *m)
1365 return (m->size1 == 1
1366 ? gsl_matrix_row (m, 0).vector
1367 : gsl_matrix_column (m, 0).vector);
1371 matrix_eval_ABS (double d)
1377 matrix_eval_ALL (gsl_matrix *m)
1379 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1386 matrix_eval_ANY (gsl_matrix *m)
1388 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1395 matrix_eval_ARSIN (double d)
1401 matrix_eval_ARTAN (double d)
1407 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
1411 for (size_t i = 0; i < n; i++)
1416 gsl_matrix *block = gsl_matrix_calloc (r, c);
1418 for (size_t i = 0; i < n; i++)
1420 for (size_t y = 0; y < m[i]->size1; y++)
1421 for (size_t x = 0; x < m[i]->size2; x++)
1422 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
1430 matrix_eval_CHOL (gsl_matrix *m, const struct matrix_expr *e)
1432 if (!gsl_linalg_cholesky_decomp1 (m))
1434 for (size_t y = 0; y < m->size1; y++)
1435 for (size_t x = y + 1; x < m->size2; x++)
1436 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
1438 for (size_t y = 0; y < m->size1; y++)
1439 for (size_t x = 0; x < y; x++)
1440 gsl_matrix_set (m, y, x, 0);
1445 msg_at (SE, e->subs[0]->location,
1446 _("Input to CHOL function is not positive-definite."));
1452 matrix_eval_col_extremum (gsl_matrix *m, bool min)
1457 return gsl_matrix_alloc (1, 0);
1459 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
1460 for (size_t x = 0; x < m->size2; x++)
1462 double ext = gsl_matrix_get (m, 0, x);
1463 for (size_t y = 1; y < m->size1; y++)
1465 double value = gsl_matrix_get (m, y, x);
1466 if (min ? value < ext : value > ext)
1469 gsl_matrix_set (cext, 0, x, ext);
1475 matrix_eval_CMAX (gsl_matrix *m)
1477 return matrix_eval_col_extremum (m, false);
1481 matrix_eval_CMIN (gsl_matrix *m)
1483 return matrix_eval_col_extremum (m, true);
1487 matrix_eval_COS (double d)
1493 matrix_eval_col_sum (gsl_matrix *m, bool square)
1498 return gsl_matrix_alloc (1, 0);
1500 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
1501 for (size_t x = 0; x < m->size2; x++)
1504 for (size_t y = 0; y < m->size1; y++)
1506 double d = gsl_matrix_get (m, y, x);
1507 sum += square ? pow2 (d) : d;
1509 gsl_matrix_set (result, 0, x, sum);
1515 matrix_eval_CSSQ (gsl_matrix *m)
1517 return matrix_eval_col_sum (m, true);
1521 matrix_eval_CSUM (gsl_matrix *m)
1523 return matrix_eval_col_sum (m, false);
1527 compare_double_3way (const void *a_, const void *b_)
1529 const double *a = a_;
1530 const double *b = b_;
1531 return *a < *b ? -1 : *a > *b;
1535 matrix_eval_DESIGN (gsl_matrix *m, const struct matrix_expr *e)
1537 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
1538 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
1539 gsl_matrix_transpose_memcpy (&m2, m);
1541 for (size_t y = 0; y < m2.size1; y++)
1542 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
1544 size_t *n = xcalloc (m2.size1, sizeof *n);
1546 for (size_t i = 0; i < m2.size1; i++)
1548 double *row = tmp + m2.size2 * i;
1549 for (size_t j = 0; j < m2.size2; )
1552 for (k = j + 1; k < m2.size2; k++)
1553 if (row[j] != row[k])
1555 row[n[i]++] = row[j];
1560 msg_at (MW, e->subs[0]->location,
1561 _("Column %zu in DESIGN argument has constant value."), i + 1);
1566 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
1568 for (size_t i = 0; i < m->size2; i++)
1573 const double *unique = tmp + m2.size2 * i;
1574 for (size_t j = 0; j < n[i]; j++, x++)
1576 double value = unique[j];
1577 for (size_t y = 0; y < m->size1; y++)
1578 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
1589 matrix_eval_DET (gsl_matrix *m)
1591 gsl_permutation *p = gsl_permutation_alloc (m->size1);
1593 gsl_linalg_LU_decomp (m, p, &signum);
1594 gsl_permutation_free (p);
1595 return gsl_linalg_LU_det (m, signum);
1599 matrix_eval_DIAG (gsl_matrix *m)
1601 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
1602 for (size_t i = 0; i < diag->size1; i++)
1603 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
1608 is_symmetric (const gsl_matrix *m)
1610 if (m->size1 != m->size2)
1613 for (size_t y = 0; y < m->size1; y++)
1614 for (size_t x = 0; x < y; x++)
1615 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
1622 compare_double_desc (const void *a_, const void *b_)
1624 const double *a = a_;
1625 const double *b = b_;
1626 return *a > *b ? -1 : *a < *b;
1630 matrix_eval_EVAL (gsl_matrix *m, const struct matrix_expr *e)
1632 if (!is_symmetric (m))
1634 msg_at (SE, e->subs[0]->location,
1635 _("Argument of EVAL must be symmetric."));
1639 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
1640 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
1641 gsl_vector v_eval = to_vector (eval);
1642 gsl_eigen_symm (m, &v_eval, w);
1643 gsl_eigen_symm_free (w);
1645 assert (v_eval.stride == 1);
1646 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
1652 matrix_eval_EXP (double d)
1657 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
1660 Charl Linssen <charl@itfromb.it>
1664 matrix_eval_GINV (gsl_matrix *A)
1666 size_t n = A->size1;
1667 size_t m = A->size2;
1669 gsl_matrix *tmp_mat = NULL;
1672 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
1673 tmp_mat = gsl_matrix_alloc (m, n);
1674 gsl_matrix_transpose_memcpy (tmp_mat, A);
1682 gsl_matrix *V = gsl_matrix_alloc (m, m);
1683 gsl_vector *u = gsl_vector_alloc (m);
1685 gsl_vector *tmp_vec = gsl_vector_alloc (m);
1686 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
1687 gsl_vector_free (tmp_vec);
1690 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
1691 gsl_matrix_set_zero (Sigma_pinv);
1692 double cutoff = 1e-15 * gsl_vector_max (u);
1694 for (size_t i = 0; i < m; ++i)
1696 double x = gsl_vector_get (u, i);
1697 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
1700 /* libgsl SVD yields "thin" SVD. Pad to full matrix by adding zeros. */
1701 gsl_matrix *U = gsl_matrix_calloc (n, n);
1702 for (size_t i = 0; i < n; i++)
1703 for (size_t j = 0; j < m; j++)
1704 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
1706 /* Two dot products to obtain pseudoinverse. */
1707 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
1708 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
1713 A_pinv = gsl_matrix_alloc (n, m);
1714 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
1718 A_pinv = gsl_matrix_alloc (m, n);
1719 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
1722 gsl_matrix_free (tmp_mat);
1723 gsl_matrix_free (tmp_mat2);
1724 gsl_matrix_free (U);
1725 gsl_matrix_free (Sigma_pinv);
1726 gsl_vector_free (u);
1727 gsl_matrix_free (V);
1739 grade_compare_3way (const void *a_, const void *b_)
1741 const struct grade *a = a_;
1742 const struct grade *b = b_;
1744 return (a->value < b->value ? -1
1745 : a->value > b->value ? 1
1753 matrix_eval_GRADE (gsl_matrix *m)
1755 size_t n = m->size1 * m->size2;
1756 struct grade *grades = xmalloc (n * sizeof *grades);
1759 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1760 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
1761 qsort (grades, n, sizeof *grades, grade_compare_3way);
1763 for (size_t i = 0; i < n; i++)
1764 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
1772 dot (gsl_vector *a, gsl_vector *b)
1774 double result = 0.0;
1775 for (size_t i = 0; i < a->size; i++)
1776 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
1781 norm2 (gsl_vector *v)
1783 double result = 0.0;
1784 for (size_t i = 0; i < v->size; i++)
1785 result += pow2 (gsl_vector_get (v, i));
1790 norm (gsl_vector *v)
1792 return sqrt (norm2 (v));
1796 matrix_eval_GSCH (gsl_matrix *v, const struct matrix_expr *e)
1798 if (v->size2 < v->size1)
1800 msg_at (SE, e->subs[0]->location,
1801 _("GSCH requires its argument to have at least as many columns "
1802 "as rows, but it has dimensions %zu×%zu."),
1803 v->size1, v->size2);
1806 if (!v->size1 || !v->size2)
1809 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
1811 for (size_t vx = 0; vx < v->size2; vx++)
1813 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
1814 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
1816 gsl_vector_memcpy (&u_i, &v_i);
1817 for (size_t j = 0; j < ux; j++)
1819 gsl_vector u_j = gsl_matrix_column (u, j).vector;
1820 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
1821 for (size_t k = 0; k < u_i.size; k++)
1822 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
1823 - scale * gsl_vector_get (&u_j, k)));
1826 double len = norm (&u_i);
1829 gsl_vector_scale (&u_i, 1.0 / len);
1830 if (++ux >= v->size1)
1837 msg_at (SE, e->subs[0]->location,
1838 _("%zu×%zu argument to GSCH contains only "
1839 "%zu linearly independent columns."),
1840 v->size1, v->size2, ux);
1841 gsl_matrix_free (u);
1845 u->size2 = v->size1;
1850 matrix_eval_IDENT (double s1, double s2)
1852 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
1853 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1858 /* Inverts X, storing the inverse into INVERSE. As a side effect, replaces X
1859 by its LU decomposition. */
1861 invert_matrix (gsl_matrix *x, gsl_matrix *inverse)
1863 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1865 gsl_linalg_LU_decomp (x, p, &signum);
1866 gsl_linalg_LU_invert (x, p, inverse);
1867 gsl_permutation_free (p);
1871 matrix_eval_INV (gsl_matrix *src)
1873 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
1874 invert_matrix (src, dst);
1879 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1881 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1882 a->size2 * b->size2);
1884 for (size_t ar = 0; ar < a->size1; ar++)
1885 for (size_t br = 0; br < b->size1; br++, y++)
1888 for (size_t ac = 0; ac < a->size2; ac++)
1889 for (size_t bc = 0; bc < b->size2; bc++, x++)
1891 double av = gsl_matrix_get (a, ar, ac);
1892 double bv = gsl_matrix_get (b, br, bc);
1893 gsl_matrix_set (k, y, x, av * bv);
1900 matrix_eval_LG10 (double d)
1906 matrix_eval_LN (double d)
1912 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1914 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1917 for (size_t i = 1; i <= n * n; i++)
1919 gsl_matrix_set (m, y, x, i);
1921 size_t y1 = !y ? n - 1 : y - 1;
1922 size_t x1 = x + 1 >= n ? 0 : x + 1;
1923 if (gsl_matrix_get (m, y1, x1) == 0)
1929 y = y + 1 >= n ? 0 : y + 1;
1934 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1936 double a = gsl_matrix_get (m, y1, x1);
1937 double b = gsl_matrix_get (m, y2, x2);
1938 gsl_matrix_set (m, y1, x1, b);
1939 gsl_matrix_set (m, y2, x2, a);
1943 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1947 /* A. Umar, "On the Construction of Even Order Magic Squares",
1948 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1950 for (size_t i = 1; i <= n * n / 2; i++)
1952 gsl_matrix_set (m, y, x, i);
1962 for (size_t i = n * n; i > n * n / 2; i--)
1964 gsl_matrix_set (m, y, x, i);
1972 for (size_t y = 0; y < n; y++)
1973 for (size_t x = 0; x < n / 2; x++)
1975 unsigned int d = gsl_matrix_get (m, y, x);
1976 if (d % 2 != (y < n / 2))
1977 magic_exchange (m, y, x, y, n - x - 1);
1982 size_t x1 = n / 2 - 1;
1984 magic_exchange (m, y1, x1, y2, x1);
1985 magic_exchange (m, y1, x2, y2, x2);
1989 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1991 /* A. Umar, "On the Construction of Even Order Magic Squares",
1992 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1996 for (size_t i = 1; ; i++)
1998 gsl_matrix_set (m, y, x, i);
1999 if (++y == n / 2 - 1)
2011 for (size_t i = n * n; ; i--)
2013 gsl_matrix_set (m, y, x, i);
2014 if (++y == n / 2 - 1)
2023 for (size_t y = 0; y < n; y++)
2024 if (y != n / 2 - 1 && y != n / 2)
2025 for (size_t x = 0; x < n / 2; x++)
2027 unsigned int d = gsl_matrix_get (m, y, x);
2028 if (d % 2 != (y < n / 2))
2029 magic_exchange (m, y, x, y, n - x - 1);
2032 size_t a0 = (n * n - 2 * n) / 2 + 1;
2033 for (size_t i = 0; i < n / 2; i++)
2036 gsl_matrix_set (m, n / 2, i, a);
2037 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
2039 for (size_t i = 0; i < n / 2; i++)
2041 size_t a = a0 + i + n / 2;
2042 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
2043 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
2045 for (size_t x = 1; x < n / 2; x += 2)
2046 magic_exchange (m, n / 2, x, n / 2 - 1, x);
2047 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
2048 magic_exchange (m, n / 2, x, n / 2 - 1, x);
2049 size_t x1 = n / 2 - 2;
2050 size_t x2 = n / 2 + 1;
2051 size_t y1 = n / 2 - 2;
2052 size_t y2 = n / 2 + 1;
2053 magic_exchange (m, y1, x1, y2, x1);
2054 magic_exchange (m, y1, x2, y2, x2);
2058 matrix_eval_MAGIC (double n_)
2062 gsl_matrix *m = gsl_matrix_calloc (n, n);
2064 matrix_eval_MAGIC_odd (m, n);
2066 matrix_eval_MAGIC_singly_even (m, n);
2068 matrix_eval_MAGIC_doubly_even (m, n);
2073 matrix_eval_MAKE (double r, double c, double value)
2075 gsl_matrix *m = gsl_matrix_alloc (r, c);
2076 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2082 matrix_eval_MDIAG (gsl_vector *v)
2084 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
2085 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
2086 gsl_vector_memcpy (&diagonal, v);
2091 matrix_eval_MMAX (gsl_matrix *m)
2093 return gsl_matrix_max (m);
2097 matrix_eval_MMIN (gsl_matrix *m)
2099 return gsl_matrix_min (m);
2103 matrix_eval_MOD (gsl_matrix *m, double divisor)
2105 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2106 *d = fmod (*d, divisor);
2111 matrix_eval_MSSQ (gsl_matrix *m)
2114 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2120 matrix_eval_MSUM (gsl_matrix *m)
2123 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2129 matrix_eval_NCOL (gsl_matrix *m)
2135 matrix_eval_NROW (gsl_matrix *m)
2141 matrix_eval_RANK (gsl_matrix *m)
2143 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
2144 gsl_linalg_QR_decomp (m, tau);
2145 gsl_vector_free (tau);
2147 return gsl_linalg_QRPT_rank (m, -1);
2151 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_,
2152 const struct matrix_expr *e)
2154 bool r_ok = r_ >= 0 && r_ < SIZE_MAX;
2155 bool c_ok = c_ >= 0 && c_ < SIZE_MAX;
2159 !r_ok ? e->subs[1]->location : e->subs[2]->location,
2160 _("Arguments 2 and 3 to RESHAPE must be integers."));
2165 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
2167 struct msg_location *loc = msg_location_dup (e->subs[1]->location);
2168 loc->end = e->subs[2]->location->end;
2169 msg_at (SE, loc, _("Product of RESHAPE size arguments (%zu×%zu = %zu) "
2170 "differs from product of matrix dimensions "
2171 "(%zu×%zu = %zu)."),
2173 m->size1, m->size2, m->size1 * m->size2);
2174 msg_location_destroy (loc);
2178 gsl_matrix *dst = gsl_matrix_alloc (r, c);
2181 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
2183 gsl_matrix_set (dst, y1, x1, *d);
2194 matrix_eval_row_extremum (gsl_matrix *m, bool min)
2199 return gsl_matrix_alloc (0, 1);
2201 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
2202 for (size_t y = 0; y < m->size1; y++)
2204 double ext = gsl_matrix_get (m, y, 0);
2205 for (size_t x = 1; x < m->size2; x++)
2207 double value = gsl_matrix_get (m, y, x);
2208 if (min ? value < ext : value > ext)
2211 gsl_matrix_set (rext, y, 0, ext);
2217 matrix_eval_RMAX (gsl_matrix *m)
2219 return matrix_eval_row_extremum (m, false);
2223 matrix_eval_RMIN (gsl_matrix *m)
2225 return matrix_eval_row_extremum (m, true);
2229 matrix_eval_RND (double d)
2241 rank_compare_3way (const void *a_, const void *b_)
2243 const struct rank *a = a_;
2244 const struct rank *b = b_;
2246 return a->value < b->value ? -1 : a->value > b->value;
2250 matrix_eval_RNKORDER (gsl_matrix *m)
2252 size_t n = m->size1 * m->size2;
2253 struct rank *ranks = xmalloc (n * sizeof *ranks);
2255 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2256 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
2257 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
2259 for (size_t i = 0; i < n; )
2262 for (j = i + 1; j < n; j++)
2263 if (ranks[i].value != ranks[j].value)
2266 double rank = (i + j + 1.0) / 2.0;
2267 for (size_t k = i; k < j; k++)
2268 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
2279 matrix_eval_row_sum (gsl_matrix *m, bool square)
2284 return gsl_matrix_alloc (0, 1);
2286 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
2287 for (size_t y = 0; y < m->size1; y++)
2290 for (size_t x = 0; x < m->size2; x++)
2292 double d = gsl_matrix_get (m, y, x);
2293 sum += square ? pow2 (d) : d;
2295 gsl_matrix_set (result, y, 0, sum);
2301 matrix_eval_RSSQ (gsl_matrix *m)
2303 return matrix_eval_row_sum (m, true);
2307 matrix_eval_RSUM (gsl_matrix *m)
2309 return matrix_eval_row_sum (m, false);
2313 matrix_eval_SIN (double d)
2319 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2, const struct matrix_expr *e)
2321 if (m1->size1 != m2->size1)
2323 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
2324 loc->end = e->subs[1]->location->end;
2326 msg_at (SE, e->location,
2327 _("SOLVE arguments must have the same number of rows."));
2328 msg_at (SN, e->subs[0]->location,
2329 _("Argument 1 has dimensions %zu×%zu."), m1->size1, m1->size2);
2330 msg_at (SN, e->subs[1]->location,
2331 _("Argument 2 has dimensions %zu×%zu."), m2->size1, m2->size2);
2333 msg_location_destroy (loc);
2337 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
2338 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
2340 gsl_linalg_LU_decomp (m1, p, &signum);
2341 for (size_t i = 0; i < m2->size2; i++)
2343 gsl_vector bi = gsl_matrix_column (m2, i).vector;
2344 gsl_vector xi = gsl_matrix_column (x, i).vector;
2345 gsl_linalg_LU_solve (m1, p, &bi, &xi);
2347 gsl_permutation_free (p);
2352 matrix_eval_SQRT (double d)
2358 matrix_eval_SSCP (gsl_matrix *m)
2360 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
2361 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
2366 matrix_eval_SVAL (gsl_matrix *m)
2368 gsl_matrix *tmp_mat = NULL;
2369 if (m->size2 > m->size1)
2371 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
2372 gsl_matrix_transpose_memcpy (tmp_mat, m);
2377 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
2378 gsl_vector *S = gsl_vector_alloc (m->size2);
2379 gsl_vector *work = gsl_vector_alloc (m->size2);
2380 gsl_linalg_SV_decomp (m, V, S, work);
2382 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
2383 for (size_t i = 0; i < m->size2; i++)
2384 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
2386 gsl_matrix_free (V);
2387 gsl_vector_free (S);
2388 gsl_vector_free (work);
2389 gsl_matrix_free (tmp_mat);
2395 matrix_eval_SWEEP (gsl_matrix *m, double d, const struct matrix_expr *e)
2397 if (d < 1 || d > SIZE_MAX)
2399 msg_at (SE, e->subs[1]->location,
2400 _("Scalar argument to SWEEP must be integer."));
2404 if (k >= MIN (m->size1, m->size2))
2406 msg_at (SE, e->subs[1]->location,
2407 _("Scalar argument to SWEEP must be integer less than or "
2408 "equal to the smaller of the matrix argument's rows and "
2413 double m_kk = gsl_matrix_get (m, k, k);
2414 if (fabs (m_kk) > 1e-19)
2416 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
2417 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
2419 double m_ij = gsl_matrix_get (m, i, j);
2420 double m_ik = gsl_matrix_get (m, i, k);
2421 double m_kj = gsl_matrix_get (m, k, j);
2422 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
2431 for (size_t i = 0; i < m->size1; i++)
2433 gsl_matrix_set (m, i, k, 0);
2434 gsl_matrix_set (m, k, i, 0);
2441 matrix_eval_TRACE (gsl_matrix *m)
2444 size_t n = MIN (m->size1, m->size2);
2445 for (size_t i = 0; i < n; i++)
2446 sum += gsl_matrix_get (m, i, i);
2451 matrix_eval_T (gsl_matrix *m)
2453 return matrix_eval_TRANSPOS (m);
2457 matrix_eval_TRANSPOS (gsl_matrix *m)
2459 if (m->size1 == m->size2)
2461 gsl_matrix_transpose (m);
2466 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
2467 gsl_matrix_transpose_memcpy (t, m);
2473 matrix_eval_TRUNC (double d)
2479 matrix_eval_UNIFORM (double r_, double c_, const struct matrix_expr *e)
2483 if (size_overflow_p (xtimes (r, xmax (c, 1))))
2485 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
2486 loc->end = e->subs[1]->location->end;
2489 _("Product of arguments to UNIFORM exceeds memory size."));
2491 msg_location_destroy (loc);
2495 gsl_matrix *m = gsl_matrix_alloc (r, c);
2496 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2497 *d = gsl_ran_flat (get_rng (), 0, 1);
2502 matrix_eval_PDF_BETA (double x, double a, double b)
2504 return gsl_ran_beta_pdf (x, a, b);
2508 matrix_eval_CDF_BETA (double x, double a, double b)
2510 return gsl_cdf_beta_P (x, a, b);
2514 matrix_eval_IDF_BETA (double P, double a, double b)
2516 return gsl_cdf_beta_Pinv (P, a, b);
2520 matrix_eval_RV_BETA (double a, double b)
2522 return gsl_ran_beta (get_rng (), a, b);
2526 matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
2528 return ncdf_beta (x, a, b, lambda);
2532 matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
2534 return npdf_beta (x, a, b, lambda);
2538 matrix_eval_CDF_BVNOR (double x0, double x1, double r)
2540 return cdf_bvnor (x0, x1, r);
2544 matrix_eval_PDF_BVNOR (double x0, double x1, double r)
2546 return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
2550 matrix_eval_CDF_CAUCHY (double x, double a, double b)
2552 return gsl_cdf_cauchy_P ((x - a) / b, 1);
2556 matrix_eval_IDF_CAUCHY (double P, double a, double b)
2558 return a + b * gsl_cdf_cauchy_Pinv (P, 1);
2562 matrix_eval_PDF_CAUCHY (double x, double a, double b)
2564 return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
2568 matrix_eval_RV_CAUCHY (double a, double b)
2570 return a + b * gsl_ran_cauchy (get_rng (), 1);
2574 matrix_eval_CDF_CHISQ (double x, double df)
2576 return gsl_cdf_chisq_P (x, df);
2580 matrix_eval_CHICDF (double x, double df)
2582 return matrix_eval_CDF_CHISQ (x, df);
2586 matrix_eval_IDF_CHISQ (double P, double df)
2588 return gsl_cdf_chisq_Pinv (P, df);
2592 matrix_eval_PDF_CHISQ (double x, double df)
2594 return gsl_ran_chisq_pdf (x, df);
2598 matrix_eval_RV_CHISQ (double df)
2600 return gsl_ran_chisq (get_rng (), df);
2604 matrix_eval_SIG_CHISQ (double x, double df)
2606 return gsl_cdf_chisq_Q (x, df);
2610 matrix_eval_CDF_EXP (double x, double a)
2612 return gsl_cdf_exponential_P (x, 1. / a);
2616 matrix_eval_IDF_EXP (double P, double a)
2618 return gsl_cdf_exponential_Pinv (P, 1. / a);
2622 matrix_eval_PDF_EXP (double x, double a)
2624 return gsl_ran_exponential_pdf (x, 1. / a);
2628 matrix_eval_RV_EXP (double a)
2630 return gsl_ran_exponential (get_rng (), 1. / a);
2634 matrix_eval_PDF_XPOWER (double x, double a, double b)
2636 return gsl_ran_exppow_pdf (x, a, b);
2640 matrix_eval_RV_XPOWER (double a, double b)
2642 return gsl_ran_exppow (get_rng (), a, b);
2646 matrix_eval_CDF_F (double x, double df1, double df2)
2648 return gsl_cdf_fdist_P (x, df1, df2);
2652 matrix_eval_FCDF (double x, double df1, double df2)
2654 return matrix_eval_CDF_F (x, df1, df2);
2658 matrix_eval_IDF_F (double P, double df1, double df2)
2660 return idf_fdist (P, df1, df2);
2664 matrix_eval_RV_F (double df1, double df2)
2666 return gsl_ran_fdist (get_rng (), df1, df2);
2670 matrix_eval_PDF_F (double x, double df1, double df2)
2672 return gsl_ran_fdist_pdf (x, df1, df2);
2676 matrix_eval_SIG_F (double x, double df1, double df2)
2678 return gsl_cdf_fdist_Q (x, df1, df2);
2682 matrix_eval_CDF_GAMMA (double x, double a, double b)
2684 return gsl_cdf_gamma_P (x, a, 1. / b);
2688 matrix_eval_IDF_GAMMA (double P, double a, double b)
2690 return gsl_cdf_gamma_Pinv (P, a, 1. / b);
2694 matrix_eval_PDF_GAMMA (double x, double a, double b)
2696 return gsl_ran_gamma_pdf (x, a, 1. / b);
2700 matrix_eval_RV_GAMMA (double a, double b)
2702 return gsl_ran_gamma (get_rng (), a, 1. / b);
2706 matrix_eval_PDF_LANDAU (double x)
2708 return gsl_ran_landau_pdf (x);
2712 matrix_eval_RV_LANDAU (void)
2714 return gsl_ran_landau (get_rng ());
2718 matrix_eval_CDF_LAPLACE (double x, double a, double b)
2720 return gsl_cdf_laplace_P ((x - a) / b, 1);
2724 matrix_eval_IDF_LAPLACE (double P, double a, double b)
2726 return a + b * gsl_cdf_laplace_Pinv (P, 1);
2730 matrix_eval_PDF_LAPLACE (double x, double a, double b)
2732 return gsl_ran_laplace_pdf ((x - a) / b, 1);
2736 matrix_eval_RV_LAPLACE (double a, double b)
2738 return a + b * gsl_ran_laplace (get_rng (), 1);
2742 matrix_eval_RV_LEVY (double c, double alpha)
2744 return gsl_ran_levy (get_rng (), c, alpha);
2748 matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
2750 return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
2754 matrix_eval_CDF_LOGISTIC (double x, double a, double b)
2756 return gsl_cdf_logistic_P ((x - a) / b, 1);
2760 matrix_eval_IDF_LOGISTIC (double P, double a, double b)
2762 return a + b * gsl_cdf_logistic_Pinv (P, 1);
2766 matrix_eval_PDF_LOGISTIC (double x, double a, double b)
2768 return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
2772 matrix_eval_RV_LOGISTIC (double a, double b)
2774 return a + b * gsl_ran_logistic (get_rng (), 1);
2778 matrix_eval_CDF_LNORMAL (double x, double m, double s)
2780 return gsl_cdf_lognormal_P (x, log (m), s);
2784 matrix_eval_IDF_LNORMAL (double P, double m, double s)
2786 return gsl_cdf_lognormal_Pinv (P, log (m), s);;
2790 matrix_eval_PDF_LNORMAL (double x, double m, double s)
2792 return gsl_ran_lognormal_pdf (x, log (m), s);
2796 matrix_eval_RV_LNORMAL (double m, double s)
2798 return gsl_ran_lognormal (get_rng (), log (m), s);
2802 matrix_eval_CDF_NORMAL (double x, double u, double s)
2804 return gsl_cdf_gaussian_P (x - u, s);
2808 matrix_eval_IDF_NORMAL (double P, double u, double s)
2810 return u + gsl_cdf_gaussian_Pinv (P, s);
2814 matrix_eval_PDF_NORMAL (double x, double u, double s)
2816 return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
2820 matrix_eval_RV_NORMAL (double u, double s)
2822 return u + gsl_ran_gaussian (get_rng (), s);
2826 matrix_eval_CDFNORM (double x)
2828 return gsl_cdf_ugaussian_P (x);
2832 matrix_eval_PROBIT (double P)
2834 return gsl_cdf_ugaussian_Pinv (P);
2838 matrix_eval_NORMAL (double s)
2840 return gsl_ran_gaussian (get_rng (), s);
2844 matrix_eval_PDF_NTAIL (double x, double a, double sigma)
2846 return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
2850 matrix_eval_RV_NTAIL (double a, double sigma)
2852 return gsl_ran_gaussian_tail (get_rng (), a, sigma);
2856 matrix_eval_CDF_PARETO (double x, double a, double b)
2858 return gsl_cdf_pareto_P (x, b, a);
2862 matrix_eval_IDF_PARETO (double P, double a, double b)
2864 return gsl_cdf_pareto_Pinv (P, b, a);
2868 matrix_eval_PDF_PARETO (double x, double a, double b)
2870 return gsl_ran_pareto_pdf (x, b, a);
2874 matrix_eval_RV_PARETO (double a, double b)
2876 return gsl_ran_pareto (get_rng (), b, a);
2880 matrix_eval_CDF_RAYLEIGH (double x, double sigma)
2882 return gsl_cdf_rayleigh_P (x, sigma);
2886 matrix_eval_IDF_RAYLEIGH (double P, double sigma)
2888 return gsl_cdf_rayleigh_Pinv (P, sigma);
2892 matrix_eval_PDF_RAYLEIGH (double x, double sigma)
2894 return gsl_ran_rayleigh_pdf (x, sigma);
2898 matrix_eval_RV_RAYLEIGH (double sigma)
2900 return gsl_ran_rayleigh (get_rng (), sigma);
2904 matrix_eval_PDF_RTAIL (double x, double a, double sigma)
2906 return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
2910 matrix_eval_RV_RTAIL (double a, double sigma)
2912 return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
2916 matrix_eval_CDF_T (double x, double df)
2918 return gsl_cdf_tdist_P (x, df);
2922 matrix_eval_TCDF (double x, double df)
2924 return matrix_eval_CDF_T (x, df);
2928 matrix_eval_IDF_T (double P, double df)
2930 return gsl_cdf_tdist_Pinv (P, df);
2934 matrix_eval_PDF_T (double x, double df)
2936 return gsl_ran_tdist_pdf (x, df);
2940 matrix_eval_RV_T (double df)
2942 return gsl_ran_tdist (get_rng (), df);
2946 matrix_eval_CDF_T1G (double x, double a, double b)
2948 return gsl_cdf_gumbel1_P (x, a, b);
2952 matrix_eval_IDF_T1G (double P, double a, double b)
2954 return gsl_cdf_gumbel1_Pinv (P, a, b);
2958 matrix_eval_PDF_T1G (double x, double a, double b)
2960 return gsl_ran_gumbel1_pdf (x, a, b);
2964 matrix_eval_RV_T1G (double a, double b)
2966 return gsl_ran_gumbel1 (get_rng (), a, b);
2970 matrix_eval_CDF_T2G (double x, double a, double b)
2972 return gsl_cdf_gumbel1_P (x, a, b);
2976 matrix_eval_IDF_T2G (double P, double a, double b)
2978 return gsl_cdf_gumbel1_Pinv (P, a, b);
2982 matrix_eval_PDF_T2G (double x, double a, double b)
2984 return gsl_ran_gumbel1_pdf (x, a, b);
2988 matrix_eval_RV_T2G (double a, double b)
2990 return gsl_ran_gumbel1 (get_rng (), a, b);
2994 matrix_eval_CDF_UNIFORM (double x, double a, double b)
2996 return gsl_cdf_flat_P (x, a, b);
3000 matrix_eval_IDF_UNIFORM (double P, double a, double b)
3002 return gsl_cdf_flat_Pinv (P, a, b);
3006 matrix_eval_PDF_UNIFORM (double x, double a, double b)
3008 return gsl_ran_flat_pdf (x, a, b);
3012 matrix_eval_RV_UNIFORM (double a, double b)
3014 return gsl_ran_flat (get_rng (), a, b);
3018 matrix_eval_CDF_WEIBULL (double x, double a, double b)
3020 return gsl_cdf_weibull_P (x, a, b);
3024 matrix_eval_IDF_WEIBULL (double P, double a, double b)
3026 return gsl_cdf_weibull_Pinv (P, a, b);
3030 matrix_eval_PDF_WEIBULL (double x, double a, double b)
3032 return gsl_ran_weibull_pdf (x, a, b);
3036 matrix_eval_RV_WEIBULL (double a, double b)
3038 return gsl_ran_weibull (get_rng (), a, b);
3042 matrix_eval_CDF_BERNOULLI (double k, double p)
3044 return k ? 1 : 1 - p;
3048 matrix_eval_PDF_BERNOULLI (double k, double p)
3050 return gsl_ran_bernoulli_pdf (k, p);
3054 matrix_eval_RV_BERNOULLI (double p)
3056 return gsl_ran_bernoulli (get_rng (), p);
3060 matrix_eval_CDF_BINOM (double k, double n, double p)
3062 return gsl_cdf_binomial_P (k, p, n);
3066 matrix_eval_PDF_BINOM (double k, double n, double p)
3068 return gsl_ran_binomial_pdf (k, p, n);
3072 matrix_eval_RV_BINOM (double n, double p)
3074 return gsl_ran_binomial (get_rng (), p, n);
3078 matrix_eval_CDF_GEOM (double k, double p)
3080 return gsl_cdf_geometric_P (k, p);
3084 matrix_eval_PDF_GEOM (double k, double p)
3086 return gsl_ran_geometric_pdf (k, p);
3090 matrix_eval_RV_GEOM (double p)
3092 return gsl_ran_geometric (get_rng (), p);
3096 matrix_eval_CDF_HYPER (double k, double a, double b, double c)
3098 return gsl_cdf_hypergeometric_P (k, c, a - c, b);
3102 matrix_eval_PDF_HYPER (double k, double a, double b, double c)
3104 return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
3108 matrix_eval_RV_HYPER (double a, double b, double c)
3110 return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
3114 matrix_eval_PDF_LOG (double k, double p)
3116 return gsl_ran_logarithmic_pdf (k, p);
3120 matrix_eval_RV_LOG (double p)
3122 return gsl_ran_logarithmic (get_rng (), p);
3126 matrix_eval_CDF_NEGBIN (double k, double n, double p)
3128 return gsl_cdf_negative_binomial_P (k, p, n);
3132 matrix_eval_PDF_NEGBIN (double k, double n, double p)
3134 return gsl_ran_negative_binomial_pdf (k, p, n);
3138 matrix_eval_RV_NEGBIN (double n, double p)
3140 return gsl_ran_negative_binomial (get_rng (), p, n);
3144 matrix_eval_CDF_POISSON (double k, double mu)
3146 return gsl_cdf_poisson_P (k, mu);
3150 matrix_eval_PDF_POISSON (double k, double mu)
3152 return gsl_ran_poisson_pdf (k, mu);
3156 matrix_eval_RV_POISSON (double mu)
3158 return gsl_ran_poisson (get_rng (), mu);
3162 matrix_op_eval (enum matrix_op op, double a, double b)
3166 case MOP_ADD_ELEMS: return a + b;
3167 case MOP_SUB_ELEMS: return a - b;
3168 case MOP_MUL_ELEMS: return a * b;
3169 case MOP_DIV_ELEMS: return a / b;
3170 case MOP_EXP_ELEMS: return pow (a, b);
3171 case MOP_GT: return a > b;
3172 case MOP_GE: return a >= b;
3173 case MOP_LT: return a < b;
3174 case MOP_LE: return a <= b;
3175 case MOP_EQ: return a == b;
3176 case MOP_NE: return a != b;
3177 case MOP_AND: return (a > 0) && (b > 0);
3178 case MOP_OR: return (a > 0) || (b > 0);
3179 case MOP_XOR: return (a > 0) != (b > 0);
3181 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3190 case MOP_PASTE_HORZ:
3191 case MOP_PASTE_VERT:
3207 matrix_op_name (enum matrix_op op)
3211 case MOP_ADD_ELEMS: return "+";
3212 case MOP_SUB_ELEMS: return "-";
3213 case MOP_MUL_ELEMS: return "&*";
3214 case MOP_DIV_ELEMS: return "&/";
3215 case MOP_EXP_ELEMS: return "&**";
3216 case MOP_GT: return ">";
3217 case MOP_GE: return ">=";
3218 case MOP_LT: return "<";
3219 case MOP_LE: return "<=";
3220 case MOP_EQ: return "=";
3221 case MOP_NE: return "<>";
3222 case MOP_AND: return "AND";
3223 case MOP_OR: return "OR";
3224 case MOP_XOR: return "XOR";
3226 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3235 case MOP_PASTE_HORZ:
3236 case MOP_PASTE_VERT:
3252 is_scalar (const gsl_matrix *m)
3254 return m->size1 == 1 && m->size2 == 1;
3258 to_scalar (const gsl_matrix *m)
3260 assert (is_scalar (m));
3261 return gsl_matrix_get (m, 0, 0);
3265 matrix_expr_evaluate_elementwise (const struct matrix_expr *e,
3267 gsl_matrix *a, gsl_matrix *b)
3271 double be = to_scalar (b);
3272 for (size_t r = 0; r < a->size1; r++)
3273 for (size_t c = 0; c < a->size2; c++)
3275 double *ae = gsl_matrix_ptr (a, r, c);
3276 *ae = matrix_op_eval (op, *ae, be);
3280 else if (is_scalar (a))
3282 double ae = to_scalar (a);
3283 for (size_t r = 0; r < b->size1; r++)
3284 for (size_t c = 0; c < b->size2; c++)
3286 double *be = gsl_matrix_ptr (b, r, c);
3287 *be = matrix_op_eval (op, ae, *be);
3291 else if (a->size1 == b->size1 && a->size2 == b->size2)
3293 for (size_t r = 0; r < a->size1; r++)
3294 for (size_t c = 0; c < a->size2; c++)
3296 double *ae = gsl_matrix_ptr (a, r, c);
3297 double be = gsl_matrix_get (b, r, c);
3298 *ae = matrix_op_eval (op, *ae, be);
3304 msg_at (SE, matrix_expr_location (e),
3305 _("The operands of %s must have the same dimensions or one "
3306 "must be a scalar."),
3307 matrix_op_name (op));
3308 msg_at (SN, matrix_expr_location (e->subs[0]),
3309 _("The left-hand operand is a %zu×%zu matrix."),
3310 a->size1, a->size2);
3311 msg_at (SN, matrix_expr_location (e->subs[1]),
3312 _("The right-hand operand is a %zu×%zu matrix."),
3313 b->size1, b->size2);
3319 matrix_expr_evaluate_mul_mat (const struct matrix_expr *e,
3320 gsl_matrix *a, gsl_matrix *b)
3322 if (is_scalar (a) || is_scalar (b))
3323 return matrix_expr_evaluate_elementwise (e, MOP_MUL_ELEMS, a, b);
3325 if (a->size2 != b->size1)
3327 msg_at (SE, e->location,
3328 _("Matrices not conformable for multiplication."));
3329 msg_at (SN, matrix_expr_location (e->subs[0]),
3330 _("The left-hand operand is a %zu×%zu matrix."),
3331 a->size1, a->size2);
3332 msg_at (SN, matrix_expr_location (e->subs[1]),
3333 _("The right-hand operand is a %zu×%zu matrix."),
3334 b->size1, b->size2);
3338 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3339 if (a->size1 && b->size2)
3340 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3345 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3347 gsl_matrix *tmp = *a;
3353 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3356 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3357 swap_matrix (z, tmp);
3361 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3363 mul_matrix (x, *x, *x, tmp);
3367 matrix_expr_evaluate_exp_mat (const struct matrix_expr *e,
3368 gsl_matrix *x_, gsl_matrix *b)
3371 if (x->size1 != x->size2)
3373 msg_at (SE, matrix_expr_location (e->subs[0]),
3374 _("Matrix exponentation with ** requires a square matrix on "
3375 "the left-hand size, not one with dimensions %zu×%zu."),
3376 x->size1, x->size2);
3381 msg_at (SE, matrix_expr_location (e->subs[1]),
3382 _("Matrix exponentiation with ** requires a scalar on the "
3383 "right-hand side, not a matrix with dimensions %zu×%zu."),
3384 b->size1, b->size2);
3387 double bf = to_scalar (b);
3388 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3390 msg_at (SE, matrix_expr_location (e->subs[1]),
3391 _("Exponent %.1f in matrix exponentiation is non-integer "
3392 "or outside the valid range."), bf);
3397 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3399 gsl_matrix_set_identity (y);
3403 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3405 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3408 mul_matrix (&y, x, y, &t);
3409 square_matrix (&x, &t);
3412 square_matrix (&x, &t);
3414 mul_matrix (&y, x, y, &t);
3417 invert_matrix (y, x);
3418 swap_matrix (&x, &y);
3421 /* Garbage collection.
3423 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3424 a permutation of them. We are returning one of them; that one must not be
3425 destroyed. We must not destroy 'x_' because the caller owns it. */
3427 gsl_matrix_free (y_);
3429 gsl_matrix_free (t_);
3435 note_operand_size (const gsl_matrix *m, const struct matrix_expr *e)
3437 msg_at (SN, matrix_expr_location (e),
3438 _("This operand is a %zu×%zu matrix."), m->size1, m->size2);
3442 note_nonscalar (const gsl_matrix *m, const struct matrix_expr *e)
3445 note_operand_size (m, e);
3449 matrix_expr_evaluate_seq (const struct matrix_expr *e,
3450 gsl_matrix *start_, gsl_matrix *end_,
3453 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3455 msg_at (SE, matrix_expr_location (e),
3456 _("All operands of : operator must be scalars."));
3458 note_nonscalar (start_, e->subs[0]);
3459 note_nonscalar (end_, e->subs[1]);
3461 note_nonscalar (by_, e->subs[2]);
3465 long int start = to_scalar (start_);
3466 long int end = to_scalar (end_);
3467 long int by = by_ ? to_scalar (by_) : 1;
3471 msg_at (SE, matrix_expr_location (e->subs[2]),
3472 _("The increment operand to : must be nonzero."));
3476 long int n = (end >= start && by > 0 ? (end - start + by) / by
3477 : end <= start && by < 0 ? (start - end - by) / -by
3479 gsl_matrix *m = gsl_matrix_alloc (1, n);
3480 for (long int i = 0; i < n; i++)
3481 gsl_matrix_set (m, 0, i, start + i * by);
3486 matrix_expr_evaluate_not (gsl_matrix *a)
3488 MATRIX_FOR_ALL_ELEMENTS (d, y, x, a)
3494 matrix_expr_evaluate_paste_horz (const struct matrix_expr *e,
3495 gsl_matrix *a, gsl_matrix *b)
3497 if (a->size1 != b->size1)
3499 if (!a->size1 || !a->size2)
3501 else if (!b->size1 || !b->size2)
3504 msg_at (SE, matrix_expr_location (e),
3505 _("This expression tries to horizontally join matrices with "
3506 "differing numbers of rows."));
3507 note_operand_size (a, e->subs[0]);
3508 note_operand_size (b, e->subs[1]);
3512 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3513 for (size_t y = 0; y < a->size1; y++)
3515 for (size_t x = 0; x < a->size2; x++)
3516 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3517 for (size_t x = 0; x < b->size2; x++)
3518 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3524 matrix_expr_evaluate_paste_vert (const struct matrix_expr *e,
3525 gsl_matrix *a, gsl_matrix *b)
3527 if (a->size2 != b->size2)
3529 if (!a->size1 || !a->size2)
3531 else if (!b->size1 || !b->size2)
3534 msg_at (SE, matrix_expr_location (e),
3535 _("This expression tries to vertically join matrices with "
3536 "differing numbers of columns."));
3537 note_operand_size (a, e->subs[0]);
3538 note_operand_size (b, e->subs[1]);
3542 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3543 for (size_t x = 0; x < a->size2; x++)
3545 for (size_t y = 0; y < a->size1; y++)
3546 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3547 for (size_t y = 0; y < b->size1; y++)
3548 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3554 matrix_to_vector (gsl_matrix *m)
3557 gsl_vector v = to_vector (m);
3558 assert (v.block == m->block || !v.block);
3562 gsl_matrix_free (m);
3563 return xmemdup (&v, sizeof v);
3577 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3580 index_vector_uninit (struct index_vector *iv)
3587 matrix_normalize_index_vector (const gsl_matrix *m,
3588 const struct matrix_expr *me, size_t size,
3589 enum index_type index_type, size_t other_size,
3590 struct index_vector *iv)
3599 msg_at (SE, matrix_expr_location (me),
3600 _("Vector index must be scalar or vector, not a "
3602 m->size1, m->size2);
3606 msg_at (SE, matrix_expr_location (me),
3607 _("Matrix row index must be scalar or vector, not a "
3609 m->size1, m->size2);
3613 msg_at (SE, matrix_expr_location (me),
3614 _("Matrix column index must be scalar or vector, not a "
3616 m->size1, m->size2);
3622 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3623 *iv = (struct index_vector) {
3624 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3627 for (size_t i = 0; i < v.size; i++)
3629 double index = gsl_vector_get (&v, i);
3630 if (index < 1 || index >= size + 1)
3635 msg_at (SE, matrix_expr_location (me),
3636 _("Index %g is out of range for vector "
3637 "with %zu elements."), index, size);
3641 msg_at (SE, matrix_expr_location (me),
3642 _("%g is not a valid row index for "
3643 "a %zu×%zu matrix."),
3644 index, size, other_size);
3648 msg_at (SE, matrix_expr_location (me),
3649 _("%g is not a valid column index for "
3650 "a %zu×%zu matrix."),
3651 index, other_size, size);
3655 index_vector_uninit (iv);
3658 iv->indexes[i] = index - 1;
3664 *iv = (struct index_vector) {
3665 .indexes = xnmalloc (size, sizeof *iv->indexes),
3668 for (size_t i = 0; i < size; i++)
3675 matrix_expr_evaluate_vec_all (const struct matrix_expr *e,
3678 if (!is_vector (sm))
3680 msg_at (SE, matrix_expr_location (e->subs[0]),
3681 _("Vector index operator may not be applied to "
3682 "a %zu×%zu matrix."),
3683 sm->size1, sm->size2);
3691 matrix_expr_evaluate_vec_index (const struct matrix_expr *e,
3692 gsl_matrix *sm, gsl_matrix *im)
3694 if (!matrix_expr_evaluate_vec_all (e, sm))
3697 gsl_vector sv = to_vector (sm);
3698 struct index_vector iv;
3699 if (!matrix_normalize_index_vector (im, e->subs[1],
3700 sv.size, IV_VECTOR, 0, &iv))
3703 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3704 sm->size1 == 1 ? iv.n : 1);
3705 gsl_vector dv = to_vector (dm);
3706 for (size_t dx = 0; dx < iv.n; dx++)
3708 size_t sx = iv.indexes[dx];
3709 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3711 index_vector_uninit (&iv);
3717 matrix_expr_evaluate_mat_index (gsl_matrix *sm,
3718 gsl_matrix *im0, const struct matrix_expr *eim0,
3719 gsl_matrix *im1, const struct matrix_expr *eim1)
3721 struct index_vector iv0;
3722 if (!matrix_normalize_index_vector (im0, eim0, sm->size1,
3723 IV_ROW, sm->size2, &iv0))
3726 struct index_vector iv1;
3727 if (!matrix_normalize_index_vector (im1, eim1, sm->size2,
3728 IV_COLUMN, sm->size1, &iv1))
3730 index_vector_uninit (&iv0);
3734 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3735 for (size_t dy = 0; dy < iv0.n; dy++)
3737 size_t sy = iv0.indexes[dy];
3739 for (size_t dx = 0; dx < iv1.n; dx++)
3741 size_t sx = iv1.indexes[dx];
3742 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3745 index_vector_uninit (&iv0);
3746 index_vector_uninit (&iv1);
3750 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3751 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3752 const struct matrix_function_properties *, gsl_matrix *[], \
3753 const struct matrix_expr *, matrix_proto_##PROTO *);
3758 check_scalar_arg (const char *name, gsl_matrix *subs[],
3759 const struct matrix_expr *e, size_t index)
3761 if (!is_scalar (subs[index]))
3763 msg_at (SE, matrix_expr_location (e->subs[index]),
3764 _("Function %s argument %zu must be a scalar, "
3765 "not a %zu×%zu matrix."),
3766 name, index + 1, subs[index]->size1, subs[index]->size2);
3773 check_vector_arg (const char *name, gsl_matrix *subs[],
3774 const struct matrix_expr *e, size_t index)
3776 if (!is_vector (subs[index]))
3778 msg_at (SE, matrix_expr_location (e->subs[index]),
3779 _("Function %s argument %zu must be a vector, "
3780 "not a %zu×%zu matrix."),
3781 name, index + 1, subs[index]->size1, subs[index]->size2);
3788 to_scalar_args (const char *name, gsl_matrix *subs[],
3789 const struct matrix_expr *e, double d[])
3791 for (size_t i = 0; i < e->n_subs; i++)
3793 if (!check_scalar_arg (name, subs, e, i))
3795 d[i] = to_scalar (subs[i]);
3801 parse_constraint_value (const char **constraintsp)
3804 long retval = strtol (*constraintsp, &tail, 10);
3805 assert (tail > *constraintsp);
3806 *constraintsp = tail;
3810 enum matrix_argument_relop
3820 argument_inequality_error (
3821 const struct matrix_function_properties *props, const struct matrix_expr *e,
3822 size_t ai, gsl_matrix *a, size_t y, size_t x,
3823 size_t bi, double b,
3824 enum matrix_argument_relop relop)
3826 const struct msg_location *loc = matrix_expr_location (e);
3830 msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
3831 "than or equal to argument %zu."),
3832 ai + 1, props->name, bi + 1);
3836 msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
3837 "than argument %zu."),
3838 ai + 1, props->name, bi + 1);
3842 msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
3843 "or equal to argument %zu."),
3844 ai + 1, props->name, bi + 1);
3848 msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
3850 ai + 1, props->name, bi + 1);
3854 msg_at (ME, loc, _("Argument %zu to matrix function %s must not be equal "
3855 "to argument %zu."),
3856 ai + 1, props->name, bi + 1);
3860 const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
3862 msg_at (SN, a_loc, _("Argument %zu is %g."),
3863 ai + 1, gsl_matrix_get (a, y, x));
3865 msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
3866 y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
3868 msg_at (SN, matrix_expr_location (e->subs[bi]),
3869 _("Argument %zu is %g."), bi + 1, b);
3873 argument_value_error (
3874 const struct matrix_function_properties *props, const struct matrix_expr *e,
3875 size_t ai, gsl_matrix *a, size_t y, size_t x,
3877 enum matrix_argument_relop relop)
3879 const struct msg_location *loc = matrix_expr_location (e);
3883 msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
3884 "than or equal to %g."),
3885 ai + 1, props->name, b);
3889 msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
3891 ai + 1, props->name, b);
3895 msg_at (SE, loc, _("Argument %zu to matrix function %s must be less than "
3897 ai + 1, props->name, b);
3901 msg_at (SE, loc, _("Argument %zu to matrix function %s must be less than "
3903 ai + 1, props->name, b);
3907 msg_at (SE, loc, _("Argument %zu to matrix function %s must not be equal "
3909 ai + 1, props->name, b);
3913 const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
3916 if (relop != MRR_NE)
3917 msg_at (SN, a_loc, _("Argument %zu is %g."),
3918 ai + 1, gsl_matrix_get (a, y, x));
3921 msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
3922 y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
3926 matrix_argument_relop_is_satisfied (double a, double b,
3927 enum matrix_argument_relop relop)
3931 case MRR_GE: return a >= b;
3932 case MRR_GT: return a > b;
3933 case MRR_LE: return a <= b;
3934 case MRR_LT: return a < b;
3935 case MRR_NE: return a != b;
3941 static enum matrix_argument_relop
3942 matrix_argument_relop_flip (enum matrix_argument_relop relop)
3946 case MRR_GE: return MRR_LE;
3947 case MRR_GT: return MRR_LT;
3948 case MRR_LE: return MRR_GE;
3949 case MRR_LT: return MRR_GT;
3950 case MRR_NE: return MRR_NE;
3957 check_constraints (const struct matrix_function_properties *props,
3958 gsl_matrix *args[], const struct matrix_expr *e)
3960 size_t n_args = e->n_subs;
3961 const char *constraints = props->constraints;
3965 size_t arg_index = SIZE_MAX;
3966 while (*constraints)
3968 if (*constraints >= 'a' && *constraints <= 'd')
3970 arg_index = *constraints++ - 'a';
3971 assert (arg_index < n_args);
3973 else if (*constraints == '[' || *constraints == '(')
3975 assert (arg_index < n_args);
3976 bool open_lower = *constraints++ == '(';
3977 int minimum = parse_constraint_value (&constraints);
3978 assert (*constraints == ',');
3980 int maximum = parse_constraint_value (&constraints);
3981 assert (*constraints == ']' || *constraints == ')');
3982 bool open_upper = *constraints++ == ')';
3984 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3985 if ((open_lower ? *d <= minimum : *d < minimum)
3986 || (open_upper ? *d >= maximum : *d > maximum))
3988 if (!is_scalar (args[arg_index]))
3989 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
3990 _("Row %zu, column %zu of argument %zu to matrix "
3991 "function %s is %g, which is outside "
3992 "the valid range %c%d,%d%c."),
3993 y + 1, x + 1, arg_index + 1, props->name, *d,
3994 open_lower ? '(' : '[',
3996 open_upper ? ')' : ']');
3998 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
3999 _("Argument %zu to matrix function %s is %g, "
4000 "which is outside the valid range %c%d,%d%c."),
4001 arg_index + 1, props->name, *d,
4002 open_lower ? '(' : '[',
4004 open_upper ? ')' : ']');
4008 else if (*constraints == 'i')
4011 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
4012 if (*d != floor (*d))
4014 if (!is_scalar (args[arg_index]))
4015 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
4016 _("Argument %zu to matrix function %s, which must be "
4017 "integer, contains non-integer value %g in "
4018 "row %zu, column %zu."),
4019 arg_index + 1, props->name, *d, y + 1, x + 1);
4021 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
4022 _("Argument %zu to matrix function %s, which must be "
4023 "integer, has non-integer value %g."),
4024 arg_index + 1, props->name, *d);
4028 else if (*constraints == '>'
4029 || *constraints == '<'
4030 || *constraints == '!')
4032 enum matrix_argument_relop relop;
4033 switch (*constraints++)
4036 if (*constraints == '=')
4046 if (*constraints == '=')
4056 assert (*constraints == '=');
4065 if (*constraints >= 'a' && *constraints <= 'd')
4067 size_t a_index = arg_index;
4068 size_t b_index = *constraints - 'a';
4069 assert (a_index < n_args);
4070 assert (b_index < n_args);
4072 /* We only support one of the two arguments being non-scalar.
4073 It's easier to support only the first one being non-scalar, so
4074 flip things around if it's the other way. */
4075 if (!is_scalar (args[b_index]))
4077 assert (is_scalar (args[a_index]));
4078 size_t tmp_index = a_index;
4080 b_index = tmp_index;
4081 relop = matrix_argument_relop_flip (relop);
4084 double b = to_scalar (args[b_index]);
4085 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
4086 if (!matrix_argument_relop_is_satisfied (*a, b, relop))
4088 argument_inequality_error (
4090 a_index, args[a_index], y, x,
4098 int comparand = parse_constraint_value (&constraints);
4100 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
4101 if (!matrix_argument_relop_is_satisfied (*d, comparand, relop))
4103 argument_value_error (
4105 arg_index, args[arg_index], y, x,
4114 assert (*constraints == ' ');
4116 arg_index = SIZE_MAX;
4123 matrix_expr_evaluate_d_none (const struct matrix_function_properties *props,
4124 gsl_matrix *subs[], const struct matrix_expr *e,
4125 matrix_proto_d_none *f)
4127 assert (e->n_subs == 0);
4129 if (!check_constraints (props, subs, e))
4132 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4133 gsl_matrix_set (m, 0, 0, f ());
4138 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
4139 gsl_matrix *subs[], const struct matrix_expr *e,
4140 matrix_proto_d_d *f)
4142 assert (e->n_subs == 1);
4145 if (!to_scalar_args (props->name, subs, e, &d)
4146 || !check_constraints (props, subs, e))
4149 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4150 gsl_matrix_set (m, 0, 0, f (d));
4155 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
4156 gsl_matrix *subs[], const struct matrix_expr *e,
4157 matrix_proto_d_dd *f)
4159 assert (e->n_subs == 2);
4162 if (!to_scalar_args (props->name, subs, e, d)
4163 && !check_constraints (props, subs, e))
4166 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4167 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
4172 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
4173 gsl_matrix *subs[], const struct matrix_expr *e,
4174 matrix_proto_d_ddd *f)
4176 assert (e->n_subs == 3);
4179 if (!to_scalar_args (props->name, subs, e, d)
4180 || !check_constraints (props, subs, e))
4183 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4184 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
4189 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
4190 gsl_matrix *subs[], const struct matrix_expr *e,
4191 matrix_proto_m_d *f)
4193 assert (e->n_subs == 1);
4196 return (to_scalar_args (props->name, subs, e, &d)
4197 && check_constraints (props, subs, e)
4203 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
4204 gsl_matrix *subs[], const struct matrix_expr *e,
4205 matrix_proto_m_ddd *f)
4207 assert (e->n_subs == 3);
4210 return (to_scalar_args (props->name, subs, e, d)
4211 && check_constraints (props, subs, e)
4212 ? f(d[0], d[1], d[2])
4217 matrix_expr_evaluate_m_ddn (const struct matrix_function_properties *props,
4218 gsl_matrix *subs[], const struct matrix_expr *e,
4219 matrix_proto_m_ddn *f)
4221 assert (e->n_subs == 2);
4224 return (to_scalar_args (props->name, subs, e, d)
4225 && check_constraints (props, subs, e)
4231 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props,
4232 gsl_matrix *subs[], const struct matrix_expr *e,
4233 matrix_proto_m_m *f)
4235 assert (e->n_subs == 1);
4236 return check_constraints (props, subs, e) ? f (subs[0]) : NULL;
4240 matrix_expr_evaluate_m_mn (const struct matrix_function_properties *props,
4241 gsl_matrix *subs[], const struct matrix_expr *e,
4242 matrix_proto_m_mn *f)
4244 assert (e->n_subs == 1);
4245 return check_constraints (props, subs, e) ? f (subs[0], e) : NULL;
4249 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
4250 gsl_matrix *subs[], const struct matrix_expr *e,
4251 matrix_proto_m_e *f)
4253 assert (e->n_subs == 1);
4255 if (!check_constraints (props, subs, e))
4258 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4264 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props,
4265 gsl_matrix *subs[], const struct matrix_expr *e,
4266 matrix_proto_m_md *f)
4268 assert (e->n_subs == 2);
4269 return (check_scalar_arg (props->name, subs, e, 1)
4270 && check_constraints (props, subs, e)
4271 ? f (subs[0], to_scalar (subs[1]))
4276 matrix_expr_evaluate_m_mdn (const struct matrix_function_properties *props,
4277 gsl_matrix *subs[], const struct matrix_expr *e,
4278 matrix_proto_m_mdn *f)
4280 assert (e->n_subs == 2);
4281 return (check_scalar_arg (props->name, subs, e, 1)
4282 && check_constraints (props, subs, e)
4283 ? f (subs[0], to_scalar (subs[1]), e)
4288 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
4289 gsl_matrix *subs[], const struct matrix_expr *e,
4290 matrix_proto_m_ed *f)
4292 assert (e->n_subs == 2);
4293 if (!check_scalar_arg (props->name, subs, e, 1)
4294 || !check_constraints (props, subs, e))
4297 double b = to_scalar (subs[1]);
4298 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4304 matrix_expr_evaluate_m_mddn (const struct matrix_function_properties *props,
4305 gsl_matrix *subs[], const struct matrix_expr *e,
4306 matrix_proto_m_mddn *f)
4308 assert (e->n_subs == 3);
4309 if (!check_scalar_arg (props->name, subs, e, 1)
4310 || !check_scalar_arg (props->name, subs, e, 2)
4311 || !check_constraints (props, subs, e))
4313 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]), e);
4317 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
4318 gsl_matrix *subs[], const struct matrix_expr *e,
4319 matrix_proto_m_edd *f)
4321 assert (e->n_subs == 3);
4322 if (!check_scalar_arg (props->name, subs, e, 1)
4323 || !check_scalar_arg (props->name, subs, e, 2)
4324 || !check_constraints (props, subs, e))
4327 double b = to_scalar (subs[1]);
4328 double c = to_scalar (subs[2]);
4329 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4335 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4336 gsl_matrix *subs[], const struct matrix_expr *e,
4337 matrix_proto_m_eddd *f)
4339 assert (e->n_subs == 4);
4340 for (size_t i = 1; i < 4; i++)
4341 if (!check_scalar_arg (props->name, subs, e, i))
4344 if (!check_constraints (props, subs, e))
4347 double b = to_scalar (subs[1]);
4348 double c = to_scalar (subs[2]);
4349 double d = to_scalar (subs[3]);
4350 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4351 *a = f (*a, b, c, d);
4356 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4357 gsl_matrix *subs[], const struct matrix_expr *e,
4358 matrix_proto_m_eed *f)
4360 assert (e->n_subs == 3);
4361 if (!check_scalar_arg (props->name, subs, e, 2))
4364 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4365 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4367 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
4368 loc->end = e->subs[1]->location->end;
4371 _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4372 "%zu×%zu, but %s requires these arguments either to have "
4373 "the same dimensions or for one of them to be a scalar."),
4375 subs[0]->size1, subs[0]->size2,
4376 subs[1]->size1, subs[1]->size2,
4379 msg_location_destroy (loc);
4383 if (!check_constraints (props, subs, e))
4386 double c = to_scalar (subs[2]);
4388 if (is_scalar (subs[0]))
4390 double a = to_scalar (subs[0]);
4391 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4397 double b = to_scalar (subs[1]);
4398 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4405 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props,
4406 gsl_matrix *subs[], const struct matrix_expr *e,
4407 matrix_proto_m_mm *f)
4409 assert (e->n_subs == 2);
4410 return check_constraints (props, subs, e) ? f (subs[0], subs[1]) : NULL;
4414 matrix_expr_evaluate_m_mmn (const struct matrix_function_properties *props,
4415 gsl_matrix *subs[], const struct matrix_expr *e,
4416 matrix_proto_m_mmn *f)
4418 assert (e->n_subs == 2);
4419 return check_constraints (props, subs, e) ? f (subs[0], subs[1], e) : NULL;
4423 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4424 gsl_matrix *subs[], const struct matrix_expr *e,
4425 matrix_proto_m_v *f)
4427 assert (e->n_subs == 1);
4428 if (!check_vector_arg (props->name, subs, e, 0)
4429 || !check_constraints (props, subs, e))
4431 gsl_vector v = to_vector (subs[0]);
4436 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props,
4437 gsl_matrix *subs[], const struct matrix_expr *e,
4438 matrix_proto_d_m *f)
4440 assert (e->n_subs == 1);
4442 if (!check_constraints (props, subs, e))
4445 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4446 gsl_matrix_set (m, 0, 0, f (subs[0]));
4451 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props,
4452 gsl_matrix *subs[], const struct matrix_expr *e,
4453 matrix_proto_m_any *f)
4455 return check_constraints (props, subs, e) ? f (subs, e->n_subs) : NULL;
4459 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props_ UNUSED,
4460 gsl_matrix *subs[], const struct matrix_expr *e,
4461 matrix_proto_IDENT *f)
4463 static const struct matrix_function_properties p1 = {
4465 .constraints = "ai>=0"
4467 static const struct matrix_function_properties p2 = {
4469 .constraints = "ai>=0 bi>=0"
4471 const struct matrix_function_properties *props = e->n_subs == 1 ? &p1 : &p2;
4473 assert (e->n_subs <= 2);
4476 return (to_scalar_args (props->name, subs, e, d)
4477 && check_constraints (props, subs, e)
4478 ? f (d[0], d[e->n_subs - 1])
4483 matrix_expr_evaluate (const struct matrix_expr *e)
4485 if (e->op == MOP_NUMBER)
4487 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4488 gsl_matrix_set (m, 0, 0, e->number);
4491 else if (e->op == MOP_VARIABLE)
4493 const gsl_matrix *src = e->variable->value;
4496 msg_at (SE, e->location,
4497 _("Uninitialized variable %s used in expression."),
4502 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4503 gsl_matrix_memcpy (dst, src);
4506 else if (e->op == MOP_EOF)
4508 struct dfm_reader *reader = read_file_open (e->eof);
4509 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4510 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4514 enum { N_LOCAL = 3 };
4515 gsl_matrix *local_subs[N_LOCAL];
4516 gsl_matrix **subs = (e->n_subs < N_LOCAL
4518 : xmalloc (e->n_subs * sizeof *subs));
4520 for (size_t i = 0; i < e->n_subs; i++)
4522 subs[i] = matrix_expr_evaluate (e->subs[i]);
4525 for (size_t j = 0; j < i; j++)
4526 gsl_matrix_free (subs[j]);
4527 if (subs != local_subs)
4533 gsl_matrix *result = NULL;
4536 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4537 case MOP_F_##ENUM: \
4539 static const struct matrix_function_properties props = { \
4541 .constraints = CONSTRAINTS, \
4543 result = matrix_expr_evaluate_##PROTO (&props, subs, e, \
4544 matrix_eval_##ENUM); \
4551 gsl_matrix_scale (subs[0], -1.0);
4569 result = matrix_expr_evaluate_elementwise (e, e->op, subs[0], subs[1]);
4573 result = matrix_expr_evaluate_not (subs[0]);
4577 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], NULL);
4581 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], subs[2]);
4585 result = matrix_expr_evaluate_mul_mat (e, subs[0], subs[1]);
4589 result = matrix_expr_evaluate_exp_mat (e, subs[0], subs[1]);
4592 case MOP_PASTE_HORZ:
4593 result = matrix_expr_evaluate_paste_horz (e, subs[0], subs[1]);
4596 case MOP_PASTE_VERT:
4597 result = matrix_expr_evaluate_paste_vert (e, subs[0], subs[1]);
4601 result = gsl_matrix_alloc (0, 0);
4605 result = matrix_expr_evaluate_vec_index (e, subs[0], subs[1]);
4609 result = matrix_expr_evaluate_vec_all (e, subs[0]);
4613 result = matrix_expr_evaluate_mat_index (subs[0],
4614 subs[1], e->subs[1],
4615 subs[2], e->subs[2]);
4619 result = matrix_expr_evaluate_mat_index (subs[0],
4620 subs[1], e->subs[1],
4625 result = matrix_expr_evaluate_mat_index (subs[0],
4627 subs[1], e->subs[1]);
4636 for (size_t i = 0; i < e->n_subs; i++)
4637 if (subs[i] != result)
4638 gsl_matrix_free (subs[i]);
4639 if (subs != local_subs)
4645 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4648 gsl_matrix *m = matrix_expr_evaluate (e);
4654 msg_at (SE, matrix_expr_location (e),
4655 _("Expression for %s must evaluate to scalar, "
4656 "not a %zu×%zu matrix."),
4657 context, m->size1, m->size2);
4658 gsl_matrix_free (m);
4663 gsl_matrix_free (m);
4668 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4672 if (!matrix_expr_evaluate_scalar (e, context, &d))
4676 if (d < LONG_MIN || d > LONG_MAX)
4678 msg_at (SE, matrix_expr_location (e),
4679 _("Expression for %s is outside the integer range."), context);
4688 An lvalue is an expression that can appear on the left side of a COMPUTE
4689 command and in other contexts that assign values.
4691 An lvalue is parsed once, with matrix_lvalue_parse(). It can then be
4692 evaluated (with matrix_lvalue_evaluate()) and assigned (with
4693 matrix_lvalue_assign()).
4695 There are three kinds of lvalues:
4697 - A variable name. A variable used as an lvalue need not be initialized,
4698 since the assignment will initialize.
4700 - A subvector, e.g. "var(index0)". The variable must be initialized and
4701 must have the form of a vector (it must have 1 column or 1 row).
4703 - A submatrix, e.g. "var(index0, index1)". The variable must be
4705 struct matrix_lvalue
4707 struct matrix_var *var; /* Destination variable. */
4708 struct matrix_expr *indexes[2]; /* Index expressions, if any. */
4709 size_t n_indexes; /* Number of indexes. */
4711 struct msg_location *var_location; /* Variable name. */
4712 struct msg_location *full_location; /* Variable name plus indexing. */
4713 struct msg_location *index_locations[2]; /* Index expressions. */
4718 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4722 msg_location_destroy (lvalue->var_location);
4723 msg_location_destroy (lvalue->full_location);
4724 for (size_t i = 0; i < lvalue->n_indexes; i++)
4726 matrix_expr_destroy (lvalue->indexes[i]);
4727 msg_location_destroy (lvalue->index_locations[i]);
4733 /* Parses and returns an lvalue at the current position in S's lexer. Returns
4734 null on parse failure. On success, the caller must eventually free the
4736 static struct matrix_lvalue *
4737 matrix_lvalue_parse (struct matrix_state *s)
4739 if (!lex_force_id (s->lexer))
4742 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4743 int start_ofs = lex_ofs (s->lexer);
4744 lvalue->var_location = lex_get_location (s->lexer, 0, 0);
4745 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4746 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4750 lex_error (s->lexer, _("Undefined variable %s."),
4751 lex_tokcstr (s->lexer));
4755 lex_get_n (s->lexer, 2);
4757 if (!matrix_parse_index_expr (s, &lvalue->indexes[0],
4758 &lvalue->index_locations[0]))
4760 lvalue->n_indexes++;
4762 if (lex_match (s->lexer, T_COMMA))
4764 if (!matrix_parse_index_expr (s, &lvalue->indexes[1],
4765 &lvalue->index_locations[1]))
4767 lvalue->n_indexes++;
4769 if (!lex_force_match (s->lexer, T_RPAREN))
4772 lvalue->full_location = lex_ofs_location (s->lexer, start_ofs,
4773 lex_ofs (s->lexer) - 1);
4778 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4784 matrix_lvalue_destroy (lvalue);
4789 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4790 enum index_type index_type, size_t other_size,
4791 struct index_vector *iv)
4796 m = matrix_expr_evaluate (e);
4803 bool ok = matrix_normalize_index_vector (m, e, size, index_type,
4805 gsl_matrix_free (m);
4809 /* Evaluates the indexes in LVALUE into IV0 and IV1, owned by the caller.
4810 Returns true if successful, false if evaluating the expressions failed or if
4811 LVALUE otherwise can't be used for an assignment.
4813 On success, the caller retains ownership of the index vectors, which are
4814 suitable for passing to matrix_lvalue_assign(). If not used for that
4815 purpose then they need to eventually be freed (with
4816 index_vector_uninit()). */
4818 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4819 struct index_vector *iv0,
4820 struct index_vector *iv1)
4822 *iv0 = INDEX_VECTOR_INIT;
4823 *iv1 = INDEX_VECTOR_INIT;
4824 if (!lvalue->n_indexes)
4827 /* Validate destination matrix exists and has the right shape. */
4828 gsl_matrix *dm = lvalue->var->value;
4831 msg_at (SE, lvalue->var_location,
4832 _("Undefined variable %s."), lvalue->var->name);
4835 else if (dm->size1 == 0 || dm->size2 == 0)
4837 msg_at (SE, lvalue->full_location, _("Cannot index %zu×%zu matrix %s."),
4838 dm->size1, dm->size2, lvalue->var->name);
4841 else if (lvalue->n_indexes == 1)
4843 if (!is_vector (dm))
4845 msg_at (SE, lvalue->full_location,
4846 _("Can't use vector indexing on %zu×%zu matrix %s."),
4847 dm->size1, dm->size2, lvalue->var->name);
4850 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4851 MAX (dm->size1, dm->size2),
4856 assert (lvalue->n_indexes == 2);
4857 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4858 IV_ROW, dm->size2, iv0))
4861 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4862 IV_COLUMN, dm->size1, iv1))
4864 index_vector_uninit (iv0);
4872 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4873 struct index_vector *iv,
4874 gsl_matrix *sm, const struct msg_location *lsm)
4876 /* Convert source matrix 'sm' to source vector 'sv'. */
4877 if (!is_vector (sm))
4879 msg_at (SE, lvalue->full_location,
4880 _("Only an %zu-element vector may be assigned to this "
4881 "%zu-element subvector of %s."),
4882 iv->n, iv->n, lvalue->var->name);
4884 _("The source is an %zu×%zu matrix."),
4885 sm->size1, sm->size2);
4888 gsl_vector sv = to_vector (sm);
4889 if (iv->n != sv.size)
4891 msg_at (SE, lvalue->full_location,
4892 _("Only an %zu-element vector may be assigned to this "
4893 "%zu-element subvector of %s."),
4894 iv->n, iv->n, lvalue->var->name);
4895 msg_at (SE, lsm, ngettext ("The source vector has %zu element.",
4896 "The source vector has %zu elements.",
4902 /* Assign elements. */
4903 gsl_vector dv = to_vector (lvalue->var->value);
4904 for (size_t x = 0; x < iv->n; x++)
4905 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4910 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4911 struct index_vector *iv0,
4912 struct index_vector *iv1,
4913 gsl_matrix *sm, const struct msg_location *lsm)
4915 gsl_matrix *dm = lvalue->var->value;
4917 /* Convert source matrix 'sm' to source vector 'sv'. */
4918 bool wrong_rows = iv0->n != sm->size1;
4919 bool wrong_columns = iv1->n != sm->size2;
4920 if (wrong_rows || wrong_columns)
4922 if (wrong_rows && wrong_columns)
4923 msg_at (SE, lvalue->full_location,
4924 _("Numbers of indexes for assigning to %s differ from the "
4925 "size of the source matrix."),
4927 else if (wrong_rows)
4928 msg_at (SE, lvalue->full_location,
4929 _("Number of row indexes for assigning to %s differs from "
4930 "number of rows in source matrix."),
4933 msg_at (SE, lvalue->full_location,
4934 _("Number of column indexes for assigning to %s differs from "
4935 "number of columns in source matrix."),
4940 if (lvalue->indexes[0])
4941 msg_at (SN, lvalue->index_locations[0],
4942 ngettext ("There is %zu row index.",
4943 "There are %zu row indexes.",
4947 msg_at (SN, lvalue->index_locations[0],
4948 ngettext ("Destination matrix %s has %zu row.",
4949 "Destination matrix %s has %zu rows.",
4951 lvalue->var->name, iv0->n);
4956 if (lvalue->indexes[1])
4957 msg_at (SN, lvalue->index_locations[1],
4958 ngettext ("There is %zu column index.",
4959 "There are %zu column indexes.",
4963 msg_at (SN, lvalue->index_locations[1],
4964 ngettext ("Destination matrix %s has %zu column.",
4965 "Destination matrix %s has %zu columns.",
4967 lvalue->var->name, iv1->n);
4970 msg_at (SN, lsm, _("The source matrix is %zu×%zu."),
4971 sm->size1, sm->size2);
4975 /* Assign elements. */
4976 for (size_t y = 0; y < iv0->n; y++)
4978 size_t f0 = iv0->indexes[y];
4979 for (size_t x = 0; x < iv1->n; x++)
4981 size_t f1 = iv1->indexes[x];
4982 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4988 /* Assigns rvalue SM to LVALUE, whose index vectors IV0 and IV1 were previously
4989 obtained with matrix_lvalue_evaluate(). Returns true if successful, false
4990 on error. Always takes ownership of M. LSM should be the source location
4991 to use for errors related to SM. */
4993 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4994 struct index_vector *iv0, struct index_vector *iv1,
4995 gsl_matrix *sm, const struct msg_location *lsm)
4997 if (!lvalue->n_indexes)
4999 gsl_matrix_free (lvalue->var->value);
5000 lvalue->var->value = sm;
5005 bool ok = (lvalue->n_indexes == 1
5006 ? matrix_lvalue_assign_vector (lvalue, iv0, sm, lsm)
5007 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm, lsm));
5008 index_vector_uninit (iv0);
5009 index_vector_uninit (iv1);
5010 gsl_matrix_free (sm);
5015 /* Evaluates and then assigns SM to LVALUE. Always takes ownership of M. LSM
5016 should be the source location to use for errors related to SM.*/
5018 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue,
5020 const struct msg_location *lsm)
5022 struct index_vector iv0, iv1;
5023 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
5025 gsl_matrix_free (sm);
5029 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm, lsm);
5032 /* Matrix command data structure. */
5034 /* An array of matrix commands. */
5035 struct matrix_commands
5037 struct matrix_command **commands;
5041 static bool matrix_commands_parse (struct matrix_state *,
5042 struct matrix_commands *,
5043 const char *command_name,
5044 const char *stop1, const char *stop2);
5045 static void matrix_commands_uninit (struct matrix_commands *);
5047 /* A single matrix command. */
5048 struct matrix_command
5050 /* The type of command. */
5051 enum matrix_command_type
5072 /* Source lines for this command. */
5073 struct msg_location *location;
5077 struct matrix_compute
5079 struct matrix_lvalue *lvalue;
5080 struct matrix_expr *rvalue;
5086 struct matrix_expr *expression;
5087 bool use_default_format;
5088 struct fmt_spec format;
5090 int space; /* -1 means NEWPAGE. */
5094 struct string_array literals; /* CLABELS/RLABELS. */
5095 struct matrix_expr *expr; /* CNAMES/RNAMES. */
5105 struct matrix_expr *condition;
5106 struct matrix_commands commands;
5116 struct matrix_var *var;
5117 struct matrix_expr *start;
5118 struct matrix_expr *end;
5119 struct matrix_expr *increment;
5121 /* Loop conditions. */
5122 struct matrix_expr *top_condition;
5123 struct matrix_expr *bottom_condition;
5126 struct matrix_commands commands;
5130 struct matrix_display
5132 struct matrix_state *state;
5136 struct matrix_release
5138 struct matrix_var **vars;
5145 struct matrix_expr *expression;
5146 struct save_file *sf;
5152 struct read_file *rf;
5153 struct matrix_lvalue *dst;
5154 struct matrix_expr *size;
5156 enum fmt_type format;
5165 struct write_file *wf;
5166 struct matrix_expr *expression;
5169 /* If this is nonnull, WRITE uses this format.
5171 If this is NULL, WRITE uses free-field format with as many
5172 digits of precision as needed. */
5173 struct fmt_spec *format;
5182 struct lexer *lexer;
5183 struct matrix_lvalue *dst;
5184 struct dataset *dataset;
5185 struct file_handle *file;
5187 struct var_syntax *vars;
5189 struct matrix_var *names;
5191 /* Treatment of missing values. */
5196 MGET_ERROR, /* Flag error on command. */
5197 MGET_ACCEPT, /* Accept without change, user-missing only. */
5198 MGET_OMIT, /* Drop this case. */
5199 MGET_RECODE /* Recode to 'substitute'. */
5202 double substitute; /* MGET_RECODE only. */
5210 struct msave_common *common;
5211 struct matrix_expr *expr;
5212 const char *rowtype;
5213 const struct matrix_expr *factors;
5214 const struct matrix_expr *splits;
5220 struct matrix_state *state;
5221 struct file_handle *file;
5223 struct stringi_set rowtypes;
5229 struct matrix_expr *expr;
5230 struct matrix_var *evec;
5231 struct matrix_var *eval;
5235 struct matrix_setdiag
5237 struct matrix_var *dst;
5238 struct matrix_expr *expr;
5244 struct matrix_expr *expr;
5245 struct matrix_var *u;
5246 struct matrix_var *s;
5247 struct matrix_var *v;
5253 static struct matrix_command *matrix_command_parse (struct matrix_state *);
5254 static bool matrix_command_execute (struct matrix_command *);
5255 static void matrix_command_destroy (struct matrix_command *);
5259 static struct matrix_command *
5260 matrix_compute_parse (struct matrix_state *s)
5262 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5263 *cmd = (struct matrix_command) {
5264 .type = MCMD_COMPUTE,
5265 .compute = { .lvalue = NULL },
5268 struct matrix_compute *compute = &cmd->compute;
5269 compute->lvalue = matrix_lvalue_parse (s);
5270 if (!compute->lvalue)
5273 if (!lex_force_match (s->lexer, T_EQUALS))
5276 compute->rvalue = matrix_expr_parse (s);
5277 if (!compute->rvalue)
5283 matrix_command_destroy (cmd);
5288 matrix_compute_execute (struct matrix_command *cmd)
5290 struct matrix_compute *compute = &cmd->compute;
5291 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
5295 matrix_lvalue_evaluate_and_assign (compute->lvalue, value,
5296 matrix_expr_location (compute->rvalue));
5302 print_labels_destroy (struct print_labels *labels)
5306 string_array_destroy (&labels->literals);
5307 matrix_expr_destroy (labels->expr);
5313 parse_literal_print_labels (struct matrix_state *s,
5314 struct print_labels **labelsp)
5316 lex_match (s->lexer, T_EQUALS);
5317 print_labels_destroy (*labelsp);
5318 *labelsp = xzalloc (sizeof **labelsp);
5319 while (lex_token (s->lexer) != T_SLASH
5320 && lex_token (s->lexer) != T_ENDCMD
5321 && lex_token (s->lexer) != T_STOP)
5323 struct string label = DS_EMPTY_INITIALIZER;
5324 while (lex_token (s->lexer) != T_COMMA
5325 && lex_token (s->lexer) != T_SLASH
5326 && lex_token (s->lexer) != T_ENDCMD
5327 && lex_token (s->lexer) != T_STOP)
5329 if (!ds_is_empty (&label))
5330 ds_put_byte (&label, ' ');
5332 if (lex_token (s->lexer) == T_STRING)
5333 ds_put_cstr (&label, lex_tokcstr (s->lexer));
5336 char *rep = lex_next_representation (s->lexer, 0, 0);
5337 ds_put_cstr (&label, rep);
5342 string_array_append_nocopy (&(*labelsp)->literals,
5343 ds_steal_cstr (&label));
5345 if (!lex_match (s->lexer, T_COMMA))
5351 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
5353 lex_match (s->lexer, T_EQUALS);
5354 struct matrix_expr *e = matrix_parse_exp (s);
5358 print_labels_destroy (*labelsp);
5359 *labelsp = xzalloc (sizeof **labelsp);
5360 (*labelsp)->expr = e;
5364 static struct matrix_command *
5365 matrix_print_parse (struct matrix_state *s)
5367 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5368 *cmd = (struct matrix_command) {
5371 .use_default_format = true,
5375 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
5377 int start_ofs = lex_ofs (s->lexer);
5378 cmd->print.expression = matrix_parse_exp (s);
5379 if (!cmd->print.expression)
5381 cmd->print.title = lex_ofs_representation (s->lexer, start_ofs,
5382 lex_ofs (s->lexer) - 1);
5385 while (lex_match (s->lexer, T_SLASH))
5387 if (lex_match_id (s->lexer, "FORMAT"))
5389 lex_match (s->lexer, T_EQUALS);
5390 if (!parse_format_specifier (s->lexer, &cmd->print.format))
5392 cmd->print.use_default_format = false;
5394 else if (lex_match_id (s->lexer, "TITLE"))
5396 lex_match (s->lexer, T_EQUALS);
5397 if (!lex_force_string (s->lexer))
5399 free (cmd->print.title);
5400 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
5403 else if (lex_match_id (s->lexer, "SPACE"))
5405 lex_match (s->lexer, T_EQUALS);
5406 if (lex_match_id (s->lexer, "NEWPAGE"))
5407 cmd->print.space = -1;
5408 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
5410 cmd->print.space = lex_integer (s->lexer);
5416 else if (lex_match_id (s->lexer, "RLABELS"))
5417 parse_literal_print_labels (s, &cmd->print.rlabels);
5418 else if (lex_match_id (s->lexer, "CLABELS"))
5419 parse_literal_print_labels (s, &cmd->print.clabels);
5420 else if (lex_match_id (s->lexer, "RNAMES"))
5422 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
5425 else if (lex_match_id (s->lexer, "CNAMES"))
5427 if (!parse_expr_print_labels (s, &cmd->print.clabels))
5432 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
5433 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
5441 matrix_command_destroy (cmd);
5446 matrix_is_integer (const gsl_matrix *m)
5448 for (size_t y = 0; y < m->size1; y++)
5449 for (size_t x = 0; x < m->size2; x++)
5451 double d = gsl_matrix_get (m, y, x);
5459 matrix_max_magnitude (const gsl_matrix *m)
5462 for (size_t y = 0; y < m->size1; y++)
5463 for (size_t x = 0; x < m->size2; x++)
5465 double d = fabs (gsl_matrix_get (m, y, x));
5473 format_fits (struct fmt_spec format, double x)
5475 char *s = data_out (&(union value) { .f = x }, NULL,
5476 format, settings_get_fmt_settings ());
5477 bool fits = *s != '*' && !strchr (s, 'E');
5482 static struct fmt_spec
5483 default_format (const gsl_matrix *m, int *log_scale)
5487 double max = matrix_max_magnitude (m);
5489 if (matrix_is_integer (m))
5490 for (int w = 1; w <= 12; w++)
5492 struct fmt_spec format = { .type = FMT_F, .w = w };
5493 if (format_fits (format, -max))
5497 if (max >= 1e9 || max <= 1e-4)
5500 snprintf (s, sizeof s, "%.1e", max);
5502 const char *e = strchr (s, 'e');
5504 *log_scale = atoi (e + 1);
5507 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5511 trimmed_string (double d)
5513 struct substring s = ss_buffer ((char *) &d, sizeof d);
5514 ss_rtrim (&s, ss_cstr (" "));
5515 return ss_xstrdup (s);
5518 static struct string_array *
5519 print_labels_get (const struct print_labels *labels, size_t n,
5520 const char *prefix, bool truncate)
5525 struct string_array *out = xzalloc (sizeof *out);
5526 if (labels->literals.n)
5527 string_array_clone (out, &labels->literals);
5528 else if (labels->expr)
5530 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5531 if (m && is_vector (m))
5533 gsl_vector v = to_vector (m);
5534 for (size_t i = 0; i < v.size; i++)
5535 string_array_append_nocopy (out, trimmed_string (
5536 gsl_vector_get (&v, i)));
5538 gsl_matrix_free (m);
5544 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5547 string_array_append (out, "");
5550 string_array_delete (out, out->n - 1);
5553 for (size_t i = 0; i < out->n; i++)
5555 char *s = out->strings[i];
5556 s[strnlen (s, 8)] = '\0';
5563 matrix_print_space (int space)
5566 output_item_submit (page_break_item_create ());
5567 for (int i = 0; i < space; i++)
5568 output_log ("%s", "");
5572 matrix_print_text (const struct matrix_print *print, const gsl_matrix *m,
5573 struct fmt_spec format, int log_scale)
5575 matrix_print_space (print->space);
5577 output_log ("%s", print->title);
5579 output_log (" 10 ** %d X", log_scale);
5581 struct string_array *clabels = print_labels_get (print->clabels,
5582 m->size2, "col", true);
5583 if (clabels && format.w < 8)
5585 struct string_array *rlabels = print_labels_get (print->rlabels,
5586 m->size1, "row", true);
5590 struct string line = DS_EMPTY_INITIALIZER;
5592 ds_put_byte_multiple (&line, ' ', 8);
5593 for (size_t x = 0; x < m->size2; x++)
5594 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5595 output_log_nocopy (ds_steal_cstr (&line));
5598 double scale = pow (10.0, log_scale);
5599 bool numeric = fmt_is_numeric (format.type);
5600 for (size_t y = 0; y < m->size1; y++)
5602 struct string line = DS_EMPTY_INITIALIZER;
5604 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5606 for (size_t x = 0; x < m->size2; x++)
5608 double f = gsl_matrix_get (m, y, x);
5610 ? data_out (&(union value) { .f = f / scale}, NULL,
5611 format, settings_get_fmt_settings ())
5612 : trimmed_string (f));
5613 ds_put_format (&line, " %s", s);
5616 output_log_nocopy (ds_steal_cstr (&line));
5619 string_array_destroy (rlabels);
5621 string_array_destroy (clabels);
5626 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5627 const struct print_labels *print_labels, size_t n,
5628 const char *name, const char *prefix)
5630 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5632 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5633 for (size_t i = 0; i < n; i++)
5634 pivot_category_create_leaf (
5636 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5637 : pivot_value_new_integer (i + 1)));
5639 d->hide_all_labels = true;
5640 string_array_destroy (labels);
5645 matrix_print_tables (const struct matrix_print *print, const gsl_matrix *m,
5646 struct fmt_spec format, int log_scale)
5648 struct pivot_table *table = pivot_table_create__ (
5649 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5652 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5654 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5655 N_("Columns"), "col");
5657 struct pivot_footnote *footnote = NULL;
5660 char *s = xasprintf ("× 10**%d", log_scale);
5661 footnote = pivot_table_create_footnote (
5662 table, pivot_value_new_user_text_nocopy (s));
5665 double scale = pow (10.0, log_scale);
5666 bool numeric = fmt_is_numeric (format.type);
5667 for (size_t y = 0; y < m->size1; y++)
5668 for (size_t x = 0; x < m->size2; x++)
5670 double f = gsl_matrix_get (m, y, x);
5671 struct pivot_value *v;
5674 v = pivot_value_new_number (f / scale);
5675 v->numeric.format = format;
5678 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5680 pivot_value_add_footnote (v, footnote);
5681 pivot_table_put2 (table, y, x, v);
5684 pivot_table_submit (table);
5688 matrix_print_execute (const struct matrix_print *print)
5690 if (print->expression)
5692 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5697 struct fmt_spec format = (print->use_default_format
5698 ? default_format (m, &log_scale)
5701 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5702 matrix_print_text (print, m, format, log_scale);
5704 matrix_print_tables (print, m, format, log_scale);
5706 gsl_matrix_free (m);
5710 matrix_print_space (print->space);
5712 output_log ("%s", print->title);
5719 matrix_do_if_clause_parse (struct matrix_state *s,
5720 struct matrix_do_if *ifc,
5721 bool parse_condition,
5722 size_t *allocated_clauses)
5724 if (ifc->n_clauses >= *allocated_clauses)
5725 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5726 sizeof *ifc->clauses);
5727 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5728 *c = (struct do_if_clause) { .condition = NULL };
5730 if (parse_condition)
5732 c->condition = matrix_expr_parse (s);
5737 return matrix_commands_parse (s, &c->commands, "DO IF", "ELSE", "END IF");
5740 static struct matrix_command *
5741 matrix_do_if_parse (struct matrix_state *s)
5743 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5744 *cmd = (struct matrix_command) {
5746 .do_if = { .n_clauses = 0 }
5749 size_t allocated_clauses = 0;
5752 if (!matrix_do_if_clause_parse (s, &cmd->do_if, true, &allocated_clauses))
5755 while (lex_match_phrase (s->lexer, "ELSE IF"));
5757 if (lex_match_id (s->lexer, "ELSE")
5758 && !matrix_do_if_clause_parse (s, &cmd->do_if, false, &allocated_clauses))
5761 if (!lex_match_phrase (s->lexer, "END IF"))
5766 matrix_command_destroy (cmd);
5771 matrix_do_if_execute (struct matrix_do_if *cmd)
5773 for (size_t i = 0; i < cmd->n_clauses; i++)
5775 struct do_if_clause *c = &cmd->clauses[i];
5779 if (!matrix_expr_evaluate_scalar (c->condition,
5780 i ? "ELSE IF" : "DO IF",
5785 for (size_t j = 0; j < c->commands.n; j++)
5786 if (!matrix_command_execute (c->commands.commands[j]))
5795 static struct matrix_command *
5796 matrix_loop_parse (struct matrix_state *s)
5798 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5799 *cmd = (struct matrix_command) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5801 struct matrix_loop *loop = &cmd->loop;
5802 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5804 struct substring name = lex_tokss (s->lexer);
5805 loop->var = matrix_var_lookup (s, name);
5807 loop->var = matrix_var_create (s, name);
5812 loop->start = matrix_expr_parse (s);
5813 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5816 loop->end = matrix_expr_parse (s);
5820 if (lex_match (s->lexer, T_BY))
5822 loop->increment = matrix_expr_parse (s);
5823 if (!loop->increment)
5828 if (lex_match_id (s->lexer, "IF"))
5830 loop->top_condition = matrix_expr_parse (s);
5831 if (!loop->top_condition)
5835 bool was_in_loop = s->in_loop;
5837 bool ok = matrix_commands_parse (s, &loop->commands, "LOOP",
5839 s->in_loop = was_in_loop;
5843 if (!lex_match_phrase (s->lexer, "END LOOP"))
5846 if (lex_match_id (s->lexer, "IF"))
5848 loop->bottom_condition = matrix_expr_parse (s);
5849 if (!loop->bottom_condition)
5856 matrix_command_destroy (cmd);
5861 matrix_loop_execute (struct matrix_loop *cmd)
5865 long int increment = 1;
5868 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5869 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5871 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5875 if (increment > 0 ? value > end
5876 : increment < 0 ? value < end
5881 int mxloops = settings_get_mxloops ();
5882 for (int i = 0; i < mxloops; i++)
5887 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5889 gsl_matrix_free (cmd->var->value);
5890 cmd->var->value = NULL;
5892 if (!cmd->var->value)
5893 cmd->var->value = gsl_matrix_alloc (1, 1);
5895 gsl_matrix_set (cmd->var->value, 0, 0, value);
5898 if (cmd->top_condition)
5901 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5906 for (size_t j = 0; j < cmd->commands.n; j++)
5907 if (!matrix_command_execute (cmd->commands.commands[j]))
5910 if (cmd->bottom_condition)
5913 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5914 "END LOOP IF", &d) || d > 0)
5921 if (increment > 0 ? value > end : value < end)
5929 static struct matrix_command *
5930 matrix_break_parse (struct matrix_state *s)
5934 lex_next_error (s->lexer, -1, -1, _("BREAK not inside LOOP."));
5938 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5939 *cmd = (struct matrix_command) { .type = MCMD_BREAK };
5945 static struct matrix_command *
5946 matrix_display_parse (struct matrix_state *s)
5948 while (lex_token (s->lexer) != T_ENDCMD)
5950 if (!lex_match_id (s->lexer, "DICTIONARY")
5951 && !lex_match_id (s->lexer, "STATUS"))
5953 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5958 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5959 *cmd = (struct matrix_command) { .type = MCMD_DISPLAY, .display = { s } };
5964 compare_matrix_var_pointers (const void *a_, const void *b_)
5966 const struct matrix_var *const *ap = a_;
5967 const struct matrix_var *const *bp = b_;
5968 const struct matrix_var *a = *ap;
5969 const struct matrix_var *b = *bp;
5970 return strcmp (a->name, b->name);
5974 matrix_display_execute (struct matrix_display *cmd)
5976 const struct matrix_state *s = cmd->state;
5978 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5979 struct pivot_dimension *attr_dimension
5980 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5981 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5982 N_("Rows"), N_("Columns"));
5983 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5985 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5988 const struct matrix_var *var;
5989 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5991 vars[n_vars++] = var;
5992 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5994 struct pivot_dimension *rows = pivot_dimension_create (
5995 table, PIVOT_AXIS_ROW, N_("Variable"));
5996 for (size_t i = 0; i < n_vars; i++)
5998 const struct matrix_var *var = vars[i];
5999 pivot_category_create_leaf (
6000 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
6002 size_t r = var->value->size1;
6003 size_t c = var->value->size2;
6004 double values[] = { r, c, r * c * sizeof (double) / 1024 };
6005 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
6006 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
6009 pivot_table_submit (table);
6014 static struct matrix_command *
6015 matrix_release_parse (struct matrix_state *s)
6017 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6018 *cmd = (struct matrix_command) { .type = MCMD_RELEASE };
6020 size_t allocated_vars = 0;
6021 while (lex_token (s->lexer) == T_ID)
6023 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
6026 if (cmd->release.n_vars >= allocated_vars)
6027 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
6028 sizeof *cmd->release.vars);
6029 cmd->release.vars[cmd->release.n_vars++] = var;
6032 lex_error (s->lexer, _("Syntax error expecting variable name."));
6035 if (!lex_match (s->lexer, T_COMMA))
6043 matrix_release_execute (struct matrix_release *cmd)
6045 for (size_t i = 0; i < cmd->n_vars; i++)
6047 struct matrix_var *var = cmd->vars[i];
6048 gsl_matrix_free (var->value);
6055 static struct save_file *
6056 save_file_create (struct matrix_state *s, struct file_handle *fh,
6057 struct string_array *variables,
6058 struct matrix_expr *names,
6059 struct stringi_set *strings)
6061 for (size_t i = 0; i < s->n_save_files; i++)
6063 struct save_file *sf = s->save_files[i];
6064 if (fh_equal (sf->file, fh))
6068 string_array_destroy (variables);
6069 matrix_expr_destroy (names);
6070 stringi_set_destroy (strings);
6076 struct save_file *sf = xmalloc (sizeof *sf);
6077 *sf = (struct save_file) {
6079 .dataset = s->dataset,
6080 .variables = *variables,
6082 .strings = STRINGI_SET_INITIALIZER (sf->strings),
6085 stringi_set_swap (&sf->strings, strings);
6086 stringi_set_destroy (strings);
6088 s->save_files = xrealloc (s->save_files,
6089 (s->n_save_files + 1) * sizeof *s->save_files);
6090 s->save_files[s->n_save_files++] = sf;
6094 static struct casewriter *
6095 save_file_open (struct save_file *sf, gsl_matrix *m,
6096 const struct msg_location *save_location)
6098 if (sf->writer || sf->error)
6102 size_t n_variables = caseproto_get_n_widths (
6103 casewriter_get_proto (sf->writer));
6104 if (m->size2 != n_variables)
6106 const char *file_name = (sf->file == fh_inline_file () ? "*"
6107 : fh_get_name (sf->file));
6108 msg_at (SE, save_location,
6109 _("Cannot save %zu×%zu matrix to %s because the "
6110 "first SAVE to %s in this matrix program wrote a "
6111 "%zu-column matrix."),
6112 m->size1, m->size2, file_name, file_name, n_variables);
6113 msg_at (SE, sf->location,
6114 _("This is the location of the first SAVE to %s."),
6123 struct dictionary *dict = dict_create (get_default_encoding ());
6125 /* Fill 'names' with user-specified names if there were any, either from
6126 sf->variables or sf->names. */
6127 struct string_array names = { .n = 0 };
6128 if (sf->variables.n)
6129 string_array_clone (&names, &sf->variables);
6132 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
6133 if (nm && is_vector (nm))
6135 gsl_vector nv = to_vector (nm);
6136 for (size_t i = 0; i < nv.size; i++)
6138 char *name = trimmed_string (gsl_vector_get (&nv, i));
6139 char *error = dict_id_is_valid__ (dict, name);
6141 string_array_append_nocopy (&names, name);
6144 msg_at (SE, save_location, "%s", error);
6150 gsl_matrix_free (nm);
6153 struct stringi_set strings;
6154 stringi_set_clone (&strings, &sf->strings);
6156 for (size_t i = 0; dict_get_n_vars (dict) < m->size2; i++)
6161 name = names.strings[i];
6164 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
6168 int width = stringi_set_delete (&strings, name) ? 8 : 0;
6169 struct variable *var = dict_create_var (dict, name, width);
6172 msg_at (ME, save_location,
6173 _("Duplicate variable name %s in SAVE statement."), name);
6178 if (!stringi_set_is_empty (&strings))
6180 size_t count = stringi_set_count (&strings);
6181 const char *example = stringi_set_node_get_string (
6182 stringi_set_first (&strings));
6185 msg_at (ME, save_location,
6186 _("The SAVE command STRINGS subcommand specifies an unknown "
6187 "variable %s."), example);
6189 msg_at (ME, save_location,
6190 ngettext ("The SAVE command STRINGS subcommand specifies %zu "
6191 "unknown variable: %s.",
6192 "The SAVE command STRINGS subcommand specifies %zu "
6193 "unknown variables, including %s.",
6198 stringi_set_destroy (&strings);
6199 string_array_destroy (&names);
6208 if (sf->file == fh_inline_file ())
6209 sf->writer = autopaging_writer_create (dict_get_proto (dict));
6211 sf->writer = any_writer_open (sf->file, dict);
6215 sf->location = msg_location_dup (save_location);
6226 save_file_destroy (struct save_file *sf)
6230 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
6232 dataset_set_dict (sf->dataset, sf->dict);
6233 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
6237 casewriter_destroy (sf->writer);
6238 dict_unref (sf->dict);
6240 fh_unref (sf->file);
6241 string_array_destroy (&sf->variables);
6242 matrix_expr_destroy (sf->names);
6243 stringi_set_destroy (&sf->strings);
6244 msg_location_destroy (sf->location);
6249 static struct matrix_command *
6250 matrix_save_parse (struct matrix_state *s)
6252 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6253 *cmd = (struct matrix_command) {
6255 .save = { .expression = NULL },
6258 struct file_handle *fh = NULL;
6259 struct matrix_save *save = &cmd->save;
6261 struct string_array variables = STRING_ARRAY_INITIALIZER;
6262 struct matrix_expr *names = NULL;
6263 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
6265 save->expression = matrix_parse_exp (s);
6266 if (!save->expression)
6269 int names_start = 0;
6271 while (lex_match (s->lexer, T_SLASH))
6273 if (lex_match_id (s->lexer, "OUTFILE"))
6275 lex_match (s->lexer, T_EQUALS);
6278 fh = (lex_match (s->lexer, T_ASTERISK)
6280 : fh_parse (s->lexer, FH_REF_FILE, s->session));
6284 else if (lex_match_id (s->lexer, "VARIABLES"))
6286 lex_match (s->lexer, T_EQUALS);
6290 struct dictionary *d = dict_create (get_default_encoding ());
6291 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
6292 PV_NO_SCRATCH | PV_NO_DUPLICATE);
6297 string_array_clear (&variables);
6298 variables = (struct string_array) {
6304 else if (lex_match_id (s->lexer, "NAMES"))
6306 lex_match (s->lexer, T_EQUALS);
6307 matrix_expr_destroy (names);
6308 names_start = lex_ofs (s->lexer);
6309 names = matrix_parse_exp (s);
6310 names_end = lex_ofs (s->lexer) - 1;
6314 else if (lex_match_id (s->lexer, "STRINGS"))
6316 lex_match (s->lexer, T_EQUALS);
6317 while (lex_token (s->lexer) == T_ID)
6319 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
6321 if (!lex_match (s->lexer, T_COMMA))
6327 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
6335 if (s->prev_save_file)
6336 fh = fh_ref (s->prev_save_file);
6339 lex_sbc_missing (s->lexer, "OUTFILE");
6343 fh_unref (s->prev_save_file);
6344 s->prev_save_file = fh_ref (fh);
6346 if (variables.n && names)
6348 lex_ofs_msg (s->lexer, SW, names_start, names_end,
6349 _("Ignoring NAMES because VARIABLES was also specified."));
6350 matrix_expr_destroy (names);
6354 save->sf = save_file_create (s, fh, &variables, names, &strings);
6358 string_array_destroy (&variables);
6359 matrix_expr_destroy (names);
6360 stringi_set_destroy (&strings);
6362 matrix_command_destroy (cmd);
6367 matrix_save_execute (const struct matrix_command *cmd)
6369 const struct matrix_save *save = &cmd->save;
6371 gsl_matrix *m = matrix_expr_evaluate (save->expression);
6375 struct casewriter *writer = save_file_open (save->sf, m, cmd->location);
6378 gsl_matrix_free (m);
6382 const struct caseproto *proto = casewriter_get_proto (writer);
6383 for (size_t y = 0; y < m->size1; y++)
6385 struct ccase *c = case_create (proto);
6386 for (size_t x = 0; x < m->size2; x++)
6388 double d = gsl_matrix_get (m, y, x);
6389 int width = caseproto_get_width (proto, x);
6390 union value *value = case_data_rw_idx (c, x);
6394 memcpy (value->s, &d, width);
6396 casewriter_write (writer, c);
6398 gsl_matrix_free (m);
6403 static struct read_file *
6404 read_file_create (struct matrix_state *s, struct file_handle *fh)
6406 for (size_t i = 0; i < s->n_read_files; i++)
6408 struct read_file *rf = s->read_files[i];
6416 struct read_file *rf = xmalloc (sizeof *rf);
6417 *rf = (struct read_file) { .file = fh };
6419 s->read_files = xrealloc (s->read_files,
6420 (s->n_read_files + 1) * sizeof *s->read_files);
6421 s->read_files[s->n_read_files++] = rf;
6425 static struct dfm_reader *
6426 read_file_open (struct read_file *rf)
6429 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
6434 read_file_destroy (struct read_file *rf)
6438 fh_unref (rf->file);
6439 dfm_close_reader (rf->reader);
6440 free (rf->encoding);
6445 static struct matrix_command *
6446 matrix_read_parse (struct matrix_state *s)
6448 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6449 *cmd = (struct matrix_command) {
6451 .read = { .format = FMT_F },
6454 struct file_handle *fh = NULL;
6455 char *encoding = NULL;
6456 struct matrix_read *read = &cmd->read;
6457 read->dst = matrix_lvalue_parse (s);
6463 int record_width_start = 0, record_width_end = 0;
6466 int repetitions = 0;
6467 int record_width = 0;
6468 bool seen_format = false;
6469 while (lex_match (s->lexer, T_SLASH))
6471 if (lex_match_id (s->lexer, "FILE"))
6473 lex_match (s->lexer, T_EQUALS);
6476 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6480 else if (lex_match_id (s->lexer, "ENCODING"))
6482 lex_match (s->lexer, T_EQUALS);
6483 if (!lex_force_string (s->lexer))
6487 encoding = ss_xstrdup (lex_tokss (s->lexer));
6491 else if (lex_match_id (s->lexer, "FIELD"))
6493 lex_match (s->lexer, T_EQUALS);
6495 record_width_start = lex_ofs (s->lexer);
6496 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6498 read->c1 = lex_integer (s->lexer);
6500 if (!lex_force_match (s->lexer, T_TO)
6501 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
6503 read->c2 = lex_integer (s->lexer) + 1;
6504 record_width_end = lex_ofs (s->lexer);
6507 record_width = read->c2 - read->c1;
6508 if (lex_match (s->lexer, T_BY))
6510 if (!lex_force_int_range (s->lexer, "BY", 1,
6511 read->c2 - read->c1))
6513 by = lex_integer (s->lexer);
6514 by_ofs = lex_ofs (s->lexer);
6515 int field_end = lex_ofs (s->lexer);
6518 if (record_width % by)
6521 s->lexer, record_width_start, field_end,
6522 _("Field width %d does not evenly divide record width %d."),
6524 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
6525 _("This syntax designates the record width."));
6526 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
6527 _("This syntax specifies the field width."));
6534 else if (lex_match_id (s->lexer, "SIZE"))
6536 lex_match (s->lexer, T_EQUALS);
6537 matrix_expr_destroy (read->size);
6538 read->size = matrix_parse_exp (s);
6542 else if (lex_match_id (s->lexer, "MODE"))
6544 lex_match (s->lexer, T_EQUALS);
6545 if (lex_match_id (s->lexer, "RECTANGULAR"))
6546 read->symmetric = false;
6547 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6548 read->symmetric = true;
6551 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6555 else if (lex_match_id (s->lexer, "REREAD"))
6556 read->reread = true;
6557 else if (lex_match_id (s->lexer, "FORMAT"))
6561 lex_sbc_only_once (s->lexer, "FORMAT");
6566 lex_match (s->lexer, T_EQUALS);
6568 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6571 format_ofs = lex_ofs (s->lexer);
6572 const char *p = lex_tokcstr (s->lexer);
6573 if (c_isdigit (p[0]))
6575 repetitions = atoi (p);
6576 p += strspn (p, "0123456789");
6577 if (!fmt_from_name (p, &read->format))
6579 lex_error (s->lexer, _("Unknown format %s."), p);
6584 else if (fmt_from_name (p, &read->format))
6588 struct fmt_spec format;
6589 if (!parse_format_specifier (s->lexer, &format))
6591 read->format = format.type;
6597 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6598 "REREAD", "FORMAT");
6605 lex_sbc_missing (s->lexer, "FIELD");
6609 if (!read->dst->n_indexes && !read->size)
6611 msg (SE, _("SIZE is required for reading data into a full matrix "
6612 "(as opposed to a submatrix)."));
6613 msg_at (SN, read->dst->var_location,
6614 _("This expression designates a full matrix."));
6620 if (s->prev_read_file)
6621 fh = fh_ref (s->prev_read_file);
6624 lex_sbc_missing (s->lexer, "FILE");
6628 fh_unref (s->prev_read_file);
6629 s->prev_read_file = fh_ref (fh);
6631 read->rf = read_file_create (s, fh);
6635 free (read->rf->encoding);
6636 read->rf->encoding = encoding;
6640 /* Field width may be specified in multiple ways:
6643 2. The format on FORMAT.
6644 3. The repetition factor on FORMAT.
6646 (2) and (3) are mutually exclusive.
6648 If more than one of these is present, they must agree. If none of them is
6649 present, then free-field format is used.
6651 if (repetitions > record_width)
6653 msg (SE, _("%d repetitions cannot fit in record width %d."),
6654 repetitions, record_width);
6655 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
6656 _("This syntax designates the number of repetitions."));
6657 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
6658 _("This syntax designates the record width."));
6661 int w = (repetitions ? record_width / repetitions
6666 msg (SE, _("This command specifies two different field widths."));
6669 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
6670 ngettext ("This syntax specifies %d repetition.",
6671 "This syntax specifies %d repetitions.",
6674 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
6675 _("This syntax designates record width %d, "
6676 "which divided by %d repetitions implies "
6678 record_width, repetitions, w);
6681 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
6682 _("This syntax specifies field width %d."), w);
6684 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
6685 _("This syntax specifies field width %d."), by);
6693 matrix_command_destroy (cmd);
6699 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6700 struct substring data, size_t y, size_t x,
6701 int first_column, int last_column, char *error)
6703 int line_number = dfm_get_line_number (reader);
6704 struct msg_location location = {
6705 .file_name = intern_new (dfm_get_file_name (reader)),
6706 .start = { .line = line_number, .column = first_column },
6707 .end = { .line = line_number, .column = last_column },
6709 msg_at (DW, &location, _("Error reading \"%.*s\" as format %s "
6710 "for matrix row %zu, column %zu: %s"),
6711 (int) data.length, data.string, fmt_name (format),
6712 y + 1, x + 1, error);
6713 msg_location_uninit (&location);
6718 matrix_read_set_field (struct matrix_read *read, struct dfm_reader *reader,
6719 gsl_matrix *m, struct substring p, size_t y, size_t x,
6720 const char *line_start)
6722 const char *input_encoding = dfm_reader_get_encoding (reader);
6725 if (fmt_is_numeric (read->format))
6728 error = data_in (p, input_encoding, read->format,
6729 settings_get_fmt_settings (), &v, 0, NULL);
6730 if (!error && v.f == SYSMIS)
6731 error = xstrdup (_("Matrix data may not contain missing value."));
6736 uint8_t s[sizeof (double)];
6737 union value v = { .s = s };
6738 error = data_in (p, input_encoding, read->format,
6739 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6740 memcpy (&f, s, sizeof f);
6745 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6746 int nc = ss_utf8_count_columns (p);
6747 int c2 = c1 + MAX (1, nc) - 1;
6748 parse_error (reader, read->format, p, y, x, c1, c2, error);
6752 gsl_matrix_set (m, y, x, f);
6753 if (read->symmetric && x != y)
6754 gsl_matrix_set (m, x, y, f);
6759 matrix_read_line (struct matrix_command *cmd, struct dfm_reader *reader,
6760 struct substring *line, const char **startp)
6762 struct matrix_read *read = &cmd->read;
6763 if (dfm_eof (reader))
6765 msg_at (SE, cmd->location,
6766 _("Unexpected end of file reading matrix data."));
6769 dfm_expand_tabs (reader);
6770 struct substring record = dfm_get_record (reader);
6771 /* XXX need to recode record into UTF-8 */
6772 *startp = record.string;
6773 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6778 matrix_read (struct matrix_command *cmd, struct dfm_reader *reader,
6781 struct matrix_read *read = &cmd->read;
6782 for (size_t y = 0; y < m->size1; y++)
6784 size_t nx = read->symmetric ? y + 1 : m->size2;
6786 struct substring line = ss_empty ();
6787 const char *line_start = line.string;
6788 for (size_t x = 0; x < nx; x++)
6795 ss_ltrim (&line, ss_cstr (" ,"));
6796 if (!ss_is_empty (line))
6798 if (!matrix_read_line (cmd, reader, &line, &line_start))
6800 dfm_forward_record (reader);
6803 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6807 if (!matrix_read_line (cmd, reader, &line, &line_start))
6809 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6810 int f = x % fields_per_line;
6811 if (f == fields_per_line - 1)
6812 dfm_forward_record (reader);
6814 p = ss_substr (line, read->w * f, read->w);
6817 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6821 dfm_forward_record (reader);
6824 ss_ltrim (&line, ss_cstr (" ,"));
6825 if (!ss_is_empty (line))
6827 int line_number = dfm_get_line_number (reader);
6828 int c1 = utf8_count_columns (line_start,
6829 line.string - line_start) + 1;
6830 int c2 = c1 + ss_utf8_count_columns (line) - 1;
6831 struct msg_location location = {
6832 .file_name = intern_new (dfm_get_file_name (reader)),
6833 .start = { .line = line_number, .column = c1 },
6834 .end = { .line = line_number, .column = c2 },
6836 msg_at (DW, &location,
6837 _("Trailing garbage following data for matrix row %zu."),
6839 msg_location_uninit (&location);
6846 matrix_read_execute (struct matrix_command *cmd)
6848 struct matrix_read *read = &cmd->read;
6849 struct index_vector iv0, iv1;
6850 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6853 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6856 gsl_matrix *m = matrix_expr_evaluate (read->size);
6862 msg_at (SE, matrix_expr_location (read->size),
6863 _("SIZE must evaluate to a scalar or a 2-element vector, "
6864 "not a %zu×%zu matrix."), m->size1, m->size2);
6865 gsl_matrix_free (m);
6866 index_vector_uninit (&iv0);
6867 index_vector_uninit (&iv1);
6871 gsl_vector v = to_vector (m);
6875 d[0] = gsl_vector_get (&v, 0);
6878 else if (v.size == 2)
6880 d[0] = gsl_vector_get (&v, 0);
6881 d[1] = gsl_vector_get (&v, 1);
6885 msg_at (SE, matrix_expr_location (read->size),
6886 _("SIZE must evaluate to a scalar or a 2-element vector, "
6887 "not a %zu×%zu matrix."),
6888 m->size1, m->size2),
6889 gsl_matrix_free (m);
6890 index_vector_uninit (&iv0);
6891 index_vector_uninit (&iv1);
6894 gsl_matrix_free (m);
6896 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6898 msg_at (SE, matrix_expr_location (read->size),
6899 _("Matrix dimensions %g×%g specified on SIZE "
6900 "are outside valid range."),
6902 index_vector_uninit (&iv0);
6903 index_vector_uninit (&iv1);
6911 if (read->dst->n_indexes)
6913 size_t submatrix_size[2];
6914 if (read->dst->n_indexes == 2)
6916 submatrix_size[0] = iv0.n;
6917 submatrix_size[1] = iv1.n;
6919 else if (read->dst->var->value->size1 == 1)
6921 submatrix_size[0] = 1;
6922 submatrix_size[1] = iv0.n;
6926 submatrix_size[0] = iv0.n;
6927 submatrix_size[1] = 1;
6932 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6934 msg_at (SE, cmd->location,
6935 _("Dimensions specified on SIZE differ from dimensions "
6936 "of destination submatrix."));
6937 msg_at (SN, matrix_expr_location (read->size),
6938 _("SIZE specifies dimensions %zu×%zu."),
6940 msg_at (SN, read->dst->full_location,
6941 _("Destination submatrix has dimensions %zu×%zu."),
6942 submatrix_size[0], submatrix_size[1]);
6943 index_vector_uninit (&iv0);
6944 index_vector_uninit (&iv1);
6950 size[0] = submatrix_size[0];
6951 size[1] = submatrix_size[1];
6955 struct dfm_reader *reader = read_file_open (read->rf);
6957 dfm_reread_record (reader, 1);
6959 if (read->symmetric && size[0] != size[1])
6961 msg_at (SE, cmd->location,
6962 _("Cannot read non-square %zu×%zu matrix "
6963 "using READ with MODE=SYMMETRIC."),
6965 index_vector_uninit (&iv0);
6966 index_vector_uninit (&iv1);
6969 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6970 matrix_read (cmd, reader, tmp);
6971 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp, cmd->location);
6976 static struct write_file *
6977 write_file_create (struct matrix_state *s, struct file_handle *fh)
6979 for (size_t i = 0; i < s->n_write_files; i++)
6981 struct write_file *wf = s->write_files[i];
6989 struct write_file *wf = xmalloc (sizeof *wf);
6990 *wf = (struct write_file) { .file = fh };
6992 s->write_files = xrealloc (s->write_files,
6993 (s->n_write_files + 1) * sizeof *s->write_files);
6994 s->write_files[s->n_write_files++] = wf;
6998 static struct dfm_writer *
6999 write_file_open (struct write_file *wf)
7002 wf->writer = dfm_open_writer (wf->file, wf->encoding);
7007 write_file_destroy (struct write_file *wf)
7013 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
7014 wf->held->s.ss.length);
7015 u8_line_destroy (wf->held);
7019 fh_unref (wf->file);
7020 dfm_close_writer (wf->writer);
7021 free (wf->encoding);
7026 static struct matrix_command *
7027 matrix_write_parse (struct matrix_state *s)
7029 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7030 *cmd = (struct matrix_command) {
7034 struct file_handle *fh = NULL;
7035 char *encoding = NULL;
7036 struct matrix_write *write = &cmd->write;
7037 write->expression = matrix_parse_exp (s);
7038 if (!write->expression)
7043 int record_width_start = 0, record_width_end = 0;
7046 int repetitions = 0;
7047 int record_width = 0;
7048 enum fmt_type format = FMT_F;
7049 bool has_format = false;
7050 while (lex_match (s->lexer, T_SLASH))
7052 if (lex_match_id (s->lexer, "OUTFILE"))
7054 lex_match (s->lexer, T_EQUALS);
7057 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
7061 else if (lex_match_id (s->lexer, "ENCODING"))
7063 lex_match (s->lexer, T_EQUALS);
7064 if (!lex_force_string (s->lexer))
7068 encoding = ss_xstrdup (lex_tokss (s->lexer));
7072 else if (lex_match_id (s->lexer, "FIELD"))
7074 lex_match (s->lexer, T_EQUALS);
7076 record_width_start = lex_ofs (s->lexer);
7078 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
7080 write->c1 = lex_integer (s->lexer);
7082 if (!lex_force_match (s->lexer, T_TO)
7083 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
7085 write->c2 = lex_integer (s->lexer) + 1;
7086 record_width_end = lex_ofs (s->lexer);
7089 record_width = write->c2 - write->c1;
7090 if (lex_match (s->lexer, T_BY))
7092 if (!lex_force_int_range (s->lexer, "BY", 1,
7093 write->c2 - write->c1))
7095 by_ofs = lex_ofs (s->lexer);
7096 int field_end = lex_ofs (s->lexer);
7097 by = lex_integer (s->lexer);
7100 if (record_width % by)
7103 s->lexer, record_width_start, field_end,
7104 _("Field width %d does not evenly divide record width %d."),
7106 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7107 _("This syntax designates the record width."));
7108 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
7109 _("This syntax specifies the field width."));
7116 else if (lex_match_id (s->lexer, "MODE"))
7118 lex_match (s->lexer, T_EQUALS);
7119 if (lex_match_id (s->lexer, "RECTANGULAR"))
7120 write->triangular = false;
7121 else if (lex_match_id (s->lexer, "TRIANGULAR"))
7122 write->triangular = true;
7125 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
7129 else if (lex_match_id (s->lexer, "HOLD"))
7131 else if (lex_match_id (s->lexer, "FORMAT"))
7133 if (has_format || write->format)
7135 lex_sbc_only_once (s->lexer, "FORMAT");
7139 lex_match (s->lexer, T_EQUALS);
7141 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
7144 format_ofs = lex_ofs (s->lexer);
7145 const char *p = lex_tokcstr (s->lexer);
7146 if (c_isdigit (p[0]))
7148 repetitions = atoi (p);
7149 p += strspn (p, "0123456789");
7150 if (!fmt_from_name (p, &format))
7152 lex_error (s->lexer, _("Unknown format %s."), p);
7158 else if (fmt_from_name (p, &format))
7165 struct fmt_spec spec;
7166 if (!parse_format_specifier (s->lexer, &spec))
7168 write->format = xmemdup (&spec, sizeof spec);
7173 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
7181 lex_sbc_missing (s->lexer, "FIELD");
7187 if (s->prev_write_file)
7188 fh = fh_ref (s->prev_write_file);
7191 lex_sbc_missing (s->lexer, "OUTFILE");
7195 fh_unref (s->prev_write_file);
7196 s->prev_write_file = fh_ref (fh);
7198 write->wf = write_file_create (s, fh);
7202 free (write->wf->encoding);
7203 write->wf->encoding = encoding;
7207 /* Field width may be specified in multiple ways:
7210 2. The format on FORMAT.
7211 3. The repetition factor on FORMAT.
7213 (2) and (3) are mutually exclusive.
7215 If more than one of these is present, they must agree. If none of them is
7216 present, then free-field format is used.
7218 if (repetitions > record_width)
7220 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7221 _("This syntax designates the number of repetitions."));
7222 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7223 _("This syntax designates the record width."));
7226 int w = (repetitions ? record_width / repetitions
7227 : write->format ? write->format->w
7231 msg (SE, _("This command specifies two different field widths."));
7234 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7235 ngettext ("This syntax specifies %d repetition.",
7236 "This syntax specifies %d repetitions.",
7239 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7240 _("This syntax designates record width %d, "
7241 "which divided by %d repetitions implies "
7243 record_width, repetitions, w);
7246 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7247 _("This syntax specifies field width %d."), w);
7249 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
7250 _("This syntax specifies field width %d."), by);
7253 if (w && !write->format)
7255 write->format = xmalloc (sizeof *write->format);
7256 *write->format = (struct fmt_spec) { .type = format, .w = w };
7258 char *error = fmt_check_output__ (*write->format);
7261 msg (SE, "%s", error);
7265 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7266 _("This syntax specifies format %s."),
7271 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7272 ngettext ("This syntax specifies %d repetition.",
7273 "This syntax specifies %d repetitions.",
7276 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7277 _("This syntax designates record width %d, "
7278 "which divided by %d repetitions implies "
7280 record_width, repetitions, w);
7284 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
7285 _("This syntax specifies field width %d."), by);
7291 if (write->format && fmt_var_width (*write->format) > sizeof (double))
7293 char fs[FMT_STRING_LEN_MAX + 1];
7294 fmt_to_string (*write->format, fs);
7295 lex_ofs_error (s->lexer, format_ofs, format_ofs,
7296 _("Format %s is too wide for %zu-byte matrix elements."),
7297 fs, sizeof (double));
7305 matrix_command_destroy (cmd);
7310 matrix_write_execute (struct matrix_write *write)
7312 gsl_matrix *m = matrix_expr_evaluate (write->expression);
7316 if (write->triangular && m->size1 != m->size2)
7318 msg_at (SE, matrix_expr_location (write->expression),
7319 _("WRITE with MODE=TRIANGULAR requires a square matrix but "
7320 "the matrix to be written has dimensions %zu×%zu."),
7321 m->size1, m->size2);
7322 gsl_matrix_free (m);
7326 struct dfm_writer *writer = write_file_open (write->wf);
7327 if (!writer || !m->size1)
7329 gsl_matrix_free (m);
7333 const struct fmt_settings *settings = settings_get_fmt_settings ();
7334 struct u8_line *line = write->wf->held;
7335 for (size_t y = 0; y < m->size1; y++)
7339 line = xmalloc (sizeof *line);
7340 u8_line_init (line);
7342 size_t nx = write->triangular ? y + 1 : m->size2;
7344 for (size_t x = 0; x < nx; x++)
7347 double f = gsl_matrix_get (m, y, x);
7351 if (fmt_is_numeric (write->format->type))
7354 v.s = (uint8_t *) &f;
7355 s = data_out (&v, NULL, *write->format, settings);
7359 s = xmalloc (DBL_BUFSIZE_BOUND);
7360 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
7361 >= DBL_BUFSIZE_BOUND)
7364 size_t len = strlen (s);
7365 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
7366 if (width + x0 > write->c2)
7368 dfm_put_record_utf8 (writer, line->s.ss.string,
7370 u8_line_clear (line);
7373 u8_line_put (line, x0, x0 + width, s, len);
7376 x0 += write->format ? write->format->w : width + 1;
7379 if (y + 1 >= m->size1 && write->hold)
7381 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
7382 u8_line_clear (line);
7386 u8_line_destroy (line);
7390 write->wf->held = line;
7392 gsl_matrix_free (m);
7397 static struct matrix_command *
7398 matrix_get_parse (struct matrix_state *s)
7400 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7401 *cmd = (struct matrix_command) {
7405 .dataset = s->dataset,
7406 .user = { .treatment = MGET_ERROR },
7407 .system = { .treatment = MGET_ERROR },
7411 struct matrix_get *get = &cmd->get;
7412 get->dst = matrix_lvalue_parse (s);
7416 while (lex_match (s->lexer, T_SLASH))
7418 if (lex_match_id (s->lexer, "FILE"))
7420 lex_match (s->lexer, T_EQUALS);
7422 fh_unref (get->file);
7423 if (lex_match (s->lexer, T_ASTERISK))
7427 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7432 else if (lex_match_id (s->lexer, "ENCODING"))
7434 lex_match (s->lexer, T_EQUALS);
7435 if (!lex_force_string (s->lexer))
7438 free (get->encoding);
7439 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
7443 else if (lex_match_id (s->lexer, "VARIABLES"))
7445 lex_match (s->lexer, T_EQUALS);
7449 lex_sbc_only_once (s->lexer, "VARIABLES");
7453 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
7456 else if (lex_match_id (s->lexer, "NAMES"))
7458 lex_match (s->lexer, T_EQUALS);
7459 if (!lex_force_id (s->lexer))
7462 struct substring name = lex_tokss (s->lexer);
7463 get->names = matrix_var_lookup (s, name);
7465 get->names = matrix_var_create (s, name);
7468 else if (lex_match_id (s->lexer, "MISSING"))
7470 lex_match (s->lexer, T_EQUALS);
7471 if (lex_match_id (s->lexer, "ACCEPT"))
7472 get->user.treatment = MGET_ACCEPT;
7473 else if (lex_match_id (s->lexer, "OMIT"))
7474 get->user.treatment = MGET_OMIT;
7475 else if (lex_is_number (s->lexer))
7477 get->user.treatment = MGET_RECODE;
7478 get->user.substitute = lex_number (s->lexer);
7483 lex_error (s->lexer, _("Syntax error expecting ACCEPT or OMIT or "
7484 "a number for MISSING."));
7488 else if (lex_match_id (s->lexer, "SYSMIS"))
7490 lex_match (s->lexer, T_EQUALS);
7491 if (lex_match_id (s->lexer, "OMIT"))
7492 get->system.treatment = MGET_OMIT;
7493 else if (lex_is_number (s->lexer))
7495 get->system.treatment = MGET_RECODE;
7496 get->system.substitute = lex_number (s->lexer);
7501 lex_error (s->lexer, _("Syntax error expecting OMIT or a number "
7508 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
7509 "MISSING", "SYSMIS");
7514 if (get->user.treatment != MGET_ACCEPT)
7515 get->system.treatment = MGET_ERROR;
7520 matrix_command_destroy (cmd);
7525 matrix_get_execute__ (struct matrix_command *cmd, struct casereader *reader,
7526 const struct dictionary *dict)
7528 struct matrix_get *get = &cmd->get;
7529 struct variable **vars;
7534 if (!var_syntax_evaluate (get->lexer, get->vars, get->n_vars, dict,
7535 &vars, &n_vars, PV_NUMERIC))
7540 n_vars = dict_get_n_vars (dict);
7541 vars = xnmalloc (n_vars, sizeof *vars);
7542 for (size_t i = 0; i < n_vars; i++)
7544 struct variable *var = dict_get_var (dict, i);
7545 if (!var_is_numeric (var))
7547 msg_at (SE, cmd->location, _("Variable %s is not numeric."),
7548 var_get_name (var));
7558 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
7559 for (size_t i = 0; i < n_vars; i++)
7561 char s[sizeof (double)];
7563 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
7564 memcpy (&f, s, sizeof f);
7565 gsl_matrix_set (names, i, 0, f);
7568 gsl_matrix_free (get->names->value);
7569 get->names->value = names;
7573 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
7574 long long int casenum = 1;
7576 for (struct ccase *c = casereader_read (reader); c;
7577 c = casereader_read (reader), casenum++)
7579 if (n_rows >= m->size1)
7581 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
7582 for (size_t y = 0; y < n_rows; y++)
7583 for (size_t x = 0; x < n_vars; x++)
7584 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
7585 gsl_matrix_free (m);
7590 for (size_t x = 0; x < n_vars; x++)
7592 const struct variable *var = vars[x];
7593 double d = case_num (c, var);
7596 if (get->system.treatment == MGET_RECODE)
7597 d = get->system.substitute;
7598 else if (get->system.treatment == MGET_OMIT)
7602 msg_at (SE, cmd->location, _("Variable %s in case %lld "
7603 "is system-missing."),
7604 var_get_name (var), casenum);
7608 else if (var_is_num_missing (var, d) == MV_USER)
7610 if (get->user.treatment == MGET_RECODE)
7611 d = get->user.substitute;
7612 else if (get->user.treatment == MGET_OMIT)
7614 else if (get->user.treatment != MGET_ACCEPT)
7616 msg_at (SE, cmd->location,
7617 _("Variable %s in case %lld has user-missing "
7619 var_get_name (var), casenum, d);
7623 gsl_matrix_set (m, n_rows, x, d);
7634 matrix_lvalue_evaluate_and_assign (get->dst, m, cmd->location);
7637 gsl_matrix_free (m);
7642 matrix_open_casereader (const struct matrix_command *cmd,
7643 const char *command_name, struct file_handle *file,
7644 const char *encoding, struct dataset *dataset,
7645 struct casereader **readerp, struct dictionary **dictp)
7649 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
7650 return *readerp != NULL;
7654 if (dict_get_n_vars (dataset_dict (dataset)) == 0)
7656 msg_at (SE, cmd->location,
7657 _("The %s command cannot read an empty active file."),
7661 *readerp = proc_open (dataset);
7662 *dictp = dict_ref (dataset_dict (dataset));
7668 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
7669 struct casereader *reader, struct dictionary *dict)
7672 casereader_destroy (reader);
7674 proc_commit (dataset);
7678 matrix_get_execute (struct matrix_command *cmd)
7680 struct matrix_get *get = &cmd->get;
7681 struct casereader *r;
7682 struct dictionary *d;
7683 if (matrix_open_casereader (cmd, "GET", get->file, get->encoding,
7684 get->dataset, &r, &d))
7686 matrix_get_execute__ (cmd, r, d);
7687 matrix_close_casereader (get->file, get->dataset, r, d);
7694 variables_changed (const char *keyword,
7695 const struct string_array *new_vars,
7696 const struct msg_location *new_vars_location,
7697 const struct msg_location *new_location,
7698 const struct string_array *old_vars,
7699 const struct msg_location *old_vars_location,
7700 const struct msg_location *old_location)
7706 msg_at (SE, new_location,
7707 _("%s may only be specified on MSAVE if it was specified "
7708 "on the first MSAVE within MATRIX."), keyword);
7709 msg_at (SN, old_location,
7710 _("The first MSAVE in MATRIX did not specify %s."),
7712 msg_at (SN, new_vars_location,
7713 _("This is the specification of %s on a later MSAVE."),
7717 if (!string_array_equal_case (old_vars, new_vars))
7719 msg_at (SE, new_location,
7720 _("%s must specify the same variables on each MSAVE "
7721 "within a given MATRIX."), keyword);
7722 msg_at (SE, old_vars_location,
7723 _("This is the specification of %s on the first MSAVE."),
7725 msg_at (SE, new_vars_location,
7726 _("This is a different specification of %s on a later MSAVE."),
7735 msave_common_changed (const struct msave_common *old,
7736 const struct msave_common *new)
7738 if (new->outfile && !fh_equal (old->outfile, new->outfile))
7740 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7741 "within a single MATRIX command."));
7742 msg_at (SN, old->outfile_location,
7743 _("This is the OUTFILE on the first MSAVE command."));
7744 msg_at (SN, new->outfile_location,
7745 _("This is the OUTFILE on a later MSAVE command."));
7749 if (!variables_changed ("VARIABLES",
7750 &new->variables, new->variables_location, new->location,
7751 &old->variables, old->variables_location, old->location)
7752 && !variables_changed ("FNAMES",
7753 &new->fnames, new->fnames_location, new->location,
7754 &old->fnames, old->fnames_location, old->location)
7755 && !variables_changed ("SNAMES",
7756 &new->snames, new->snames_location, new->location,
7757 &old->snames, old->snames_location, old->location))
7764 msave_common_destroy (struct msave_common *common)
7768 msg_location_destroy (common->location);
7769 fh_unref (common->outfile);
7770 msg_location_destroy (common->outfile_location);
7771 string_array_destroy (&common->variables);
7772 msg_location_destroy (common->variables_location);
7773 string_array_destroy (&common->fnames);
7774 msg_location_destroy (common->fnames_location);
7775 string_array_destroy (&common->snames);
7776 msg_location_destroy (common->snames_location);
7778 for (size_t i = 0; i < common->n_factors; i++)
7779 matrix_expr_destroy (common->factors[i]);
7780 free (common->factors);
7782 for (size_t i = 0; i < common->n_splits; i++)
7783 matrix_expr_destroy (common->splits[i]);
7784 free (common->splits);
7786 dict_unref (common->dict);
7787 casewriter_destroy (common->writer);
7794 match_rowtype (struct lexer *lexer)
7796 static const char *rowtypes[] = {
7797 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7799 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7801 for (size_t i = 0; i < n_rowtypes; i++)
7802 if (lex_match_id (lexer, rowtypes[i]))
7805 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7810 parse_var_names (struct lexer *lexer, struct string_array *sa,
7811 struct msg_location **locationp)
7813 lex_match (lexer, T_EQUALS);
7815 string_array_clear (sa);
7816 msg_location_destroy (*locationp);
7819 struct dictionary *dict = dict_create (get_default_encoding ());
7822 int start_ofs = lex_ofs (lexer);
7823 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7824 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7825 int end_ofs = lex_ofs (lexer) - 1;
7830 for (size_t i = 0; i < n_names; i++)
7831 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7832 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7834 lex_ofs_error (lexer, start_ofs, end_ofs,
7835 _("Variable name %s is reserved."), names[i]);
7836 for (size_t j = 0; j < n_names; j++)
7842 sa->strings = names;
7843 sa->n = sa->allocated = n_names;
7844 *locationp = lex_ofs_location (lexer, start_ofs, end_ofs);
7849 static struct matrix_command *
7850 matrix_msave_parse (struct matrix_state *s)
7852 int start_ofs = lex_ofs (s->lexer);
7854 struct msave_common *common = xmalloc (sizeof *common);
7855 *common = (struct msave_common) { .outfile = NULL };
7857 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7858 *cmd = (struct matrix_command) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7860 struct matrix_expr *splits = NULL;
7861 struct matrix_expr *factors = NULL;
7863 struct matrix_msave *msave = &cmd->msave;
7864 msave->expr = matrix_parse_exp (s);
7868 while (lex_match (s->lexer, T_SLASH))
7870 if (lex_match_id (s->lexer, "TYPE"))
7872 lex_match (s->lexer, T_EQUALS);
7874 msave->rowtype = match_rowtype (s->lexer);
7875 if (!msave->rowtype)
7878 else if (lex_match_id (s->lexer, "OUTFILE"))
7880 lex_match (s->lexer, T_EQUALS);
7882 fh_unref (common->outfile);
7883 int start_ofs = lex_ofs (s->lexer);
7884 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7885 if (!common->outfile)
7887 msg_location_destroy (common->outfile_location);
7888 common->outfile_location = lex_ofs_location (s->lexer, start_ofs,
7889 lex_ofs (s->lexer) - 1);
7891 else if (lex_match_id (s->lexer, "VARIABLES"))
7893 if (!parse_var_names (s->lexer, &common->variables,
7894 &common->variables_location))
7897 else if (lex_match_id (s->lexer, "FNAMES"))
7899 if (!parse_var_names (s->lexer, &common->fnames,
7900 &common->fnames_location))
7903 else if (lex_match_id (s->lexer, "SNAMES"))
7905 if (!parse_var_names (s->lexer, &common->snames,
7906 &common->snames_location))
7909 else if (lex_match_id (s->lexer, "SPLIT"))
7911 lex_match (s->lexer, T_EQUALS);
7913 matrix_expr_destroy (splits);
7914 splits = matrix_parse_exp (s);
7918 else if (lex_match_id (s->lexer, "FACTOR"))
7920 lex_match (s->lexer, T_EQUALS);
7922 matrix_expr_destroy (factors);
7923 factors = matrix_parse_exp (s);
7929 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7930 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7934 if (!msave->rowtype)
7936 lex_sbc_missing (s->lexer, "TYPE");
7940 if (!s->msave_common)
7942 if (common->fnames.n && !factors)
7944 msg_at (SE, common->fnames_location, _("FNAMES requires FACTOR."));
7947 if (common->snames.n && !splits)
7949 msg_at (SE, common->snames_location, _("SNAMES requires SPLIT."));
7952 if (!common->outfile)
7954 lex_sbc_missing (s->lexer, "OUTFILE");
7957 common->location = lex_ofs_location (s->lexer, start_ofs,
7958 lex_ofs (s->lexer));
7959 msg_location_remove_columns (common->location);
7960 s->msave_common = common;
7964 if (msave_common_changed (s->msave_common, common))
7966 msave_common_destroy (common);
7968 msave->common = s->msave_common;
7970 struct msave_common *c = s->msave_common;
7973 if (c->n_factors >= c->allocated_factors)
7974 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7975 sizeof *c->factors);
7976 c->factors[c->n_factors++] = factors;
7978 if (c->n_factors > 0)
7979 msave->factors = c->factors[c->n_factors - 1];
7983 if (c->n_splits >= c->allocated_splits)
7984 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7986 c->splits[c->n_splits++] = splits;
7988 if (c->n_splits > 0)
7989 msave->splits = c->splits[c->n_splits - 1];
7994 matrix_expr_destroy (splits);
7995 matrix_expr_destroy (factors);
7996 msave_common_destroy (common);
7997 matrix_command_destroy (cmd);
8002 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
8004 gsl_matrix *m = matrix_expr_evaluate (e);
8010 msg_at (SE, matrix_expr_location (e),
8011 _("%s expression must evaluate to vector, "
8012 "not a %zu×%zu matrix."),
8013 name, m->size1, m->size2);
8014 gsl_matrix_free (m);
8018 return matrix_to_vector (m);
8022 msave_add_vars (struct dictionary *d, const struct string_array *vars)
8024 for (size_t i = 0; i < vars->n; i++)
8025 if (!dict_create_var (d, vars->strings[i], 0))
8026 return vars->strings[i];
8030 static struct dictionary *
8031 msave_create_dict (const struct msave_common *common)
8033 struct dictionary *dict = dict_create (get_default_encoding ());
8035 const char *dup_split = msave_add_vars (dict, &common->snames);
8038 /* Should not be possible because the parser ensures that the names are
8043 dict_create_var_assert (dict, "ROWTYPE_", 8);
8045 const char *dup_factor = msave_add_vars (dict, &common->fnames);
8048 msg_at (SE, common->fnames_location,
8049 _("Duplicate or invalid FACTOR variable name %s."),
8054 dict_create_var_assert (dict, "VARNAME_", 8);
8056 const char *dup_var = msave_add_vars (dict, &common->variables);
8059 msg_at (SE, common->variables_location,
8060 _("Duplicate or invalid variable name %s."),
8073 matrix_msave_execute (struct matrix_command *cmd)
8075 struct matrix_msave *msave = &cmd->msave;
8076 struct msave_common *common = msave->common;
8077 gsl_matrix *m = NULL;
8078 gsl_vector *factors = NULL;
8079 gsl_vector *splits = NULL;
8081 m = matrix_expr_evaluate (msave->expr);
8085 if (!common->variables.n)
8086 for (size_t i = 0; i < m->size2; i++)
8087 string_array_append_nocopy (&common->variables,
8088 xasprintf ("COL%zu", i + 1));
8089 else if (m->size2 != common->variables.n)
8091 msg_at (SE, matrix_expr_location (msave->expr),
8092 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
8093 m->size2, common->variables.n);
8099 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
8103 if (!common->fnames.n)
8104 for (size_t i = 0; i < factors->size; i++)
8105 string_array_append_nocopy (&common->fnames,
8106 xasprintf ("FAC%zu", i + 1));
8107 else if (factors->size != common->fnames.n)
8109 msg_at (SE, matrix_expr_location (msave->factors),
8110 _("There are %zu factor variables, "
8111 "but %zu factor values were supplied."),
8112 common->fnames.n, factors->size);
8119 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
8123 if (!common->snames.n)
8124 for (size_t i = 0; i < splits->size; i++)
8125 string_array_append_nocopy (&common->snames,
8126 xasprintf ("SPL%zu", i + 1));
8127 else if (splits->size != common->snames.n)
8129 msg_at (SE, matrix_expr_location (msave->splits),
8130 _("There are %zu split variables, "
8131 "but %zu split values were supplied."),
8132 common->snames.n, splits->size);
8137 if (!common->writer)
8139 struct dictionary *dict = msave_create_dict (common);
8143 common->writer = any_writer_open (common->outfile, dict);
8144 if (!common->writer)
8150 common->dict = dict;
8153 bool matrix = (!strcmp (msave->rowtype, "COV")
8154 || !strcmp (msave->rowtype, "CORR"));
8155 for (size_t y = 0; y < m->size1; y++)
8157 struct ccase *c = case_create (dict_get_proto (common->dict));
8160 /* Split variables */
8162 for (size_t i = 0; i < splits->size; i++)
8163 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
8166 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8167 msave->rowtype, ' ');
8171 for (size_t i = 0; i < factors->size; i++)
8172 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
8175 const char *varname_ = (matrix && y < common->variables.n
8176 ? common->variables.strings[y]
8178 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8181 /* Continuous variables. */
8182 for (size_t x = 0; x < m->size2; x++)
8183 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
8184 casewriter_write (common->writer, c);
8188 gsl_matrix_free (m);
8189 gsl_vector_free (factors);
8190 gsl_vector_free (splits);
8195 static struct matrix_command *
8196 matrix_mget_parse (struct matrix_state *s)
8198 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8199 *cmd = (struct matrix_command) {
8203 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
8207 struct matrix_mget *mget = &cmd->mget;
8209 lex_match (s->lexer, T_SLASH);
8210 while (lex_token (s->lexer) != T_ENDCMD)
8212 if (lex_match_id (s->lexer, "FILE"))
8214 lex_match (s->lexer, T_EQUALS);
8216 fh_unref (mget->file);
8217 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
8221 else if (lex_match_id (s->lexer, "ENCODING"))
8223 lex_match (s->lexer, T_EQUALS);
8224 if (!lex_force_string (s->lexer))
8227 free (mget->encoding);
8228 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
8232 else if (lex_match_id (s->lexer, "TYPE"))
8234 lex_match (s->lexer, T_EQUALS);
8235 while (lex_token (s->lexer) != T_SLASH
8236 && lex_token (s->lexer) != T_ENDCMD)
8238 const char *rowtype = match_rowtype (s->lexer);
8242 stringi_set_insert (&mget->rowtypes, rowtype);
8247 lex_error_expecting (s->lexer, "FILE", "TYPE");
8250 lex_match (s->lexer, T_SLASH);
8255 matrix_command_destroy (cmd);
8259 static const struct variable *
8260 get_a8_var (const struct msg_location *loc,
8261 const struct dictionary *d, const char *name)
8263 const struct variable *v = dict_lookup_var (d, name);
8266 msg_at (SE, loc, _("Matrix data file lacks %s variable."), name);
8269 if (var_get_width (v) != 8)
8271 msg_at (SE, loc, _("%s variable in matrix data file must be "
8272 "8-byte string, but it has width %d."),
8273 name, var_get_width (v));
8280 var_changed (const struct ccase *ca, const struct ccase *cb,
8281 const struct variable *var)
8284 ? !value_equal (case_data (ca, var), case_data (cb, var),
8285 var_get_width (var))
8290 vars_changed (const struct ccase *ca, const struct ccase *cb,
8291 const struct dictionary *d,
8292 size_t first_var, size_t n_vars)
8294 for (size_t i = 0; i < n_vars; i++)
8296 const struct variable *v = dict_get_var (d, first_var + i);
8297 if (var_changed (ca, cb, v))
8304 vars_all_missing (const struct ccase *c, const struct dictionary *d,
8305 size_t first_var, size_t n_vars)
8307 for (size_t i = 0; i < n_vars; i++)
8308 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
8314 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
8315 const struct dictionary *d,
8316 const struct variable *rowtype_var,
8317 const struct stringi_set *accepted_rowtypes,
8318 struct matrix_state *s,
8319 size_t ss, size_t sn, size_t si,
8320 size_t fs, size_t fn, size_t fi,
8321 size_t cs, size_t cn,
8322 struct pivot_table *pt,
8323 struct pivot_dimension *var_dimension)
8328 /* Is this a matrix for pooled data, either where there are no factor
8329 variables or the factor variables are missing? */
8330 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
8332 struct substring rowtype = case_ss (rows[0], rowtype_var);
8333 ss_rtrim (&rowtype, ss_cstr (" "));
8334 if (!stringi_set_is_empty (accepted_rowtypes)
8335 && !stringi_set_contains_len (accepted_rowtypes,
8336 rowtype.string, rowtype.length))
8339 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
8340 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
8341 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
8342 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
8343 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
8344 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
8348 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
8349 (int) rowtype.length, rowtype.string);
8353 struct string name = DS_EMPTY_INITIALIZER;
8354 ds_put_cstr (&name, prefix);
8356 ds_put_format (&name, "F%zu", fi);
8358 ds_put_format (&name, "S%zu", si);
8360 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
8362 mv = matrix_var_create (s, ds_ss (&name));
8365 msg (SW, _("Matrix data file contains variable with existing name %s."),
8367 goto exit_free_name;
8370 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
8371 size_t n_missing = 0;
8372 for (size_t y = 0; y < n_rows; y++)
8374 for (size_t x = 0; x < cn; x++)
8376 struct variable *var = dict_get_var (d, cs + x);
8377 double value = case_num (rows[y], var);
8378 if (var_is_num_missing (var, value))
8383 gsl_matrix_set (m, y, x, value);
8387 int var_index = pivot_category_create_leaf (
8388 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
8389 double values[] = { n_rows, cn };
8390 for (size_t j = 0; j < sn; j++)
8392 struct variable *var = dict_get_var (d, ss + j);
8393 const union value *value = case_data (rows[0], var);
8394 pivot_table_put2 (pt, j, var_index,
8395 pivot_value_new_var_value (var, value));
8397 for (size_t j = 0; j < fn; j++)
8399 struct variable *var = dict_get_var (d, fs + j);
8400 const union value sysmis = { .f = SYSMIS };
8401 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
8402 pivot_table_put2 (pt, j + sn, var_index,
8403 pivot_value_new_var_value (var, value));
8405 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
8406 pivot_table_put2 (pt, j + sn + fn, var_index,
8407 pivot_value_new_integer (values[j]));
8410 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
8411 "value, which was treated as zero.",
8412 "Matrix data file variable %s contains %zu missing "
8413 "values, which were treated as zero.", n_missing),
8414 ds_cstr (&name), n_missing);
8421 for (size_t y = 0; y < n_rows; y++)
8422 case_unref (rows[y]);
8426 matrix_mget_execute__ (struct matrix_command *cmd, struct casereader *r,
8427 const struct dictionary *d)
8429 struct matrix_mget *mget = &cmd->mget;
8430 const struct msg_location *loc = cmd->location;
8431 const struct variable *rowtype_ = get_a8_var (loc, d, "ROWTYPE_");
8432 const struct variable *varname_ = get_a8_var (loc, d, "VARNAME_");
8433 if (!rowtype_ || !varname_)
8436 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
8439 _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
8442 if (var_get_dict_index (varname_) + 1 >= dict_get_n_vars (d))
8444 msg_at (SE, loc, _("Matrix data file contains no continuous variables."));
8448 for (size_t i = 0; i < dict_get_n_vars (d); i++)
8450 const struct variable *v = dict_get_var (d, i);
8451 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
8454 _("Matrix data file contains unexpected string variable %s."),
8460 /* SPLIT variables. */
8462 size_t sn = var_get_dict_index (rowtype_);
8463 struct ccase *sc = NULL;
8466 /* FACTOR variables. */
8467 size_t fs = var_get_dict_index (rowtype_) + 1;
8468 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
8469 struct ccase *fc = NULL;
8472 /* Continuous variables. */
8473 size_t cs = var_get_dict_index (varname_) + 1;
8474 size_t cn = dict_get_n_vars (d) - cs;
8475 struct ccase *cc = NULL;
8478 struct pivot_table *pt = pivot_table_create (
8479 N_("Matrix Variables Created by MGET"));
8480 struct pivot_dimension *attr_dimension = pivot_dimension_create (
8481 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
8482 struct pivot_dimension *var_dimension = pivot_dimension_create (
8483 pt, PIVOT_AXIS_ROW, N_("Variable"));
8486 struct pivot_category *splits = pivot_category_create_group (
8487 attr_dimension->root, N_("Split Values"));
8488 for (size_t i = 0; i < sn; i++)
8489 pivot_category_create_leaf (splits, pivot_value_new_variable (
8490 dict_get_var (d, ss + i)));
8494 struct pivot_category *factors = pivot_category_create_group (
8495 attr_dimension->root, N_("Factors"));
8496 for (size_t i = 0; i < fn; i++)
8497 pivot_category_create_leaf (factors, pivot_value_new_variable (
8498 dict_get_var (d, fs + i)));
8500 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
8501 N_("Rows"), N_("Columns"));
8504 struct ccase **rows = NULL;
8505 size_t allocated_rows = 0;
8509 while ((c = casereader_read (r)) != NULL)
8511 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
8521 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
8522 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
8523 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
8526 if (change != NOTHING_CHANGED)
8528 matrix_mget_commit_var (rows, n_rows, d,
8529 rowtype_, &mget->rowtypes,
8540 if (n_rows >= allocated_rows)
8541 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
8544 if (change == SPLITS_CHANGED)
8550 /* Reset the factor number, if there are factors. */
8554 if (row_has_factors)
8560 else if (change == FACTORS_CHANGED)
8562 if (row_has_factors)
8568 matrix_mget_commit_var (rows, n_rows, d,
8569 rowtype_, &mget->rowtypes,
8581 if (var_dimension->n_leaves)
8582 pivot_table_submit (pt);
8584 pivot_table_unref (pt);
8588 matrix_mget_execute (struct matrix_command *cmd)
8590 struct matrix_mget *mget = &cmd->mget;
8591 struct casereader *r;
8592 struct dictionary *d;
8593 if (matrix_open_casereader (cmd, "MGET", mget->file, mget->encoding,
8594 mget->state->dataset, &r, &d))
8596 matrix_mget_execute__ (cmd, r, d);
8597 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
8604 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
8606 if (!lex_force_id (s->lexer))
8609 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
8611 *varp = matrix_var_create (s, lex_tokss (s->lexer));
8616 static struct matrix_command *
8617 matrix_eigen_parse (struct matrix_state *s)
8619 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8620 *cmd = (struct matrix_command) {
8622 .eigen = { .expr = NULL }
8625 struct matrix_eigen *eigen = &cmd->eigen;
8626 if (!lex_force_match (s->lexer, T_LPAREN))
8628 eigen->expr = matrix_expr_parse (s);
8630 || !lex_force_match (s->lexer, T_COMMA)
8631 || !matrix_parse_dst_var (s, &eigen->evec)
8632 || !lex_force_match (s->lexer, T_COMMA)
8633 || !matrix_parse_dst_var (s, &eigen->eval)
8634 || !lex_force_match (s->lexer, T_RPAREN))
8640 matrix_command_destroy (cmd);
8645 matrix_eigen_execute (struct matrix_command *cmd)
8647 struct matrix_eigen *eigen = &cmd->eigen;
8648 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
8651 if (!is_symmetric (A))
8653 msg_at (SE, cmd->location, _("Argument of EIGEN must be symmetric."));
8654 gsl_matrix_free (A);
8658 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
8659 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
8660 gsl_vector v_eval = to_vector (eval);
8661 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
8662 gsl_eigen_symmv (A, &v_eval, evec, w);
8663 gsl_eigen_symmv_free (w);
8665 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
8667 gsl_matrix_free (A);
8669 gsl_matrix_free (eigen->eval->value);
8670 eigen->eval->value = eval;
8672 gsl_matrix_free (eigen->evec->value);
8673 eigen->evec->value = evec;
8678 static struct matrix_command *
8679 matrix_setdiag_parse (struct matrix_state *s)
8681 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8682 *cmd = (struct matrix_command) {
8683 .type = MCMD_SETDIAG,
8684 .setdiag = { .dst = NULL }
8687 struct matrix_setdiag *setdiag = &cmd->setdiag;
8688 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
8691 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
8694 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
8699 if (!lex_force_match (s->lexer, T_COMMA))
8702 setdiag->expr = matrix_expr_parse (s);
8706 if (!lex_force_match (s->lexer, T_RPAREN))
8712 matrix_command_destroy (cmd);
8717 matrix_setdiag_execute (struct matrix_command *cmd)
8719 struct matrix_setdiag *setdiag = &cmd->setdiag;
8720 gsl_matrix *dst = setdiag->dst->value;
8723 msg_at (SE, cmd->location,
8724 _("SETDIAG destination matrix %s is uninitialized."),
8725 setdiag->dst->name);
8729 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
8733 size_t n = MIN (dst->size1, dst->size2);
8734 if (is_scalar (src))
8736 double d = to_scalar (src);
8737 for (size_t i = 0; i < n; i++)
8738 gsl_matrix_set (dst, i, i, d);
8740 else if (is_vector (src))
8742 gsl_vector v = to_vector (src);
8743 for (size_t i = 0; i < n && i < v.size; i++)
8744 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
8747 msg_at (SE, matrix_expr_location (setdiag->expr),
8748 _("SETDIAG argument 2 must be a scalar or a vector, "
8749 "not a %zu×%zu matrix."),
8750 src->size1, src->size2);
8751 gsl_matrix_free (src);
8756 static struct matrix_command *
8757 matrix_svd_parse (struct matrix_state *s)
8759 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8760 *cmd = (struct matrix_command) {
8762 .svd = { .expr = NULL }
8765 struct matrix_svd *svd = &cmd->svd;
8766 if (!lex_force_match (s->lexer, T_LPAREN))
8768 svd->expr = matrix_expr_parse (s);
8770 || !lex_force_match (s->lexer, T_COMMA)
8771 || !matrix_parse_dst_var (s, &svd->u)
8772 || !lex_force_match (s->lexer, T_COMMA)
8773 || !matrix_parse_dst_var (s, &svd->s)
8774 || !lex_force_match (s->lexer, T_COMMA)
8775 || !matrix_parse_dst_var (s, &svd->v)
8776 || !lex_force_match (s->lexer, T_RPAREN))
8782 matrix_command_destroy (cmd);
8787 matrix_svd_execute (struct matrix_svd *svd)
8789 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8793 if (m->size1 >= m->size2)
8796 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8797 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8798 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8799 gsl_vector *work = gsl_vector_alloc (A->size2);
8800 gsl_linalg_SV_decomp (A, V, &Sv, work);
8801 gsl_vector_free (work);
8803 matrix_var_set (svd->u, A);
8804 matrix_var_set (svd->s, S);
8805 matrix_var_set (svd->v, V);
8809 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8810 gsl_matrix_transpose_memcpy (At, m);
8811 gsl_matrix_free (m);
8813 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8814 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8815 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8816 gsl_vector *work = gsl_vector_alloc (At->size2);
8817 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8818 gsl_vector_free (work);
8820 matrix_var_set (svd->v, At);
8821 matrix_var_set (svd->s, St);
8822 matrix_var_set (svd->u, Vt);
8826 /* The main MATRIX command logic. */
8829 matrix_command_execute (struct matrix_command *cmd)
8834 matrix_compute_execute (cmd);
8838 matrix_print_execute (&cmd->print);
8842 return matrix_do_if_execute (&cmd->do_if);
8845 matrix_loop_execute (&cmd->loop);
8852 matrix_display_execute (&cmd->display);
8856 matrix_release_execute (&cmd->release);
8860 matrix_save_execute (cmd);
8864 matrix_read_execute (cmd);
8868 matrix_write_execute (&cmd->write);
8872 matrix_get_execute (cmd);
8876 matrix_msave_execute (cmd);
8880 matrix_mget_execute (cmd);
8884 matrix_eigen_execute (cmd);
8888 matrix_setdiag_execute (cmd);
8892 matrix_svd_execute (&cmd->svd);
8900 matrix_command_destroy (struct matrix_command *cmd)
8905 msg_location_destroy (cmd->location);
8910 matrix_lvalue_destroy (cmd->compute.lvalue);
8911 matrix_expr_destroy (cmd->compute.rvalue);
8915 matrix_expr_destroy (cmd->print.expression);
8916 free (cmd->print.title);
8917 print_labels_destroy (cmd->print.rlabels);
8918 print_labels_destroy (cmd->print.clabels);
8922 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8924 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8925 matrix_commands_uninit (&cmd->do_if.clauses[i].commands);
8927 free (cmd->do_if.clauses);
8931 matrix_expr_destroy (cmd->loop.start);
8932 matrix_expr_destroy (cmd->loop.end);
8933 matrix_expr_destroy (cmd->loop.increment);
8934 matrix_expr_destroy (cmd->loop.top_condition);
8935 matrix_expr_destroy (cmd->loop.bottom_condition);
8936 matrix_commands_uninit (&cmd->loop.commands);
8946 free (cmd->release.vars);
8950 matrix_expr_destroy (cmd->save.expression);
8954 matrix_lvalue_destroy (cmd->read.dst);
8955 matrix_expr_destroy (cmd->read.size);
8959 matrix_expr_destroy (cmd->write.expression);
8960 free (cmd->write.format);
8964 matrix_lvalue_destroy (cmd->get.dst);
8965 fh_unref (cmd->get.file);
8966 free (cmd->get.encoding);
8967 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8971 matrix_expr_destroy (cmd->msave.expr);
8975 fh_unref (cmd->mget.file);
8976 stringi_set_destroy (&cmd->mget.rowtypes);
8980 matrix_expr_destroy (cmd->eigen.expr);
8984 matrix_expr_destroy (cmd->setdiag.expr);
8988 matrix_expr_destroy (cmd->svd.expr);
8995 matrix_commands_parse (struct matrix_state *s, struct matrix_commands *c,
8996 const char *command_name,
8997 const char *stop1, const char *stop2)
8999 lex_end_of_command (s->lexer);
9000 lex_discard_rest_of_command (s->lexer);
9002 size_t allocated = 0;
9005 while (lex_token (s->lexer) == T_ENDCMD)
9008 if (lex_at_phrase (s->lexer, stop1)
9009 || (stop2 && lex_at_phrase (s->lexer, stop2)))
9012 if (lex_at_phrase (s->lexer, "END MATRIX"))
9014 lex_next_error (s->lexer, 0, 1,
9015 _("Premature END MATRIX within %s."), command_name);
9019 struct matrix_command *cmd = matrix_command_parse (s);
9023 if (c->n >= allocated)
9024 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
9025 c->commands[c->n++] = cmd;
9030 matrix_commands_uninit (struct matrix_commands *cmds)
9032 for (size_t i = 0; i < cmds->n; i++)
9033 matrix_command_destroy (cmds->commands[i]);
9034 free (cmds->commands);
9037 struct matrix_command_name
9040 struct matrix_command *(*parse) (struct matrix_state *);
9043 static const struct matrix_command_name *
9044 matrix_command_name_parse (struct lexer *lexer)
9046 static const struct matrix_command_name commands[] = {
9047 { "COMPUTE", matrix_compute_parse },
9048 { "CALL EIGEN", matrix_eigen_parse },
9049 { "CALL SETDIAG", matrix_setdiag_parse },
9050 { "CALL SVD", matrix_svd_parse },
9051 { "PRINT", matrix_print_parse },
9052 { "DO IF", matrix_do_if_parse },
9053 { "LOOP", matrix_loop_parse },
9054 { "BREAK", matrix_break_parse },
9055 { "READ", matrix_read_parse },
9056 { "WRITE", matrix_write_parse },
9057 { "GET", matrix_get_parse },
9058 { "SAVE", matrix_save_parse },
9059 { "MGET", matrix_mget_parse },
9060 { "MSAVE", matrix_msave_parse },
9061 { "DISPLAY", matrix_display_parse },
9062 { "RELEASE", matrix_release_parse },
9064 static size_t n = sizeof commands / sizeof *commands;
9066 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
9067 if (lex_match_phrase (lexer, c->name))
9072 static struct matrix_command *
9073 matrix_command_parse (struct matrix_state *s)
9075 int start_ofs = lex_ofs (s->lexer);
9076 size_t nesting_level = SIZE_MAX;
9078 struct matrix_command *c = NULL;
9079 const struct matrix_command_name *cmd = matrix_command_name_parse (s->lexer);
9081 lex_error (s->lexer, _("Unknown matrix command."));
9082 else if (!cmd->parse)
9083 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
9087 nesting_level = output_open_group (
9088 group_item_create_nocopy (utf8_to_title (cmd->name),
9089 utf8_to_title (cmd->name)));
9095 c->location = lex_ofs_location (s->lexer, start_ofs, lex_ofs (s->lexer));
9096 msg_location_remove_columns (c->location);
9097 lex_end_of_command (s->lexer);
9099 lex_discard_rest_of_command (s->lexer);
9100 if (nesting_level != SIZE_MAX)
9101 output_close_groups (nesting_level);
9107 cmd_matrix (struct lexer *lexer, struct dataset *ds)
9109 if (!lex_force_match (lexer, T_ENDCMD))
9112 struct matrix_state state = {
9114 .session = dataset_session (ds),
9116 .vars = HMAP_INITIALIZER (state.vars),
9121 while (lex_match (lexer, T_ENDCMD))
9123 if (lex_token (lexer) == T_STOP)
9125 msg (SE, _("Unexpected end of input expecting matrix command."));
9129 if (lex_match_phrase (lexer, "END MATRIX"))
9132 struct matrix_command *c = matrix_command_parse (&state);
9135 matrix_command_execute (c);
9136 matrix_command_destroy (c);
9140 struct matrix_var *var, *next;
9141 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
9144 gsl_matrix_free (var->value);
9145 hmap_delete (&state.vars, &var->hmap_node);
9148 hmap_destroy (&state.vars);
9149 msave_common_destroy (state.msave_common);
9150 fh_unref (state.prev_read_file);
9151 for (size_t i = 0; i < state.n_read_files; i++)
9152 read_file_destroy (state.read_files[i]);
9153 free (state.read_files);
9154 fh_unref (state.prev_write_file);
9155 for (size_t i = 0; i < state.n_write_files; i++)
9156 write_file_destroy (state.write_files[i]);
9157 free (state.write_files);
9158 fh_unref (state.prev_save_file);
9159 for (size_t i = 0; i < state.n_save_files; i++)
9160 save_file_destroy (state.save_files[i]);
9161 free (state.save_files);