1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2021 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
19 #include <gsl/gsl_blas.h>
20 #include <gsl/gsl_cdf.h>
21 #include <gsl/gsl_eigen.h>
22 #include <gsl/gsl_linalg.h>
23 #include <gsl/gsl_matrix.h>
24 #include <gsl/gsl_permutation.h>
25 #include <gsl/gsl_randist.h>
26 #include <gsl/gsl_vector.h>
31 #include "data/any-reader.h"
32 #include "data/any-writer.h"
33 #include "data/casereader.h"
34 #include "data/casewriter.h"
35 #include "data/data-in.h"
36 #include "data/data-out.h"
37 #include "data/dataset.h"
38 #include "data/dictionary.h"
39 #include "data/file-handle-def.h"
40 #include "language/command.h"
41 #include "language/data-io/data-reader.h"
42 #include "language/data-io/data-writer.h"
43 #include "language/data-io/file-handle.h"
44 #include "language/lexer/format-parser.h"
45 #include "language/lexer/lexer.h"
46 #include "language/lexer/variable-parser.h"
47 #include "libpspp/array.h"
48 #include "libpspp/assertion.h"
49 #include "libpspp/compiler.h"
50 #include "libpspp/hmap.h"
51 #include "libpspp/i18n.h"
52 #include "libpspp/intern.h"
53 #include "libpspp/misc.h"
54 #include "libpspp/str.h"
55 #include "libpspp/string-array.h"
56 #include "libpspp/stringi-set.h"
57 #include "libpspp/u8-line.h"
58 #include "math/distributions.h"
59 #include "math/random.h"
60 #include "output/driver.h"
61 #include "output/output-item.h"
62 #include "output/pivot-table.h"
64 #include "gl/c-ctype.h"
65 #include "gl/c-strcase.h"
66 #include "gl/ftoastr.h"
67 #include "gl/minmax.h"
71 #define _(msgid) gettext (msgid)
72 #define N_(msgid) (msgid)
76 /* 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 string_array variables; /* VARIABLES subcommand. */
93 struct string_array fnames; /* FNAMES subcommand. */
94 struct string_array snames; /* SNAMES subcommand. */
96 /* Collects and owns factors and splits. The individual msave_command
97 structs point to these but do not own them. (This is because factors
98 and splits can be carried over from one MSAVE to the next, so it's
99 easiest to just take the most recent.) */
100 struct matrix_expr **factors;
101 size_t n_factors, allocated_factors;
102 struct matrix_expr **splits;
103 size_t n_splits, allocated_splits;
105 /* Execution state. */
106 struct dictionary *dict;
107 struct casewriter *writer;
110 /* A file used by one or more READ commands. */
114 struct file_handle *file;
116 /* Execution state. */
117 struct dfm_reader *reader;
121 static struct read_file *read_file_create (struct matrix_state *,
122 struct file_handle *);
123 static struct dfm_reader *read_file_open (struct read_file *);
125 /* A file used by one or more WRITE comamnds. */
129 struct file_handle *file;
131 /* Execution state. */
132 struct dfm_writer *writer;
134 struct u8_line *held; /* Output held by a previous WRITE with /HOLD. */
137 static struct write_file *write_file_create (struct matrix_state *,
138 struct file_handle *);
139 static struct dfm_writer *write_file_open (struct write_file *);
140 static void write_file_destroy (struct write_file *);
142 /* A file used by one or more SAVE commands. */
146 struct file_handle *file;
147 struct dataset *dataset;
148 struct string_array variables;
149 struct matrix_expr *names;
150 struct stringi_set strings;
152 /* Execution state. */
154 struct casewriter *writer;
155 struct dictionary *dict;
156 struct msg_location *location;
159 /* State of an entire matrix program. */
162 /* State passed into MATRIX from outside. */
163 struct dataset *dataset;
164 struct session *session;
167 /* Matrix program's own state. */
168 struct hmap vars; /* Dictionary of matrix variables. */
169 bool in_loop; /* True if parsing within a LOOP. */
172 struct msave_common *msave_common;
175 struct file_handle *prev_read_file;
176 struct read_file **read_files;
180 struct file_handle *prev_write_file;
181 struct write_file **write_files;
182 size_t n_write_files;
185 struct file_handle *prev_save_file;
186 struct save_file **save_files;
190 /* Finds and returns the variable with the given NAME (case-insensitive) within
191 S, if there is one, or a null pointer if there is not. */
192 static struct matrix_var *
193 matrix_var_lookup (struct matrix_state *s, struct substring name)
195 struct matrix_var *var;
197 HMAP_FOR_EACH_WITH_HASH (var, struct matrix_var, hmap_node,
198 utf8_hash_case_substring (name, 0), &s->vars)
199 if (!utf8_sscasecmp (ss_cstr (var->name), name))
205 /* Creates and returns a new variable named NAME within S. There must not
206 already be a variable with the same (case-insensitive) name. The variable
207 is created uninitialized. */
208 static struct matrix_var *
209 matrix_var_create (struct matrix_state *s, struct substring name)
211 struct matrix_var *var = xmalloc (sizeof *var);
212 *var = (struct matrix_var) { .name = ss_xstrdup (name) };
213 hmap_insert (&s->vars, &var->hmap_node, utf8_hash_case_substring (name, 0));
217 /* Replaces VAR's value by VALUE. Takes ownership of VALUE. */
219 matrix_var_set (struct matrix_var *var, gsl_matrix *value)
221 gsl_matrix_free (var->value);
225 /* Matrix function catalog. */
227 /* The third argument to F() is a "prototype". For most prototypes, the first
228 letter (before the _) represents the return type and each other letter
229 (after the _) is an argument type. The types are:
231 - "m": A matrix of unrestricted dimensions.
235 - "v": A row or column vector.
237 - "e": Primarily for the first argument, this is a matrix with
238 unrestricted dimensions treated elementwise. Each element in the matrix
239 is passed to the implementation function separately.
241 - "n": This gets passed the "const struct matrix_expr *" that represents
242 the expression. This allows the evaluation function to grab the source
243 location of arguments so that it can report accurate error locations.
244 This type doesn't correspond to an argument passed in by the user.
246 The fourth argument is an optional constraints string. For this purpose the
247 first argument is named "a", the second "b", and so on. The following kinds
248 of constraints are supported. For matrix arguments, the constraints are
249 applied to each value in the matrix separately:
251 - "a(0,1)" or "a[0,1]": 0 < a < 1 or 0 <= a <= 1, respectively. Any
252 integer may substitute for 0 and 1. Half-open constraints (] and [) are
255 - "ai": Restrict a to integer values.
257 - "a>0", "a<0", "a>=0", "a<=0", "a!=0".
259 - "a<b", "a>b", "a<=b", "a>=b", "b!=0".
261 #define MATRIX_FUNCTIONS \
262 F(ABS, "ABS", m_e, NULL) \
263 F(ALL, "ALL", d_m, NULL) \
264 F(ANY, "ANY", d_m, NULL) \
265 F(ARSIN, "ARSIN", m_e, "a[-1,1]") \
266 F(ARTAN, "ARTAN", m_e, NULL) \
267 F(BLOCK, "BLOCK", m_any, NULL) \
268 F(CHOL, "CHOL", m_mn, NULL) \
269 F(CMIN, "CMIN", m_m, NULL) \
270 F(CMAX, "CMAX", m_m, NULL) \
271 F(COS, "COS", m_e, NULL) \
272 F(CSSQ, "CSSQ", m_m, NULL) \
273 F(CSUM, "CSUM", m_m, NULL) \
274 F(DESIGN, "DESIGN", m_mn, NULL) \
275 F(DET, "DET", d_m, NULL) \
276 F(DIAG, "DIAG", m_m, NULL) \
277 F(EVAL, "EVAL", m_mn, NULL) \
278 F(EXP, "EXP", m_e, NULL) \
279 F(GINV, "GINV", m_m, NULL) \
280 F(GRADE, "GRADE", m_m, NULL) \
281 F(GSCH, "GSCH", m_mn, NULL) \
282 F(IDENT, "IDENT", IDENT, NULL) \
283 F(INV, "INV", m_m, NULL) \
284 F(KRONEKER, "KRONEKER", m_mm, NULL) \
285 F(LG10, "LG10", m_e, "a>0") \
286 F(LN, "LN", m_e, "a>0") \
287 F(MAGIC, "MAGIC", m_d, "ai>=3") \
288 F(MAKE, "MAKE", m_ddd, "ai>=0 bi>=0") \
289 F(MDIAG, "MDIAG", m_v, NULL) \
290 F(MMAX, "MMAX", d_m, NULL) \
291 F(MMIN, "MMIN", d_m, NULL) \
292 F(MOD, "MOD", m_md, "b!=0") \
293 F(MSSQ, "MSSQ", d_m, NULL) \
294 F(MSUM, "MSUM", d_m, NULL) \
295 F(NCOL, "NCOL", d_m, NULL) \
296 F(NROW, "NROW", d_m, NULL) \
297 F(RANK, "RANK", d_m, NULL) \
298 F(RESHAPE, "RESHAPE", m_mddn, NULL) \
299 F(RMAX, "RMAX", m_m, NULL) \
300 F(RMIN, "RMIN", m_m, NULL) \
301 F(RND, "RND", m_e, NULL) \
302 F(RNKORDER, "RNKORDER", m_m, NULL) \
303 F(RSSQ, "RSSQ", m_m, NULL) \
304 F(RSUM, "RSUM", m_m, NULL) \
305 F(SIN, "SIN", m_e, NULL) \
306 F(SOLVE, "SOLVE", m_mmn, NULL) \
307 F(SQRT, "SQRT", m_e, "a>=0") \
308 F(SSCP, "SSCP", m_m, NULL) \
309 F(SVAL, "SVAL", m_m, NULL) \
310 F(SWEEP, "SWEEP", m_mdn, NULL) \
311 F(T, "T", m_m, NULL) \
312 F(TRACE, "TRACE", d_m, NULL) \
313 F(TRANSPOS, "TRANSPOS", m_m, NULL) \
314 F(TRUNC, "TRUNC", m_e, NULL) \
315 F(UNIFORM, "UNIFORM", m_ddn, "ai>=0 bi>=0") \
316 F(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0") \
317 F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0") \
318 F(IDF_BETA, "IDF.BETA", m_edd, "a[0,1] b>0 c>0") \
319 F(RV_BETA, "RV.BETA", d_dd, "a>0 b>0") \
320 F(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
321 F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
322 F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]") \
323 F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]") \
324 F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0") \
325 F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0") \
326 F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "c>0") \
327 F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \
328 F(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0") \
329 F(CHICDF, "CHICDF", m_ed, "a>=0 b>0") \
330 F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0") \
331 F(PDF_CHISQ, "PDF.CHISQ", m_ed, "a>=0 b>0") \
332 F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \
333 F(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0") \
334 F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0") \
335 F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0") \
336 F(PDF_EXP, "PDF.EXP", m_ed, "a>=0 b>0") \
337 F(RV_EXP, "RV.EXP", d_d, "a>0") \
338 F(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0") \
339 F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \
340 F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0") \
341 F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0") \
342 F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0") \
343 F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0") \
344 F(RV_F, "RV.F", d_dd, "a>0 b>0") \
345 F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0") \
346 F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
347 F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0") \
348 F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
349 F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0") \
350 F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL) \
351 F(RV_LANDAU, "RV.LANDAU", d_none, NULL) \
352 F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0") \
353 F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0") \
354 F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0") \
355 F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0") \
356 F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]") \
357 F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]") \
358 F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0") \
359 F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0") \
360 F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0") \
361 F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0") \
362 F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
363 F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0") \
364 F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
365 F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0") \
366 F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0") \
367 F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0") \
368 F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0") \
369 F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0") \
370 F(CDFNORM, "CDFNORM", m_e, NULL) \
371 F(PROBIT, "PROBIT", m_e, "a(0,1)") \
372 F(NORMAL, "NORMAL", m_e, "a>0") \
373 F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0") \
374 F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0") \
375 F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0") \
376 F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0") \
377 F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0") \
378 F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0") \
379 F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0") \
380 F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0") \
381 F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0") \
382 F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0") \
383 F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL) \
384 F(RV_RTAIL, "RV.RTAIL", d_dd, NULL) \
385 F(CDF_T, "CDF.T", m_ed, "b>0") \
386 F(TCDF, "TCDF", m_ed, "b>0") \
387 F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0") \
388 F(PDF_T, "PDF.T", m_ed, "b>0") \
389 F(RV_T, "RV.T", d_d, "a>0") \
390 F(CDF_T1G, "CDF.T1G", m_edd, NULL) \
391 F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)") \
392 F(PDF_T1G, "PDF.T1G", m_edd, NULL) \
393 F(RV_T1G, "RV.T1G", d_dd, NULL) \
394 F(CDF_T2G, "CDF.T2G", m_edd, NULL) \
395 F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)") \
396 F(PDF_T2G, "PDF.T2G", m_edd, NULL) \
397 F(RV_T2G, "RV.T2G", d_dd, NULL) \
398 F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c") \
399 F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c") \
400 F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c") \
401 F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b") \
402 F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
403 F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0") \
404 F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
405 F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0") \
406 F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
407 F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
408 F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]") \
409 F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]") \
410 F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]") \
411 F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]") \
412 F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]") \
413 F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]") \
414 F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]") \
415 F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
416 F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
417 F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a") \
418 F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]") \
419 F(RV_LOG, "RV.LOG", d_d, "a(0,1]") \
420 F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
421 F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
422 F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]") \
423 F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0") \
424 F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0") \
425 F(RV_POISSON, "RV.POISSON", d_d, "a>0")
427 /* Properties of a matrix function.
429 These come straight from the macro invocations above. */
430 struct matrix_function_properties
433 const char *constraints;
436 /* Minimum and maximum argument counts for each matrix function prototype. */
437 enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 };
438 enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
439 enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
440 enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 };
441 enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 };
442 enum { d_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 };
443 enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX };
444 enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 };
445 enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 };
446 enum { m_ddn_MIN_ARGS = 2, m_ddn_MAX_ARGS = 2 };
447 enum { m_e_MIN_ARGS = 1, m_e_MAX_ARGS = 1 };
448 enum { m_ed_MIN_ARGS = 2, m_ed_MAX_ARGS = 2 };
449 enum { m_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
450 enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
451 enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
452 enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 };
453 enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 };
454 enum { m_mddn_MIN_ARGS = 3, m_mddn_MAX_ARGS = 3 };
455 enum { m_mdn_MIN_ARGS = 2, m_mdn_MAX_ARGS = 2 };
456 enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 };
457 enum { m_mmn_MIN_ARGS = 2, m_mmn_MAX_ARGS = 2 };
458 enum { m_mn_MIN_ARGS = 1, m_mn_MAX_ARGS = 1 };
459 enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 };
461 /* C function prototype for each matrix function prototype. */
462 typedef double matrix_proto_d_none (void);
463 typedef double matrix_proto_d_d (double);
464 typedef double matrix_proto_d_dd (double, double);
465 typedef double matrix_proto_d_dd (double, double);
466 typedef double matrix_proto_d_ddd (double, double, double);
467 typedef gsl_matrix *matrix_proto_m_d (double);
468 typedef gsl_matrix *matrix_proto_m_ddd (double, double, double);
469 typedef gsl_matrix *matrix_proto_m_ddn (double, double,
470 const struct matrix_expr *);
471 typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *);
472 typedef gsl_matrix *matrix_proto_m_mn (gsl_matrix *,
473 const struct matrix_expr *);
474 typedef double matrix_proto_m_e (double);
475 typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
476 typedef gsl_matrix *matrix_proto_m_mdn (gsl_matrix *, double,
477 const struct matrix_expr *);
478 typedef double matrix_proto_m_ed (double, double);
479 typedef gsl_matrix *matrix_proto_m_mddn (gsl_matrix *, double, double,
480 const struct matrix_expr *);
481 typedef double matrix_proto_m_edd (double, double, double);
482 typedef double matrix_proto_m_eddd (double, double, double, double);
483 typedef double matrix_proto_m_eed (double, double, double);
484 typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *);
485 typedef gsl_matrix *matrix_proto_m_mmn (gsl_matrix *, gsl_matrix *,
486 const struct matrix_expr *);
487 typedef gsl_matrix *matrix_proto_m_v (gsl_vector *);
488 typedef double matrix_proto_d_m (gsl_matrix *);
489 typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n);
490 typedef gsl_matrix *matrix_proto_IDENT (double, double);
492 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
493 static matrix_proto_##PROTO matrix_eval_##ENUM;
497 /* Matrix expression data structure and parsing. */
499 /* A node in a matrix expression. */
505 #define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM,
509 /* Elementwise and scalar arithmetic. */
510 MOP_NEGATE, /* unary - */
511 MOP_ADD_ELEMS, /* + */
512 MOP_SUB_ELEMS, /* - */
513 MOP_MUL_ELEMS, /* &* */
514 MOP_DIV_ELEMS, /* / and &/ */
515 MOP_EXP_ELEMS, /* &** */
517 MOP_SEQ_BY, /* a:b:c */
519 /* Matrix arithmetic. */
521 MOP_EXP_MAT, /* ** */
538 MOP_PASTE_HORZ, /* a, b, c, ... */
539 MOP_PASTE_VERT, /* a; b; c; ... */
543 MOP_VEC_INDEX, /* x(y) */
544 MOP_VEC_ALL, /* x(:) */
545 MOP_MAT_INDEX, /* x(y,z) */
546 MOP_ROW_INDEX, /* x(y,:) */
547 MOP_COL_INDEX, /* x(:,z) */
554 MOP_EOF, /* EOF('file') */
560 /* Nonterminal expression nodes. */
563 struct matrix_expr **subs;
567 /* Terminal expression nodes. */
568 double number; /* MOP_NUMBER. */
569 struct matrix_var *variable; /* MOP_VARIABLE. */
570 struct read_file *eof; /* MOP_EOF. */
573 /* The syntax location corresponding to this expression node, for use in
574 error messages. This is always nonnull for terminal expression nodes.
575 For most others, it is null because it can be computed lazily if and
578 Use matrix_expr_location() instead of using this member directly, so
579 that it gets computed lazily if needed. */
580 struct msg_location *location;
584 matrix_expr_location__ (const struct matrix_expr *e,
585 const struct msg_location **minp,
586 const struct msg_location **maxp)
588 struct msg_location *loc = e->location;
591 const struct msg_location *min = *minp;
594 || loc->start.line < min->start.line
595 || (loc->start.line == min->start.line
596 && loc->start.column < min->start.column)))
599 const struct msg_location *max = *maxp;
602 || loc->end.line > max->end.line
603 || (loc->end.line == max->end.line
604 && loc->end.column > max->end.column)))
610 assert (e->op != MOP_NUMBER && e->op != MOP_VARIABLE && e->op != MOP_EOF);
611 for (size_t i = 0; i < e->n_subs; i++)
612 matrix_expr_location__ (e->subs[i], minp, maxp);
615 /* Returns the source code location corresponding to expression E, computing it
617 static const struct msg_location *
618 matrix_expr_location (const struct matrix_expr *e_)
620 struct matrix_expr *e = CONST_CAST (struct matrix_expr *, e_);
626 const struct msg_location *min = NULL;
627 const struct msg_location *max = NULL;
628 matrix_expr_location__ (e, &min, &max);
631 e->location = msg_location_dup (min);
632 e->location->end = max->end;
638 /* Sets e->location to the tokens in S's lexer from offset START_OFS to the
639 token before the current one. Has no effect if E already has a location or
642 matrix_expr_add_location (struct matrix_state *s, int start_ofs,
643 struct matrix_expr *e)
645 if (e && !e->location)
646 e->location = lex_ofs_location (s->lexer, start_ofs,
647 lex_ofs (s->lexer) - 1);
650 /* Frees E and all the data and sub-expressions that it references. */
652 matrix_expr_destroy (struct matrix_expr *e)
659 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
690 for (size_t i = 0; i < e->n_subs; i++)
691 matrix_expr_destroy (e->subs[i]);
700 msg_location_destroy (e->location);
704 /* Creates and returns a new matrix_expr with type OP, which must be a
705 nonterminal type. Initializes the new matrix_expr with the N_SUBS
706 expressions in SUBS as subexpressions. */
707 static struct matrix_expr *
708 matrix_expr_create_subs (enum matrix_op op, struct matrix_expr **subs,
711 struct matrix_expr *e = xmalloc (sizeof *e);
712 *e = (struct matrix_expr) {
714 .subs = xmemdup (subs, n_subs * sizeof *subs),
720 static struct matrix_expr *
721 matrix_expr_create_0 (enum matrix_op op)
723 struct matrix_expr *sub;
724 return matrix_expr_create_subs (op, &sub, 0);
727 static struct matrix_expr *
728 matrix_expr_create_1 (enum matrix_op op, struct matrix_expr *sub)
730 return matrix_expr_create_subs (op, &sub, 1);
733 static struct matrix_expr *
734 matrix_expr_create_2 (enum matrix_op op,
735 struct matrix_expr *sub0, struct matrix_expr *sub1)
737 struct matrix_expr *subs[] = { sub0, sub1 };
738 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
741 static struct matrix_expr *
742 matrix_expr_create_3 (enum matrix_op op, struct matrix_expr *sub0,
743 struct matrix_expr *sub1, struct matrix_expr *sub2)
745 struct matrix_expr *subs[] = { sub0, sub1, sub2 };
746 return matrix_expr_create_subs (op, subs, sizeof subs / sizeof *subs);
749 /* Creates and returns a new MOP_NUMBER expression node to contain NUMBER. */
750 static struct matrix_expr *
751 matrix_expr_create_number (double number)
753 struct matrix_expr *e = xmalloc (sizeof *e);
754 *e = (struct matrix_expr) {
761 static struct matrix_expr *matrix_expr_parse (struct matrix_state *);
763 /* A binary operator for matrix_parse_binary_operator(). */
764 struct matrix_operator_syntax
766 /* Exactly one of these specifies the operator syntax. */
767 enum token_type token; /* A token, e.g. T_ASTERISK. */
768 const char *id; /* An identifier, e.g. "XOR". */
769 const char *phrase; /* A token phrase, e.g. "&**". */
771 /* The matrix operator corresponding to the syntax. */
776 matrix_operator_syntax_match (struct lexer *lexer,
777 const struct matrix_operator_syntax *syntax,
778 size_t n_syntax, enum matrix_op *op)
780 const struct matrix_operator_syntax *end = &syntax[n_syntax];
781 for (const struct matrix_operator_syntax *syn = syntax; syn < end; syn++)
782 if (syn->id ? lex_match_id (lexer, syn->id)
783 : syn->phrase ? lex_match_phrase (lexer, syn->phrase)
784 : lex_match (lexer, syn->token))
792 /* Parses a binary operator level in the recursive descent parser, returning a
793 matrix expression if successful or a null pointer otherwise. PARSE_NEXT
794 must be the function to parse the next level of precedence. The N_SYNTAX
795 elements of SYNTAX must specify the syntax and matrix_expr node type to
796 parse at this level. */
797 static struct matrix_expr *
798 matrix_parse_binary_operator (
799 struct matrix_state *s,
800 struct matrix_expr *(*parse_next) (struct matrix_state *),
801 const struct matrix_operator_syntax *syntax, size_t n_syntax)
803 struct matrix_expr *lhs = parse_next (s);
810 if (!matrix_operator_syntax_match (s->lexer, syntax, n_syntax, &op))
813 struct matrix_expr *rhs = parse_next (s);
816 matrix_expr_destroy (lhs);
819 lhs = matrix_expr_create_2 (op, lhs, rhs);
823 /* Parses a comma-separated list of expressions within {}, transforming them
824 into MOP_PASTE_HORZ operators. Returns the new expression or NULL on
826 static struct matrix_expr *
827 matrix_parse_curly_comma (struct matrix_state *s)
829 static const struct matrix_operator_syntax op = {
830 .token = T_COMMA, .op = MOP_PASTE_HORZ
832 return matrix_parse_binary_operator (s, matrix_expr_parse, &op, 1);
835 /* Parses a semicolon-separated list of expressions within {}, transforming
836 them into MOP_PASTE_VERT operators. Returns the new expression or NULL on
838 static struct matrix_expr *
839 matrix_parse_curly_semi (struct matrix_state *s)
841 if (lex_token (s->lexer) == T_RCURLY)
843 /* {} is a special case for a 0×0 matrix. */
844 return matrix_expr_create_0 (MOP_EMPTY);
847 static const struct matrix_operator_syntax op = {
848 .token = T_SEMICOLON, .op = MOP_PASTE_VERT
850 return matrix_parse_binary_operator (s, matrix_parse_curly_comma, &op, 1);
853 struct matrix_function
857 size_t min_args, max_args;
860 static struct matrix_expr *matrix_expr_parse (struct matrix_state *);
863 word_matches (const char **test, const char **name)
865 size_t test_len = strcspn (*test, ".");
866 size_t name_len = strcspn (*name, ".");
867 if (test_len == name_len)
869 if (buf_compare_case (*test, *name, test_len))
872 else if (test_len < 3 || test_len > name_len)
876 if (buf_compare_case (*test, *name, test_len))
882 if (**test != **name)
893 /* Returns 0 if TOKEN and FUNC do not match,
894 1 if TOKEN is an acceptable abbreviation for FUNC,
895 2 if TOKEN equals FUNC. */
897 compare_function_names (const char *token_, const char *func_)
899 const char *token = token_;
900 const char *func = func_;
901 while (*token || *func)
902 if (!word_matches (&token, &func))
904 return !c_strcasecmp (token_, func_) ? 2 : 1;
907 static const struct matrix_function *
908 matrix_parse_function_name (const char *token)
910 static const struct matrix_function functions[] =
912 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
913 { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS },
917 enum { N_FUNCTIONS = sizeof functions / sizeof *functions };
919 for (size_t i = 0; i < N_FUNCTIONS; i++)
921 if (compare_function_names (token, functions[i].name) > 0)
922 return &functions[i];
928 matrix_parse_function (struct matrix_state *s, const char *token,
929 struct matrix_expr **exprp)
932 if (lex_next_token (s->lexer, 1) != T_LPAREN)
935 int start_ofs = lex_ofs (s->lexer);
936 if (lex_match_id (s->lexer, "EOF"))
939 struct file_handle *fh = fh_parse (s->lexer, FH_REF_FILE, s->session);
943 if (!lex_force_match (s->lexer, T_RPAREN))
949 struct read_file *rf = read_file_create (s, fh);
951 struct matrix_expr *e = xmalloc (sizeof *e);
952 *e = (struct matrix_expr) { .op = MOP_EOF, .eof = rf };
953 matrix_expr_add_location (s, start_ofs, e);
958 const struct matrix_function *f = matrix_parse_function_name (token);
962 struct matrix_expr *e = xmalloc (sizeof *e);
963 *e = (struct matrix_expr) { .op = f->op };
965 lex_get_n (s->lexer, 2);
966 if (lex_token (s->lexer) != T_RPAREN)
968 size_t allocated_subs = 0;
971 struct matrix_expr *sub = matrix_expr_parse (s);
975 if (e->n_subs >= allocated_subs)
976 e->subs = x2nrealloc (e->subs, &allocated_subs, sizeof *e->subs);
977 e->subs[e->n_subs++] = sub;
979 while (lex_match (s->lexer, T_COMMA));
981 if (!lex_force_match (s->lexer, T_RPAREN))
984 if (e->n_subs < f->min_args || e->n_subs > f->max_args)
986 if (f->min_args == f->max_args)
987 msg_at (SE, e->location,
988 ngettext ("Matrix function %s requires %zu argument.",
989 "Matrix function %s requires %zu arguments.",
991 f->name, f->min_args);
992 else if (f->min_args == 1 && f->max_args == 2)
993 msg_at (SE, e->location,
994 ngettext ("Matrix function %s requires 1 or 2 arguments, "
995 "but %zu was provided.",
996 "Matrix function %s requires 1 or 2 arguments, "
997 "but %zu were provided.",
1000 else if (f->min_args == 1 && f->max_args == INT_MAX)
1001 msg_at (SE, e->location,
1002 _("Matrix function %s requires at least one argument."),
1010 matrix_expr_add_location (s, start_ofs, e);
1016 matrix_expr_destroy (e);
1020 static struct matrix_expr *
1021 matrix_parse_primary__ (struct matrix_state *s)
1023 if (lex_is_number (s->lexer))
1025 double number = lex_number (s->lexer);
1028 return matrix_expr_create_number (number);
1030 else if (lex_is_string (s->lexer))
1032 char string[sizeof (double)];
1033 buf_copy_str_rpad (string, sizeof string, lex_tokcstr (s->lexer), ' ');
1037 memcpy (&number, string, sizeof number);
1039 return matrix_expr_create_number (number);
1041 else if (lex_match (s->lexer, T_LPAREN))
1043 struct matrix_expr *e = matrix_expr_parse (s);
1044 if (!e || !lex_force_match (s->lexer, T_RPAREN))
1046 matrix_expr_destroy (e);
1051 else if (lex_match (s->lexer, T_LCURLY))
1053 struct matrix_expr *e = matrix_parse_curly_semi (s);
1054 if (!e || !lex_force_match (s->lexer, T_RCURLY))
1056 matrix_expr_destroy (e);
1061 else if (lex_token (s->lexer) == T_ID)
1063 struct matrix_expr *retval;
1064 if (matrix_parse_function (s, lex_tokcstr (s->lexer), &retval))
1067 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
1070 lex_error (s->lexer, _("Unknown variable %s."),
1071 lex_tokcstr (s->lexer));
1076 struct matrix_expr *e = xmalloc (sizeof *e);
1077 *e = (struct matrix_expr) { .op = MOP_VARIABLE, .variable = var };
1080 else if (lex_token (s->lexer) == T_ALL)
1082 struct matrix_expr *retval;
1083 if (matrix_parse_function (s, "ALL", &retval))
1087 lex_error (s->lexer, NULL);
1091 static struct matrix_expr *
1092 matrix_parse_primary (struct matrix_state *s)
1094 int start_ofs = lex_ofs (s->lexer);
1095 struct matrix_expr *e = matrix_parse_primary__ (s);
1096 matrix_expr_add_location (s, start_ofs, e);
1100 static struct matrix_expr *matrix_parse_postfix (struct matrix_state *);
1103 matrix_parse_index_expr (struct matrix_state *s,
1104 struct matrix_expr **indexp,
1105 struct msg_location **locationp)
1107 if (lex_match (s->lexer, T_COLON))
1110 *locationp = lex_get_location (s->lexer, -1, -1);
1116 *indexp = matrix_expr_parse (s);
1117 if (locationp && *indexp)
1118 *locationp = msg_location_dup (matrix_expr_location (*indexp));
1119 return *indexp != NULL;
1123 static struct matrix_expr *
1124 matrix_parse_postfix (struct matrix_state *s)
1126 struct matrix_expr *lhs = matrix_parse_primary (s);
1127 if (!lhs || !lex_match (s->lexer, T_LPAREN))
1130 struct matrix_expr *i0;
1131 if (!matrix_parse_index_expr (s, &i0, NULL))
1133 matrix_expr_destroy (lhs);
1136 if (lex_match (s->lexer, T_RPAREN))
1138 ? matrix_expr_create_2 (MOP_VEC_INDEX, lhs, i0)
1139 : matrix_expr_create_1 (MOP_VEC_ALL, lhs));
1140 else if (lex_match (s->lexer, T_COMMA))
1142 struct matrix_expr *i1;
1143 if (!matrix_parse_index_expr (s, &i1, NULL)
1144 || !lex_force_match (s->lexer, T_RPAREN))
1146 matrix_expr_destroy (lhs);
1147 matrix_expr_destroy (i0);
1148 matrix_expr_destroy (i1);
1151 return (i0 && i1 ? matrix_expr_create_3 (MOP_MAT_INDEX, lhs, i0, i1)
1152 : i0 ? matrix_expr_create_2 (MOP_ROW_INDEX, lhs, i0)
1153 : i1 ? matrix_expr_create_2 (MOP_COL_INDEX, lhs, i1)
1158 lex_error_expecting (s->lexer, "`)'", "`,'");
1163 static struct matrix_expr *
1164 matrix_parse_unary (struct matrix_state *s)
1166 int start_ofs = lex_ofs (s->lexer);
1168 struct matrix_expr *e;
1169 if (lex_match (s->lexer, T_DASH))
1171 struct matrix_expr *sub = matrix_parse_unary (s);
1174 e = matrix_expr_create_1 (MOP_NEGATE, sub);
1176 else if (lex_match (s->lexer, T_PLUS))
1178 e = matrix_parse_unary (s);
1183 return matrix_parse_postfix (s);
1185 matrix_expr_add_location (s, start_ofs, e);
1186 e->location->start = lex_ofs_start_point (s->lexer, start_ofs);
1190 static struct matrix_expr *
1191 matrix_parse_seq (struct matrix_state *s)
1193 struct matrix_expr *start = matrix_parse_unary (s);
1194 if (!start || !lex_match (s->lexer, T_COLON))
1197 struct matrix_expr *end = matrix_parse_unary (s);
1200 matrix_expr_destroy (start);
1204 if (lex_match (s->lexer, T_COLON))
1206 struct matrix_expr *increment = matrix_parse_unary (s);
1209 matrix_expr_destroy (start);
1210 matrix_expr_destroy (end);
1213 return matrix_expr_create_3 (MOP_SEQ_BY, start, end, increment);
1216 return matrix_expr_create_2 (MOP_SEQ, start, end);
1219 static struct matrix_expr *
1220 matrix_parse_exp (struct matrix_state *s)
1222 static const struct matrix_operator_syntax syntax[] = {
1223 { .token = T_EXP, .op = MOP_EXP_MAT },
1224 { .phrase = "&**", .op = MOP_EXP_ELEMS },
1226 size_t n_syntax = sizeof syntax / sizeof *syntax;
1228 return matrix_parse_binary_operator (s, matrix_parse_seq, syntax, n_syntax);
1231 static struct matrix_expr *
1232 matrix_parse_mul_div (struct matrix_state *s)
1234 static const struct matrix_operator_syntax syntax[] = {
1235 { .token = T_ASTERISK, .op = MOP_MUL_MAT },
1236 { .token = T_SLASH, .op = MOP_DIV_ELEMS },
1237 { .phrase = "&*", .op = MOP_MUL_ELEMS },
1238 { .phrase = "&/", .op = MOP_DIV_ELEMS },
1240 size_t n_syntax = sizeof syntax / sizeof *syntax;
1242 return matrix_parse_binary_operator (s, matrix_parse_exp, syntax, n_syntax);
1245 static struct matrix_expr *
1246 matrix_parse_add_sub (struct matrix_state *s)
1248 struct matrix_expr *lhs = matrix_parse_mul_div (s);
1255 if (lex_match (s->lexer, T_PLUS))
1257 else if (lex_match (s->lexer, T_DASH))
1259 else if (lex_token (s->lexer) == T_NEG_NUM)
1264 struct matrix_expr *rhs = matrix_parse_mul_div (s);
1267 matrix_expr_destroy (lhs);
1270 lhs = matrix_expr_create_2 (op, lhs, rhs);
1274 static struct matrix_expr *
1275 matrix_parse_relational (struct matrix_state *s)
1277 static const struct matrix_operator_syntax syntax[] = {
1278 { .token = T_GT, .op = MOP_GT },
1279 { .token = T_GE, .op = MOP_GE },
1280 { .token = T_LT, .op = MOP_LT },
1281 { .token = T_LE, .op = MOP_LE },
1282 { .token = T_EQUALS, .op = MOP_EQ },
1283 { .token = T_EQ, .op = MOP_EQ },
1284 { .token = T_NE, .op = MOP_NE },
1286 size_t n_syntax = sizeof syntax / sizeof *syntax;
1288 return matrix_parse_binary_operator (s, matrix_parse_add_sub,
1292 static struct matrix_expr *
1293 matrix_parse_not (struct matrix_state *s)
1295 int start_ofs = lex_ofs (s->lexer);
1296 if (lex_match (s->lexer, T_NOT))
1298 struct matrix_expr *sub = matrix_parse_not (s);
1302 struct matrix_expr *e = matrix_expr_create_1 (MOP_NOT, sub);
1303 matrix_expr_add_location (s, start_ofs, e);
1304 e->location->start = lex_ofs_start_point (s->lexer, start_ofs);
1308 return matrix_parse_relational (s);
1311 static struct matrix_expr *
1312 matrix_parse_and (struct matrix_state *s)
1314 static const struct matrix_operator_syntax op = {
1315 .token = T_AND, .op = MOP_AND
1318 return matrix_parse_binary_operator (s, matrix_parse_not, &op, 1);
1321 static struct matrix_expr *
1322 matrix_expr_parse__ (struct matrix_state *s)
1324 static const struct matrix_operator_syntax syntax[] = {
1325 { .token = T_OR, .op = MOP_OR },
1326 { .id = "XOR", .op = MOP_XOR },
1328 size_t n_syntax = sizeof syntax / sizeof *syntax;
1330 return matrix_parse_binary_operator (s, matrix_parse_and, syntax, n_syntax);
1333 static struct matrix_expr *
1334 matrix_expr_parse (struct matrix_state *s)
1336 int start_ofs = lex_ofs (s->lexer);
1337 struct matrix_expr *e = matrix_expr_parse__ (s);
1338 matrix_expr_add_location (s, start_ofs, e);
1342 /* Matrix expression evaluation. */
1344 /* Iterates over all the elements in matrix M, setting Y and X to the row and
1345 column indexes, respectively, and pointing D to the entry at each
1347 #define MATRIX_FOR_ALL_ELEMENTS(D, Y, X, M) \
1348 for (size_t Y = 0; Y < (M)->size1; Y++) \
1349 for (size_t X = 0; X < (M)->size2; X++) \
1350 for (double *D = gsl_matrix_ptr ((M), Y, X); D; D = NULL)
1353 is_vector (const gsl_matrix *m)
1355 return m->size1 <= 1 || m->size2 <= 1;
1359 to_vector (gsl_matrix *m)
1361 return (m->size1 == 1
1362 ? gsl_matrix_row (m, 0).vector
1363 : gsl_matrix_column (m, 0).vector);
1367 matrix_eval_ABS (double d)
1373 matrix_eval_ALL (gsl_matrix *m)
1375 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1382 matrix_eval_ANY (gsl_matrix *m)
1384 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1391 matrix_eval_ARSIN (double d)
1397 matrix_eval_ARTAN (double d)
1403 matrix_eval_BLOCK (gsl_matrix *m[], size_t n)
1407 for (size_t i = 0; i < n; i++)
1412 gsl_matrix *block = gsl_matrix_calloc (r, c);
1414 for (size_t i = 0; i < n; i++)
1416 for (size_t y = 0; y < m[i]->size1; y++)
1417 for (size_t x = 0; x < m[i]->size2; x++)
1418 gsl_matrix_set (block, r + y, c + x, gsl_matrix_get (m[i], y, x));
1426 matrix_eval_CHOL (gsl_matrix *m, const struct matrix_expr *e)
1428 if (!gsl_linalg_cholesky_decomp1 (m))
1430 for (size_t y = 0; y < m->size1; y++)
1431 for (size_t x = y + 1; x < m->size2; x++)
1432 gsl_matrix_set (m, y, x, gsl_matrix_get (m, x, y));
1434 for (size_t y = 0; y < m->size1; y++)
1435 for (size_t x = 0; x < y; x++)
1436 gsl_matrix_set (m, y, x, 0);
1441 msg_at (SE, e->subs[0]->location,
1442 _("Input to CHOL function is not positive-definite."));
1448 matrix_eval_col_extremum (gsl_matrix *m, bool min)
1453 return gsl_matrix_alloc (1, 0);
1455 gsl_matrix *cext = gsl_matrix_alloc (1, m->size2);
1456 for (size_t x = 0; x < m->size2; x++)
1458 double ext = gsl_matrix_get (m, 0, x);
1459 for (size_t y = 1; y < m->size1; y++)
1461 double value = gsl_matrix_get (m, y, x);
1462 if (min ? value < ext : value > ext)
1465 gsl_matrix_set (cext, 0, x, ext);
1471 matrix_eval_CMAX (gsl_matrix *m)
1473 return matrix_eval_col_extremum (m, false);
1477 matrix_eval_CMIN (gsl_matrix *m)
1479 return matrix_eval_col_extremum (m, true);
1483 matrix_eval_COS (double d)
1489 matrix_eval_col_sum (gsl_matrix *m, bool square)
1494 return gsl_matrix_alloc (1, 0);
1496 gsl_matrix *result = gsl_matrix_alloc (1, m->size2);
1497 for (size_t x = 0; x < m->size2; x++)
1500 for (size_t y = 0; y < m->size1; y++)
1502 double d = gsl_matrix_get (m, y, x);
1503 sum += square ? pow2 (d) : d;
1505 gsl_matrix_set (result, 0, x, sum);
1511 matrix_eval_CSSQ (gsl_matrix *m)
1513 return matrix_eval_col_sum (m, true);
1517 matrix_eval_CSUM (gsl_matrix *m)
1519 return matrix_eval_col_sum (m, false);
1523 compare_double_3way (const void *a_, const void *b_)
1525 const double *a = a_;
1526 const double *b = b_;
1527 return *a < *b ? -1 : *a > *b;
1531 matrix_eval_DESIGN (gsl_matrix *m, const struct matrix_expr *e)
1533 double *tmp = xmalloc (m->size1 * m->size2 * sizeof *tmp);
1534 gsl_matrix m2 = gsl_matrix_view_array (tmp, m->size2, m->size1).matrix;
1535 gsl_matrix_transpose_memcpy (&m2, m);
1537 for (size_t y = 0; y < m2.size1; y++)
1538 qsort (tmp + y * m2.size2, m2.size2, sizeof *tmp, compare_double_3way);
1540 size_t *n = xcalloc (m2.size1, sizeof *n);
1542 for (size_t i = 0; i < m2.size1; i++)
1544 double *row = tmp + m2.size2 * i;
1545 for (size_t j = 0; j < m2.size2; )
1548 for (k = j + 1; k < m2.size2; k++)
1549 if (row[j] != row[k])
1551 row[n[i]++] = row[j];
1556 msg_at (MW, e->subs[0]->location,
1557 _("Column %zu in DESIGN argument has constant value."), i + 1);
1562 gsl_matrix *result = gsl_matrix_alloc (m->size1, n_total);
1564 for (size_t i = 0; i < m->size2; i++)
1569 const double *unique = tmp + m2.size2 * i;
1570 for (size_t j = 0; j < n[i]; j++, x++)
1572 double value = unique[j];
1573 for (size_t y = 0; y < m->size1; y++)
1574 gsl_matrix_set (result, y, x, gsl_matrix_get (m, y, i) == value);
1585 matrix_eval_DET (gsl_matrix *m)
1587 gsl_permutation *p = gsl_permutation_alloc (m->size1);
1589 gsl_linalg_LU_decomp (m, p, &signum);
1590 gsl_permutation_free (p);
1591 return gsl_linalg_LU_det (m, signum);
1595 matrix_eval_DIAG (gsl_matrix *m)
1597 gsl_matrix *diag = gsl_matrix_alloc (MIN (m->size1, m->size2), 1);
1598 for (size_t i = 0; i < diag->size1; i++)
1599 gsl_matrix_set (diag, i, 0, gsl_matrix_get (m, i, i));
1604 is_symmetric (const gsl_matrix *m)
1606 if (m->size1 != m->size2)
1609 for (size_t y = 0; y < m->size1; y++)
1610 for (size_t x = 0; x < y; x++)
1611 if (gsl_matrix_get (m, y, x) != gsl_matrix_get (m, x, y))
1618 compare_double_desc (const void *a_, const void *b_)
1620 const double *a = a_;
1621 const double *b = b_;
1622 return *a > *b ? -1 : *a < *b;
1626 matrix_eval_EVAL (gsl_matrix *m, const struct matrix_expr *e)
1628 if (!is_symmetric (m))
1630 msg_at (SE, e->subs[0]->location,
1631 _("Argument of EVAL must be symmetric."));
1635 gsl_eigen_symm_workspace *w = gsl_eigen_symm_alloc (m->size1);
1636 gsl_matrix *eval = gsl_matrix_alloc (m->size1, 1);
1637 gsl_vector v_eval = to_vector (eval);
1638 gsl_eigen_symm (m, &v_eval, w);
1639 gsl_eigen_symm_free (w);
1641 assert (v_eval.stride == 1);
1642 qsort (v_eval.data, v_eval.size, sizeof *v_eval.data, compare_double_desc);
1648 matrix_eval_EXP (double d)
1653 /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was
1656 Charl Linssen <charl@itfromb.it>
1660 matrix_eval_GINV (gsl_matrix *A)
1662 size_t n = A->size1;
1663 size_t m = A->size2;
1665 gsl_matrix *tmp_mat = NULL;
1668 /* libgsl SVD can only handle the case m <= n, so transpose matrix. */
1669 tmp_mat = gsl_matrix_alloc (m, n);
1670 gsl_matrix_transpose_memcpy (tmp_mat, A);
1678 gsl_matrix *V = gsl_matrix_alloc (m, m);
1679 gsl_vector *u = gsl_vector_alloc (m);
1681 gsl_vector *tmp_vec = gsl_vector_alloc (m);
1682 gsl_linalg_SV_decomp (A, V, u, tmp_vec);
1683 gsl_vector_free (tmp_vec);
1686 gsl_matrix *Sigma_pinv = gsl_matrix_alloc (m, n);
1687 gsl_matrix_set_zero (Sigma_pinv);
1688 double cutoff = 1e-15 * gsl_vector_max (u);
1690 for (size_t i = 0; i < m; ++i)
1692 double x = gsl_vector_get (u, i);
1693 gsl_matrix_set (Sigma_pinv, i, i, x > cutoff ? 1.0 / x : 0);
1696 /* libgsl SVD yields "thin" SVD. Pad to full matrix by adding zeros. */
1697 gsl_matrix *U = gsl_matrix_calloc (n, n);
1698 for (size_t i = 0; i < n; i++)
1699 for (size_t j = 0; j < m; j++)
1700 gsl_matrix_set (U, i, j, gsl_matrix_get (A, i, j));
1702 /* Two dot products to obtain pseudoinverse. */
1703 gsl_matrix *tmp_mat2 = gsl_matrix_alloc (m, n);
1704 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., tmp_mat2);
1709 A_pinv = gsl_matrix_alloc (n, m);
1710 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., U, tmp_mat2, 0., A_pinv);
1714 A_pinv = gsl_matrix_alloc (m, n);
1715 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., tmp_mat2, U, 0., A_pinv);
1718 gsl_matrix_free (tmp_mat);
1719 gsl_matrix_free (tmp_mat2);
1720 gsl_matrix_free (U);
1721 gsl_matrix_free (Sigma_pinv);
1722 gsl_vector_free (u);
1723 gsl_matrix_free (V);
1735 grade_compare_3way (const void *a_, const void *b_)
1737 const struct grade *a = a_;
1738 const struct grade *b = b_;
1740 return (a->value < b->value ? -1
1741 : a->value > b->value ? 1
1749 matrix_eval_GRADE (gsl_matrix *m)
1751 size_t n = m->size1 * m->size2;
1752 struct grade *grades = xmalloc (n * sizeof *grades);
1755 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1756 grades[i++] = (struct grade) { .y = y, .x = x, .value = *d };
1757 qsort (grades, n, sizeof *grades, grade_compare_3way);
1759 for (size_t i = 0; i < n; i++)
1760 gsl_matrix_set (m, grades[i].y, grades[i].x, i + 1);
1768 dot (gsl_vector *a, gsl_vector *b)
1770 double result = 0.0;
1771 for (size_t i = 0; i < a->size; i++)
1772 result += gsl_vector_get (a, i) * gsl_vector_get (b, i);
1777 norm2 (gsl_vector *v)
1779 double result = 0.0;
1780 for (size_t i = 0; i < v->size; i++)
1781 result += pow2 (gsl_vector_get (v, i));
1786 norm (gsl_vector *v)
1788 return sqrt (norm2 (v));
1792 matrix_eval_GSCH (gsl_matrix *v, const struct matrix_expr *e)
1794 if (v->size2 < v->size1)
1796 msg_at (SE, e->subs[0]->location,
1797 _("GSCH requires its argument to have at least as many columns "
1798 "as rows, but it has dimensions %zu×%zu."),
1799 v->size1, v->size2);
1802 if (!v->size1 || !v->size2)
1805 gsl_matrix *u = gsl_matrix_calloc (v->size1, v->size2);
1807 for (size_t vx = 0; vx < v->size2; vx++)
1809 gsl_vector u_i = gsl_matrix_column (u, ux).vector;
1810 gsl_vector v_i = gsl_matrix_column (v, vx).vector;
1812 gsl_vector_memcpy (&u_i, &v_i);
1813 for (size_t j = 0; j < ux; j++)
1815 gsl_vector u_j = gsl_matrix_column (u, j).vector;
1816 double scale = dot (&u_j, &u_i) / norm2 (&u_j);
1817 for (size_t k = 0; k < u_i.size; k++)
1818 gsl_vector_set (&u_i, k, (gsl_vector_get (&u_i, k)
1819 - scale * gsl_vector_get (&u_j, k)));
1822 double len = norm (&u_i);
1825 gsl_vector_scale (&u_i, 1.0 / len);
1826 if (++ux >= v->size1)
1833 msg_at (SE, e->subs[0]->location,
1834 _("%zu×%zu argument to GSCH contains only "
1835 "%zu linearly independent columns."),
1836 v->size1, v->size2, ux);
1837 gsl_matrix_free (u);
1841 u->size2 = v->size1;
1846 matrix_eval_IDENT (double s1, double s2)
1848 gsl_matrix *m = gsl_matrix_alloc (s1, s2);
1849 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
1854 /* Inverts X, storing the inverse into INVERSE. As a side effect, replaces X
1855 by its LU decomposition. */
1857 invert_matrix (gsl_matrix *x, gsl_matrix *inverse)
1859 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1861 gsl_linalg_LU_decomp (x, p, &signum);
1862 gsl_linalg_LU_invert (x, p, inverse);
1863 gsl_permutation_free (p);
1867 matrix_eval_INV (gsl_matrix *src)
1869 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
1870 invert_matrix (src, dst);
1875 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1877 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1878 a->size2 * b->size2);
1880 for (size_t ar = 0; ar < a->size1; ar++)
1881 for (size_t br = 0; br < b->size1; br++, y++)
1884 for (size_t ac = 0; ac < a->size2; ac++)
1885 for (size_t bc = 0; bc < b->size2; bc++, x++)
1887 double av = gsl_matrix_get (a, ar, ac);
1888 double bv = gsl_matrix_get (b, br, bc);
1889 gsl_matrix_set (k, y, x, av * bv);
1896 matrix_eval_LG10 (double d)
1902 matrix_eval_LN (double d)
1908 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1910 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1913 for (size_t i = 1; i <= n * n; i++)
1915 gsl_matrix_set (m, y, x, i);
1917 size_t y1 = !y ? n - 1 : y - 1;
1918 size_t x1 = x + 1 >= n ? 0 : x + 1;
1919 if (gsl_matrix_get (m, y1, x1) == 0)
1925 y = y + 1 >= n ? 0 : y + 1;
1930 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1932 double a = gsl_matrix_get (m, y1, x1);
1933 double b = gsl_matrix_get (m, y2, x2);
1934 gsl_matrix_set (m, y1, x1, b);
1935 gsl_matrix_set (m, y2, x2, a);
1939 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1943 /* A. Umar, "On the Construction of Even Order Magic Squares",
1944 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1946 for (size_t i = 1; i <= n * n / 2; i++)
1948 gsl_matrix_set (m, y, x, i);
1958 for (size_t i = n * n; i > n * n / 2; i--)
1960 gsl_matrix_set (m, y, x, i);
1968 for (size_t y = 0; y < n; y++)
1969 for (size_t x = 0; x < n / 2; x++)
1971 unsigned int d = gsl_matrix_get (m, y, x);
1972 if (d % 2 != (y < n / 2))
1973 magic_exchange (m, y, x, y, n - x - 1);
1978 size_t x1 = n / 2 - 1;
1980 magic_exchange (m, y1, x1, y2, x1);
1981 magic_exchange (m, y1, x2, y2, x2);
1985 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1987 /* A. Umar, "On the Construction of Even Order Magic Squares",
1988 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1992 for (size_t i = 1; ; i++)
1994 gsl_matrix_set (m, y, x, i);
1995 if (++y == n / 2 - 1)
2007 for (size_t i = n * n; ; i--)
2009 gsl_matrix_set (m, y, x, i);
2010 if (++y == n / 2 - 1)
2019 for (size_t y = 0; y < n; y++)
2020 if (y != n / 2 - 1 && y != n / 2)
2021 for (size_t x = 0; x < n / 2; x++)
2023 unsigned int d = gsl_matrix_get (m, y, x);
2024 if (d % 2 != (y < n / 2))
2025 magic_exchange (m, y, x, y, n - x - 1);
2028 size_t a0 = (n * n - 2 * n) / 2 + 1;
2029 for (size_t i = 0; i < n / 2; i++)
2032 gsl_matrix_set (m, n / 2, i, a);
2033 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
2035 for (size_t i = 0; i < n / 2; i++)
2037 size_t a = a0 + i + n / 2;
2038 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
2039 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
2041 for (size_t x = 1; x < n / 2; x += 2)
2042 magic_exchange (m, n / 2, x, n / 2 - 1, x);
2043 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
2044 magic_exchange (m, n / 2, x, n / 2 - 1, x);
2045 size_t x1 = n / 2 - 2;
2046 size_t x2 = n / 2 + 1;
2047 size_t y1 = n / 2 - 2;
2048 size_t y2 = n / 2 + 1;
2049 magic_exchange (m, y1, x1, y2, x1);
2050 magic_exchange (m, y1, x2, y2, x2);
2054 matrix_eval_MAGIC (double n_)
2058 gsl_matrix *m = gsl_matrix_calloc (n, n);
2060 matrix_eval_MAGIC_odd (m, n);
2062 matrix_eval_MAGIC_singly_even (m, n);
2064 matrix_eval_MAGIC_doubly_even (m, n);
2069 matrix_eval_MAKE (double r, double c, double value)
2071 gsl_matrix *m = gsl_matrix_alloc (r, c);
2072 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2078 matrix_eval_MDIAG (gsl_vector *v)
2080 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
2081 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
2082 gsl_vector_memcpy (&diagonal, v);
2087 matrix_eval_MMAX (gsl_matrix *m)
2089 return gsl_matrix_max (m);
2093 matrix_eval_MMIN (gsl_matrix *m)
2095 return gsl_matrix_min (m);
2099 matrix_eval_MOD (gsl_matrix *m, double divisor)
2101 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2102 *d = fmod (*d, divisor);
2107 matrix_eval_MSSQ (gsl_matrix *m)
2110 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2116 matrix_eval_MSUM (gsl_matrix *m)
2119 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2125 matrix_eval_NCOL (gsl_matrix *m)
2131 matrix_eval_NROW (gsl_matrix *m)
2137 matrix_eval_RANK (gsl_matrix *m)
2139 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
2140 gsl_linalg_QR_decomp (m, tau);
2141 gsl_vector_free (tau);
2143 return gsl_linalg_QRPT_rank (m, -1);
2147 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_,
2148 const struct matrix_expr *e)
2150 bool r_ok = r_ >= 0 && r_ < SIZE_MAX;
2151 bool c_ok = c_ >= 0 && c_ < SIZE_MAX;
2155 !r_ok ? e->subs[1]->location : e->subs[2]->location,
2156 _("Arguments 2 and 3 to RESHAPE must be integers."));
2161 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
2163 struct msg_location *loc = msg_location_dup (e->subs[1]->location);
2164 loc->end = e->subs[2]->location->end;
2165 msg_at (SE, loc, _("Product of RESHAPE size arguments (%zu×%zu = %zu) "
2166 "differs from product of matrix dimensions "
2167 "(%zu×%zu = %zu)."),
2169 m->size1, m->size2, m->size1 * m->size2);
2170 msg_location_destroy (loc);
2174 gsl_matrix *dst = gsl_matrix_alloc (r, c);
2177 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
2179 gsl_matrix_set (dst, y1, x1, *d);
2190 matrix_eval_row_extremum (gsl_matrix *m, bool min)
2195 return gsl_matrix_alloc (0, 1);
2197 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
2198 for (size_t y = 0; y < m->size1; y++)
2200 double ext = gsl_matrix_get (m, y, 0);
2201 for (size_t x = 1; x < m->size2; x++)
2203 double value = gsl_matrix_get (m, y, x);
2204 if (min ? value < ext : value > ext)
2207 gsl_matrix_set (rext, y, 0, ext);
2213 matrix_eval_RMAX (gsl_matrix *m)
2215 return matrix_eval_row_extremum (m, false);
2219 matrix_eval_RMIN (gsl_matrix *m)
2221 return matrix_eval_row_extremum (m, true);
2225 matrix_eval_RND (double d)
2237 rank_compare_3way (const void *a_, const void *b_)
2239 const struct rank *a = a_;
2240 const struct rank *b = b_;
2242 return a->value < b->value ? -1 : a->value > b->value;
2246 matrix_eval_RNKORDER (gsl_matrix *m)
2248 size_t n = m->size1 * m->size2;
2249 struct rank *ranks = xmalloc (n * sizeof *ranks);
2251 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2252 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
2253 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
2255 for (size_t i = 0; i < n; )
2258 for (j = i + 1; j < n; j++)
2259 if (ranks[i].value != ranks[j].value)
2262 double rank = (i + j + 1.0) / 2.0;
2263 for (size_t k = i; k < j; k++)
2264 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
2275 matrix_eval_row_sum (gsl_matrix *m, bool square)
2280 return gsl_matrix_alloc (0, 1);
2282 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
2283 for (size_t y = 0; y < m->size1; y++)
2286 for (size_t x = 0; x < m->size2; x++)
2288 double d = gsl_matrix_get (m, y, x);
2289 sum += square ? pow2 (d) : d;
2291 gsl_matrix_set (result, y, 0, sum);
2297 matrix_eval_RSSQ (gsl_matrix *m)
2299 return matrix_eval_row_sum (m, true);
2303 matrix_eval_RSUM (gsl_matrix *m)
2305 return matrix_eval_row_sum (m, false);
2309 matrix_eval_SIN (double d)
2315 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2, const struct matrix_expr *e)
2317 if (m1->size1 != m2->size1)
2319 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
2320 loc->end = e->subs[1]->location->end;
2322 msg_at (SE, e->location,
2323 _("SOLVE arguments must have the same number of rows."));
2324 msg_at (SN, e->subs[0]->location,
2325 _("Argument 1 has dimensions %zu×%zu."), m1->size1, m1->size2);
2326 msg_at (SN, e->subs[1]->location,
2327 _("Argument 2 has dimensions %zu×%zu."), m2->size1, m2->size2);
2329 msg_location_destroy (loc);
2333 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
2334 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
2336 gsl_linalg_LU_decomp (m1, p, &signum);
2337 for (size_t i = 0; i < m2->size2; i++)
2339 gsl_vector bi = gsl_matrix_column (m2, i).vector;
2340 gsl_vector xi = gsl_matrix_column (x, i).vector;
2341 gsl_linalg_LU_solve (m1, p, &bi, &xi);
2343 gsl_permutation_free (p);
2348 matrix_eval_SQRT (double d)
2354 matrix_eval_SSCP (gsl_matrix *m)
2356 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
2357 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
2362 matrix_eval_SVAL (gsl_matrix *m)
2364 gsl_matrix *tmp_mat = NULL;
2365 if (m->size2 > m->size1)
2367 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
2368 gsl_matrix_transpose_memcpy (tmp_mat, m);
2373 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
2374 gsl_vector *S = gsl_vector_alloc (m->size2);
2375 gsl_vector *work = gsl_vector_alloc (m->size2);
2376 gsl_linalg_SV_decomp (m, V, S, work);
2378 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
2379 for (size_t i = 0; i < m->size2; i++)
2380 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
2382 gsl_matrix_free (V);
2383 gsl_vector_free (S);
2384 gsl_vector_free (work);
2385 gsl_matrix_free (tmp_mat);
2391 matrix_eval_SWEEP (gsl_matrix *m, double d, const struct matrix_expr *e)
2393 if (d < 1 || d > SIZE_MAX)
2395 msg_at (SE, e->subs[1]->location,
2396 _("Scalar argument to SWEEP must be integer."));
2400 if (k >= MIN (m->size1, m->size2))
2402 msg_at (SE, e->subs[1]->location,
2403 _("Scalar argument to SWEEP must be integer less than or "
2404 "equal to the smaller of the matrix argument's rows and "
2409 double m_kk = gsl_matrix_get (m, k, k);
2410 if (fabs (m_kk) > 1e-19)
2412 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
2413 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
2415 double m_ij = gsl_matrix_get (m, i, j);
2416 double m_ik = gsl_matrix_get (m, i, k);
2417 double m_kj = gsl_matrix_get (m, k, j);
2418 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
2427 for (size_t i = 0; i < m->size1; i++)
2429 gsl_matrix_set (m, i, k, 0);
2430 gsl_matrix_set (m, k, i, 0);
2437 matrix_eval_TRACE (gsl_matrix *m)
2440 size_t n = MIN (m->size1, m->size2);
2441 for (size_t i = 0; i < n; i++)
2442 sum += gsl_matrix_get (m, i, i);
2447 matrix_eval_T (gsl_matrix *m)
2449 return matrix_eval_TRANSPOS (m);
2453 matrix_eval_TRANSPOS (gsl_matrix *m)
2455 if (m->size1 == m->size2)
2457 gsl_matrix_transpose (m);
2462 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
2463 gsl_matrix_transpose_memcpy (t, m);
2469 matrix_eval_TRUNC (double d)
2475 matrix_eval_UNIFORM (double r_, double c_, const struct matrix_expr *e)
2479 if (size_overflow_p (xtimes (r, xmax (c, 1))))
2481 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
2482 loc->end = e->subs[1]->location->end;
2485 _("Product of arguments to UNIFORM exceeds memory size."));
2487 msg_location_destroy (loc);
2491 gsl_matrix *m = gsl_matrix_alloc (r, c);
2492 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2493 *d = gsl_ran_flat (get_rng (), 0, 1);
2498 matrix_eval_PDF_BETA (double x, double a, double b)
2500 return gsl_ran_beta_pdf (x, a, b);
2504 matrix_eval_CDF_BETA (double x, double a, double b)
2506 return gsl_cdf_beta_P (x, a, b);
2510 matrix_eval_IDF_BETA (double P, double a, double b)
2512 return gsl_cdf_beta_Pinv (P, a, b);
2516 matrix_eval_RV_BETA (double a, double b)
2518 return gsl_ran_beta (get_rng (), a, b);
2522 matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
2524 return ncdf_beta (x, a, b, lambda);
2528 matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
2530 return npdf_beta (x, a, b, lambda);
2534 matrix_eval_CDF_BVNOR (double x0, double x1, double r)
2536 return cdf_bvnor (x0, x1, r);
2540 matrix_eval_PDF_BVNOR (double x0, double x1, double r)
2542 return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
2546 matrix_eval_CDF_CAUCHY (double x, double a, double b)
2548 return gsl_cdf_cauchy_P ((x - a) / b, 1);
2552 matrix_eval_IDF_CAUCHY (double P, double a, double b)
2554 return a + b * gsl_cdf_cauchy_Pinv (P, 1);
2558 matrix_eval_PDF_CAUCHY (double x, double a, double b)
2560 return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
2564 matrix_eval_RV_CAUCHY (double a, double b)
2566 return a + b * gsl_ran_cauchy (get_rng (), 1);
2570 matrix_eval_CDF_CHISQ (double x, double df)
2572 return gsl_cdf_chisq_P (x, df);
2576 matrix_eval_CHICDF (double x, double df)
2578 return matrix_eval_CDF_CHISQ (x, df);
2582 matrix_eval_IDF_CHISQ (double P, double df)
2584 return gsl_cdf_chisq_Pinv (P, df);
2588 matrix_eval_PDF_CHISQ (double x, double df)
2590 return gsl_ran_chisq_pdf (x, df);
2594 matrix_eval_RV_CHISQ (double df)
2596 return gsl_ran_chisq (get_rng (), df);
2600 matrix_eval_SIG_CHISQ (double x, double df)
2602 return gsl_cdf_chisq_Q (x, df);
2606 matrix_eval_CDF_EXP (double x, double a)
2608 return gsl_cdf_exponential_P (x, 1. / a);
2612 matrix_eval_IDF_EXP (double P, double a)
2614 return gsl_cdf_exponential_Pinv (P, 1. / a);
2618 matrix_eval_PDF_EXP (double x, double a)
2620 return gsl_ran_exponential_pdf (x, 1. / a);
2624 matrix_eval_RV_EXP (double a)
2626 return gsl_ran_exponential (get_rng (), 1. / a);
2630 matrix_eval_PDF_XPOWER (double x, double a, double b)
2632 return gsl_ran_exppow_pdf (x, a, b);
2636 matrix_eval_RV_XPOWER (double a, double b)
2638 return gsl_ran_exppow (get_rng (), a, b);
2642 matrix_eval_CDF_F (double x, double df1, double df2)
2644 return gsl_cdf_fdist_P (x, df1, df2);
2648 matrix_eval_FCDF (double x, double df1, double df2)
2650 return matrix_eval_CDF_F (x, df1, df2);
2654 matrix_eval_IDF_F (double P, double df1, double df2)
2656 return idf_fdist (P, df1, df2);
2660 matrix_eval_RV_F (double df1, double df2)
2662 return gsl_ran_fdist (get_rng (), df1, df2);
2666 matrix_eval_PDF_F (double x, double df1, double df2)
2668 return gsl_ran_fdist_pdf (x, df1, df2);
2672 matrix_eval_SIG_F (double x, double df1, double df2)
2674 return gsl_cdf_fdist_Q (x, df1, df2);
2678 matrix_eval_CDF_GAMMA (double x, double a, double b)
2680 return gsl_cdf_gamma_P (x, a, 1. / b);
2684 matrix_eval_IDF_GAMMA (double P, double a, double b)
2686 return gsl_cdf_gamma_Pinv (P, a, 1. / b);
2690 matrix_eval_PDF_GAMMA (double x, double a, double b)
2692 return gsl_ran_gamma_pdf (x, a, 1. / b);
2696 matrix_eval_RV_GAMMA (double a, double b)
2698 return gsl_ran_gamma (get_rng (), a, 1. / b);
2702 matrix_eval_PDF_LANDAU (double x)
2704 return gsl_ran_landau_pdf (x);
2708 matrix_eval_RV_LANDAU (void)
2710 return gsl_ran_landau (get_rng ());
2714 matrix_eval_CDF_LAPLACE (double x, double a, double b)
2716 return gsl_cdf_laplace_P ((x - a) / b, 1);
2720 matrix_eval_IDF_LAPLACE (double P, double a, double b)
2722 return a + b * gsl_cdf_laplace_Pinv (P, 1);
2726 matrix_eval_PDF_LAPLACE (double x, double a, double b)
2728 return gsl_ran_laplace_pdf ((x - a) / b, 1);
2732 matrix_eval_RV_LAPLACE (double a, double b)
2734 return a + b * gsl_ran_laplace (get_rng (), 1);
2738 matrix_eval_RV_LEVY (double c, double alpha)
2740 return gsl_ran_levy (get_rng (), c, alpha);
2744 matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
2746 return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
2750 matrix_eval_CDF_LOGISTIC (double x, double a, double b)
2752 return gsl_cdf_logistic_P ((x - a) / b, 1);
2756 matrix_eval_IDF_LOGISTIC (double P, double a, double b)
2758 return a + b * gsl_cdf_logistic_Pinv (P, 1);
2762 matrix_eval_PDF_LOGISTIC (double x, double a, double b)
2764 return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
2768 matrix_eval_RV_LOGISTIC (double a, double b)
2770 return a + b * gsl_ran_logistic (get_rng (), 1);
2774 matrix_eval_CDF_LNORMAL (double x, double m, double s)
2776 return gsl_cdf_lognormal_P (x, log (m), s);
2780 matrix_eval_IDF_LNORMAL (double P, double m, double s)
2782 return gsl_cdf_lognormal_Pinv (P, log (m), s);;
2786 matrix_eval_PDF_LNORMAL (double x, double m, double s)
2788 return gsl_ran_lognormal_pdf (x, log (m), s);
2792 matrix_eval_RV_LNORMAL (double m, double s)
2794 return gsl_ran_lognormal (get_rng (), log (m), s);
2798 matrix_eval_CDF_NORMAL (double x, double u, double s)
2800 return gsl_cdf_gaussian_P (x - u, s);
2804 matrix_eval_IDF_NORMAL (double P, double u, double s)
2806 return u + gsl_cdf_gaussian_Pinv (P, s);
2810 matrix_eval_PDF_NORMAL (double x, double u, double s)
2812 return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
2816 matrix_eval_RV_NORMAL (double u, double s)
2818 return u + gsl_ran_gaussian (get_rng (), s);
2822 matrix_eval_CDFNORM (double x)
2824 return gsl_cdf_ugaussian_P (x);
2828 matrix_eval_PROBIT (double P)
2830 return gsl_cdf_ugaussian_Pinv (P);
2834 matrix_eval_NORMAL (double s)
2836 return gsl_ran_gaussian (get_rng (), s);
2840 matrix_eval_PDF_NTAIL (double x, double a, double sigma)
2842 return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
2846 matrix_eval_RV_NTAIL (double a, double sigma)
2848 return gsl_ran_gaussian_tail (get_rng (), a, sigma);
2852 matrix_eval_CDF_PARETO (double x, double a, double b)
2854 return gsl_cdf_pareto_P (x, b, a);
2858 matrix_eval_IDF_PARETO (double P, double a, double b)
2860 return gsl_cdf_pareto_Pinv (P, b, a);
2864 matrix_eval_PDF_PARETO (double x, double a, double b)
2866 return gsl_ran_pareto_pdf (x, b, a);
2870 matrix_eval_RV_PARETO (double a, double b)
2872 return gsl_ran_pareto (get_rng (), b, a);
2876 matrix_eval_CDF_RAYLEIGH (double x, double sigma)
2878 return gsl_cdf_rayleigh_P (x, sigma);
2882 matrix_eval_IDF_RAYLEIGH (double P, double sigma)
2884 return gsl_cdf_rayleigh_Pinv (P, sigma);
2888 matrix_eval_PDF_RAYLEIGH (double x, double sigma)
2890 return gsl_ran_rayleigh_pdf (x, sigma);
2894 matrix_eval_RV_RAYLEIGH (double sigma)
2896 return gsl_ran_rayleigh (get_rng (), sigma);
2900 matrix_eval_PDF_RTAIL (double x, double a, double sigma)
2902 return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
2906 matrix_eval_RV_RTAIL (double a, double sigma)
2908 return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
2912 matrix_eval_CDF_T (double x, double df)
2914 return gsl_cdf_tdist_P (x, df);
2918 matrix_eval_TCDF (double x, double df)
2920 return matrix_eval_CDF_T (x, df);
2924 matrix_eval_IDF_T (double P, double df)
2926 return gsl_cdf_tdist_Pinv (P, df);
2930 matrix_eval_PDF_T (double x, double df)
2932 return gsl_ran_tdist_pdf (x, df);
2936 matrix_eval_RV_T (double df)
2938 return gsl_ran_tdist (get_rng (), df);
2942 matrix_eval_CDF_T1G (double x, double a, double b)
2944 return gsl_cdf_gumbel1_P (x, a, b);
2948 matrix_eval_IDF_T1G (double P, double a, double b)
2950 return gsl_cdf_gumbel1_Pinv (P, a, b);
2954 matrix_eval_PDF_T1G (double x, double a, double b)
2956 return gsl_ran_gumbel1_pdf (x, a, b);
2960 matrix_eval_RV_T1G (double a, double b)
2962 return gsl_ran_gumbel1 (get_rng (), a, b);
2966 matrix_eval_CDF_T2G (double x, double a, double b)
2968 return gsl_cdf_gumbel1_P (x, a, b);
2972 matrix_eval_IDF_T2G (double P, double a, double b)
2974 return gsl_cdf_gumbel1_Pinv (P, a, b);
2978 matrix_eval_PDF_T2G (double x, double a, double b)
2980 return gsl_ran_gumbel1_pdf (x, a, b);
2984 matrix_eval_RV_T2G (double a, double b)
2986 return gsl_ran_gumbel1 (get_rng (), a, b);
2990 matrix_eval_CDF_UNIFORM (double x, double a, double b)
2992 return gsl_cdf_flat_P (x, a, b);
2996 matrix_eval_IDF_UNIFORM (double P, double a, double b)
2998 return gsl_cdf_flat_Pinv (P, a, b);
3002 matrix_eval_PDF_UNIFORM (double x, double a, double b)
3004 return gsl_ran_flat_pdf (x, a, b);
3008 matrix_eval_RV_UNIFORM (double a, double b)
3010 return gsl_ran_flat (get_rng (), a, b);
3014 matrix_eval_CDF_WEIBULL (double x, double a, double b)
3016 return gsl_cdf_weibull_P (x, a, b);
3020 matrix_eval_IDF_WEIBULL (double P, double a, double b)
3022 return gsl_cdf_weibull_Pinv (P, a, b);
3026 matrix_eval_PDF_WEIBULL (double x, double a, double b)
3028 return gsl_ran_weibull_pdf (x, a, b);
3032 matrix_eval_RV_WEIBULL (double a, double b)
3034 return gsl_ran_weibull (get_rng (), a, b);
3038 matrix_eval_CDF_BERNOULLI (double k, double p)
3040 return k ? 1 : 1 - p;
3044 matrix_eval_PDF_BERNOULLI (double k, double p)
3046 return gsl_ran_bernoulli_pdf (k, p);
3050 matrix_eval_RV_BERNOULLI (double p)
3052 return gsl_ran_bernoulli (get_rng (), p);
3056 matrix_eval_CDF_BINOM (double k, double n, double p)
3058 return gsl_cdf_binomial_P (k, p, n);
3062 matrix_eval_PDF_BINOM (double k, double n, double p)
3064 return gsl_ran_binomial_pdf (k, p, n);
3068 matrix_eval_RV_BINOM (double n, double p)
3070 return gsl_ran_binomial (get_rng (), p, n);
3074 matrix_eval_CDF_GEOM (double k, double p)
3076 return gsl_cdf_geometric_P (k, p);
3080 matrix_eval_PDF_GEOM (double k, double p)
3082 return gsl_ran_geometric_pdf (k, p);
3086 matrix_eval_RV_GEOM (double p)
3088 return gsl_ran_geometric (get_rng (), p);
3092 matrix_eval_CDF_HYPER (double k, double a, double b, double c)
3094 return gsl_cdf_hypergeometric_P (k, c, a - c, b);
3098 matrix_eval_PDF_HYPER (double k, double a, double b, double c)
3100 return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
3104 matrix_eval_RV_HYPER (double a, double b, double c)
3106 return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
3110 matrix_eval_PDF_LOG (double k, double p)
3112 return gsl_ran_logarithmic_pdf (k, p);
3116 matrix_eval_RV_LOG (double p)
3118 return gsl_ran_logarithmic (get_rng (), p);
3122 matrix_eval_CDF_NEGBIN (double k, double n, double p)
3124 return gsl_cdf_negative_binomial_P (k, p, n);
3128 matrix_eval_PDF_NEGBIN (double k, double n, double p)
3130 return gsl_ran_negative_binomial_pdf (k, p, n);
3134 matrix_eval_RV_NEGBIN (double n, double p)
3136 return gsl_ran_negative_binomial (get_rng (), p, n);
3140 matrix_eval_CDF_POISSON (double k, double mu)
3142 return gsl_cdf_poisson_P (k, mu);
3146 matrix_eval_PDF_POISSON (double k, double mu)
3148 return gsl_ran_poisson_pdf (k, mu);
3152 matrix_eval_RV_POISSON (double mu)
3154 return gsl_ran_poisson (get_rng (), mu);
3158 matrix_op_eval (enum matrix_op op, double a, double b)
3162 case MOP_ADD_ELEMS: return a + b;
3163 case MOP_SUB_ELEMS: return a - b;
3164 case MOP_MUL_ELEMS: return a * b;
3165 case MOP_DIV_ELEMS: return a / b;
3166 case MOP_EXP_ELEMS: return pow (a, b);
3167 case MOP_GT: return a > b;
3168 case MOP_GE: return a >= b;
3169 case MOP_LT: return a < b;
3170 case MOP_LE: return a <= b;
3171 case MOP_EQ: return a == b;
3172 case MOP_NE: return a != b;
3173 case MOP_AND: return (a > 0) && (b > 0);
3174 case MOP_OR: return (a > 0) || (b > 0);
3175 case MOP_XOR: return (a > 0) != (b > 0);
3177 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3186 case MOP_PASTE_HORZ:
3187 case MOP_PASTE_VERT:
3203 matrix_op_name (enum matrix_op op)
3207 case MOP_ADD_ELEMS: return "+";
3208 case MOP_SUB_ELEMS: return "-";
3209 case MOP_MUL_ELEMS: return "&*";
3210 case MOP_DIV_ELEMS: return "&/";
3211 case MOP_EXP_ELEMS: return "&**";
3212 case MOP_GT: return ">";
3213 case MOP_GE: return ">=";
3214 case MOP_LT: return "<";
3215 case MOP_LE: return "<=";
3216 case MOP_EQ: return "=";
3217 case MOP_NE: return "<>";
3218 case MOP_AND: return "AND";
3219 case MOP_OR: return "OR";
3220 case MOP_XOR: return "XOR";
3222 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3231 case MOP_PASTE_HORZ:
3232 case MOP_PASTE_VERT:
3248 is_scalar (const gsl_matrix *m)
3250 return m->size1 == 1 && m->size2 == 1;
3254 to_scalar (const gsl_matrix *m)
3256 assert (is_scalar (m));
3257 return gsl_matrix_get (m, 0, 0);
3261 matrix_expr_evaluate_elementwise (const struct matrix_expr *e,
3263 gsl_matrix *a, gsl_matrix *b)
3267 double be = to_scalar (b);
3268 for (size_t r = 0; r < a->size1; r++)
3269 for (size_t c = 0; c < a->size2; c++)
3271 double *ae = gsl_matrix_ptr (a, r, c);
3272 *ae = matrix_op_eval (op, *ae, be);
3276 else if (is_scalar (a))
3278 double ae = to_scalar (a);
3279 for (size_t r = 0; r < b->size1; r++)
3280 for (size_t c = 0; c < b->size2; c++)
3282 double *be = gsl_matrix_ptr (b, r, c);
3283 *be = matrix_op_eval (op, ae, *be);
3287 else if (a->size1 == b->size1 && a->size2 == b->size2)
3289 for (size_t r = 0; r < a->size1; r++)
3290 for (size_t c = 0; c < a->size2; c++)
3292 double *ae = gsl_matrix_ptr (a, r, c);
3293 double be = gsl_matrix_get (b, r, c);
3294 *ae = matrix_op_eval (op, *ae, be);
3300 msg_at (SE, matrix_expr_location (e),
3301 _("The operands of %s must have the same dimensions or one "
3302 "must be a scalar."),
3303 matrix_op_name (op));
3304 msg_at (SN, matrix_expr_location (e->subs[0]),
3305 _("The left-hand operand is a %zu×%zu matrix."),
3306 a->size1, a->size2);
3307 msg_at (SN, matrix_expr_location (e->subs[1]),
3308 _("The right-hand operand is a %zu×%zu matrix."),
3309 b->size1, b->size2);
3315 matrix_expr_evaluate_mul_mat (const struct matrix_expr *e,
3316 gsl_matrix *a, gsl_matrix *b)
3318 if (is_scalar (a) || is_scalar (b))
3319 return matrix_expr_evaluate_elementwise (e, MOP_MUL_ELEMS, a, b);
3321 if (a->size2 != b->size1)
3323 msg_at (SE, e->location,
3324 _("Matrices not conformable for multiplication."));
3325 msg_at (SN, matrix_expr_location (e->subs[0]),
3326 _("The left-hand operand is a %zu×%zu matrix."),
3327 a->size1, a->size2);
3328 msg_at (SN, matrix_expr_location (e->subs[1]),
3329 _("The right-hand operand is a %zu×%zu matrix."),
3330 b->size1, b->size2);
3334 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3335 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3340 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3342 gsl_matrix *tmp = *a;
3348 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3351 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3352 swap_matrix (z, tmp);
3356 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3358 mul_matrix (x, *x, *x, tmp);
3362 matrix_expr_evaluate_exp_mat (const struct matrix_expr *e,
3363 gsl_matrix *x_, gsl_matrix *b)
3366 if (x->size1 != x->size2)
3368 msg_at (SE, matrix_expr_location (e->subs[0]),
3369 _("Matrix exponentation with ** requires a square matrix on "
3370 "the left-hand size, not one with dimensions %zu×%zu."),
3371 x->size1, x->size2);
3376 msg_at (SE, matrix_expr_location (e->subs[1]),
3377 _("Matrix exponentiation with ** requires a scalar on the "
3378 "right-hand side, not a matrix with dimensions %zu×%zu."),
3379 b->size1, b->size2);
3382 double bf = to_scalar (b);
3383 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3385 msg_at (SE, matrix_expr_location (e->subs[1]),
3386 _("Exponent %.1f in matrix multiplication is non-integer "
3387 "or outside the valid range."), bf);
3392 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3394 gsl_matrix_set_identity (y);
3398 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3400 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3403 mul_matrix (&y, x, y, &t);
3404 square_matrix (&x, &t);
3407 square_matrix (&x, &t);
3409 mul_matrix (&y, x, y, &t);
3412 invert_matrix (y, x);
3413 swap_matrix (&x, &y);
3416 /* Garbage collection.
3418 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3419 a permutation of them. We are returning one of them; that one must not be
3420 destroyed. We must not destroy 'x_' because the caller owns it. */
3422 gsl_matrix_free (y_);
3424 gsl_matrix_free (t_);
3430 note_operand_size (const gsl_matrix *m, const struct matrix_expr *e)
3432 msg_at (SN, matrix_expr_location (e),
3433 _("This operand is a %zu×%zu matrix."), m->size1, m->size2);
3437 note_nonscalar (const gsl_matrix *m, const struct matrix_expr *e)
3440 note_operand_size (m, e);
3444 matrix_expr_evaluate_seq (const struct matrix_expr *e,
3445 gsl_matrix *start_, gsl_matrix *end_,
3448 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3450 msg_at (SE, matrix_expr_location (e),
3451 _("All operands of : operator must be scalars."));
3453 note_nonscalar (start_, e->subs[0]);
3454 note_nonscalar (end_, e->subs[1]);
3456 note_nonscalar (by_, e->subs[2]);
3460 long int start = to_scalar (start_);
3461 long int end = to_scalar (end_);
3462 long int by = by_ ? to_scalar (by_) : 1;
3466 msg_at (SE, matrix_expr_location (e->subs[2]),
3467 _("The increment operand to : must be nonzero."));
3471 long int n = (end >= start && by > 0 ? (end - start + by) / by
3472 : end <= start && by < 0 ? (start - end - by) / -by
3474 gsl_matrix *m = gsl_matrix_alloc (1, n);
3475 for (long int i = 0; i < n; i++)
3476 gsl_matrix_set (m, 0, i, start + i * by);
3481 matrix_expr_evaluate_not (gsl_matrix *a)
3483 MATRIX_FOR_ALL_ELEMENTS (d, y, x, a)
3489 matrix_expr_evaluate_paste_horz (const struct matrix_expr *e,
3490 gsl_matrix *a, gsl_matrix *b)
3492 if (a->size1 != b->size1)
3494 if (!a->size1 || !a->size2)
3496 else if (!b->size1 || !b->size2)
3499 msg_at (SE, matrix_expr_location (e),
3500 _("This expression tries to horizontally join matrices with "
3501 "differing numbers of rows."));
3502 note_operand_size (a, e->subs[0]);
3503 note_operand_size (b, e->subs[1]);
3507 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3508 for (size_t y = 0; y < a->size1; y++)
3510 for (size_t x = 0; x < a->size2; x++)
3511 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3512 for (size_t x = 0; x < b->size2; x++)
3513 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3519 matrix_expr_evaluate_paste_vert (const struct matrix_expr *e,
3520 gsl_matrix *a, gsl_matrix *b)
3522 if (a->size2 != b->size2)
3524 if (!a->size1 || !a->size2)
3526 else if (!b->size1 || !b->size2)
3529 msg_at (SE, matrix_expr_location (e),
3530 _("This expression tries to vertically join matrices with "
3531 "differing numbers of columns."));
3532 note_operand_size (a, e->subs[0]);
3533 note_operand_size (b, e->subs[1]);
3537 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3538 for (size_t x = 0; x < a->size2; x++)
3540 for (size_t y = 0; y < a->size1; y++)
3541 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3542 for (size_t y = 0; y < b->size1; y++)
3543 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3549 matrix_to_vector (gsl_matrix *m)
3552 gsl_vector v = to_vector (m);
3553 assert (v.block == m->block || !v.block);
3557 gsl_matrix_free (m);
3558 return xmemdup (&v, sizeof v);
3572 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3575 index_vector_uninit (struct index_vector *iv)
3582 matrix_normalize_index_vector (const gsl_matrix *m,
3583 const struct matrix_expr *me, size_t size,
3584 enum index_type index_type, size_t other_size,
3585 struct index_vector *iv)
3594 msg_at (SE, matrix_expr_location (me),
3595 _("Vector index must be scalar or vector, not a "
3597 m->size1, m->size2);
3601 msg_at (SE, matrix_expr_location (me),
3602 _("Matrix row index must be scalar or vector, not a "
3604 m->size1, m->size2);
3608 msg_at (SE, matrix_expr_location (me),
3609 _("Matrix column index must be scalar or vector, not a "
3611 m->size1, m->size2);
3617 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3618 *iv = (struct index_vector) {
3619 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3622 for (size_t i = 0; i < v.size; i++)
3624 double index = gsl_vector_get (&v, i);
3625 if (index < 1 || index >= size + 1)
3630 msg_at (SE, matrix_expr_location (me),
3631 _("Index %g is out of range for vector "
3632 "with %zu elements."), index, size);
3636 msg_at (SE, matrix_expr_location (me),
3637 _("%g is not a valid row index for "
3638 "a %zu×%zu matrix."),
3639 index, size, other_size);
3643 msg_at (SE, matrix_expr_location (me),
3644 _("%g is not a valid column index for "
3645 "a %zu×%zu matrix."),
3646 index, other_size, size);
3650 index_vector_uninit (iv);
3653 iv->indexes[i] = index - 1;
3659 *iv = (struct index_vector) {
3660 .indexes = xnmalloc (size, sizeof *iv->indexes),
3663 for (size_t i = 0; i < size; i++)
3670 matrix_expr_evaluate_vec_all (const struct matrix_expr *e,
3673 if (!is_vector (sm))
3675 msg_at (SE, matrix_expr_location (e->subs[0]),
3676 _("Vector index operator may not be applied to "
3677 "a %zu×%zu matrix."),
3678 sm->size1, sm->size2);
3686 matrix_expr_evaluate_vec_index (const struct matrix_expr *e,
3687 gsl_matrix *sm, gsl_matrix *im)
3689 if (!matrix_expr_evaluate_vec_all (e, sm))
3692 gsl_vector sv = to_vector (sm);
3693 struct index_vector iv;
3694 if (!matrix_normalize_index_vector (im, e->subs[1],
3695 sv.size, IV_VECTOR, 0, &iv))
3698 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3699 sm->size1 == 1 ? iv.n : 1);
3700 gsl_vector dv = to_vector (dm);
3701 for (size_t dx = 0; dx < iv.n; dx++)
3703 size_t sx = iv.indexes[dx];
3704 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3706 index_vector_uninit (&iv);
3712 matrix_expr_evaluate_mat_index (gsl_matrix *sm,
3713 gsl_matrix *im0, const struct matrix_expr *eim0,
3714 gsl_matrix *im1, const struct matrix_expr *eim1)
3716 struct index_vector iv0;
3717 if (!matrix_normalize_index_vector (im0, eim0, sm->size1,
3718 IV_ROW, sm->size2, &iv0))
3721 struct index_vector iv1;
3722 if (!matrix_normalize_index_vector (im1, eim1, sm->size2,
3723 IV_COLUMN, sm->size1, &iv1))
3725 index_vector_uninit (&iv0);
3729 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3730 for (size_t dy = 0; dy < iv0.n; dy++)
3732 size_t sy = iv0.indexes[dy];
3734 for (size_t dx = 0; dx < iv1.n; dx++)
3736 size_t sx = iv1.indexes[dx];
3737 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3740 index_vector_uninit (&iv0);
3741 index_vector_uninit (&iv1);
3745 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3746 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3747 const struct matrix_function_properties *, gsl_matrix *[], \
3748 const struct matrix_expr *, matrix_proto_##PROTO *);
3753 check_scalar_arg (const char *name, gsl_matrix *subs[],
3754 const struct matrix_expr *e, size_t index)
3756 if (!is_scalar (subs[index]))
3758 msg_at (SE, matrix_expr_location (e->subs[index]),
3759 _("Function %s argument %zu must be a scalar, "
3760 "not a %zu×%zu matrix."),
3761 name, index + 1, subs[index]->size1, subs[index]->size2);
3768 check_vector_arg (const char *name, gsl_matrix *subs[],
3769 const struct matrix_expr *e, size_t index)
3771 if (!is_vector (subs[index]))
3773 msg_at (SE, matrix_expr_location (e->subs[index]),
3774 _("Function %s argument %zu must be a vector, "
3775 "not a %zu×%zu matrix."),
3776 name, index + 1, subs[index]->size1, subs[index]->size2);
3783 to_scalar_args (const char *name, gsl_matrix *subs[],
3784 const struct matrix_expr *e, double d[])
3786 for (size_t i = 0; i < e->n_subs; i++)
3788 if (!check_scalar_arg (name, subs, e, i))
3790 d[i] = to_scalar (subs[i]);
3796 parse_constraint_value (const char **constraintsp)
3799 long retval = strtol (*constraintsp, &tail, 10);
3800 assert (tail > *constraintsp);
3801 *constraintsp = tail;
3805 enum matrix_argument_relop
3815 argument_inequality_error (
3816 const struct matrix_function_properties *props, const struct matrix_expr *e,
3817 size_t ai, gsl_matrix *a, size_t y, size_t x,
3818 size_t bi, double b,
3819 enum matrix_argument_relop relop)
3821 const struct msg_location *loc = matrix_expr_location (e);
3825 msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
3826 "than or equal to argument %zu."),
3827 ai + 1, props->name, bi + 1);
3831 msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
3832 "than argument %zu."),
3833 ai + 1, props->name, bi + 1);
3837 msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
3838 "or equal to argument %zu."),
3839 ai + 1, props->name, bi + 1);
3843 msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
3845 ai + 1, props->name, bi + 1);
3849 msg_at (ME, loc, _("Argument %zu to matrix function %s must not be equal "
3850 "to argument %zu."),
3851 ai + 1, props->name, bi + 1);
3855 const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
3857 msg_at (SN, a_loc, _("Argument %zu is %g."),
3858 ai + 1, gsl_matrix_get (a, y, x));
3860 msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
3861 y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
3863 msg_at (SN, matrix_expr_location (e->subs[bi]),
3864 _("Argument %zu is %g."), bi + 1, b);
3868 argument_value_error (
3869 const struct matrix_function_properties *props, const struct matrix_expr *e,
3870 size_t ai, gsl_matrix *a, size_t y, size_t x,
3872 enum matrix_argument_relop relop)
3874 const struct msg_location *loc = matrix_expr_location (e);
3878 msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
3879 "than or equal to %g."),
3880 ai + 1, props->name, b);
3884 msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
3886 ai + 1, props->name, b);
3890 msg_at (SE, loc, _("Argument %zu to matrix function %s must be less than "
3892 ai + 1, props->name, b);
3896 msg_at (SE, loc, _("Argument %zu to matrix function %s must be less than "
3898 ai + 1, props->name, b);
3902 msg_at (SE, loc, _("Argument %zu to matrix function %s must not be equal "
3904 ai + 1, props->name, b);
3908 const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
3911 if (relop != MRR_NE)
3912 msg_at (SN, a_loc, _("Argument %zu is %g."),
3913 ai + 1, gsl_matrix_get (a, y, x));
3916 msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
3917 y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
3921 matrix_argument_relop_is_satisfied (double a, double b,
3922 enum matrix_argument_relop relop)
3926 case MRR_GE: return a >= b;
3927 case MRR_GT: return a > b;
3928 case MRR_LE: return a <= b;
3929 case MRR_LT: return a < b;
3930 case MRR_NE: return a != b;
3936 static enum matrix_argument_relop
3937 matrix_argument_relop_flip (enum matrix_argument_relop relop)
3941 case MRR_GE: return MRR_LE;
3942 case MRR_GT: return MRR_LT;
3943 case MRR_LE: return MRR_GE;
3944 case MRR_LT: return MRR_GT;
3945 case MRR_NE: return MRR_NE;
3952 check_constraints (const struct matrix_function_properties *props,
3953 gsl_matrix *args[], const struct matrix_expr *e)
3955 size_t n_args = e->n_subs;
3956 const char *constraints = props->constraints;
3960 size_t arg_index = SIZE_MAX;
3961 while (*constraints)
3963 if (*constraints >= 'a' && *constraints <= 'd')
3965 arg_index = *constraints++ - 'a';
3966 assert (arg_index < n_args);
3968 else if (*constraints == '[' || *constraints == '(')
3970 assert (arg_index < n_args);
3971 bool open_lower = *constraints++ == '(';
3972 int minimum = parse_constraint_value (&constraints);
3973 assert (*constraints == ',');
3975 int maximum = parse_constraint_value (&constraints);
3976 assert (*constraints == ']' || *constraints == ')');
3977 bool open_upper = *constraints++ == ')';
3979 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3980 if ((open_lower ? *d <= minimum : *d < minimum)
3981 || (open_upper ? *d >= maximum : *d > maximum))
3983 if (!is_scalar (args[arg_index]))
3984 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
3985 _("Row %zu, column %zu of argument %zu to matrix "
3986 "function %s is %g, which is outside "
3987 "the valid range %c%d,%d%c."),
3988 y + 1, x + 1, arg_index + 1, props->name, *d,
3989 open_lower ? '(' : '[',
3991 open_upper ? ')' : ']');
3993 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
3994 _("Argument %zu to matrix function %s is %g, "
3995 "which is outside the valid range %c%d,%d%c."),
3996 arg_index + 1, props->name, *d,
3997 open_lower ? '(' : '[',
3999 open_upper ? ')' : ']');
4003 else if (*constraints == 'i')
4006 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
4007 if (*d != floor (*d))
4009 if (!is_scalar (args[arg_index]))
4010 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
4011 _("Argument %zu to matrix function %s, which must be "
4012 "integer, contains non-integer value %g in "
4013 "row %zu, column %zu."),
4014 arg_index + 1, props->name, *d, y + 1, x + 1);
4016 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
4017 _("Argument %zu to matrix function %s, which must be "
4018 "integer, has non-integer value %g."),
4019 arg_index + 1, props->name, *d);
4023 else if (*constraints == '>'
4024 || *constraints == '<'
4025 || *constraints == '!')
4027 enum matrix_argument_relop relop;
4028 switch (*constraints++)
4031 if (*constraints == '=')
4041 if (*constraints == '=')
4051 assert (*constraints == '=');
4060 if (*constraints >= 'a' && *constraints <= 'd')
4062 size_t a_index = arg_index;
4063 size_t b_index = *constraints - 'a';
4064 assert (a_index < n_args);
4065 assert (b_index < n_args);
4067 /* We only support one of the two arguments being non-scalar.
4068 It's easier to support only the first one being non-scalar, so
4069 flip things around if it's the other way. */
4070 if (!is_scalar (args[b_index]))
4072 assert (is_scalar (args[a_index]));
4073 size_t tmp_index = a_index;
4075 b_index = tmp_index;
4076 relop = matrix_argument_relop_flip (relop);
4079 double b = to_scalar (args[b_index]);
4080 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
4081 if (!matrix_argument_relop_is_satisfied (*a, b, relop))
4083 argument_inequality_error (
4085 a_index, args[a_index], y, x,
4093 int comparand = parse_constraint_value (&constraints);
4095 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
4096 if (!matrix_argument_relop_is_satisfied (*d, comparand, relop))
4098 argument_value_error (
4100 arg_index, args[arg_index], y, x,
4109 assert (*constraints == ' ');
4111 arg_index = SIZE_MAX;
4118 matrix_expr_evaluate_d_none (const struct matrix_function_properties *props,
4119 gsl_matrix *subs[], const struct matrix_expr *e,
4120 matrix_proto_d_none *f)
4122 assert (e->n_subs == 0);
4124 if (!check_constraints (props, subs, e))
4127 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4128 gsl_matrix_set (m, 0, 0, f ());
4133 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
4134 gsl_matrix *subs[], const struct matrix_expr *e,
4135 matrix_proto_d_d *f)
4137 assert (e->n_subs == 1);
4140 if (!to_scalar_args (props->name, subs, e, &d)
4141 || !check_constraints (props, subs, e))
4144 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4145 gsl_matrix_set (m, 0, 0, f (d));
4150 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
4151 gsl_matrix *subs[], const struct matrix_expr *e,
4152 matrix_proto_d_dd *f)
4154 assert (e->n_subs == 2);
4157 if (!to_scalar_args (props->name, subs, e, d)
4158 && !check_constraints (props, subs, e))
4161 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4162 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
4167 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
4168 gsl_matrix *subs[], const struct matrix_expr *e,
4169 matrix_proto_d_ddd *f)
4171 assert (e->n_subs == 3);
4174 if (!to_scalar_args (props->name, subs, e, d)
4175 || !check_constraints (props, subs, e))
4178 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4179 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
4184 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
4185 gsl_matrix *subs[], const struct matrix_expr *e,
4186 matrix_proto_m_d *f)
4188 assert (e->n_subs == 1);
4191 return (to_scalar_args (props->name, subs, e, &d)
4192 && check_constraints (props, subs, e)
4198 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
4199 gsl_matrix *subs[], const struct matrix_expr *e,
4200 matrix_proto_m_ddd *f)
4202 assert (e->n_subs == 3);
4205 return (to_scalar_args (props->name, subs, e, d)
4206 && check_constraints (props, subs, e)
4207 ? f(d[0], d[1], d[2])
4212 matrix_expr_evaluate_m_ddn (const struct matrix_function_properties *props,
4213 gsl_matrix *subs[], const struct matrix_expr *e,
4214 matrix_proto_m_ddn *f)
4216 assert (e->n_subs == 2);
4219 return (to_scalar_args (props->name, subs, e, d)
4220 && check_constraints (props, subs, e)
4226 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props,
4227 gsl_matrix *subs[], const struct matrix_expr *e,
4228 matrix_proto_m_m *f)
4230 assert (e->n_subs == 1);
4231 return check_constraints (props, subs, e) ? f (subs[0]) : NULL;
4235 matrix_expr_evaluate_m_mn (const struct matrix_function_properties *props,
4236 gsl_matrix *subs[], const struct matrix_expr *e,
4237 matrix_proto_m_mn *f)
4239 assert (e->n_subs == 1);
4240 return check_constraints (props, subs, e) ? f (subs[0], e) : NULL;
4244 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
4245 gsl_matrix *subs[], const struct matrix_expr *e,
4246 matrix_proto_m_e *f)
4248 assert (e->n_subs == 1);
4250 if (!check_constraints (props, subs, e))
4253 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4259 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props,
4260 gsl_matrix *subs[], const struct matrix_expr *e,
4261 matrix_proto_m_md *f)
4263 assert (e->n_subs == 2);
4264 return (check_scalar_arg (props->name, subs, e, 1)
4265 && check_constraints (props, subs, e)
4266 ? f (subs[0], to_scalar (subs[1]))
4271 matrix_expr_evaluate_m_mdn (const struct matrix_function_properties *props,
4272 gsl_matrix *subs[], const struct matrix_expr *e,
4273 matrix_proto_m_mdn *f)
4275 assert (e->n_subs == 2);
4276 return (check_scalar_arg (props->name, subs, e, 1)
4277 && check_constraints (props, subs, e)
4278 ? f (subs[0], to_scalar (subs[1]), e)
4283 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
4284 gsl_matrix *subs[], const struct matrix_expr *e,
4285 matrix_proto_m_ed *f)
4287 assert (e->n_subs == 2);
4288 if (!check_scalar_arg (props->name, subs, e, 1)
4289 || !check_constraints (props, subs, e))
4292 double b = to_scalar (subs[1]);
4293 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4299 matrix_expr_evaluate_m_mddn (const struct matrix_function_properties *props,
4300 gsl_matrix *subs[], const struct matrix_expr *e,
4301 matrix_proto_m_mddn *f)
4303 assert (e->n_subs == 3);
4304 if (!check_scalar_arg (props->name, subs, e, 1)
4305 || !check_scalar_arg (props->name, subs, e, 2)
4306 || !check_constraints (props, subs, e))
4308 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]), e);
4312 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
4313 gsl_matrix *subs[], const struct matrix_expr *e,
4314 matrix_proto_m_edd *f)
4316 assert (e->n_subs == 3);
4317 if (!check_scalar_arg (props->name, subs, e, 1)
4318 || !check_scalar_arg (props->name, subs, e, 2)
4319 || !check_constraints (props, subs, e))
4322 double b = to_scalar (subs[1]);
4323 double c = to_scalar (subs[2]);
4324 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4330 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4331 gsl_matrix *subs[], const struct matrix_expr *e,
4332 matrix_proto_m_eddd *f)
4334 assert (e->n_subs == 4);
4335 for (size_t i = 1; i < 4; i++)
4336 if (!check_scalar_arg (props->name, subs, e, i))
4339 if (!check_constraints (props, subs, e))
4342 double b = to_scalar (subs[1]);
4343 double c = to_scalar (subs[2]);
4344 double d = to_scalar (subs[3]);
4345 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4346 *a = f (*a, b, c, d);
4351 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4352 gsl_matrix *subs[], const struct matrix_expr *e,
4353 matrix_proto_m_eed *f)
4355 assert (e->n_subs == 3);
4356 if (!check_scalar_arg (props->name, subs, e, 2))
4359 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4360 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4362 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
4363 loc->end = e->subs[1]->location->end;
4366 _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4367 "%zu×%zu, but %s requires these arguments either to have "
4368 "the same dimensions or for one of them to be a scalar."),
4370 subs[0]->size1, subs[0]->size2,
4371 subs[1]->size1, subs[1]->size2,
4374 msg_location_destroy (loc);
4378 if (!check_constraints (props, subs, e))
4381 double c = to_scalar (subs[2]);
4383 if (is_scalar (subs[0]))
4385 double a = to_scalar (subs[0]);
4386 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4392 double b = to_scalar (subs[1]);
4393 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4400 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props,
4401 gsl_matrix *subs[], const struct matrix_expr *e,
4402 matrix_proto_m_mm *f)
4404 assert (e->n_subs == 2);
4405 return check_constraints (props, subs, e) ? f (subs[0], subs[1]) : NULL;
4409 matrix_expr_evaluate_m_mmn (const struct matrix_function_properties *props,
4410 gsl_matrix *subs[], const struct matrix_expr *e,
4411 matrix_proto_m_mmn *f)
4413 assert (e->n_subs == 2);
4414 return check_constraints (props, subs, e) ? f (subs[0], subs[1], e) : NULL;
4418 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4419 gsl_matrix *subs[], const struct matrix_expr *e,
4420 matrix_proto_m_v *f)
4422 assert (e->n_subs == 1);
4423 if (!check_vector_arg (props->name, subs, e, 0)
4424 || !check_constraints (props, subs, e))
4426 gsl_vector v = to_vector (subs[0]);
4431 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props,
4432 gsl_matrix *subs[], const struct matrix_expr *e,
4433 matrix_proto_d_m *f)
4435 assert (e->n_subs == 1);
4437 if (!check_constraints (props, subs, e))
4440 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4441 gsl_matrix_set (m, 0, 0, f (subs[0]));
4446 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props,
4447 gsl_matrix *subs[], const struct matrix_expr *e,
4448 matrix_proto_m_any *f)
4450 return check_constraints (props, subs, e) ? f (subs, e->n_subs) : NULL;
4454 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props_ UNUSED,
4455 gsl_matrix *subs[], const struct matrix_expr *e,
4456 matrix_proto_IDENT *f)
4458 static const struct matrix_function_properties p1 = {
4460 .constraints = "ai>=0"
4462 static const struct matrix_function_properties p2 = {
4464 .constraints = "ai>=0 bi>=0"
4466 const struct matrix_function_properties *props = e->n_subs == 1 ? &p1 : &p2;
4468 assert (e->n_subs <= 2);
4471 return (to_scalar_args (props->name, subs, e, d)
4472 && check_constraints (props, subs, e)
4473 ? f (d[0], d[e->n_subs - 1])
4478 matrix_expr_evaluate (const struct matrix_expr *e)
4480 if (e->op == MOP_NUMBER)
4482 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4483 gsl_matrix_set (m, 0, 0, e->number);
4486 else if (e->op == MOP_VARIABLE)
4488 const gsl_matrix *src = e->variable->value;
4491 msg_at (SE, e->location,
4492 _("Uninitialized variable %s used in expression."),
4497 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4498 gsl_matrix_memcpy (dst, src);
4501 else if (e->op == MOP_EOF)
4503 struct dfm_reader *reader = read_file_open (e->eof);
4504 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4505 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4509 enum { N_LOCAL = 3 };
4510 gsl_matrix *local_subs[N_LOCAL];
4511 gsl_matrix **subs = (e->n_subs < N_LOCAL
4513 : xmalloc (e->n_subs * sizeof *subs));
4515 for (size_t i = 0; i < e->n_subs; i++)
4517 subs[i] = matrix_expr_evaluate (e->subs[i]);
4520 for (size_t j = 0; j < i; j++)
4521 gsl_matrix_free (subs[j]);
4522 if (subs != local_subs)
4528 gsl_matrix *result = NULL;
4531 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4532 case MOP_F_##ENUM: \
4534 static const struct matrix_function_properties props = { \
4536 .constraints = CONSTRAINTS, \
4538 result = matrix_expr_evaluate_##PROTO (&props, subs, e, \
4539 matrix_eval_##ENUM); \
4546 gsl_matrix_scale (subs[0], -1.0);
4564 result = matrix_expr_evaluate_elementwise (e, e->op, subs[0], subs[1]);
4568 result = matrix_expr_evaluate_not (subs[0]);
4572 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], NULL);
4576 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], subs[2]);
4580 result = matrix_expr_evaluate_mul_mat (e, subs[0], subs[1]);
4584 result = matrix_expr_evaluate_exp_mat (e, subs[0], subs[1]);
4587 case MOP_PASTE_HORZ:
4588 result = matrix_expr_evaluate_paste_horz (e, subs[0], subs[1]);
4591 case MOP_PASTE_VERT:
4592 result = matrix_expr_evaluate_paste_vert (e, subs[0], subs[1]);
4596 result = gsl_matrix_alloc (0, 0);
4600 result = matrix_expr_evaluate_vec_index (e, subs[0], subs[1]);
4604 result = matrix_expr_evaluate_vec_all (e, subs[0]);
4608 result = matrix_expr_evaluate_mat_index (subs[0],
4609 subs[1], e->subs[1],
4610 subs[2], e->subs[2]);
4614 result = matrix_expr_evaluate_mat_index (subs[0],
4615 subs[1], e->subs[1],
4620 result = matrix_expr_evaluate_mat_index (subs[0],
4622 subs[1], e->subs[1]);
4631 for (size_t i = 0; i < e->n_subs; i++)
4632 if (subs[i] != result)
4633 gsl_matrix_free (subs[i]);
4634 if (subs != local_subs)
4640 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4643 gsl_matrix *m = matrix_expr_evaluate (e);
4649 msg_at (SE, matrix_expr_location (e),
4650 _("Expression for %s must evaluate to scalar, "
4651 "not a %zu×%zu matrix."),
4652 context, m->size1, m->size2);
4653 gsl_matrix_free (m);
4658 gsl_matrix_free (m);
4663 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4667 if (!matrix_expr_evaluate_scalar (e, context, &d))
4671 if (d < LONG_MIN || d > LONG_MAX)
4673 msg_at (SE, matrix_expr_location (e),
4674 _("Expression for %s is outside the integer range."), context);
4683 An lvalue is an expression that can appear on the left side of a COMPUTE
4684 command and in other contexts that assign values.
4686 An lvalue is parsed once, with matrix_lvalue_parse(). It can then be
4687 evaluated (with matrix_lvalue_evaluate()) and assigned (with
4688 matrix_lvalue_assign()).
4690 There are three kinds of lvalues:
4692 - A variable name. A variable used as an lvalue need not be initialized,
4693 since the assignment will initialize.
4695 - A subvector, e.g. "var(index0)". The variable must be initialized and
4696 must have the form of a vector (it must have 1 column or 1 row).
4698 - A submatrix, e.g. "var(index0, index1)". The variable must be
4700 struct matrix_lvalue
4702 struct matrix_var *var; /* Destination variable. */
4703 struct matrix_expr *indexes[2]; /* Index expressions, if any. */
4704 size_t n_indexes; /* Number of indexes. */
4706 struct msg_location *var_location; /* Variable name. */
4707 struct msg_location *full_location; /* Variable name plus indexing. */
4708 struct msg_location *index_locations[2]; /* Index expressions. */
4713 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4717 msg_location_destroy (lvalue->var_location);
4718 msg_location_destroy (lvalue->full_location);
4719 for (size_t i = 0; i < lvalue->n_indexes; i++)
4721 matrix_expr_destroy (lvalue->indexes[i]);
4722 msg_location_destroy (lvalue->index_locations[i]);
4728 /* Parses and returns an lvalue at the current position in S's lexer. Returns
4729 null on parse failure. On success, the caller must eventually free the
4731 static struct matrix_lvalue *
4732 matrix_lvalue_parse (struct matrix_state *s)
4734 if (!lex_force_id (s->lexer))
4737 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4738 int start_ofs = lex_ofs (s->lexer);
4739 lvalue->var_location = lex_get_location (s->lexer, 0, 0);
4740 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4741 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4745 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
4749 lex_get_n (s->lexer, 2);
4751 if (!matrix_parse_index_expr (s, &lvalue->indexes[0],
4752 &lvalue->index_locations[0]))
4754 lvalue->n_indexes++;
4756 if (lex_match (s->lexer, T_COMMA))
4758 if (!matrix_parse_index_expr (s, &lvalue->indexes[1],
4759 &lvalue->index_locations[1]))
4761 lvalue->n_indexes++;
4763 if (!lex_force_match (s->lexer, T_RPAREN))
4766 lvalue->full_location = lex_ofs_location (s->lexer, start_ofs,
4767 lex_ofs (s->lexer) - 1);
4772 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4778 matrix_lvalue_destroy (lvalue);
4783 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4784 enum index_type index_type, size_t other_size,
4785 struct index_vector *iv)
4790 m = matrix_expr_evaluate (e);
4797 bool ok = matrix_normalize_index_vector (m, e, size, index_type,
4799 gsl_matrix_free (m);
4803 /* Evaluates the indexes in LVALUE into IV0 and IV1, owned by the caller.
4804 Returns true if successful, false if evaluating the expressions failed or if
4805 LVALUE otherwise can't be used for an assignment.
4807 On success, the caller retains ownership of the index vectors, which are
4808 suitable for passing to matrix_lvalue_assign(). If not used for that
4809 purpose then they need to eventually be freed (with
4810 index_vector_uninit()). */
4812 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4813 struct index_vector *iv0,
4814 struct index_vector *iv1)
4816 *iv0 = INDEX_VECTOR_INIT;
4817 *iv1 = INDEX_VECTOR_INIT;
4818 if (!lvalue->n_indexes)
4821 /* Validate destination matrix exists and has the right shape. */
4822 gsl_matrix *dm = lvalue->var->value;
4825 msg_at (SE, lvalue->var_location,
4826 _("Undefined variable %s."), lvalue->var->name);
4829 else if (dm->size1 == 0 || dm->size2 == 0)
4831 msg_at (SE, lvalue->full_location, _("Cannot index %zu×%zu matrix %s."),
4832 dm->size1, dm->size2, lvalue->var->name);
4835 else if (lvalue->n_indexes == 1)
4837 if (!is_vector (dm))
4839 msg_at (SE, lvalue->full_location,
4840 _("Can't use vector indexing on %zu×%zu matrix %s."),
4841 dm->size1, dm->size2, lvalue->var->name);
4844 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4845 MAX (dm->size1, dm->size2),
4850 assert (lvalue->n_indexes == 2);
4851 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4852 IV_ROW, dm->size2, iv0))
4855 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4856 IV_COLUMN, dm->size1, iv1))
4858 index_vector_uninit (iv0);
4866 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4867 struct index_vector *iv,
4868 gsl_matrix *sm, const struct msg_location *lsm)
4870 /* Convert source matrix 'sm' to source vector 'sv'. */
4871 if (!is_vector (sm))
4873 msg_at (SE, lvalue->full_location,
4874 _("Only an %zu-element vector may be assigned to this "
4875 "%zu-element subvector of %s."),
4876 iv->n, iv->n, lvalue->var->name);
4878 _("The source is an %zu×%zu matrix."),
4879 sm->size1, sm->size2);
4882 gsl_vector sv = to_vector (sm);
4883 if (iv->n != sv.size)
4885 msg_at (SE, lvalue->full_location,
4886 _("Only an %zu-element vector may be assigned to this "
4887 "%zu-element subvector of %s."),
4888 iv->n, iv->n, lvalue->var->name);
4889 msg_at (SE, lsm, ngettext ("The source vector has %zu element.",
4890 "The source vector has %zu elements.",
4896 /* Assign elements. */
4897 gsl_vector dv = to_vector (lvalue->var->value);
4898 for (size_t x = 0; x < iv->n; x++)
4899 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4904 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4905 struct index_vector *iv0,
4906 struct index_vector *iv1,
4907 gsl_matrix *sm, const struct msg_location *lsm)
4909 gsl_matrix *dm = lvalue->var->value;
4911 /* Convert source matrix 'sm' to source vector 'sv'. */
4912 bool wrong_rows = iv0->n != sm->size1;
4913 bool wrong_columns = iv1->n != sm->size2;
4914 if (wrong_rows || wrong_columns)
4916 if (wrong_rows && wrong_columns)
4917 msg_at (SE, lvalue->full_location,
4918 _("Numbers of indexes for assigning to %s differ from the "
4919 "size of the source matrix."),
4921 else if (wrong_rows)
4922 msg_at (SE, lvalue->full_location,
4923 _("Number of row indexes for assigning to %s differs from "
4924 "number of rows in source matrix."),
4927 msg_at (SE, lvalue->full_location,
4928 _("Number of column indexes for assigning to %s differs from "
4929 "number of columns in source matrix."),
4934 if (lvalue->indexes[0])
4935 msg_at (SN, lvalue->index_locations[0],
4936 ngettext ("There is %zu row index.",
4937 "There are %zu row indexes.",
4941 msg_at (SN, lvalue->index_locations[0],
4942 ngettext ("Destination matrix %s has %zu row.",
4943 "Destination matrix %s has %zu rows.",
4945 lvalue->var->name, iv0->n);
4950 if (lvalue->indexes[1])
4951 msg_at (SN, lvalue->index_locations[1],
4952 ngettext ("There is %zu column index.",
4953 "There are %zu column indexes.",
4957 msg_at (SN, lvalue->index_locations[1],
4958 ngettext ("Destination matrix %s has %zu column.",
4959 "Destination matrix %s has %zu columns.",
4961 lvalue->var->name, iv1->n);
4964 msg_at (SN, lsm, _("The source matrix is %zu×%zu."),
4965 sm->size1, sm->size2);
4969 /* Assign elements. */
4970 for (size_t y = 0; y < iv0->n; y++)
4972 size_t f0 = iv0->indexes[y];
4973 for (size_t x = 0; x < iv1->n; x++)
4975 size_t f1 = iv1->indexes[x];
4976 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4982 /* Assigns rvalue SM to LVALUE, whose index vectors IV0 and IV1 were previously
4983 obtained with matrix_lvalue_evaluate(). Returns true if successful, false
4984 on error. Always takes ownership of M. LSM should be the source location
4985 to use for errors related to SM. */
4987 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4988 struct index_vector *iv0, struct index_vector *iv1,
4989 gsl_matrix *sm, const struct msg_location *lsm)
4991 if (!lvalue->n_indexes)
4993 gsl_matrix_free (lvalue->var->value);
4994 lvalue->var->value = sm;
4999 bool ok = (lvalue->n_indexes == 1
5000 ? matrix_lvalue_assign_vector (lvalue, iv0, sm, lsm)
5001 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm, lsm));
5002 index_vector_uninit (iv0);
5003 index_vector_uninit (iv1);
5004 gsl_matrix_free (sm);
5009 /* Evaluates and then assigns SM to LVALUE. Always takes ownership of M. LSM
5010 should be the source location to use for errors related to SM.*/
5012 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue,
5014 const struct msg_location *lsm)
5016 struct index_vector iv0, iv1;
5017 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
5019 gsl_matrix_free (sm);
5023 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm, lsm);
5026 /* Matrix command data structure. */
5028 /* An array of matrix commands. */
5029 struct matrix_commands
5031 struct matrix_command **commands;
5035 static bool matrix_commands_parse (struct matrix_state *,
5036 struct matrix_commands *,
5037 const char *command_name,
5038 const char *stop1, const char *stop2);
5039 static void matrix_commands_uninit (struct matrix_commands *);
5041 /* A single matrix command. */
5042 struct matrix_command
5044 /* The type of command. */
5045 enum matrix_command_type
5066 /* Source lines for this command. */
5067 struct msg_location *location;
5071 struct matrix_compute
5073 struct matrix_lvalue *lvalue;
5074 struct matrix_expr *rvalue;
5080 struct matrix_expr *expression;
5081 bool use_default_format;
5082 struct fmt_spec format;
5084 int space; /* -1 means NEWPAGE. */
5088 struct string_array literals; /* CLABELS/RLABELS. */
5089 struct matrix_expr *expr; /* CNAMES/RNAMES. */
5099 struct matrix_expr *condition;
5100 struct matrix_commands commands;
5110 struct matrix_var *var;
5111 struct matrix_expr *start;
5112 struct matrix_expr *end;
5113 struct matrix_expr *increment;
5115 /* Loop conditions. */
5116 struct matrix_expr *top_condition;
5117 struct matrix_expr *bottom_condition;
5120 struct matrix_commands commands;
5124 struct matrix_display
5126 struct matrix_state *state;
5130 struct matrix_release
5132 struct matrix_var **vars;
5139 struct matrix_expr *expression;
5140 struct save_file *sf;
5146 struct read_file *rf;
5147 struct matrix_lvalue *dst;
5148 struct matrix_expr *size;
5150 enum fmt_type format;
5159 struct write_file *wf;
5160 struct matrix_expr *expression;
5163 /* If this is nonnull, WRITE uses this format.
5165 If this is NULL, WRITE uses free-field format with as many
5166 digits of precision as needed. */
5167 struct fmt_spec *format;
5176 struct matrix_lvalue *dst;
5177 struct dataset *dataset;
5178 struct file_handle *file;
5180 struct var_syntax *vars;
5182 struct matrix_var *names;
5184 /* Treatment of missing values. */
5189 MGET_ERROR, /* Flag error on command. */
5190 MGET_ACCEPT, /* Accept without change, user-missing only. */
5191 MGET_OMIT, /* Drop this case. */
5192 MGET_RECODE /* Recode to 'substitute'. */
5195 double substitute; /* MGET_RECODE only. */
5203 struct msave_common *common;
5204 struct matrix_expr *expr;
5205 const char *rowtype;
5206 const struct matrix_expr *factors;
5207 const struct matrix_expr *splits;
5213 struct matrix_state *state;
5214 struct file_handle *file;
5216 struct stringi_set rowtypes;
5222 struct matrix_expr *expr;
5223 struct matrix_var *evec;
5224 struct matrix_var *eval;
5228 struct matrix_setdiag
5230 struct matrix_var *dst;
5231 struct matrix_expr *expr;
5237 struct matrix_expr *expr;
5238 struct matrix_var *u;
5239 struct matrix_var *s;
5240 struct matrix_var *v;
5246 static struct matrix_command *matrix_command_parse (struct matrix_state *);
5247 static bool matrix_command_execute (struct matrix_command *);
5248 static void matrix_command_destroy (struct matrix_command *);
5252 static struct matrix_command *
5253 matrix_compute_parse (struct matrix_state *s)
5255 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5256 *cmd = (struct matrix_command) {
5257 .type = MCMD_COMPUTE,
5258 .compute = { .lvalue = NULL },
5261 struct matrix_compute *compute = &cmd->compute;
5262 compute->lvalue = matrix_lvalue_parse (s);
5263 if (!compute->lvalue)
5266 if (!lex_force_match (s->lexer, T_EQUALS))
5269 compute->rvalue = matrix_expr_parse (s);
5270 if (!compute->rvalue)
5276 matrix_command_destroy (cmd);
5281 matrix_compute_execute (struct matrix_command *cmd)
5283 struct matrix_compute *compute = &cmd->compute;
5284 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
5288 matrix_lvalue_evaluate_and_assign (compute->lvalue, value,
5289 matrix_expr_location (compute->rvalue));
5295 print_labels_destroy (struct print_labels *labels)
5299 string_array_destroy (&labels->literals);
5300 matrix_expr_destroy (labels->expr);
5306 parse_literal_print_labels (struct matrix_state *s,
5307 struct print_labels **labelsp)
5309 lex_match (s->lexer, T_EQUALS);
5310 print_labels_destroy (*labelsp);
5311 *labelsp = xzalloc (sizeof **labelsp);
5312 while (lex_token (s->lexer) != T_SLASH
5313 && lex_token (s->lexer) != T_ENDCMD
5314 && lex_token (s->lexer) != T_STOP)
5316 struct string label = DS_EMPTY_INITIALIZER;
5317 while (lex_token (s->lexer) != T_COMMA
5318 && lex_token (s->lexer) != T_SLASH
5319 && lex_token (s->lexer) != T_ENDCMD
5320 && lex_token (s->lexer) != T_STOP)
5322 if (!ds_is_empty (&label))
5323 ds_put_byte (&label, ' ');
5325 if (lex_token (s->lexer) == T_STRING)
5326 ds_put_cstr (&label, lex_tokcstr (s->lexer));
5329 char *rep = lex_next_representation (s->lexer, 0, 0);
5330 ds_put_cstr (&label, rep);
5335 string_array_append_nocopy (&(*labelsp)->literals,
5336 ds_steal_cstr (&label));
5338 if (!lex_match (s->lexer, T_COMMA))
5344 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
5346 lex_match (s->lexer, T_EQUALS);
5347 struct matrix_expr *e = matrix_parse_exp (s);
5351 print_labels_destroy (*labelsp);
5352 *labelsp = xzalloc (sizeof **labelsp);
5353 (*labelsp)->expr = e;
5357 static struct matrix_command *
5358 matrix_print_parse (struct matrix_state *s)
5360 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5361 *cmd = (struct matrix_command) {
5364 .use_default_format = true,
5368 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
5371 for (size_t i = 0; ; i++)
5373 enum token_type t = lex_next_token (s->lexer, i);
5374 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
5376 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
5378 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
5381 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
5386 cmd->print.expression = matrix_parse_exp (s);
5387 if (!cmd->print.expression)
5391 while (lex_match (s->lexer, T_SLASH))
5393 if (lex_match_id (s->lexer, "FORMAT"))
5395 lex_match (s->lexer, T_EQUALS);
5396 if (!parse_format_specifier (s->lexer, &cmd->print.format))
5398 cmd->print.use_default_format = false;
5400 else if (lex_match_id (s->lexer, "TITLE"))
5402 lex_match (s->lexer, T_EQUALS);
5403 if (!lex_force_string (s->lexer))
5405 free (cmd->print.title);
5406 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
5409 else if (lex_match_id (s->lexer, "SPACE"))
5411 lex_match (s->lexer, T_EQUALS);
5412 if (lex_match_id (s->lexer, "NEWPAGE"))
5413 cmd->print.space = -1;
5414 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
5416 cmd->print.space = lex_integer (s->lexer);
5422 else if (lex_match_id (s->lexer, "RLABELS"))
5423 parse_literal_print_labels (s, &cmd->print.rlabels);
5424 else if (lex_match_id (s->lexer, "CLABELS"))
5425 parse_literal_print_labels (s, &cmd->print.clabels);
5426 else if (lex_match_id (s->lexer, "RNAMES"))
5428 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
5431 else if (lex_match_id (s->lexer, "CNAMES"))
5433 if (!parse_expr_print_labels (s, &cmd->print.clabels))
5438 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
5439 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
5447 matrix_command_destroy (cmd);
5452 matrix_is_integer (const gsl_matrix *m)
5454 for (size_t y = 0; y < m->size1; y++)
5455 for (size_t x = 0; x < m->size2; x++)
5457 double d = gsl_matrix_get (m, y, x);
5465 matrix_max_magnitude (const gsl_matrix *m)
5468 for (size_t y = 0; y < m->size1; y++)
5469 for (size_t x = 0; x < m->size2; x++)
5471 double d = fabs (gsl_matrix_get (m, y, x));
5479 format_fits (struct fmt_spec format, double x)
5481 char *s = data_out (&(union value) { .f = x }, NULL,
5482 &format, settings_get_fmt_settings ());
5483 bool fits = *s != '*' && !strchr (s, 'E');
5488 static struct fmt_spec
5489 default_format (const gsl_matrix *m, int *log_scale)
5493 double max = matrix_max_magnitude (m);
5495 if (matrix_is_integer (m))
5496 for (int w = 1; w <= 12; w++)
5498 struct fmt_spec format = { .type = FMT_F, .w = w };
5499 if (format_fits (format, -max))
5503 if (max >= 1e9 || max <= 1e-4)
5506 snprintf (s, sizeof s, "%.1e", max);
5508 const char *e = strchr (s, 'e');
5510 *log_scale = atoi (e + 1);
5513 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5517 trimmed_string (double d)
5519 struct substring s = ss_buffer ((char *) &d, sizeof d);
5520 ss_rtrim (&s, ss_cstr (" "));
5521 return ss_xstrdup (s);
5524 static struct string_array *
5525 print_labels_get (const struct print_labels *labels, size_t n,
5526 const char *prefix, bool truncate)
5531 struct string_array *out = xzalloc (sizeof *out);
5532 if (labels->literals.n)
5533 string_array_clone (out, &labels->literals);
5534 else if (labels->expr)
5536 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5537 if (m && is_vector (m))
5539 gsl_vector v = to_vector (m);
5540 for (size_t i = 0; i < v.size; i++)
5541 string_array_append_nocopy (out, trimmed_string (
5542 gsl_vector_get (&v, i)));
5544 gsl_matrix_free (m);
5550 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5553 string_array_append (out, "");
5556 string_array_delete (out, out->n - 1);
5559 for (size_t i = 0; i < out->n; i++)
5561 char *s = out->strings[i];
5562 s[strnlen (s, 8)] = '\0';
5569 matrix_print_space (int space)
5572 output_item_submit (page_break_item_create ());
5573 for (int i = 0; i < space; i++)
5574 output_log ("%s", "");
5578 matrix_print_text (const struct matrix_print *print, const gsl_matrix *m,
5579 struct fmt_spec format, int log_scale)
5581 matrix_print_space (print->space);
5583 output_log ("%s", print->title);
5585 output_log (" 10 ** %d X", log_scale);
5587 struct string_array *clabels = print_labels_get (print->clabels,
5588 m->size2, "col", true);
5589 if (clabels && format.w < 8)
5591 struct string_array *rlabels = print_labels_get (print->rlabels,
5592 m->size1, "row", true);
5596 struct string line = DS_EMPTY_INITIALIZER;
5598 ds_put_byte_multiple (&line, ' ', 8);
5599 for (size_t x = 0; x < m->size2; x++)
5600 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5601 output_log_nocopy (ds_steal_cstr (&line));
5604 double scale = pow (10.0, log_scale);
5605 bool numeric = fmt_is_numeric (format.type);
5606 for (size_t y = 0; y < m->size1; y++)
5608 struct string line = DS_EMPTY_INITIALIZER;
5610 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5612 for (size_t x = 0; x < m->size2; x++)
5614 double f = gsl_matrix_get (m, y, x);
5616 ? data_out (&(union value) { .f = f / scale}, NULL,
5617 &format, settings_get_fmt_settings ())
5618 : trimmed_string (f));
5619 ds_put_format (&line, " %s", s);
5622 output_log_nocopy (ds_steal_cstr (&line));
5625 string_array_destroy (rlabels);
5627 string_array_destroy (clabels);
5632 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5633 const struct print_labels *print_labels, size_t n,
5634 const char *name, const char *prefix)
5636 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5638 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5639 for (size_t i = 0; i < n; i++)
5640 pivot_category_create_leaf (
5642 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5643 : pivot_value_new_integer (i + 1)));
5645 d->hide_all_labels = true;
5646 string_array_destroy (labels);
5651 matrix_print_tables (const struct matrix_print *print, const gsl_matrix *m,
5652 struct fmt_spec format, int log_scale)
5654 struct pivot_table *table = pivot_table_create__ (
5655 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5658 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5660 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5661 N_("Columns"), "col");
5663 struct pivot_footnote *footnote = NULL;
5666 char *s = xasprintf ("× 10**%d", log_scale);
5667 footnote = pivot_table_create_footnote (
5668 table, pivot_value_new_user_text_nocopy (s));
5671 double scale = pow (10.0, log_scale);
5672 bool numeric = fmt_is_numeric (format.type);
5673 for (size_t y = 0; y < m->size1; y++)
5674 for (size_t x = 0; x < m->size2; x++)
5676 double f = gsl_matrix_get (m, y, x);
5677 struct pivot_value *v;
5680 v = pivot_value_new_number (f / scale);
5681 v->numeric.format = format;
5684 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5686 pivot_value_add_footnote (v, footnote);
5687 pivot_table_put2 (table, y, x, v);
5690 pivot_table_submit (table);
5694 matrix_print_execute (const struct matrix_print *print)
5696 if (print->expression)
5698 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5703 struct fmt_spec format = (print->use_default_format
5704 ? default_format (m, &log_scale)
5707 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5708 matrix_print_text (print, m, format, log_scale);
5710 matrix_print_tables (print, m, format, log_scale);
5712 gsl_matrix_free (m);
5716 matrix_print_space (print->space);
5718 output_log ("%s", print->title);
5725 matrix_do_if_clause_parse (struct matrix_state *s,
5726 struct matrix_do_if *ifc,
5727 bool parse_condition,
5728 size_t *allocated_clauses)
5730 if (ifc->n_clauses >= *allocated_clauses)
5731 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5732 sizeof *ifc->clauses);
5733 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5734 *c = (struct do_if_clause) { .condition = NULL };
5736 if (parse_condition)
5738 c->condition = matrix_expr_parse (s);
5743 return matrix_commands_parse (s, &c->commands, "DO IF", "ELSE", "END IF");
5746 static struct matrix_command *
5747 matrix_do_if_parse (struct matrix_state *s)
5749 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5750 *cmd = (struct matrix_command) {
5752 .do_if = { .n_clauses = 0 }
5755 size_t allocated_clauses = 0;
5758 if (!matrix_do_if_clause_parse (s, &cmd->do_if, true, &allocated_clauses))
5761 while (lex_match_phrase (s->lexer, "ELSE IF"));
5763 if (lex_match_id (s->lexer, "ELSE")
5764 && !matrix_do_if_clause_parse (s, &cmd->do_if, false, &allocated_clauses))
5767 if (!lex_match_phrase (s->lexer, "END IF"))
5772 matrix_command_destroy (cmd);
5777 matrix_do_if_execute (struct matrix_do_if *cmd)
5779 for (size_t i = 0; i < cmd->n_clauses; i++)
5781 struct do_if_clause *c = &cmd->clauses[i];
5785 if (!matrix_expr_evaluate_scalar (c->condition,
5786 i ? "ELSE IF" : "DO IF",
5791 for (size_t j = 0; j < c->commands.n; j++)
5792 if (!matrix_command_execute (c->commands.commands[j]))
5801 static struct matrix_command *
5802 matrix_loop_parse (struct matrix_state *s)
5804 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5805 *cmd = (struct matrix_command) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5807 struct matrix_loop *loop = &cmd->loop;
5808 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5810 struct substring name = lex_tokss (s->lexer);
5811 loop->var = matrix_var_lookup (s, name);
5813 loop->var = matrix_var_create (s, name);
5818 loop->start = matrix_expr_parse (s);
5819 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5822 loop->end = matrix_expr_parse (s);
5826 if (lex_match (s->lexer, T_BY))
5828 loop->increment = matrix_expr_parse (s);
5829 if (!loop->increment)
5834 if (lex_match_id (s->lexer, "IF"))
5836 loop->top_condition = matrix_expr_parse (s);
5837 if (!loop->top_condition)
5841 bool was_in_loop = s->in_loop;
5843 bool ok = matrix_commands_parse (s, &loop->commands, "LOOP",
5845 s->in_loop = was_in_loop;
5849 if (!lex_match_phrase (s->lexer, "END LOOP"))
5852 if (lex_match_id (s->lexer, "IF"))
5854 loop->bottom_condition = matrix_expr_parse (s);
5855 if (!loop->bottom_condition)
5862 matrix_command_destroy (cmd);
5867 matrix_loop_execute (struct matrix_loop *cmd)
5871 long int increment = 1;
5874 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5875 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5877 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5881 if (increment > 0 ? value > end
5882 : increment < 0 ? value < end
5887 int mxloops = settings_get_mxloops ();
5888 for (int i = 0; i < mxloops; i++)
5893 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5895 gsl_matrix_free (cmd->var->value);
5896 cmd->var->value = NULL;
5898 if (!cmd->var->value)
5899 cmd->var->value = gsl_matrix_alloc (1, 1);
5901 gsl_matrix_set (cmd->var->value, 0, 0, value);
5904 if (cmd->top_condition)
5907 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5912 for (size_t j = 0; j < cmd->commands.n; j++)
5913 if (!matrix_command_execute (cmd->commands.commands[j]))
5916 if (cmd->bottom_condition)
5919 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5920 "END LOOP IF", &d) || d > 0)
5927 if (increment > 0 ? value > end : value < end)
5935 static struct matrix_command *
5936 matrix_break_parse (struct matrix_state *s)
5940 msg (SE, _("BREAK not inside LOOP."));
5944 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5945 *cmd = (struct matrix_command) { .type = MCMD_BREAK };
5951 static struct matrix_command *
5952 matrix_display_parse (struct matrix_state *s)
5954 while (lex_token (s->lexer) != T_ENDCMD)
5956 if (!lex_match_id (s->lexer, "DICTIONARY")
5957 && !lex_match_id (s->lexer, "STATUS"))
5959 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5964 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5965 *cmd = (struct matrix_command) { .type = MCMD_DISPLAY, .display = { s } };
5970 compare_matrix_var_pointers (const void *a_, const void *b_)
5972 const struct matrix_var *const *ap = a_;
5973 const struct matrix_var *const *bp = b_;
5974 const struct matrix_var *a = *ap;
5975 const struct matrix_var *b = *bp;
5976 return strcmp (a->name, b->name);
5980 matrix_display_execute (struct matrix_display *cmd)
5982 const struct matrix_state *s = cmd->state;
5984 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5985 struct pivot_dimension *attr_dimension
5986 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5987 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5988 N_("Rows"), N_("Columns"));
5989 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5991 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5994 const struct matrix_var *var;
5995 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5997 vars[n_vars++] = var;
5998 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
6000 struct pivot_dimension *rows = pivot_dimension_create (
6001 table, PIVOT_AXIS_ROW, N_("Variable"));
6002 for (size_t i = 0; i < n_vars; i++)
6004 const struct matrix_var *var = vars[i];
6005 pivot_category_create_leaf (
6006 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
6008 size_t r = var->value->size1;
6009 size_t c = var->value->size2;
6010 double values[] = { r, c, r * c * sizeof (double) / 1024 };
6011 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
6012 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
6015 pivot_table_submit (table);
6020 static struct matrix_command *
6021 matrix_release_parse (struct matrix_state *s)
6023 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6024 *cmd = (struct matrix_command) { .type = MCMD_RELEASE };
6026 size_t allocated_vars = 0;
6027 while (lex_token (s->lexer) == T_ID)
6029 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
6032 if (cmd->release.n_vars >= allocated_vars)
6033 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
6034 sizeof *cmd->release.vars);
6035 cmd->release.vars[cmd->release.n_vars++] = var;
6038 lex_error (s->lexer, _("Variable name expected."));
6041 if (!lex_match (s->lexer, T_COMMA))
6049 matrix_release_execute (struct matrix_release *cmd)
6051 for (size_t i = 0; i < cmd->n_vars; i++)
6053 struct matrix_var *var = cmd->vars[i];
6054 gsl_matrix_free (var->value);
6061 static struct save_file *
6062 save_file_create (struct matrix_state *s, struct file_handle *fh,
6063 struct string_array *variables,
6064 struct matrix_expr *names,
6065 struct stringi_set *strings)
6067 for (size_t i = 0; i < s->n_save_files; i++)
6069 struct save_file *sf = s->save_files[i];
6070 if (fh_equal (sf->file, fh))
6074 string_array_destroy (variables);
6075 matrix_expr_destroy (names);
6076 stringi_set_destroy (strings);
6082 struct save_file *sf = xmalloc (sizeof *sf);
6083 *sf = (struct save_file) {
6085 .dataset = s->dataset,
6086 .variables = *variables,
6088 .strings = STRINGI_SET_INITIALIZER (sf->strings),
6091 stringi_set_swap (&sf->strings, strings);
6092 stringi_set_destroy (strings);
6094 s->save_files = xrealloc (s->save_files,
6095 (s->n_save_files + 1) * sizeof *s->save_files);
6096 s->save_files[s->n_save_files++] = sf;
6100 static struct casewriter *
6101 save_file_open (struct save_file *sf, gsl_matrix *m,
6102 const struct msg_location *save_location)
6104 if (sf->writer || sf->error)
6108 size_t n_variables = caseproto_get_n_widths (
6109 casewriter_get_proto (sf->writer));
6110 if (m->size2 != n_variables)
6112 const char *file_name = (sf->file == fh_inline_file () ? "*"
6113 : fh_get_name (sf->file));
6114 msg_at (SE, save_location,
6115 _("Cannot save %zu×%zu matrix to %s because the "
6116 "first SAVE to %s in this matrix program wrote a "
6117 "%zu-column matrix."),
6118 m->size1, m->size2, file_name, file_name, n_variables);
6119 msg_at (SE, sf->location,
6120 _("This is the location of the first SAVE to %s."),
6129 struct dictionary *dict = dict_create (get_default_encoding ());
6131 /* Fill 'names' with user-specified names if there were any, either from
6132 sf->variables or sf->names. */
6133 struct string_array names = { .n = 0 };
6134 if (sf->variables.n)
6135 string_array_clone (&names, &sf->variables);
6138 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
6139 if (nm && is_vector (nm))
6141 gsl_vector nv = to_vector (nm);
6142 for (size_t i = 0; i < nv.size; i++)
6144 char *name = trimmed_string (gsl_vector_get (&nv, i));
6145 if (dict_id_is_valid (dict, name, true))
6146 string_array_append_nocopy (&names, name);
6151 gsl_matrix_free (nm);
6154 struct stringi_set strings;
6155 stringi_set_clone (&strings, &sf->strings);
6157 for (size_t i = 0; dict_get_n_vars (dict) < m->size2; i++)
6162 name = names.strings[i];
6165 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
6169 int width = stringi_set_delete (&strings, name) ? 8 : 0;
6170 struct variable *var = dict_create_var (dict, name, width);
6173 msg_at (ME, save_location,
6174 _("Duplicate variable name %s in SAVE statement."), name);
6179 if (!stringi_set_is_empty (&strings))
6181 size_t count = stringi_set_count (&strings);
6182 const char *example = stringi_set_node_get_string (
6183 stringi_set_first (&strings));
6186 msg_at (ME, save_location,
6187 _("The SAVE command STRINGS subcommand specifies an unknown "
6188 "variable %s."), example);
6190 msg_at (ME, save_location,
6191 ngettext ("The SAVE command STRINGS subcommand specifies %zu "
6192 "unknown variable: %s.",
6193 "The SAVE command STRINGS subcommand specifies %zu "
6194 "unknown variables, including %s.",
6199 stringi_set_destroy (&strings);
6200 string_array_destroy (&names);
6209 if (sf->file == fh_inline_file ())
6210 sf->writer = autopaging_writer_create (dict_get_proto (dict));
6212 sf->writer = any_writer_open (sf->file, dict);
6216 sf->location = msg_location_dup (save_location);
6227 save_file_destroy (struct save_file *sf)
6231 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
6233 dataset_set_dict (sf->dataset, sf->dict);
6234 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
6238 casewriter_destroy (sf->writer);
6239 dict_unref (sf->dict);
6241 fh_unref (sf->file);
6242 string_array_destroy (&sf->variables);
6243 matrix_expr_destroy (sf->names);
6244 stringi_set_destroy (&sf->strings);
6245 msg_location_destroy (sf->location);
6250 static struct matrix_command *
6251 matrix_save_parse (struct matrix_state *s)
6253 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6254 *cmd = (struct matrix_command) {
6256 .save = { .expression = NULL },
6259 struct file_handle *fh = NULL;
6260 struct matrix_save *save = &cmd->save;
6262 struct string_array variables = STRING_ARRAY_INITIALIZER;
6263 struct matrix_expr *names = NULL;
6264 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
6266 save->expression = matrix_parse_exp (s);
6267 if (!save->expression)
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 = matrix_parse_exp (s);
6311 else if (lex_match_id (s->lexer, "STRINGS"))
6313 lex_match (s->lexer, T_EQUALS);
6314 while (lex_token (s->lexer) == T_ID)
6316 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
6318 if (!lex_match (s->lexer, T_COMMA))
6324 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
6332 if (s->prev_save_file)
6333 fh = fh_ref (s->prev_save_file);
6336 lex_sbc_missing ("OUTFILE");
6340 fh_unref (s->prev_save_file);
6341 s->prev_save_file = fh_ref (fh);
6343 if (variables.n && names)
6345 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
6346 matrix_expr_destroy (names);
6350 save->sf = save_file_create (s, fh, &variables, names, &strings);
6354 string_array_destroy (&variables);
6355 matrix_expr_destroy (names);
6356 stringi_set_destroy (&strings);
6358 matrix_command_destroy (cmd);
6363 matrix_save_execute (const struct matrix_command *cmd)
6365 const struct matrix_save *save = &cmd->save;
6367 gsl_matrix *m = matrix_expr_evaluate (save->expression);
6371 struct casewriter *writer = save_file_open (save->sf, m, cmd->location);
6374 gsl_matrix_free (m);
6378 const struct caseproto *proto = casewriter_get_proto (writer);
6379 for (size_t y = 0; y < m->size1; y++)
6381 struct ccase *c = case_create (proto);
6382 for (size_t x = 0; x < m->size2; x++)
6384 double d = gsl_matrix_get (m, y, x);
6385 int width = caseproto_get_width (proto, x);
6386 union value *value = case_data_rw_idx (c, x);
6390 memcpy (value->s, &d, width);
6392 casewriter_write (writer, c);
6394 gsl_matrix_free (m);
6399 static struct read_file *
6400 read_file_create (struct matrix_state *s, struct file_handle *fh)
6402 for (size_t i = 0; i < s->n_read_files; i++)
6404 struct read_file *rf = s->read_files[i];
6412 struct read_file *rf = xmalloc (sizeof *rf);
6413 *rf = (struct read_file) { .file = fh };
6415 s->read_files = xrealloc (s->read_files,
6416 (s->n_read_files + 1) * sizeof *s->read_files);
6417 s->read_files[s->n_read_files++] = rf;
6421 static struct dfm_reader *
6422 read_file_open (struct read_file *rf)
6425 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
6430 read_file_destroy (struct read_file *rf)
6434 fh_unref (rf->file);
6435 dfm_close_reader (rf->reader);
6436 free (rf->encoding);
6441 static struct matrix_command *
6442 matrix_read_parse (struct matrix_state *s)
6444 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6445 *cmd = (struct matrix_command) {
6447 .read = { .format = FMT_F },
6450 struct file_handle *fh = NULL;
6451 char *encoding = NULL;
6452 struct matrix_read *read = &cmd->read;
6453 read->dst = matrix_lvalue_parse (s);
6458 int repetitions = 0;
6459 int record_width = 0;
6460 bool seen_format = false;
6461 while (lex_match (s->lexer, T_SLASH))
6463 if (lex_match_id (s->lexer, "FILE"))
6465 lex_match (s->lexer, T_EQUALS);
6468 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6472 else if (lex_match_id (s->lexer, "ENCODING"))
6474 lex_match (s->lexer, T_EQUALS);
6475 if (!lex_force_string (s->lexer))
6479 encoding = ss_xstrdup (lex_tokss (s->lexer));
6483 else if (lex_match_id (s->lexer, "FIELD"))
6485 lex_match (s->lexer, T_EQUALS);
6487 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6489 read->c1 = lex_integer (s->lexer);
6491 if (!lex_force_match (s->lexer, T_TO)
6492 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
6494 read->c2 = lex_integer (s->lexer) + 1;
6497 record_width = read->c2 - read->c1;
6498 if (lex_match (s->lexer, T_BY))
6500 if (!lex_force_int_range (s->lexer, "BY", 1,
6501 read->c2 - read->c1))
6503 by = lex_integer (s->lexer);
6506 if (record_width % by)
6508 msg (SE, _("BY %d does not evenly divide record width %d."),
6516 else if (lex_match_id (s->lexer, "SIZE"))
6518 lex_match (s->lexer, T_EQUALS);
6519 matrix_expr_destroy (read->size);
6520 read->size = matrix_parse_exp (s);
6524 else if (lex_match_id (s->lexer, "MODE"))
6526 lex_match (s->lexer, T_EQUALS);
6527 if (lex_match_id (s->lexer, "RECTANGULAR"))
6528 read->symmetric = false;
6529 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6530 read->symmetric = true;
6533 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6537 else if (lex_match_id (s->lexer, "REREAD"))
6538 read->reread = true;
6539 else if (lex_match_id (s->lexer, "FORMAT"))
6543 lex_sbc_only_once ("FORMAT");
6548 lex_match (s->lexer, T_EQUALS);
6550 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6553 const char *p = lex_tokcstr (s->lexer);
6554 if (c_isdigit (p[0]))
6556 repetitions = atoi (p);
6557 p += strspn (p, "0123456789");
6558 if (!fmt_from_name (p, &read->format))
6560 lex_error (s->lexer, _("Unknown format %s."), p);
6565 else if (fmt_from_name (p, &read->format))
6569 struct fmt_spec format;
6570 if (!parse_format_specifier (s->lexer, &format))
6572 read->format = format.type;
6578 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6579 "REREAD", "FORMAT");
6586 lex_sbc_missing ("FIELD");
6590 if (!read->dst->n_indexes && !read->size)
6592 msg (SE, _("SIZE is required for reading data into a full matrix "
6593 "(as opposed to a submatrix)."));
6599 if (s->prev_read_file)
6600 fh = fh_ref (s->prev_read_file);
6603 lex_sbc_missing ("FILE");
6607 fh_unref (s->prev_read_file);
6608 s->prev_read_file = fh_ref (fh);
6610 read->rf = read_file_create (s, fh);
6614 free (read->rf->encoding);
6615 read->rf->encoding = encoding;
6619 /* Field width may be specified in multiple ways:
6622 2. The format on FORMAT.
6623 3. The repetition factor on FORMAT.
6625 (2) and (3) are mutually exclusive.
6627 If more than one of these is present, they must agree. If none of them is
6628 present, then free-field format is used.
6630 if (repetitions > record_width)
6632 msg (SE, _("%d repetitions cannot fit in record width %d."),
6633 repetitions, record_width);
6636 int w = (repetitions ? record_width / repetitions
6642 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6643 "which implies field width %d, "
6644 "but BY specifies field width %d."),
6645 repetitions, record_width, w, by);
6647 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6656 matrix_command_destroy (cmd);
6662 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6663 struct substring data, size_t y, size_t x,
6664 int first_column, int last_column, char *error)
6666 int line_number = dfm_get_line_number (reader);
6667 struct msg_location location = {
6668 .file_name = intern_new (dfm_get_file_name (reader)),
6669 .start = { .line = line_number, .column = first_column },
6670 .end = { .line = line_number, .column = last_column },
6672 msg_at (DW, &location, _("Error reading \"%.*s\" as format %s "
6673 "for matrix row %zu, column %zu: %s"),
6674 (int) data.length, data.string, fmt_name (format),
6675 y + 1, x + 1, error);
6676 msg_location_uninit (&location);
6681 matrix_read_set_field (struct matrix_read *read, struct dfm_reader *reader,
6682 gsl_matrix *m, struct substring p, size_t y, size_t x,
6683 const char *line_start)
6685 const char *input_encoding = dfm_reader_get_encoding (reader);
6688 if (fmt_is_numeric (read->format))
6691 error = data_in (p, input_encoding, read->format,
6692 settings_get_fmt_settings (), &v, 0, NULL);
6693 if (!error && v.f == SYSMIS)
6694 error = xstrdup (_("Matrix data may not contain missing value."));
6699 uint8_t s[sizeof (double)];
6700 union value v = { .s = s };
6701 error = data_in (p, input_encoding, read->format,
6702 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6703 memcpy (&f, s, sizeof f);
6708 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6709 int nc = ss_utf8_count_columns (p);
6710 int c2 = c1 + MAX (1, nc) - 1;
6711 parse_error (reader, read->format, p, y, x, c1, c2, error);
6715 gsl_matrix_set (m, y, x, f);
6716 if (read->symmetric && x != y)
6717 gsl_matrix_set (m, x, y, f);
6722 matrix_read_line (struct matrix_command *cmd, struct dfm_reader *reader,
6723 struct substring *line, const char **startp)
6725 struct matrix_read *read = &cmd->read;
6726 if (dfm_eof (reader))
6728 msg_at (SE, cmd->location,
6729 _("Unexpected end of file reading matrix data."));
6732 dfm_expand_tabs (reader);
6733 struct substring record = dfm_get_record (reader);
6734 /* XXX need to recode record into UTF-8 */
6735 *startp = record.string;
6736 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6741 matrix_read (struct matrix_command *cmd, struct dfm_reader *reader,
6744 struct matrix_read *read = &cmd->read;
6745 for (size_t y = 0; y < m->size1; y++)
6747 size_t nx = read->symmetric ? y + 1 : m->size2;
6749 struct substring line = ss_empty ();
6750 const char *line_start = line.string;
6751 for (size_t x = 0; x < nx; x++)
6758 ss_ltrim (&line, ss_cstr (" ,"));
6759 if (!ss_is_empty (line))
6761 if (!matrix_read_line (cmd, reader, &line, &line_start))
6763 dfm_forward_record (reader);
6766 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6770 if (!matrix_read_line (cmd, reader, &line, &line_start))
6772 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6773 int f = x % fields_per_line;
6774 if (f == fields_per_line - 1)
6775 dfm_forward_record (reader);
6777 p = ss_substr (line, read->w * f, read->w);
6780 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6784 dfm_forward_record (reader);
6787 ss_ltrim (&line, ss_cstr (" ,"));
6788 if (!ss_is_empty (line))
6790 int line_number = dfm_get_line_number (reader);
6791 int c1 = utf8_count_columns (line_start,
6792 line.string - line_start) + 1;
6793 int c2 = c1 + ss_utf8_count_columns (line) - 1;
6794 struct msg_location location = {
6795 .file_name = intern_new (dfm_get_file_name (reader)),
6796 .start = { .line = line_number, .column = c1 },
6797 .end = { .line = line_number, .column = c2 },
6799 msg_at (DW, &location,
6800 _("Trailing garbage following data for matrix row %zu."),
6802 msg_location_uninit (&location);
6809 matrix_read_execute (struct matrix_command *cmd)
6811 struct matrix_read *read = &cmd->read;
6812 struct index_vector iv0, iv1;
6813 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6816 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6819 gsl_matrix *m = matrix_expr_evaluate (read->size);
6825 msg_at (SE, matrix_expr_location (read->size),
6826 _("SIZE must evaluate to a scalar or a 2-element vector, "
6827 "not a %zu×%zu matrix."), m->size1, m->size2);
6828 gsl_matrix_free (m);
6829 index_vector_uninit (&iv0);
6830 index_vector_uninit (&iv1);
6834 gsl_vector v = to_vector (m);
6838 d[0] = gsl_vector_get (&v, 0);
6841 else if (v.size == 2)
6843 d[0] = gsl_vector_get (&v, 0);
6844 d[1] = gsl_vector_get (&v, 1);
6848 msg_at (SE, matrix_expr_location (read->size),
6849 _("SIZE must evaluate to a scalar or a 2-element vector, "
6850 "not a %zu×%zu matrix."),
6851 m->size1, m->size2),
6852 gsl_matrix_free (m);
6853 index_vector_uninit (&iv0);
6854 index_vector_uninit (&iv1);
6857 gsl_matrix_free (m);
6859 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6861 msg_at (SE, matrix_expr_location (read->size),
6862 _("Matrix dimensions %g×%g specified on SIZE "
6863 "are outside valid range."),
6865 index_vector_uninit (&iv0);
6866 index_vector_uninit (&iv1);
6874 if (read->dst->n_indexes)
6876 size_t submatrix_size[2];
6877 if (read->dst->n_indexes == 2)
6879 submatrix_size[0] = iv0.n;
6880 submatrix_size[1] = iv1.n;
6882 else if (read->dst->var->value->size1 == 1)
6884 submatrix_size[0] = 1;
6885 submatrix_size[1] = iv0.n;
6889 submatrix_size[0] = iv0.n;
6890 submatrix_size[1] = 1;
6895 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6897 msg_at (SE, cmd->location,
6898 _("Dimensions specified on SIZE differ from dimensions "
6899 "of destination submatrix."));
6900 msg_at (SN, matrix_expr_location (read->size),
6901 _("SIZE specifies dimensions %zu×%zu."),
6903 msg_at (SN, read->dst->full_location,
6904 _("Destination submatrix has dimensions %zu×%zu."),
6905 submatrix_size[0], submatrix_size[1]);
6906 index_vector_uninit (&iv0);
6907 index_vector_uninit (&iv1);
6913 size[0] = submatrix_size[0];
6914 size[1] = submatrix_size[1];
6918 struct dfm_reader *reader = read_file_open (read->rf);
6920 dfm_reread_record (reader, 1);
6922 if (read->symmetric && size[0] != size[1])
6924 msg_at (SE, cmd->location,
6925 _("Cannot read non-square %zu×%zu matrix "
6926 "using READ with MODE=SYMMETRIC."),
6928 index_vector_uninit (&iv0);
6929 index_vector_uninit (&iv1);
6932 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6933 matrix_read (cmd, reader, tmp);
6934 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp, cmd->location);
6939 static struct write_file *
6940 write_file_create (struct matrix_state *s, struct file_handle *fh)
6942 for (size_t i = 0; i < s->n_write_files; i++)
6944 struct write_file *wf = s->write_files[i];
6952 struct write_file *wf = xmalloc (sizeof *wf);
6953 *wf = (struct write_file) { .file = fh };
6955 s->write_files = xrealloc (s->write_files,
6956 (s->n_write_files + 1) * sizeof *s->write_files);
6957 s->write_files[s->n_write_files++] = wf;
6961 static struct dfm_writer *
6962 write_file_open (struct write_file *wf)
6965 wf->writer = dfm_open_writer (wf->file, wf->encoding);
6970 write_file_destroy (struct write_file *wf)
6976 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
6977 wf->held->s.ss.length);
6978 u8_line_destroy (wf->held);
6982 fh_unref (wf->file);
6983 dfm_close_writer (wf->writer);
6984 free (wf->encoding);
6989 static struct matrix_command *
6990 matrix_write_parse (struct matrix_state *s)
6992 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6993 *cmd = (struct matrix_command) {
6997 struct file_handle *fh = NULL;
6998 char *encoding = NULL;
6999 struct matrix_write *write = &cmd->write;
7000 write->expression = matrix_parse_exp (s);
7001 if (!write->expression)
7005 int repetitions = 0;
7006 int record_width = 0;
7007 enum fmt_type format = FMT_F;
7008 bool has_format = false;
7009 while (lex_match (s->lexer, T_SLASH))
7011 if (lex_match_id (s->lexer, "OUTFILE"))
7013 lex_match (s->lexer, T_EQUALS);
7016 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
7020 else if (lex_match_id (s->lexer, "ENCODING"))
7022 lex_match (s->lexer, T_EQUALS);
7023 if (!lex_force_string (s->lexer))
7027 encoding = ss_xstrdup (lex_tokss (s->lexer));
7031 else if (lex_match_id (s->lexer, "FIELD"))
7033 lex_match (s->lexer, T_EQUALS);
7035 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
7037 write->c1 = lex_integer (s->lexer);
7039 if (!lex_force_match (s->lexer, T_TO)
7040 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
7042 write->c2 = lex_integer (s->lexer) + 1;
7045 record_width = write->c2 - write->c1;
7046 if (lex_match (s->lexer, T_BY))
7048 if (!lex_force_int_range (s->lexer, "BY", 1,
7049 write->c2 - write->c1))
7051 by = lex_integer (s->lexer);
7054 if (record_width % by)
7056 msg (SE, _("BY %d does not evenly divide record width %d."),
7064 else if (lex_match_id (s->lexer, "MODE"))
7066 lex_match (s->lexer, T_EQUALS);
7067 if (lex_match_id (s->lexer, "RECTANGULAR"))
7068 write->triangular = false;
7069 else if (lex_match_id (s->lexer, "TRIANGULAR"))
7070 write->triangular = true;
7073 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
7077 else if (lex_match_id (s->lexer, "HOLD"))
7079 else if (lex_match_id (s->lexer, "FORMAT"))
7081 if (has_format || write->format)
7083 lex_sbc_only_once ("FORMAT");
7087 lex_match (s->lexer, T_EQUALS);
7089 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
7092 const char *p = lex_tokcstr (s->lexer);
7093 if (c_isdigit (p[0]))
7095 repetitions = atoi (p);
7096 p += strspn (p, "0123456789");
7097 if (!fmt_from_name (p, &format))
7099 lex_error (s->lexer, _("Unknown format %s."), p);
7105 else if (fmt_from_name (p, &format))
7112 struct fmt_spec spec;
7113 if (!parse_format_specifier (s->lexer, &spec))
7115 write->format = xmemdup (&spec, sizeof spec);
7120 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
7128 lex_sbc_missing ("FIELD");
7134 if (s->prev_write_file)
7135 fh = fh_ref (s->prev_write_file);
7138 lex_sbc_missing ("OUTFILE");
7142 fh_unref (s->prev_write_file);
7143 s->prev_write_file = fh_ref (fh);
7145 write->wf = write_file_create (s, fh);
7149 free (write->wf->encoding);
7150 write->wf->encoding = encoding;
7154 /* Field width may be specified in multiple ways:
7157 2. The format on FORMAT.
7158 3. The repetition factor on FORMAT.
7160 (2) and (3) are mutually exclusive.
7162 If more than one of these is present, they must agree. If none of them is
7163 present, then free-field format is used.
7165 if (repetitions > record_width)
7167 msg (SE, _("%d repetitions cannot fit in record width %d."),
7168 repetitions, record_width);
7171 int w = (repetitions ? record_width / repetitions
7172 : write->format ? write->format->w
7177 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
7178 "which implies field width %d, "
7179 "but BY specifies field width %d."),
7180 repetitions, record_width, w, by);
7182 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
7186 if (w && !write->format)
7188 write->format = xmalloc (sizeof *write->format);
7189 *write->format = (struct fmt_spec) { .type = format, .w = w };
7191 if (!fmt_check_output (write->format))
7195 if (write->format && fmt_var_width (write->format) > sizeof (double))
7197 char s[FMT_STRING_LEN_MAX + 1];
7198 fmt_to_string (write->format, s);
7199 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
7200 s, sizeof (double));
7208 matrix_command_destroy (cmd);
7213 matrix_write_execute (struct matrix_write *write)
7215 gsl_matrix *m = matrix_expr_evaluate (write->expression);
7219 if (write->triangular && m->size1 != m->size2)
7221 msg_at (SE, matrix_expr_location (write->expression),
7222 _("WRITE with MODE=TRIANGULAR requires a square matrix but "
7223 "the matrix to be written has dimensions %zu×%zu."),
7224 m->size1, m->size2);
7225 gsl_matrix_free (m);
7229 struct dfm_writer *writer = write_file_open (write->wf);
7230 if (!writer || !m->size1)
7232 gsl_matrix_free (m);
7236 const struct fmt_settings *settings = settings_get_fmt_settings ();
7237 struct u8_line *line = write->wf->held;
7238 for (size_t y = 0; y < m->size1; y++)
7242 line = xmalloc (sizeof *line);
7243 u8_line_init (line);
7245 size_t nx = write->triangular ? y + 1 : m->size2;
7247 for (size_t x = 0; x < nx; x++)
7250 double f = gsl_matrix_get (m, y, x);
7254 if (fmt_is_numeric (write->format->type))
7257 v.s = (uint8_t *) &f;
7258 s = data_out (&v, NULL, write->format, settings);
7262 s = xmalloc (DBL_BUFSIZE_BOUND);
7263 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
7264 >= DBL_BUFSIZE_BOUND)
7267 size_t len = strlen (s);
7268 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
7269 if (width + x0 > write->c2)
7271 dfm_put_record_utf8 (writer, line->s.ss.string,
7273 u8_line_clear (line);
7276 u8_line_put (line, x0, x0 + width, s, len);
7279 x0 += write->format ? write->format->w : width + 1;
7282 if (y + 1 >= m->size1 && write->hold)
7284 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
7285 u8_line_clear (line);
7289 u8_line_destroy (line);
7293 write->wf->held = line;
7295 gsl_matrix_free (m);
7300 static struct matrix_command *
7301 matrix_get_parse (struct matrix_state *s)
7303 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7304 *cmd = (struct matrix_command) {
7307 .dataset = s->dataset,
7308 .user = { .treatment = MGET_ERROR },
7309 .system = { .treatment = MGET_ERROR },
7313 struct matrix_get *get = &cmd->get;
7314 get->dst = matrix_lvalue_parse (s);
7318 while (lex_match (s->lexer, T_SLASH))
7320 if (lex_match_id (s->lexer, "FILE"))
7322 lex_match (s->lexer, T_EQUALS);
7324 fh_unref (get->file);
7325 if (lex_match (s->lexer, T_ASTERISK))
7329 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7334 else if (lex_match_id (s->lexer, "ENCODING"))
7336 lex_match (s->lexer, T_EQUALS);
7337 if (!lex_force_string (s->lexer))
7340 free (get->encoding);
7341 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
7345 else if (lex_match_id (s->lexer, "VARIABLES"))
7347 lex_match (s->lexer, T_EQUALS);
7351 lex_sbc_only_once ("VARIABLES");
7355 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
7358 else if (lex_match_id (s->lexer, "NAMES"))
7360 lex_match (s->lexer, T_EQUALS);
7361 if (!lex_force_id (s->lexer))
7364 struct substring name = lex_tokss (s->lexer);
7365 get->names = matrix_var_lookup (s, name);
7367 get->names = matrix_var_create (s, name);
7370 else if (lex_match_id (s->lexer, "MISSING"))
7372 lex_match (s->lexer, T_EQUALS);
7373 if (lex_match_id (s->lexer, "ACCEPT"))
7374 get->user.treatment = MGET_ACCEPT;
7375 else if (lex_match_id (s->lexer, "OMIT"))
7376 get->user.treatment = MGET_OMIT;
7377 else if (lex_is_number (s->lexer))
7379 get->user.treatment = MGET_RECODE;
7380 get->user.substitute = lex_number (s->lexer);
7385 lex_error (s->lexer, NULL);
7389 else if (lex_match_id (s->lexer, "SYSMIS"))
7391 lex_match (s->lexer, T_EQUALS);
7392 if (lex_match_id (s->lexer, "OMIT"))
7393 get->system.treatment = MGET_OMIT;
7394 else if (lex_is_number (s->lexer))
7396 get->system.treatment = MGET_RECODE;
7397 get->system.substitute = lex_number (s->lexer);
7402 lex_error (s->lexer, NULL);
7408 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
7409 "MISSING", "SYSMIS");
7414 if (get->user.treatment != MGET_ACCEPT)
7415 get->system.treatment = MGET_ERROR;
7420 matrix_command_destroy (cmd);
7425 matrix_get_execute__ (struct matrix_command *cmd, struct casereader *reader,
7426 const struct dictionary *dict)
7428 struct matrix_get *get = &cmd->get;
7429 struct variable **vars;
7434 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
7435 &vars, &n_vars, PV_NUMERIC))
7440 n_vars = dict_get_n_vars (dict);
7441 vars = xnmalloc (n_vars, sizeof *vars);
7442 for (size_t i = 0; i < n_vars; i++)
7444 struct variable *var = dict_get_var (dict, i);
7445 if (!var_is_numeric (var))
7447 msg_at (SE, cmd->location, _("Variable %s is not numeric."),
7448 var_get_name (var));
7458 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
7459 for (size_t i = 0; i < n_vars; i++)
7461 char s[sizeof (double)];
7463 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
7464 memcpy (&f, s, sizeof f);
7465 gsl_matrix_set (names, i, 0, f);
7468 gsl_matrix_free (get->names->value);
7469 get->names->value = names;
7473 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
7474 long long int casenum = 1;
7476 for (struct ccase *c = casereader_read (reader); c;
7477 c = casereader_read (reader), casenum++)
7479 if (n_rows >= m->size1)
7481 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
7482 for (size_t y = 0; y < n_rows; y++)
7483 for (size_t x = 0; x < n_vars; x++)
7484 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
7485 gsl_matrix_free (m);
7490 for (size_t x = 0; x < n_vars; x++)
7492 const struct variable *var = vars[x];
7493 double d = case_num (c, var);
7496 if (get->system.treatment == MGET_RECODE)
7497 d = get->system.substitute;
7498 else if (get->system.treatment == MGET_OMIT)
7502 msg_at (SE, cmd->location, _("Variable %s in case %lld "
7503 "is system-missing."),
7504 var_get_name (var), casenum);
7508 else if (var_is_num_missing (var, d) == MV_USER)
7510 if (get->user.treatment == MGET_RECODE)
7511 d = get->user.substitute;
7512 else if (get->user.treatment == MGET_OMIT)
7514 else if (get->user.treatment != MGET_ACCEPT)
7516 msg_at (SE, cmd->location,
7517 _("Variable %s in case %lld has user-missing "
7519 var_get_name (var), casenum, d);
7523 gsl_matrix_set (m, n_rows, x, d);
7534 matrix_lvalue_evaluate_and_assign (get->dst, m, cmd->location);
7537 gsl_matrix_free (m);
7542 matrix_open_casereader (const struct matrix_command *cmd,
7543 const char *command_name, struct file_handle *file,
7544 const char *encoding, struct dataset *dataset,
7545 struct casereader **readerp, struct dictionary **dictp)
7549 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
7550 return *readerp != NULL;
7554 if (dict_get_n_vars (dataset_dict (dataset)) == 0)
7556 msg_at (ME, cmd->location,
7557 _("The %s command cannot read an empty active file."),
7561 *readerp = proc_open (dataset);
7562 *dictp = dict_ref (dataset_dict (dataset));
7568 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
7569 struct casereader *reader, struct dictionary *dict)
7572 casereader_destroy (reader);
7574 proc_commit (dataset);
7578 matrix_get_execute (struct matrix_command *cmd)
7580 struct matrix_get *get = &cmd->get;
7581 struct casereader *r;
7582 struct dictionary *d;
7583 if (matrix_open_casereader (cmd, "GET", get->file, get->encoding,
7584 get->dataset, &r, &d))
7586 matrix_get_execute__ (cmd, r, d);
7587 matrix_close_casereader (get->file, get->dataset, r, d);
7594 variables_changed (const char *keyword,
7595 const struct string_array *new,
7596 const struct string_array *old)
7602 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7603 "on the first MSAVE within MATRIX."), keyword);
7606 else if (!string_array_equal_case (old, new))
7608 msg (SE, _("%s must specify the same variables each time within "
7609 "a given MATRIX."), keyword);
7617 msave_common_changed (const struct msave_common *old,
7618 const struct msave_common *new)
7620 if (new->outfile && !fh_equal (old->outfile, new->outfile))
7621 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7622 "within a single MATRIX command."));
7623 else if (variables_changed ("VARIABLES", &new->variables, &old->variables)
7624 || variables_changed ("FNAMES", &new->fnames, &old->fnames)
7625 || variables_changed ("SNAMES", &new->snames, &old->snames))
7626 msg_at (SN, old->location,
7627 _("This is the location of the first MSAVE command."));
7635 msave_common_destroy (struct msave_common *common)
7639 msg_location_destroy (common->location);
7640 fh_unref (common->outfile);
7641 string_array_destroy (&common->variables);
7642 string_array_destroy (&common->fnames);
7643 string_array_destroy (&common->snames);
7645 for (size_t i = 0; i < common->n_factors; i++)
7646 matrix_expr_destroy (common->factors[i]);
7647 free (common->factors);
7649 for (size_t i = 0; i < common->n_splits; i++)
7650 matrix_expr_destroy (common->splits[i]);
7651 free (common->splits);
7653 dict_unref (common->dict);
7654 casewriter_destroy (common->writer);
7661 match_rowtype (struct lexer *lexer)
7663 static const char *rowtypes[] = {
7664 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7666 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7668 for (size_t i = 0; i < n_rowtypes; i++)
7669 if (lex_match_id (lexer, rowtypes[i]))
7672 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7677 parse_var_names (struct lexer *lexer, struct string_array *sa)
7679 lex_match (lexer, T_EQUALS);
7681 string_array_clear (sa);
7683 struct dictionary *dict = dict_create (get_default_encoding ());
7686 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7687 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7692 for (size_t i = 0; i < n_names; i++)
7693 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7694 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7696 msg (SE, _("Variable name %s is reserved."), names[i]);
7697 for (size_t j = 0; j < n_names; j++)
7703 string_array_clear (sa);
7704 sa->strings = names;
7705 sa->n = sa->allocated = n_names;
7710 static struct matrix_command *
7711 matrix_msave_parse (struct matrix_state *s)
7713 int start_ofs = lex_ofs (s->lexer);
7715 struct msave_common *common = xmalloc (sizeof *common);
7716 *common = (struct msave_common) { .outfile = NULL };
7718 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7719 *cmd = (struct matrix_command) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7721 struct matrix_expr *splits = NULL;
7722 struct matrix_expr *factors = NULL;
7724 struct matrix_msave *msave = &cmd->msave;
7725 msave->expr = matrix_parse_exp (s);
7729 while (lex_match (s->lexer, T_SLASH))
7731 if (lex_match_id (s->lexer, "TYPE"))
7733 lex_match (s->lexer, T_EQUALS);
7735 msave->rowtype = match_rowtype (s->lexer);
7736 if (!msave->rowtype)
7739 else if (lex_match_id (s->lexer, "OUTFILE"))
7741 lex_match (s->lexer, T_EQUALS);
7743 fh_unref (common->outfile);
7744 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7745 if (!common->outfile)
7748 else if (lex_match_id (s->lexer, "VARIABLES"))
7750 if (!parse_var_names (s->lexer, &common->variables))
7753 else if (lex_match_id (s->lexer, "FNAMES"))
7755 if (!parse_var_names (s->lexer, &common->fnames))
7758 else if (lex_match_id (s->lexer, "SNAMES"))
7760 if (!parse_var_names (s->lexer, &common->snames))
7763 else if (lex_match_id (s->lexer, "SPLIT"))
7765 lex_match (s->lexer, T_EQUALS);
7767 matrix_expr_destroy (splits);
7768 splits = matrix_parse_exp (s);
7772 else if (lex_match_id (s->lexer, "FACTOR"))
7774 lex_match (s->lexer, T_EQUALS);
7776 matrix_expr_destroy (factors);
7777 factors = matrix_parse_exp (s);
7783 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7784 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7788 if (!msave->rowtype)
7790 lex_sbc_missing ("TYPE");
7794 if (!s->msave_common)
7796 if (common->fnames.n && !factors)
7798 msg (SE, _("FNAMES requires FACTOR."));
7801 if (common->snames.n && !splits)
7803 msg (SE, _("SNAMES requires SPLIT."));
7806 if (!common->outfile)
7808 lex_sbc_missing ("OUTFILE");
7811 common->location = lex_ofs_location (s->lexer, start_ofs,
7812 lex_ofs (s->lexer));
7813 msg_location_remove_columns (common->location);
7814 s->msave_common = common;
7818 if (msave_common_changed (s->msave_common, common))
7820 msave_common_destroy (common);
7822 msave->common = s->msave_common;
7824 struct msave_common *c = s->msave_common;
7827 if (c->n_factors >= c->allocated_factors)
7828 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7829 sizeof *c->factors);
7830 c->factors[c->n_factors++] = factors;
7832 if (c->n_factors > 0)
7833 msave->factors = c->factors[c->n_factors - 1];
7837 if (c->n_splits >= c->allocated_splits)
7838 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7840 c->splits[c->n_splits++] = splits;
7842 if (c->n_splits > 0)
7843 msave->splits = c->splits[c->n_splits - 1];
7848 matrix_expr_destroy (splits);
7849 matrix_expr_destroy (factors);
7850 msave_common_destroy (common);
7851 matrix_command_destroy (cmd);
7856 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7858 gsl_matrix *m = matrix_expr_evaluate (e);
7864 msg_at (SE, matrix_expr_location (e),
7865 _("%s expression must evaluate to vector, "
7866 "not a %zu×%zu matrix."),
7867 name, m->size1, m->size2);
7868 gsl_matrix_free (m);
7872 return matrix_to_vector (m);
7876 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7878 for (size_t i = 0; i < vars->n; i++)
7879 if (!dict_create_var (d, vars->strings[i], 0))
7880 return vars->strings[i];
7884 static struct dictionary *
7885 msave_create_dict (const struct msave_common *common,
7886 const struct msg_location *location)
7888 struct dictionary *dict = dict_create (get_default_encoding ());
7890 const char *dup_split = msave_add_vars (dict, &common->snames);
7893 /* Should not be possible because the parser ensures that the names are
7898 dict_create_var_assert (dict, "ROWTYPE_", 8);
7900 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7903 msg_at (SE, location, _("Duplicate or invalid FACTOR variable name %s."),
7908 dict_create_var_assert (dict, "VARNAME_", 8);
7910 const char *dup_var = msave_add_vars (dict, &common->variables);
7913 msg_at (SE, location, _("Duplicate or invalid variable name %s."),
7926 matrix_msave_execute (struct matrix_command *cmd)
7928 struct matrix_msave *msave = &cmd->msave;
7929 struct msave_common *common = msave->common;
7930 gsl_matrix *m = NULL;
7931 gsl_vector *factors = NULL;
7932 gsl_vector *splits = NULL;
7934 m = matrix_expr_evaluate (msave->expr);
7938 if (!common->variables.n)
7939 for (size_t i = 0; i < m->size2; i++)
7940 string_array_append_nocopy (&common->variables,
7941 xasprintf ("COL%zu", i + 1));
7942 else if (m->size2 != common->variables.n)
7944 msg_at (SE, matrix_expr_location (msave->expr),
7945 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7946 m->size2, common->variables.n);
7952 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7956 if (!common->fnames.n)
7957 for (size_t i = 0; i < factors->size; i++)
7958 string_array_append_nocopy (&common->fnames,
7959 xasprintf ("FAC%zu", i + 1));
7960 else if (factors->size != common->fnames.n)
7962 msg_at (SE, matrix_expr_location (msave->factors),
7963 _("There are %zu factor variables, "
7964 "but %zu factor values were supplied."),
7965 common->fnames.n, factors->size);
7972 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7976 if (!common->snames.n)
7977 for (size_t i = 0; i < splits->size; i++)
7978 string_array_append_nocopy (&common->snames,
7979 xasprintf ("SPL%zu", i + 1));
7980 else if (splits->size != common->snames.n)
7982 msg_at (SE, matrix_expr_location (msave->splits),
7983 _("There are %zu split variables, "
7984 "but %zu split values were supplied."),
7985 common->snames.n, splits->size);
7990 if (!common->writer)
7992 struct dictionary *dict = msave_create_dict (common, cmd->location);
7996 common->writer = any_writer_open (common->outfile, dict);
7997 if (!common->writer)
8003 common->dict = dict;
8006 bool matrix = (!strcmp (msave->rowtype, "COV")
8007 || !strcmp (msave->rowtype, "CORR"));
8008 for (size_t y = 0; y < m->size1; y++)
8010 struct ccase *c = case_create (dict_get_proto (common->dict));
8013 /* Split variables */
8015 for (size_t i = 0; i < splits->size; i++)
8016 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
8019 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8020 msave->rowtype, ' ');
8024 for (size_t i = 0; i < factors->size; i++)
8025 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
8028 const char *varname_ = (matrix && y < common->variables.n
8029 ? common->variables.strings[y]
8031 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8034 /* Continuous variables. */
8035 for (size_t x = 0; x < m->size2; x++)
8036 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
8037 casewriter_write (common->writer, c);
8041 gsl_matrix_free (m);
8042 gsl_vector_free (factors);
8043 gsl_vector_free (splits);
8048 static struct matrix_command *
8049 matrix_mget_parse (struct matrix_state *s)
8051 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8052 *cmd = (struct matrix_command) {
8056 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
8060 struct matrix_mget *mget = &cmd->mget;
8062 lex_match (s->lexer, T_SLASH);
8063 while (lex_token (s->lexer) != T_ENDCMD)
8065 if (lex_match_id (s->lexer, "FILE"))
8067 lex_match (s->lexer, T_EQUALS);
8069 fh_unref (mget->file);
8070 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
8074 else if (lex_match_id (s->lexer, "ENCODING"))
8076 lex_match (s->lexer, T_EQUALS);
8077 if (!lex_force_string (s->lexer))
8080 free (mget->encoding);
8081 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
8085 else if (lex_match_id (s->lexer, "TYPE"))
8087 lex_match (s->lexer, T_EQUALS);
8088 while (lex_token (s->lexer) != T_SLASH
8089 && lex_token (s->lexer) != T_ENDCMD)
8091 const char *rowtype = match_rowtype (s->lexer);
8095 stringi_set_insert (&mget->rowtypes, rowtype);
8100 lex_error_expecting (s->lexer, "FILE", "TYPE");
8103 lex_match (s->lexer, T_SLASH);
8108 matrix_command_destroy (cmd);
8112 static const struct variable *
8113 get_a8_var (const struct msg_location *loc,
8114 const struct dictionary *d, const char *name)
8116 const struct variable *v = dict_lookup_var (d, name);
8119 msg_at (SE, loc, _("Matrix data file lacks %s variable."), name);
8122 if (var_get_width (v) != 8)
8124 msg_at (SE, loc, _("%s variable in matrix data file must be "
8125 "8-byte string, but it has width %d."),
8126 name, var_get_width (v));
8133 var_changed (const struct ccase *ca, const struct ccase *cb,
8134 const struct variable *var)
8137 ? !value_equal (case_data (ca, var), case_data (cb, var),
8138 var_get_width (var))
8143 vars_changed (const struct ccase *ca, const struct ccase *cb,
8144 const struct dictionary *d,
8145 size_t first_var, size_t n_vars)
8147 for (size_t i = 0; i < n_vars; i++)
8149 const struct variable *v = dict_get_var (d, first_var + i);
8150 if (var_changed (ca, cb, v))
8157 vars_all_missing (const struct ccase *c, const struct dictionary *d,
8158 size_t first_var, size_t n_vars)
8160 for (size_t i = 0; i < n_vars; i++)
8161 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
8167 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
8168 const struct dictionary *d,
8169 const struct variable *rowtype_var,
8170 const struct stringi_set *accepted_rowtypes,
8171 struct matrix_state *s,
8172 size_t ss, size_t sn, size_t si,
8173 size_t fs, size_t fn, size_t fi,
8174 size_t cs, size_t cn,
8175 struct pivot_table *pt,
8176 struct pivot_dimension *var_dimension)
8181 /* Is this a matrix for pooled data, either where there are no factor
8182 variables or the factor variables are missing? */
8183 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
8185 struct substring rowtype = case_ss (rows[0], rowtype_var);
8186 ss_rtrim (&rowtype, ss_cstr (" "));
8187 if (!stringi_set_is_empty (accepted_rowtypes)
8188 && !stringi_set_contains_len (accepted_rowtypes,
8189 rowtype.string, rowtype.length))
8192 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
8193 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
8194 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
8195 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
8196 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
8197 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
8201 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
8202 (int) rowtype.length, rowtype.string);
8206 struct string name = DS_EMPTY_INITIALIZER;
8207 ds_put_cstr (&name, prefix);
8209 ds_put_format (&name, "F%zu", fi);
8211 ds_put_format (&name, "S%zu", si);
8213 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
8215 mv = matrix_var_create (s, ds_ss (&name));
8218 msg (SW, _("Matrix data file contains variable with existing name %s."),
8220 goto exit_free_name;
8223 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
8224 size_t n_missing = 0;
8225 for (size_t y = 0; y < n_rows; y++)
8227 for (size_t x = 0; x < cn; x++)
8229 struct variable *var = dict_get_var (d, cs + x);
8230 double value = case_num (rows[y], var);
8231 if (var_is_num_missing (var, value))
8236 gsl_matrix_set (m, y, x, value);
8240 int var_index = pivot_category_create_leaf (
8241 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
8242 double values[] = { n_rows, cn };
8243 for (size_t j = 0; j < sn; j++)
8245 struct variable *var = dict_get_var (d, ss + j);
8246 const union value *value = case_data (rows[0], var);
8247 pivot_table_put2 (pt, j, var_index,
8248 pivot_value_new_var_value (var, value));
8250 for (size_t j = 0; j < fn; j++)
8252 struct variable *var = dict_get_var (d, fs + j);
8253 const union value sysmis = { .f = SYSMIS };
8254 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
8255 pivot_table_put2 (pt, j + sn, var_index,
8256 pivot_value_new_var_value (var, value));
8258 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
8259 pivot_table_put2 (pt, j + sn + fn, var_index,
8260 pivot_value_new_integer (values[j]));
8263 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
8264 "value, which was treated as zero.",
8265 "Matrix data file variable %s contains %zu missing "
8266 "values, which were treated as zero.", n_missing),
8267 ds_cstr (&name), n_missing);
8274 for (size_t y = 0; y < n_rows; y++)
8275 case_unref (rows[y]);
8279 matrix_mget_execute__ (struct matrix_command *cmd, struct casereader *r,
8280 const struct dictionary *d)
8282 struct matrix_mget *mget = &cmd->mget;
8283 const struct msg_location *loc = cmd->location;
8284 const struct variable *rowtype_ = get_a8_var (loc, d, "ROWTYPE_");
8285 const struct variable *varname_ = get_a8_var (loc, d, "VARNAME_");
8286 if (!rowtype_ || !varname_)
8289 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
8292 _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
8295 if (var_get_dict_index (varname_) + 1 >= dict_get_n_vars (d))
8297 msg_at (SE, loc, _("Matrix data file contains no continuous variables."));
8301 for (size_t i = 0; i < dict_get_n_vars (d); i++)
8303 const struct variable *v = dict_get_var (d, i);
8304 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
8307 _("Matrix data file contains unexpected string variable %s."),
8313 /* SPLIT variables. */
8315 size_t sn = var_get_dict_index (rowtype_);
8316 struct ccase *sc = NULL;
8319 /* FACTOR variables. */
8320 size_t fs = var_get_dict_index (rowtype_) + 1;
8321 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
8322 struct ccase *fc = NULL;
8325 /* Continuous variables. */
8326 size_t cs = var_get_dict_index (varname_) + 1;
8327 size_t cn = dict_get_n_vars (d) - cs;
8328 struct ccase *cc = NULL;
8331 struct pivot_table *pt = pivot_table_create (
8332 N_("Matrix Variables Created by MGET"));
8333 struct pivot_dimension *attr_dimension = pivot_dimension_create (
8334 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
8335 struct pivot_dimension *var_dimension = pivot_dimension_create (
8336 pt, PIVOT_AXIS_ROW, N_("Variable"));
8339 struct pivot_category *splits = pivot_category_create_group (
8340 attr_dimension->root, N_("Split Values"));
8341 for (size_t i = 0; i < sn; i++)
8342 pivot_category_create_leaf (splits, pivot_value_new_variable (
8343 dict_get_var (d, ss + i)));
8347 struct pivot_category *factors = pivot_category_create_group (
8348 attr_dimension->root, N_("Factors"));
8349 for (size_t i = 0; i < fn; i++)
8350 pivot_category_create_leaf (factors, pivot_value_new_variable (
8351 dict_get_var (d, fs + i)));
8353 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
8354 N_("Rows"), N_("Columns"));
8357 struct ccase **rows = NULL;
8358 size_t allocated_rows = 0;
8362 while ((c = casereader_read (r)) != NULL)
8364 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
8374 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
8375 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
8376 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
8379 if (change != NOTHING_CHANGED)
8381 matrix_mget_commit_var (rows, n_rows, d,
8382 rowtype_, &mget->rowtypes,
8393 if (n_rows >= allocated_rows)
8394 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
8397 if (change == SPLITS_CHANGED)
8403 /* Reset the factor number, if there are factors. */
8407 if (row_has_factors)
8413 else if (change == FACTORS_CHANGED)
8415 if (row_has_factors)
8421 matrix_mget_commit_var (rows, n_rows, d,
8422 rowtype_, &mget->rowtypes,
8434 if (var_dimension->n_leaves)
8435 pivot_table_submit (pt);
8437 pivot_table_unref (pt);
8441 matrix_mget_execute (struct matrix_command *cmd)
8443 struct matrix_mget *mget = &cmd->mget;
8444 struct casereader *r;
8445 struct dictionary *d;
8446 if (matrix_open_casereader (cmd, "MGET", mget->file, mget->encoding,
8447 mget->state->dataset, &r, &d))
8449 matrix_mget_execute__ (cmd, r, d);
8450 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
8457 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
8459 if (!lex_force_id (s->lexer))
8462 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
8464 *varp = matrix_var_create (s, lex_tokss (s->lexer));
8469 static struct matrix_command *
8470 matrix_eigen_parse (struct matrix_state *s)
8472 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8473 *cmd = (struct matrix_command) {
8475 .eigen = { .expr = NULL }
8478 struct matrix_eigen *eigen = &cmd->eigen;
8479 if (!lex_force_match (s->lexer, T_LPAREN))
8481 eigen->expr = matrix_expr_parse (s);
8483 || !lex_force_match (s->lexer, T_COMMA)
8484 || !matrix_parse_dst_var (s, &eigen->evec)
8485 || !lex_force_match (s->lexer, T_COMMA)
8486 || !matrix_parse_dst_var (s, &eigen->eval)
8487 || !lex_force_match (s->lexer, T_RPAREN))
8493 matrix_command_destroy (cmd);
8498 matrix_eigen_execute (struct matrix_command *cmd)
8500 struct matrix_eigen *eigen = &cmd->eigen;
8501 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
8504 if (!is_symmetric (A))
8506 msg_at (SE, cmd->location, _("Argument of EIGEN must be symmetric."));
8507 gsl_matrix_free (A);
8511 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
8512 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
8513 gsl_vector v_eval = to_vector (eval);
8514 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
8515 gsl_eigen_symmv (A, &v_eval, evec, w);
8516 gsl_eigen_symmv_free (w);
8518 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
8520 gsl_matrix_free (A);
8522 gsl_matrix_free (eigen->eval->value);
8523 eigen->eval->value = eval;
8525 gsl_matrix_free (eigen->evec->value);
8526 eigen->evec->value = evec;
8531 static struct matrix_command *
8532 matrix_setdiag_parse (struct matrix_state *s)
8534 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8535 *cmd = (struct matrix_command) {
8536 .type = MCMD_SETDIAG,
8537 .setdiag = { .dst = NULL }
8540 struct matrix_setdiag *setdiag = &cmd->setdiag;
8541 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
8544 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
8547 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
8552 if (!lex_force_match (s->lexer, T_COMMA))
8555 setdiag->expr = matrix_expr_parse (s);
8559 if (!lex_force_match (s->lexer, T_RPAREN))
8565 matrix_command_destroy (cmd);
8570 matrix_setdiag_execute (struct matrix_command *cmd)
8572 struct matrix_setdiag *setdiag = &cmd->setdiag;
8573 gsl_matrix *dst = setdiag->dst->value;
8576 msg_at (SE, cmd->location,
8577 _("SETDIAG destination matrix %s is uninitialized."),
8578 setdiag->dst->name);
8582 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
8586 size_t n = MIN (dst->size1, dst->size2);
8587 if (is_scalar (src))
8589 double d = to_scalar (src);
8590 for (size_t i = 0; i < n; i++)
8591 gsl_matrix_set (dst, i, i, d);
8593 else if (is_vector (src))
8595 gsl_vector v = to_vector (src);
8596 for (size_t i = 0; i < n && i < v.size; i++)
8597 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
8600 msg_at (SE, matrix_expr_location (setdiag->expr),
8601 _("SETDIAG argument 2 must be a scalar or a vector, "
8602 "not a %zu×%zu matrix."),
8603 src->size1, src->size2);
8604 gsl_matrix_free (src);
8609 static struct matrix_command *
8610 matrix_svd_parse (struct matrix_state *s)
8612 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8613 *cmd = (struct matrix_command) {
8615 .svd = { .expr = NULL }
8618 struct matrix_svd *svd = &cmd->svd;
8619 if (!lex_force_match (s->lexer, T_LPAREN))
8621 svd->expr = matrix_expr_parse (s);
8623 || !lex_force_match (s->lexer, T_COMMA)
8624 || !matrix_parse_dst_var (s, &svd->u)
8625 || !lex_force_match (s->lexer, T_COMMA)
8626 || !matrix_parse_dst_var (s, &svd->s)
8627 || !lex_force_match (s->lexer, T_COMMA)
8628 || !matrix_parse_dst_var (s, &svd->v)
8629 || !lex_force_match (s->lexer, T_RPAREN))
8635 matrix_command_destroy (cmd);
8640 matrix_svd_execute (struct matrix_svd *svd)
8642 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8646 if (m->size1 >= m->size2)
8649 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8650 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8651 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8652 gsl_vector *work = gsl_vector_alloc (A->size2);
8653 gsl_linalg_SV_decomp (A, V, &Sv, work);
8654 gsl_vector_free (work);
8656 matrix_var_set (svd->u, A);
8657 matrix_var_set (svd->s, S);
8658 matrix_var_set (svd->v, V);
8662 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8663 gsl_matrix_transpose_memcpy (At, m);
8664 gsl_matrix_free (m);
8666 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8667 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8668 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8669 gsl_vector *work = gsl_vector_alloc (At->size2);
8670 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8671 gsl_vector_free (work);
8673 matrix_var_set (svd->v, At);
8674 matrix_var_set (svd->s, St);
8675 matrix_var_set (svd->u, Vt);
8679 /* The main MATRIX command logic. */
8682 matrix_command_execute (struct matrix_command *cmd)
8687 matrix_compute_execute (cmd);
8691 matrix_print_execute (&cmd->print);
8695 return matrix_do_if_execute (&cmd->do_if);
8698 matrix_loop_execute (&cmd->loop);
8705 matrix_display_execute (&cmd->display);
8709 matrix_release_execute (&cmd->release);
8713 matrix_save_execute (cmd);
8717 matrix_read_execute (cmd);
8721 matrix_write_execute (&cmd->write);
8725 matrix_get_execute (cmd);
8729 matrix_msave_execute (cmd);
8733 matrix_mget_execute (cmd);
8737 matrix_eigen_execute (cmd);
8741 matrix_setdiag_execute (cmd);
8745 matrix_svd_execute (&cmd->svd);
8753 matrix_command_destroy (struct matrix_command *cmd)
8758 msg_location_destroy (cmd->location);
8763 matrix_lvalue_destroy (cmd->compute.lvalue);
8764 matrix_expr_destroy (cmd->compute.rvalue);
8768 matrix_expr_destroy (cmd->print.expression);
8769 free (cmd->print.title);
8770 print_labels_destroy (cmd->print.rlabels);
8771 print_labels_destroy (cmd->print.clabels);
8775 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8777 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8778 matrix_commands_uninit (&cmd->do_if.clauses[i].commands);
8780 free (cmd->do_if.clauses);
8784 matrix_expr_destroy (cmd->loop.start);
8785 matrix_expr_destroy (cmd->loop.end);
8786 matrix_expr_destroy (cmd->loop.increment);
8787 matrix_expr_destroy (cmd->loop.top_condition);
8788 matrix_expr_destroy (cmd->loop.bottom_condition);
8789 matrix_commands_uninit (&cmd->loop.commands);
8799 free (cmd->release.vars);
8803 matrix_expr_destroy (cmd->save.expression);
8807 matrix_lvalue_destroy (cmd->read.dst);
8808 matrix_expr_destroy (cmd->read.size);
8812 matrix_expr_destroy (cmd->write.expression);
8813 free (cmd->write.format);
8817 matrix_lvalue_destroy (cmd->get.dst);
8818 fh_unref (cmd->get.file);
8819 free (cmd->get.encoding);
8820 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8824 matrix_expr_destroy (cmd->msave.expr);
8828 fh_unref (cmd->mget.file);
8829 stringi_set_destroy (&cmd->mget.rowtypes);
8833 matrix_expr_destroy (cmd->eigen.expr);
8837 matrix_expr_destroy (cmd->setdiag.expr);
8841 matrix_expr_destroy (cmd->svd.expr);
8848 matrix_commands_parse (struct matrix_state *s, struct matrix_commands *c,
8849 const char *command_name,
8850 const char *stop1, const char *stop2)
8852 lex_end_of_command (s->lexer);
8853 lex_discard_rest_of_command (s->lexer);
8855 size_t allocated = 0;
8858 while (lex_token (s->lexer) == T_ENDCMD)
8861 if (lex_at_phrase (s->lexer, stop1)
8862 || (stop2 && lex_at_phrase (s->lexer, stop2)))
8865 if (lex_at_phrase (s->lexer, "END MATRIX"))
8867 msg (SE, _("Premature END MATRIX within %s."), command_name);
8871 struct matrix_command *cmd = matrix_command_parse (s);
8875 if (c->n >= allocated)
8876 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
8877 c->commands[c->n++] = cmd;
8882 matrix_commands_uninit (struct matrix_commands *cmds)
8884 for (size_t i = 0; i < cmds->n; i++)
8885 matrix_command_destroy (cmds->commands[i]);
8886 free (cmds->commands);
8889 struct matrix_command_name
8892 struct matrix_command *(*parse) (struct matrix_state *);
8895 static const struct matrix_command_name *
8896 matrix_command_name_parse (struct lexer *lexer)
8898 static const struct matrix_command_name commands[] = {
8899 { "COMPUTE", matrix_compute_parse },
8900 { "CALL EIGEN", matrix_eigen_parse },
8901 { "CALL SETDIAG", matrix_setdiag_parse },
8902 { "CALL SVD", matrix_svd_parse },
8903 { "PRINT", matrix_print_parse },
8904 { "DO IF", matrix_do_if_parse },
8905 { "LOOP", matrix_loop_parse },
8906 { "BREAK", matrix_break_parse },
8907 { "READ", matrix_read_parse },
8908 { "WRITE", matrix_write_parse },
8909 { "GET", matrix_get_parse },
8910 { "SAVE", matrix_save_parse },
8911 { "MGET", matrix_mget_parse },
8912 { "MSAVE", matrix_msave_parse },
8913 { "DISPLAY", matrix_display_parse },
8914 { "RELEASE", matrix_release_parse },
8916 static size_t n = sizeof commands / sizeof *commands;
8918 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8919 if (lex_match_phrase (lexer, c->name))
8924 static struct matrix_command *
8925 matrix_command_parse (struct matrix_state *s)
8927 int start_ofs = lex_ofs (s->lexer);
8928 size_t nesting_level = SIZE_MAX;
8930 struct matrix_command *c = NULL;
8931 const struct matrix_command_name *cmd = matrix_command_name_parse (s->lexer);
8933 lex_error (s->lexer, _("Unknown matrix command."));
8934 else if (!cmd->parse)
8935 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8939 nesting_level = output_open_group (
8940 group_item_create_nocopy (utf8_to_title (cmd->name),
8941 utf8_to_title (cmd->name)));
8947 c->location = lex_ofs_location (s->lexer, start_ofs, lex_ofs (s->lexer));
8948 msg_location_remove_columns (c->location);
8949 lex_end_of_command (s->lexer);
8951 lex_discard_rest_of_command (s->lexer);
8952 if (nesting_level != SIZE_MAX)
8953 output_close_groups (nesting_level);
8959 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8961 if (!lex_force_match (lexer, T_ENDCMD))
8964 struct matrix_state state = {
8966 .session = dataset_session (ds),
8968 .vars = HMAP_INITIALIZER (state.vars),
8973 while (lex_match (lexer, T_ENDCMD))
8975 if (lex_token (lexer) == T_STOP)
8977 msg (SE, _("Unexpected end of input expecting matrix command."));
8981 if (lex_match_phrase (lexer, "END MATRIX"))
8984 struct matrix_command *c = matrix_command_parse (&state);
8987 matrix_command_execute (c);
8988 matrix_command_destroy (c);
8992 struct matrix_var *var, *next;
8993 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8996 gsl_matrix_free (var->value);
8997 hmap_delete (&state.vars, &var->hmap_node);
9000 hmap_destroy (&state.vars);
9001 msave_common_destroy (state.msave_common);
9002 fh_unref (state.prev_read_file);
9003 for (size_t i = 0; i < state.n_read_files; i++)
9004 read_file_destroy (state.read_files[i]);
9005 free (state.read_files);
9006 fh_unref (state.prev_write_file);
9007 for (size_t i = 0; i < state.n_write_files; i++)
9008 write_file_destroy (state.write_files[i]);
9009 free (state.write_files);
9010 fh_unref (state.prev_save_file);
9011 for (size_t i = 0; i < state.n_save_files; i++)
9012 save_file_destroy (state.save_files[i]);
9013 free (state.save_files);