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))
5393 char *error = fmt_check_output__ (cmd->print.format);
5396 lex_next_error (s->lexer, -1, -1, "%s", error);
5401 cmd->print.use_default_format = false;
5403 else if (lex_match_id (s->lexer, "TITLE"))
5405 lex_match (s->lexer, T_EQUALS);
5406 if (!lex_force_string (s->lexer))
5408 free (cmd->print.title);
5409 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
5412 else if (lex_match_id (s->lexer, "SPACE"))
5414 lex_match (s->lexer, T_EQUALS);
5415 if (lex_match_id (s->lexer, "NEWPAGE"))
5416 cmd->print.space = -1;
5417 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
5419 cmd->print.space = lex_integer (s->lexer);
5425 else if (lex_match_id (s->lexer, "RLABELS"))
5426 parse_literal_print_labels (s, &cmd->print.rlabels);
5427 else if (lex_match_id (s->lexer, "CLABELS"))
5428 parse_literal_print_labels (s, &cmd->print.clabels);
5429 else if (lex_match_id (s->lexer, "RNAMES"))
5431 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
5434 else if (lex_match_id (s->lexer, "CNAMES"))
5436 if (!parse_expr_print_labels (s, &cmd->print.clabels))
5441 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
5442 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
5450 matrix_command_destroy (cmd);
5455 matrix_is_integer (const gsl_matrix *m)
5457 for (size_t y = 0; y < m->size1; y++)
5458 for (size_t x = 0; x < m->size2; x++)
5460 double d = gsl_matrix_get (m, y, x);
5468 matrix_max_magnitude (const gsl_matrix *m)
5471 for (size_t y = 0; y < m->size1; y++)
5472 for (size_t x = 0; x < m->size2; x++)
5474 double d = fabs (gsl_matrix_get (m, y, x));
5482 format_fits (struct fmt_spec format, double x)
5484 char *s = data_out (&(union value) { .f = x }, NULL,
5485 format, settings_get_fmt_settings ());
5486 bool fits = *s != '*' && !strchr (s, 'E');
5491 static struct fmt_spec
5492 default_format (const gsl_matrix *m, int *log_scale)
5496 double max = matrix_max_magnitude (m);
5498 if (matrix_is_integer (m))
5499 for (int w = 1; w <= 12; w++)
5501 struct fmt_spec format = { .type = FMT_F, .w = w };
5502 if (format_fits (format, -max))
5506 if (max >= 1e9 || max <= 1e-4)
5509 snprintf (s, sizeof s, "%.1e", max);
5511 const char *e = strchr (s, 'e');
5513 *log_scale = atoi (e + 1);
5516 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5520 trimmed_string (double d)
5522 struct substring s = ss_buffer ((char *) &d, sizeof d);
5523 ss_rtrim (&s, ss_cstr (" "));
5524 return ss_xstrdup (s);
5527 static struct string_array *
5528 print_labels_get (const struct print_labels *labels, size_t n,
5529 const char *prefix, bool truncate)
5534 struct string_array *out = xzalloc (sizeof *out);
5535 if (labels->literals.n)
5536 string_array_clone (out, &labels->literals);
5537 else if (labels->expr)
5539 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5540 if (m && is_vector (m))
5542 gsl_vector v = to_vector (m);
5543 for (size_t i = 0; i < v.size; i++)
5544 string_array_append_nocopy (out, trimmed_string (
5545 gsl_vector_get (&v, i)));
5547 gsl_matrix_free (m);
5553 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5556 string_array_append (out, "");
5559 string_array_delete (out, out->n - 1);
5562 for (size_t i = 0; i < out->n; i++)
5564 char *s = out->strings[i];
5565 s[strnlen (s, 8)] = '\0';
5572 matrix_print_space (int space)
5575 output_item_submit (page_break_item_create ());
5576 for (int i = 0; i < space; i++)
5577 output_log ("%s", "");
5581 matrix_print_text (const struct matrix_print *print, const gsl_matrix *m,
5582 struct fmt_spec format, int log_scale)
5584 matrix_print_space (print->space);
5586 output_log ("%s", print->title);
5588 output_log (" 10 ** %d X", log_scale);
5590 struct string_array *clabels = print_labels_get (print->clabels,
5591 m->size2, "col", true);
5592 if (clabels && format.w < 8)
5594 struct string_array *rlabels = print_labels_get (print->rlabels,
5595 m->size1, "row", true);
5599 struct string line = DS_EMPTY_INITIALIZER;
5601 ds_put_byte_multiple (&line, ' ', 8);
5602 for (size_t x = 0; x < m->size2; x++)
5603 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5604 output_log_nocopy (ds_steal_cstr (&line));
5607 double scale = pow (10.0, log_scale);
5608 bool numeric = fmt_is_numeric (format.type);
5609 for (size_t y = 0; y < m->size1; y++)
5611 struct string line = DS_EMPTY_INITIALIZER;
5613 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5615 for (size_t x = 0; x < m->size2; x++)
5617 double f = gsl_matrix_get (m, y, x);
5619 ? data_out (&(union value) { .f = f / scale}, NULL,
5620 format, settings_get_fmt_settings ())
5621 : trimmed_string (f));
5622 ds_put_format (&line, " %s", s);
5625 output_log_nocopy (ds_steal_cstr (&line));
5628 string_array_destroy (rlabels);
5630 string_array_destroy (clabels);
5635 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5636 const struct print_labels *print_labels, size_t n,
5637 const char *name, const char *prefix)
5639 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5641 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5642 for (size_t i = 0; i < n; i++)
5643 pivot_category_create_leaf (
5645 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5646 : pivot_value_new_integer (i + 1)));
5648 d->hide_all_labels = true;
5649 string_array_destroy (labels);
5654 matrix_print_tables (const struct matrix_print *print, const gsl_matrix *m,
5655 struct fmt_spec format, int log_scale)
5657 struct pivot_table *table = pivot_table_create__ (
5658 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5661 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5663 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5664 N_("Columns"), "col");
5666 struct pivot_footnote *footnote = NULL;
5669 char *s = xasprintf ("× 10**%d", log_scale);
5670 footnote = pivot_table_create_footnote (
5671 table, pivot_value_new_user_text_nocopy (s));
5674 double scale = pow (10.0, log_scale);
5675 bool numeric = fmt_is_numeric (format.type);
5676 for (size_t y = 0; y < m->size1; y++)
5677 for (size_t x = 0; x < m->size2; x++)
5679 double f = gsl_matrix_get (m, y, x);
5680 struct pivot_value *v;
5683 v = pivot_value_new_number (f / scale);
5684 v->numeric.format = format;
5687 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5689 pivot_value_add_footnote (v, footnote);
5690 pivot_table_put2 (table, y, x, v);
5693 pivot_table_submit (table);
5697 matrix_print_execute (const struct matrix_print *print)
5699 if (print->expression)
5701 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5706 struct fmt_spec format = (print->use_default_format
5707 ? default_format (m, &log_scale)
5710 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5711 matrix_print_text (print, m, format, log_scale);
5713 matrix_print_tables (print, m, format, log_scale);
5715 gsl_matrix_free (m);
5719 matrix_print_space (print->space);
5721 output_log ("%s", print->title);
5728 matrix_do_if_clause_parse (struct matrix_state *s,
5729 struct matrix_do_if *ifc,
5730 bool parse_condition,
5731 size_t *allocated_clauses)
5733 if (ifc->n_clauses >= *allocated_clauses)
5734 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5735 sizeof *ifc->clauses);
5736 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5737 *c = (struct do_if_clause) { .condition = NULL };
5739 if (parse_condition)
5741 c->condition = matrix_expr_parse (s);
5746 return matrix_commands_parse (s, &c->commands, "DO IF", "ELSE", "END IF");
5749 static struct matrix_command *
5750 matrix_do_if_parse (struct matrix_state *s)
5752 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5753 *cmd = (struct matrix_command) {
5755 .do_if = { .n_clauses = 0 }
5758 size_t allocated_clauses = 0;
5761 if (!matrix_do_if_clause_parse (s, &cmd->do_if, true, &allocated_clauses))
5764 while (lex_match_phrase (s->lexer, "ELSE IF"));
5766 if (lex_match_id (s->lexer, "ELSE")
5767 && !matrix_do_if_clause_parse (s, &cmd->do_if, false, &allocated_clauses))
5770 if (!lex_match_phrase (s->lexer, "END IF"))
5775 matrix_command_destroy (cmd);
5780 matrix_do_if_execute (struct matrix_do_if *cmd)
5782 for (size_t i = 0; i < cmd->n_clauses; i++)
5784 struct do_if_clause *c = &cmd->clauses[i];
5788 if (!matrix_expr_evaluate_scalar (c->condition,
5789 i ? "ELSE IF" : "DO IF",
5794 for (size_t j = 0; j < c->commands.n; j++)
5795 if (!matrix_command_execute (c->commands.commands[j]))
5804 static struct matrix_command *
5805 matrix_loop_parse (struct matrix_state *s)
5807 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5808 *cmd = (struct matrix_command) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5810 struct matrix_loop *loop = &cmd->loop;
5811 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5813 struct substring name = lex_tokss (s->lexer);
5814 loop->var = matrix_var_lookup (s, name);
5816 loop->var = matrix_var_create (s, name);
5821 loop->start = matrix_expr_parse (s);
5822 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5825 loop->end = matrix_expr_parse (s);
5829 if (lex_match (s->lexer, T_BY))
5831 loop->increment = matrix_expr_parse (s);
5832 if (!loop->increment)
5837 if (lex_match_id (s->lexer, "IF"))
5839 loop->top_condition = matrix_expr_parse (s);
5840 if (!loop->top_condition)
5844 bool was_in_loop = s->in_loop;
5846 bool ok = matrix_commands_parse (s, &loop->commands, "LOOP",
5848 s->in_loop = was_in_loop;
5852 if (!lex_match_phrase (s->lexer, "END LOOP"))
5855 if (lex_match_id (s->lexer, "IF"))
5857 loop->bottom_condition = matrix_expr_parse (s);
5858 if (!loop->bottom_condition)
5865 matrix_command_destroy (cmd);
5870 matrix_loop_execute (struct matrix_loop *cmd)
5874 long int increment = 1;
5877 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5878 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5880 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5884 if (increment > 0 ? value > end
5885 : increment < 0 ? value < end
5890 int mxloops = settings_get_mxloops ();
5891 for (int i = 0; i < mxloops; i++)
5896 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5898 gsl_matrix_free (cmd->var->value);
5899 cmd->var->value = NULL;
5901 if (!cmd->var->value)
5902 cmd->var->value = gsl_matrix_alloc (1, 1);
5904 gsl_matrix_set (cmd->var->value, 0, 0, value);
5907 if (cmd->top_condition)
5910 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5915 for (size_t j = 0; j < cmd->commands.n; j++)
5916 if (!matrix_command_execute (cmd->commands.commands[j]))
5919 if (cmd->bottom_condition)
5922 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5923 "END LOOP IF", &d) || d > 0)
5930 if (increment > 0 ? value > end : value < end)
5938 static struct matrix_command *
5939 matrix_break_parse (struct matrix_state *s)
5943 lex_next_error (s->lexer, -1, -1, _("BREAK not inside LOOP."));
5947 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5948 *cmd = (struct matrix_command) { .type = MCMD_BREAK };
5954 static struct matrix_command *
5955 matrix_display_parse (struct matrix_state *s)
5957 while (lex_token (s->lexer) != T_ENDCMD)
5959 if (!lex_match_id (s->lexer, "DICTIONARY")
5960 && !lex_match_id (s->lexer, "STATUS"))
5962 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5967 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5968 *cmd = (struct matrix_command) { .type = MCMD_DISPLAY, .display = { s } };
5973 compare_matrix_var_pointers (const void *a_, const void *b_)
5975 const struct matrix_var *const *ap = a_;
5976 const struct matrix_var *const *bp = b_;
5977 const struct matrix_var *a = *ap;
5978 const struct matrix_var *b = *bp;
5979 return strcmp (a->name, b->name);
5983 matrix_display_execute (struct matrix_display *cmd)
5985 const struct matrix_state *s = cmd->state;
5987 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5988 struct pivot_dimension *attr_dimension
5989 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5990 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5991 N_("Rows"), N_("Columns"));
5992 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5994 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5997 const struct matrix_var *var;
5998 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
6000 vars[n_vars++] = var;
6001 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
6003 struct pivot_dimension *rows = pivot_dimension_create (
6004 table, PIVOT_AXIS_ROW, N_("Variable"));
6005 for (size_t i = 0; i < n_vars; i++)
6007 const struct matrix_var *var = vars[i];
6008 pivot_category_create_leaf (
6009 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
6011 size_t r = var->value->size1;
6012 size_t c = var->value->size2;
6013 double values[] = { r, c, r * c * sizeof (double) / 1024 };
6014 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
6015 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
6018 pivot_table_submit (table);
6023 static struct matrix_command *
6024 matrix_release_parse (struct matrix_state *s)
6026 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6027 *cmd = (struct matrix_command) { .type = MCMD_RELEASE };
6029 size_t allocated_vars = 0;
6030 while (lex_token (s->lexer) == T_ID)
6032 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
6035 if (cmd->release.n_vars >= allocated_vars)
6036 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
6037 sizeof *cmd->release.vars);
6038 cmd->release.vars[cmd->release.n_vars++] = var;
6041 lex_error (s->lexer, _("Syntax error expecting variable name."));
6044 if (!lex_match (s->lexer, T_COMMA))
6052 matrix_release_execute (struct matrix_release *cmd)
6054 for (size_t i = 0; i < cmd->n_vars; i++)
6056 struct matrix_var *var = cmd->vars[i];
6057 gsl_matrix_free (var->value);
6064 static struct save_file *
6065 save_file_create (struct matrix_state *s, struct file_handle *fh,
6066 struct string_array *variables,
6067 struct matrix_expr *names,
6068 struct stringi_set *strings)
6070 for (size_t i = 0; i < s->n_save_files; i++)
6072 struct save_file *sf = s->save_files[i];
6073 if (fh_equal (sf->file, fh))
6077 string_array_destroy (variables);
6078 matrix_expr_destroy (names);
6079 stringi_set_destroy (strings);
6085 struct save_file *sf = xmalloc (sizeof *sf);
6086 *sf = (struct save_file) {
6088 .dataset = s->dataset,
6089 .variables = *variables,
6091 .strings = STRINGI_SET_INITIALIZER (sf->strings),
6094 stringi_set_swap (&sf->strings, strings);
6095 stringi_set_destroy (strings);
6097 s->save_files = xrealloc (s->save_files,
6098 (s->n_save_files + 1) * sizeof *s->save_files);
6099 s->save_files[s->n_save_files++] = sf;
6103 static struct casewriter *
6104 save_file_open (struct save_file *sf, gsl_matrix *m,
6105 const struct msg_location *save_location)
6107 if (sf->writer || sf->error)
6111 size_t n_variables = caseproto_get_n_widths (
6112 casewriter_get_proto (sf->writer));
6113 if (m->size2 != n_variables)
6115 const char *file_name = (sf->file == fh_inline_file () ? "*"
6116 : fh_get_name (sf->file));
6117 msg_at (SE, save_location,
6118 _("Cannot save %zu×%zu matrix to %s because the "
6119 "first SAVE to %s in this matrix program wrote a "
6120 "%zu-column matrix."),
6121 m->size1, m->size2, file_name, file_name, n_variables);
6122 msg_at (SE, sf->location,
6123 _("This is the location of the first SAVE to %s."),
6132 struct dictionary *dict = dict_create (get_default_encoding ());
6134 /* Fill 'names' with user-specified names if there were any, either from
6135 sf->variables or sf->names. */
6136 struct string_array names = { .n = 0 };
6137 if (sf->variables.n)
6138 string_array_clone (&names, &sf->variables);
6141 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
6142 if (nm && is_vector (nm))
6144 gsl_vector nv = to_vector (nm);
6145 for (size_t i = 0; i < nv.size; i++)
6147 char *name = trimmed_string (gsl_vector_get (&nv, i));
6148 char *error = dict_id_is_valid__ (dict, name);
6150 string_array_append_nocopy (&names, name);
6153 msg_at (SE, save_location, "%s", error);
6159 gsl_matrix_free (nm);
6162 struct stringi_set strings;
6163 stringi_set_clone (&strings, &sf->strings);
6165 for (size_t i = 0; dict_get_n_vars (dict) < m->size2; i++)
6170 name = names.strings[i];
6173 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
6177 int width = stringi_set_delete (&strings, name) ? 8 : 0;
6178 struct variable *var = dict_create_var (dict, name, width);
6181 msg_at (ME, save_location,
6182 _("Duplicate variable name %s in SAVE statement."), name);
6187 if (!stringi_set_is_empty (&strings))
6189 size_t count = stringi_set_count (&strings);
6190 const char *example = stringi_set_node_get_string (
6191 stringi_set_first (&strings));
6194 msg_at (ME, save_location,
6195 _("The SAVE command STRINGS subcommand specifies an unknown "
6196 "variable %s."), example);
6198 msg_at (ME, save_location,
6199 ngettext ("The SAVE command STRINGS subcommand specifies %zu "
6200 "unknown variable: %s.",
6201 "The SAVE command STRINGS subcommand specifies %zu "
6202 "unknown variables, including %s.",
6207 stringi_set_destroy (&strings);
6208 string_array_destroy (&names);
6217 if (sf->file == fh_inline_file ())
6218 sf->writer = autopaging_writer_create (dict_get_proto (dict));
6220 sf->writer = any_writer_open (sf->file, dict);
6224 sf->location = msg_location_dup (save_location);
6235 save_file_destroy (struct save_file *sf)
6239 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
6241 dataset_set_dict (sf->dataset, sf->dict);
6242 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
6246 casewriter_destroy (sf->writer);
6247 dict_unref (sf->dict);
6249 fh_unref (sf->file);
6250 string_array_destroy (&sf->variables);
6251 matrix_expr_destroy (sf->names);
6252 stringi_set_destroy (&sf->strings);
6253 msg_location_destroy (sf->location);
6258 static struct matrix_command *
6259 matrix_save_parse (struct matrix_state *s)
6261 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6262 *cmd = (struct matrix_command) {
6264 .save = { .expression = NULL },
6267 struct file_handle *fh = NULL;
6268 struct matrix_save *save = &cmd->save;
6270 struct string_array variables = STRING_ARRAY_INITIALIZER;
6271 struct matrix_expr *names = NULL;
6272 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
6274 save->expression = matrix_parse_exp (s);
6275 if (!save->expression)
6278 int names_start = 0;
6280 while (lex_match (s->lexer, T_SLASH))
6282 if (lex_match_id (s->lexer, "OUTFILE"))
6284 lex_match (s->lexer, T_EQUALS);
6287 fh = (lex_match (s->lexer, T_ASTERISK)
6289 : fh_parse (s->lexer, FH_REF_FILE, s->session));
6293 else if (lex_match_id (s->lexer, "VARIABLES"))
6295 lex_match (s->lexer, T_EQUALS);
6299 struct dictionary *d = dict_create (get_default_encoding ());
6300 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
6301 PV_NO_SCRATCH | PV_NO_DUPLICATE);
6306 string_array_clear (&variables);
6307 variables = (struct string_array) {
6313 else if (lex_match_id (s->lexer, "NAMES"))
6315 lex_match (s->lexer, T_EQUALS);
6316 matrix_expr_destroy (names);
6317 names_start = lex_ofs (s->lexer);
6318 names = matrix_parse_exp (s);
6319 names_end = lex_ofs (s->lexer) - 1;
6323 else if (lex_match_id (s->lexer, "STRINGS"))
6325 lex_match (s->lexer, T_EQUALS);
6326 while (lex_token (s->lexer) == T_ID)
6328 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
6330 if (!lex_match (s->lexer, T_COMMA))
6336 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
6344 if (s->prev_save_file)
6345 fh = fh_ref (s->prev_save_file);
6348 lex_sbc_missing (s->lexer, "OUTFILE");
6352 fh_unref (s->prev_save_file);
6353 s->prev_save_file = fh_ref (fh);
6355 if (variables.n && names)
6357 lex_ofs_msg (s->lexer, SW, names_start, names_end,
6358 _("Ignoring NAMES because VARIABLES was also specified."));
6359 matrix_expr_destroy (names);
6363 save->sf = save_file_create (s, fh, &variables, names, &strings);
6367 string_array_destroy (&variables);
6368 matrix_expr_destroy (names);
6369 stringi_set_destroy (&strings);
6371 matrix_command_destroy (cmd);
6376 matrix_save_execute (const struct matrix_command *cmd)
6378 const struct matrix_save *save = &cmd->save;
6380 gsl_matrix *m = matrix_expr_evaluate (save->expression);
6384 struct casewriter *writer = save_file_open (save->sf, m, cmd->location);
6387 gsl_matrix_free (m);
6391 const struct caseproto *proto = casewriter_get_proto (writer);
6392 for (size_t y = 0; y < m->size1; y++)
6394 struct ccase *c = case_create (proto);
6395 for (size_t x = 0; x < m->size2; x++)
6397 double d = gsl_matrix_get (m, y, x);
6398 int width = caseproto_get_width (proto, x);
6399 union value *value = case_data_rw_idx (c, x);
6403 memcpy (value->s, &d, width);
6405 casewriter_write (writer, c);
6407 gsl_matrix_free (m);
6412 static struct read_file *
6413 read_file_create (struct matrix_state *s, struct file_handle *fh)
6415 for (size_t i = 0; i < s->n_read_files; i++)
6417 struct read_file *rf = s->read_files[i];
6425 struct read_file *rf = xmalloc (sizeof *rf);
6426 *rf = (struct read_file) { .file = fh };
6428 s->read_files = xrealloc (s->read_files,
6429 (s->n_read_files + 1) * sizeof *s->read_files);
6430 s->read_files[s->n_read_files++] = rf;
6434 static struct dfm_reader *
6435 read_file_open (struct read_file *rf)
6438 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
6443 read_file_destroy (struct read_file *rf)
6447 fh_unref (rf->file);
6448 dfm_close_reader (rf->reader);
6449 free (rf->encoding);
6454 static struct matrix_command *
6455 matrix_read_parse (struct matrix_state *s)
6457 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6458 *cmd = (struct matrix_command) {
6460 .read = { .format = FMT_F },
6463 struct file_handle *fh = NULL;
6464 char *encoding = NULL;
6465 struct matrix_read *read = &cmd->read;
6466 read->dst = matrix_lvalue_parse (s);
6472 int record_width_start = 0, record_width_end = 0;
6475 int repetitions = 0;
6476 int record_width = 0;
6477 bool seen_format = false;
6478 while (lex_match (s->lexer, T_SLASH))
6480 if (lex_match_id (s->lexer, "FILE"))
6482 lex_match (s->lexer, T_EQUALS);
6485 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6489 else if (lex_match_id (s->lexer, "ENCODING"))
6491 lex_match (s->lexer, T_EQUALS);
6492 if (!lex_force_string (s->lexer))
6496 encoding = ss_xstrdup (lex_tokss (s->lexer));
6500 else if (lex_match_id (s->lexer, "FIELD"))
6502 lex_match (s->lexer, T_EQUALS);
6504 record_width_start = lex_ofs (s->lexer);
6505 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6507 read->c1 = lex_integer (s->lexer);
6509 if (!lex_force_match (s->lexer, T_TO)
6510 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
6512 read->c2 = lex_integer (s->lexer) + 1;
6513 record_width_end = lex_ofs (s->lexer);
6516 record_width = read->c2 - read->c1;
6517 if (lex_match (s->lexer, T_BY))
6519 if (!lex_force_int_range (s->lexer, "BY", 1,
6520 read->c2 - read->c1))
6522 by = lex_integer (s->lexer);
6523 by_ofs = lex_ofs (s->lexer);
6524 int field_end = lex_ofs (s->lexer);
6527 if (record_width % by)
6530 s->lexer, record_width_start, field_end,
6531 _("Field width %d does not evenly divide record width %d."),
6533 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
6534 _("This syntax designates the record width."));
6535 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
6536 _("This syntax specifies the field width."));
6543 else if (lex_match_id (s->lexer, "SIZE"))
6545 lex_match (s->lexer, T_EQUALS);
6546 matrix_expr_destroy (read->size);
6547 read->size = matrix_parse_exp (s);
6551 else if (lex_match_id (s->lexer, "MODE"))
6553 lex_match (s->lexer, T_EQUALS);
6554 if (lex_match_id (s->lexer, "RECTANGULAR"))
6555 read->symmetric = false;
6556 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6557 read->symmetric = true;
6560 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6564 else if (lex_match_id (s->lexer, "REREAD"))
6565 read->reread = true;
6566 else if (lex_match_id (s->lexer, "FORMAT"))
6570 lex_sbc_only_once (s->lexer, "FORMAT");
6575 lex_match (s->lexer, T_EQUALS);
6577 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6580 format_ofs = lex_ofs (s->lexer);
6581 const char *p = lex_tokcstr (s->lexer);
6582 if (c_isdigit (p[0]))
6584 repetitions = atoi (p);
6585 p += strspn (p, "0123456789");
6586 if (!fmt_from_name (p, &read->format))
6588 lex_error (s->lexer, _("Unknown format %s."), p);
6593 else if (fmt_from_name (p, &read->format))
6597 struct fmt_spec format;
6598 if (!parse_format_specifier (s->lexer, &format))
6600 read->format = format.type;
6606 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6607 "REREAD", "FORMAT");
6614 lex_sbc_missing (s->lexer, "FIELD");
6618 if (!read->dst->n_indexes && !read->size)
6620 msg (SE, _("SIZE is required for reading data into a full matrix "
6621 "(as opposed to a submatrix)."));
6622 msg_at (SN, read->dst->var_location,
6623 _("This expression designates a full matrix."));
6629 if (s->prev_read_file)
6630 fh = fh_ref (s->prev_read_file);
6633 lex_sbc_missing (s->lexer, "FILE");
6637 fh_unref (s->prev_read_file);
6638 s->prev_read_file = fh_ref (fh);
6640 read->rf = read_file_create (s, fh);
6644 free (read->rf->encoding);
6645 read->rf->encoding = encoding;
6649 /* Field width may be specified in multiple ways:
6652 2. The format on FORMAT.
6653 3. The repetition factor on FORMAT.
6655 (2) and (3) are mutually exclusive.
6657 If more than one of these is present, they must agree. If none of them is
6658 present, then free-field format is used.
6660 if (repetitions > record_width)
6662 msg (SE, _("%d repetitions cannot fit in record width %d."),
6663 repetitions, record_width);
6664 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
6665 _("This syntax designates the number of repetitions."));
6666 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
6667 _("This syntax designates the record width."));
6670 int w = (repetitions ? record_width / repetitions
6675 msg (SE, _("This command specifies two different field widths."));
6678 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
6679 ngettext ("This syntax specifies %d repetition.",
6680 "This syntax specifies %d repetitions.",
6683 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
6684 _("This syntax designates record width %d, "
6685 "which divided by %d repetitions implies "
6687 record_width, repetitions, w);
6690 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
6691 _("This syntax specifies field width %d."), w);
6693 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
6694 _("This syntax specifies field width %d."), by);
6702 matrix_command_destroy (cmd);
6708 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6709 struct substring data, size_t y, size_t x,
6710 int first_column, int last_column, char *error)
6712 int line_number = dfm_get_line_number (reader);
6713 struct msg_location location = {
6714 .file_name = intern_new (dfm_get_file_name (reader)),
6715 .start = { .line = line_number, .column = first_column },
6716 .end = { .line = line_number, .column = last_column },
6718 msg_at (DW, &location, _("Error reading \"%.*s\" as format %s "
6719 "for matrix row %zu, column %zu: %s"),
6720 (int) data.length, data.string, fmt_name (format),
6721 y + 1, x + 1, error);
6722 msg_location_uninit (&location);
6727 matrix_read_set_field (struct matrix_read *read, struct dfm_reader *reader,
6728 gsl_matrix *m, struct substring p, size_t y, size_t x,
6729 const char *line_start)
6731 const char *input_encoding = dfm_reader_get_encoding (reader);
6734 if (fmt_is_numeric (read->format))
6737 error = data_in (p, input_encoding, read->format,
6738 settings_get_fmt_settings (), &v, 0, NULL);
6739 if (!error && v.f == SYSMIS)
6740 error = xstrdup (_("Matrix data may not contain missing value."));
6745 uint8_t s[sizeof (double)];
6746 union value v = { .s = s };
6747 error = data_in (p, input_encoding, read->format,
6748 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6749 memcpy (&f, s, sizeof f);
6754 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6755 int nc = ss_utf8_count_columns (p);
6756 int c2 = c1 + MAX (1, nc) - 1;
6757 parse_error (reader, read->format, p, y, x, c1, c2, error);
6761 gsl_matrix_set (m, y, x, f);
6762 if (read->symmetric && x != y)
6763 gsl_matrix_set (m, x, y, f);
6768 matrix_read_line (struct matrix_command *cmd, struct dfm_reader *reader,
6769 struct substring *line, const char **startp)
6771 struct matrix_read *read = &cmd->read;
6772 if (dfm_eof (reader))
6774 msg_at (SE, cmd->location,
6775 _("Unexpected end of file reading matrix data."));
6778 dfm_expand_tabs (reader);
6779 struct substring record = dfm_get_record (reader);
6780 /* XXX need to recode record into UTF-8 */
6781 *startp = record.string;
6782 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6787 matrix_read (struct matrix_command *cmd, struct dfm_reader *reader,
6790 struct matrix_read *read = &cmd->read;
6791 for (size_t y = 0; y < m->size1; y++)
6793 size_t nx = read->symmetric ? y + 1 : m->size2;
6795 struct substring line = ss_empty ();
6796 const char *line_start = line.string;
6797 for (size_t x = 0; x < nx; x++)
6804 ss_ltrim (&line, ss_cstr (" ,"));
6805 if (!ss_is_empty (line))
6807 if (!matrix_read_line (cmd, reader, &line, &line_start))
6809 dfm_forward_record (reader);
6812 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6816 if (!matrix_read_line (cmd, reader, &line, &line_start))
6818 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6819 int f = x % fields_per_line;
6820 if (f == fields_per_line - 1)
6821 dfm_forward_record (reader);
6823 p = ss_substr (line, read->w * f, read->w);
6826 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6830 dfm_forward_record (reader);
6833 ss_ltrim (&line, ss_cstr (" ,"));
6834 if (!ss_is_empty (line))
6836 int line_number = dfm_get_line_number (reader);
6837 int c1 = utf8_count_columns (line_start,
6838 line.string - line_start) + 1;
6839 int c2 = c1 + ss_utf8_count_columns (line) - 1;
6840 struct msg_location location = {
6841 .file_name = intern_new (dfm_get_file_name (reader)),
6842 .start = { .line = line_number, .column = c1 },
6843 .end = { .line = line_number, .column = c2 },
6845 msg_at (DW, &location,
6846 _("Trailing garbage following data for matrix row %zu."),
6848 msg_location_uninit (&location);
6855 matrix_read_execute (struct matrix_command *cmd)
6857 struct matrix_read *read = &cmd->read;
6858 struct index_vector iv0, iv1;
6859 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6862 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6865 gsl_matrix *m = matrix_expr_evaluate (read->size);
6871 msg_at (SE, matrix_expr_location (read->size),
6872 _("SIZE must evaluate to a scalar or a 2-element vector, "
6873 "not a %zu×%zu matrix."), m->size1, m->size2);
6874 gsl_matrix_free (m);
6875 index_vector_uninit (&iv0);
6876 index_vector_uninit (&iv1);
6880 gsl_vector v = to_vector (m);
6884 d[0] = gsl_vector_get (&v, 0);
6887 else if (v.size == 2)
6889 d[0] = gsl_vector_get (&v, 0);
6890 d[1] = gsl_vector_get (&v, 1);
6894 msg_at (SE, matrix_expr_location (read->size),
6895 _("SIZE must evaluate to a scalar or a 2-element vector, "
6896 "not a %zu×%zu matrix."),
6897 m->size1, m->size2),
6898 gsl_matrix_free (m);
6899 index_vector_uninit (&iv0);
6900 index_vector_uninit (&iv1);
6903 gsl_matrix_free (m);
6905 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6907 msg_at (SE, matrix_expr_location (read->size),
6908 _("Matrix dimensions %g×%g specified on SIZE "
6909 "are outside valid range."),
6911 index_vector_uninit (&iv0);
6912 index_vector_uninit (&iv1);
6920 if (read->dst->n_indexes)
6922 size_t submatrix_size[2];
6923 if (read->dst->n_indexes == 2)
6925 submatrix_size[0] = iv0.n;
6926 submatrix_size[1] = iv1.n;
6928 else if (read->dst->var->value->size1 == 1)
6930 submatrix_size[0] = 1;
6931 submatrix_size[1] = iv0.n;
6935 submatrix_size[0] = iv0.n;
6936 submatrix_size[1] = 1;
6941 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6943 msg_at (SE, cmd->location,
6944 _("Dimensions specified on SIZE differ from dimensions "
6945 "of destination submatrix."));
6946 msg_at (SN, matrix_expr_location (read->size),
6947 _("SIZE specifies dimensions %zu×%zu."),
6949 msg_at (SN, read->dst->full_location,
6950 _("Destination submatrix has dimensions %zu×%zu."),
6951 submatrix_size[0], submatrix_size[1]);
6952 index_vector_uninit (&iv0);
6953 index_vector_uninit (&iv1);
6959 size[0] = submatrix_size[0];
6960 size[1] = submatrix_size[1];
6964 struct dfm_reader *reader = read_file_open (read->rf);
6966 dfm_reread_record (reader, 1);
6968 if (read->symmetric && size[0] != size[1])
6970 msg_at (SE, cmd->location,
6971 _("Cannot read non-square %zu×%zu matrix "
6972 "using READ with MODE=SYMMETRIC."),
6974 index_vector_uninit (&iv0);
6975 index_vector_uninit (&iv1);
6978 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6979 matrix_read (cmd, reader, tmp);
6980 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp, cmd->location);
6985 static struct write_file *
6986 write_file_create (struct matrix_state *s, struct file_handle *fh)
6988 for (size_t i = 0; i < s->n_write_files; i++)
6990 struct write_file *wf = s->write_files[i];
6998 struct write_file *wf = xmalloc (sizeof *wf);
6999 *wf = (struct write_file) { .file = fh };
7001 s->write_files = xrealloc (s->write_files,
7002 (s->n_write_files + 1) * sizeof *s->write_files);
7003 s->write_files[s->n_write_files++] = wf;
7007 static struct dfm_writer *
7008 write_file_open (struct write_file *wf)
7011 wf->writer = dfm_open_writer (wf->file, wf->encoding);
7016 write_file_destroy (struct write_file *wf)
7022 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
7023 wf->held->s.ss.length);
7024 u8_line_destroy (wf->held);
7028 fh_unref (wf->file);
7029 dfm_close_writer (wf->writer);
7030 free (wf->encoding);
7035 static struct matrix_command *
7036 matrix_write_parse (struct matrix_state *s)
7038 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7039 *cmd = (struct matrix_command) {
7043 struct file_handle *fh = NULL;
7044 char *encoding = NULL;
7045 struct matrix_write *write = &cmd->write;
7046 write->expression = matrix_parse_exp (s);
7047 if (!write->expression)
7052 int record_width_start = 0, record_width_end = 0;
7055 int repetitions = 0;
7056 int record_width = 0;
7057 enum fmt_type format = FMT_F;
7058 bool has_format = false;
7059 while (lex_match (s->lexer, T_SLASH))
7061 if (lex_match_id (s->lexer, "OUTFILE"))
7063 lex_match (s->lexer, T_EQUALS);
7066 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
7070 else if (lex_match_id (s->lexer, "ENCODING"))
7072 lex_match (s->lexer, T_EQUALS);
7073 if (!lex_force_string (s->lexer))
7077 encoding = ss_xstrdup (lex_tokss (s->lexer));
7081 else if (lex_match_id (s->lexer, "FIELD"))
7083 lex_match (s->lexer, T_EQUALS);
7085 record_width_start = lex_ofs (s->lexer);
7087 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
7089 write->c1 = lex_integer (s->lexer);
7091 if (!lex_force_match (s->lexer, T_TO)
7092 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
7094 write->c2 = lex_integer (s->lexer) + 1;
7095 record_width_end = lex_ofs (s->lexer);
7098 record_width = write->c2 - write->c1;
7099 if (lex_match (s->lexer, T_BY))
7101 if (!lex_force_int_range (s->lexer, "BY", 1,
7102 write->c2 - write->c1))
7104 by_ofs = lex_ofs (s->lexer);
7105 int field_end = lex_ofs (s->lexer);
7106 by = lex_integer (s->lexer);
7109 if (record_width % by)
7112 s->lexer, record_width_start, field_end,
7113 _("Field width %d does not evenly divide record width %d."),
7115 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7116 _("This syntax designates the record width."));
7117 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
7118 _("This syntax specifies the field width."));
7125 else if (lex_match_id (s->lexer, "MODE"))
7127 lex_match (s->lexer, T_EQUALS);
7128 if (lex_match_id (s->lexer, "RECTANGULAR"))
7129 write->triangular = false;
7130 else if (lex_match_id (s->lexer, "TRIANGULAR"))
7131 write->triangular = true;
7134 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
7138 else if (lex_match_id (s->lexer, "HOLD"))
7140 else if (lex_match_id (s->lexer, "FORMAT"))
7142 if (has_format || write->format)
7144 lex_sbc_only_once (s->lexer, "FORMAT");
7148 lex_match (s->lexer, T_EQUALS);
7150 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
7153 format_ofs = lex_ofs (s->lexer);
7154 const char *p = lex_tokcstr (s->lexer);
7155 if (c_isdigit (p[0]))
7157 repetitions = atoi (p);
7158 p += strspn (p, "0123456789");
7159 if (!fmt_from_name (p, &format))
7161 lex_error (s->lexer, _("Unknown format %s."), p);
7167 else if (fmt_from_name (p, &format))
7174 struct fmt_spec spec;
7175 if (!parse_format_specifier (s->lexer, &spec))
7177 write->format = xmemdup (&spec, sizeof spec);
7182 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
7190 lex_sbc_missing (s->lexer, "FIELD");
7196 if (s->prev_write_file)
7197 fh = fh_ref (s->prev_write_file);
7200 lex_sbc_missing (s->lexer, "OUTFILE");
7204 fh_unref (s->prev_write_file);
7205 s->prev_write_file = fh_ref (fh);
7207 write->wf = write_file_create (s, fh);
7211 free (write->wf->encoding);
7212 write->wf->encoding = encoding;
7216 /* Field width may be specified in multiple ways:
7219 2. The format on FORMAT.
7220 3. The repetition factor on FORMAT.
7222 (2) and (3) are mutually exclusive.
7224 If more than one of these is present, they must agree. If none of them is
7225 present, then free-field format is used.
7227 if (repetitions > record_width)
7229 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7230 _("This syntax designates the number of repetitions."));
7231 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7232 _("This syntax designates the record width."));
7235 int w = (repetitions ? record_width / repetitions
7236 : write->format ? write->format->w
7240 msg (SE, _("This command specifies two different field widths."));
7243 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7244 ngettext ("This syntax specifies %d repetition.",
7245 "This syntax specifies %d repetitions.",
7248 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7249 _("This syntax designates record width %d, "
7250 "which divided by %d repetitions implies "
7252 record_width, repetitions, w);
7255 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7256 _("This syntax specifies field width %d."), w);
7258 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
7259 _("This syntax specifies field width %d."), by);
7262 if (w && !write->format)
7264 write->format = xmalloc (sizeof *write->format);
7265 *write->format = (struct fmt_spec) { .type = format, .w = w };
7267 char *error = fmt_check_output__ (*write->format);
7270 msg (SE, "%s", error);
7274 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7275 _("This syntax specifies format %s."),
7280 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7281 ngettext ("This syntax specifies %d repetition.",
7282 "This syntax specifies %d repetitions.",
7285 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7286 _("This syntax designates record width %d, "
7287 "which divided by %d repetitions implies "
7289 record_width, repetitions, w);
7293 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
7294 _("This syntax specifies field width %d."), by);
7300 if (write->format && fmt_var_width (*write->format) > sizeof (double))
7302 char fs[FMT_STRING_LEN_MAX + 1];
7303 fmt_to_string (*write->format, fs);
7304 lex_ofs_error (s->lexer, format_ofs, format_ofs,
7305 _("Format %s is too wide for %zu-byte matrix elements."),
7306 fs, sizeof (double));
7314 matrix_command_destroy (cmd);
7319 matrix_write_execute (struct matrix_write *write)
7321 gsl_matrix *m = matrix_expr_evaluate (write->expression);
7325 if (write->triangular && m->size1 != m->size2)
7327 msg_at (SE, matrix_expr_location (write->expression),
7328 _("WRITE with MODE=TRIANGULAR requires a square matrix but "
7329 "the matrix to be written has dimensions %zu×%zu."),
7330 m->size1, m->size2);
7331 gsl_matrix_free (m);
7335 struct dfm_writer *writer = write_file_open (write->wf);
7336 if (!writer || !m->size1)
7338 gsl_matrix_free (m);
7342 const struct fmt_settings *settings = settings_get_fmt_settings ();
7343 struct u8_line *line = write->wf->held;
7344 for (size_t y = 0; y < m->size1; y++)
7348 line = xmalloc (sizeof *line);
7349 u8_line_init (line);
7351 size_t nx = write->triangular ? y + 1 : m->size2;
7353 for (size_t x = 0; x < nx; x++)
7356 double f = gsl_matrix_get (m, y, x);
7360 if (fmt_is_numeric (write->format->type))
7363 v.s = (uint8_t *) &f;
7364 s = data_out (&v, NULL, *write->format, settings);
7368 s = xmalloc (DBL_BUFSIZE_BOUND);
7369 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
7370 >= DBL_BUFSIZE_BOUND)
7373 size_t len = strlen (s);
7374 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
7375 if (width + x0 > write->c2)
7377 dfm_put_record_utf8 (writer, line->s.ss.string,
7379 u8_line_clear (line);
7382 u8_line_put (line, x0, x0 + width, s, len);
7385 x0 += write->format ? write->format->w : width + 1;
7388 if (y + 1 >= m->size1 && write->hold)
7390 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
7391 u8_line_clear (line);
7395 u8_line_destroy (line);
7399 write->wf->held = line;
7401 gsl_matrix_free (m);
7406 static struct matrix_command *
7407 matrix_get_parse (struct matrix_state *s)
7409 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7410 *cmd = (struct matrix_command) {
7414 .dataset = s->dataset,
7415 .user = { .treatment = MGET_ERROR },
7416 .system = { .treatment = MGET_ERROR },
7420 struct matrix_get *get = &cmd->get;
7421 get->dst = matrix_lvalue_parse (s);
7425 while (lex_match (s->lexer, T_SLASH))
7427 if (lex_match_id (s->lexer, "FILE"))
7429 lex_match (s->lexer, T_EQUALS);
7431 fh_unref (get->file);
7432 if (lex_match (s->lexer, T_ASTERISK))
7436 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7441 else if (lex_match_id (s->lexer, "ENCODING"))
7443 lex_match (s->lexer, T_EQUALS);
7444 if (!lex_force_string (s->lexer))
7447 free (get->encoding);
7448 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
7452 else if (lex_match_id (s->lexer, "VARIABLES"))
7454 lex_match (s->lexer, T_EQUALS);
7458 lex_sbc_only_once (s->lexer, "VARIABLES");
7462 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
7465 else if (lex_match_id (s->lexer, "NAMES"))
7467 lex_match (s->lexer, T_EQUALS);
7468 if (!lex_force_id (s->lexer))
7471 struct substring name = lex_tokss (s->lexer);
7472 get->names = matrix_var_lookup (s, name);
7474 get->names = matrix_var_create (s, name);
7477 else if (lex_match_id (s->lexer, "MISSING"))
7479 lex_match (s->lexer, T_EQUALS);
7480 if (lex_match_id (s->lexer, "ACCEPT"))
7481 get->user.treatment = MGET_ACCEPT;
7482 else if (lex_match_id (s->lexer, "OMIT"))
7483 get->user.treatment = MGET_OMIT;
7484 else if (lex_is_number (s->lexer))
7486 get->user.treatment = MGET_RECODE;
7487 get->user.substitute = lex_number (s->lexer);
7492 lex_error (s->lexer, _("Syntax error expecting ACCEPT or OMIT or "
7493 "a number for MISSING."));
7497 else if (lex_match_id (s->lexer, "SYSMIS"))
7499 lex_match (s->lexer, T_EQUALS);
7500 if (lex_match_id (s->lexer, "OMIT"))
7501 get->system.treatment = MGET_OMIT;
7502 else if (lex_is_number (s->lexer))
7504 get->system.treatment = MGET_RECODE;
7505 get->system.substitute = lex_number (s->lexer);
7510 lex_error (s->lexer, _("Syntax error expecting OMIT or a number "
7517 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
7518 "MISSING", "SYSMIS");
7523 if (get->user.treatment != MGET_ACCEPT)
7524 get->system.treatment = MGET_ERROR;
7529 matrix_command_destroy (cmd);
7534 matrix_get_execute__ (struct matrix_command *cmd, struct casereader *reader,
7535 const struct dictionary *dict)
7537 struct matrix_get *get = &cmd->get;
7538 struct variable **vars;
7543 if (!var_syntax_evaluate (get->lexer, get->vars, get->n_vars, dict,
7544 &vars, &n_vars, PV_NUMERIC))
7549 n_vars = dict_get_n_vars (dict);
7550 vars = xnmalloc (n_vars, sizeof *vars);
7551 for (size_t i = 0; i < n_vars; i++)
7553 struct variable *var = dict_get_var (dict, i);
7554 if (!var_is_numeric (var))
7556 msg_at (SE, cmd->location, _("Variable %s is not numeric."),
7557 var_get_name (var));
7567 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
7568 for (size_t i = 0; i < n_vars; i++)
7570 char s[sizeof (double)];
7572 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
7573 memcpy (&f, s, sizeof f);
7574 gsl_matrix_set (names, i, 0, f);
7577 gsl_matrix_free (get->names->value);
7578 get->names->value = names;
7582 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
7583 long long int casenum = 1;
7585 for (struct ccase *c = casereader_read (reader); c;
7586 c = casereader_read (reader), casenum++)
7588 if (n_rows >= m->size1)
7590 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
7591 for (size_t y = 0; y < n_rows; y++)
7592 for (size_t x = 0; x < n_vars; x++)
7593 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
7594 gsl_matrix_free (m);
7599 for (size_t x = 0; x < n_vars; x++)
7601 const struct variable *var = vars[x];
7602 double d = case_num (c, var);
7605 if (get->system.treatment == MGET_RECODE)
7606 d = get->system.substitute;
7607 else if (get->system.treatment == MGET_OMIT)
7611 msg_at (SE, cmd->location, _("Variable %s in case %lld "
7612 "is system-missing."),
7613 var_get_name (var), casenum);
7617 else if (var_is_num_missing (var, d) == MV_USER)
7619 if (get->user.treatment == MGET_RECODE)
7620 d = get->user.substitute;
7621 else if (get->user.treatment == MGET_OMIT)
7623 else if (get->user.treatment != MGET_ACCEPT)
7625 msg_at (SE, cmd->location,
7626 _("Variable %s in case %lld has user-missing "
7628 var_get_name (var), casenum, d);
7632 gsl_matrix_set (m, n_rows, x, d);
7643 matrix_lvalue_evaluate_and_assign (get->dst, m, cmd->location);
7646 gsl_matrix_free (m);
7651 matrix_open_casereader (const struct matrix_command *cmd,
7652 const char *command_name, struct file_handle *file,
7653 const char *encoding, struct dataset *dataset,
7654 struct casereader **readerp, struct dictionary **dictp)
7658 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
7659 return *readerp != NULL;
7663 if (dict_get_n_vars (dataset_dict (dataset)) == 0)
7665 msg_at (SE, cmd->location,
7666 _("The %s command cannot read an empty active file."),
7670 *readerp = proc_open (dataset);
7671 *dictp = dict_ref (dataset_dict (dataset));
7677 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
7678 struct casereader *reader, struct dictionary *dict)
7681 casereader_destroy (reader);
7683 proc_commit (dataset);
7687 matrix_get_execute (struct matrix_command *cmd)
7689 struct matrix_get *get = &cmd->get;
7690 struct casereader *r;
7691 struct dictionary *d;
7692 if (matrix_open_casereader (cmd, "GET", get->file, get->encoding,
7693 get->dataset, &r, &d))
7695 matrix_get_execute__ (cmd, r, d);
7696 matrix_close_casereader (get->file, get->dataset, r, d);
7703 variables_changed (const char *keyword,
7704 const struct string_array *new_vars,
7705 const struct msg_location *new_vars_location,
7706 const struct msg_location *new_location,
7707 const struct string_array *old_vars,
7708 const struct msg_location *old_vars_location,
7709 const struct msg_location *old_location)
7715 msg_at (SE, new_location,
7716 _("%s may only be specified on MSAVE if it was specified "
7717 "on the first MSAVE within MATRIX."), keyword);
7718 msg_at (SN, old_location,
7719 _("The first MSAVE in MATRIX did not specify %s."),
7721 msg_at (SN, new_vars_location,
7722 _("This is the specification of %s on a later MSAVE."),
7726 if (!string_array_equal_case (old_vars, new_vars))
7728 msg_at (SE, new_location,
7729 _("%s must specify the same variables on each MSAVE "
7730 "within a given MATRIX."), keyword);
7731 msg_at (SE, old_vars_location,
7732 _("This is the specification of %s on the first MSAVE."),
7734 msg_at (SE, new_vars_location,
7735 _("This is a different specification of %s on a later MSAVE."),
7744 msave_common_changed (const struct msave_common *old,
7745 const struct msave_common *new)
7747 if (new->outfile && !fh_equal (old->outfile, new->outfile))
7749 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7750 "within a single MATRIX command."));
7751 msg_at (SN, old->outfile_location,
7752 _("This is the OUTFILE on the first MSAVE command."));
7753 msg_at (SN, new->outfile_location,
7754 _("This is the OUTFILE on a later MSAVE command."));
7758 if (!variables_changed ("VARIABLES",
7759 &new->variables, new->variables_location, new->location,
7760 &old->variables, old->variables_location, old->location)
7761 && !variables_changed ("FNAMES",
7762 &new->fnames, new->fnames_location, new->location,
7763 &old->fnames, old->fnames_location, old->location)
7764 && !variables_changed ("SNAMES",
7765 &new->snames, new->snames_location, new->location,
7766 &old->snames, old->snames_location, old->location))
7773 msave_common_destroy (struct msave_common *common)
7777 msg_location_destroy (common->location);
7778 fh_unref (common->outfile);
7779 msg_location_destroy (common->outfile_location);
7780 string_array_destroy (&common->variables);
7781 msg_location_destroy (common->variables_location);
7782 string_array_destroy (&common->fnames);
7783 msg_location_destroy (common->fnames_location);
7784 string_array_destroy (&common->snames);
7785 msg_location_destroy (common->snames_location);
7787 for (size_t i = 0; i < common->n_factors; i++)
7788 matrix_expr_destroy (common->factors[i]);
7789 free (common->factors);
7791 for (size_t i = 0; i < common->n_splits; i++)
7792 matrix_expr_destroy (common->splits[i]);
7793 free (common->splits);
7795 dict_unref (common->dict);
7796 casewriter_destroy (common->writer);
7803 match_rowtype (struct lexer *lexer)
7805 static const char *rowtypes[] = {
7806 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7808 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7810 for (size_t i = 0; i < n_rowtypes; i++)
7811 if (lex_match_id (lexer, rowtypes[i]))
7814 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7819 parse_var_names (struct lexer *lexer, struct string_array *sa,
7820 struct msg_location **locationp)
7822 lex_match (lexer, T_EQUALS);
7824 string_array_clear (sa);
7825 msg_location_destroy (*locationp);
7828 struct dictionary *dict = dict_create (get_default_encoding ());
7831 int start_ofs = lex_ofs (lexer);
7832 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7833 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7834 int end_ofs = lex_ofs (lexer) - 1;
7839 for (size_t i = 0; i < n_names; i++)
7840 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7841 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7843 lex_ofs_error (lexer, start_ofs, end_ofs,
7844 _("Variable name %s is reserved."), names[i]);
7845 for (size_t j = 0; j < n_names; j++)
7851 sa->strings = names;
7852 sa->n = sa->allocated = n_names;
7853 *locationp = lex_ofs_location (lexer, start_ofs, end_ofs);
7858 static struct matrix_command *
7859 matrix_msave_parse (struct matrix_state *s)
7861 int start_ofs = lex_ofs (s->lexer);
7863 struct msave_common *common = xmalloc (sizeof *common);
7864 *common = (struct msave_common) { .outfile = NULL };
7866 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7867 *cmd = (struct matrix_command) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7869 struct matrix_expr *splits = NULL;
7870 struct matrix_expr *factors = NULL;
7872 struct matrix_msave *msave = &cmd->msave;
7873 msave->expr = matrix_parse_exp (s);
7877 while (lex_match (s->lexer, T_SLASH))
7879 if (lex_match_id (s->lexer, "TYPE"))
7881 lex_match (s->lexer, T_EQUALS);
7883 msave->rowtype = match_rowtype (s->lexer);
7884 if (!msave->rowtype)
7887 else if (lex_match_id (s->lexer, "OUTFILE"))
7889 lex_match (s->lexer, T_EQUALS);
7891 fh_unref (common->outfile);
7892 int start_ofs = lex_ofs (s->lexer);
7893 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7894 if (!common->outfile)
7896 msg_location_destroy (common->outfile_location);
7897 common->outfile_location = lex_ofs_location (s->lexer, start_ofs,
7898 lex_ofs (s->lexer) - 1);
7900 else if (lex_match_id (s->lexer, "VARIABLES"))
7902 if (!parse_var_names (s->lexer, &common->variables,
7903 &common->variables_location))
7906 else if (lex_match_id (s->lexer, "FNAMES"))
7908 if (!parse_var_names (s->lexer, &common->fnames,
7909 &common->fnames_location))
7912 else if (lex_match_id (s->lexer, "SNAMES"))
7914 if (!parse_var_names (s->lexer, &common->snames,
7915 &common->snames_location))
7918 else if (lex_match_id (s->lexer, "SPLIT"))
7920 lex_match (s->lexer, T_EQUALS);
7922 matrix_expr_destroy (splits);
7923 splits = matrix_parse_exp (s);
7927 else if (lex_match_id (s->lexer, "FACTOR"))
7929 lex_match (s->lexer, T_EQUALS);
7931 matrix_expr_destroy (factors);
7932 factors = matrix_parse_exp (s);
7938 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7939 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7943 if (!msave->rowtype)
7945 lex_sbc_missing (s->lexer, "TYPE");
7949 if (!s->msave_common)
7951 if (common->fnames.n && !factors)
7953 msg_at (SE, common->fnames_location, _("FNAMES requires FACTOR."));
7956 if (common->snames.n && !splits)
7958 msg_at (SE, common->snames_location, _("SNAMES requires SPLIT."));
7961 if (!common->outfile)
7963 lex_sbc_missing (s->lexer, "OUTFILE");
7966 common->location = lex_ofs_location (s->lexer, start_ofs,
7967 lex_ofs (s->lexer));
7968 msg_location_remove_columns (common->location);
7969 s->msave_common = common;
7973 if (msave_common_changed (s->msave_common, common))
7975 msave_common_destroy (common);
7977 msave->common = s->msave_common;
7979 struct msave_common *c = s->msave_common;
7982 if (c->n_factors >= c->allocated_factors)
7983 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7984 sizeof *c->factors);
7985 c->factors[c->n_factors++] = factors;
7987 if (c->n_factors > 0)
7988 msave->factors = c->factors[c->n_factors - 1];
7992 if (c->n_splits >= c->allocated_splits)
7993 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7995 c->splits[c->n_splits++] = splits;
7997 if (c->n_splits > 0)
7998 msave->splits = c->splits[c->n_splits - 1];
8003 matrix_expr_destroy (splits);
8004 matrix_expr_destroy (factors);
8005 msave_common_destroy (common);
8006 matrix_command_destroy (cmd);
8011 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
8013 gsl_matrix *m = matrix_expr_evaluate (e);
8019 msg_at (SE, matrix_expr_location (e),
8020 _("%s expression must evaluate to vector, "
8021 "not a %zu×%zu matrix."),
8022 name, m->size1, m->size2);
8023 gsl_matrix_free (m);
8027 return matrix_to_vector (m);
8031 msave_add_vars (struct dictionary *d, const struct string_array *vars)
8033 for (size_t i = 0; i < vars->n; i++)
8034 if (!dict_create_var (d, vars->strings[i], 0))
8035 return vars->strings[i];
8039 static struct dictionary *
8040 msave_create_dict (const struct msave_common *common)
8042 struct dictionary *dict = dict_create (get_default_encoding ());
8044 const char *dup_split = msave_add_vars (dict, &common->snames);
8047 /* Should not be possible because the parser ensures that the names are
8052 dict_create_var_assert (dict, "ROWTYPE_", 8);
8054 const char *dup_factor = msave_add_vars (dict, &common->fnames);
8057 msg_at (SE, common->fnames_location,
8058 _("Duplicate or invalid FACTOR variable name %s."),
8063 dict_create_var_assert (dict, "VARNAME_", 8);
8065 const char *dup_var = msave_add_vars (dict, &common->variables);
8068 msg_at (SE, common->variables_location,
8069 _("Duplicate or invalid variable name %s."),
8082 matrix_msave_execute (struct matrix_command *cmd)
8084 struct matrix_msave *msave = &cmd->msave;
8085 struct msave_common *common = msave->common;
8086 gsl_matrix *m = NULL;
8087 gsl_vector *factors = NULL;
8088 gsl_vector *splits = NULL;
8090 m = matrix_expr_evaluate (msave->expr);
8094 if (!common->variables.n)
8095 for (size_t i = 0; i < m->size2; i++)
8096 string_array_append_nocopy (&common->variables,
8097 xasprintf ("COL%zu", i + 1));
8098 else if (m->size2 != common->variables.n)
8100 msg_at (SE, matrix_expr_location (msave->expr),
8101 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
8102 m->size2, common->variables.n);
8108 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
8112 if (!common->fnames.n)
8113 for (size_t i = 0; i < factors->size; i++)
8114 string_array_append_nocopy (&common->fnames,
8115 xasprintf ("FAC%zu", i + 1));
8116 else if (factors->size != common->fnames.n)
8118 msg_at (SE, matrix_expr_location (msave->factors),
8119 _("There are %zu factor variables, "
8120 "but %zu factor values were supplied."),
8121 common->fnames.n, factors->size);
8128 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
8132 if (!common->snames.n)
8133 for (size_t i = 0; i < splits->size; i++)
8134 string_array_append_nocopy (&common->snames,
8135 xasprintf ("SPL%zu", i + 1));
8136 else if (splits->size != common->snames.n)
8138 msg_at (SE, matrix_expr_location (msave->splits),
8139 _("There are %zu split variables, "
8140 "but %zu split values were supplied."),
8141 common->snames.n, splits->size);
8146 if (!common->writer)
8148 struct dictionary *dict = msave_create_dict (common);
8152 common->writer = any_writer_open (common->outfile, dict);
8153 if (!common->writer)
8159 common->dict = dict;
8162 bool matrix = (!strcmp (msave->rowtype, "COV")
8163 || !strcmp (msave->rowtype, "CORR"));
8164 for (size_t y = 0; y < m->size1; y++)
8166 struct ccase *c = case_create (dict_get_proto (common->dict));
8169 /* Split variables */
8171 for (size_t i = 0; i < splits->size; i++)
8172 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
8175 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8176 msave->rowtype, ' ');
8180 for (size_t i = 0; i < factors->size; i++)
8181 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
8184 const char *varname_ = (matrix && y < common->variables.n
8185 ? common->variables.strings[y]
8187 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8190 /* Continuous variables. */
8191 for (size_t x = 0; x < m->size2; x++)
8192 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
8193 casewriter_write (common->writer, c);
8197 gsl_matrix_free (m);
8198 gsl_vector_free (factors);
8199 gsl_vector_free (splits);
8204 static struct matrix_command *
8205 matrix_mget_parse (struct matrix_state *s)
8207 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8208 *cmd = (struct matrix_command) {
8212 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
8216 struct matrix_mget *mget = &cmd->mget;
8218 lex_match (s->lexer, T_SLASH);
8219 while (lex_token (s->lexer) != T_ENDCMD)
8221 if (lex_match_id (s->lexer, "FILE"))
8223 lex_match (s->lexer, T_EQUALS);
8225 fh_unref (mget->file);
8226 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
8230 else if (lex_match_id (s->lexer, "ENCODING"))
8232 lex_match (s->lexer, T_EQUALS);
8233 if (!lex_force_string (s->lexer))
8236 free (mget->encoding);
8237 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
8241 else if (lex_match_id (s->lexer, "TYPE"))
8243 lex_match (s->lexer, T_EQUALS);
8244 while (lex_token (s->lexer) != T_SLASH
8245 && lex_token (s->lexer) != T_ENDCMD)
8247 const char *rowtype = match_rowtype (s->lexer);
8251 stringi_set_insert (&mget->rowtypes, rowtype);
8256 lex_error_expecting (s->lexer, "FILE", "TYPE");
8259 lex_match (s->lexer, T_SLASH);
8264 matrix_command_destroy (cmd);
8268 static const struct variable *
8269 get_a8_var (const struct msg_location *loc,
8270 const struct dictionary *d, const char *name)
8272 const struct variable *v = dict_lookup_var (d, name);
8275 msg_at (SE, loc, _("Matrix data file lacks %s variable."), name);
8278 if (var_get_width (v) != 8)
8280 msg_at (SE, loc, _("%s variable in matrix data file must be "
8281 "8-byte string, but it has width %d."),
8282 name, var_get_width (v));
8289 var_changed (const struct ccase *ca, const struct ccase *cb,
8290 const struct variable *var)
8293 ? !value_equal (case_data (ca, var), case_data (cb, var),
8294 var_get_width (var))
8299 vars_changed (const struct ccase *ca, const struct ccase *cb,
8300 const struct dictionary *d,
8301 size_t first_var, size_t n_vars)
8303 for (size_t i = 0; i < n_vars; i++)
8305 const struct variable *v = dict_get_var (d, first_var + i);
8306 if (var_changed (ca, cb, v))
8313 vars_all_missing (const struct ccase *c, const struct dictionary *d,
8314 size_t first_var, size_t n_vars)
8316 for (size_t i = 0; i < n_vars; i++)
8317 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
8323 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
8324 const struct dictionary *d,
8325 const struct variable *rowtype_var,
8326 const struct stringi_set *accepted_rowtypes,
8327 struct matrix_state *s,
8328 size_t ss, size_t sn, size_t si,
8329 size_t fs, size_t fn, size_t fi,
8330 size_t cs, size_t cn,
8331 struct pivot_table *pt,
8332 struct pivot_dimension *var_dimension)
8337 /* Is this a matrix for pooled data, either where there are no factor
8338 variables or the factor variables are missing? */
8339 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
8341 struct substring rowtype = case_ss (rows[0], rowtype_var);
8342 ss_rtrim (&rowtype, ss_cstr (" "));
8343 if (!stringi_set_is_empty (accepted_rowtypes)
8344 && !stringi_set_contains_len (accepted_rowtypes,
8345 rowtype.string, rowtype.length))
8348 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
8349 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
8350 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
8351 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
8352 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
8353 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
8357 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
8358 (int) rowtype.length, rowtype.string);
8362 struct string name = DS_EMPTY_INITIALIZER;
8363 ds_put_cstr (&name, prefix);
8365 ds_put_format (&name, "F%zu", fi);
8367 ds_put_format (&name, "S%zu", si);
8369 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
8371 mv = matrix_var_create (s, ds_ss (&name));
8374 msg (SW, _("Matrix data file contains variable with existing name %s."),
8376 goto exit_free_name;
8379 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
8380 size_t n_missing = 0;
8381 for (size_t y = 0; y < n_rows; y++)
8383 for (size_t x = 0; x < cn; x++)
8385 struct variable *var = dict_get_var (d, cs + x);
8386 double value = case_num (rows[y], var);
8387 if (var_is_num_missing (var, value))
8392 gsl_matrix_set (m, y, x, value);
8396 int var_index = pivot_category_create_leaf (
8397 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
8398 double values[] = { n_rows, cn };
8399 for (size_t j = 0; j < sn; j++)
8401 struct variable *var = dict_get_var (d, ss + j);
8402 const union value *value = case_data (rows[0], var);
8403 pivot_table_put2 (pt, j, var_index,
8404 pivot_value_new_var_value (var, value));
8406 for (size_t j = 0; j < fn; j++)
8408 struct variable *var = dict_get_var (d, fs + j);
8409 const union value sysmis = { .f = SYSMIS };
8410 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
8411 pivot_table_put2 (pt, j + sn, var_index,
8412 pivot_value_new_var_value (var, value));
8414 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
8415 pivot_table_put2 (pt, j + sn + fn, var_index,
8416 pivot_value_new_integer (values[j]));
8419 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
8420 "value, which was treated as zero.",
8421 "Matrix data file variable %s contains %zu missing "
8422 "values, which were treated as zero.", n_missing),
8423 ds_cstr (&name), n_missing);
8430 for (size_t y = 0; y < n_rows; y++)
8431 case_unref (rows[y]);
8435 matrix_mget_execute__ (struct matrix_command *cmd, struct casereader *r,
8436 const struct dictionary *d)
8438 struct matrix_mget *mget = &cmd->mget;
8439 const struct msg_location *loc = cmd->location;
8440 const struct variable *rowtype_ = get_a8_var (loc, d, "ROWTYPE_");
8441 const struct variable *varname_ = get_a8_var (loc, d, "VARNAME_");
8442 if (!rowtype_ || !varname_)
8445 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
8448 _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
8451 if (var_get_dict_index (varname_) + 1 >= dict_get_n_vars (d))
8453 msg_at (SE, loc, _("Matrix data file contains no continuous variables."));
8457 for (size_t i = 0; i < dict_get_n_vars (d); i++)
8459 const struct variable *v = dict_get_var (d, i);
8460 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
8463 _("Matrix data file contains unexpected string variable %s."),
8469 /* SPLIT variables. */
8471 size_t sn = var_get_dict_index (rowtype_);
8472 struct ccase *sc = NULL;
8475 /* FACTOR variables. */
8476 size_t fs = var_get_dict_index (rowtype_) + 1;
8477 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
8478 struct ccase *fc = NULL;
8481 /* Continuous variables. */
8482 size_t cs = var_get_dict_index (varname_) + 1;
8483 size_t cn = dict_get_n_vars (d) - cs;
8484 struct ccase *cc = NULL;
8487 struct pivot_table *pt = pivot_table_create (
8488 N_("Matrix Variables Created by MGET"));
8489 struct pivot_dimension *attr_dimension = pivot_dimension_create (
8490 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
8491 struct pivot_dimension *var_dimension = pivot_dimension_create (
8492 pt, PIVOT_AXIS_ROW, N_("Variable"));
8495 struct pivot_category *splits = pivot_category_create_group (
8496 attr_dimension->root, N_("Split Values"));
8497 for (size_t i = 0; i < sn; i++)
8498 pivot_category_create_leaf (splits, pivot_value_new_variable (
8499 dict_get_var (d, ss + i)));
8503 struct pivot_category *factors = pivot_category_create_group (
8504 attr_dimension->root, N_("Factors"));
8505 for (size_t i = 0; i < fn; i++)
8506 pivot_category_create_leaf (factors, pivot_value_new_variable (
8507 dict_get_var (d, fs + i)));
8509 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
8510 N_("Rows"), N_("Columns"));
8513 struct ccase **rows = NULL;
8514 size_t allocated_rows = 0;
8518 while ((c = casereader_read (r)) != NULL)
8520 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
8530 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
8531 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
8532 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
8535 if (change != NOTHING_CHANGED)
8537 matrix_mget_commit_var (rows, n_rows, d,
8538 rowtype_, &mget->rowtypes,
8549 if (n_rows >= allocated_rows)
8550 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
8553 if (change == SPLITS_CHANGED)
8559 /* Reset the factor number, if there are factors. */
8563 if (row_has_factors)
8569 else if (change == FACTORS_CHANGED)
8571 if (row_has_factors)
8577 matrix_mget_commit_var (rows, n_rows, d,
8578 rowtype_, &mget->rowtypes,
8590 if (var_dimension->n_leaves)
8591 pivot_table_submit (pt);
8593 pivot_table_unref (pt);
8597 matrix_mget_execute (struct matrix_command *cmd)
8599 struct matrix_mget *mget = &cmd->mget;
8600 struct casereader *r;
8601 struct dictionary *d;
8602 if (matrix_open_casereader (cmd, "MGET", mget->file, mget->encoding,
8603 mget->state->dataset, &r, &d))
8605 matrix_mget_execute__ (cmd, r, d);
8606 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
8613 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
8615 if (!lex_force_id (s->lexer))
8618 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
8620 *varp = matrix_var_create (s, lex_tokss (s->lexer));
8625 static struct matrix_command *
8626 matrix_eigen_parse (struct matrix_state *s)
8628 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8629 *cmd = (struct matrix_command) {
8631 .eigen = { .expr = NULL }
8634 struct matrix_eigen *eigen = &cmd->eigen;
8635 if (!lex_force_match (s->lexer, T_LPAREN))
8637 eigen->expr = matrix_expr_parse (s);
8639 || !lex_force_match (s->lexer, T_COMMA)
8640 || !matrix_parse_dst_var (s, &eigen->evec)
8641 || !lex_force_match (s->lexer, T_COMMA)
8642 || !matrix_parse_dst_var (s, &eigen->eval)
8643 || !lex_force_match (s->lexer, T_RPAREN))
8649 matrix_command_destroy (cmd);
8654 matrix_eigen_execute (struct matrix_command *cmd)
8656 struct matrix_eigen *eigen = &cmd->eigen;
8657 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
8660 if (!is_symmetric (A))
8662 msg_at (SE, cmd->location, _("Argument of EIGEN must be symmetric."));
8663 gsl_matrix_free (A);
8667 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
8668 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
8669 gsl_vector v_eval = to_vector (eval);
8670 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
8671 gsl_eigen_symmv (A, &v_eval, evec, w);
8672 gsl_eigen_symmv_free (w);
8674 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
8676 gsl_matrix_free (A);
8678 gsl_matrix_free (eigen->eval->value);
8679 eigen->eval->value = eval;
8681 gsl_matrix_free (eigen->evec->value);
8682 eigen->evec->value = evec;
8687 static struct matrix_command *
8688 matrix_setdiag_parse (struct matrix_state *s)
8690 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8691 *cmd = (struct matrix_command) {
8692 .type = MCMD_SETDIAG,
8693 .setdiag = { .dst = NULL }
8696 struct matrix_setdiag *setdiag = &cmd->setdiag;
8697 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
8700 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
8703 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
8708 if (!lex_force_match (s->lexer, T_COMMA))
8711 setdiag->expr = matrix_expr_parse (s);
8715 if (!lex_force_match (s->lexer, T_RPAREN))
8721 matrix_command_destroy (cmd);
8726 matrix_setdiag_execute (struct matrix_command *cmd)
8728 struct matrix_setdiag *setdiag = &cmd->setdiag;
8729 gsl_matrix *dst = setdiag->dst->value;
8732 msg_at (SE, cmd->location,
8733 _("SETDIAG destination matrix %s is uninitialized."),
8734 setdiag->dst->name);
8738 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
8742 size_t n = MIN (dst->size1, dst->size2);
8743 if (is_scalar (src))
8745 double d = to_scalar (src);
8746 for (size_t i = 0; i < n; i++)
8747 gsl_matrix_set (dst, i, i, d);
8749 else if (is_vector (src))
8751 gsl_vector v = to_vector (src);
8752 for (size_t i = 0; i < n && i < v.size; i++)
8753 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
8756 msg_at (SE, matrix_expr_location (setdiag->expr),
8757 _("SETDIAG argument 2 must be a scalar or a vector, "
8758 "not a %zu×%zu matrix."),
8759 src->size1, src->size2);
8760 gsl_matrix_free (src);
8765 static struct matrix_command *
8766 matrix_svd_parse (struct matrix_state *s)
8768 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8769 *cmd = (struct matrix_command) {
8771 .svd = { .expr = NULL }
8774 struct matrix_svd *svd = &cmd->svd;
8775 if (!lex_force_match (s->lexer, T_LPAREN))
8777 svd->expr = matrix_expr_parse (s);
8779 || !lex_force_match (s->lexer, T_COMMA)
8780 || !matrix_parse_dst_var (s, &svd->u)
8781 || !lex_force_match (s->lexer, T_COMMA)
8782 || !matrix_parse_dst_var (s, &svd->s)
8783 || !lex_force_match (s->lexer, T_COMMA)
8784 || !matrix_parse_dst_var (s, &svd->v)
8785 || !lex_force_match (s->lexer, T_RPAREN))
8791 matrix_command_destroy (cmd);
8796 matrix_svd_execute (struct matrix_svd *svd)
8798 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8802 if (m->size1 >= m->size2)
8805 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8806 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8807 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8808 gsl_vector *work = gsl_vector_alloc (A->size2);
8809 gsl_linalg_SV_decomp (A, V, &Sv, work);
8810 gsl_vector_free (work);
8812 matrix_var_set (svd->u, A);
8813 matrix_var_set (svd->s, S);
8814 matrix_var_set (svd->v, V);
8818 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8819 gsl_matrix_transpose_memcpy (At, m);
8820 gsl_matrix_free (m);
8822 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8823 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8824 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8825 gsl_vector *work = gsl_vector_alloc (At->size2);
8826 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8827 gsl_vector_free (work);
8829 matrix_var_set (svd->v, At);
8830 matrix_var_set (svd->s, St);
8831 matrix_var_set (svd->u, Vt);
8835 /* The main MATRIX command logic. */
8838 matrix_command_execute (struct matrix_command *cmd)
8843 matrix_compute_execute (cmd);
8847 matrix_print_execute (&cmd->print);
8851 return matrix_do_if_execute (&cmd->do_if);
8854 matrix_loop_execute (&cmd->loop);
8861 matrix_display_execute (&cmd->display);
8865 matrix_release_execute (&cmd->release);
8869 matrix_save_execute (cmd);
8873 matrix_read_execute (cmd);
8877 matrix_write_execute (&cmd->write);
8881 matrix_get_execute (cmd);
8885 matrix_msave_execute (cmd);
8889 matrix_mget_execute (cmd);
8893 matrix_eigen_execute (cmd);
8897 matrix_setdiag_execute (cmd);
8901 matrix_svd_execute (&cmd->svd);
8909 matrix_command_destroy (struct matrix_command *cmd)
8914 msg_location_destroy (cmd->location);
8919 matrix_lvalue_destroy (cmd->compute.lvalue);
8920 matrix_expr_destroy (cmd->compute.rvalue);
8924 matrix_expr_destroy (cmd->print.expression);
8925 free (cmd->print.title);
8926 print_labels_destroy (cmd->print.rlabels);
8927 print_labels_destroy (cmd->print.clabels);
8931 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8933 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8934 matrix_commands_uninit (&cmd->do_if.clauses[i].commands);
8936 free (cmd->do_if.clauses);
8940 matrix_expr_destroy (cmd->loop.start);
8941 matrix_expr_destroy (cmd->loop.end);
8942 matrix_expr_destroy (cmd->loop.increment);
8943 matrix_expr_destroy (cmd->loop.top_condition);
8944 matrix_expr_destroy (cmd->loop.bottom_condition);
8945 matrix_commands_uninit (&cmd->loop.commands);
8955 free (cmd->release.vars);
8959 matrix_expr_destroy (cmd->save.expression);
8963 matrix_lvalue_destroy (cmd->read.dst);
8964 matrix_expr_destroy (cmd->read.size);
8968 matrix_expr_destroy (cmd->write.expression);
8969 free (cmd->write.format);
8973 matrix_lvalue_destroy (cmd->get.dst);
8974 fh_unref (cmd->get.file);
8975 free (cmd->get.encoding);
8976 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8980 matrix_expr_destroy (cmd->msave.expr);
8984 fh_unref (cmd->mget.file);
8985 stringi_set_destroy (&cmd->mget.rowtypes);
8989 matrix_expr_destroy (cmd->eigen.expr);
8993 matrix_expr_destroy (cmd->setdiag.expr);
8997 matrix_expr_destroy (cmd->svd.expr);
9004 matrix_commands_parse (struct matrix_state *s, struct matrix_commands *c,
9005 const char *command_name,
9006 const char *stop1, const char *stop2)
9008 lex_end_of_command (s->lexer);
9009 lex_discard_rest_of_command (s->lexer);
9011 size_t allocated = 0;
9014 while (lex_token (s->lexer) == T_ENDCMD)
9017 if (lex_at_phrase (s->lexer, stop1)
9018 || (stop2 && lex_at_phrase (s->lexer, stop2)))
9021 if (lex_at_phrase (s->lexer, "END MATRIX"))
9023 lex_next_error (s->lexer, 0, 1,
9024 _("Premature END MATRIX within %s."), command_name);
9028 struct matrix_command *cmd = matrix_command_parse (s);
9032 if (c->n >= allocated)
9033 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
9034 c->commands[c->n++] = cmd;
9039 matrix_commands_uninit (struct matrix_commands *cmds)
9041 for (size_t i = 0; i < cmds->n; i++)
9042 matrix_command_destroy (cmds->commands[i]);
9043 free (cmds->commands);
9046 struct matrix_command_name
9049 struct matrix_command *(*parse) (struct matrix_state *);
9052 static const struct matrix_command_name *
9053 matrix_command_name_parse (struct lexer *lexer)
9055 static const struct matrix_command_name commands[] = {
9056 { "COMPUTE", matrix_compute_parse },
9057 { "CALL EIGEN", matrix_eigen_parse },
9058 { "CALL SETDIAG", matrix_setdiag_parse },
9059 { "CALL SVD", matrix_svd_parse },
9060 { "PRINT", matrix_print_parse },
9061 { "DO IF", matrix_do_if_parse },
9062 { "LOOP", matrix_loop_parse },
9063 { "BREAK", matrix_break_parse },
9064 { "READ", matrix_read_parse },
9065 { "WRITE", matrix_write_parse },
9066 { "GET", matrix_get_parse },
9067 { "SAVE", matrix_save_parse },
9068 { "MGET", matrix_mget_parse },
9069 { "MSAVE", matrix_msave_parse },
9070 { "DISPLAY", matrix_display_parse },
9071 { "RELEASE", matrix_release_parse },
9073 static size_t n = sizeof commands / sizeof *commands;
9075 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
9076 if (lex_match_phrase (lexer, c->name))
9081 static struct matrix_command *
9082 matrix_command_parse (struct matrix_state *s)
9084 int start_ofs = lex_ofs (s->lexer);
9085 size_t nesting_level = SIZE_MAX;
9087 struct matrix_command *c = NULL;
9088 const struct matrix_command_name *cmd = matrix_command_name_parse (s->lexer);
9090 lex_error (s->lexer, _("Unknown matrix command."));
9091 else if (!cmd->parse)
9092 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
9096 nesting_level = output_open_group (
9097 group_item_create_nocopy (utf8_to_title (cmd->name),
9098 utf8_to_title (cmd->name)));
9104 c->location = lex_ofs_location (s->lexer, start_ofs, lex_ofs (s->lexer));
9105 msg_location_remove_columns (c->location);
9106 lex_end_of_command (s->lexer);
9108 lex_discard_rest_of_command (s->lexer);
9109 if (nesting_level != SIZE_MAX)
9110 output_close_groups (nesting_level);
9116 cmd_matrix (struct lexer *lexer, struct dataset *ds)
9118 if (!lex_force_match (lexer, T_ENDCMD))
9121 struct matrix_state state = {
9123 .session = dataset_session (ds),
9125 .vars = HMAP_INITIALIZER (state.vars),
9130 while (lex_match (lexer, T_ENDCMD))
9132 if (lex_token (lexer) == T_STOP)
9134 msg (SE, _("Unexpected end of input expecting matrix command."));
9138 if (lex_match_phrase (lexer, "END MATRIX"))
9141 struct matrix_command *c = matrix_command_parse (&state);
9144 matrix_command_execute (c);
9145 matrix_command_destroy (c);
9149 struct matrix_var *var, *next;
9150 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
9153 gsl_matrix_free (var->value);
9154 hmap_delete (&state.vars, &var->hmap_node);
9157 hmap_destroy (&state.vars);
9158 msave_common_destroy (state.msave_common);
9159 fh_unref (state.prev_read_file);
9160 for (size_t i = 0; i < state.n_read_files; i++)
9161 read_file_destroy (state.read_files[i]);
9162 free (state.read_files);
9163 fh_unref (state.prev_write_file);
9164 for (size_t i = 0; i < state.n_write_files; i++)
9165 write_file_destroy (state.write_files[i]);
9166 free (state.write_files);
9167 fh_unref (state.prev_save_file);
9168 for (size_t i = 0; i < state.n_save_files; i++)
9169 save_file_destroy (state.save_files[i]);
9170 free (state.save_files);