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)
5370 int start_ofs = lex_ofs (s->lexer);
5371 cmd->print.expression = matrix_parse_exp (s);
5372 if (!cmd->print.expression)
5374 cmd->print.title = lex_ofs_representation (s->lexer, start_ofs,
5375 lex_ofs (s->lexer) - 1);
5378 while (lex_match (s->lexer, T_SLASH))
5380 if (lex_match_id (s->lexer, "FORMAT"))
5382 lex_match (s->lexer, T_EQUALS);
5383 if (!parse_format_specifier (s->lexer, &cmd->print.format))
5385 cmd->print.use_default_format = false;
5387 else if (lex_match_id (s->lexer, "TITLE"))
5389 lex_match (s->lexer, T_EQUALS);
5390 if (!lex_force_string (s->lexer))
5392 free (cmd->print.title);
5393 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
5396 else if (lex_match_id (s->lexer, "SPACE"))
5398 lex_match (s->lexer, T_EQUALS);
5399 if (lex_match_id (s->lexer, "NEWPAGE"))
5400 cmd->print.space = -1;
5401 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
5403 cmd->print.space = lex_integer (s->lexer);
5409 else if (lex_match_id (s->lexer, "RLABELS"))
5410 parse_literal_print_labels (s, &cmd->print.rlabels);
5411 else if (lex_match_id (s->lexer, "CLABELS"))
5412 parse_literal_print_labels (s, &cmd->print.clabels);
5413 else if (lex_match_id (s->lexer, "RNAMES"))
5415 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
5418 else if (lex_match_id (s->lexer, "CNAMES"))
5420 if (!parse_expr_print_labels (s, &cmd->print.clabels))
5425 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
5426 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
5434 matrix_command_destroy (cmd);
5439 matrix_is_integer (const gsl_matrix *m)
5441 for (size_t y = 0; y < m->size1; y++)
5442 for (size_t x = 0; x < m->size2; x++)
5444 double d = gsl_matrix_get (m, y, x);
5452 matrix_max_magnitude (const gsl_matrix *m)
5455 for (size_t y = 0; y < m->size1; y++)
5456 for (size_t x = 0; x < m->size2; x++)
5458 double d = fabs (gsl_matrix_get (m, y, x));
5466 format_fits (struct fmt_spec format, double x)
5468 char *s = data_out (&(union value) { .f = x }, NULL,
5469 &format, settings_get_fmt_settings ());
5470 bool fits = *s != '*' && !strchr (s, 'E');
5475 static struct fmt_spec
5476 default_format (const gsl_matrix *m, int *log_scale)
5480 double max = matrix_max_magnitude (m);
5482 if (matrix_is_integer (m))
5483 for (int w = 1; w <= 12; w++)
5485 struct fmt_spec format = { .type = FMT_F, .w = w };
5486 if (format_fits (format, -max))
5490 if (max >= 1e9 || max <= 1e-4)
5493 snprintf (s, sizeof s, "%.1e", max);
5495 const char *e = strchr (s, 'e');
5497 *log_scale = atoi (e + 1);
5500 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5504 trimmed_string (double d)
5506 struct substring s = ss_buffer ((char *) &d, sizeof d);
5507 ss_rtrim (&s, ss_cstr (" "));
5508 return ss_xstrdup (s);
5511 static struct string_array *
5512 print_labels_get (const struct print_labels *labels, size_t n,
5513 const char *prefix, bool truncate)
5518 struct string_array *out = xzalloc (sizeof *out);
5519 if (labels->literals.n)
5520 string_array_clone (out, &labels->literals);
5521 else if (labels->expr)
5523 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5524 if (m && is_vector (m))
5526 gsl_vector v = to_vector (m);
5527 for (size_t i = 0; i < v.size; i++)
5528 string_array_append_nocopy (out, trimmed_string (
5529 gsl_vector_get (&v, i)));
5531 gsl_matrix_free (m);
5537 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5540 string_array_append (out, "");
5543 string_array_delete (out, out->n - 1);
5546 for (size_t i = 0; i < out->n; i++)
5548 char *s = out->strings[i];
5549 s[strnlen (s, 8)] = '\0';
5556 matrix_print_space (int space)
5559 output_item_submit (page_break_item_create ());
5560 for (int i = 0; i < space; i++)
5561 output_log ("%s", "");
5565 matrix_print_text (const struct matrix_print *print, const gsl_matrix *m,
5566 struct fmt_spec format, int log_scale)
5568 matrix_print_space (print->space);
5570 output_log ("%s", print->title);
5572 output_log (" 10 ** %d X", log_scale);
5574 struct string_array *clabels = print_labels_get (print->clabels,
5575 m->size2, "col", true);
5576 if (clabels && format.w < 8)
5578 struct string_array *rlabels = print_labels_get (print->rlabels,
5579 m->size1, "row", true);
5583 struct string line = DS_EMPTY_INITIALIZER;
5585 ds_put_byte_multiple (&line, ' ', 8);
5586 for (size_t x = 0; x < m->size2; x++)
5587 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5588 output_log_nocopy (ds_steal_cstr (&line));
5591 double scale = pow (10.0, log_scale);
5592 bool numeric = fmt_is_numeric (format.type);
5593 for (size_t y = 0; y < m->size1; y++)
5595 struct string line = DS_EMPTY_INITIALIZER;
5597 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5599 for (size_t x = 0; x < m->size2; x++)
5601 double f = gsl_matrix_get (m, y, x);
5603 ? data_out (&(union value) { .f = f / scale}, NULL,
5604 &format, settings_get_fmt_settings ())
5605 : trimmed_string (f));
5606 ds_put_format (&line, " %s", s);
5609 output_log_nocopy (ds_steal_cstr (&line));
5612 string_array_destroy (rlabels);
5614 string_array_destroy (clabels);
5619 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5620 const struct print_labels *print_labels, size_t n,
5621 const char *name, const char *prefix)
5623 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5625 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5626 for (size_t i = 0; i < n; i++)
5627 pivot_category_create_leaf (
5629 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5630 : pivot_value_new_integer (i + 1)));
5632 d->hide_all_labels = true;
5633 string_array_destroy (labels);
5638 matrix_print_tables (const struct matrix_print *print, const gsl_matrix *m,
5639 struct fmt_spec format, int log_scale)
5641 struct pivot_table *table = pivot_table_create__ (
5642 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5645 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5647 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5648 N_("Columns"), "col");
5650 struct pivot_footnote *footnote = NULL;
5653 char *s = xasprintf ("× 10**%d", log_scale);
5654 footnote = pivot_table_create_footnote (
5655 table, pivot_value_new_user_text_nocopy (s));
5658 double scale = pow (10.0, log_scale);
5659 bool numeric = fmt_is_numeric (format.type);
5660 for (size_t y = 0; y < m->size1; y++)
5661 for (size_t x = 0; x < m->size2; x++)
5663 double f = gsl_matrix_get (m, y, x);
5664 struct pivot_value *v;
5667 v = pivot_value_new_number (f / scale);
5668 v->numeric.format = format;
5671 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5673 pivot_value_add_footnote (v, footnote);
5674 pivot_table_put2 (table, y, x, v);
5677 pivot_table_submit (table);
5681 matrix_print_execute (const struct matrix_print *print)
5683 if (print->expression)
5685 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5690 struct fmt_spec format = (print->use_default_format
5691 ? default_format (m, &log_scale)
5694 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5695 matrix_print_text (print, m, format, log_scale);
5697 matrix_print_tables (print, m, format, log_scale);
5699 gsl_matrix_free (m);
5703 matrix_print_space (print->space);
5705 output_log ("%s", print->title);
5712 matrix_do_if_clause_parse (struct matrix_state *s,
5713 struct matrix_do_if *ifc,
5714 bool parse_condition,
5715 size_t *allocated_clauses)
5717 if (ifc->n_clauses >= *allocated_clauses)
5718 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5719 sizeof *ifc->clauses);
5720 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5721 *c = (struct do_if_clause) { .condition = NULL };
5723 if (parse_condition)
5725 c->condition = matrix_expr_parse (s);
5730 return matrix_commands_parse (s, &c->commands, "DO IF", "ELSE", "END IF");
5733 static struct matrix_command *
5734 matrix_do_if_parse (struct matrix_state *s)
5736 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5737 *cmd = (struct matrix_command) {
5739 .do_if = { .n_clauses = 0 }
5742 size_t allocated_clauses = 0;
5745 if (!matrix_do_if_clause_parse (s, &cmd->do_if, true, &allocated_clauses))
5748 while (lex_match_phrase (s->lexer, "ELSE IF"));
5750 if (lex_match_id (s->lexer, "ELSE")
5751 && !matrix_do_if_clause_parse (s, &cmd->do_if, false, &allocated_clauses))
5754 if (!lex_match_phrase (s->lexer, "END IF"))
5759 matrix_command_destroy (cmd);
5764 matrix_do_if_execute (struct matrix_do_if *cmd)
5766 for (size_t i = 0; i < cmd->n_clauses; i++)
5768 struct do_if_clause *c = &cmd->clauses[i];
5772 if (!matrix_expr_evaluate_scalar (c->condition,
5773 i ? "ELSE IF" : "DO IF",
5778 for (size_t j = 0; j < c->commands.n; j++)
5779 if (!matrix_command_execute (c->commands.commands[j]))
5788 static struct matrix_command *
5789 matrix_loop_parse (struct matrix_state *s)
5791 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5792 *cmd = (struct matrix_command) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5794 struct matrix_loop *loop = &cmd->loop;
5795 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5797 struct substring name = lex_tokss (s->lexer);
5798 loop->var = matrix_var_lookup (s, name);
5800 loop->var = matrix_var_create (s, name);
5805 loop->start = matrix_expr_parse (s);
5806 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5809 loop->end = matrix_expr_parse (s);
5813 if (lex_match (s->lexer, T_BY))
5815 loop->increment = matrix_expr_parse (s);
5816 if (!loop->increment)
5821 if (lex_match_id (s->lexer, "IF"))
5823 loop->top_condition = matrix_expr_parse (s);
5824 if (!loop->top_condition)
5828 bool was_in_loop = s->in_loop;
5830 bool ok = matrix_commands_parse (s, &loop->commands, "LOOP",
5832 s->in_loop = was_in_loop;
5836 if (!lex_match_phrase (s->lexer, "END LOOP"))
5839 if (lex_match_id (s->lexer, "IF"))
5841 loop->bottom_condition = matrix_expr_parse (s);
5842 if (!loop->bottom_condition)
5849 matrix_command_destroy (cmd);
5854 matrix_loop_execute (struct matrix_loop *cmd)
5858 long int increment = 1;
5861 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5862 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5864 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5868 if (increment > 0 ? value > end
5869 : increment < 0 ? value < end
5874 int mxloops = settings_get_mxloops ();
5875 for (int i = 0; i < mxloops; i++)
5880 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5882 gsl_matrix_free (cmd->var->value);
5883 cmd->var->value = NULL;
5885 if (!cmd->var->value)
5886 cmd->var->value = gsl_matrix_alloc (1, 1);
5888 gsl_matrix_set (cmd->var->value, 0, 0, value);
5891 if (cmd->top_condition)
5894 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5899 for (size_t j = 0; j < cmd->commands.n; j++)
5900 if (!matrix_command_execute (cmd->commands.commands[j]))
5903 if (cmd->bottom_condition)
5906 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5907 "END LOOP IF", &d) || d > 0)
5914 if (increment > 0 ? value > end : value < end)
5922 static struct matrix_command *
5923 matrix_break_parse (struct matrix_state *s)
5927 msg (SE, _("BREAK not inside LOOP."));
5931 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5932 *cmd = (struct matrix_command) { .type = MCMD_BREAK };
5938 static struct matrix_command *
5939 matrix_display_parse (struct matrix_state *s)
5941 while (lex_token (s->lexer) != T_ENDCMD)
5943 if (!lex_match_id (s->lexer, "DICTIONARY")
5944 && !lex_match_id (s->lexer, "STATUS"))
5946 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5951 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5952 *cmd = (struct matrix_command) { .type = MCMD_DISPLAY, .display = { s } };
5957 compare_matrix_var_pointers (const void *a_, const void *b_)
5959 const struct matrix_var *const *ap = a_;
5960 const struct matrix_var *const *bp = b_;
5961 const struct matrix_var *a = *ap;
5962 const struct matrix_var *b = *bp;
5963 return strcmp (a->name, b->name);
5967 matrix_display_execute (struct matrix_display *cmd)
5969 const struct matrix_state *s = cmd->state;
5971 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5972 struct pivot_dimension *attr_dimension
5973 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5974 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5975 N_("Rows"), N_("Columns"));
5976 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5978 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5981 const struct matrix_var *var;
5982 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5984 vars[n_vars++] = var;
5985 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5987 struct pivot_dimension *rows = pivot_dimension_create (
5988 table, PIVOT_AXIS_ROW, N_("Variable"));
5989 for (size_t i = 0; i < n_vars; i++)
5991 const struct matrix_var *var = vars[i];
5992 pivot_category_create_leaf (
5993 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
5995 size_t r = var->value->size1;
5996 size_t c = var->value->size2;
5997 double values[] = { r, c, r * c * sizeof (double) / 1024 };
5998 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
5999 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
6002 pivot_table_submit (table);
6007 static struct matrix_command *
6008 matrix_release_parse (struct matrix_state *s)
6010 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6011 *cmd = (struct matrix_command) { .type = MCMD_RELEASE };
6013 size_t allocated_vars = 0;
6014 while (lex_token (s->lexer) == T_ID)
6016 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
6019 if (cmd->release.n_vars >= allocated_vars)
6020 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
6021 sizeof *cmd->release.vars);
6022 cmd->release.vars[cmd->release.n_vars++] = var;
6025 lex_error (s->lexer, _("Variable name expected."));
6028 if (!lex_match (s->lexer, T_COMMA))
6036 matrix_release_execute (struct matrix_release *cmd)
6038 for (size_t i = 0; i < cmd->n_vars; i++)
6040 struct matrix_var *var = cmd->vars[i];
6041 gsl_matrix_free (var->value);
6048 static struct save_file *
6049 save_file_create (struct matrix_state *s, struct file_handle *fh,
6050 struct string_array *variables,
6051 struct matrix_expr *names,
6052 struct stringi_set *strings)
6054 for (size_t i = 0; i < s->n_save_files; i++)
6056 struct save_file *sf = s->save_files[i];
6057 if (fh_equal (sf->file, fh))
6061 string_array_destroy (variables);
6062 matrix_expr_destroy (names);
6063 stringi_set_destroy (strings);
6069 struct save_file *sf = xmalloc (sizeof *sf);
6070 *sf = (struct save_file) {
6072 .dataset = s->dataset,
6073 .variables = *variables,
6075 .strings = STRINGI_SET_INITIALIZER (sf->strings),
6078 stringi_set_swap (&sf->strings, strings);
6079 stringi_set_destroy (strings);
6081 s->save_files = xrealloc (s->save_files,
6082 (s->n_save_files + 1) * sizeof *s->save_files);
6083 s->save_files[s->n_save_files++] = sf;
6087 static struct casewriter *
6088 save_file_open (struct save_file *sf, gsl_matrix *m,
6089 const struct msg_location *save_location)
6091 if (sf->writer || sf->error)
6095 size_t n_variables = caseproto_get_n_widths (
6096 casewriter_get_proto (sf->writer));
6097 if (m->size2 != n_variables)
6099 const char *file_name = (sf->file == fh_inline_file () ? "*"
6100 : fh_get_name (sf->file));
6101 msg_at (SE, save_location,
6102 _("Cannot save %zu×%zu matrix to %s because the "
6103 "first SAVE to %s in this matrix program wrote a "
6104 "%zu-column matrix."),
6105 m->size1, m->size2, file_name, file_name, n_variables);
6106 msg_at (SE, sf->location,
6107 _("This is the location of the first SAVE to %s."),
6116 struct dictionary *dict = dict_create (get_default_encoding ());
6118 /* Fill 'names' with user-specified names if there were any, either from
6119 sf->variables or sf->names. */
6120 struct string_array names = { .n = 0 };
6121 if (sf->variables.n)
6122 string_array_clone (&names, &sf->variables);
6125 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
6126 if (nm && is_vector (nm))
6128 gsl_vector nv = to_vector (nm);
6129 for (size_t i = 0; i < nv.size; i++)
6131 char *name = trimmed_string (gsl_vector_get (&nv, i));
6132 if (dict_id_is_valid (dict, name, true))
6133 string_array_append_nocopy (&names, name);
6138 gsl_matrix_free (nm);
6141 struct stringi_set strings;
6142 stringi_set_clone (&strings, &sf->strings);
6144 for (size_t i = 0; dict_get_n_vars (dict) < m->size2; i++)
6149 name = names.strings[i];
6152 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
6156 int width = stringi_set_delete (&strings, name) ? 8 : 0;
6157 struct variable *var = dict_create_var (dict, name, width);
6160 msg_at (ME, save_location,
6161 _("Duplicate variable name %s in SAVE statement."), name);
6166 if (!stringi_set_is_empty (&strings))
6168 size_t count = stringi_set_count (&strings);
6169 const char *example = stringi_set_node_get_string (
6170 stringi_set_first (&strings));
6173 msg_at (ME, save_location,
6174 _("The SAVE command STRINGS subcommand specifies an unknown "
6175 "variable %s."), example);
6177 msg_at (ME, save_location,
6178 ngettext ("The SAVE command STRINGS subcommand specifies %zu "
6179 "unknown variable: %s.",
6180 "The SAVE command STRINGS subcommand specifies %zu "
6181 "unknown variables, including %s.",
6186 stringi_set_destroy (&strings);
6187 string_array_destroy (&names);
6196 if (sf->file == fh_inline_file ())
6197 sf->writer = autopaging_writer_create (dict_get_proto (dict));
6199 sf->writer = any_writer_open (sf->file, dict);
6203 sf->location = msg_location_dup (save_location);
6214 save_file_destroy (struct save_file *sf)
6218 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
6220 dataset_set_dict (sf->dataset, sf->dict);
6221 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
6225 casewriter_destroy (sf->writer);
6226 dict_unref (sf->dict);
6228 fh_unref (sf->file);
6229 string_array_destroy (&sf->variables);
6230 matrix_expr_destroy (sf->names);
6231 stringi_set_destroy (&sf->strings);
6232 msg_location_destroy (sf->location);
6237 static struct matrix_command *
6238 matrix_save_parse (struct matrix_state *s)
6240 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6241 *cmd = (struct matrix_command) {
6243 .save = { .expression = NULL },
6246 struct file_handle *fh = NULL;
6247 struct matrix_save *save = &cmd->save;
6249 struct string_array variables = STRING_ARRAY_INITIALIZER;
6250 struct matrix_expr *names = NULL;
6251 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
6253 save->expression = matrix_parse_exp (s);
6254 if (!save->expression)
6257 while (lex_match (s->lexer, T_SLASH))
6259 if (lex_match_id (s->lexer, "OUTFILE"))
6261 lex_match (s->lexer, T_EQUALS);
6264 fh = (lex_match (s->lexer, T_ASTERISK)
6266 : fh_parse (s->lexer, FH_REF_FILE, s->session));
6270 else if (lex_match_id (s->lexer, "VARIABLES"))
6272 lex_match (s->lexer, T_EQUALS);
6276 struct dictionary *d = dict_create (get_default_encoding ());
6277 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
6278 PV_NO_SCRATCH | PV_NO_DUPLICATE);
6283 string_array_clear (&variables);
6284 variables = (struct string_array) {
6290 else if (lex_match_id (s->lexer, "NAMES"))
6292 lex_match (s->lexer, T_EQUALS);
6293 matrix_expr_destroy (names);
6294 names = matrix_parse_exp (s);
6298 else if (lex_match_id (s->lexer, "STRINGS"))
6300 lex_match (s->lexer, T_EQUALS);
6301 while (lex_token (s->lexer) == T_ID)
6303 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
6305 if (!lex_match (s->lexer, T_COMMA))
6311 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
6319 if (s->prev_save_file)
6320 fh = fh_ref (s->prev_save_file);
6323 lex_sbc_missing ("OUTFILE");
6327 fh_unref (s->prev_save_file);
6328 s->prev_save_file = fh_ref (fh);
6330 if (variables.n && names)
6332 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
6333 matrix_expr_destroy (names);
6337 save->sf = save_file_create (s, fh, &variables, names, &strings);
6341 string_array_destroy (&variables);
6342 matrix_expr_destroy (names);
6343 stringi_set_destroy (&strings);
6345 matrix_command_destroy (cmd);
6350 matrix_save_execute (const struct matrix_command *cmd)
6352 const struct matrix_save *save = &cmd->save;
6354 gsl_matrix *m = matrix_expr_evaluate (save->expression);
6358 struct casewriter *writer = save_file_open (save->sf, m, cmd->location);
6361 gsl_matrix_free (m);
6365 const struct caseproto *proto = casewriter_get_proto (writer);
6366 for (size_t y = 0; y < m->size1; y++)
6368 struct ccase *c = case_create (proto);
6369 for (size_t x = 0; x < m->size2; x++)
6371 double d = gsl_matrix_get (m, y, x);
6372 int width = caseproto_get_width (proto, x);
6373 union value *value = case_data_rw_idx (c, x);
6377 memcpy (value->s, &d, width);
6379 casewriter_write (writer, c);
6381 gsl_matrix_free (m);
6386 static struct read_file *
6387 read_file_create (struct matrix_state *s, struct file_handle *fh)
6389 for (size_t i = 0; i < s->n_read_files; i++)
6391 struct read_file *rf = s->read_files[i];
6399 struct read_file *rf = xmalloc (sizeof *rf);
6400 *rf = (struct read_file) { .file = fh };
6402 s->read_files = xrealloc (s->read_files,
6403 (s->n_read_files + 1) * sizeof *s->read_files);
6404 s->read_files[s->n_read_files++] = rf;
6408 static struct dfm_reader *
6409 read_file_open (struct read_file *rf)
6412 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
6417 read_file_destroy (struct read_file *rf)
6421 fh_unref (rf->file);
6422 dfm_close_reader (rf->reader);
6423 free (rf->encoding);
6428 static struct matrix_command *
6429 matrix_read_parse (struct matrix_state *s)
6431 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6432 *cmd = (struct matrix_command) {
6434 .read = { .format = FMT_F },
6437 struct file_handle *fh = NULL;
6438 char *encoding = NULL;
6439 struct matrix_read *read = &cmd->read;
6440 read->dst = matrix_lvalue_parse (s);
6445 int repetitions = 0;
6446 int record_width = 0;
6447 bool seen_format = false;
6448 while (lex_match (s->lexer, T_SLASH))
6450 if (lex_match_id (s->lexer, "FILE"))
6452 lex_match (s->lexer, T_EQUALS);
6455 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6459 else if (lex_match_id (s->lexer, "ENCODING"))
6461 lex_match (s->lexer, T_EQUALS);
6462 if (!lex_force_string (s->lexer))
6466 encoding = ss_xstrdup (lex_tokss (s->lexer));
6470 else if (lex_match_id (s->lexer, "FIELD"))
6472 lex_match (s->lexer, T_EQUALS);
6474 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6476 read->c1 = lex_integer (s->lexer);
6478 if (!lex_force_match (s->lexer, T_TO)
6479 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
6481 read->c2 = lex_integer (s->lexer) + 1;
6484 record_width = read->c2 - read->c1;
6485 if (lex_match (s->lexer, T_BY))
6487 if (!lex_force_int_range (s->lexer, "BY", 1,
6488 read->c2 - read->c1))
6490 by = lex_integer (s->lexer);
6493 if (record_width % by)
6495 msg (SE, _("BY %d does not evenly divide record width %d."),
6503 else if (lex_match_id (s->lexer, "SIZE"))
6505 lex_match (s->lexer, T_EQUALS);
6506 matrix_expr_destroy (read->size);
6507 read->size = matrix_parse_exp (s);
6511 else if (lex_match_id (s->lexer, "MODE"))
6513 lex_match (s->lexer, T_EQUALS);
6514 if (lex_match_id (s->lexer, "RECTANGULAR"))
6515 read->symmetric = false;
6516 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6517 read->symmetric = true;
6520 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6524 else if (lex_match_id (s->lexer, "REREAD"))
6525 read->reread = true;
6526 else if (lex_match_id (s->lexer, "FORMAT"))
6530 lex_sbc_only_once ("FORMAT");
6535 lex_match (s->lexer, T_EQUALS);
6537 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6540 const char *p = lex_tokcstr (s->lexer);
6541 if (c_isdigit (p[0]))
6543 repetitions = atoi (p);
6544 p += strspn (p, "0123456789");
6545 if (!fmt_from_name (p, &read->format))
6547 lex_error (s->lexer, _("Unknown format %s."), p);
6552 else if (fmt_from_name (p, &read->format))
6556 struct fmt_spec format;
6557 if (!parse_format_specifier (s->lexer, &format))
6559 read->format = format.type;
6565 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6566 "REREAD", "FORMAT");
6573 lex_sbc_missing ("FIELD");
6577 if (!read->dst->n_indexes && !read->size)
6579 msg (SE, _("SIZE is required for reading data into a full matrix "
6580 "(as opposed to a submatrix)."));
6586 if (s->prev_read_file)
6587 fh = fh_ref (s->prev_read_file);
6590 lex_sbc_missing ("FILE");
6594 fh_unref (s->prev_read_file);
6595 s->prev_read_file = fh_ref (fh);
6597 read->rf = read_file_create (s, fh);
6601 free (read->rf->encoding);
6602 read->rf->encoding = encoding;
6606 /* Field width may be specified in multiple ways:
6609 2. The format on FORMAT.
6610 3. The repetition factor on FORMAT.
6612 (2) and (3) are mutually exclusive.
6614 If more than one of these is present, they must agree. If none of them is
6615 present, then free-field format is used.
6617 if (repetitions > record_width)
6619 msg (SE, _("%d repetitions cannot fit in record width %d."),
6620 repetitions, record_width);
6623 int w = (repetitions ? record_width / repetitions
6629 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6630 "which implies field width %d, "
6631 "but BY specifies field width %d."),
6632 repetitions, record_width, w, by);
6634 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6643 matrix_command_destroy (cmd);
6649 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6650 struct substring data, size_t y, size_t x,
6651 int first_column, int last_column, char *error)
6653 int line_number = dfm_get_line_number (reader);
6654 struct msg_location location = {
6655 .file_name = intern_new (dfm_get_file_name (reader)),
6656 .start = { .line = line_number, .column = first_column },
6657 .end = { .line = line_number, .column = last_column },
6659 msg_at (DW, &location, _("Error reading \"%.*s\" as format %s "
6660 "for matrix row %zu, column %zu: %s"),
6661 (int) data.length, data.string, fmt_name (format),
6662 y + 1, x + 1, error);
6663 msg_location_uninit (&location);
6668 matrix_read_set_field (struct matrix_read *read, struct dfm_reader *reader,
6669 gsl_matrix *m, struct substring p, size_t y, size_t x,
6670 const char *line_start)
6672 const char *input_encoding = dfm_reader_get_encoding (reader);
6675 if (fmt_is_numeric (read->format))
6678 error = data_in (p, input_encoding, read->format,
6679 settings_get_fmt_settings (), &v, 0, NULL);
6680 if (!error && v.f == SYSMIS)
6681 error = xstrdup (_("Matrix data may not contain missing value."));
6686 uint8_t s[sizeof (double)];
6687 union value v = { .s = s };
6688 error = data_in (p, input_encoding, read->format,
6689 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6690 memcpy (&f, s, sizeof f);
6695 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6696 int nc = ss_utf8_count_columns (p);
6697 int c2 = c1 + MAX (1, nc) - 1;
6698 parse_error (reader, read->format, p, y, x, c1, c2, error);
6702 gsl_matrix_set (m, y, x, f);
6703 if (read->symmetric && x != y)
6704 gsl_matrix_set (m, x, y, f);
6709 matrix_read_line (struct matrix_command *cmd, struct dfm_reader *reader,
6710 struct substring *line, const char **startp)
6712 struct matrix_read *read = &cmd->read;
6713 if (dfm_eof (reader))
6715 msg_at (SE, cmd->location,
6716 _("Unexpected end of file reading matrix data."));
6719 dfm_expand_tabs (reader);
6720 struct substring record = dfm_get_record (reader);
6721 /* XXX need to recode record into UTF-8 */
6722 *startp = record.string;
6723 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6728 matrix_read (struct matrix_command *cmd, struct dfm_reader *reader,
6731 struct matrix_read *read = &cmd->read;
6732 for (size_t y = 0; y < m->size1; y++)
6734 size_t nx = read->symmetric ? y + 1 : m->size2;
6736 struct substring line = ss_empty ();
6737 const char *line_start = line.string;
6738 for (size_t x = 0; x < nx; x++)
6745 ss_ltrim (&line, ss_cstr (" ,"));
6746 if (!ss_is_empty (line))
6748 if (!matrix_read_line (cmd, reader, &line, &line_start))
6750 dfm_forward_record (reader);
6753 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6757 if (!matrix_read_line (cmd, reader, &line, &line_start))
6759 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6760 int f = x % fields_per_line;
6761 if (f == fields_per_line - 1)
6762 dfm_forward_record (reader);
6764 p = ss_substr (line, read->w * f, read->w);
6767 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6771 dfm_forward_record (reader);
6774 ss_ltrim (&line, ss_cstr (" ,"));
6775 if (!ss_is_empty (line))
6777 int line_number = dfm_get_line_number (reader);
6778 int c1 = utf8_count_columns (line_start,
6779 line.string - line_start) + 1;
6780 int c2 = c1 + ss_utf8_count_columns (line) - 1;
6781 struct msg_location location = {
6782 .file_name = intern_new (dfm_get_file_name (reader)),
6783 .start = { .line = line_number, .column = c1 },
6784 .end = { .line = line_number, .column = c2 },
6786 msg_at (DW, &location,
6787 _("Trailing garbage following data for matrix row %zu."),
6789 msg_location_uninit (&location);
6796 matrix_read_execute (struct matrix_command *cmd)
6798 struct matrix_read *read = &cmd->read;
6799 struct index_vector iv0, iv1;
6800 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6803 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6806 gsl_matrix *m = matrix_expr_evaluate (read->size);
6812 msg_at (SE, matrix_expr_location (read->size),
6813 _("SIZE must evaluate to a scalar or a 2-element vector, "
6814 "not a %zu×%zu matrix."), m->size1, m->size2);
6815 gsl_matrix_free (m);
6816 index_vector_uninit (&iv0);
6817 index_vector_uninit (&iv1);
6821 gsl_vector v = to_vector (m);
6825 d[0] = gsl_vector_get (&v, 0);
6828 else if (v.size == 2)
6830 d[0] = gsl_vector_get (&v, 0);
6831 d[1] = gsl_vector_get (&v, 1);
6835 msg_at (SE, matrix_expr_location (read->size),
6836 _("SIZE must evaluate to a scalar or a 2-element vector, "
6837 "not a %zu×%zu matrix."),
6838 m->size1, m->size2),
6839 gsl_matrix_free (m);
6840 index_vector_uninit (&iv0);
6841 index_vector_uninit (&iv1);
6844 gsl_matrix_free (m);
6846 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6848 msg_at (SE, matrix_expr_location (read->size),
6849 _("Matrix dimensions %g×%g specified on SIZE "
6850 "are outside valid range."),
6852 index_vector_uninit (&iv0);
6853 index_vector_uninit (&iv1);
6861 if (read->dst->n_indexes)
6863 size_t submatrix_size[2];
6864 if (read->dst->n_indexes == 2)
6866 submatrix_size[0] = iv0.n;
6867 submatrix_size[1] = iv1.n;
6869 else if (read->dst->var->value->size1 == 1)
6871 submatrix_size[0] = 1;
6872 submatrix_size[1] = iv0.n;
6876 submatrix_size[0] = iv0.n;
6877 submatrix_size[1] = 1;
6882 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6884 msg_at (SE, cmd->location,
6885 _("Dimensions specified on SIZE differ from dimensions "
6886 "of destination submatrix."));
6887 msg_at (SN, matrix_expr_location (read->size),
6888 _("SIZE specifies dimensions %zu×%zu."),
6890 msg_at (SN, read->dst->full_location,
6891 _("Destination submatrix has dimensions %zu×%zu."),
6892 submatrix_size[0], submatrix_size[1]);
6893 index_vector_uninit (&iv0);
6894 index_vector_uninit (&iv1);
6900 size[0] = submatrix_size[0];
6901 size[1] = submatrix_size[1];
6905 struct dfm_reader *reader = read_file_open (read->rf);
6907 dfm_reread_record (reader, 1);
6909 if (read->symmetric && size[0] != size[1])
6911 msg_at (SE, cmd->location,
6912 _("Cannot read non-square %zu×%zu matrix "
6913 "using READ with MODE=SYMMETRIC."),
6915 index_vector_uninit (&iv0);
6916 index_vector_uninit (&iv1);
6919 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6920 matrix_read (cmd, reader, tmp);
6921 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp, cmd->location);
6926 static struct write_file *
6927 write_file_create (struct matrix_state *s, struct file_handle *fh)
6929 for (size_t i = 0; i < s->n_write_files; i++)
6931 struct write_file *wf = s->write_files[i];
6939 struct write_file *wf = xmalloc (sizeof *wf);
6940 *wf = (struct write_file) { .file = fh };
6942 s->write_files = xrealloc (s->write_files,
6943 (s->n_write_files + 1) * sizeof *s->write_files);
6944 s->write_files[s->n_write_files++] = wf;
6948 static struct dfm_writer *
6949 write_file_open (struct write_file *wf)
6952 wf->writer = dfm_open_writer (wf->file, wf->encoding);
6957 write_file_destroy (struct write_file *wf)
6963 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
6964 wf->held->s.ss.length);
6965 u8_line_destroy (wf->held);
6969 fh_unref (wf->file);
6970 dfm_close_writer (wf->writer);
6971 free (wf->encoding);
6976 static struct matrix_command *
6977 matrix_write_parse (struct matrix_state *s)
6979 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6980 *cmd = (struct matrix_command) {
6984 struct file_handle *fh = NULL;
6985 char *encoding = NULL;
6986 struct matrix_write *write = &cmd->write;
6987 write->expression = matrix_parse_exp (s);
6988 if (!write->expression)
6992 int repetitions = 0;
6993 int record_width = 0;
6994 enum fmt_type format = FMT_F;
6995 bool has_format = false;
6996 while (lex_match (s->lexer, T_SLASH))
6998 if (lex_match_id (s->lexer, "OUTFILE"))
7000 lex_match (s->lexer, T_EQUALS);
7003 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
7007 else if (lex_match_id (s->lexer, "ENCODING"))
7009 lex_match (s->lexer, T_EQUALS);
7010 if (!lex_force_string (s->lexer))
7014 encoding = ss_xstrdup (lex_tokss (s->lexer));
7018 else if (lex_match_id (s->lexer, "FIELD"))
7020 lex_match (s->lexer, T_EQUALS);
7022 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
7024 write->c1 = lex_integer (s->lexer);
7026 if (!lex_force_match (s->lexer, T_TO)
7027 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
7029 write->c2 = lex_integer (s->lexer) + 1;
7032 record_width = write->c2 - write->c1;
7033 if (lex_match (s->lexer, T_BY))
7035 if (!lex_force_int_range (s->lexer, "BY", 1,
7036 write->c2 - write->c1))
7038 by = lex_integer (s->lexer);
7041 if (record_width % by)
7043 msg (SE, _("BY %d does not evenly divide record width %d."),
7051 else if (lex_match_id (s->lexer, "MODE"))
7053 lex_match (s->lexer, T_EQUALS);
7054 if (lex_match_id (s->lexer, "RECTANGULAR"))
7055 write->triangular = false;
7056 else if (lex_match_id (s->lexer, "TRIANGULAR"))
7057 write->triangular = true;
7060 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
7064 else if (lex_match_id (s->lexer, "HOLD"))
7066 else if (lex_match_id (s->lexer, "FORMAT"))
7068 if (has_format || write->format)
7070 lex_sbc_only_once ("FORMAT");
7074 lex_match (s->lexer, T_EQUALS);
7076 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
7079 const char *p = lex_tokcstr (s->lexer);
7080 if (c_isdigit (p[0]))
7082 repetitions = atoi (p);
7083 p += strspn (p, "0123456789");
7084 if (!fmt_from_name (p, &format))
7086 lex_error (s->lexer, _("Unknown format %s."), p);
7092 else if (fmt_from_name (p, &format))
7099 struct fmt_spec spec;
7100 if (!parse_format_specifier (s->lexer, &spec))
7102 write->format = xmemdup (&spec, sizeof spec);
7107 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
7115 lex_sbc_missing ("FIELD");
7121 if (s->prev_write_file)
7122 fh = fh_ref (s->prev_write_file);
7125 lex_sbc_missing ("OUTFILE");
7129 fh_unref (s->prev_write_file);
7130 s->prev_write_file = fh_ref (fh);
7132 write->wf = write_file_create (s, fh);
7136 free (write->wf->encoding);
7137 write->wf->encoding = encoding;
7141 /* Field width may be specified in multiple ways:
7144 2. The format on FORMAT.
7145 3. The repetition factor on FORMAT.
7147 (2) and (3) are mutually exclusive.
7149 If more than one of these is present, they must agree. If none of them is
7150 present, then free-field format is used.
7152 if (repetitions > record_width)
7154 msg (SE, _("%d repetitions cannot fit in record width %d."),
7155 repetitions, record_width);
7158 int w = (repetitions ? record_width / repetitions
7159 : write->format ? write->format->w
7164 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
7165 "which implies field width %d, "
7166 "but BY specifies field width %d."),
7167 repetitions, record_width, w, by);
7169 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
7173 if (w && !write->format)
7175 write->format = xmalloc (sizeof *write->format);
7176 *write->format = (struct fmt_spec) { .type = format, .w = w };
7178 if (!fmt_check_output (write->format))
7182 if (write->format && fmt_var_width (write->format) > sizeof (double))
7184 char s[FMT_STRING_LEN_MAX + 1];
7185 fmt_to_string (write->format, s);
7186 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
7187 s, sizeof (double));
7195 matrix_command_destroy (cmd);
7200 matrix_write_execute (struct matrix_write *write)
7202 gsl_matrix *m = matrix_expr_evaluate (write->expression);
7206 if (write->triangular && m->size1 != m->size2)
7208 msg_at (SE, matrix_expr_location (write->expression),
7209 _("WRITE with MODE=TRIANGULAR requires a square matrix but "
7210 "the matrix to be written has dimensions %zu×%zu."),
7211 m->size1, m->size2);
7212 gsl_matrix_free (m);
7216 struct dfm_writer *writer = write_file_open (write->wf);
7217 if (!writer || !m->size1)
7219 gsl_matrix_free (m);
7223 const struct fmt_settings *settings = settings_get_fmt_settings ();
7224 struct u8_line *line = write->wf->held;
7225 for (size_t y = 0; y < m->size1; y++)
7229 line = xmalloc (sizeof *line);
7230 u8_line_init (line);
7232 size_t nx = write->triangular ? y + 1 : m->size2;
7234 for (size_t x = 0; x < nx; x++)
7237 double f = gsl_matrix_get (m, y, x);
7241 if (fmt_is_numeric (write->format->type))
7244 v.s = (uint8_t *) &f;
7245 s = data_out (&v, NULL, write->format, settings);
7249 s = xmalloc (DBL_BUFSIZE_BOUND);
7250 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
7251 >= DBL_BUFSIZE_BOUND)
7254 size_t len = strlen (s);
7255 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
7256 if (width + x0 > write->c2)
7258 dfm_put_record_utf8 (writer, line->s.ss.string,
7260 u8_line_clear (line);
7263 u8_line_put (line, x0, x0 + width, s, len);
7266 x0 += write->format ? write->format->w : width + 1;
7269 if (y + 1 >= m->size1 && write->hold)
7271 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
7272 u8_line_clear (line);
7276 u8_line_destroy (line);
7280 write->wf->held = line;
7282 gsl_matrix_free (m);
7287 static struct matrix_command *
7288 matrix_get_parse (struct matrix_state *s)
7290 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7291 *cmd = (struct matrix_command) {
7294 .dataset = s->dataset,
7295 .user = { .treatment = MGET_ERROR },
7296 .system = { .treatment = MGET_ERROR },
7300 struct matrix_get *get = &cmd->get;
7301 get->dst = matrix_lvalue_parse (s);
7305 while (lex_match (s->lexer, T_SLASH))
7307 if (lex_match_id (s->lexer, "FILE"))
7309 lex_match (s->lexer, T_EQUALS);
7311 fh_unref (get->file);
7312 if (lex_match (s->lexer, T_ASTERISK))
7316 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7321 else if (lex_match_id (s->lexer, "ENCODING"))
7323 lex_match (s->lexer, T_EQUALS);
7324 if (!lex_force_string (s->lexer))
7327 free (get->encoding);
7328 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
7332 else if (lex_match_id (s->lexer, "VARIABLES"))
7334 lex_match (s->lexer, T_EQUALS);
7338 lex_sbc_only_once ("VARIABLES");
7342 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
7345 else if (lex_match_id (s->lexer, "NAMES"))
7347 lex_match (s->lexer, T_EQUALS);
7348 if (!lex_force_id (s->lexer))
7351 struct substring name = lex_tokss (s->lexer);
7352 get->names = matrix_var_lookup (s, name);
7354 get->names = matrix_var_create (s, name);
7357 else if (lex_match_id (s->lexer, "MISSING"))
7359 lex_match (s->lexer, T_EQUALS);
7360 if (lex_match_id (s->lexer, "ACCEPT"))
7361 get->user.treatment = MGET_ACCEPT;
7362 else if (lex_match_id (s->lexer, "OMIT"))
7363 get->user.treatment = MGET_OMIT;
7364 else if (lex_is_number (s->lexer))
7366 get->user.treatment = MGET_RECODE;
7367 get->user.substitute = lex_number (s->lexer);
7372 lex_error (s->lexer, NULL);
7376 else if (lex_match_id (s->lexer, "SYSMIS"))
7378 lex_match (s->lexer, T_EQUALS);
7379 if (lex_match_id (s->lexer, "OMIT"))
7380 get->system.treatment = MGET_OMIT;
7381 else if (lex_is_number (s->lexer))
7383 get->system.treatment = MGET_RECODE;
7384 get->system.substitute = lex_number (s->lexer);
7389 lex_error (s->lexer, NULL);
7395 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
7396 "MISSING", "SYSMIS");
7401 if (get->user.treatment != MGET_ACCEPT)
7402 get->system.treatment = MGET_ERROR;
7407 matrix_command_destroy (cmd);
7412 matrix_get_execute__ (struct matrix_command *cmd, struct casereader *reader,
7413 const struct dictionary *dict)
7415 struct matrix_get *get = &cmd->get;
7416 struct variable **vars;
7421 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
7422 &vars, &n_vars, PV_NUMERIC))
7427 n_vars = dict_get_n_vars (dict);
7428 vars = xnmalloc (n_vars, sizeof *vars);
7429 for (size_t i = 0; i < n_vars; i++)
7431 struct variable *var = dict_get_var (dict, i);
7432 if (!var_is_numeric (var))
7434 msg_at (SE, cmd->location, _("Variable %s is not numeric."),
7435 var_get_name (var));
7445 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
7446 for (size_t i = 0; i < n_vars; i++)
7448 char s[sizeof (double)];
7450 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
7451 memcpy (&f, s, sizeof f);
7452 gsl_matrix_set (names, i, 0, f);
7455 gsl_matrix_free (get->names->value);
7456 get->names->value = names;
7460 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
7461 long long int casenum = 1;
7463 for (struct ccase *c = casereader_read (reader); c;
7464 c = casereader_read (reader), casenum++)
7466 if (n_rows >= m->size1)
7468 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
7469 for (size_t y = 0; y < n_rows; y++)
7470 for (size_t x = 0; x < n_vars; x++)
7471 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
7472 gsl_matrix_free (m);
7477 for (size_t x = 0; x < n_vars; x++)
7479 const struct variable *var = vars[x];
7480 double d = case_num (c, var);
7483 if (get->system.treatment == MGET_RECODE)
7484 d = get->system.substitute;
7485 else if (get->system.treatment == MGET_OMIT)
7489 msg_at (SE, cmd->location, _("Variable %s in case %lld "
7490 "is system-missing."),
7491 var_get_name (var), casenum);
7495 else if (var_is_num_missing (var, d) == MV_USER)
7497 if (get->user.treatment == MGET_RECODE)
7498 d = get->user.substitute;
7499 else if (get->user.treatment == MGET_OMIT)
7501 else if (get->user.treatment != MGET_ACCEPT)
7503 msg_at (SE, cmd->location,
7504 _("Variable %s in case %lld has user-missing "
7506 var_get_name (var), casenum, d);
7510 gsl_matrix_set (m, n_rows, x, d);
7521 matrix_lvalue_evaluate_and_assign (get->dst, m, cmd->location);
7524 gsl_matrix_free (m);
7529 matrix_open_casereader (const struct matrix_command *cmd,
7530 const char *command_name, struct file_handle *file,
7531 const char *encoding, struct dataset *dataset,
7532 struct casereader **readerp, struct dictionary **dictp)
7536 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
7537 return *readerp != NULL;
7541 if (dict_get_n_vars (dataset_dict (dataset)) == 0)
7543 msg_at (ME, cmd->location,
7544 _("The %s command cannot read an empty active file."),
7548 *readerp = proc_open (dataset);
7549 *dictp = dict_ref (dataset_dict (dataset));
7555 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
7556 struct casereader *reader, struct dictionary *dict)
7559 casereader_destroy (reader);
7561 proc_commit (dataset);
7565 matrix_get_execute (struct matrix_command *cmd)
7567 struct matrix_get *get = &cmd->get;
7568 struct casereader *r;
7569 struct dictionary *d;
7570 if (matrix_open_casereader (cmd, "GET", get->file, get->encoding,
7571 get->dataset, &r, &d))
7573 matrix_get_execute__ (cmd, r, d);
7574 matrix_close_casereader (get->file, get->dataset, r, d);
7581 variables_changed (const char *keyword,
7582 const struct string_array *new,
7583 const struct string_array *old)
7589 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7590 "on the first MSAVE within MATRIX."), keyword);
7593 else if (!string_array_equal_case (old, new))
7595 msg (SE, _("%s must specify the same variables each time within "
7596 "a given MATRIX."), keyword);
7604 msave_common_changed (const struct msave_common *old,
7605 const struct msave_common *new)
7607 if (new->outfile && !fh_equal (old->outfile, new->outfile))
7608 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7609 "within a single MATRIX command."));
7610 else if (variables_changed ("VARIABLES", &new->variables, &old->variables)
7611 || variables_changed ("FNAMES", &new->fnames, &old->fnames)
7612 || variables_changed ("SNAMES", &new->snames, &old->snames))
7613 msg_at (SN, old->location,
7614 _("This is the location of the first MSAVE command."));
7622 msave_common_destroy (struct msave_common *common)
7626 msg_location_destroy (common->location);
7627 fh_unref (common->outfile);
7628 string_array_destroy (&common->variables);
7629 string_array_destroy (&common->fnames);
7630 string_array_destroy (&common->snames);
7632 for (size_t i = 0; i < common->n_factors; i++)
7633 matrix_expr_destroy (common->factors[i]);
7634 free (common->factors);
7636 for (size_t i = 0; i < common->n_splits; i++)
7637 matrix_expr_destroy (common->splits[i]);
7638 free (common->splits);
7640 dict_unref (common->dict);
7641 casewriter_destroy (common->writer);
7648 match_rowtype (struct lexer *lexer)
7650 static const char *rowtypes[] = {
7651 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7653 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7655 for (size_t i = 0; i < n_rowtypes; i++)
7656 if (lex_match_id (lexer, rowtypes[i]))
7659 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7664 parse_var_names (struct lexer *lexer, struct string_array *sa)
7666 lex_match (lexer, T_EQUALS);
7668 string_array_clear (sa);
7670 struct dictionary *dict = dict_create (get_default_encoding ());
7673 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7674 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7679 for (size_t i = 0; i < n_names; i++)
7680 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7681 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7683 msg (SE, _("Variable name %s is reserved."), names[i]);
7684 for (size_t j = 0; j < n_names; j++)
7690 string_array_clear (sa);
7691 sa->strings = names;
7692 sa->n = sa->allocated = n_names;
7697 static struct matrix_command *
7698 matrix_msave_parse (struct matrix_state *s)
7700 int start_ofs = lex_ofs (s->lexer);
7702 struct msave_common *common = xmalloc (sizeof *common);
7703 *common = (struct msave_common) { .outfile = NULL };
7705 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7706 *cmd = (struct matrix_command) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7708 struct matrix_expr *splits = NULL;
7709 struct matrix_expr *factors = NULL;
7711 struct matrix_msave *msave = &cmd->msave;
7712 msave->expr = matrix_parse_exp (s);
7716 while (lex_match (s->lexer, T_SLASH))
7718 if (lex_match_id (s->lexer, "TYPE"))
7720 lex_match (s->lexer, T_EQUALS);
7722 msave->rowtype = match_rowtype (s->lexer);
7723 if (!msave->rowtype)
7726 else if (lex_match_id (s->lexer, "OUTFILE"))
7728 lex_match (s->lexer, T_EQUALS);
7730 fh_unref (common->outfile);
7731 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7732 if (!common->outfile)
7735 else if (lex_match_id (s->lexer, "VARIABLES"))
7737 if (!parse_var_names (s->lexer, &common->variables))
7740 else if (lex_match_id (s->lexer, "FNAMES"))
7742 if (!parse_var_names (s->lexer, &common->fnames))
7745 else if (lex_match_id (s->lexer, "SNAMES"))
7747 if (!parse_var_names (s->lexer, &common->snames))
7750 else if (lex_match_id (s->lexer, "SPLIT"))
7752 lex_match (s->lexer, T_EQUALS);
7754 matrix_expr_destroy (splits);
7755 splits = matrix_parse_exp (s);
7759 else if (lex_match_id (s->lexer, "FACTOR"))
7761 lex_match (s->lexer, T_EQUALS);
7763 matrix_expr_destroy (factors);
7764 factors = matrix_parse_exp (s);
7770 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7771 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7775 if (!msave->rowtype)
7777 lex_sbc_missing ("TYPE");
7781 if (!s->msave_common)
7783 if (common->fnames.n && !factors)
7785 msg (SE, _("FNAMES requires FACTOR."));
7788 if (common->snames.n && !splits)
7790 msg (SE, _("SNAMES requires SPLIT."));
7793 if (!common->outfile)
7795 lex_sbc_missing ("OUTFILE");
7798 common->location = lex_ofs_location (s->lexer, start_ofs,
7799 lex_ofs (s->lexer));
7800 msg_location_remove_columns (common->location);
7801 s->msave_common = common;
7805 if (msave_common_changed (s->msave_common, common))
7807 msave_common_destroy (common);
7809 msave->common = s->msave_common;
7811 struct msave_common *c = s->msave_common;
7814 if (c->n_factors >= c->allocated_factors)
7815 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7816 sizeof *c->factors);
7817 c->factors[c->n_factors++] = factors;
7819 if (c->n_factors > 0)
7820 msave->factors = c->factors[c->n_factors - 1];
7824 if (c->n_splits >= c->allocated_splits)
7825 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7827 c->splits[c->n_splits++] = splits;
7829 if (c->n_splits > 0)
7830 msave->splits = c->splits[c->n_splits - 1];
7835 matrix_expr_destroy (splits);
7836 matrix_expr_destroy (factors);
7837 msave_common_destroy (common);
7838 matrix_command_destroy (cmd);
7843 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7845 gsl_matrix *m = matrix_expr_evaluate (e);
7851 msg_at (SE, matrix_expr_location (e),
7852 _("%s expression must evaluate to vector, "
7853 "not a %zu×%zu matrix."),
7854 name, m->size1, m->size2);
7855 gsl_matrix_free (m);
7859 return matrix_to_vector (m);
7863 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7865 for (size_t i = 0; i < vars->n; i++)
7866 if (!dict_create_var (d, vars->strings[i], 0))
7867 return vars->strings[i];
7871 static struct dictionary *
7872 msave_create_dict (const struct msave_common *common,
7873 const struct msg_location *location)
7875 struct dictionary *dict = dict_create (get_default_encoding ());
7877 const char *dup_split = msave_add_vars (dict, &common->snames);
7880 /* Should not be possible because the parser ensures that the names are
7885 dict_create_var_assert (dict, "ROWTYPE_", 8);
7887 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7890 msg_at (SE, location, _("Duplicate or invalid FACTOR variable name %s."),
7895 dict_create_var_assert (dict, "VARNAME_", 8);
7897 const char *dup_var = msave_add_vars (dict, &common->variables);
7900 msg_at (SE, location, _("Duplicate or invalid variable name %s."),
7913 matrix_msave_execute (struct matrix_command *cmd)
7915 struct matrix_msave *msave = &cmd->msave;
7916 struct msave_common *common = msave->common;
7917 gsl_matrix *m = NULL;
7918 gsl_vector *factors = NULL;
7919 gsl_vector *splits = NULL;
7921 m = matrix_expr_evaluate (msave->expr);
7925 if (!common->variables.n)
7926 for (size_t i = 0; i < m->size2; i++)
7927 string_array_append_nocopy (&common->variables,
7928 xasprintf ("COL%zu", i + 1));
7929 else if (m->size2 != common->variables.n)
7931 msg_at (SE, matrix_expr_location (msave->expr),
7932 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7933 m->size2, common->variables.n);
7939 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7943 if (!common->fnames.n)
7944 for (size_t i = 0; i < factors->size; i++)
7945 string_array_append_nocopy (&common->fnames,
7946 xasprintf ("FAC%zu", i + 1));
7947 else if (factors->size != common->fnames.n)
7949 msg_at (SE, matrix_expr_location (msave->factors),
7950 _("There are %zu factor variables, "
7951 "but %zu factor values were supplied."),
7952 common->fnames.n, factors->size);
7959 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7963 if (!common->snames.n)
7964 for (size_t i = 0; i < splits->size; i++)
7965 string_array_append_nocopy (&common->snames,
7966 xasprintf ("SPL%zu", i + 1));
7967 else if (splits->size != common->snames.n)
7969 msg_at (SE, matrix_expr_location (msave->splits),
7970 _("There are %zu split variables, "
7971 "but %zu split values were supplied."),
7972 common->snames.n, splits->size);
7977 if (!common->writer)
7979 struct dictionary *dict = msave_create_dict (common, cmd->location);
7983 common->writer = any_writer_open (common->outfile, dict);
7984 if (!common->writer)
7990 common->dict = dict;
7993 bool matrix = (!strcmp (msave->rowtype, "COV")
7994 || !strcmp (msave->rowtype, "CORR"));
7995 for (size_t y = 0; y < m->size1; y++)
7997 struct ccase *c = case_create (dict_get_proto (common->dict));
8000 /* Split variables */
8002 for (size_t i = 0; i < splits->size; i++)
8003 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
8006 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8007 msave->rowtype, ' ');
8011 for (size_t i = 0; i < factors->size; i++)
8012 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
8015 const char *varname_ = (matrix && y < common->variables.n
8016 ? common->variables.strings[y]
8018 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8021 /* Continuous variables. */
8022 for (size_t x = 0; x < m->size2; x++)
8023 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
8024 casewriter_write (common->writer, c);
8028 gsl_matrix_free (m);
8029 gsl_vector_free (factors);
8030 gsl_vector_free (splits);
8035 static struct matrix_command *
8036 matrix_mget_parse (struct matrix_state *s)
8038 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8039 *cmd = (struct matrix_command) {
8043 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
8047 struct matrix_mget *mget = &cmd->mget;
8049 lex_match (s->lexer, T_SLASH);
8050 while (lex_token (s->lexer) != T_ENDCMD)
8052 if (lex_match_id (s->lexer, "FILE"))
8054 lex_match (s->lexer, T_EQUALS);
8056 fh_unref (mget->file);
8057 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
8061 else if (lex_match_id (s->lexer, "ENCODING"))
8063 lex_match (s->lexer, T_EQUALS);
8064 if (!lex_force_string (s->lexer))
8067 free (mget->encoding);
8068 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
8072 else if (lex_match_id (s->lexer, "TYPE"))
8074 lex_match (s->lexer, T_EQUALS);
8075 while (lex_token (s->lexer) != T_SLASH
8076 && lex_token (s->lexer) != T_ENDCMD)
8078 const char *rowtype = match_rowtype (s->lexer);
8082 stringi_set_insert (&mget->rowtypes, rowtype);
8087 lex_error_expecting (s->lexer, "FILE", "TYPE");
8090 lex_match (s->lexer, T_SLASH);
8095 matrix_command_destroy (cmd);
8099 static const struct variable *
8100 get_a8_var (const struct msg_location *loc,
8101 const struct dictionary *d, const char *name)
8103 const struct variable *v = dict_lookup_var (d, name);
8106 msg_at (SE, loc, _("Matrix data file lacks %s variable."), name);
8109 if (var_get_width (v) != 8)
8111 msg_at (SE, loc, _("%s variable in matrix data file must be "
8112 "8-byte string, but it has width %d."),
8113 name, var_get_width (v));
8120 var_changed (const struct ccase *ca, const struct ccase *cb,
8121 const struct variable *var)
8124 ? !value_equal (case_data (ca, var), case_data (cb, var),
8125 var_get_width (var))
8130 vars_changed (const struct ccase *ca, const struct ccase *cb,
8131 const struct dictionary *d,
8132 size_t first_var, size_t n_vars)
8134 for (size_t i = 0; i < n_vars; i++)
8136 const struct variable *v = dict_get_var (d, first_var + i);
8137 if (var_changed (ca, cb, v))
8144 vars_all_missing (const struct ccase *c, const struct dictionary *d,
8145 size_t first_var, size_t n_vars)
8147 for (size_t i = 0; i < n_vars; i++)
8148 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
8154 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
8155 const struct dictionary *d,
8156 const struct variable *rowtype_var,
8157 const struct stringi_set *accepted_rowtypes,
8158 struct matrix_state *s,
8159 size_t ss, size_t sn, size_t si,
8160 size_t fs, size_t fn, size_t fi,
8161 size_t cs, size_t cn,
8162 struct pivot_table *pt,
8163 struct pivot_dimension *var_dimension)
8168 /* Is this a matrix for pooled data, either where there are no factor
8169 variables or the factor variables are missing? */
8170 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
8172 struct substring rowtype = case_ss (rows[0], rowtype_var);
8173 ss_rtrim (&rowtype, ss_cstr (" "));
8174 if (!stringi_set_is_empty (accepted_rowtypes)
8175 && !stringi_set_contains_len (accepted_rowtypes,
8176 rowtype.string, rowtype.length))
8179 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
8180 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
8181 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
8182 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
8183 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
8184 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
8188 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
8189 (int) rowtype.length, rowtype.string);
8193 struct string name = DS_EMPTY_INITIALIZER;
8194 ds_put_cstr (&name, prefix);
8196 ds_put_format (&name, "F%zu", fi);
8198 ds_put_format (&name, "S%zu", si);
8200 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
8202 mv = matrix_var_create (s, ds_ss (&name));
8205 msg (SW, _("Matrix data file contains variable with existing name %s."),
8207 goto exit_free_name;
8210 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
8211 size_t n_missing = 0;
8212 for (size_t y = 0; y < n_rows; y++)
8214 for (size_t x = 0; x < cn; x++)
8216 struct variable *var = dict_get_var (d, cs + x);
8217 double value = case_num (rows[y], var);
8218 if (var_is_num_missing (var, value))
8223 gsl_matrix_set (m, y, x, value);
8227 int var_index = pivot_category_create_leaf (
8228 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
8229 double values[] = { n_rows, cn };
8230 for (size_t j = 0; j < sn; j++)
8232 struct variable *var = dict_get_var (d, ss + j);
8233 const union value *value = case_data (rows[0], var);
8234 pivot_table_put2 (pt, j, var_index,
8235 pivot_value_new_var_value (var, value));
8237 for (size_t j = 0; j < fn; j++)
8239 struct variable *var = dict_get_var (d, fs + j);
8240 const union value sysmis = { .f = SYSMIS };
8241 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
8242 pivot_table_put2 (pt, j + sn, var_index,
8243 pivot_value_new_var_value (var, value));
8245 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
8246 pivot_table_put2 (pt, j + sn + fn, var_index,
8247 pivot_value_new_integer (values[j]));
8250 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
8251 "value, which was treated as zero.",
8252 "Matrix data file variable %s contains %zu missing "
8253 "values, which were treated as zero.", n_missing),
8254 ds_cstr (&name), n_missing);
8261 for (size_t y = 0; y < n_rows; y++)
8262 case_unref (rows[y]);
8266 matrix_mget_execute__ (struct matrix_command *cmd, struct casereader *r,
8267 const struct dictionary *d)
8269 struct matrix_mget *mget = &cmd->mget;
8270 const struct msg_location *loc = cmd->location;
8271 const struct variable *rowtype_ = get_a8_var (loc, d, "ROWTYPE_");
8272 const struct variable *varname_ = get_a8_var (loc, d, "VARNAME_");
8273 if (!rowtype_ || !varname_)
8276 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
8279 _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
8282 if (var_get_dict_index (varname_) + 1 >= dict_get_n_vars (d))
8284 msg_at (SE, loc, _("Matrix data file contains no continuous variables."));
8288 for (size_t i = 0; i < dict_get_n_vars (d); i++)
8290 const struct variable *v = dict_get_var (d, i);
8291 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
8294 _("Matrix data file contains unexpected string variable %s."),
8300 /* SPLIT variables. */
8302 size_t sn = var_get_dict_index (rowtype_);
8303 struct ccase *sc = NULL;
8306 /* FACTOR variables. */
8307 size_t fs = var_get_dict_index (rowtype_) + 1;
8308 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
8309 struct ccase *fc = NULL;
8312 /* Continuous variables. */
8313 size_t cs = var_get_dict_index (varname_) + 1;
8314 size_t cn = dict_get_n_vars (d) - cs;
8315 struct ccase *cc = NULL;
8318 struct pivot_table *pt = pivot_table_create (
8319 N_("Matrix Variables Created by MGET"));
8320 struct pivot_dimension *attr_dimension = pivot_dimension_create (
8321 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
8322 struct pivot_dimension *var_dimension = pivot_dimension_create (
8323 pt, PIVOT_AXIS_ROW, N_("Variable"));
8326 struct pivot_category *splits = pivot_category_create_group (
8327 attr_dimension->root, N_("Split Values"));
8328 for (size_t i = 0; i < sn; i++)
8329 pivot_category_create_leaf (splits, pivot_value_new_variable (
8330 dict_get_var (d, ss + i)));
8334 struct pivot_category *factors = pivot_category_create_group (
8335 attr_dimension->root, N_("Factors"));
8336 for (size_t i = 0; i < fn; i++)
8337 pivot_category_create_leaf (factors, pivot_value_new_variable (
8338 dict_get_var (d, fs + i)));
8340 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
8341 N_("Rows"), N_("Columns"));
8344 struct ccase **rows = NULL;
8345 size_t allocated_rows = 0;
8349 while ((c = casereader_read (r)) != NULL)
8351 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
8361 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
8362 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
8363 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
8366 if (change != NOTHING_CHANGED)
8368 matrix_mget_commit_var (rows, n_rows, d,
8369 rowtype_, &mget->rowtypes,
8380 if (n_rows >= allocated_rows)
8381 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
8384 if (change == SPLITS_CHANGED)
8390 /* Reset the factor number, if there are factors. */
8394 if (row_has_factors)
8400 else if (change == FACTORS_CHANGED)
8402 if (row_has_factors)
8408 matrix_mget_commit_var (rows, n_rows, d,
8409 rowtype_, &mget->rowtypes,
8421 if (var_dimension->n_leaves)
8422 pivot_table_submit (pt);
8424 pivot_table_unref (pt);
8428 matrix_mget_execute (struct matrix_command *cmd)
8430 struct matrix_mget *mget = &cmd->mget;
8431 struct casereader *r;
8432 struct dictionary *d;
8433 if (matrix_open_casereader (cmd, "MGET", mget->file, mget->encoding,
8434 mget->state->dataset, &r, &d))
8436 matrix_mget_execute__ (cmd, r, d);
8437 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
8444 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
8446 if (!lex_force_id (s->lexer))
8449 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
8451 *varp = matrix_var_create (s, lex_tokss (s->lexer));
8456 static struct matrix_command *
8457 matrix_eigen_parse (struct matrix_state *s)
8459 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8460 *cmd = (struct matrix_command) {
8462 .eigen = { .expr = NULL }
8465 struct matrix_eigen *eigen = &cmd->eigen;
8466 if (!lex_force_match (s->lexer, T_LPAREN))
8468 eigen->expr = matrix_expr_parse (s);
8470 || !lex_force_match (s->lexer, T_COMMA)
8471 || !matrix_parse_dst_var (s, &eigen->evec)
8472 || !lex_force_match (s->lexer, T_COMMA)
8473 || !matrix_parse_dst_var (s, &eigen->eval)
8474 || !lex_force_match (s->lexer, T_RPAREN))
8480 matrix_command_destroy (cmd);
8485 matrix_eigen_execute (struct matrix_command *cmd)
8487 struct matrix_eigen *eigen = &cmd->eigen;
8488 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
8491 if (!is_symmetric (A))
8493 msg_at (SE, cmd->location, _("Argument of EIGEN must be symmetric."));
8494 gsl_matrix_free (A);
8498 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
8499 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
8500 gsl_vector v_eval = to_vector (eval);
8501 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
8502 gsl_eigen_symmv (A, &v_eval, evec, w);
8503 gsl_eigen_symmv_free (w);
8505 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
8507 gsl_matrix_free (A);
8509 gsl_matrix_free (eigen->eval->value);
8510 eigen->eval->value = eval;
8512 gsl_matrix_free (eigen->evec->value);
8513 eigen->evec->value = evec;
8518 static struct matrix_command *
8519 matrix_setdiag_parse (struct matrix_state *s)
8521 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8522 *cmd = (struct matrix_command) {
8523 .type = MCMD_SETDIAG,
8524 .setdiag = { .dst = NULL }
8527 struct matrix_setdiag *setdiag = &cmd->setdiag;
8528 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
8531 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
8534 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
8539 if (!lex_force_match (s->lexer, T_COMMA))
8542 setdiag->expr = matrix_expr_parse (s);
8546 if (!lex_force_match (s->lexer, T_RPAREN))
8552 matrix_command_destroy (cmd);
8557 matrix_setdiag_execute (struct matrix_command *cmd)
8559 struct matrix_setdiag *setdiag = &cmd->setdiag;
8560 gsl_matrix *dst = setdiag->dst->value;
8563 msg_at (SE, cmd->location,
8564 _("SETDIAG destination matrix %s is uninitialized."),
8565 setdiag->dst->name);
8569 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
8573 size_t n = MIN (dst->size1, dst->size2);
8574 if (is_scalar (src))
8576 double d = to_scalar (src);
8577 for (size_t i = 0; i < n; i++)
8578 gsl_matrix_set (dst, i, i, d);
8580 else if (is_vector (src))
8582 gsl_vector v = to_vector (src);
8583 for (size_t i = 0; i < n && i < v.size; i++)
8584 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
8587 msg_at (SE, matrix_expr_location (setdiag->expr),
8588 _("SETDIAG argument 2 must be a scalar or a vector, "
8589 "not a %zu×%zu matrix."),
8590 src->size1, src->size2);
8591 gsl_matrix_free (src);
8596 static struct matrix_command *
8597 matrix_svd_parse (struct matrix_state *s)
8599 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8600 *cmd = (struct matrix_command) {
8602 .svd = { .expr = NULL }
8605 struct matrix_svd *svd = &cmd->svd;
8606 if (!lex_force_match (s->lexer, T_LPAREN))
8608 svd->expr = matrix_expr_parse (s);
8610 || !lex_force_match (s->lexer, T_COMMA)
8611 || !matrix_parse_dst_var (s, &svd->u)
8612 || !lex_force_match (s->lexer, T_COMMA)
8613 || !matrix_parse_dst_var (s, &svd->s)
8614 || !lex_force_match (s->lexer, T_COMMA)
8615 || !matrix_parse_dst_var (s, &svd->v)
8616 || !lex_force_match (s->lexer, T_RPAREN))
8622 matrix_command_destroy (cmd);
8627 matrix_svd_execute (struct matrix_svd *svd)
8629 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8633 if (m->size1 >= m->size2)
8636 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8637 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8638 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8639 gsl_vector *work = gsl_vector_alloc (A->size2);
8640 gsl_linalg_SV_decomp (A, V, &Sv, work);
8641 gsl_vector_free (work);
8643 matrix_var_set (svd->u, A);
8644 matrix_var_set (svd->s, S);
8645 matrix_var_set (svd->v, V);
8649 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8650 gsl_matrix_transpose_memcpy (At, m);
8651 gsl_matrix_free (m);
8653 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8654 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8655 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8656 gsl_vector *work = gsl_vector_alloc (At->size2);
8657 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8658 gsl_vector_free (work);
8660 matrix_var_set (svd->v, At);
8661 matrix_var_set (svd->s, St);
8662 matrix_var_set (svd->u, Vt);
8666 /* The main MATRIX command logic. */
8669 matrix_command_execute (struct matrix_command *cmd)
8674 matrix_compute_execute (cmd);
8678 matrix_print_execute (&cmd->print);
8682 return matrix_do_if_execute (&cmd->do_if);
8685 matrix_loop_execute (&cmd->loop);
8692 matrix_display_execute (&cmd->display);
8696 matrix_release_execute (&cmd->release);
8700 matrix_save_execute (cmd);
8704 matrix_read_execute (cmd);
8708 matrix_write_execute (&cmd->write);
8712 matrix_get_execute (cmd);
8716 matrix_msave_execute (cmd);
8720 matrix_mget_execute (cmd);
8724 matrix_eigen_execute (cmd);
8728 matrix_setdiag_execute (cmd);
8732 matrix_svd_execute (&cmd->svd);
8740 matrix_command_destroy (struct matrix_command *cmd)
8745 msg_location_destroy (cmd->location);
8750 matrix_lvalue_destroy (cmd->compute.lvalue);
8751 matrix_expr_destroy (cmd->compute.rvalue);
8755 matrix_expr_destroy (cmd->print.expression);
8756 free (cmd->print.title);
8757 print_labels_destroy (cmd->print.rlabels);
8758 print_labels_destroy (cmd->print.clabels);
8762 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8764 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8765 matrix_commands_uninit (&cmd->do_if.clauses[i].commands);
8767 free (cmd->do_if.clauses);
8771 matrix_expr_destroy (cmd->loop.start);
8772 matrix_expr_destroy (cmd->loop.end);
8773 matrix_expr_destroy (cmd->loop.increment);
8774 matrix_expr_destroy (cmd->loop.top_condition);
8775 matrix_expr_destroy (cmd->loop.bottom_condition);
8776 matrix_commands_uninit (&cmd->loop.commands);
8786 free (cmd->release.vars);
8790 matrix_expr_destroy (cmd->save.expression);
8794 matrix_lvalue_destroy (cmd->read.dst);
8795 matrix_expr_destroy (cmd->read.size);
8799 matrix_expr_destroy (cmd->write.expression);
8800 free (cmd->write.format);
8804 matrix_lvalue_destroy (cmd->get.dst);
8805 fh_unref (cmd->get.file);
8806 free (cmd->get.encoding);
8807 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8811 matrix_expr_destroy (cmd->msave.expr);
8815 fh_unref (cmd->mget.file);
8816 stringi_set_destroy (&cmd->mget.rowtypes);
8820 matrix_expr_destroy (cmd->eigen.expr);
8824 matrix_expr_destroy (cmd->setdiag.expr);
8828 matrix_expr_destroy (cmd->svd.expr);
8835 matrix_commands_parse (struct matrix_state *s, struct matrix_commands *c,
8836 const char *command_name,
8837 const char *stop1, const char *stop2)
8839 lex_end_of_command (s->lexer);
8840 lex_discard_rest_of_command (s->lexer);
8842 size_t allocated = 0;
8845 while (lex_token (s->lexer) == T_ENDCMD)
8848 if (lex_at_phrase (s->lexer, stop1)
8849 || (stop2 && lex_at_phrase (s->lexer, stop2)))
8852 if (lex_at_phrase (s->lexer, "END MATRIX"))
8854 msg (SE, _("Premature END MATRIX within %s."), command_name);
8858 struct matrix_command *cmd = matrix_command_parse (s);
8862 if (c->n >= allocated)
8863 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
8864 c->commands[c->n++] = cmd;
8869 matrix_commands_uninit (struct matrix_commands *cmds)
8871 for (size_t i = 0; i < cmds->n; i++)
8872 matrix_command_destroy (cmds->commands[i]);
8873 free (cmds->commands);
8876 struct matrix_command_name
8879 struct matrix_command *(*parse) (struct matrix_state *);
8882 static const struct matrix_command_name *
8883 matrix_command_name_parse (struct lexer *lexer)
8885 static const struct matrix_command_name commands[] = {
8886 { "COMPUTE", matrix_compute_parse },
8887 { "CALL EIGEN", matrix_eigen_parse },
8888 { "CALL SETDIAG", matrix_setdiag_parse },
8889 { "CALL SVD", matrix_svd_parse },
8890 { "PRINT", matrix_print_parse },
8891 { "DO IF", matrix_do_if_parse },
8892 { "LOOP", matrix_loop_parse },
8893 { "BREAK", matrix_break_parse },
8894 { "READ", matrix_read_parse },
8895 { "WRITE", matrix_write_parse },
8896 { "GET", matrix_get_parse },
8897 { "SAVE", matrix_save_parse },
8898 { "MGET", matrix_mget_parse },
8899 { "MSAVE", matrix_msave_parse },
8900 { "DISPLAY", matrix_display_parse },
8901 { "RELEASE", matrix_release_parse },
8903 static size_t n = sizeof commands / sizeof *commands;
8905 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8906 if (lex_match_phrase (lexer, c->name))
8911 static struct matrix_command *
8912 matrix_command_parse (struct matrix_state *s)
8914 int start_ofs = lex_ofs (s->lexer);
8915 size_t nesting_level = SIZE_MAX;
8917 struct matrix_command *c = NULL;
8918 const struct matrix_command_name *cmd = matrix_command_name_parse (s->lexer);
8920 lex_error (s->lexer, _("Unknown matrix command."));
8921 else if (!cmd->parse)
8922 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8926 nesting_level = output_open_group (
8927 group_item_create_nocopy (utf8_to_title (cmd->name),
8928 utf8_to_title (cmd->name)));
8934 c->location = lex_ofs_location (s->lexer, start_ofs, lex_ofs (s->lexer));
8935 msg_location_remove_columns (c->location);
8936 lex_end_of_command (s->lexer);
8938 lex_discard_rest_of_command (s->lexer);
8939 if (nesting_level != SIZE_MAX)
8940 output_close_groups (nesting_level);
8946 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8948 if (!lex_force_match (lexer, T_ENDCMD))
8951 struct matrix_state state = {
8953 .session = dataset_session (ds),
8955 .vars = HMAP_INITIALIZER (state.vars),
8960 while (lex_match (lexer, T_ENDCMD))
8962 if (lex_token (lexer) == T_STOP)
8964 msg (SE, _("Unexpected end of input expecting matrix command."));
8968 if (lex_match_phrase (lexer, "END MATRIX"))
8971 struct matrix_command *c = matrix_command_parse (&state);
8974 matrix_command_execute (c);
8975 matrix_command_destroy (c);
8979 struct matrix_var *var, *next;
8980 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8983 gsl_matrix_free (var->value);
8984 hmap_delete (&state.vars, &var->hmap_node);
8987 hmap_destroy (&state.vars);
8988 msave_common_destroy (state.msave_common);
8989 fh_unref (state.prev_read_file);
8990 for (size_t i = 0; i < state.n_read_files; i++)
8991 read_file_destroy (state.read_files[i]);
8992 free (state.read_files);
8993 fh_unref (state.prev_write_file);
8994 for (size_t i = 0; i < state.n_write_files; i++)
8995 write_file_destroy (state.write_files[i]);
8996 free (state.write_files);
8997 fh_unref (state.prev_save_file);
8998 for (size_t i = 0; i < state.n_save_files; i++)
8999 save_file_destroy (state.save_files[i]);
9000 free (state.save_files);