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 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3344 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3346 gsl_matrix *tmp = *a;
3352 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3355 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3356 swap_matrix (z, tmp);
3360 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3362 mul_matrix (x, *x, *x, tmp);
3366 matrix_expr_evaluate_exp_mat (const struct matrix_expr *e,
3367 gsl_matrix *x_, gsl_matrix *b)
3370 if (x->size1 != x->size2)
3372 msg_at (SE, matrix_expr_location (e->subs[0]),
3373 _("Matrix exponentation with ** requires a square matrix on "
3374 "the left-hand size, not one with dimensions %zu×%zu."),
3375 x->size1, x->size2);
3380 msg_at (SE, matrix_expr_location (e->subs[1]),
3381 _("Matrix exponentiation with ** requires a scalar on the "
3382 "right-hand side, not a matrix with dimensions %zu×%zu."),
3383 b->size1, b->size2);
3386 double bf = to_scalar (b);
3387 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3389 msg_at (SE, matrix_expr_location (e->subs[1]),
3390 _("Exponent %.1f in matrix exponentiation is non-integer "
3391 "or outside the valid range."), bf);
3396 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3398 gsl_matrix_set_identity (y);
3402 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3404 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3407 mul_matrix (&y, x, y, &t);
3408 square_matrix (&x, &t);
3411 square_matrix (&x, &t);
3413 mul_matrix (&y, x, y, &t);
3416 invert_matrix (y, x);
3417 swap_matrix (&x, &y);
3420 /* Garbage collection.
3422 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3423 a permutation of them. We are returning one of them; that one must not be
3424 destroyed. We must not destroy 'x_' because the caller owns it. */
3426 gsl_matrix_free (y_);
3428 gsl_matrix_free (t_);
3434 note_operand_size (const gsl_matrix *m, const struct matrix_expr *e)
3436 msg_at (SN, matrix_expr_location (e),
3437 _("This operand is a %zu×%zu matrix."), m->size1, m->size2);
3441 note_nonscalar (const gsl_matrix *m, const struct matrix_expr *e)
3444 note_operand_size (m, e);
3448 matrix_expr_evaluate_seq (const struct matrix_expr *e,
3449 gsl_matrix *start_, gsl_matrix *end_,
3452 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3454 msg_at (SE, matrix_expr_location (e),
3455 _("All operands of : operator must be scalars."));
3457 note_nonscalar (start_, e->subs[0]);
3458 note_nonscalar (end_, e->subs[1]);
3460 note_nonscalar (by_, e->subs[2]);
3464 long int start = to_scalar (start_);
3465 long int end = to_scalar (end_);
3466 long int by = by_ ? to_scalar (by_) : 1;
3470 msg_at (SE, matrix_expr_location (e->subs[2]),
3471 _("The increment operand to : must be nonzero."));
3475 long int n = (end >= start && by > 0 ? (end - start + by) / by
3476 : end <= start && by < 0 ? (start - end - by) / -by
3478 gsl_matrix *m = gsl_matrix_alloc (1, n);
3479 for (long int i = 0; i < n; i++)
3480 gsl_matrix_set (m, 0, i, start + i * by);
3485 matrix_expr_evaluate_not (gsl_matrix *a)
3487 MATRIX_FOR_ALL_ELEMENTS (d, y, x, a)
3493 matrix_expr_evaluate_paste_horz (const struct matrix_expr *e,
3494 gsl_matrix *a, gsl_matrix *b)
3496 if (a->size1 != b->size1)
3498 if (!a->size1 || !a->size2)
3500 else if (!b->size1 || !b->size2)
3503 msg_at (SE, matrix_expr_location (e),
3504 _("This expression tries to horizontally join matrices with "
3505 "differing numbers of rows."));
3506 note_operand_size (a, e->subs[0]);
3507 note_operand_size (b, e->subs[1]);
3511 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3512 for (size_t y = 0; y < a->size1; y++)
3514 for (size_t x = 0; x < a->size2; x++)
3515 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3516 for (size_t x = 0; x < b->size2; x++)
3517 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3523 matrix_expr_evaluate_paste_vert (const struct matrix_expr *e,
3524 gsl_matrix *a, gsl_matrix *b)
3526 if (a->size2 != b->size2)
3528 if (!a->size1 || !a->size2)
3530 else if (!b->size1 || !b->size2)
3533 msg_at (SE, matrix_expr_location (e),
3534 _("This expression tries to vertically join matrices with "
3535 "differing numbers of columns."));
3536 note_operand_size (a, e->subs[0]);
3537 note_operand_size (b, e->subs[1]);
3541 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3542 for (size_t x = 0; x < a->size2; x++)
3544 for (size_t y = 0; y < a->size1; y++)
3545 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3546 for (size_t y = 0; y < b->size1; y++)
3547 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3553 matrix_to_vector (gsl_matrix *m)
3556 gsl_vector v = to_vector (m);
3557 assert (v.block == m->block || !v.block);
3561 gsl_matrix_free (m);
3562 return xmemdup (&v, sizeof v);
3576 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3579 index_vector_uninit (struct index_vector *iv)
3586 matrix_normalize_index_vector (const gsl_matrix *m,
3587 const struct matrix_expr *me, size_t size,
3588 enum index_type index_type, size_t other_size,
3589 struct index_vector *iv)
3598 msg_at (SE, matrix_expr_location (me),
3599 _("Vector index must be scalar or vector, not a "
3601 m->size1, m->size2);
3605 msg_at (SE, matrix_expr_location (me),
3606 _("Matrix row index must be scalar or vector, not a "
3608 m->size1, m->size2);
3612 msg_at (SE, matrix_expr_location (me),
3613 _("Matrix column index must be scalar or vector, not a "
3615 m->size1, m->size2);
3621 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3622 *iv = (struct index_vector) {
3623 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3626 for (size_t i = 0; i < v.size; i++)
3628 double index = gsl_vector_get (&v, i);
3629 if (index < 1 || index >= size + 1)
3634 msg_at (SE, matrix_expr_location (me),
3635 _("Index %g is out of range for vector "
3636 "with %zu elements."), index, size);
3640 msg_at (SE, matrix_expr_location (me),
3641 _("%g is not a valid row index for "
3642 "a %zu×%zu matrix."),
3643 index, size, other_size);
3647 msg_at (SE, matrix_expr_location (me),
3648 _("%g is not a valid column index for "
3649 "a %zu×%zu matrix."),
3650 index, other_size, size);
3654 index_vector_uninit (iv);
3657 iv->indexes[i] = index - 1;
3663 *iv = (struct index_vector) {
3664 .indexes = xnmalloc (size, sizeof *iv->indexes),
3667 for (size_t i = 0; i < size; i++)
3674 matrix_expr_evaluate_vec_all (const struct matrix_expr *e,
3677 if (!is_vector (sm))
3679 msg_at (SE, matrix_expr_location (e->subs[0]),
3680 _("Vector index operator may not be applied to "
3681 "a %zu×%zu matrix."),
3682 sm->size1, sm->size2);
3690 matrix_expr_evaluate_vec_index (const struct matrix_expr *e,
3691 gsl_matrix *sm, gsl_matrix *im)
3693 if (!matrix_expr_evaluate_vec_all (e, sm))
3696 gsl_vector sv = to_vector (sm);
3697 struct index_vector iv;
3698 if (!matrix_normalize_index_vector (im, e->subs[1],
3699 sv.size, IV_VECTOR, 0, &iv))
3702 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3703 sm->size1 == 1 ? iv.n : 1);
3704 gsl_vector dv = to_vector (dm);
3705 for (size_t dx = 0; dx < iv.n; dx++)
3707 size_t sx = iv.indexes[dx];
3708 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3710 index_vector_uninit (&iv);
3716 matrix_expr_evaluate_mat_index (gsl_matrix *sm,
3717 gsl_matrix *im0, const struct matrix_expr *eim0,
3718 gsl_matrix *im1, const struct matrix_expr *eim1)
3720 struct index_vector iv0;
3721 if (!matrix_normalize_index_vector (im0, eim0, sm->size1,
3722 IV_ROW, sm->size2, &iv0))
3725 struct index_vector iv1;
3726 if (!matrix_normalize_index_vector (im1, eim1, sm->size2,
3727 IV_COLUMN, sm->size1, &iv1))
3729 index_vector_uninit (&iv0);
3733 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3734 for (size_t dy = 0; dy < iv0.n; dy++)
3736 size_t sy = iv0.indexes[dy];
3738 for (size_t dx = 0; dx < iv1.n; dx++)
3740 size_t sx = iv1.indexes[dx];
3741 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3744 index_vector_uninit (&iv0);
3745 index_vector_uninit (&iv1);
3749 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3750 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3751 const struct matrix_function_properties *, gsl_matrix *[], \
3752 const struct matrix_expr *, matrix_proto_##PROTO *);
3757 check_scalar_arg (const char *name, gsl_matrix *subs[],
3758 const struct matrix_expr *e, size_t index)
3760 if (!is_scalar (subs[index]))
3762 msg_at (SE, matrix_expr_location (e->subs[index]),
3763 _("Function %s argument %zu must be a scalar, "
3764 "not a %zu×%zu matrix."),
3765 name, index + 1, subs[index]->size1, subs[index]->size2);
3772 check_vector_arg (const char *name, gsl_matrix *subs[],
3773 const struct matrix_expr *e, size_t index)
3775 if (!is_vector (subs[index]))
3777 msg_at (SE, matrix_expr_location (e->subs[index]),
3778 _("Function %s argument %zu must be a vector, "
3779 "not a %zu×%zu matrix."),
3780 name, index + 1, subs[index]->size1, subs[index]->size2);
3787 to_scalar_args (const char *name, gsl_matrix *subs[],
3788 const struct matrix_expr *e, double d[])
3790 for (size_t i = 0; i < e->n_subs; i++)
3792 if (!check_scalar_arg (name, subs, e, i))
3794 d[i] = to_scalar (subs[i]);
3800 parse_constraint_value (const char **constraintsp)
3803 long retval = strtol (*constraintsp, &tail, 10);
3804 assert (tail > *constraintsp);
3805 *constraintsp = tail;
3809 enum matrix_argument_relop
3819 argument_inequality_error (
3820 const struct matrix_function_properties *props, const struct matrix_expr *e,
3821 size_t ai, gsl_matrix *a, size_t y, size_t x,
3822 size_t bi, double b,
3823 enum matrix_argument_relop relop)
3825 const struct msg_location *loc = matrix_expr_location (e);
3829 msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
3830 "than or equal to argument %zu."),
3831 ai + 1, props->name, bi + 1);
3835 msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
3836 "than argument %zu."),
3837 ai + 1, props->name, bi + 1);
3841 msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
3842 "or equal to argument %zu."),
3843 ai + 1, props->name, bi + 1);
3847 msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
3849 ai + 1, props->name, bi + 1);
3853 msg_at (ME, loc, _("Argument %zu to matrix function %s must not be equal "
3854 "to argument %zu."),
3855 ai + 1, props->name, bi + 1);
3859 const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
3861 msg_at (SN, a_loc, _("Argument %zu is %g."),
3862 ai + 1, gsl_matrix_get (a, y, x));
3864 msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
3865 y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
3867 msg_at (SN, matrix_expr_location (e->subs[bi]),
3868 _("Argument %zu is %g."), bi + 1, b);
3872 argument_value_error (
3873 const struct matrix_function_properties *props, const struct matrix_expr *e,
3874 size_t ai, gsl_matrix *a, size_t y, size_t x,
3876 enum matrix_argument_relop relop)
3878 const struct msg_location *loc = matrix_expr_location (e);
3882 msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
3883 "than or equal to %g."),
3884 ai + 1, props->name, b);
3888 msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
3890 ai + 1, props->name, b);
3894 msg_at (SE, loc, _("Argument %zu to matrix function %s must be less than "
3896 ai + 1, props->name, b);
3900 msg_at (SE, loc, _("Argument %zu to matrix function %s must be less than "
3902 ai + 1, props->name, b);
3906 msg_at (SE, loc, _("Argument %zu to matrix function %s must not be equal "
3908 ai + 1, props->name, b);
3912 const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
3915 if (relop != MRR_NE)
3916 msg_at (SN, a_loc, _("Argument %zu is %g."),
3917 ai + 1, gsl_matrix_get (a, y, x));
3920 msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
3921 y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
3925 matrix_argument_relop_is_satisfied (double a, double b,
3926 enum matrix_argument_relop relop)
3930 case MRR_GE: return a >= b;
3931 case MRR_GT: return a > b;
3932 case MRR_LE: return a <= b;
3933 case MRR_LT: return a < b;
3934 case MRR_NE: return a != b;
3940 static enum matrix_argument_relop
3941 matrix_argument_relop_flip (enum matrix_argument_relop relop)
3945 case MRR_GE: return MRR_LE;
3946 case MRR_GT: return MRR_LT;
3947 case MRR_LE: return MRR_GE;
3948 case MRR_LT: return MRR_GT;
3949 case MRR_NE: return MRR_NE;
3956 check_constraints (const struct matrix_function_properties *props,
3957 gsl_matrix *args[], const struct matrix_expr *e)
3959 size_t n_args = e->n_subs;
3960 const char *constraints = props->constraints;
3964 size_t arg_index = SIZE_MAX;
3965 while (*constraints)
3967 if (*constraints >= 'a' && *constraints <= 'd')
3969 arg_index = *constraints++ - 'a';
3970 assert (arg_index < n_args);
3972 else if (*constraints == '[' || *constraints == '(')
3974 assert (arg_index < n_args);
3975 bool open_lower = *constraints++ == '(';
3976 int minimum = parse_constraint_value (&constraints);
3977 assert (*constraints == ',');
3979 int maximum = parse_constraint_value (&constraints);
3980 assert (*constraints == ']' || *constraints == ')');
3981 bool open_upper = *constraints++ == ')';
3983 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3984 if ((open_lower ? *d <= minimum : *d < minimum)
3985 || (open_upper ? *d >= maximum : *d > maximum))
3987 if (!is_scalar (args[arg_index]))
3988 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
3989 _("Row %zu, column %zu of argument %zu to matrix "
3990 "function %s is %g, which is outside "
3991 "the valid range %c%d,%d%c."),
3992 y + 1, x + 1, arg_index + 1, props->name, *d,
3993 open_lower ? '(' : '[',
3995 open_upper ? ')' : ']');
3997 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
3998 _("Argument %zu to matrix function %s is %g, "
3999 "which is outside the valid range %c%d,%d%c."),
4000 arg_index + 1, props->name, *d,
4001 open_lower ? '(' : '[',
4003 open_upper ? ')' : ']');
4007 else if (*constraints == 'i')
4010 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
4011 if (*d != floor (*d))
4013 if (!is_scalar (args[arg_index]))
4014 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
4015 _("Argument %zu to matrix function %s, which must be "
4016 "integer, contains non-integer value %g in "
4017 "row %zu, column %zu."),
4018 arg_index + 1, props->name, *d, y + 1, x + 1);
4020 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
4021 _("Argument %zu to matrix function %s, which must be "
4022 "integer, has non-integer value %g."),
4023 arg_index + 1, props->name, *d);
4027 else if (*constraints == '>'
4028 || *constraints == '<'
4029 || *constraints == '!')
4031 enum matrix_argument_relop relop;
4032 switch (*constraints++)
4035 if (*constraints == '=')
4045 if (*constraints == '=')
4055 assert (*constraints == '=');
4064 if (*constraints >= 'a' && *constraints <= 'd')
4066 size_t a_index = arg_index;
4067 size_t b_index = *constraints - 'a';
4068 assert (a_index < n_args);
4069 assert (b_index < n_args);
4071 /* We only support one of the two arguments being non-scalar.
4072 It's easier to support only the first one being non-scalar, so
4073 flip things around if it's the other way. */
4074 if (!is_scalar (args[b_index]))
4076 assert (is_scalar (args[a_index]));
4077 size_t tmp_index = a_index;
4079 b_index = tmp_index;
4080 relop = matrix_argument_relop_flip (relop);
4083 double b = to_scalar (args[b_index]);
4084 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
4085 if (!matrix_argument_relop_is_satisfied (*a, b, relop))
4087 argument_inequality_error (
4089 a_index, args[a_index], y, x,
4097 int comparand = parse_constraint_value (&constraints);
4099 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
4100 if (!matrix_argument_relop_is_satisfied (*d, comparand, relop))
4102 argument_value_error (
4104 arg_index, args[arg_index], y, x,
4113 assert (*constraints == ' ');
4115 arg_index = SIZE_MAX;
4122 matrix_expr_evaluate_d_none (const struct matrix_function_properties *props,
4123 gsl_matrix *subs[], const struct matrix_expr *e,
4124 matrix_proto_d_none *f)
4126 assert (e->n_subs == 0);
4128 if (!check_constraints (props, subs, e))
4131 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4132 gsl_matrix_set (m, 0, 0, f ());
4137 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
4138 gsl_matrix *subs[], const struct matrix_expr *e,
4139 matrix_proto_d_d *f)
4141 assert (e->n_subs == 1);
4144 if (!to_scalar_args (props->name, subs, e, &d)
4145 || !check_constraints (props, subs, e))
4148 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4149 gsl_matrix_set (m, 0, 0, f (d));
4154 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
4155 gsl_matrix *subs[], const struct matrix_expr *e,
4156 matrix_proto_d_dd *f)
4158 assert (e->n_subs == 2);
4161 if (!to_scalar_args (props->name, subs, e, d)
4162 && !check_constraints (props, subs, e))
4165 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4166 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
4171 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
4172 gsl_matrix *subs[], const struct matrix_expr *e,
4173 matrix_proto_d_ddd *f)
4175 assert (e->n_subs == 3);
4178 if (!to_scalar_args (props->name, subs, e, d)
4179 || !check_constraints (props, subs, e))
4182 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4183 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
4188 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
4189 gsl_matrix *subs[], const struct matrix_expr *e,
4190 matrix_proto_m_d *f)
4192 assert (e->n_subs == 1);
4195 return (to_scalar_args (props->name, subs, e, &d)
4196 && check_constraints (props, subs, e)
4202 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
4203 gsl_matrix *subs[], const struct matrix_expr *e,
4204 matrix_proto_m_ddd *f)
4206 assert (e->n_subs == 3);
4209 return (to_scalar_args (props->name, subs, e, d)
4210 && check_constraints (props, subs, e)
4211 ? f(d[0], d[1], d[2])
4216 matrix_expr_evaluate_m_ddn (const struct matrix_function_properties *props,
4217 gsl_matrix *subs[], const struct matrix_expr *e,
4218 matrix_proto_m_ddn *f)
4220 assert (e->n_subs == 2);
4223 return (to_scalar_args (props->name, subs, e, d)
4224 && check_constraints (props, subs, e)
4230 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props,
4231 gsl_matrix *subs[], const struct matrix_expr *e,
4232 matrix_proto_m_m *f)
4234 assert (e->n_subs == 1);
4235 return check_constraints (props, subs, e) ? f (subs[0]) : NULL;
4239 matrix_expr_evaluate_m_mn (const struct matrix_function_properties *props,
4240 gsl_matrix *subs[], const struct matrix_expr *e,
4241 matrix_proto_m_mn *f)
4243 assert (e->n_subs == 1);
4244 return check_constraints (props, subs, e) ? f (subs[0], e) : NULL;
4248 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
4249 gsl_matrix *subs[], const struct matrix_expr *e,
4250 matrix_proto_m_e *f)
4252 assert (e->n_subs == 1);
4254 if (!check_constraints (props, subs, e))
4257 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4263 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props,
4264 gsl_matrix *subs[], const struct matrix_expr *e,
4265 matrix_proto_m_md *f)
4267 assert (e->n_subs == 2);
4268 return (check_scalar_arg (props->name, subs, e, 1)
4269 && check_constraints (props, subs, e)
4270 ? f (subs[0], to_scalar (subs[1]))
4275 matrix_expr_evaluate_m_mdn (const struct matrix_function_properties *props,
4276 gsl_matrix *subs[], const struct matrix_expr *e,
4277 matrix_proto_m_mdn *f)
4279 assert (e->n_subs == 2);
4280 return (check_scalar_arg (props->name, subs, e, 1)
4281 && check_constraints (props, subs, e)
4282 ? f (subs[0], to_scalar (subs[1]), e)
4287 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
4288 gsl_matrix *subs[], const struct matrix_expr *e,
4289 matrix_proto_m_ed *f)
4291 assert (e->n_subs == 2);
4292 if (!check_scalar_arg (props->name, subs, e, 1)
4293 || !check_constraints (props, subs, e))
4296 double b = to_scalar (subs[1]);
4297 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4303 matrix_expr_evaluate_m_mddn (const struct matrix_function_properties *props,
4304 gsl_matrix *subs[], const struct matrix_expr *e,
4305 matrix_proto_m_mddn *f)
4307 assert (e->n_subs == 3);
4308 if (!check_scalar_arg (props->name, subs, e, 1)
4309 || !check_scalar_arg (props->name, subs, e, 2)
4310 || !check_constraints (props, subs, e))
4312 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]), e);
4316 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
4317 gsl_matrix *subs[], const struct matrix_expr *e,
4318 matrix_proto_m_edd *f)
4320 assert (e->n_subs == 3);
4321 if (!check_scalar_arg (props->name, subs, e, 1)
4322 || !check_scalar_arg (props->name, subs, e, 2)
4323 || !check_constraints (props, subs, e))
4326 double b = to_scalar (subs[1]);
4327 double c = to_scalar (subs[2]);
4328 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4334 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4335 gsl_matrix *subs[], const struct matrix_expr *e,
4336 matrix_proto_m_eddd *f)
4338 assert (e->n_subs == 4);
4339 for (size_t i = 1; i < 4; i++)
4340 if (!check_scalar_arg (props->name, subs, e, i))
4343 if (!check_constraints (props, subs, e))
4346 double b = to_scalar (subs[1]);
4347 double c = to_scalar (subs[2]);
4348 double d = to_scalar (subs[3]);
4349 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4350 *a = f (*a, b, c, d);
4355 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4356 gsl_matrix *subs[], const struct matrix_expr *e,
4357 matrix_proto_m_eed *f)
4359 assert (e->n_subs == 3);
4360 if (!check_scalar_arg (props->name, subs, e, 2))
4363 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4364 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4366 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
4367 loc->end = e->subs[1]->location->end;
4370 _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4371 "%zu×%zu, but %s requires these arguments either to have "
4372 "the same dimensions or for one of them to be a scalar."),
4374 subs[0]->size1, subs[0]->size2,
4375 subs[1]->size1, subs[1]->size2,
4378 msg_location_destroy (loc);
4382 if (!check_constraints (props, subs, e))
4385 double c = to_scalar (subs[2]);
4387 if (is_scalar (subs[0]))
4389 double a = to_scalar (subs[0]);
4390 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4396 double b = to_scalar (subs[1]);
4397 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4404 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props,
4405 gsl_matrix *subs[], const struct matrix_expr *e,
4406 matrix_proto_m_mm *f)
4408 assert (e->n_subs == 2);
4409 return check_constraints (props, subs, e) ? f (subs[0], subs[1]) : NULL;
4413 matrix_expr_evaluate_m_mmn (const struct matrix_function_properties *props,
4414 gsl_matrix *subs[], const struct matrix_expr *e,
4415 matrix_proto_m_mmn *f)
4417 assert (e->n_subs == 2);
4418 return check_constraints (props, subs, e) ? f (subs[0], subs[1], e) : NULL;
4422 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4423 gsl_matrix *subs[], const struct matrix_expr *e,
4424 matrix_proto_m_v *f)
4426 assert (e->n_subs == 1);
4427 if (!check_vector_arg (props->name, subs, e, 0)
4428 || !check_constraints (props, subs, e))
4430 gsl_vector v = to_vector (subs[0]);
4435 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props,
4436 gsl_matrix *subs[], const struct matrix_expr *e,
4437 matrix_proto_d_m *f)
4439 assert (e->n_subs == 1);
4441 if (!check_constraints (props, subs, e))
4444 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4445 gsl_matrix_set (m, 0, 0, f (subs[0]));
4450 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props,
4451 gsl_matrix *subs[], const struct matrix_expr *e,
4452 matrix_proto_m_any *f)
4454 return check_constraints (props, subs, e) ? f (subs, e->n_subs) : NULL;
4458 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props_ UNUSED,
4459 gsl_matrix *subs[], const struct matrix_expr *e,
4460 matrix_proto_IDENT *f)
4462 static const struct matrix_function_properties p1 = {
4464 .constraints = "ai>=0"
4466 static const struct matrix_function_properties p2 = {
4468 .constraints = "ai>=0 bi>=0"
4470 const struct matrix_function_properties *props = e->n_subs == 1 ? &p1 : &p2;
4472 assert (e->n_subs <= 2);
4475 return (to_scalar_args (props->name, subs, e, d)
4476 && check_constraints (props, subs, e)
4477 ? f (d[0], d[e->n_subs - 1])
4482 matrix_expr_evaluate (const struct matrix_expr *e)
4484 if (e->op == MOP_NUMBER)
4486 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4487 gsl_matrix_set (m, 0, 0, e->number);
4490 else if (e->op == MOP_VARIABLE)
4492 const gsl_matrix *src = e->variable->value;
4495 msg_at (SE, e->location,
4496 _("Uninitialized variable %s used in expression."),
4501 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4502 gsl_matrix_memcpy (dst, src);
4505 else if (e->op == MOP_EOF)
4507 struct dfm_reader *reader = read_file_open (e->eof);
4508 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4509 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4513 enum { N_LOCAL = 3 };
4514 gsl_matrix *local_subs[N_LOCAL];
4515 gsl_matrix **subs = (e->n_subs < N_LOCAL
4517 : xmalloc (e->n_subs * sizeof *subs));
4519 for (size_t i = 0; i < e->n_subs; i++)
4521 subs[i] = matrix_expr_evaluate (e->subs[i]);
4524 for (size_t j = 0; j < i; j++)
4525 gsl_matrix_free (subs[j]);
4526 if (subs != local_subs)
4532 gsl_matrix *result = NULL;
4535 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4536 case MOP_F_##ENUM: \
4538 static const struct matrix_function_properties props = { \
4540 .constraints = CONSTRAINTS, \
4542 result = matrix_expr_evaluate_##PROTO (&props, subs, e, \
4543 matrix_eval_##ENUM); \
4550 gsl_matrix_scale (subs[0], -1.0);
4568 result = matrix_expr_evaluate_elementwise (e, e->op, subs[0], subs[1]);
4572 result = matrix_expr_evaluate_not (subs[0]);
4576 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], NULL);
4580 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], subs[2]);
4584 result = matrix_expr_evaluate_mul_mat (e, subs[0], subs[1]);
4588 result = matrix_expr_evaluate_exp_mat (e, subs[0], subs[1]);
4591 case MOP_PASTE_HORZ:
4592 result = matrix_expr_evaluate_paste_horz (e, subs[0], subs[1]);
4595 case MOP_PASTE_VERT:
4596 result = matrix_expr_evaluate_paste_vert (e, subs[0], subs[1]);
4600 result = gsl_matrix_alloc (0, 0);
4604 result = matrix_expr_evaluate_vec_index (e, subs[0], subs[1]);
4608 result = matrix_expr_evaluate_vec_all (e, subs[0]);
4612 result = matrix_expr_evaluate_mat_index (subs[0],
4613 subs[1], e->subs[1],
4614 subs[2], e->subs[2]);
4618 result = matrix_expr_evaluate_mat_index (subs[0],
4619 subs[1], e->subs[1],
4624 result = matrix_expr_evaluate_mat_index (subs[0],
4626 subs[1], e->subs[1]);
4635 for (size_t i = 0; i < e->n_subs; i++)
4636 if (subs[i] != result)
4637 gsl_matrix_free (subs[i]);
4638 if (subs != local_subs)
4644 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4647 gsl_matrix *m = matrix_expr_evaluate (e);
4653 msg_at (SE, matrix_expr_location (e),
4654 _("Expression for %s must evaluate to scalar, "
4655 "not a %zu×%zu matrix."),
4656 context, m->size1, m->size2);
4657 gsl_matrix_free (m);
4662 gsl_matrix_free (m);
4667 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4671 if (!matrix_expr_evaluate_scalar (e, context, &d))
4675 if (d < LONG_MIN || d > LONG_MAX)
4677 msg_at (SE, matrix_expr_location (e),
4678 _("Expression for %s is outside the integer range."), context);
4687 An lvalue is an expression that can appear on the left side of a COMPUTE
4688 command and in other contexts that assign values.
4690 An lvalue is parsed once, with matrix_lvalue_parse(). It can then be
4691 evaluated (with matrix_lvalue_evaluate()) and assigned (with
4692 matrix_lvalue_assign()).
4694 There are three kinds of lvalues:
4696 - A variable name. A variable used as an lvalue need not be initialized,
4697 since the assignment will initialize.
4699 - A subvector, e.g. "var(index0)". The variable must be initialized and
4700 must have the form of a vector (it must have 1 column or 1 row).
4702 - A submatrix, e.g. "var(index0, index1)". The variable must be
4704 struct matrix_lvalue
4706 struct matrix_var *var; /* Destination variable. */
4707 struct matrix_expr *indexes[2]; /* Index expressions, if any. */
4708 size_t n_indexes; /* Number of indexes. */
4710 struct msg_location *var_location; /* Variable name. */
4711 struct msg_location *full_location; /* Variable name plus indexing. */
4712 struct msg_location *index_locations[2]; /* Index expressions. */
4717 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4721 msg_location_destroy (lvalue->var_location);
4722 msg_location_destroy (lvalue->full_location);
4723 for (size_t i = 0; i < lvalue->n_indexes; i++)
4725 matrix_expr_destroy (lvalue->indexes[i]);
4726 msg_location_destroy (lvalue->index_locations[i]);
4732 /* Parses and returns an lvalue at the current position in S's lexer. Returns
4733 null on parse failure. On success, the caller must eventually free the
4735 static struct matrix_lvalue *
4736 matrix_lvalue_parse (struct matrix_state *s)
4738 if (!lex_force_id (s->lexer))
4741 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4742 int start_ofs = lex_ofs (s->lexer);
4743 lvalue->var_location = lex_get_location (s->lexer, 0, 0);
4744 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4745 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4749 lex_error (s->lexer, _("Undefined variable %s."),
4750 lex_tokcstr (s->lexer));
4754 lex_get_n (s->lexer, 2);
4756 if (!matrix_parse_index_expr (s, &lvalue->indexes[0],
4757 &lvalue->index_locations[0]))
4759 lvalue->n_indexes++;
4761 if (lex_match (s->lexer, T_COMMA))
4763 if (!matrix_parse_index_expr (s, &lvalue->indexes[1],
4764 &lvalue->index_locations[1]))
4766 lvalue->n_indexes++;
4768 if (!lex_force_match (s->lexer, T_RPAREN))
4771 lvalue->full_location = lex_ofs_location (s->lexer, start_ofs,
4772 lex_ofs (s->lexer) - 1);
4777 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4783 matrix_lvalue_destroy (lvalue);
4788 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4789 enum index_type index_type, size_t other_size,
4790 struct index_vector *iv)
4795 m = matrix_expr_evaluate (e);
4802 bool ok = matrix_normalize_index_vector (m, e, size, index_type,
4804 gsl_matrix_free (m);
4808 /* Evaluates the indexes in LVALUE into IV0 and IV1, owned by the caller.
4809 Returns true if successful, false if evaluating the expressions failed or if
4810 LVALUE otherwise can't be used for an assignment.
4812 On success, the caller retains ownership of the index vectors, which are
4813 suitable for passing to matrix_lvalue_assign(). If not used for that
4814 purpose then they need to eventually be freed (with
4815 index_vector_uninit()). */
4817 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4818 struct index_vector *iv0,
4819 struct index_vector *iv1)
4821 *iv0 = INDEX_VECTOR_INIT;
4822 *iv1 = INDEX_VECTOR_INIT;
4823 if (!lvalue->n_indexes)
4826 /* Validate destination matrix exists and has the right shape. */
4827 gsl_matrix *dm = lvalue->var->value;
4830 msg_at (SE, lvalue->var_location,
4831 _("Undefined variable %s."), lvalue->var->name);
4834 else if (dm->size1 == 0 || dm->size2 == 0)
4836 msg_at (SE, lvalue->full_location, _("Cannot index %zu×%zu matrix %s."),
4837 dm->size1, dm->size2, lvalue->var->name);
4840 else if (lvalue->n_indexes == 1)
4842 if (!is_vector (dm))
4844 msg_at (SE, lvalue->full_location,
4845 _("Can't use vector indexing on %zu×%zu matrix %s."),
4846 dm->size1, dm->size2, lvalue->var->name);
4849 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4850 MAX (dm->size1, dm->size2),
4855 assert (lvalue->n_indexes == 2);
4856 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4857 IV_ROW, dm->size2, iv0))
4860 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4861 IV_COLUMN, dm->size1, iv1))
4863 index_vector_uninit (iv0);
4871 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4872 struct index_vector *iv,
4873 gsl_matrix *sm, const struct msg_location *lsm)
4875 /* Convert source matrix 'sm' to source vector 'sv'. */
4876 if (!is_vector (sm))
4878 msg_at (SE, lvalue->full_location,
4879 _("Only an %zu-element vector may be assigned to this "
4880 "%zu-element subvector of %s."),
4881 iv->n, iv->n, lvalue->var->name);
4883 _("The source is an %zu×%zu matrix."),
4884 sm->size1, sm->size2);
4887 gsl_vector sv = to_vector (sm);
4888 if (iv->n != sv.size)
4890 msg_at (SE, lvalue->full_location,
4891 _("Only an %zu-element vector may be assigned to this "
4892 "%zu-element subvector of %s."),
4893 iv->n, iv->n, lvalue->var->name);
4894 msg_at (SE, lsm, ngettext ("The source vector has %zu element.",
4895 "The source vector has %zu elements.",
4901 /* Assign elements. */
4902 gsl_vector dv = to_vector (lvalue->var->value);
4903 for (size_t x = 0; x < iv->n; x++)
4904 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4909 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4910 struct index_vector *iv0,
4911 struct index_vector *iv1,
4912 gsl_matrix *sm, const struct msg_location *lsm)
4914 gsl_matrix *dm = lvalue->var->value;
4916 /* Convert source matrix 'sm' to source vector 'sv'. */
4917 bool wrong_rows = iv0->n != sm->size1;
4918 bool wrong_columns = iv1->n != sm->size2;
4919 if (wrong_rows || wrong_columns)
4921 if (wrong_rows && wrong_columns)
4922 msg_at (SE, lvalue->full_location,
4923 _("Numbers of indexes for assigning to %s differ from the "
4924 "size of the source matrix."),
4926 else if (wrong_rows)
4927 msg_at (SE, lvalue->full_location,
4928 _("Number of row indexes for assigning to %s differs from "
4929 "number of rows in source matrix."),
4932 msg_at (SE, lvalue->full_location,
4933 _("Number of column indexes for assigning to %s differs from "
4934 "number of columns in source matrix."),
4939 if (lvalue->indexes[0])
4940 msg_at (SN, lvalue->index_locations[0],
4941 ngettext ("There is %zu row index.",
4942 "There are %zu row indexes.",
4946 msg_at (SN, lvalue->index_locations[0],
4947 ngettext ("Destination matrix %s has %zu row.",
4948 "Destination matrix %s has %zu rows.",
4950 lvalue->var->name, iv0->n);
4955 if (lvalue->indexes[1])
4956 msg_at (SN, lvalue->index_locations[1],
4957 ngettext ("There is %zu column index.",
4958 "There are %zu column indexes.",
4962 msg_at (SN, lvalue->index_locations[1],
4963 ngettext ("Destination matrix %s has %zu column.",
4964 "Destination matrix %s has %zu columns.",
4966 lvalue->var->name, iv1->n);
4969 msg_at (SN, lsm, _("The source matrix is %zu×%zu."),
4970 sm->size1, sm->size2);
4974 /* Assign elements. */
4975 for (size_t y = 0; y < iv0->n; y++)
4977 size_t f0 = iv0->indexes[y];
4978 for (size_t x = 0; x < iv1->n; x++)
4980 size_t f1 = iv1->indexes[x];
4981 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4987 /* Assigns rvalue SM to LVALUE, whose index vectors IV0 and IV1 were previously
4988 obtained with matrix_lvalue_evaluate(). Returns true if successful, false
4989 on error. Always takes ownership of M. LSM should be the source location
4990 to use for errors related to SM. */
4992 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4993 struct index_vector *iv0, struct index_vector *iv1,
4994 gsl_matrix *sm, const struct msg_location *lsm)
4996 if (!lvalue->n_indexes)
4998 gsl_matrix_free (lvalue->var->value);
4999 lvalue->var->value = sm;
5004 bool ok = (lvalue->n_indexes == 1
5005 ? matrix_lvalue_assign_vector (lvalue, iv0, sm, lsm)
5006 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm, lsm));
5007 index_vector_uninit (iv0);
5008 index_vector_uninit (iv1);
5009 gsl_matrix_free (sm);
5014 /* Evaluates and then assigns SM to LVALUE. Always takes ownership of M. LSM
5015 should be the source location to use for errors related to SM.*/
5017 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue,
5019 const struct msg_location *lsm)
5021 struct index_vector iv0, iv1;
5022 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
5024 gsl_matrix_free (sm);
5028 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm, lsm);
5031 /* Matrix command data structure. */
5033 /* An array of matrix commands. */
5034 struct matrix_commands
5036 struct matrix_command **commands;
5040 static bool matrix_commands_parse (struct matrix_state *,
5041 struct matrix_commands *,
5042 const char *command_name,
5043 const char *stop1, const char *stop2);
5044 static void matrix_commands_uninit (struct matrix_commands *);
5046 /* A single matrix command. */
5047 struct matrix_command
5049 /* The type of command. */
5050 enum matrix_command_type
5071 /* Source lines for this command. */
5072 struct msg_location *location;
5076 struct matrix_compute
5078 struct matrix_lvalue *lvalue;
5079 struct matrix_expr *rvalue;
5085 struct matrix_expr *expression;
5086 bool use_default_format;
5087 struct fmt_spec format;
5089 int space; /* -1 means NEWPAGE. */
5093 struct string_array literals; /* CLABELS/RLABELS. */
5094 struct matrix_expr *expr; /* CNAMES/RNAMES. */
5104 struct matrix_expr *condition;
5105 struct matrix_commands commands;
5115 struct matrix_var *var;
5116 struct matrix_expr *start;
5117 struct matrix_expr *end;
5118 struct matrix_expr *increment;
5120 /* Loop conditions. */
5121 struct matrix_expr *top_condition;
5122 struct matrix_expr *bottom_condition;
5125 struct matrix_commands commands;
5129 struct matrix_display
5131 struct matrix_state *state;
5135 struct matrix_release
5137 struct matrix_var **vars;
5144 struct matrix_expr *expression;
5145 struct save_file *sf;
5151 struct read_file *rf;
5152 struct matrix_lvalue *dst;
5153 struct matrix_expr *size;
5155 enum fmt_type format;
5164 struct write_file *wf;
5165 struct matrix_expr *expression;
5168 /* If this is nonnull, WRITE uses this format.
5170 If this is NULL, WRITE uses free-field format with as many
5171 digits of precision as needed. */
5172 struct fmt_spec *format;
5181 struct lexer *lexer;
5182 struct matrix_lvalue *dst;
5183 struct dataset *dataset;
5184 struct file_handle *file;
5186 struct var_syntax *vars;
5188 struct matrix_var *names;
5190 /* Treatment of missing values. */
5195 MGET_ERROR, /* Flag error on command. */
5196 MGET_ACCEPT, /* Accept without change, user-missing only. */
5197 MGET_OMIT, /* Drop this case. */
5198 MGET_RECODE /* Recode to 'substitute'. */
5201 double substitute; /* MGET_RECODE only. */
5209 struct msave_common *common;
5210 struct matrix_expr *expr;
5211 const char *rowtype;
5212 const struct matrix_expr *factors;
5213 const struct matrix_expr *splits;
5219 struct matrix_state *state;
5220 struct file_handle *file;
5222 struct stringi_set rowtypes;
5228 struct matrix_expr *expr;
5229 struct matrix_var *evec;
5230 struct matrix_var *eval;
5234 struct matrix_setdiag
5236 struct matrix_var *dst;
5237 struct matrix_expr *expr;
5243 struct matrix_expr *expr;
5244 struct matrix_var *u;
5245 struct matrix_var *s;
5246 struct matrix_var *v;
5252 static struct matrix_command *matrix_command_parse (struct matrix_state *);
5253 static bool matrix_command_execute (struct matrix_command *);
5254 static void matrix_command_destroy (struct matrix_command *);
5258 static struct matrix_command *
5259 matrix_compute_parse (struct matrix_state *s)
5261 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5262 *cmd = (struct matrix_command) {
5263 .type = MCMD_COMPUTE,
5264 .compute = { .lvalue = NULL },
5267 struct matrix_compute *compute = &cmd->compute;
5268 compute->lvalue = matrix_lvalue_parse (s);
5269 if (!compute->lvalue)
5272 if (!lex_force_match (s->lexer, T_EQUALS))
5275 compute->rvalue = matrix_expr_parse (s);
5276 if (!compute->rvalue)
5282 matrix_command_destroy (cmd);
5287 matrix_compute_execute (struct matrix_command *cmd)
5289 struct matrix_compute *compute = &cmd->compute;
5290 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
5294 matrix_lvalue_evaluate_and_assign (compute->lvalue, value,
5295 matrix_expr_location (compute->rvalue));
5301 print_labels_destroy (struct print_labels *labels)
5305 string_array_destroy (&labels->literals);
5306 matrix_expr_destroy (labels->expr);
5312 parse_literal_print_labels (struct matrix_state *s,
5313 struct print_labels **labelsp)
5315 lex_match (s->lexer, T_EQUALS);
5316 print_labels_destroy (*labelsp);
5317 *labelsp = xzalloc (sizeof **labelsp);
5318 while (lex_token (s->lexer) != T_SLASH
5319 && lex_token (s->lexer) != T_ENDCMD
5320 && lex_token (s->lexer) != T_STOP)
5322 struct string label = DS_EMPTY_INITIALIZER;
5323 while (lex_token (s->lexer) != T_COMMA
5324 && lex_token (s->lexer) != T_SLASH
5325 && lex_token (s->lexer) != T_ENDCMD
5326 && lex_token (s->lexer) != T_STOP)
5328 if (!ds_is_empty (&label))
5329 ds_put_byte (&label, ' ');
5331 if (lex_token (s->lexer) == T_STRING)
5332 ds_put_cstr (&label, lex_tokcstr (s->lexer));
5335 char *rep = lex_next_representation (s->lexer, 0, 0);
5336 ds_put_cstr (&label, rep);
5341 string_array_append_nocopy (&(*labelsp)->literals,
5342 ds_steal_cstr (&label));
5344 if (!lex_match (s->lexer, T_COMMA))
5350 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
5352 lex_match (s->lexer, T_EQUALS);
5353 struct matrix_expr *e = matrix_parse_exp (s);
5357 print_labels_destroy (*labelsp);
5358 *labelsp = xzalloc (sizeof **labelsp);
5359 (*labelsp)->expr = e;
5363 static struct matrix_command *
5364 matrix_print_parse (struct matrix_state *s)
5366 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5367 *cmd = (struct matrix_command) {
5370 .use_default_format = true,
5374 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
5376 int start_ofs = lex_ofs (s->lexer);
5377 cmd->print.expression = matrix_parse_exp (s);
5378 if (!cmd->print.expression)
5380 cmd->print.title = lex_ofs_representation (s->lexer, start_ofs,
5381 lex_ofs (s->lexer) - 1);
5384 while (lex_match (s->lexer, T_SLASH))
5386 if (lex_match_id (s->lexer, "FORMAT"))
5388 lex_match (s->lexer, T_EQUALS);
5389 if (!parse_format_specifier (s->lexer, &cmd->print.format))
5391 cmd->print.use_default_format = false;
5393 else if (lex_match_id (s->lexer, "TITLE"))
5395 lex_match (s->lexer, T_EQUALS);
5396 if (!lex_force_string (s->lexer))
5398 free (cmd->print.title);
5399 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
5402 else if (lex_match_id (s->lexer, "SPACE"))
5404 lex_match (s->lexer, T_EQUALS);
5405 if (lex_match_id (s->lexer, "NEWPAGE"))
5406 cmd->print.space = -1;
5407 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
5409 cmd->print.space = lex_integer (s->lexer);
5415 else if (lex_match_id (s->lexer, "RLABELS"))
5416 parse_literal_print_labels (s, &cmd->print.rlabels);
5417 else if (lex_match_id (s->lexer, "CLABELS"))
5418 parse_literal_print_labels (s, &cmd->print.clabels);
5419 else if (lex_match_id (s->lexer, "RNAMES"))
5421 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
5424 else if (lex_match_id (s->lexer, "CNAMES"))
5426 if (!parse_expr_print_labels (s, &cmd->print.clabels))
5431 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
5432 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
5440 matrix_command_destroy (cmd);
5445 matrix_is_integer (const gsl_matrix *m)
5447 for (size_t y = 0; y < m->size1; y++)
5448 for (size_t x = 0; x < m->size2; x++)
5450 double d = gsl_matrix_get (m, y, x);
5458 matrix_max_magnitude (const gsl_matrix *m)
5461 for (size_t y = 0; y < m->size1; y++)
5462 for (size_t x = 0; x < m->size2; x++)
5464 double d = fabs (gsl_matrix_get (m, y, x));
5472 format_fits (struct fmt_spec format, double x)
5474 char *s = data_out (&(union value) { .f = x }, NULL,
5475 &format, settings_get_fmt_settings ());
5476 bool fits = *s != '*' && !strchr (s, 'E');
5481 static struct fmt_spec
5482 default_format (const gsl_matrix *m, int *log_scale)
5486 double max = matrix_max_magnitude (m);
5488 if (matrix_is_integer (m))
5489 for (int w = 1; w <= 12; w++)
5491 struct fmt_spec format = { .type = FMT_F, .w = w };
5492 if (format_fits (format, -max))
5496 if (max >= 1e9 || max <= 1e-4)
5499 snprintf (s, sizeof s, "%.1e", max);
5501 const char *e = strchr (s, 'e');
5503 *log_scale = atoi (e + 1);
5506 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5510 trimmed_string (double d)
5512 struct substring s = ss_buffer ((char *) &d, sizeof d);
5513 ss_rtrim (&s, ss_cstr (" "));
5514 return ss_xstrdup (s);
5517 static struct string_array *
5518 print_labels_get (const struct print_labels *labels, size_t n,
5519 const char *prefix, bool truncate)
5524 struct string_array *out = xzalloc (sizeof *out);
5525 if (labels->literals.n)
5526 string_array_clone (out, &labels->literals);
5527 else if (labels->expr)
5529 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5530 if (m && is_vector (m))
5532 gsl_vector v = to_vector (m);
5533 for (size_t i = 0; i < v.size; i++)
5534 string_array_append_nocopy (out, trimmed_string (
5535 gsl_vector_get (&v, i)));
5537 gsl_matrix_free (m);
5543 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5546 string_array_append (out, "");
5549 string_array_delete (out, out->n - 1);
5552 for (size_t i = 0; i < out->n; i++)
5554 char *s = out->strings[i];
5555 s[strnlen (s, 8)] = '\0';
5562 matrix_print_space (int space)
5565 output_item_submit (page_break_item_create ());
5566 for (int i = 0; i < space; i++)
5567 output_log ("%s", "");
5571 matrix_print_text (const struct matrix_print *print, const gsl_matrix *m,
5572 struct fmt_spec format, int log_scale)
5574 matrix_print_space (print->space);
5576 output_log ("%s", print->title);
5578 output_log (" 10 ** %d X", log_scale);
5580 struct string_array *clabels = print_labels_get (print->clabels,
5581 m->size2, "col", true);
5582 if (clabels && format.w < 8)
5584 struct string_array *rlabels = print_labels_get (print->rlabels,
5585 m->size1, "row", true);
5589 struct string line = DS_EMPTY_INITIALIZER;
5591 ds_put_byte_multiple (&line, ' ', 8);
5592 for (size_t x = 0; x < m->size2; x++)
5593 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5594 output_log_nocopy (ds_steal_cstr (&line));
5597 double scale = pow (10.0, log_scale);
5598 bool numeric = fmt_is_numeric (format.type);
5599 for (size_t y = 0; y < m->size1; y++)
5601 struct string line = DS_EMPTY_INITIALIZER;
5603 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5605 for (size_t x = 0; x < m->size2; x++)
5607 double f = gsl_matrix_get (m, y, x);
5609 ? data_out (&(union value) { .f = f / scale}, NULL,
5610 &format, settings_get_fmt_settings ())
5611 : trimmed_string (f));
5612 ds_put_format (&line, " %s", s);
5615 output_log_nocopy (ds_steal_cstr (&line));
5618 string_array_destroy (rlabels);
5620 string_array_destroy (clabels);
5625 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5626 const struct print_labels *print_labels, size_t n,
5627 const char *name, const char *prefix)
5629 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5631 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5632 for (size_t i = 0; i < n; i++)
5633 pivot_category_create_leaf (
5635 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5636 : pivot_value_new_integer (i + 1)));
5638 d->hide_all_labels = true;
5639 string_array_destroy (labels);
5644 matrix_print_tables (const struct matrix_print *print, const gsl_matrix *m,
5645 struct fmt_spec format, int log_scale)
5647 struct pivot_table *table = pivot_table_create__ (
5648 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5651 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5653 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5654 N_("Columns"), "col");
5656 struct pivot_footnote *footnote = NULL;
5659 char *s = xasprintf ("× 10**%d", log_scale);
5660 footnote = pivot_table_create_footnote (
5661 table, pivot_value_new_user_text_nocopy (s));
5664 double scale = pow (10.0, log_scale);
5665 bool numeric = fmt_is_numeric (format.type);
5666 for (size_t y = 0; y < m->size1; y++)
5667 for (size_t x = 0; x < m->size2; x++)
5669 double f = gsl_matrix_get (m, y, x);
5670 struct pivot_value *v;
5673 v = pivot_value_new_number (f / scale);
5674 v->numeric.format = format;
5677 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5679 pivot_value_add_footnote (v, footnote);
5680 pivot_table_put2 (table, y, x, v);
5683 pivot_table_submit (table);
5687 matrix_print_execute (const struct matrix_print *print)
5689 if (print->expression)
5691 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5696 struct fmt_spec format = (print->use_default_format
5697 ? default_format (m, &log_scale)
5700 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5701 matrix_print_text (print, m, format, log_scale);
5703 matrix_print_tables (print, m, format, log_scale);
5705 gsl_matrix_free (m);
5709 matrix_print_space (print->space);
5711 output_log ("%s", print->title);
5718 matrix_do_if_clause_parse (struct matrix_state *s,
5719 struct matrix_do_if *ifc,
5720 bool parse_condition,
5721 size_t *allocated_clauses)
5723 if (ifc->n_clauses >= *allocated_clauses)
5724 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5725 sizeof *ifc->clauses);
5726 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5727 *c = (struct do_if_clause) { .condition = NULL };
5729 if (parse_condition)
5731 c->condition = matrix_expr_parse (s);
5736 return matrix_commands_parse (s, &c->commands, "DO IF", "ELSE", "END IF");
5739 static struct matrix_command *
5740 matrix_do_if_parse (struct matrix_state *s)
5742 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5743 *cmd = (struct matrix_command) {
5745 .do_if = { .n_clauses = 0 }
5748 size_t allocated_clauses = 0;
5751 if (!matrix_do_if_clause_parse (s, &cmd->do_if, true, &allocated_clauses))
5754 while (lex_match_phrase (s->lexer, "ELSE IF"));
5756 if (lex_match_id (s->lexer, "ELSE")
5757 && !matrix_do_if_clause_parse (s, &cmd->do_if, false, &allocated_clauses))
5760 if (!lex_match_phrase (s->lexer, "END IF"))
5765 matrix_command_destroy (cmd);
5770 matrix_do_if_execute (struct matrix_do_if *cmd)
5772 for (size_t i = 0; i < cmd->n_clauses; i++)
5774 struct do_if_clause *c = &cmd->clauses[i];
5778 if (!matrix_expr_evaluate_scalar (c->condition,
5779 i ? "ELSE IF" : "DO IF",
5784 for (size_t j = 0; j < c->commands.n; j++)
5785 if (!matrix_command_execute (c->commands.commands[j]))
5794 static struct matrix_command *
5795 matrix_loop_parse (struct matrix_state *s)
5797 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5798 *cmd = (struct matrix_command) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5800 struct matrix_loop *loop = &cmd->loop;
5801 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5803 struct substring name = lex_tokss (s->lexer);
5804 loop->var = matrix_var_lookup (s, name);
5806 loop->var = matrix_var_create (s, name);
5811 loop->start = matrix_expr_parse (s);
5812 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5815 loop->end = matrix_expr_parse (s);
5819 if (lex_match (s->lexer, T_BY))
5821 loop->increment = matrix_expr_parse (s);
5822 if (!loop->increment)
5827 if (lex_match_id (s->lexer, "IF"))
5829 loop->top_condition = matrix_expr_parse (s);
5830 if (!loop->top_condition)
5834 bool was_in_loop = s->in_loop;
5836 bool ok = matrix_commands_parse (s, &loop->commands, "LOOP",
5838 s->in_loop = was_in_loop;
5842 if (!lex_match_phrase (s->lexer, "END LOOP"))
5845 if (lex_match_id (s->lexer, "IF"))
5847 loop->bottom_condition = matrix_expr_parse (s);
5848 if (!loop->bottom_condition)
5855 matrix_command_destroy (cmd);
5860 matrix_loop_execute (struct matrix_loop *cmd)
5864 long int increment = 1;
5867 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5868 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5870 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5874 if (increment > 0 ? value > end
5875 : increment < 0 ? value < end
5880 int mxloops = settings_get_mxloops ();
5881 for (int i = 0; i < mxloops; i++)
5886 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5888 gsl_matrix_free (cmd->var->value);
5889 cmd->var->value = NULL;
5891 if (!cmd->var->value)
5892 cmd->var->value = gsl_matrix_alloc (1, 1);
5894 gsl_matrix_set (cmd->var->value, 0, 0, value);
5897 if (cmd->top_condition)
5900 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5905 for (size_t j = 0; j < cmd->commands.n; j++)
5906 if (!matrix_command_execute (cmd->commands.commands[j]))
5909 if (cmd->bottom_condition)
5912 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5913 "END LOOP IF", &d) || d > 0)
5920 if (increment > 0 ? value > end : value < end)
5928 static struct matrix_command *
5929 matrix_break_parse (struct matrix_state *s)
5933 lex_next_error (s->lexer, -1, -1, _("BREAK not inside LOOP."));
5937 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5938 *cmd = (struct matrix_command) { .type = MCMD_BREAK };
5944 static struct matrix_command *
5945 matrix_display_parse (struct matrix_state *s)
5947 while (lex_token (s->lexer) != T_ENDCMD)
5949 if (!lex_match_id (s->lexer, "DICTIONARY")
5950 && !lex_match_id (s->lexer, "STATUS"))
5952 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5957 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5958 *cmd = (struct matrix_command) { .type = MCMD_DISPLAY, .display = { s } };
5963 compare_matrix_var_pointers (const void *a_, const void *b_)
5965 const struct matrix_var *const *ap = a_;
5966 const struct matrix_var *const *bp = b_;
5967 const struct matrix_var *a = *ap;
5968 const struct matrix_var *b = *bp;
5969 return strcmp (a->name, b->name);
5973 matrix_display_execute (struct matrix_display *cmd)
5975 const struct matrix_state *s = cmd->state;
5977 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5978 struct pivot_dimension *attr_dimension
5979 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5980 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5981 N_("Rows"), N_("Columns"));
5982 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5984 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5987 const struct matrix_var *var;
5988 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5990 vars[n_vars++] = var;
5991 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5993 struct pivot_dimension *rows = pivot_dimension_create (
5994 table, PIVOT_AXIS_ROW, N_("Variable"));
5995 for (size_t i = 0; i < n_vars; i++)
5997 const struct matrix_var *var = vars[i];
5998 pivot_category_create_leaf (
5999 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
6001 size_t r = var->value->size1;
6002 size_t c = var->value->size2;
6003 double values[] = { r, c, r * c * sizeof (double) / 1024 };
6004 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
6005 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
6008 pivot_table_submit (table);
6013 static struct matrix_command *
6014 matrix_release_parse (struct matrix_state *s)
6016 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6017 *cmd = (struct matrix_command) { .type = MCMD_RELEASE };
6019 size_t allocated_vars = 0;
6020 while (lex_token (s->lexer) == T_ID)
6022 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
6025 if (cmd->release.n_vars >= allocated_vars)
6026 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
6027 sizeof *cmd->release.vars);
6028 cmd->release.vars[cmd->release.n_vars++] = var;
6031 lex_error (s->lexer, _("Syntax error expecting variable name."));
6034 if (!lex_match (s->lexer, T_COMMA))
6042 matrix_release_execute (struct matrix_release *cmd)
6044 for (size_t i = 0; i < cmd->n_vars; i++)
6046 struct matrix_var *var = cmd->vars[i];
6047 gsl_matrix_free (var->value);
6054 static struct save_file *
6055 save_file_create (struct matrix_state *s, struct file_handle *fh,
6056 struct string_array *variables,
6057 struct matrix_expr *names,
6058 struct stringi_set *strings)
6060 for (size_t i = 0; i < s->n_save_files; i++)
6062 struct save_file *sf = s->save_files[i];
6063 if (fh_equal (sf->file, fh))
6067 string_array_destroy (variables);
6068 matrix_expr_destroy (names);
6069 stringi_set_destroy (strings);
6075 struct save_file *sf = xmalloc (sizeof *sf);
6076 *sf = (struct save_file) {
6078 .dataset = s->dataset,
6079 .variables = *variables,
6081 .strings = STRINGI_SET_INITIALIZER (sf->strings),
6084 stringi_set_swap (&sf->strings, strings);
6085 stringi_set_destroy (strings);
6087 s->save_files = xrealloc (s->save_files,
6088 (s->n_save_files + 1) * sizeof *s->save_files);
6089 s->save_files[s->n_save_files++] = sf;
6093 static struct casewriter *
6094 save_file_open (struct save_file *sf, gsl_matrix *m,
6095 const struct msg_location *save_location)
6097 if (sf->writer || sf->error)
6101 size_t n_variables = caseproto_get_n_widths (
6102 casewriter_get_proto (sf->writer));
6103 if (m->size2 != n_variables)
6105 const char *file_name = (sf->file == fh_inline_file () ? "*"
6106 : fh_get_name (sf->file));
6107 msg_at (SE, save_location,
6108 _("Cannot save %zu×%zu matrix to %s because the "
6109 "first SAVE to %s in this matrix program wrote a "
6110 "%zu-column matrix."),
6111 m->size1, m->size2, file_name, file_name, n_variables);
6112 msg_at (SE, sf->location,
6113 _("This is the location of the first SAVE to %s."),
6122 struct dictionary *dict = dict_create (get_default_encoding ());
6124 /* Fill 'names' with user-specified names if there were any, either from
6125 sf->variables or sf->names. */
6126 struct string_array names = { .n = 0 };
6127 if (sf->variables.n)
6128 string_array_clone (&names, &sf->variables);
6131 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
6132 if (nm && is_vector (nm))
6134 gsl_vector nv = to_vector (nm);
6135 for (size_t i = 0; i < nv.size; i++)
6137 char *name = trimmed_string (gsl_vector_get (&nv, i));
6138 char *error = dict_id_is_valid__ (dict, name);
6140 string_array_append_nocopy (&names, name);
6143 msg_at (SE, save_location, "%s", error);
6149 gsl_matrix_free (nm);
6152 struct stringi_set strings;
6153 stringi_set_clone (&strings, &sf->strings);
6155 for (size_t i = 0; dict_get_n_vars (dict) < m->size2; i++)
6160 name = names.strings[i];
6163 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
6167 int width = stringi_set_delete (&strings, name) ? 8 : 0;
6168 struct variable *var = dict_create_var (dict, name, width);
6171 msg_at (ME, save_location,
6172 _("Duplicate variable name %s in SAVE statement."), name);
6177 if (!stringi_set_is_empty (&strings))
6179 size_t count = stringi_set_count (&strings);
6180 const char *example = stringi_set_node_get_string (
6181 stringi_set_first (&strings));
6184 msg_at (ME, save_location,
6185 _("The SAVE command STRINGS subcommand specifies an unknown "
6186 "variable %s."), example);
6188 msg_at (ME, save_location,
6189 ngettext ("The SAVE command STRINGS subcommand specifies %zu "
6190 "unknown variable: %s.",
6191 "The SAVE command STRINGS subcommand specifies %zu "
6192 "unknown variables, including %s.",
6197 stringi_set_destroy (&strings);
6198 string_array_destroy (&names);
6207 if (sf->file == fh_inline_file ())
6208 sf->writer = autopaging_writer_create (dict_get_proto (dict));
6210 sf->writer = any_writer_open (sf->file, dict);
6214 sf->location = msg_location_dup (save_location);
6225 save_file_destroy (struct save_file *sf)
6229 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
6231 dataset_set_dict (sf->dataset, sf->dict);
6232 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
6236 casewriter_destroy (sf->writer);
6237 dict_unref (sf->dict);
6239 fh_unref (sf->file);
6240 string_array_destroy (&sf->variables);
6241 matrix_expr_destroy (sf->names);
6242 stringi_set_destroy (&sf->strings);
6243 msg_location_destroy (sf->location);
6248 static struct matrix_command *
6249 matrix_save_parse (struct matrix_state *s)
6251 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6252 *cmd = (struct matrix_command) {
6254 .save = { .expression = NULL },
6257 struct file_handle *fh = NULL;
6258 struct matrix_save *save = &cmd->save;
6260 struct string_array variables = STRING_ARRAY_INITIALIZER;
6261 struct matrix_expr *names = NULL;
6262 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
6264 save->expression = matrix_parse_exp (s);
6265 if (!save->expression)
6268 int names_start = 0;
6270 while (lex_match (s->lexer, T_SLASH))
6272 if (lex_match_id (s->lexer, "OUTFILE"))
6274 lex_match (s->lexer, T_EQUALS);
6277 fh = (lex_match (s->lexer, T_ASTERISK)
6279 : fh_parse (s->lexer, FH_REF_FILE, s->session));
6283 else if (lex_match_id (s->lexer, "VARIABLES"))
6285 lex_match (s->lexer, T_EQUALS);
6289 struct dictionary *d = dict_create (get_default_encoding ());
6290 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
6291 PV_NO_SCRATCH | PV_NO_DUPLICATE);
6296 string_array_clear (&variables);
6297 variables = (struct string_array) {
6303 else if (lex_match_id (s->lexer, "NAMES"))
6305 lex_match (s->lexer, T_EQUALS);
6306 matrix_expr_destroy (names);
6307 names_start = lex_ofs (s->lexer);
6308 names = matrix_parse_exp (s);
6309 names_end = lex_ofs (s->lexer) - 1;
6313 else if (lex_match_id (s->lexer, "STRINGS"))
6315 lex_match (s->lexer, T_EQUALS);
6316 while (lex_token (s->lexer) == T_ID)
6318 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
6320 if (!lex_match (s->lexer, T_COMMA))
6326 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
6334 if (s->prev_save_file)
6335 fh = fh_ref (s->prev_save_file);
6338 lex_sbc_missing (s->lexer, "OUTFILE");
6342 fh_unref (s->prev_save_file);
6343 s->prev_save_file = fh_ref (fh);
6345 if (variables.n && names)
6347 lex_ofs_msg (s->lexer, SW, names_start, names_end,
6348 _("Ignoring NAMES because VARIABLES was also specified."));
6349 matrix_expr_destroy (names);
6353 save->sf = save_file_create (s, fh, &variables, names, &strings);
6357 string_array_destroy (&variables);
6358 matrix_expr_destroy (names);
6359 stringi_set_destroy (&strings);
6361 matrix_command_destroy (cmd);
6366 matrix_save_execute (const struct matrix_command *cmd)
6368 const struct matrix_save *save = &cmd->save;
6370 gsl_matrix *m = matrix_expr_evaluate (save->expression);
6374 struct casewriter *writer = save_file_open (save->sf, m, cmd->location);
6377 gsl_matrix_free (m);
6381 const struct caseproto *proto = casewriter_get_proto (writer);
6382 for (size_t y = 0; y < m->size1; y++)
6384 struct ccase *c = case_create (proto);
6385 for (size_t x = 0; x < m->size2; x++)
6387 double d = gsl_matrix_get (m, y, x);
6388 int width = caseproto_get_width (proto, x);
6389 union value *value = case_data_rw_idx (c, x);
6393 memcpy (value->s, &d, width);
6395 casewriter_write (writer, c);
6397 gsl_matrix_free (m);
6402 static struct read_file *
6403 read_file_create (struct matrix_state *s, struct file_handle *fh)
6405 for (size_t i = 0; i < s->n_read_files; i++)
6407 struct read_file *rf = s->read_files[i];
6415 struct read_file *rf = xmalloc (sizeof *rf);
6416 *rf = (struct read_file) { .file = fh };
6418 s->read_files = xrealloc (s->read_files,
6419 (s->n_read_files + 1) * sizeof *s->read_files);
6420 s->read_files[s->n_read_files++] = rf;
6424 static struct dfm_reader *
6425 read_file_open (struct read_file *rf)
6428 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
6433 read_file_destroy (struct read_file *rf)
6437 fh_unref (rf->file);
6438 dfm_close_reader (rf->reader);
6439 free (rf->encoding);
6444 static struct matrix_command *
6445 matrix_read_parse (struct matrix_state *s)
6447 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6448 *cmd = (struct matrix_command) {
6450 .read = { .format = FMT_F },
6453 struct file_handle *fh = NULL;
6454 char *encoding = NULL;
6455 struct matrix_read *read = &cmd->read;
6456 read->dst = matrix_lvalue_parse (s);
6462 int record_width_start = 0, record_width_end = 0;
6465 int repetitions = 0;
6466 int record_width = 0;
6467 bool seen_format = false;
6468 while (lex_match (s->lexer, T_SLASH))
6470 if (lex_match_id (s->lexer, "FILE"))
6472 lex_match (s->lexer, T_EQUALS);
6475 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6479 else if (lex_match_id (s->lexer, "ENCODING"))
6481 lex_match (s->lexer, T_EQUALS);
6482 if (!lex_force_string (s->lexer))
6486 encoding = ss_xstrdup (lex_tokss (s->lexer));
6490 else if (lex_match_id (s->lexer, "FIELD"))
6492 lex_match (s->lexer, T_EQUALS);
6494 record_width_start = lex_ofs (s->lexer);
6495 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6497 read->c1 = lex_integer (s->lexer);
6499 if (!lex_force_match (s->lexer, T_TO)
6500 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
6502 read->c2 = lex_integer (s->lexer) + 1;
6503 record_width_end = lex_ofs (s->lexer);
6506 record_width = read->c2 - read->c1;
6507 if (lex_match (s->lexer, T_BY))
6509 if (!lex_force_int_range (s->lexer, "BY", 1,
6510 read->c2 - read->c1))
6512 by = lex_integer (s->lexer);
6513 by_ofs = lex_ofs (s->lexer);
6514 int field_end = lex_ofs (s->lexer);
6517 if (record_width % by)
6520 s->lexer, record_width_start, field_end,
6521 _("Field width %d does not evenly divide record width %d."),
6523 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
6524 _("This syntax designates the record width."));
6525 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
6526 _("This syntax specifies the field width."));
6533 else if (lex_match_id (s->lexer, "SIZE"))
6535 lex_match (s->lexer, T_EQUALS);
6536 matrix_expr_destroy (read->size);
6537 read->size = matrix_parse_exp (s);
6541 else if (lex_match_id (s->lexer, "MODE"))
6543 lex_match (s->lexer, T_EQUALS);
6544 if (lex_match_id (s->lexer, "RECTANGULAR"))
6545 read->symmetric = false;
6546 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6547 read->symmetric = true;
6550 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6554 else if (lex_match_id (s->lexer, "REREAD"))
6555 read->reread = true;
6556 else if (lex_match_id (s->lexer, "FORMAT"))
6560 lex_sbc_only_once (s->lexer, "FORMAT");
6565 lex_match (s->lexer, T_EQUALS);
6567 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6570 format_ofs = lex_ofs (s->lexer);
6571 const char *p = lex_tokcstr (s->lexer);
6572 if (c_isdigit (p[0]))
6574 repetitions = atoi (p);
6575 p += strspn (p, "0123456789");
6576 if (!fmt_from_name (p, &read->format))
6578 lex_error (s->lexer, _("Unknown format %s."), p);
6583 else if (fmt_from_name (p, &read->format))
6587 struct fmt_spec format;
6588 if (!parse_format_specifier (s->lexer, &format))
6590 read->format = format.type;
6596 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6597 "REREAD", "FORMAT");
6604 lex_sbc_missing (s->lexer, "FIELD");
6608 if (!read->dst->n_indexes && !read->size)
6610 msg (SE, _("SIZE is required for reading data into a full matrix "
6611 "(as opposed to a submatrix)."));
6612 msg_at (SN, read->dst->var_location,
6613 _("This expression designates a full matrix."));
6619 if (s->prev_read_file)
6620 fh = fh_ref (s->prev_read_file);
6623 lex_sbc_missing (s->lexer, "FILE");
6627 fh_unref (s->prev_read_file);
6628 s->prev_read_file = fh_ref (fh);
6630 read->rf = read_file_create (s, fh);
6634 free (read->rf->encoding);
6635 read->rf->encoding = encoding;
6639 /* Field width may be specified in multiple ways:
6642 2. The format on FORMAT.
6643 3. The repetition factor on FORMAT.
6645 (2) and (3) are mutually exclusive.
6647 If more than one of these is present, they must agree. If none of them is
6648 present, then free-field format is used.
6650 if (repetitions > record_width)
6652 msg (SE, _("%d repetitions cannot fit in record width %d."),
6653 repetitions, record_width);
6654 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
6655 _("This syntax designates the number of repetitions."));
6656 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
6657 _("This syntax designates the record width."));
6660 int w = (repetitions ? record_width / repetitions
6665 msg (SE, _("This command specifies two different field widths."));
6668 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
6669 ngettext ("This syntax specifies %d repetition.",
6670 "This syntax specifies %d repetitions.",
6673 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
6674 _("This syntax designates record width %d, "
6675 "which divided by %d repetitions implies "
6677 record_width, repetitions, w);
6680 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
6681 _("This syntax specifies field width %d."), w);
6683 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
6684 _("This syntax specifies field width %d."), by);
6692 matrix_command_destroy (cmd);
6698 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6699 struct substring data, size_t y, size_t x,
6700 int first_column, int last_column, char *error)
6702 int line_number = dfm_get_line_number (reader);
6703 struct msg_location location = {
6704 .file_name = intern_new (dfm_get_file_name (reader)),
6705 .start = { .line = line_number, .column = first_column },
6706 .end = { .line = line_number, .column = last_column },
6708 msg_at (DW, &location, _("Error reading \"%.*s\" as format %s "
6709 "for matrix row %zu, column %zu: %s"),
6710 (int) data.length, data.string, fmt_name (format),
6711 y + 1, x + 1, error);
6712 msg_location_uninit (&location);
6717 matrix_read_set_field (struct matrix_read *read, struct dfm_reader *reader,
6718 gsl_matrix *m, struct substring p, size_t y, size_t x,
6719 const char *line_start)
6721 const char *input_encoding = dfm_reader_get_encoding (reader);
6724 if (fmt_is_numeric (read->format))
6727 error = data_in (p, input_encoding, read->format,
6728 settings_get_fmt_settings (), &v, 0, NULL);
6729 if (!error && v.f == SYSMIS)
6730 error = xstrdup (_("Matrix data may not contain missing value."));
6735 uint8_t s[sizeof (double)];
6736 union value v = { .s = s };
6737 error = data_in (p, input_encoding, read->format,
6738 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6739 memcpy (&f, s, sizeof f);
6744 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6745 int nc = ss_utf8_count_columns (p);
6746 int c2 = c1 + MAX (1, nc) - 1;
6747 parse_error (reader, read->format, p, y, x, c1, c2, error);
6751 gsl_matrix_set (m, y, x, f);
6752 if (read->symmetric && x != y)
6753 gsl_matrix_set (m, x, y, f);
6758 matrix_read_line (struct matrix_command *cmd, struct dfm_reader *reader,
6759 struct substring *line, const char **startp)
6761 struct matrix_read *read = &cmd->read;
6762 if (dfm_eof (reader))
6764 msg_at (SE, cmd->location,
6765 _("Unexpected end of file reading matrix data."));
6768 dfm_expand_tabs (reader);
6769 struct substring record = dfm_get_record (reader);
6770 /* XXX need to recode record into UTF-8 */
6771 *startp = record.string;
6772 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6777 matrix_read (struct matrix_command *cmd, struct dfm_reader *reader,
6780 struct matrix_read *read = &cmd->read;
6781 for (size_t y = 0; y < m->size1; y++)
6783 size_t nx = read->symmetric ? y + 1 : m->size2;
6785 struct substring line = ss_empty ();
6786 const char *line_start = line.string;
6787 for (size_t x = 0; x < nx; x++)
6794 ss_ltrim (&line, ss_cstr (" ,"));
6795 if (!ss_is_empty (line))
6797 if (!matrix_read_line (cmd, reader, &line, &line_start))
6799 dfm_forward_record (reader);
6802 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6806 if (!matrix_read_line (cmd, reader, &line, &line_start))
6808 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6809 int f = x % fields_per_line;
6810 if (f == fields_per_line - 1)
6811 dfm_forward_record (reader);
6813 p = ss_substr (line, read->w * f, read->w);
6816 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6820 dfm_forward_record (reader);
6823 ss_ltrim (&line, ss_cstr (" ,"));
6824 if (!ss_is_empty (line))
6826 int line_number = dfm_get_line_number (reader);
6827 int c1 = utf8_count_columns (line_start,
6828 line.string - line_start) + 1;
6829 int c2 = c1 + ss_utf8_count_columns (line) - 1;
6830 struct msg_location location = {
6831 .file_name = intern_new (dfm_get_file_name (reader)),
6832 .start = { .line = line_number, .column = c1 },
6833 .end = { .line = line_number, .column = c2 },
6835 msg_at (DW, &location,
6836 _("Trailing garbage following data for matrix row %zu."),
6838 msg_location_uninit (&location);
6845 matrix_read_execute (struct matrix_command *cmd)
6847 struct matrix_read *read = &cmd->read;
6848 struct index_vector iv0, iv1;
6849 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6852 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6855 gsl_matrix *m = matrix_expr_evaluate (read->size);
6861 msg_at (SE, matrix_expr_location (read->size),
6862 _("SIZE must evaluate to a scalar or a 2-element vector, "
6863 "not a %zu×%zu matrix."), m->size1, m->size2);
6864 gsl_matrix_free (m);
6865 index_vector_uninit (&iv0);
6866 index_vector_uninit (&iv1);
6870 gsl_vector v = to_vector (m);
6874 d[0] = gsl_vector_get (&v, 0);
6877 else if (v.size == 2)
6879 d[0] = gsl_vector_get (&v, 0);
6880 d[1] = gsl_vector_get (&v, 1);
6884 msg_at (SE, matrix_expr_location (read->size),
6885 _("SIZE must evaluate to a scalar or a 2-element vector, "
6886 "not a %zu×%zu matrix."),
6887 m->size1, m->size2),
6888 gsl_matrix_free (m);
6889 index_vector_uninit (&iv0);
6890 index_vector_uninit (&iv1);
6893 gsl_matrix_free (m);
6895 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6897 msg_at (SE, matrix_expr_location (read->size),
6898 _("Matrix dimensions %g×%g specified on SIZE "
6899 "are outside valid range."),
6901 index_vector_uninit (&iv0);
6902 index_vector_uninit (&iv1);
6910 if (read->dst->n_indexes)
6912 size_t submatrix_size[2];
6913 if (read->dst->n_indexes == 2)
6915 submatrix_size[0] = iv0.n;
6916 submatrix_size[1] = iv1.n;
6918 else if (read->dst->var->value->size1 == 1)
6920 submatrix_size[0] = 1;
6921 submatrix_size[1] = iv0.n;
6925 submatrix_size[0] = iv0.n;
6926 submatrix_size[1] = 1;
6931 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6933 msg_at (SE, cmd->location,
6934 _("Dimensions specified on SIZE differ from dimensions "
6935 "of destination submatrix."));
6936 msg_at (SN, matrix_expr_location (read->size),
6937 _("SIZE specifies dimensions %zu×%zu."),
6939 msg_at (SN, read->dst->full_location,
6940 _("Destination submatrix has dimensions %zu×%zu."),
6941 submatrix_size[0], submatrix_size[1]);
6942 index_vector_uninit (&iv0);
6943 index_vector_uninit (&iv1);
6949 size[0] = submatrix_size[0];
6950 size[1] = submatrix_size[1];
6954 struct dfm_reader *reader = read_file_open (read->rf);
6956 dfm_reread_record (reader, 1);
6958 if (read->symmetric && size[0] != size[1])
6960 msg_at (SE, cmd->location,
6961 _("Cannot read non-square %zu×%zu matrix "
6962 "using READ with MODE=SYMMETRIC."),
6964 index_vector_uninit (&iv0);
6965 index_vector_uninit (&iv1);
6968 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6969 matrix_read (cmd, reader, tmp);
6970 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp, cmd->location);
6975 static struct write_file *
6976 write_file_create (struct matrix_state *s, struct file_handle *fh)
6978 for (size_t i = 0; i < s->n_write_files; i++)
6980 struct write_file *wf = s->write_files[i];
6988 struct write_file *wf = xmalloc (sizeof *wf);
6989 *wf = (struct write_file) { .file = fh };
6991 s->write_files = xrealloc (s->write_files,
6992 (s->n_write_files + 1) * sizeof *s->write_files);
6993 s->write_files[s->n_write_files++] = wf;
6997 static struct dfm_writer *
6998 write_file_open (struct write_file *wf)
7001 wf->writer = dfm_open_writer (wf->file, wf->encoding);
7006 write_file_destroy (struct write_file *wf)
7012 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
7013 wf->held->s.ss.length);
7014 u8_line_destroy (wf->held);
7018 fh_unref (wf->file);
7019 dfm_close_writer (wf->writer);
7020 free (wf->encoding);
7025 static struct matrix_command *
7026 matrix_write_parse (struct matrix_state *s)
7028 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7029 *cmd = (struct matrix_command) {
7033 struct file_handle *fh = NULL;
7034 char *encoding = NULL;
7035 struct matrix_write *write = &cmd->write;
7036 write->expression = matrix_parse_exp (s);
7037 if (!write->expression)
7042 int record_width_start = 0, record_width_end = 0;
7045 int repetitions = 0;
7046 int record_width = 0;
7047 enum fmt_type format = FMT_F;
7048 bool has_format = false;
7049 while (lex_match (s->lexer, T_SLASH))
7051 if (lex_match_id (s->lexer, "OUTFILE"))
7053 lex_match (s->lexer, T_EQUALS);
7056 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
7060 else if (lex_match_id (s->lexer, "ENCODING"))
7062 lex_match (s->lexer, T_EQUALS);
7063 if (!lex_force_string (s->lexer))
7067 encoding = ss_xstrdup (lex_tokss (s->lexer));
7071 else if (lex_match_id (s->lexer, "FIELD"))
7073 lex_match (s->lexer, T_EQUALS);
7075 record_width_start = lex_ofs (s->lexer);
7077 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
7079 write->c1 = lex_integer (s->lexer);
7081 if (!lex_force_match (s->lexer, T_TO)
7082 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
7084 write->c2 = lex_integer (s->lexer) + 1;
7085 record_width_end = lex_ofs (s->lexer);
7088 record_width = write->c2 - write->c1;
7089 if (lex_match (s->lexer, T_BY))
7091 if (!lex_force_int_range (s->lexer, "BY", 1,
7092 write->c2 - write->c1))
7094 by_ofs = lex_ofs (s->lexer);
7095 int field_end = lex_ofs (s->lexer);
7096 by = lex_integer (s->lexer);
7099 if (record_width % by)
7102 s->lexer, record_width_start, field_end,
7103 _("Field width %d does not evenly divide record width %d."),
7105 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7106 _("This syntax designates the record width."));
7107 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
7108 _("This syntax specifies the field width."));
7115 else if (lex_match_id (s->lexer, "MODE"))
7117 lex_match (s->lexer, T_EQUALS);
7118 if (lex_match_id (s->lexer, "RECTANGULAR"))
7119 write->triangular = false;
7120 else if (lex_match_id (s->lexer, "TRIANGULAR"))
7121 write->triangular = true;
7124 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
7128 else if (lex_match_id (s->lexer, "HOLD"))
7130 else if (lex_match_id (s->lexer, "FORMAT"))
7132 if (has_format || write->format)
7134 lex_sbc_only_once (s->lexer, "FORMAT");
7138 lex_match (s->lexer, T_EQUALS);
7140 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
7143 format_ofs = lex_ofs (s->lexer);
7144 const char *p = lex_tokcstr (s->lexer);
7145 if (c_isdigit (p[0]))
7147 repetitions = atoi (p);
7148 p += strspn (p, "0123456789");
7149 if (!fmt_from_name (p, &format))
7151 lex_error (s->lexer, _("Unknown format %s."), p);
7157 else if (fmt_from_name (p, &format))
7164 struct fmt_spec spec;
7165 if (!parse_format_specifier (s->lexer, &spec))
7167 write->format = xmemdup (&spec, sizeof spec);
7172 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
7180 lex_sbc_missing (s->lexer, "FIELD");
7186 if (s->prev_write_file)
7187 fh = fh_ref (s->prev_write_file);
7190 lex_sbc_missing (s->lexer, "OUTFILE");
7194 fh_unref (s->prev_write_file);
7195 s->prev_write_file = fh_ref (fh);
7197 write->wf = write_file_create (s, fh);
7201 free (write->wf->encoding);
7202 write->wf->encoding = encoding;
7206 /* Field width may be specified in multiple ways:
7209 2. The format on FORMAT.
7210 3. The repetition factor on FORMAT.
7212 (2) and (3) are mutually exclusive.
7214 If more than one of these is present, they must agree. If none of them is
7215 present, then free-field format is used.
7217 if (repetitions > record_width)
7219 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7220 _("This syntax designates the number of repetitions."));
7221 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7222 _("This syntax designates the record width."));
7225 int w = (repetitions ? record_width / repetitions
7226 : write->format ? write->format->w
7230 msg (SE, _("This command specifies two different field widths."));
7233 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7234 ngettext ("This syntax specifies %d repetition.",
7235 "This syntax specifies %d repetitions.",
7238 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7239 _("This syntax designates record width %d, "
7240 "which divided by %d repetitions implies "
7242 record_width, repetitions, w);
7245 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7246 _("This syntax specifies field width %d."), w);
7248 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
7249 _("This syntax specifies field width %d."), by);
7252 if (w && !write->format)
7254 write->format = xmalloc (sizeof *write->format);
7255 *write->format = (struct fmt_spec) { .type = format, .w = w };
7257 char *error = fmt_check_output__ (write->format);
7260 msg (SE, "%s", error);
7264 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7265 _("This syntax specifies format %s."),
7270 lex_ofs_msg (s->lexer, SN, format_ofs, format_ofs,
7271 ngettext ("This syntax specifies %d repetition.",
7272 "This syntax specifies %d repetitions.",
7275 lex_ofs_msg (s->lexer, SN, record_width_start, record_width_end,
7276 _("This syntax designates record width %d, "
7277 "which divided by %d repetitions implies "
7279 record_width, repetitions, w);
7283 lex_ofs_msg (s->lexer, SN, by_ofs, by_ofs,
7284 _("This syntax specifies field width %d."), by);
7290 if (write->format && fmt_var_width (write->format) > sizeof (double))
7292 char fs[FMT_STRING_LEN_MAX + 1];
7293 fmt_to_string (write->format, fs);
7294 lex_ofs_error (s->lexer, format_ofs, format_ofs,
7295 _("Format %s is too wide for %zu-byte matrix elements."),
7296 fs, sizeof (double));
7304 matrix_command_destroy (cmd);
7309 matrix_write_execute (struct matrix_write *write)
7311 gsl_matrix *m = matrix_expr_evaluate (write->expression);
7315 if (write->triangular && m->size1 != m->size2)
7317 msg_at (SE, matrix_expr_location (write->expression),
7318 _("WRITE with MODE=TRIANGULAR requires a square matrix but "
7319 "the matrix to be written has dimensions %zu×%zu."),
7320 m->size1, m->size2);
7321 gsl_matrix_free (m);
7325 struct dfm_writer *writer = write_file_open (write->wf);
7326 if (!writer || !m->size1)
7328 gsl_matrix_free (m);
7332 const struct fmt_settings *settings = settings_get_fmt_settings ();
7333 struct u8_line *line = write->wf->held;
7334 for (size_t y = 0; y < m->size1; y++)
7338 line = xmalloc (sizeof *line);
7339 u8_line_init (line);
7341 size_t nx = write->triangular ? y + 1 : m->size2;
7343 for (size_t x = 0; x < nx; x++)
7346 double f = gsl_matrix_get (m, y, x);
7350 if (fmt_is_numeric (write->format->type))
7353 v.s = (uint8_t *) &f;
7354 s = data_out (&v, NULL, write->format, settings);
7358 s = xmalloc (DBL_BUFSIZE_BOUND);
7359 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
7360 >= DBL_BUFSIZE_BOUND)
7363 size_t len = strlen (s);
7364 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
7365 if (width + x0 > write->c2)
7367 dfm_put_record_utf8 (writer, line->s.ss.string,
7369 u8_line_clear (line);
7372 u8_line_put (line, x0, x0 + width, s, len);
7375 x0 += write->format ? write->format->w : width + 1;
7378 if (y + 1 >= m->size1 && write->hold)
7380 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
7381 u8_line_clear (line);
7385 u8_line_destroy (line);
7389 write->wf->held = line;
7391 gsl_matrix_free (m);
7396 static struct matrix_command *
7397 matrix_get_parse (struct matrix_state *s)
7399 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7400 *cmd = (struct matrix_command) {
7404 .dataset = s->dataset,
7405 .user = { .treatment = MGET_ERROR },
7406 .system = { .treatment = MGET_ERROR },
7410 struct matrix_get *get = &cmd->get;
7411 get->dst = matrix_lvalue_parse (s);
7415 while (lex_match (s->lexer, T_SLASH))
7417 if (lex_match_id (s->lexer, "FILE"))
7419 lex_match (s->lexer, T_EQUALS);
7421 fh_unref (get->file);
7422 if (lex_match (s->lexer, T_ASTERISK))
7426 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7431 else if (lex_match_id (s->lexer, "ENCODING"))
7433 lex_match (s->lexer, T_EQUALS);
7434 if (!lex_force_string (s->lexer))
7437 free (get->encoding);
7438 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
7442 else if (lex_match_id (s->lexer, "VARIABLES"))
7444 lex_match (s->lexer, T_EQUALS);
7448 lex_sbc_only_once (s->lexer, "VARIABLES");
7452 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
7455 else if (lex_match_id (s->lexer, "NAMES"))
7457 lex_match (s->lexer, T_EQUALS);
7458 if (!lex_force_id (s->lexer))
7461 struct substring name = lex_tokss (s->lexer);
7462 get->names = matrix_var_lookup (s, name);
7464 get->names = matrix_var_create (s, name);
7467 else if (lex_match_id (s->lexer, "MISSING"))
7469 lex_match (s->lexer, T_EQUALS);
7470 if (lex_match_id (s->lexer, "ACCEPT"))
7471 get->user.treatment = MGET_ACCEPT;
7472 else if (lex_match_id (s->lexer, "OMIT"))
7473 get->user.treatment = MGET_OMIT;
7474 else if (lex_is_number (s->lexer))
7476 get->user.treatment = MGET_RECODE;
7477 get->user.substitute = lex_number (s->lexer);
7482 lex_error (s->lexer, _("Syntax error expecting ACCEPT or OMIT or "
7483 "a number for MISSING."));
7487 else if (lex_match_id (s->lexer, "SYSMIS"))
7489 lex_match (s->lexer, T_EQUALS);
7490 if (lex_match_id (s->lexer, "OMIT"))
7491 get->system.treatment = MGET_OMIT;
7492 else if (lex_is_number (s->lexer))
7494 get->system.treatment = MGET_RECODE;
7495 get->system.substitute = lex_number (s->lexer);
7500 lex_error (s->lexer, _("Syntax error expecting OMIT or a number "
7507 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
7508 "MISSING", "SYSMIS");
7513 if (get->user.treatment != MGET_ACCEPT)
7514 get->system.treatment = MGET_ERROR;
7519 matrix_command_destroy (cmd);
7524 matrix_get_execute__ (struct matrix_command *cmd, struct casereader *reader,
7525 const struct dictionary *dict)
7527 struct matrix_get *get = &cmd->get;
7528 struct variable **vars;
7533 if (!var_syntax_evaluate (get->lexer, get->vars, get->n_vars, dict,
7534 &vars, &n_vars, PV_NUMERIC))
7539 n_vars = dict_get_n_vars (dict);
7540 vars = xnmalloc (n_vars, sizeof *vars);
7541 for (size_t i = 0; i < n_vars; i++)
7543 struct variable *var = dict_get_var (dict, i);
7544 if (!var_is_numeric (var))
7546 msg_at (SE, cmd->location, _("Variable %s is not numeric."),
7547 var_get_name (var));
7557 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
7558 for (size_t i = 0; i < n_vars; i++)
7560 char s[sizeof (double)];
7562 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
7563 memcpy (&f, s, sizeof f);
7564 gsl_matrix_set (names, i, 0, f);
7567 gsl_matrix_free (get->names->value);
7568 get->names->value = names;
7572 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
7573 long long int casenum = 1;
7575 for (struct ccase *c = casereader_read (reader); c;
7576 c = casereader_read (reader), casenum++)
7578 if (n_rows >= m->size1)
7580 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
7581 for (size_t y = 0; y < n_rows; y++)
7582 for (size_t x = 0; x < n_vars; x++)
7583 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
7584 gsl_matrix_free (m);
7589 for (size_t x = 0; x < n_vars; x++)
7591 const struct variable *var = vars[x];
7592 double d = case_num (c, var);
7595 if (get->system.treatment == MGET_RECODE)
7596 d = get->system.substitute;
7597 else if (get->system.treatment == MGET_OMIT)
7601 msg_at (SE, cmd->location, _("Variable %s in case %lld "
7602 "is system-missing."),
7603 var_get_name (var), casenum);
7607 else if (var_is_num_missing (var, d) == MV_USER)
7609 if (get->user.treatment == MGET_RECODE)
7610 d = get->user.substitute;
7611 else if (get->user.treatment == MGET_OMIT)
7613 else if (get->user.treatment != MGET_ACCEPT)
7615 msg_at (SE, cmd->location,
7616 _("Variable %s in case %lld has user-missing "
7618 var_get_name (var), casenum, d);
7622 gsl_matrix_set (m, n_rows, x, d);
7633 matrix_lvalue_evaluate_and_assign (get->dst, m, cmd->location);
7636 gsl_matrix_free (m);
7641 matrix_open_casereader (const struct matrix_command *cmd,
7642 const char *command_name, struct file_handle *file,
7643 const char *encoding, struct dataset *dataset,
7644 struct casereader **readerp, struct dictionary **dictp)
7648 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
7649 return *readerp != NULL;
7653 if (dict_get_n_vars (dataset_dict (dataset)) == 0)
7655 msg_at (SE, cmd->location,
7656 _("The %s command cannot read an empty active file."),
7660 *readerp = proc_open (dataset);
7661 *dictp = dict_ref (dataset_dict (dataset));
7667 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
7668 struct casereader *reader, struct dictionary *dict)
7671 casereader_destroy (reader);
7673 proc_commit (dataset);
7677 matrix_get_execute (struct matrix_command *cmd)
7679 struct matrix_get *get = &cmd->get;
7680 struct casereader *r;
7681 struct dictionary *d;
7682 if (matrix_open_casereader (cmd, "GET", get->file, get->encoding,
7683 get->dataset, &r, &d))
7685 matrix_get_execute__ (cmd, r, d);
7686 matrix_close_casereader (get->file, get->dataset, r, d);
7693 variables_changed (const char *keyword,
7694 const struct string_array *new_vars,
7695 const struct msg_location *new_vars_location,
7696 const struct msg_location *new_location,
7697 const struct string_array *old_vars,
7698 const struct msg_location *old_vars_location,
7699 const struct msg_location *old_location)
7705 msg_at (SE, new_location,
7706 _("%s may only be specified on MSAVE if it was specified "
7707 "on the first MSAVE within MATRIX."), keyword);
7708 msg_at (SN, old_location,
7709 _("The first MSAVE in MATRIX did not specify %s."),
7711 msg_at (SN, new_vars_location,
7712 _("This is the specification of %s on a later MSAVE."),
7716 if (!string_array_equal_case (old_vars, new_vars))
7718 msg_at (SE, new_location,
7719 _("%s must specify the same variables on each MSAVE "
7720 "within a given MATRIX."), keyword);
7721 msg_at (SE, old_vars_location,
7722 _("This is the specification of %s on the first MSAVE."),
7724 msg_at (SE, new_vars_location,
7725 _("This is a different specification of %s on a later MSAVE."),
7734 msave_common_changed (const struct msave_common *old,
7735 const struct msave_common *new)
7737 if (new->outfile && !fh_equal (old->outfile, new->outfile))
7739 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7740 "within a single MATRIX command."));
7741 msg_at (SN, old->outfile_location,
7742 _("This is the OUTFILE on the first MSAVE command."));
7743 msg_at (SN, new->outfile_location,
7744 _("This is the OUTFILE on a later MSAVE command."));
7748 if (!variables_changed ("VARIABLES",
7749 &new->variables, new->variables_location, new->location,
7750 &old->variables, old->variables_location, old->location)
7751 && !variables_changed ("FNAMES",
7752 &new->fnames, new->fnames_location, new->location,
7753 &old->fnames, old->fnames_location, old->location)
7754 && !variables_changed ("SNAMES",
7755 &new->snames, new->snames_location, new->location,
7756 &old->snames, old->snames_location, old->location))
7763 msave_common_destroy (struct msave_common *common)
7767 msg_location_destroy (common->location);
7768 fh_unref (common->outfile);
7769 msg_location_destroy (common->outfile_location);
7770 string_array_destroy (&common->variables);
7771 msg_location_destroy (common->variables_location);
7772 string_array_destroy (&common->fnames);
7773 msg_location_destroy (common->fnames_location);
7774 string_array_destroy (&common->snames);
7775 msg_location_destroy (common->snames_location);
7777 for (size_t i = 0; i < common->n_factors; i++)
7778 matrix_expr_destroy (common->factors[i]);
7779 free (common->factors);
7781 for (size_t i = 0; i < common->n_splits; i++)
7782 matrix_expr_destroy (common->splits[i]);
7783 free (common->splits);
7785 dict_unref (common->dict);
7786 casewriter_destroy (common->writer);
7793 match_rowtype (struct lexer *lexer)
7795 static const char *rowtypes[] = {
7796 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7798 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7800 for (size_t i = 0; i < n_rowtypes; i++)
7801 if (lex_match_id (lexer, rowtypes[i]))
7804 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7809 parse_var_names (struct lexer *lexer, struct string_array *sa,
7810 struct msg_location **locationp)
7812 lex_match (lexer, T_EQUALS);
7814 string_array_clear (sa);
7815 msg_location_destroy (*locationp);
7818 struct dictionary *dict = dict_create (get_default_encoding ());
7821 int start_ofs = lex_ofs (lexer);
7822 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7823 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7824 int end_ofs = lex_ofs (lexer) - 1;
7829 for (size_t i = 0; i < n_names; i++)
7830 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7831 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7833 lex_ofs_error (lexer, start_ofs, end_ofs,
7834 _("Variable name %s is reserved."), names[i]);
7835 for (size_t j = 0; j < n_names; j++)
7841 sa->strings = names;
7842 sa->n = sa->allocated = n_names;
7843 *locationp = lex_ofs_location (lexer, start_ofs, end_ofs);
7848 static struct matrix_command *
7849 matrix_msave_parse (struct matrix_state *s)
7851 int start_ofs = lex_ofs (s->lexer);
7853 struct msave_common *common = xmalloc (sizeof *common);
7854 *common = (struct msave_common) { .outfile = NULL };
7856 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7857 *cmd = (struct matrix_command) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7859 struct matrix_expr *splits = NULL;
7860 struct matrix_expr *factors = NULL;
7862 struct matrix_msave *msave = &cmd->msave;
7863 msave->expr = matrix_parse_exp (s);
7867 while (lex_match (s->lexer, T_SLASH))
7869 if (lex_match_id (s->lexer, "TYPE"))
7871 lex_match (s->lexer, T_EQUALS);
7873 msave->rowtype = match_rowtype (s->lexer);
7874 if (!msave->rowtype)
7877 else if (lex_match_id (s->lexer, "OUTFILE"))
7879 lex_match (s->lexer, T_EQUALS);
7881 fh_unref (common->outfile);
7882 int start_ofs = lex_ofs (s->lexer);
7883 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7884 if (!common->outfile)
7886 msg_location_destroy (common->outfile_location);
7887 common->outfile_location = lex_ofs_location (s->lexer, start_ofs,
7888 lex_ofs (s->lexer) - 1);
7890 else if (lex_match_id (s->lexer, "VARIABLES"))
7892 if (!parse_var_names (s->lexer, &common->variables,
7893 &common->variables_location))
7896 else if (lex_match_id (s->lexer, "FNAMES"))
7898 if (!parse_var_names (s->lexer, &common->fnames,
7899 &common->fnames_location))
7902 else if (lex_match_id (s->lexer, "SNAMES"))
7904 if (!parse_var_names (s->lexer, &common->snames,
7905 &common->snames_location))
7908 else if (lex_match_id (s->lexer, "SPLIT"))
7910 lex_match (s->lexer, T_EQUALS);
7912 matrix_expr_destroy (splits);
7913 splits = matrix_parse_exp (s);
7917 else if (lex_match_id (s->lexer, "FACTOR"))
7919 lex_match (s->lexer, T_EQUALS);
7921 matrix_expr_destroy (factors);
7922 factors = matrix_parse_exp (s);
7928 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7929 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7933 if (!msave->rowtype)
7935 lex_sbc_missing (s->lexer, "TYPE");
7939 if (!s->msave_common)
7941 if (common->fnames.n && !factors)
7943 msg_at (SE, common->fnames_location, _("FNAMES requires FACTOR."));
7946 if (common->snames.n && !splits)
7948 msg_at (SE, common->snames_location, _("SNAMES requires SPLIT."));
7951 if (!common->outfile)
7953 lex_sbc_missing (s->lexer, "OUTFILE");
7956 common->location = lex_ofs_location (s->lexer, start_ofs,
7957 lex_ofs (s->lexer));
7958 msg_location_remove_columns (common->location);
7959 s->msave_common = common;
7963 if (msave_common_changed (s->msave_common, common))
7965 msave_common_destroy (common);
7967 msave->common = s->msave_common;
7969 struct msave_common *c = s->msave_common;
7972 if (c->n_factors >= c->allocated_factors)
7973 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7974 sizeof *c->factors);
7975 c->factors[c->n_factors++] = factors;
7977 if (c->n_factors > 0)
7978 msave->factors = c->factors[c->n_factors - 1];
7982 if (c->n_splits >= c->allocated_splits)
7983 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7985 c->splits[c->n_splits++] = splits;
7987 if (c->n_splits > 0)
7988 msave->splits = c->splits[c->n_splits - 1];
7993 matrix_expr_destroy (splits);
7994 matrix_expr_destroy (factors);
7995 msave_common_destroy (common);
7996 matrix_command_destroy (cmd);
8001 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
8003 gsl_matrix *m = matrix_expr_evaluate (e);
8009 msg_at (SE, matrix_expr_location (e),
8010 _("%s expression must evaluate to vector, "
8011 "not a %zu×%zu matrix."),
8012 name, m->size1, m->size2);
8013 gsl_matrix_free (m);
8017 return matrix_to_vector (m);
8021 msave_add_vars (struct dictionary *d, const struct string_array *vars)
8023 for (size_t i = 0; i < vars->n; i++)
8024 if (!dict_create_var (d, vars->strings[i], 0))
8025 return vars->strings[i];
8029 static struct dictionary *
8030 msave_create_dict (const struct msave_common *common)
8032 struct dictionary *dict = dict_create (get_default_encoding ());
8034 const char *dup_split = msave_add_vars (dict, &common->snames);
8037 /* Should not be possible because the parser ensures that the names are
8042 dict_create_var_assert (dict, "ROWTYPE_", 8);
8044 const char *dup_factor = msave_add_vars (dict, &common->fnames);
8047 msg_at (SE, common->fnames_location,
8048 _("Duplicate or invalid FACTOR variable name %s."),
8053 dict_create_var_assert (dict, "VARNAME_", 8);
8055 const char *dup_var = msave_add_vars (dict, &common->variables);
8058 msg_at (SE, common->variables_location,
8059 _("Duplicate or invalid variable name %s."),
8072 matrix_msave_execute (struct matrix_command *cmd)
8074 struct matrix_msave *msave = &cmd->msave;
8075 struct msave_common *common = msave->common;
8076 gsl_matrix *m = NULL;
8077 gsl_vector *factors = NULL;
8078 gsl_vector *splits = NULL;
8080 m = matrix_expr_evaluate (msave->expr);
8084 if (!common->variables.n)
8085 for (size_t i = 0; i < m->size2; i++)
8086 string_array_append_nocopy (&common->variables,
8087 xasprintf ("COL%zu", i + 1));
8088 else if (m->size2 != common->variables.n)
8090 msg_at (SE, matrix_expr_location (msave->expr),
8091 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
8092 m->size2, common->variables.n);
8098 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
8102 if (!common->fnames.n)
8103 for (size_t i = 0; i < factors->size; i++)
8104 string_array_append_nocopy (&common->fnames,
8105 xasprintf ("FAC%zu", i + 1));
8106 else if (factors->size != common->fnames.n)
8108 msg_at (SE, matrix_expr_location (msave->factors),
8109 _("There are %zu factor variables, "
8110 "but %zu factor values were supplied."),
8111 common->fnames.n, factors->size);
8118 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
8122 if (!common->snames.n)
8123 for (size_t i = 0; i < splits->size; i++)
8124 string_array_append_nocopy (&common->snames,
8125 xasprintf ("SPL%zu", i + 1));
8126 else if (splits->size != common->snames.n)
8128 msg_at (SE, matrix_expr_location (msave->splits),
8129 _("There are %zu split variables, "
8130 "but %zu split values were supplied."),
8131 common->snames.n, splits->size);
8136 if (!common->writer)
8138 struct dictionary *dict = msave_create_dict (common);
8142 common->writer = any_writer_open (common->outfile, dict);
8143 if (!common->writer)
8149 common->dict = dict;
8152 bool matrix = (!strcmp (msave->rowtype, "COV")
8153 || !strcmp (msave->rowtype, "CORR"));
8154 for (size_t y = 0; y < m->size1; y++)
8156 struct ccase *c = case_create (dict_get_proto (common->dict));
8159 /* Split variables */
8161 for (size_t i = 0; i < splits->size; i++)
8162 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
8165 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8166 msave->rowtype, ' ');
8170 for (size_t i = 0; i < factors->size; i++)
8171 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
8174 const char *varname_ = (matrix && y < common->variables.n
8175 ? common->variables.strings[y]
8177 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8180 /* Continuous variables. */
8181 for (size_t x = 0; x < m->size2; x++)
8182 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
8183 casewriter_write (common->writer, c);
8187 gsl_matrix_free (m);
8188 gsl_vector_free (factors);
8189 gsl_vector_free (splits);
8194 static struct matrix_command *
8195 matrix_mget_parse (struct matrix_state *s)
8197 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8198 *cmd = (struct matrix_command) {
8202 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
8206 struct matrix_mget *mget = &cmd->mget;
8208 lex_match (s->lexer, T_SLASH);
8209 while (lex_token (s->lexer) != T_ENDCMD)
8211 if (lex_match_id (s->lexer, "FILE"))
8213 lex_match (s->lexer, T_EQUALS);
8215 fh_unref (mget->file);
8216 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
8220 else if (lex_match_id (s->lexer, "ENCODING"))
8222 lex_match (s->lexer, T_EQUALS);
8223 if (!lex_force_string (s->lexer))
8226 free (mget->encoding);
8227 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
8231 else if (lex_match_id (s->lexer, "TYPE"))
8233 lex_match (s->lexer, T_EQUALS);
8234 while (lex_token (s->lexer) != T_SLASH
8235 && lex_token (s->lexer) != T_ENDCMD)
8237 const char *rowtype = match_rowtype (s->lexer);
8241 stringi_set_insert (&mget->rowtypes, rowtype);
8246 lex_error_expecting (s->lexer, "FILE", "TYPE");
8249 lex_match (s->lexer, T_SLASH);
8254 matrix_command_destroy (cmd);
8258 static const struct variable *
8259 get_a8_var (const struct msg_location *loc,
8260 const struct dictionary *d, const char *name)
8262 const struct variable *v = dict_lookup_var (d, name);
8265 msg_at (SE, loc, _("Matrix data file lacks %s variable."), name);
8268 if (var_get_width (v) != 8)
8270 msg_at (SE, loc, _("%s variable in matrix data file must be "
8271 "8-byte string, but it has width %d."),
8272 name, var_get_width (v));
8279 var_changed (const struct ccase *ca, const struct ccase *cb,
8280 const struct variable *var)
8283 ? !value_equal (case_data (ca, var), case_data (cb, var),
8284 var_get_width (var))
8289 vars_changed (const struct ccase *ca, const struct ccase *cb,
8290 const struct dictionary *d,
8291 size_t first_var, size_t n_vars)
8293 for (size_t i = 0; i < n_vars; i++)
8295 const struct variable *v = dict_get_var (d, first_var + i);
8296 if (var_changed (ca, cb, v))
8303 vars_all_missing (const struct ccase *c, const struct dictionary *d,
8304 size_t first_var, size_t n_vars)
8306 for (size_t i = 0; i < n_vars; i++)
8307 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
8313 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
8314 const struct dictionary *d,
8315 const struct variable *rowtype_var,
8316 const struct stringi_set *accepted_rowtypes,
8317 struct matrix_state *s,
8318 size_t ss, size_t sn, size_t si,
8319 size_t fs, size_t fn, size_t fi,
8320 size_t cs, size_t cn,
8321 struct pivot_table *pt,
8322 struct pivot_dimension *var_dimension)
8327 /* Is this a matrix for pooled data, either where there are no factor
8328 variables or the factor variables are missing? */
8329 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
8331 struct substring rowtype = case_ss (rows[0], rowtype_var);
8332 ss_rtrim (&rowtype, ss_cstr (" "));
8333 if (!stringi_set_is_empty (accepted_rowtypes)
8334 && !stringi_set_contains_len (accepted_rowtypes,
8335 rowtype.string, rowtype.length))
8338 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
8339 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
8340 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
8341 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
8342 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
8343 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
8347 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
8348 (int) rowtype.length, rowtype.string);
8352 struct string name = DS_EMPTY_INITIALIZER;
8353 ds_put_cstr (&name, prefix);
8355 ds_put_format (&name, "F%zu", fi);
8357 ds_put_format (&name, "S%zu", si);
8359 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
8361 mv = matrix_var_create (s, ds_ss (&name));
8364 msg (SW, _("Matrix data file contains variable with existing name %s."),
8366 goto exit_free_name;
8369 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
8370 size_t n_missing = 0;
8371 for (size_t y = 0; y < n_rows; y++)
8373 for (size_t x = 0; x < cn; x++)
8375 struct variable *var = dict_get_var (d, cs + x);
8376 double value = case_num (rows[y], var);
8377 if (var_is_num_missing (var, value))
8382 gsl_matrix_set (m, y, x, value);
8386 int var_index = pivot_category_create_leaf (
8387 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
8388 double values[] = { n_rows, cn };
8389 for (size_t j = 0; j < sn; j++)
8391 struct variable *var = dict_get_var (d, ss + j);
8392 const union value *value = case_data (rows[0], var);
8393 pivot_table_put2 (pt, j, var_index,
8394 pivot_value_new_var_value (var, value));
8396 for (size_t j = 0; j < fn; j++)
8398 struct variable *var = dict_get_var (d, fs + j);
8399 const union value sysmis = { .f = SYSMIS };
8400 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
8401 pivot_table_put2 (pt, j + sn, var_index,
8402 pivot_value_new_var_value (var, value));
8404 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
8405 pivot_table_put2 (pt, j + sn + fn, var_index,
8406 pivot_value_new_integer (values[j]));
8409 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
8410 "value, which was treated as zero.",
8411 "Matrix data file variable %s contains %zu missing "
8412 "values, which were treated as zero.", n_missing),
8413 ds_cstr (&name), n_missing);
8420 for (size_t y = 0; y < n_rows; y++)
8421 case_unref (rows[y]);
8425 matrix_mget_execute__ (struct matrix_command *cmd, struct casereader *r,
8426 const struct dictionary *d)
8428 struct matrix_mget *mget = &cmd->mget;
8429 const struct msg_location *loc = cmd->location;
8430 const struct variable *rowtype_ = get_a8_var (loc, d, "ROWTYPE_");
8431 const struct variable *varname_ = get_a8_var (loc, d, "VARNAME_");
8432 if (!rowtype_ || !varname_)
8435 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
8438 _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
8441 if (var_get_dict_index (varname_) + 1 >= dict_get_n_vars (d))
8443 msg_at (SE, loc, _("Matrix data file contains no continuous variables."));
8447 for (size_t i = 0; i < dict_get_n_vars (d); i++)
8449 const struct variable *v = dict_get_var (d, i);
8450 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
8453 _("Matrix data file contains unexpected string variable %s."),
8459 /* SPLIT variables. */
8461 size_t sn = var_get_dict_index (rowtype_);
8462 struct ccase *sc = NULL;
8465 /* FACTOR variables. */
8466 size_t fs = var_get_dict_index (rowtype_) + 1;
8467 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
8468 struct ccase *fc = NULL;
8471 /* Continuous variables. */
8472 size_t cs = var_get_dict_index (varname_) + 1;
8473 size_t cn = dict_get_n_vars (d) - cs;
8474 struct ccase *cc = NULL;
8477 struct pivot_table *pt = pivot_table_create (
8478 N_("Matrix Variables Created by MGET"));
8479 struct pivot_dimension *attr_dimension = pivot_dimension_create (
8480 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
8481 struct pivot_dimension *var_dimension = pivot_dimension_create (
8482 pt, PIVOT_AXIS_ROW, N_("Variable"));
8485 struct pivot_category *splits = pivot_category_create_group (
8486 attr_dimension->root, N_("Split Values"));
8487 for (size_t i = 0; i < sn; i++)
8488 pivot_category_create_leaf (splits, pivot_value_new_variable (
8489 dict_get_var (d, ss + i)));
8493 struct pivot_category *factors = pivot_category_create_group (
8494 attr_dimension->root, N_("Factors"));
8495 for (size_t i = 0; i < fn; i++)
8496 pivot_category_create_leaf (factors, pivot_value_new_variable (
8497 dict_get_var (d, fs + i)));
8499 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
8500 N_("Rows"), N_("Columns"));
8503 struct ccase **rows = NULL;
8504 size_t allocated_rows = 0;
8508 while ((c = casereader_read (r)) != NULL)
8510 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
8520 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
8521 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
8522 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
8525 if (change != NOTHING_CHANGED)
8527 matrix_mget_commit_var (rows, n_rows, d,
8528 rowtype_, &mget->rowtypes,
8539 if (n_rows >= allocated_rows)
8540 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
8543 if (change == SPLITS_CHANGED)
8549 /* Reset the factor number, if there are factors. */
8553 if (row_has_factors)
8559 else if (change == FACTORS_CHANGED)
8561 if (row_has_factors)
8567 matrix_mget_commit_var (rows, n_rows, d,
8568 rowtype_, &mget->rowtypes,
8580 if (var_dimension->n_leaves)
8581 pivot_table_submit (pt);
8583 pivot_table_unref (pt);
8587 matrix_mget_execute (struct matrix_command *cmd)
8589 struct matrix_mget *mget = &cmd->mget;
8590 struct casereader *r;
8591 struct dictionary *d;
8592 if (matrix_open_casereader (cmd, "MGET", mget->file, mget->encoding,
8593 mget->state->dataset, &r, &d))
8595 matrix_mget_execute__ (cmd, r, d);
8596 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
8603 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
8605 if (!lex_force_id (s->lexer))
8608 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
8610 *varp = matrix_var_create (s, lex_tokss (s->lexer));
8615 static struct matrix_command *
8616 matrix_eigen_parse (struct matrix_state *s)
8618 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8619 *cmd = (struct matrix_command) {
8621 .eigen = { .expr = NULL }
8624 struct matrix_eigen *eigen = &cmd->eigen;
8625 if (!lex_force_match (s->lexer, T_LPAREN))
8627 eigen->expr = matrix_expr_parse (s);
8629 || !lex_force_match (s->lexer, T_COMMA)
8630 || !matrix_parse_dst_var (s, &eigen->evec)
8631 || !lex_force_match (s->lexer, T_COMMA)
8632 || !matrix_parse_dst_var (s, &eigen->eval)
8633 || !lex_force_match (s->lexer, T_RPAREN))
8639 matrix_command_destroy (cmd);
8644 matrix_eigen_execute (struct matrix_command *cmd)
8646 struct matrix_eigen *eigen = &cmd->eigen;
8647 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
8650 if (!is_symmetric (A))
8652 msg_at (SE, cmd->location, _("Argument of EIGEN must be symmetric."));
8653 gsl_matrix_free (A);
8657 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
8658 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
8659 gsl_vector v_eval = to_vector (eval);
8660 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
8661 gsl_eigen_symmv (A, &v_eval, evec, w);
8662 gsl_eigen_symmv_free (w);
8664 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
8666 gsl_matrix_free (A);
8668 gsl_matrix_free (eigen->eval->value);
8669 eigen->eval->value = eval;
8671 gsl_matrix_free (eigen->evec->value);
8672 eigen->evec->value = evec;
8677 static struct matrix_command *
8678 matrix_setdiag_parse (struct matrix_state *s)
8680 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8681 *cmd = (struct matrix_command) {
8682 .type = MCMD_SETDIAG,
8683 .setdiag = { .dst = NULL }
8686 struct matrix_setdiag *setdiag = &cmd->setdiag;
8687 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
8690 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
8693 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
8698 if (!lex_force_match (s->lexer, T_COMMA))
8701 setdiag->expr = matrix_expr_parse (s);
8705 if (!lex_force_match (s->lexer, T_RPAREN))
8711 matrix_command_destroy (cmd);
8716 matrix_setdiag_execute (struct matrix_command *cmd)
8718 struct matrix_setdiag *setdiag = &cmd->setdiag;
8719 gsl_matrix *dst = setdiag->dst->value;
8722 msg_at (SE, cmd->location,
8723 _("SETDIAG destination matrix %s is uninitialized."),
8724 setdiag->dst->name);
8728 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
8732 size_t n = MIN (dst->size1, dst->size2);
8733 if (is_scalar (src))
8735 double d = to_scalar (src);
8736 for (size_t i = 0; i < n; i++)
8737 gsl_matrix_set (dst, i, i, d);
8739 else if (is_vector (src))
8741 gsl_vector v = to_vector (src);
8742 for (size_t i = 0; i < n && i < v.size; i++)
8743 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
8746 msg_at (SE, matrix_expr_location (setdiag->expr),
8747 _("SETDIAG argument 2 must be a scalar or a vector, "
8748 "not a %zu×%zu matrix."),
8749 src->size1, src->size2);
8750 gsl_matrix_free (src);
8755 static struct matrix_command *
8756 matrix_svd_parse (struct matrix_state *s)
8758 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8759 *cmd = (struct matrix_command) {
8761 .svd = { .expr = NULL }
8764 struct matrix_svd *svd = &cmd->svd;
8765 if (!lex_force_match (s->lexer, T_LPAREN))
8767 svd->expr = matrix_expr_parse (s);
8769 || !lex_force_match (s->lexer, T_COMMA)
8770 || !matrix_parse_dst_var (s, &svd->u)
8771 || !lex_force_match (s->lexer, T_COMMA)
8772 || !matrix_parse_dst_var (s, &svd->s)
8773 || !lex_force_match (s->lexer, T_COMMA)
8774 || !matrix_parse_dst_var (s, &svd->v)
8775 || !lex_force_match (s->lexer, T_RPAREN))
8781 matrix_command_destroy (cmd);
8786 matrix_svd_execute (struct matrix_svd *svd)
8788 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8792 if (m->size1 >= m->size2)
8795 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8796 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8797 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8798 gsl_vector *work = gsl_vector_alloc (A->size2);
8799 gsl_linalg_SV_decomp (A, V, &Sv, work);
8800 gsl_vector_free (work);
8802 matrix_var_set (svd->u, A);
8803 matrix_var_set (svd->s, S);
8804 matrix_var_set (svd->v, V);
8808 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8809 gsl_matrix_transpose_memcpy (At, m);
8810 gsl_matrix_free (m);
8812 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8813 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8814 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8815 gsl_vector *work = gsl_vector_alloc (At->size2);
8816 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8817 gsl_vector_free (work);
8819 matrix_var_set (svd->v, At);
8820 matrix_var_set (svd->s, St);
8821 matrix_var_set (svd->u, Vt);
8825 /* The main MATRIX command logic. */
8828 matrix_command_execute (struct matrix_command *cmd)
8833 matrix_compute_execute (cmd);
8837 matrix_print_execute (&cmd->print);
8841 return matrix_do_if_execute (&cmd->do_if);
8844 matrix_loop_execute (&cmd->loop);
8851 matrix_display_execute (&cmd->display);
8855 matrix_release_execute (&cmd->release);
8859 matrix_save_execute (cmd);
8863 matrix_read_execute (cmd);
8867 matrix_write_execute (&cmd->write);
8871 matrix_get_execute (cmd);
8875 matrix_msave_execute (cmd);
8879 matrix_mget_execute (cmd);
8883 matrix_eigen_execute (cmd);
8887 matrix_setdiag_execute (cmd);
8891 matrix_svd_execute (&cmd->svd);
8899 matrix_command_destroy (struct matrix_command *cmd)
8904 msg_location_destroy (cmd->location);
8909 matrix_lvalue_destroy (cmd->compute.lvalue);
8910 matrix_expr_destroy (cmd->compute.rvalue);
8914 matrix_expr_destroy (cmd->print.expression);
8915 free (cmd->print.title);
8916 print_labels_destroy (cmd->print.rlabels);
8917 print_labels_destroy (cmd->print.clabels);
8921 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8923 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8924 matrix_commands_uninit (&cmd->do_if.clauses[i].commands);
8926 free (cmd->do_if.clauses);
8930 matrix_expr_destroy (cmd->loop.start);
8931 matrix_expr_destroy (cmd->loop.end);
8932 matrix_expr_destroy (cmd->loop.increment);
8933 matrix_expr_destroy (cmd->loop.top_condition);
8934 matrix_expr_destroy (cmd->loop.bottom_condition);
8935 matrix_commands_uninit (&cmd->loop.commands);
8945 free (cmd->release.vars);
8949 matrix_expr_destroy (cmd->save.expression);
8953 matrix_lvalue_destroy (cmd->read.dst);
8954 matrix_expr_destroy (cmd->read.size);
8958 matrix_expr_destroy (cmd->write.expression);
8959 free (cmd->write.format);
8963 matrix_lvalue_destroy (cmd->get.dst);
8964 fh_unref (cmd->get.file);
8965 free (cmd->get.encoding);
8966 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8970 matrix_expr_destroy (cmd->msave.expr);
8974 fh_unref (cmd->mget.file);
8975 stringi_set_destroy (&cmd->mget.rowtypes);
8979 matrix_expr_destroy (cmd->eigen.expr);
8983 matrix_expr_destroy (cmd->setdiag.expr);
8987 matrix_expr_destroy (cmd->svd.expr);
8994 matrix_commands_parse (struct matrix_state *s, struct matrix_commands *c,
8995 const char *command_name,
8996 const char *stop1, const char *stop2)
8998 lex_end_of_command (s->lexer);
8999 lex_discard_rest_of_command (s->lexer);
9001 size_t allocated = 0;
9004 while (lex_token (s->lexer) == T_ENDCMD)
9007 if (lex_at_phrase (s->lexer, stop1)
9008 || (stop2 && lex_at_phrase (s->lexer, stop2)))
9011 if (lex_at_phrase (s->lexer, "END MATRIX"))
9013 lex_next_error (s->lexer, 0, 1,
9014 _("Premature END MATRIX within %s."), command_name);
9018 struct matrix_command *cmd = matrix_command_parse (s);
9022 if (c->n >= allocated)
9023 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
9024 c->commands[c->n++] = cmd;
9029 matrix_commands_uninit (struct matrix_commands *cmds)
9031 for (size_t i = 0; i < cmds->n; i++)
9032 matrix_command_destroy (cmds->commands[i]);
9033 free (cmds->commands);
9036 struct matrix_command_name
9039 struct matrix_command *(*parse) (struct matrix_state *);
9042 static const struct matrix_command_name *
9043 matrix_command_name_parse (struct lexer *lexer)
9045 static const struct matrix_command_name commands[] = {
9046 { "COMPUTE", matrix_compute_parse },
9047 { "CALL EIGEN", matrix_eigen_parse },
9048 { "CALL SETDIAG", matrix_setdiag_parse },
9049 { "CALL SVD", matrix_svd_parse },
9050 { "PRINT", matrix_print_parse },
9051 { "DO IF", matrix_do_if_parse },
9052 { "LOOP", matrix_loop_parse },
9053 { "BREAK", matrix_break_parse },
9054 { "READ", matrix_read_parse },
9055 { "WRITE", matrix_write_parse },
9056 { "GET", matrix_get_parse },
9057 { "SAVE", matrix_save_parse },
9058 { "MGET", matrix_mget_parse },
9059 { "MSAVE", matrix_msave_parse },
9060 { "DISPLAY", matrix_display_parse },
9061 { "RELEASE", matrix_release_parse },
9063 static size_t n = sizeof commands / sizeof *commands;
9065 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
9066 if (lex_match_phrase (lexer, c->name))
9071 static struct matrix_command *
9072 matrix_command_parse (struct matrix_state *s)
9074 int start_ofs = lex_ofs (s->lexer);
9075 size_t nesting_level = SIZE_MAX;
9077 struct matrix_command *c = NULL;
9078 const struct matrix_command_name *cmd = matrix_command_name_parse (s->lexer);
9080 lex_error (s->lexer, _("Unknown matrix command."));
9081 else if (!cmd->parse)
9082 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
9086 nesting_level = output_open_group (
9087 group_item_create_nocopy (utf8_to_title (cmd->name),
9088 utf8_to_title (cmd->name)));
9094 c->location = lex_ofs_location (s->lexer, start_ofs, lex_ofs (s->lexer));
9095 msg_location_remove_columns (c->location);
9096 lex_end_of_command (s->lexer);
9098 lex_discard_rest_of_command (s->lexer);
9099 if (nesting_level != SIZE_MAX)
9100 output_close_groups (nesting_level);
9106 cmd_matrix (struct lexer *lexer, struct dataset *ds)
9108 if (!lex_force_match (lexer, T_ENDCMD))
9111 struct matrix_state state = {
9113 .session = dataset_session (ds),
9115 .vars = HMAP_INITIALIZER (state.vars),
9120 while (lex_match (lexer, T_ENDCMD))
9122 if (lex_token (lexer) == T_STOP)
9124 msg (SE, _("Unexpected end of input expecting matrix command."));
9128 if (lex_match_phrase (lexer, "END MATRIX"))
9131 struct matrix_command *c = matrix_command_parse (&state);
9134 matrix_command_execute (c);
9135 matrix_command_destroy (c);
9139 struct matrix_var *var, *next;
9140 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
9143 gsl_matrix_free (var->value);
9144 hmap_delete (&state.vars, &var->hmap_node);
9147 hmap_destroy (&state.vars);
9148 msave_common_destroy (state.msave_common);
9149 fh_unref (state.prev_read_file);
9150 for (size_t i = 0; i < state.n_read_files; i++)
9151 read_file_destroy (state.read_files[i]);
9152 free (state.read_files);
9153 fh_unref (state.prev_write_file);
9154 for (size_t i = 0; i < state.n_write_files; i++)
9155 write_file_destroy (state.write_files[i]);
9156 free (state.write_files);
9157 fh_unref (state.prev_save_file);
9158 for (size_t i = 0; i < state.n_save_files; i++)
9159 save_file_destroy (state.save_files[i]);
9160 free (state.save_files);