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)
1855 invert_matrix (gsl_matrix *x)
1857 gsl_permutation *p = gsl_permutation_alloc (x->size1);
1859 gsl_linalg_LU_decomp (x, p, &signum);
1860 gsl_linalg_LU_invx (x, p);
1861 gsl_permutation_free (p);
1865 matrix_eval_INV (gsl_matrix *m)
1872 matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b)
1874 gsl_matrix *k = gsl_matrix_alloc (a->size1 * b->size1,
1875 a->size2 * b->size2);
1877 for (size_t ar = 0; ar < a->size1; ar++)
1878 for (size_t br = 0; br < b->size1; br++, y++)
1881 for (size_t ac = 0; ac < a->size2; ac++)
1882 for (size_t bc = 0; bc < b->size2; bc++, x++)
1884 double av = gsl_matrix_get (a, ar, ac);
1885 double bv = gsl_matrix_get (b, br, bc);
1886 gsl_matrix_set (k, y, x, av * bv);
1893 matrix_eval_LG10 (double d)
1899 matrix_eval_LN (double d)
1905 matrix_eval_MAGIC_odd (gsl_matrix *m, size_t n)
1907 /* Siamese method: https://en.wikipedia.org/wiki/Siamese_method. */
1910 for (size_t i = 1; i <= n * n; i++)
1912 gsl_matrix_set (m, y, x, i);
1914 size_t y1 = !y ? n - 1 : y - 1;
1915 size_t x1 = x + 1 >= n ? 0 : x + 1;
1916 if (gsl_matrix_get (m, y1, x1) == 0)
1922 y = y + 1 >= n ? 0 : y + 1;
1927 magic_exchange (gsl_matrix *m, size_t y1, size_t x1, size_t y2, size_t x2)
1929 double a = gsl_matrix_get (m, y1, x1);
1930 double b = gsl_matrix_get (m, y2, x2);
1931 gsl_matrix_set (m, y1, x1, b);
1932 gsl_matrix_set (m, y2, x2, a);
1936 matrix_eval_MAGIC_doubly_even (gsl_matrix *m, size_t n)
1940 /* A. Umar, "On the Construction of Even Order Magic Squares",
1941 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1943 for (size_t i = 1; i <= n * n / 2; i++)
1945 gsl_matrix_set (m, y, x, i);
1955 for (size_t i = n * n; i > n * n / 2; i--)
1957 gsl_matrix_set (m, y, x, i);
1965 for (size_t y = 0; y < n; y++)
1966 for (size_t x = 0; x < n / 2; x++)
1968 unsigned int d = gsl_matrix_get (m, y, x);
1969 if (d % 2 != (y < n / 2))
1970 magic_exchange (m, y, x, y, n - x - 1);
1975 size_t x1 = n / 2 - 1;
1977 magic_exchange (m, y1, x1, y2, x1);
1978 magic_exchange (m, y1, x2, y2, x2);
1982 matrix_eval_MAGIC_singly_even (gsl_matrix *m, size_t n)
1984 /* A. Umar, "On the Construction of Even Order Magic Squares",
1985 https://arxiv.org/ftp/arxiv/papers/1202/1202.0948.pdf. */
1989 for (size_t i = 1; ; i++)
1991 gsl_matrix_set (m, y, x, i);
1992 if (++y == n / 2 - 1)
2004 for (size_t i = n * n; ; i--)
2006 gsl_matrix_set (m, y, x, i);
2007 if (++y == n / 2 - 1)
2016 for (size_t y = 0; y < n; y++)
2017 if (y != n / 2 - 1 && y != n / 2)
2018 for (size_t x = 0; x < n / 2; x++)
2020 unsigned int d = gsl_matrix_get (m, y, x);
2021 if (d % 2 != (y < n / 2))
2022 magic_exchange (m, y, x, y, n - x - 1);
2025 size_t a0 = (n * n - 2 * n) / 2 + 1;
2026 for (size_t i = 0; i < n / 2; i++)
2029 gsl_matrix_set (m, n / 2, i, a);
2030 gsl_matrix_set (m, n / 2 - 1, i, (n * n + 1) - a);
2032 for (size_t i = 0; i < n / 2; i++)
2034 size_t a = a0 + i + n / 2;
2035 gsl_matrix_set (m, n / 2 - 1, n - i - 1, a);
2036 gsl_matrix_set (m, n / 2, n - i - 1, (n * n + 1) - a);
2038 for (size_t x = 1; x < n / 2; x += 2)
2039 magic_exchange (m, n / 2, x, n / 2 - 1, x);
2040 for (size_t x = n / 2 + 2; x <= n - 3; x += 2)
2041 magic_exchange (m, n / 2, x, n / 2 - 1, x);
2042 size_t x1 = n / 2 - 2;
2043 size_t x2 = n / 2 + 1;
2044 size_t y1 = n / 2 - 2;
2045 size_t y2 = n / 2 + 1;
2046 magic_exchange (m, y1, x1, y2, x1);
2047 magic_exchange (m, y1, x2, y2, x2);
2051 matrix_eval_MAGIC (double n_)
2055 gsl_matrix *m = gsl_matrix_calloc (n, n);
2057 matrix_eval_MAGIC_odd (m, n);
2059 matrix_eval_MAGIC_singly_even (m, n);
2061 matrix_eval_MAGIC_doubly_even (m, n);
2066 matrix_eval_MAKE (double r, double c, double value)
2068 gsl_matrix *m = gsl_matrix_alloc (r, c);
2069 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2075 matrix_eval_MDIAG (gsl_vector *v)
2077 gsl_matrix *m = gsl_matrix_calloc (v->size, v->size);
2078 gsl_vector diagonal = gsl_matrix_diagonal (m).vector;
2079 gsl_vector_memcpy (&diagonal, v);
2084 matrix_eval_MMAX (gsl_matrix *m)
2086 return gsl_matrix_max (m);
2090 matrix_eval_MMIN (gsl_matrix *m)
2092 return gsl_matrix_min (m);
2096 matrix_eval_MOD (gsl_matrix *m, double divisor)
2098 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2099 *d = fmod (*d, divisor);
2104 matrix_eval_MSSQ (gsl_matrix *m)
2107 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2113 matrix_eval_MSUM (gsl_matrix *m)
2116 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2122 matrix_eval_NCOL (gsl_matrix *m)
2128 matrix_eval_NROW (gsl_matrix *m)
2134 matrix_eval_RANK (gsl_matrix *m)
2136 gsl_vector *tau = gsl_vector_alloc (MIN (m->size1, m->size2));
2137 gsl_linalg_QR_decomp (m, tau);
2138 gsl_vector_free (tau);
2140 return gsl_linalg_QRPT_rank (m, -1);
2144 matrix_eval_RESHAPE (gsl_matrix *m, double r_, double c_,
2145 const struct matrix_expr *e)
2147 bool r_ok = r_ >= 0 && r_ < SIZE_MAX;
2148 bool c_ok = c_ >= 0 && c_ < SIZE_MAX;
2152 !r_ok ? e->subs[1]->location : e->subs[2]->location,
2153 _("Arguments 2 and 3 to RESHAPE must be integers."));
2158 if (size_overflow_p (xtimes (r, xmax (c, 1))) || c * r != m->size1 * m->size2)
2160 struct msg_location *loc = msg_location_dup (e->subs[1]->location);
2161 loc->end = e->subs[2]->location->end;
2162 msg_at (SE, loc, _("Product of RESHAPE size arguments (%zu×%zu = %zu) "
2163 "differs from product of matrix dimensions "
2164 "(%zu×%zu = %zu)."),
2166 m->size1, m->size2, m->size1 * m->size2);
2167 msg_location_destroy (loc);
2171 gsl_matrix *dst = gsl_matrix_alloc (r, c);
2174 MATRIX_FOR_ALL_ELEMENTS (d, y2, x2, m)
2176 gsl_matrix_set (dst, y1, x1, *d);
2187 matrix_eval_row_extremum (gsl_matrix *m, bool min)
2192 return gsl_matrix_alloc (0, 1);
2194 gsl_matrix *rext = gsl_matrix_alloc (m->size1, 1);
2195 for (size_t y = 0; y < m->size1; y++)
2197 double ext = gsl_matrix_get (m, y, 0);
2198 for (size_t x = 1; x < m->size2; x++)
2200 double value = gsl_matrix_get (m, y, x);
2201 if (min ? value < ext : value > ext)
2204 gsl_matrix_set (rext, y, 0, ext);
2210 matrix_eval_RMAX (gsl_matrix *m)
2212 return matrix_eval_row_extremum (m, false);
2216 matrix_eval_RMIN (gsl_matrix *m)
2218 return matrix_eval_row_extremum (m, true);
2222 matrix_eval_RND (double d)
2234 rank_compare_3way (const void *a_, const void *b_)
2236 const struct rank *a = a_;
2237 const struct rank *b = b_;
2239 return a->value < b->value ? -1 : a->value > b->value;
2243 matrix_eval_RNKORDER (gsl_matrix *m)
2245 size_t n = m->size1 * m->size2;
2246 struct rank *ranks = xmalloc (n * sizeof *ranks);
2248 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2249 ranks[i++] = (struct rank) { .y = y, .x = x, .value = *d };
2250 qsort (ranks, n, sizeof *ranks, rank_compare_3way);
2252 for (size_t i = 0; i < n; )
2255 for (j = i + 1; j < n; j++)
2256 if (ranks[i].value != ranks[j].value)
2259 double rank = (i + j + 1.0) / 2.0;
2260 for (size_t k = i; k < j; k++)
2261 gsl_matrix_set (m, ranks[k].y, ranks[k].x, rank);
2272 matrix_eval_row_sum (gsl_matrix *m, bool square)
2277 return gsl_matrix_alloc (0, 1);
2279 gsl_matrix *result = gsl_matrix_alloc (m->size1, 1);
2280 for (size_t y = 0; y < m->size1; y++)
2283 for (size_t x = 0; x < m->size2; x++)
2285 double d = gsl_matrix_get (m, y, x);
2286 sum += square ? pow2 (d) : d;
2288 gsl_matrix_set (result, y, 0, sum);
2294 matrix_eval_RSSQ (gsl_matrix *m)
2296 return matrix_eval_row_sum (m, true);
2300 matrix_eval_RSUM (gsl_matrix *m)
2302 return matrix_eval_row_sum (m, false);
2306 matrix_eval_SIN (double d)
2312 matrix_eval_SOLVE (gsl_matrix *m1, gsl_matrix *m2, const struct matrix_expr *e)
2314 if (m1->size1 != m2->size1)
2316 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
2317 loc->end = e->subs[1]->location->end;
2319 msg_at (SE, e->location,
2320 _("SOLVE arguments must have the same number of rows."));
2321 msg_at (SN, e->subs[0]->location,
2322 _("Argument 1 has dimensions %zu×%zu."), m1->size1, m1->size2);
2323 msg_at (SN, e->subs[1]->location,
2324 _("Argument 2 has dimensions %zu×%zu."), m2->size1, m2->size2);
2326 msg_location_destroy (loc);
2330 gsl_matrix *x = gsl_matrix_alloc (m2->size1, m2->size2);
2331 gsl_permutation *p = gsl_permutation_alloc (m1->size1);
2333 gsl_linalg_LU_decomp (m1, p, &signum);
2334 for (size_t i = 0; i < m2->size2; i++)
2336 gsl_vector bi = gsl_matrix_column (m2, i).vector;
2337 gsl_vector xi = gsl_matrix_column (x, i).vector;
2338 gsl_linalg_LU_solve (m1, p, &bi, &xi);
2340 gsl_permutation_free (p);
2345 matrix_eval_SQRT (double d)
2351 matrix_eval_SSCP (gsl_matrix *m)
2353 gsl_matrix *sscp = gsl_matrix_alloc (m->size2, m->size2);
2354 gsl_blas_dgemm (CblasTrans, CblasNoTrans, 1.0, m, m, 0.0, sscp);
2359 matrix_eval_SVAL (gsl_matrix *m)
2361 gsl_matrix *tmp_mat = NULL;
2362 if (m->size2 > m->size1)
2364 tmp_mat = gsl_matrix_alloc (m->size2, m->size1);
2365 gsl_matrix_transpose_memcpy (tmp_mat, m);
2370 gsl_matrix *V = gsl_matrix_alloc (m->size2, m->size2);
2371 gsl_vector *S = gsl_vector_alloc (m->size2);
2372 gsl_vector *work = gsl_vector_alloc (m->size2);
2373 gsl_linalg_SV_decomp (m, V, S, work);
2375 gsl_matrix *vals = gsl_matrix_alloc (m->size2, 1);
2376 for (size_t i = 0; i < m->size2; i++)
2377 gsl_matrix_set (vals, i, 0, gsl_vector_get (S, i));
2379 gsl_matrix_free (V);
2380 gsl_vector_free (S);
2381 gsl_vector_free (work);
2382 gsl_matrix_free (tmp_mat);
2388 matrix_eval_SWEEP (gsl_matrix *m, double d, const struct matrix_expr *e)
2390 if (d < 1 || d > SIZE_MAX)
2392 msg_at (SE, e->subs[1]->location,
2393 _("Scalar argument to SWEEP must be integer."));
2397 if (k >= MIN (m->size1, m->size2))
2399 msg_at (SE, e->subs[1]->location,
2400 _("Scalar argument to SWEEP must be integer less than or "
2401 "equal to the smaller of the matrix argument's rows and "
2406 double m_kk = gsl_matrix_get (m, k, k);
2407 if (fabs (m_kk) > 1e-19)
2409 gsl_matrix *a = gsl_matrix_alloc (m->size1, m->size2);
2410 MATRIX_FOR_ALL_ELEMENTS (a_ij, i, j, a)
2412 double m_ij = gsl_matrix_get (m, i, j);
2413 double m_ik = gsl_matrix_get (m, i, k);
2414 double m_kj = gsl_matrix_get (m, k, j);
2415 *a_ij = (i != k && j != k ? m_ij * m_kk - m_ik * m_kj
2424 for (size_t i = 0; i < m->size1; i++)
2426 gsl_matrix_set (m, i, k, 0);
2427 gsl_matrix_set (m, k, i, 0);
2434 matrix_eval_TRACE (gsl_matrix *m)
2437 size_t n = MIN (m->size1, m->size2);
2438 for (size_t i = 0; i < n; i++)
2439 sum += gsl_matrix_get (m, i, i);
2444 matrix_eval_T (gsl_matrix *m)
2446 return matrix_eval_TRANSPOS (m);
2450 matrix_eval_TRANSPOS (gsl_matrix *m)
2452 if (m->size1 == m->size2)
2454 gsl_matrix_transpose (m);
2459 gsl_matrix *t = gsl_matrix_alloc (m->size2, m->size1);
2460 gsl_matrix_transpose_memcpy (t, m);
2466 matrix_eval_TRUNC (double d)
2472 matrix_eval_UNIFORM (double r_, double c_, const struct matrix_expr *e)
2476 if (size_overflow_p (xtimes (r, xmax (c, 1))))
2478 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
2479 loc->end = e->subs[1]->location->end;
2482 _("Product of arguments to UNIFORM exceeds memory size."));
2484 msg_location_destroy (loc);
2488 gsl_matrix *m = gsl_matrix_alloc (r, c);
2489 MATRIX_FOR_ALL_ELEMENTS (d, y, x, m)
2490 *d = gsl_ran_flat (get_rng (), 0, 1);
2495 matrix_eval_PDF_BETA (double x, double a, double b)
2497 return gsl_ran_beta_pdf (x, a, b);
2501 matrix_eval_CDF_BETA (double x, double a, double b)
2503 return gsl_cdf_beta_P (x, a, b);
2507 matrix_eval_IDF_BETA (double P, double a, double b)
2509 return gsl_cdf_beta_Pinv (P, a, b);
2513 matrix_eval_RV_BETA (double a, double b)
2515 return gsl_ran_beta (get_rng (), a, b);
2519 matrix_eval_NCDF_BETA (double x, double a, double b, double lambda)
2521 return ncdf_beta (x, a, b, lambda);
2525 matrix_eval_NPDF_BETA (double x, double a, double b, double lambda)
2527 return npdf_beta (x, a, b, lambda);
2531 matrix_eval_CDF_BVNOR (double x0, double x1, double r)
2533 return cdf_bvnor (x0, x1, r);
2537 matrix_eval_PDF_BVNOR (double x0, double x1, double r)
2539 return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r);
2543 matrix_eval_CDF_CAUCHY (double x, double a, double b)
2545 return gsl_cdf_cauchy_P ((x - a) / b, 1);
2549 matrix_eval_IDF_CAUCHY (double P, double a, double b)
2551 return a + b * gsl_cdf_cauchy_Pinv (P, 1);
2555 matrix_eval_PDF_CAUCHY (double x, double a, double b)
2557 return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b;
2561 matrix_eval_RV_CAUCHY (double a, double b)
2563 return a + b * gsl_ran_cauchy (get_rng (), 1);
2567 matrix_eval_CDF_CHISQ (double x, double df)
2569 return gsl_cdf_chisq_P (x, df);
2573 matrix_eval_CHICDF (double x, double df)
2575 return matrix_eval_CDF_CHISQ (x, df);
2579 matrix_eval_IDF_CHISQ (double P, double df)
2581 return gsl_cdf_chisq_Pinv (P, df);
2585 matrix_eval_PDF_CHISQ (double x, double df)
2587 return gsl_ran_chisq_pdf (x, df);
2591 matrix_eval_RV_CHISQ (double df)
2593 return gsl_ran_chisq (get_rng (), df);
2597 matrix_eval_SIG_CHISQ (double x, double df)
2599 return gsl_cdf_chisq_Q (x, df);
2603 matrix_eval_CDF_EXP (double x, double a)
2605 return gsl_cdf_exponential_P (x, 1. / a);
2609 matrix_eval_IDF_EXP (double P, double a)
2611 return gsl_cdf_exponential_Pinv (P, 1. / a);
2615 matrix_eval_PDF_EXP (double x, double a)
2617 return gsl_ran_exponential_pdf (x, 1. / a);
2621 matrix_eval_RV_EXP (double a)
2623 return gsl_ran_exponential (get_rng (), 1. / a);
2627 matrix_eval_PDF_XPOWER (double x, double a, double b)
2629 return gsl_ran_exppow_pdf (x, a, b);
2633 matrix_eval_RV_XPOWER (double a, double b)
2635 return gsl_ran_exppow (get_rng (), a, b);
2639 matrix_eval_CDF_F (double x, double df1, double df2)
2641 return gsl_cdf_fdist_P (x, df1, df2);
2645 matrix_eval_FCDF (double x, double df1, double df2)
2647 return matrix_eval_CDF_F (x, df1, df2);
2651 matrix_eval_IDF_F (double P, double df1, double df2)
2653 return idf_fdist (P, df1, df2);
2657 matrix_eval_RV_F (double df1, double df2)
2659 return gsl_ran_fdist (get_rng (), df1, df2);
2663 matrix_eval_PDF_F (double x, double df1, double df2)
2665 return gsl_ran_fdist_pdf (x, df1, df2);
2669 matrix_eval_SIG_F (double x, double df1, double df2)
2671 return gsl_cdf_fdist_Q (x, df1, df2);
2675 matrix_eval_CDF_GAMMA (double x, double a, double b)
2677 return gsl_cdf_gamma_P (x, a, 1. / b);
2681 matrix_eval_IDF_GAMMA (double P, double a, double b)
2683 return gsl_cdf_gamma_Pinv (P, a, 1. / b);
2687 matrix_eval_PDF_GAMMA (double x, double a, double b)
2689 return gsl_ran_gamma_pdf (x, a, 1. / b);
2693 matrix_eval_RV_GAMMA (double a, double b)
2695 return gsl_ran_gamma (get_rng (), a, 1. / b);
2699 matrix_eval_PDF_LANDAU (double x)
2701 return gsl_ran_landau_pdf (x);
2705 matrix_eval_RV_LANDAU (void)
2707 return gsl_ran_landau (get_rng ());
2711 matrix_eval_CDF_LAPLACE (double x, double a, double b)
2713 return gsl_cdf_laplace_P ((x - a) / b, 1);
2717 matrix_eval_IDF_LAPLACE (double P, double a, double b)
2719 return a + b * gsl_cdf_laplace_Pinv (P, 1);
2723 matrix_eval_PDF_LAPLACE (double x, double a, double b)
2725 return gsl_ran_laplace_pdf ((x - a) / b, 1);
2729 matrix_eval_RV_LAPLACE (double a, double b)
2731 return a + b * gsl_ran_laplace (get_rng (), 1);
2735 matrix_eval_RV_LEVY (double c, double alpha)
2737 return gsl_ran_levy (get_rng (), c, alpha);
2741 matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
2743 return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
2747 matrix_eval_CDF_LOGISTIC (double x, double a, double b)
2749 return gsl_cdf_logistic_P ((x - a) / b, 1);
2753 matrix_eval_IDF_LOGISTIC (double P, double a, double b)
2755 return a + b * gsl_cdf_logistic_Pinv (P, 1);
2759 matrix_eval_PDF_LOGISTIC (double x, double a, double b)
2761 return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
2765 matrix_eval_RV_LOGISTIC (double a, double b)
2767 return a + b * gsl_ran_logistic (get_rng (), 1);
2771 matrix_eval_CDF_LNORMAL (double x, double m, double s)
2773 return gsl_cdf_lognormal_P (x, log (m), s);
2777 matrix_eval_IDF_LNORMAL (double P, double m, double s)
2779 return gsl_cdf_lognormal_Pinv (P, log (m), s);;
2783 matrix_eval_PDF_LNORMAL (double x, double m, double s)
2785 return gsl_ran_lognormal_pdf (x, log (m), s);
2789 matrix_eval_RV_LNORMAL (double m, double s)
2791 return gsl_ran_lognormal (get_rng (), log (m), s);
2795 matrix_eval_CDF_NORMAL (double x, double u, double s)
2797 return gsl_cdf_gaussian_P (x - u, s);
2801 matrix_eval_IDF_NORMAL (double P, double u, double s)
2803 return u + gsl_cdf_gaussian_Pinv (P, s);
2807 matrix_eval_PDF_NORMAL (double x, double u, double s)
2809 return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
2813 matrix_eval_RV_NORMAL (double u, double s)
2815 return u + gsl_ran_gaussian (get_rng (), s);
2819 matrix_eval_CDFNORM (double x)
2821 return gsl_cdf_ugaussian_P (x);
2825 matrix_eval_PROBIT (double P)
2827 return gsl_cdf_ugaussian_Pinv (P);
2831 matrix_eval_NORMAL (double s)
2833 return gsl_ran_gaussian (get_rng (), s);
2837 matrix_eval_PDF_NTAIL (double x, double a, double sigma)
2839 return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
2843 matrix_eval_RV_NTAIL (double a, double sigma)
2845 return gsl_ran_gaussian_tail (get_rng (), a, sigma);
2849 matrix_eval_CDF_PARETO (double x, double a, double b)
2851 return gsl_cdf_pareto_P (x, b, a);
2855 matrix_eval_IDF_PARETO (double P, double a, double b)
2857 return gsl_cdf_pareto_Pinv (P, b, a);
2861 matrix_eval_PDF_PARETO (double x, double a, double b)
2863 return gsl_ran_pareto_pdf (x, b, a);
2867 matrix_eval_RV_PARETO (double a, double b)
2869 return gsl_ran_pareto (get_rng (), b, a);
2873 matrix_eval_CDF_RAYLEIGH (double x, double sigma)
2875 return gsl_cdf_rayleigh_P (x, sigma);
2879 matrix_eval_IDF_RAYLEIGH (double P, double sigma)
2881 return gsl_cdf_rayleigh_Pinv (P, sigma);
2885 matrix_eval_PDF_RAYLEIGH (double x, double sigma)
2887 return gsl_ran_rayleigh_pdf (x, sigma);
2891 matrix_eval_RV_RAYLEIGH (double sigma)
2893 return gsl_ran_rayleigh (get_rng (), sigma);
2897 matrix_eval_PDF_RTAIL (double x, double a, double sigma)
2899 return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
2903 matrix_eval_RV_RTAIL (double a, double sigma)
2905 return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
2909 matrix_eval_CDF_T (double x, double df)
2911 return gsl_cdf_tdist_P (x, df);
2915 matrix_eval_TCDF (double x, double df)
2917 return matrix_eval_CDF_T (x, df);
2921 matrix_eval_IDF_T (double P, double df)
2923 return gsl_cdf_tdist_Pinv (P, df);
2927 matrix_eval_PDF_T (double x, double df)
2929 return gsl_ran_tdist_pdf (x, df);
2933 matrix_eval_RV_T (double df)
2935 return gsl_ran_tdist (get_rng (), df);
2939 matrix_eval_CDF_T1G (double x, double a, double b)
2941 return gsl_cdf_gumbel1_P (x, a, b);
2945 matrix_eval_IDF_T1G (double P, double a, double b)
2947 return gsl_cdf_gumbel1_Pinv (P, a, b);
2951 matrix_eval_PDF_T1G (double x, double a, double b)
2953 return gsl_ran_gumbel1_pdf (x, a, b);
2957 matrix_eval_RV_T1G (double a, double b)
2959 return gsl_ran_gumbel1 (get_rng (), a, b);
2963 matrix_eval_CDF_T2G (double x, double a, double b)
2965 return gsl_cdf_gumbel1_P (x, a, b);
2969 matrix_eval_IDF_T2G (double P, double a, double b)
2971 return gsl_cdf_gumbel1_Pinv (P, a, b);
2975 matrix_eval_PDF_T2G (double x, double a, double b)
2977 return gsl_ran_gumbel1_pdf (x, a, b);
2981 matrix_eval_RV_T2G (double a, double b)
2983 return gsl_ran_gumbel1 (get_rng (), a, b);
2987 matrix_eval_CDF_UNIFORM (double x, double a, double b)
2989 return gsl_cdf_flat_P (x, a, b);
2993 matrix_eval_IDF_UNIFORM (double P, double a, double b)
2995 return gsl_cdf_flat_Pinv (P, a, b);
2999 matrix_eval_PDF_UNIFORM (double x, double a, double b)
3001 return gsl_ran_flat_pdf (x, a, b);
3005 matrix_eval_RV_UNIFORM (double a, double b)
3007 return gsl_ran_flat (get_rng (), a, b);
3011 matrix_eval_CDF_WEIBULL (double x, double a, double b)
3013 return gsl_cdf_weibull_P (x, a, b);
3017 matrix_eval_IDF_WEIBULL (double P, double a, double b)
3019 return gsl_cdf_weibull_Pinv (P, a, b);
3023 matrix_eval_PDF_WEIBULL (double x, double a, double b)
3025 return gsl_ran_weibull_pdf (x, a, b);
3029 matrix_eval_RV_WEIBULL (double a, double b)
3031 return gsl_ran_weibull (get_rng (), a, b);
3035 matrix_eval_CDF_BERNOULLI (double k, double p)
3037 return k ? 1 : 1 - p;
3041 matrix_eval_PDF_BERNOULLI (double k, double p)
3043 return gsl_ran_bernoulli_pdf (k, p);
3047 matrix_eval_RV_BERNOULLI (double p)
3049 return gsl_ran_bernoulli (get_rng (), p);
3053 matrix_eval_CDF_BINOM (double k, double n, double p)
3055 return gsl_cdf_binomial_P (k, p, n);
3059 matrix_eval_PDF_BINOM (double k, double n, double p)
3061 return gsl_ran_binomial_pdf (k, p, n);
3065 matrix_eval_RV_BINOM (double n, double p)
3067 return gsl_ran_binomial (get_rng (), p, n);
3071 matrix_eval_CDF_GEOM (double k, double p)
3073 return gsl_cdf_geometric_P (k, p);
3077 matrix_eval_PDF_GEOM (double k, double p)
3079 return gsl_ran_geometric_pdf (k, p);
3083 matrix_eval_RV_GEOM (double p)
3085 return gsl_ran_geometric (get_rng (), p);
3089 matrix_eval_CDF_HYPER (double k, double a, double b, double c)
3091 return gsl_cdf_hypergeometric_P (k, c, a - c, b);
3095 matrix_eval_PDF_HYPER (double k, double a, double b, double c)
3097 return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
3101 matrix_eval_RV_HYPER (double a, double b, double c)
3103 return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
3107 matrix_eval_PDF_LOG (double k, double p)
3109 return gsl_ran_logarithmic_pdf (k, p);
3113 matrix_eval_RV_LOG (double p)
3115 return gsl_ran_logarithmic (get_rng (), p);
3119 matrix_eval_CDF_NEGBIN (double k, double n, double p)
3121 return gsl_cdf_negative_binomial_P (k, p, n);
3125 matrix_eval_PDF_NEGBIN (double k, double n, double p)
3127 return gsl_ran_negative_binomial_pdf (k, p, n);
3131 matrix_eval_RV_NEGBIN (double n, double p)
3133 return gsl_ran_negative_binomial (get_rng (), p, n);
3137 matrix_eval_CDF_POISSON (double k, double mu)
3139 return gsl_cdf_poisson_P (k, mu);
3143 matrix_eval_PDF_POISSON (double k, double mu)
3145 return gsl_ran_poisson_pdf (k, mu);
3149 matrix_eval_RV_POISSON (double mu)
3151 return gsl_ran_poisson (get_rng (), mu);
3155 matrix_op_eval (enum matrix_op op, double a, double b)
3159 case MOP_ADD_ELEMS: return a + b;
3160 case MOP_SUB_ELEMS: return a - b;
3161 case MOP_MUL_ELEMS: return a * b;
3162 case MOP_DIV_ELEMS: return a / b;
3163 case MOP_EXP_ELEMS: return pow (a, b);
3164 case MOP_GT: return a > b;
3165 case MOP_GE: return a >= b;
3166 case MOP_LT: return a < b;
3167 case MOP_LE: return a <= b;
3168 case MOP_EQ: return a == b;
3169 case MOP_NE: return a != b;
3170 case MOP_AND: return (a > 0) && (b > 0);
3171 case MOP_OR: return (a > 0) || (b > 0);
3172 case MOP_XOR: return (a > 0) != (b > 0);
3174 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3183 case MOP_PASTE_HORZ:
3184 case MOP_PASTE_VERT:
3200 matrix_op_name (enum matrix_op op)
3204 case MOP_ADD_ELEMS: return "+";
3205 case MOP_SUB_ELEMS: return "-";
3206 case MOP_MUL_ELEMS: return "&*";
3207 case MOP_DIV_ELEMS: return "&/";
3208 case MOP_EXP_ELEMS: return "&**";
3209 case MOP_GT: return ">";
3210 case MOP_GE: return ">=";
3211 case MOP_LT: return "<";
3212 case MOP_LE: return "<=";
3213 case MOP_EQ: return "=";
3214 case MOP_NE: return "<>";
3215 case MOP_AND: return "AND";
3216 case MOP_OR: return "OR";
3217 case MOP_XOR: return "XOR";
3219 #define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM:
3228 case MOP_PASTE_HORZ:
3229 case MOP_PASTE_VERT:
3245 is_scalar (const gsl_matrix *m)
3247 return m->size1 == 1 && m->size2 == 1;
3251 to_scalar (const gsl_matrix *m)
3253 assert (is_scalar (m));
3254 return gsl_matrix_get (m, 0, 0);
3258 matrix_expr_evaluate_elementwise (const struct matrix_expr *e,
3260 gsl_matrix *a, gsl_matrix *b)
3264 double be = to_scalar (b);
3265 for (size_t r = 0; r < a->size1; r++)
3266 for (size_t c = 0; c < a->size2; c++)
3268 double *ae = gsl_matrix_ptr (a, r, c);
3269 *ae = matrix_op_eval (op, *ae, be);
3273 else if (is_scalar (a))
3275 double ae = to_scalar (a);
3276 for (size_t r = 0; r < b->size1; r++)
3277 for (size_t c = 0; c < b->size2; c++)
3279 double *be = gsl_matrix_ptr (b, r, c);
3280 *be = matrix_op_eval (op, ae, *be);
3284 else if (a->size1 == b->size1 && a->size2 == b->size2)
3286 for (size_t r = 0; r < a->size1; r++)
3287 for (size_t c = 0; c < a->size2; c++)
3289 double *ae = gsl_matrix_ptr (a, r, c);
3290 double be = gsl_matrix_get (b, r, c);
3291 *ae = matrix_op_eval (op, *ae, be);
3297 msg_at (SE, matrix_expr_location (e),
3298 _("The operands of %s must have the same dimensions or one "
3299 "must be a scalar."),
3300 matrix_op_name (op));
3301 msg_at (SN, matrix_expr_location (e->subs[0]),
3302 _("The left-hand operand is a %zu×%zu matrix."),
3303 a->size1, a->size2);
3304 msg_at (SN, matrix_expr_location (e->subs[1]),
3305 _("The right-hand operand is a %zu×%zu matrix."),
3306 b->size1, b->size2);
3312 matrix_expr_evaluate_mul_mat (const struct matrix_expr *e,
3313 gsl_matrix *a, gsl_matrix *b)
3315 if (is_scalar (a) || is_scalar (b))
3316 return matrix_expr_evaluate_elementwise (e, MOP_MUL_ELEMS, a, b);
3318 if (a->size2 != b->size1)
3320 msg_at (SE, e->location,
3321 _("Matrices not conformable for multiplication."));
3322 msg_at (SN, matrix_expr_location (e->subs[0]),
3323 _("The left-hand operand is a %zu×%zu matrix."),
3324 a->size1, a->size2);
3325 msg_at (SN, matrix_expr_location (e->subs[1]),
3326 _("The right-hand operand is a %zu×%zu matrix."),
3327 b->size1, b->size2);
3331 gsl_matrix *c = gsl_matrix_alloc (a->size1, b->size2);
3332 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, a, b, 0.0, c);
3337 swap_matrix (gsl_matrix **a, gsl_matrix **b)
3339 gsl_matrix *tmp = *a;
3345 mul_matrix (gsl_matrix **z, const gsl_matrix *x, const gsl_matrix *y,
3348 gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, x, y, 0.0, *tmp);
3349 swap_matrix (z, tmp);
3353 square_matrix (gsl_matrix **x, gsl_matrix **tmp)
3355 mul_matrix (x, *x, *x, tmp);
3359 matrix_expr_evaluate_exp_mat (const struct matrix_expr *e,
3360 gsl_matrix *x_, gsl_matrix *b)
3363 if (x->size1 != x->size2)
3365 msg_at (SE, matrix_expr_location (e->subs[0]),
3366 _("Matrix exponentation with ** requires a square matrix on "
3367 "the left-hand size, not one with dimensions %zu×%zu."),
3368 x->size1, x->size2);
3373 msg_at (SE, matrix_expr_location (e->subs[1]),
3374 _("Matrix exponentiation with ** requires a scalar on the "
3375 "right-hand side, not a matrix with dimensions %zu×%zu."),
3376 b->size1, b->size2);
3379 double bf = to_scalar (b);
3380 if (bf != floor (bf) || bf <= LONG_MIN || bf > LONG_MAX)
3382 msg_at (SE, matrix_expr_location (e->subs[1]),
3383 _("Exponent %.1f in matrix multiplication is non-integer "
3384 "or outside the valid range."), bf);
3389 gsl_matrix *y_ = gsl_matrix_alloc (x->size1, x->size2);
3391 gsl_matrix_set_identity (y);
3395 gsl_matrix *t_ = gsl_matrix_alloc (x->size1, x->size2);
3397 for (unsigned long int n = labs (bl); n > 1; n /= 2)
3400 mul_matrix (&y, x, y, &t);
3401 square_matrix (&x, &t);
3404 square_matrix (&x, &t);
3406 mul_matrix (&y, x, y, &t);
3410 /* Garbage collection.
3412 There are three matrices: 'x_', 'y_', and 't_', and 'x', 'y', and 't' are
3413 a permutation of them. We are returning one of them; that one must not be
3414 destroyed. We must not destroy 'x_' because the caller owns it. */
3416 gsl_matrix_free (y_);
3418 gsl_matrix_free (t_);
3424 note_operand_size (const gsl_matrix *m, const struct matrix_expr *e)
3426 msg_at (SN, matrix_expr_location (e),
3427 _("This operand is a %zu×%zu matrix."), m->size1, m->size2);
3431 note_nonscalar (const gsl_matrix *m, const struct matrix_expr *e)
3434 note_operand_size (m, e);
3438 matrix_expr_evaluate_seq (const struct matrix_expr *e,
3439 gsl_matrix *start_, gsl_matrix *end_,
3442 if (!is_scalar (start_) || !is_scalar (end_) || (by_ && !is_scalar (by_)))
3444 msg_at (SE, matrix_expr_location (e),
3445 _("All operands of : operator must be scalars."));
3447 note_nonscalar (start_, e->subs[0]);
3448 note_nonscalar (end_, e->subs[1]);
3450 note_nonscalar (by_, e->subs[2]);
3454 long int start = to_scalar (start_);
3455 long int end = to_scalar (end_);
3456 long int by = by_ ? to_scalar (by_) : 1;
3460 msg_at (SE, matrix_expr_location (e->subs[2]),
3461 _("The increment operand to : must be nonzero."));
3465 long int n = (end >= start && by > 0 ? (end - start + by) / by
3466 : end <= start && by < 0 ? (start - end - by) / -by
3468 gsl_matrix *m = gsl_matrix_alloc (1, n);
3469 for (long int i = 0; i < n; i++)
3470 gsl_matrix_set (m, 0, i, start + i * by);
3475 matrix_expr_evaluate_not (gsl_matrix *a)
3477 MATRIX_FOR_ALL_ELEMENTS (d, y, x, a)
3483 matrix_expr_evaluate_paste_horz (const struct matrix_expr *e,
3484 gsl_matrix *a, gsl_matrix *b)
3486 if (a->size1 != b->size1)
3488 if (!a->size1 || !a->size2)
3490 else if (!b->size1 || !b->size2)
3493 msg_at (SE, matrix_expr_location (e),
3494 _("This expression tries to horizontally join matrices with "
3495 "differing numbers of rows."));
3496 note_operand_size (a, e->subs[0]);
3497 note_operand_size (b, e->subs[1]);
3501 gsl_matrix *c = gsl_matrix_alloc (a->size1, a->size2 + b->size2);
3502 for (size_t y = 0; y < a->size1; y++)
3504 for (size_t x = 0; x < a->size2; x++)
3505 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3506 for (size_t x = 0; x < b->size2; x++)
3507 gsl_matrix_set (c, y, x + a->size2, gsl_matrix_get (b, y, x));
3513 matrix_expr_evaluate_paste_vert (const struct matrix_expr *e,
3514 gsl_matrix *a, gsl_matrix *b)
3516 if (a->size2 != b->size2)
3518 if (!a->size1 || !a->size2)
3520 else if (!b->size1 || !b->size2)
3523 msg_at (SE, matrix_expr_location (e),
3524 _("This expression tries to vertically join matrices with "
3525 "differing numbers of columns."));
3526 note_operand_size (a, e->subs[0]);
3527 note_operand_size (b, e->subs[1]);
3531 gsl_matrix *c = gsl_matrix_alloc (a->size1 + b->size1, a->size2);
3532 for (size_t x = 0; x < a->size2; x++)
3534 for (size_t y = 0; y < a->size1; y++)
3535 gsl_matrix_set (c, y, x, gsl_matrix_get (a, y, x));
3536 for (size_t y = 0; y < b->size1; y++)
3537 gsl_matrix_set (c, y + a->size1, x, gsl_matrix_get (b, y, x));
3543 matrix_to_vector (gsl_matrix *m)
3546 gsl_vector v = to_vector (m);
3547 assert (v.block == m->block || !v.block);
3551 gsl_matrix_free (m);
3552 return xmemdup (&v, sizeof v);
3566 #define INDEX_VECTOR_INIT (struct index_vector) { .n = 0 }
3569 index_vector_uninit (struct index_vector *iv)
3576 matrix_normalize_index_vector (const gsl_matrix *m,
3577 const struct matrix_expr *me, size_t size,
3578 enum index_type index_type, size_t other_size,
3579 struct index_vector *iv)
3588 msg_at (SE, matrix_expr_location (me),
3589 _("Vector index must be scalar or vector, not a "
3591 m->size1, m->size2);
3595 msg_at (SE, matrix_expr_location (me),
3596 _("Matrix row index must be scalar or vector, not a "
3598 m->size1, m->size2);
3602 msg_at (SE, matrix_expr_location (me),
3603 _("Matrix column index must be scalar or vector, not a "
3605 m->size1, m->size2);
3611 gsl_vector v = to_vector (CONST_CAST (gsl_matrix *, m));
3612 *iv = (struct index_vector) {
3613 .indexes = xnmalloc (v.size, sizeof *iv->indexes),
3616 for (size_t i = 0; i < v.size; i++)
3618 double index = gsl_vector_get (&v, i);
3619 if (index < 1 || index >= size + 1)
3624 msg_at (SE, matrix_expr_location (me),
3625 _("Index %g is out of range for vector "
3626 "with %zu elements."), index, size);
3630 msg_at (SE, matrix_expr_location (me),
3631 _("%g is not a valid row index for "
3632 "a %zu×%zu matrix."),
3633 index, size, other_size);
3637 msg_at (SE, matrix_expr_location (me),
3638 _("%g is not a valid column index for "
3639 "a %zu×%zu matrix."),
3640 index, other_size, size);
3644 index_vector_uninit (iv);
3647 iv->indexes[i] = index - 1;
3653 *iv = (struct index_vector) {
3654 .indexes = xnmalloc (size, sizeof *iv->indexes),
3657 for (size_t i = 0; i < size; i++)
3664 matrix_expr_evaluate_vec_all (const struct matrix_expr *e,
3667 if (!is_vector (sm))
3669 msg_at (SE, matrix_expr_location (e->subs[0]),
3670 _("Vector index operator may not be applied to "
3671 "a %zu×%zu matrix."),
3672 sm->size1, sm->size2);
3680 matrix_expr_evaluate_vec_index (const struct matrix_expr *e,
3681 gsl_matrix *sm, gsl_matrix *im)
3683 if (!matrix_expr_evaluate_vec_all (e, sm))
3686 gsl_vector sv = to_vector (sm);
3687 struct index_vector iv;
3688 if (!matrix_normalize_index_vector (im, e->subs[1],
3689 sv.size, IV_VECTOR, 0, &iv))
3692 gsl_matrix *dm = gsl_matrix_alloc (sm->size1 == 1 ? 1 : iv.n,
3693 sm->size1 == 1 ? iv.n : 1);
3694 gsl_vector dv = to_vector (dm);
3695 for (size_t dx = 0; dx < iv.n; dx++)
3697 size_t sx = iv.indexes[dx];
3698 gsl_vector_set (&dv, dx, gsl_vector_get (&sv, sx));
3700 index_vector_uninit (&iv);
3706 matrix_expr_evaluate_mat_index (gsl_matrix *sm,
3707 gsl_matrix *im0, const struct matrix_expr *eim0,
3708 gsl_matrix *im1, const struct matrix_expr *eim1)
3710 struct index_vector iv0;
3711 if (!matrix_normalize_index_vector (im0, eim0, sm->size1,
3712 IV_ROW, sm->size2, &iv0))
3715 struct index_vector iv1;
3716 if (!matrix_normalize_index_vector (im1, eim1, sm->size2,
3717 IV_COLUMN, sm->size1, &iv1))
3719 index_vector_uninit (&iv0);
3723 gsl_matrix *dm = gsl_matrix_alloc (iv0.n, iv1.n);
3724 for (size_t dy = 0; dy < iv0.n; dy++)
3726 size_t sy = iv0.indexes[dy];
3728 for (size_t dx = 0; dx < iv1.n; dx++)
3730 size_t sx = iv1.indexes[dx];
3731 gsl_matrix_set (dm, dy, dx, gsl_matrix_get (sm, sy, sx));
3734 index_vector_uninit (&iv0);
3735 index_vector_uninit (&iv1);
3739 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
3740 static gsl_matrix *matrix_expr_evaluate_##PROTO ( \
3741 const struct matrix_function_properties *, gsl_matrix *[], \
3742 const struct matrix_expr *, matrix_proto_##PROTO *);
3747 check_scalar_arg (const char *name, gsl_matrix *subs[],
3748 const struct matrix_expr *e, size_t index)
3750 if (!is_scalar (subs[index]))
3752 msg_at (SE, matrix_expr_location (e->subs[index]),
3753 _("Function %s argument %zu must be a scalar, "
3754 "not a %zu×%zu matrix."),
3755 name, index + 1, subs[index]->size1, subs[index]->size2);
3762 check_vector_arg (const char *name, gsl_matrix *subs[],
3763 const struct matrix_expr *e, size_t index)
3765 if (!is_vector (subs[index]))
3767 msg_at (SE, matrix_expr_location (e->subs[index]),
3768 _("Function %s argument %zu must be a vector, "
3769 "not a %zu×%zu matrix."),
3770 name, index + 1, subs[index]->size1, subs[index]->size2);
3777 to_scalar_args (const char *name, gsl_matrix *subs[],
3778 const struct matrix_expr *e, double d[])
3780 for (size_t i = 0; i < e->n_subs; i++)
3782 if (!check_scalar_arg (name, subs, e, i))
3784 d[i] = to_scalar (subs[i]);
3790 parse_constraint_value (const char **constraintsp)
3793 long retval = strtol (*constraintsp, &tail, 10);
3794 assert (tail > *constraintsp);
3795 *constraintsp = tail;
3799 enum matrix_argument_relop
3809 argument_inequality_error (
3810 const struct matrix_function_properties *props, const struct matrix_expr *e,
3811 size_t ai, gsl_matrix *a, size_t y, size_t x,
3812 size_t bi, double b,
3813 enum matrix_argument_relop relop)
3815 const struct msg_location *loc = matrix_expr_location (e);
3819 msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
3820 "than or equal to argument %zu."),
3821 ai + 1, props->name, bi + 1);
3825 msg_at (ME, loc, _("Argument %zu to matrix function %s must be greater "
3826 "than argument %zu."),
3827 ai + 1, props->name, bi + 1);
3831 msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
3832 "or equal to argument %zu."),
3833 ai + 1, props->name, bi + 1);
3837 msg_at (ME, loc, _("Argument %zu to matrix function %s must be less than "
3839 ai + 1, props->name, bi + 1);
3843 msg_at (ME, loc, _("Argument %zu to matrix function %s must not be equal "
3844 "to argument %zu."),
3845 ai + 1, props->name, bi + 1);
3849 const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
3851 msg_at (SN, a_loc, _("Argument %zu is %g."),
3852 ai + 1, gsl_matrix_get (a, y, x));
3854 msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
3855 y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
3857 msg_at (SN, matrix_expr_location (e->subs[bi]),
3858 _("Argument %zu is %g."), bi + 1, b);
3862 argument_value_error (
3863 const struct matrix_function_properties *props, const struct matrix_expr *e,
3864 size_t ai, gsl_matrix *a, size_t y, size_t x,
3866 enum matrix_argument_relop relop)
3868 const struct msg_location *loc = matrix_expr_location (e);
3872 msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
3873 "than or equal to %g."),
3874 ai + 1, props->name, b);
3878 msg_at (SE, loc, _("Argument %zu to matrix function %s must be greater "
3880 ai + 1, props->name, b);
3884 msg_at (SE, loc, _("Argument %zu to matrix function %s must be less than "
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 not be equal "
3898 ai + 1, props->name, b);
3902 const struct msg_location *a_loc = matrix_expr_location (e->subs[ai]);
3905 if (relop != MRR_NE)
3906 msg_at (SN, a_loc, _("Argument %zu is %g."),
3907 ai + 1, gsl_matrix_get (a, y, x));
3910 msg_at (SN, a_loc, _("Row %zu, column %zu of argument %zu is %g."),
3911 y + 1, x + 1, ai + 1, gsl_matrix_get (a, y, x));
3915 matrix_argument_relop_is_satisfied (double a, double b,
3916 enum matrix_argument_relop relop)
3920 case MRR_GE: return a >= b;
3921 case MRR_GT: return a > b;
3922 case MRR_LE: return a <= b;
3923 case MRR_LT: return a < b;
3924 case MRR_NE: return a != b;
3930 static enum matrix_argument_relop
3931 matrix_argument_relop_flip (enum matrix_argument_relop relop)
3935 case MRR_GE: return MRR_LE;
3936 case MRR_GT: return MRR_LT;
3937 case MRR_LE: return MRR_GE;
3938 case MRR_LT: return MRR_GT;
3939 case MRR_NE: return MRR_NE;
3946 check_constraints (const struct matrix_function_properties *props,
3947 gsl_matrix *args[], const struct matrix_expr *e)
3949 size_t n_args = e->n_subs;
3950 const char *constraints = props->constraints;
3954 size_t arg_index = SIZE_MAX;
3955 while (*constraints)
3957 if (*constraints >= 'a' && *constraints <= 'd')
3959 arg_index = *constraints++ - 'a';
3960 assert (arg_index < n_args);
3962 else if (*constraints == '[' || *constraints == '(')
3964 assert (arg_index < n_args);
3965 bool open_lower = *constraints++ == '(';
3966 int minimum = parse_constraint_value (&constraints);
3967 assert (*constraints == ',');
3969 int maximum = parse_constraint_value (&constraints);
3970 assert (*constraints == ']' || *constraints == ')');
3971 bool open_upper = *constraints++ == ')';
3973 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
3974 if ((open_lower ? *d <= minimum : *d < minimum)
3975 || (open_upper ? *d >= maximum : *d > maximum))
3977 if (!is_scalar (args[arg_index]))
3978 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
3979 _("Row %zu, column %zu of argument %zu to matrix "
3980 "function %s is %g, which is outside "
3981 "the valid range %c%d,%d%c."),
3982 y + 1, x + 1, arg_index + 1, props->name, *d,
3983 open_lower ? '(' : '[',
3985 open_upper ? ')' : ']');
3987 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
3988 _("Argument %zu to matrix function %s is %g, "
3989 "which is outside the valid range %c%d,%d%c."),
3990 arg_index + 1, props->name, *d,
3991 open_lower ? '(' : '[',
3993 open_upper ? ')' : ']');
3997 else if (*constraints == 'i')
4000 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
4001 if (*d != floor (*d))
4003 if (!is_scalar (args[arg_index]))
4004 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
4005 _("Argument %zu to matrix function %s, which must be "
4006 "integer, contains non-integer value %g in "
4007 "row %zu, column %zu."),
4008 arg_index + 1, props->name, *d, y + 1, x + 1);
4010 msg_at (SE, matrix_expr_location (e->subs[arg_index]),
4011 _("Argument %zu to matrix function %s, which must be "
4012 "integer, has non-integer value %g."),
4013 arg_index + 1, props->name, *d);
4017 else if (*constraints == '>'
4018 || *constraints == '<'
4019 || *constraints == '!')
4021 enum matrix_argument_relop relop;
4022 switch (*constraints++)
4025 if (*constraints == '=')
4035 if (*constraints == '=')
4045 assert (*constraints == '=');
4054 if (*constraints >= 'a' && *constraints <= 'd')
4056 size_t a_index = arg_index;
4057 size_t b_index = *constraints - 'a';
4058 assert (a_index < n_args);
4059 assert (b_index < n_args);
4061 /* We only support one of the two arguments being non-scalar.
4062 It's easier to support only the first one being non-scalar, so
4063 flip things around if it's the other way. */
4064 if (!is_scalar (args[b_index]))
4066 assert (is_scalar (args[a_index]));
4067 size_t tmp_index = a_index;
4069 b_index = tmp_index;
4070 relop = matrix_argument_relop_flip (relop);
4073 double b = to_scalar (args[b_index]);
4074 MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
4075 if (!matrix_argument_relop_is_satisfied (*a, b, relop))
4077 argument_inequality_error (
4079 a_index, args[a_index], y, x,
4087 int comparand = parse_constraint_value (&constraints);
4089 MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
4090 if (!matrix_argument_relop_is_satisfied (*d, comparand, relop))
4092 argument_value_error (
4094 arg_index, args[arg_index], y, x,
4103 assert (*constraints == ' ');
4105 arg_index = SIZE_MAX;
4112 matrix_expr_evaluate_d_none (const struct matrix_function_properties *props,
4113 gsl_matrix *subs[], const struct matrix_expr *e,
4114 matrix_proto_d_none *f)
4116 assert (e->n_subs == 0);
4118 if (!check_constraints (props, subs, e))
4121 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4122 gsl_matrix_set (m, 0, 0, f ());
4127 matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
4128 gsl_matrix *subs[], const struct matrix_expr *e,
4129 matrix_proto_d_d *f)
4131 assert (e->n_subs == 1);
4134 if (!to_scalar_args (props->name, subs, e, &d)
4135 || !check_constraints (props, subs, e))
4138 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4139 gsl_matrix_set (m, 0, 0, f (d));
4144 matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props,
4145 gsl_matrix *subs[], const struct matrix_expr *e,
4146 matrix_proto_d_dd *f)
4148 assert (e->n_subs == 2);
4151 if (!to_scalar_args (props->name, subs, e, d)
4152 && !check_constraints (props, subs, e))
4155 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4156 gsl_matrix_set (m, 0, 0, f (d[0], d[1]));
4161 matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
4162 gsl_matrix *subs[], const struct matrix_expr *e,
4163 matrix_proto_d_ddd *f)
4165 assert (e->n_subs == 3);
4168 if (!to_scalar_args (props->name, subs, e, d)
4169 || !check_constraints (props, subs, e))
4172 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4173 gsl_matrix_set (m, 0, 0, f (d[0], d[1], d[2]));
4178 matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
4179 gsl_matrix *subs[], const struct matrix_expr *e,
4180 matrix_proto_m_d *f)
4182 assert (e->n_subs == 1);
4185 return (to_scalar_args (props->name, subs, e, &d)
4186 && check_constraints (props, subs, e)
4192 matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props,
4193 gsl_matrix *subs[], const struct matrix_expr *e,
4194 matrix_proto_m_ddd *f)
4196 assert (e->n_subs == 3);
4199 return (to_scalar_args (props->name, subs, e, d)
4200 && check_constraints (props, subs, e)
4201 ? f(d[0], d[1], d[2])
4206 matrix_expr_evaluate_m_ddn (const struct matrix_function_properties *props,
4207 gsl_matrix *subs[], const struct matrix_expr *e,
4208 matrix_proto_m_ddn *f)
4210 assert (e->n_subs == 2);
4213 return (to_scalar_args (props->name, subs, e, d)
4214 && check_constraints (props, subs, e)
4220 matrix_expr_evaluate_m_m (const struct matrix_function_properties *props,
4221 gsl_matrix *subs[], const struct matrix_expr *e,
4222 matrix_proto_m_m *f)
4224 assert (e->n_subs == 1);
4225 return check_constraints (props, subs, e) ? f (subs[0]) : NULL;
4229 matrix_expr_evaluate_m_mn (const struct matrix_function_properties *props,
4230 gsl_matrix *subs[], const struct matrix_expr *e,
4231 matrix_proto_m_mn *f)
4233 assert (e->n_subs == 1);
4234 return check_constraints (props, subs, e) ? f (subs[0], e) : NULL;
4238 matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
4239 gsl_matrix *subs[], const struct matrix_expr *e,
4240 matrix_proto_m_e *f)
4242 assert (e->n_subs == 1);
4244 if (!check_constraints (props, subs, e))
4247 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4253 matrix_expr_evaluate_m_md (const struct matrix_function_properties *props,
4254 gsl_matrix *subs[], const struct matrix_expr *e,
4255 matrix_proto_m_md *f)
4257 assert (e->n_subs == 2);
4258 return (check_scalar_arg (props->name, subs, e, 1)
4259 && check_constraints (props, subs, e)
4260 ? f (subs[0], to_scalar (subs[1]))
4265 matrix_expr_evaluate_m_mdn (const struct matrix_function_properties *props,
4266 gsl_matrix *subs[], const struct matrix_expr *e,
4267 matrix_proto_m_mdn *f)
4269 assert (e->n_subs == 2);
4270 return (check_scalar_arg (props->name, subs, e, 1)
4271 && check_constraints (props, subs, e)
4272 ? f (subs[0], to_scalar (subs[1]), e)
4277 matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
4278 gsl_matrix *subs[], const struct matrix_expr *e,
4279 matrix_proto_m_ed *f)
4281 assert (e->n_subs == 2);
4282 if (!check_scalar_arg (props->name, subs, e, 1)
4283 || !check_constraints (props, subs, e))
4286 double b = to_scalar (subs[1]);
4287 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4293 matrix_expr_evaluate_m_mddn (const struct matrix_function_properties *props,
4294 gsl_matrix *subs[], const struct matrix_expr *e,
4295 matrix_proto_m_mddn *f)
4297 assert (e->n_subs == 3);
4298 if (!check_scalar_arg (props->name, subs, e, 1)
4299 || !check_scalar_arg (props->name, subs, e, 2)
4300 || !check_constraints (props, subs, e))
4302 return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2]), e);
4306 matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
4307 gsl_matrix *subs[], const struct matrix_expr *e,
4308 matrix_proto_m_edd *f)
4310 assert (e->n_subs == 3);
4311 if (!check_scalar_arg (props->name, subs, e, 1)
4312 || !check_scalar_arg (props->name, subs, e, 2)
4313 || !check_constraints (props, subs, e))
4316 double b = to_scalar (subs[1]);
4317 double c = to_scalar (subs[2]);
4318 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4324 matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
4325 gsl_matrix *subs[], const struct matrix_expr *e,
4326 matrix_proto_m_eddd *f)
4328 assert (e->n_subs == 4);
4329 for (size_t i = 1; i < 4; i++)
4330 if (!check_scalar_arg (props->name, subs, e, i))
4333 if (!check_constraints (props, subs, e))
4336 double b = to_scalar (subs[1]);
4337 double c = to_scalar (subs[2]);
4338 double d = to_scalar (subs[3]);
4339 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4340 *a = f (*a, b, c, d);
4345 matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
4346 gsl_matrix *subs[], const struct matrix_expr *e,
4347 matrix_proto_m_eed *f)
4349 assert (e->n_subs == 3);
4350 if (!check_scalar_arg (props->name, subs, e, 2))
4353 if (!is_scalar (subs[0]) && !is_scalar (subs[1])
4354 && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
4356 struct msg_location *loc = msg_location_dup (e->subs[0]->location);
4357 loc->end = e->subs[1]->location->end;
4360 _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
4361 "%zu×%zu, but %s requires these arguments either to have "
4362 "the same dimensions or for one of them to be a scalar."),
4364 subs[0]->size1, subs[0]->size2,
4365 subs[1]->size1, subs[1]->size2,
4368 msg_location_destroy (loc);
4372 if (!check_constraints (props, subs, e))
4375 double c = to_scalar (subs[2]);
4377 if (is_scalar (subs[0]))
4379 double a = to_scalar (subs[0]);
4380 MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
4386 double b = to_scalar (subs[1]);
4387 MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
4394 matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props,
4395 gsl_matrix *subs[], const struct matrix_expr *e,
4396 matrix_proto_m_mm *f)
4398 assert (e->n_subs == 2);
4399 return check_constraints (props, subs, e) ? f (subs[0], subs[1]) : NULL;
4403 matrix_expr_evaluate_m_mmn (const struct matrix_function_properties *props,
4404 gsl_matrix *subs[], const struct matrix_expr *e,
4405 matrix_proto_m_mmn *f)
4407 assert (e->n_subs == 2);
4408 return check_constraints (props, subs, e) ? f (subs[0], subs[1], e) : NULL;
4412 matrix_expr_evaluate_m_v (const struct matrix_function_properties *props,
4413 gsl_matrix *subs[], const struct matrix_expr *e,
4414 matrix_proto_m_v *f)
4416 assert (e->n_subs == 1);
4417 if (!check_vector_arg (props->name, subs, e, 0)
4418 || !check_constraints (props, subs, e))
4420 gsl_vector v = to_vector (subs[0]);
4425 matrix_expr_evaluate_d_m (const struct matrix_function_properties *props,
4426 gsl_matrix *subs[], const struct matrix_expr *e,
4427 matrix_proto_d_m *f)
4429 assert (e->n_subs == 1);
4431 if (!check_constraints (props, subs, e))
4434 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4435 gsl_matrix_set (m, 0, 0, f (subs[0]));
4440 matrix_expr_evaluate_m_any (const struct matrix_function_properties *props,
4441 gsl_matrix *subs[], const struct matrix_expr *e,
4442 matrix_proto_m_any *f)
4444 return check_constraints (props, subs, e) ? f (subs, e->n_subs) : NULL;
4448 matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props_ UNUSED,
4449 gsl_matrix *subs[], const struct matrix_expr *e,
4450 matrix_proto_IDENT *f)
4452 static const struct matrix_function_properties p1 = {
4454 .constraints = "ai>=0"
4456 static const struct matrix_function_properties p2 = {
4458 .constraints = "ai>=0 bi>=0"
4460 const struct matrix_function_properties *props = e->n_subs == 1 ? &p1 : &p2;
4462 assert (e->n_subs <= 2);
4465 return (to_scalar_args (props->name, subs, e, d)
4466 && check_constraints (props, subs, e)
4467 ? f (d[0], d[e->n_subs - 1])
4472 matrix_expr_evaluate (const struct matrix_expr *e)
4474 if (e->op == MOP_NUMBER)
4476 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4477 gsl_matrix_set (m, 0, 0, e->number);
4480 else if (e->op == MOP_VARIABLE)
4482 const gsl_matrix *src = e->variable->value;
4485 msg_at (SE, e->location,
4486 _("Uninitialized variable %s used in expression."),
4491 gsl_matrix *dst = gsl_matrix_alloc (src->size1, src->size2);
4492 gsl_matrix_memcpy (dst, src);
4495 else if (e->op == MOP_EOF)
4497 struct dfm_reader *reader = read_file_open (e->eof);
4498 gsl_matrix *m = gsl_matrix_alloc (1, 1);
4499 gsl_matrix_set (m, 0, 0, !reader || dfm_eof (reader));
4503 enum { N_LOCAL = 3 };
4504 gsl_matrix *local_subs[N_LOCAL];
4505 gsl_matrix **subs = (e->n_subs < N_LOCAL
4507 : xmalloc (e->n_subs * sizeof *subs));
4509 for (size_t i = 0; i < e->n_subs; i++)
4511 subs[i] = matrix_expr_evaluate (e->subs[i]);
4514 for (size_t j = 0; j < i; j++)
4515 gsl_matrix_free (subs[j]);
4516 if (subs != local_subs)
4522 gsl_matrix *result = NULL;
4525 #define F(ENUM, STRING, PROTO, CONSTRAINTS) \
4526 case MOP_F_##ENUM: \
4528 static const struct matrix_function_properties props = { \
4530 .constraints = CONSTRAINTS, \
4532 result = matrix_expr_evaluate_##PROTO (&props, subs, e, \
4533 matrix_eval_##ENUM); \
4540 gsl_matrix_scale (subs[0], -1.0);
4558 result = matrix_expr_evaluate_elementwise (e, e->op, subs[0], subs[1]);
4562 result = matrix_expr_evaluate_not (subs[0]);
4566 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], NULL);
4570 result = matrix_expr_evaluate_seq (e, subs[0], subs[1], subs[2]);
4574 result = matrix_expr_evaluate_mul_mat (e, subs[0], subs[1]);
4578 result = matrix_expr_evaluate_exp_mat (e, subs[0], subs[1]);
4581 case MOP_PASTE_HORZ:
4582 result = matrix_expr_evaluate_paste_horz (e, subs[0], subs[1]);
4585 case MOP_PASTE_VERT:
4586 result = matrix_expr_evaluate_paste_vert (e, subs[0], subs[1]);
4590 result = gsl_matrix_alloc (0, 0);
4594 result = matrix_expr_evaluate_vec_index (e, subs[0], subs[1]);
4598 result = matrix_expr_evaluate_vec_all (e, subs[0]);
4602 result = matrix_expr_evaluate_mat_index (subs[0],
4603 subs[1], e->subs[1],
4604 subs[2], e->subs[2]);
4608 result = matrix_expr_evaluate_mat_index (subs[0],
4609 subs[1], e->subs[1],
4614 result = matrix_expr_evaluate_mat_index (subs[0],
4616 subs[1], e->subs[1]);
4625 for (size_t i = 0; i < e->n_subs; i++)
4626 if (subs[i] != result)
4627 gsl_matrix_free (subs[i]);
4628 if (subs != local_subs)
4634 matrix_expr_evaluate_scalar (const struct matrix_expr *e, const char *context,
4637 gsl_matrix *m = matrix_expr_evaluate (e);
4643 msg_at (SE, matrix_expr_location (e),
4644 _("Expression for %s must evaluate to scalar, "
4645 "not a %zu×%zu matrix."),
4646 context, m->size1, m->size2);
4647 gsl_matrix_free (m);
4652 gsl_matrix_free (m);
4657 matrix_expr_evaluate_integer (const struct matrix_expr *e, const char *context,
4661 if (!matrix_expr_evaluate_scalar (e, context, &d))
4665 if (d < LONG_MIN || d > LONG_MAX)
4667 msg_at (SE, matrix_expr_location (e),
4668 _("Expression for %s is outside the integer range."), context);
4677 An lvalue is an expression that can appear on the left side of a COMPUTE
4678 command and in other contexts that assign values.
4680 An lvalue is parsed once, with matrix_lvalue_parse(). It can then be
4681 evaluated (with matrix_lvalue_evaluate()) and assigned (with
4682 matrix_lvalue_assign()).
4684 There are three kinds of lvalues:
4686 - A variable name. A variable used as an lvalue need not be initialized,
4687 since the assignment will initialize.
4689 - A subvector, e.g. "var(index0)". The variable must be initialized and
4690 must have the form of a vector (it must have 1 column or 1 row).
4692 - A submatrix, e.g. "var(index0, index1)". The variable must be
4694 struct matrix_lvalue
4696 struct matrix_var *var; /* Destination variable. */
4697 struct matrix_expr *indexes[2]; /* Index expressions, if any. */
4698 size_t n_indexes; /* Number of indexes. */
4700 struct msg_location *var_location; /* Variable name. */
4701 struct msg_location *full_location; /* Variable name plus indexing. */
4702 struct msg_location *index_locations[2]; /* Index expressions. */
4707 matrix_lvalue_destroy (struct matrix_lvalue *lvalue)
4711 msg_location_destroy (lvalue->var_location);
4712 msg_location_destroy (lvalue->full_location);
4713 for (size_t i = 0; i < lvalue->n_indexes; i++)
4715 matrix_expr_destroy (lvalue->indexes[i]);
4716 msg_location_destroy (lvalue->index_locations[i]);
4722 /* Parses and returns an lvalue at the current position in S's lexer. Returns
4723 null on parse failure. On success, the caller must eventually free the
4725 static struct matrix_lvalue *
4726 matrix_lvalue_parse (struct matrix_state *s)
4728 if (!lex_force_id (s->lexer))
4731 struct matrix_lvalue *lvalue = xzalloc (sizeof *lvalue);
4732 int start_ofs = lex_ofs (s->lexer);
4733 lvalue->var_location = lex_get_location (s->lexer, 0, 0);
4734 lvalue->var = matrix_var_lookup (s, lex_tokss (s->lexer));
4735 if (lex_next_token (s->lexer, 1) == T_LPAREN)
4739 msg (SE, _("Undefined variable %s."), lex_tokcstr (s->lexer));
4743 lex_get_n (s->lexer, 2);
4745 if (!matrix_parse_index_expr (s, &lvalue->indexes[0],
4746 &lvalue->index_locations[0]))
4748 lvalue->n_indexes++;
4750 if (lex_match (s->lexer, T_COMMA))
4752 if (!matrix_parse_index_expr (s, &lvalue->indexes[1],
4753 &lvalue->index_locations[1]))
4755 lvalue->n_indexes++;
4757 if (!lex_force_match (s->lexer, T_RPAREN))
4760 lvalue->full_location = lex_ofs_location (s->lexer, start_ofs,
4761 lex_ofs (s->lexer) - 1);
4766 lvalue->var = matrix_var_create (s, lex_tokss (s->lexer));
4772 matrix_lvalue_destroy (lvalue);
4777 matrix_lvalue_evaluate_vector (struct matrix_expr *e, size_t size,
4778 enum index_type index_type, size_t other_size,
4779 struct index_vector *iv)
4784 m = matrix_expr_evaluate (e);
4791 bool ok = matrix_normalize_index_vector (m, e, size, index_type,
4793 gsl_matrix_free (m);
4797 /* Evaluates the indexes in LVALUE into IV0 and IV1, owned by the caller.
4798 Returns true if successful, false if evaluating the expressions failed or if
4799 LVALUE otherwise can't be used for an assignment.
4801 On success, the caller retains ownership of the index vectors, which are
4802 suitable for passing to matrix_lvalue_assign(). If not used for that
4803 purpose then they need to eventually be freed (with
4804 index_vector_uninit()). */
4806 matrix_lvalue_evaluate (struct matrix_lvalue *lvalue,
4807 struct index_vector *iv0,
4808 struct index_vector *iv1)
4810 *iv0 = INDEX_VECTOR_INIT;
4811 *iv1 = INDEX_VECTOR_INIT;
4812 if (!lvalue->n_indexes)
4815 /* Validate destination matrix exists and has the right shape. */
4816 gsl_matrix *dm = lvalue->var->value;
4819 msg_at (SE, lvalue->var_location,
4820 _("Undefined variable %s."), lvalue->var->name);
4823 else if (dm->size1 == 0 || dm->size2 == 0)
4825 msg_at (SE, lvalue->full_location, _("Cannot index %zu×%zu matrix %s."),
4826 dm->size1, dm->size2, lvalue->var->name);
4829 else if (lvalue->n_indexes == 1)
4831 if (!is_vector (dm))
4833 msg_at (SE, lvalue->full_location,
4834 _("Can't use vector indexing on %zu×%zu matrix %s."),
4835 dm->size1, dm->size2, lvalue->var->name);
4838 return matrix_lvalue_evaluate_vector (lvalue->indexes[0],
4839 MAX (dm->size1, dm->size2),
4844 assert (lvalue->n_indexes == 2);
4845 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[0], dm->size1,
4846 IV_ROW, dm->size2, iv0))
4849 if (!matrix_lvalue_evaluate_vector (lvalue->indexes[1], dm->size2,
4850 IV_COLUMN, dm->size1, iv1))
4852 index_vector_uninit (iv0);
4860 matrix_lvalue_assign_vector (struct matrix_lvalue *lvalue,
4861 struct index_vector *iv,
4862 gsl_matrix *sm, const struct msg_location *lsm)
4864 /* Convert source matrix 'sm' to source vector 'sv'. */
4865 if (!is_vector (sm))
4867 msg_at (SE, lvalue->full_location,
4868 _("Only an %zu-element vector may be assigned to this "
4869 "%zu-element subvector of %s."),
4870 iv->n, iv->n, lvalue->var->name);
4872 _("The source is an %zu×%zu matrix."),
4873 sm->size1, sm->size2);
4876 gsl_vector sv = to_vector (sm);
4877 if (iv->n != sv.size)
4879 msg_at (SE, lvalue->full_location,
4880 _("Only an %zu-element vector may be assigned to this "
4881 "%zu-element subvector of %s."),
4882 iv->n, iv->n, lvalue->var->name);
4883 msg_at (SE, lsm, ngettext ("The source vector has %zu element.",
4884 "The source vector has %zu elements.",
4890 /* Assign elements. */
4891 gsl_vector dv = to_vector (lvalue->var->value);
4892 for (size_t x = 0; x < iv->n; x++)
4893 gsl_vector_set (&dv, iv->indexes[x], gsl_vector_get (&sv, x));
4898 matrix_lvalue_assign_matrix (struct matrix_lvalue *lvalue,
4899 struct index_vector *iv0,
4900 struct index_vector *iv1,
4901 gsl_matrix *sm, const struct msg_location *lsm)
4903 gsl_matrix *dm = lvalue->var->value;
4905 /* Convert source matrix 'sm' to source vector 'sv'. */
4906 bool wrong_rows = iv0->n != sm->size1;
4907 bool wrong_columns = iv1->n != sm->size2;
4908 if (wrong_rows || wrong_columns)
4910 if (wrong_rows && wrong_columns)
4911 msg_at (SE, lvalue->full_location,
4912 _("Numbers of indexes for assigning to %s differ from the "
4913 "size of the source matrix."),
4915 else if (wrong_rows)
4916 msg_at (SE, lvalue->full_location,
4917 _("Number of row indexes for assigning to %s differs from "
4918 "number of rows in source matrix."),
4921 msg_at (SE, lvalue->full_location,
4922 _("Number of column indexes for assigning to %s differs from "
4923 "number of columns in source matrix."),
4928 if (lvalue->indexes[0])
4929 msg_at (SN, lvalue->index_locations[0],
4930 ngettext ("There is %zu row index.",
4931 "There are %zu row indexes.",
4935 msg_at (SN, lvalue->index_locations[0],
4936 ngettext ("Destination matrix %s has %zu row.",
4937 "Destination matrix %s has %zu rows.",
4939 lvalue->var->name, iv0->n);
4944 if (lvalue->indexes[1])
4945 msg_at (SN, lvalue->index_locations[1],
4946 ngettext ("There is %zu column index.",
4947 "There are %zu column indexes.",
4951 msg_at (SN, lvalue->index_locations[1],
4952 ngettext ("Destination matrix %s has %zu column.",
4953 "Destination matrix %s has %zu columns.",
4955 lvalue->var->name, iv1->n);
4958 msg_at (SN, lsm, _("The source matrix is %zu×%zu."),
4959 sm->size1, sm->size2);
4963 /* Assign elements. */
4964 for (size_t y = 0; y < iv0->n; y++)
4966 size_t f0 = iv0->indexes[y];
4967 for (size_t x = 0; x < iv1->n; x++)
4969 size_t f1 = iv1->indexes[x];
4970 gsl_matrix_set (dm, f0, f1, gsl_matrix_get (sm, y, x));
4976 /* Assigns rvalue SM to LVALUE, whose index vectors IV0 and IV1 were previously
4977 obtained with matrix_lvalue_evaluate(). Returns true if successful, false
4978 on error. Always takes ownership of M. LSM should be the source location
4979 to use for errors related to SM. */
4981 matrix_lvalue_assign (struct matrix_lvalue *lvalue,
4982 struct index_vector *iv0, struct index_vector *iv1,
4983 gsl_matrix *sm, const struct msg_location *lsm)
4985 if (!lvalue->n_indexes)
4987 gsl_matrix_free (lvalue->var->value);
4988 lvalue->var->value = sm;
4993 bool ok = (lvalue->n_indexes == 1
4994 ? matrix_lvalue_assign_vector (lvalue, iv0, sm, lsm)
4995 : matrix_lvalue_assign_matrix (lvalue, iv0, iv1, sm, lsm));
4996 index_vector_uninit (iv0);
4997 index_vector_uninit (iv1);
4998 gsl_matrix_free (sm);
5003 /* Evaluates and then assigns SM to LVALUE. Always takes ownership of M. LSM
5004 should be the source location to use for errors related to SM.*/
5006 matrix_lvalue_evaluate_and_assign (struct matrix_lvalue *lvalue,
5008 const struct msg_location *lsm)
5010 struct index_vector iv0, iv1;
5011 if (!matrix_lvalue_evaluate (lvalue, &iv0, &iv1))
5013 gsl_matrix_free (sm);
5017 return matrix_lvalue_assign (lvalue, &iv0, &iv1, sm, lsm);
5020 /* Matrix command data structure. */
5022 /* An array of matrix commands. */
5023 struct matrix_commands
5025 struct matrix_command **commands;
5029 static bool matrix_commands_parse (struct matrix_state *,
5030 struct matrix_commands *,
5031 const char *command_name,
5032 const char *stop1, const char *stop2);
5033 static void matrix_commands_uninit (struct matrix_commands *);
5035 /* A single matrix command. */
5036 struct matrix_command
5038 /* The type of command. */
5039 enum matrix_command_type
5060 /* Source lines for this command. */
5061 struct msg_location *location;
5065 struct matrix_compute
5067 struct matrix_lvalue *lvalue;
5068 struct matrix_expr *rvalue;
5074 struct matrix_expr *expression;
5075 bool use_default_format;
5076 struct fmt_spec format;
5078 int space; /* -1 means NEWPAGE. */
5082 struct string_array literals; /* CLABELS/RLABELS. */
5083 struct matrix_expr *expr; /* CNAMES/RNAMES. */
5093 struct matrix_expr *condition;
5094 struct matrix_commands commands;
5104 struct matrix_var *var;
5105 struct matrix_expr *start;
5106 struct matrix_expr *end;
5107 struct matrix_expr *increment;
5109 /* Loop conditions. */
5110 struct matrix_expr *top_condition;
5111 struct matrix_expr *bottom_condition;
5114 struct matrix_commands commands;
5118 struct matrix_display
5120 struct matrix_state *state;
5124 struct matrix_release
5126 struct matrix_var **vars;
5133 struct matrix_expr *expression;
5134 struct save_file *sf;
5140 struct read_file *rf;
5141 struct matrix_lvalue *dst;
5142 struct matrix_expr *size;
5144 enum fmt_type format;
5153 struct write_file *wf;
5154 struct matrix_expr *expression;
5157 /* If this is nonnull, WRITE uses this format.
5159 If this is NULL, WRITE uses free-field format with as many
5160 digits of precision as needed. */
5161 struct fmt_spec *format;
5170 struct matrix_lvalue *dst;
5171 struct dataset *dataset;
5172 struct file_handle *file;
5174 struct var_syntax *vars;
5176 struct matrix_var *names;
5178 /* Treatment of missing values. */
5183 MGET_ERROR, /* Flag error on command. */
5184 MGET_ACCEPT, /* Accept without change, user-missing only. */
5185 MGET_OMIT, /* Drop this case. */
5186 MGET_RECODE /* Recode to 'substitute'. */
5189 double substitute; /* MGET_RECODE only. */
5197 struct msave_common *common;
5198 struct matrix_expr *expr;
5199 const char *rowtype;
5200 const struct matrix_expr *factors;
5201 const struct matrix_expr *splits;
5207 struct matrix_state *state;
5208 struct file_handle *file;
5210 struct stringi_set rowtypes;
5216 struct matrix_expr *expr;
5217 struct matrix_var *evec;
5218 struct matrix_var *eval;
5222 struct matrix_setdiag
5224 struct matrix_var *dst;
5225 struct matrix_expr *expr;
5231 struct matrix_expr *expr;
5232 struct matrix_var *u;
5233 struct matrix_var *s;
5234 struct matrix_var *v;
5240 static struct matrix_command *matrix_command_parse (struct matrix_state *);
5241 static bool matrix_command_execute (struct matrix_command *);
5242 static void matrix_command_destroy (struct matrix_command *);
5246 static struct matrix_command *
5247 matrix_compute_parse (struct matrix_state *s)
5249 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5250 *cmd = (struct matrix_command) {
5251 .type = MCMD_COMPUTE,
5252 .compute = { .lvalue = NULL },
5255 struct matrix_compute *compute = &cmd->compute;
5256 compute->lvalue = matrix_lvalue_parse (s);
5257 if (!compute->lvalue)
5260 if (!lex_force_match (s->lexer, T_EQUALS))
5263 compute->rvalue = matrix_expr_parse (s);
5264 if (!compute->rvalue)
5270 matrix_command_destroy (cmd);
5275 matrix_compute_execute (struct matrix_command *cmd)
5277 struct matrix_compute *compute = &cmd->compute;
5278 gsl_matrix *value = matrix_expr_evaluate (compute->rvalue);
5282 matrix_lvalue_evaluate_and_assign (compute->lvalue, value,
5283 matrix_expr_location (compute->rvalue));
5289 print_labels_destroy (struct print_labels *labels)
5293 string_array_destroy (&labels->literals);
5294 matrix_expr_destroy (labels->expr);
5300 parse_literal_print_labels (struct matrix_state *s,
5301 struct print_labels **labelsp)
5303 lex_match (s->lexer, T_EQUALS);
5304 print_labels_destroy (*labelsp);
5305 *labelsp = xzalloc (sizeof **labelsp);
5306 while (lex_token (s->lexer) != T_SLASH
5307 && lex_token (s->lexer) != T_ENDCMD
5308 && lex_token (s->lexer) != T_STOP)
5310 struct string label = DS_EMPTY_INITIALIZER;
5311 while (lex_token (s->lexer) != T_COMMA
5312 && lex_token (s->lexer) != T_SLASH
5313 && lex_token (s->lexer) != T_ENDCMD
5314 && lex_token (s->lexer) != T_STOP)
5316 if (!ds_is_empty (&label))
5317 ds_put_byte (&label, ' ');
5319 if (lex_token (s->lexer) == T_STRING)
5320 ds_put_cstr (&label, lex_tokcstr (s->lexer));
5323 char *rep = lex_next_representation (s->lexer, 0, 0);
5324 ds_put_cstr (&label, rep);
5329 string_array_append_nocopy (&(*labelsp)->literals,
5330 ds_steal_cstr (&label));
5332 if (!lex_match (s->lexer, T_COMMA))
5338 parse_expr_print_labels (struct matrix_state *s, struct print_labels **labelsp)
5340 lex_match (s->lexer, T_EQUALS);
5341 struct matrix_expr *e = matrix_parse_exp (s);
5345 print_labels_destroy (*labelsp);
5346 *labelsp = xzalloc (sizeof **labelsp);
5347 (*labelsp)->expr = e;
5351 static struct matrix_command *
5352 matrix_print_parse (struct matrix_state *s)
5354 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5355 *cmd = (struct matrix_command) {
5358 .use_default_format = true,
5362 if (lex_token (s->lexer) != T_SLASH && lex_token (s->lexer) != T_ENDCMD)
5365 for (size_t i = 0; ; i++)
5367 enum token_type t = lex_next_token (s->lexer, i);
5368 if (t == T_LPAREN || t == T_LBRACK || t == T_LCURLY)
5370 else if ((t == T_RPAREN || t == T_RBRACK || t == T_RCURLY) && depth)
5372 else if ((t == T_SLASH && !depth) || t == T_ENDCMD || t == T_STOP)
5375 cmd->print.title = lex_next_representation (s->lexer, 0, i - 1);
5380 cmd->print.expression = matrix_parse_exp (s);
5381 if (!cmd->print.expression)
5385 while (lex_match (s->lexer, T_SLASH))
5387 if (lex_match_id (s->lexer, "FORMAT"))
5389 lex_match (s->lexer, T_EQUALS);
5390 if (!parse_format_specifier (s->lexer, &cmd->print.format))
5392 cmd->print.use_default_format = false;
5394 else if (lex_match_id (s->lexer, "TITLE"))
5396 lex_match (s->lexer, T_EQUALS);
5397 if (!lex_force_string (s->lexer))
5399 free (cmd->print.title);
5400 cmd->print.title = ss_xstrdup (lex_tokss (s->lexer));
5403 else if (lex_match_id (s->lexer, "SPACE"))
5405 lex_match (s->lexer, T_EQUALS);
5406 if (lex_match_id (s->lexer, "NEWPAGE"))
5407 cmd->print.space = -1;
5408 else if (lex_force_int_range (s->lexer, "SPACE", 1, INT_MAX))
5410 cmd->print.space = lex_integer (s->lexer);
5416 else if (lex_match_id (s->lexer, "RLABELS"))
5417 parse_literal_print_labels (s, &cmd->print.rlabels);
5418 else if (lex_match_id (s->lexer, "CLABELS"))
5419 parse_literal_print_labels (s, &cmd->print.clabels);
5420 else if (lex_match_id (s->lexer, "RNAMES"))
5422 if (!parse_expr_print_labels (s, &cmd->print.rlabels))
5425 else if (lex_match_id (s->lexer, "CNAMES"))
5427 if (!parse_expr_print_labels (s, &cmd->print.clabels))
5432 lex_error_expecting (s->lexer, "FORMAT", "TITLE", "SPACE",
5433 "RLABELS", "CLABELS", "RNAMES", "CNAMES");
5441 matrix_command_destroy (cmd);
5446 matrix_is_integer (const gsl_matrix *m)
5448 for (size_t y = 0; y < m->size1; y++)
5449 for (size_t x = 0; x < m->size2; x++)
5451 double d = gsl_matrix_get (m, y, x);
5459 matrix_max_magnitude (const gsl_matrix *m)
5462 for (size_t y = 0; y < m->size1; y++)
5463 for (size_t x = 0; x < m->size2; x++)
5465 double d = fabs (gsl_matrix_get (m, y, x));
5473 format_fits (struct fmt_spec format, double x)
5475 char *s = data_out (&(union value) { .f = x }, NULL,
5476 &format, settings_get_fmt_settings ());
5477 bool fits = *s != '*' && !strchr (s, 'E');
5482 static struct fmt_spec
5483 default_format (const gsl_matrix *m, int *log_scale)
5487 double max = matrix_max_magnitude (m);
5489 if (matrix_is_integer (m))
5490 for (int w = 1; w <= 12; w++)
5492 struct fmt_spec format = { .type = FMT_F, .w = w };
5493 if (format_fits (format, -max))
5497 if (max >= 1e9 || max <= 1e-4)
5500 snprintf (s, sizeof s, "%.1e", max);
5502 const char *e = strchr (s, 'e');
5504 *log_scale = atoi (e + 1);
5507 return (struct fmt_spec) { .type = FMT_F, .w = 13, .d = 10 };
5511 trimmed_string (double d)
5513 struct substring s = ss_buffer ((char *) &d, sizeof d);
5514 ss_rtrim (&s, ss_cstr (" "));
5515 return ss_xstrdup (s);
5518 static struct string_array *
5519 print_labels_get (const struct print_labels *labels, size_t n,
5520 const char *prefix, bool truncate)
5525 struct string_array *out = xzalloc (sizeof *out);
5526 if (labels->literals.n)
5527 string_array_clone (out, &labels->literals);
5528 else if (labels->expr)
5530 gsl_matrix *m = matrix_expr_evaluate (labels->expr);
5531 if (m && is_vector (m))
5533 gsl_vector v = to_vector (m);
5534 for (size_t i = 0; i < v.size; i++)
5535 string_array_append_nocopy (out, trimmed_string (
5536 gsl_vector_get (&v, i)));
5538 gsl_matrix_free (m);
5544 string_array_append_nocopy (out, xasprintf ("%s%zu", prefix,
5547 string_array_append (out, "");
5550 string_array_delete (out, out->n - 1);
5553 for (size_t i = 0; i < out->n; i++)
5555 char *s = out->strings[i];
5556 s[strnlen (s, 8)] = '\0';
5563 matrix_print_space (int space)
5566 output_item_submit (page_break_item_create ());
5567 for (int i = 0; i < space; i++)
5568 output_log ("%s", "");
5572 matrix_print_text (const struct matrix_print *print, const gsl_matrix *m,
5573 struct fmt_spec format, int log_scale)
5575 matrix_print_space (print->space);
5577 output_log ("%s", print->title);
5579 output_log (" 10 ** %d X", log_scale);
5581 struct string_array *clabels = print_labels_get (print->clabels,
5582 m->size2, "col", true);
5583 if (clabels && format.w < 8)
5585 struct string_array *rlabels = print_labels_get (print->rlabels,
5586 m->size1, "row", true);
5590 struct string line = DS_EMPTY_INITIALIZER;
5592 ds_put_byte_multiple (&line, ' ', 8);
5593 for (size_t x = 0; x < m->size2; x++)
5594 ds_put_format (&line, " %*s", format.w, clabels->strings[x]);
5595 output_log_nocopy (ds_steal_cstr (&line));
5598 double scale = pow (10.0, log_scale);
5599 bool numeric = fmt_is_numeric (format.type);
5600 for (size_t y = 0; y < m->size1; y++)
5602 struct string line = DS_EMPTY_INITIALIZER;
5604 ds_put_format (&line, "%-8s", rlabels->strings[y]);
5606 for (size_t x = 0; x < m->size2; x++)
5608 double f = gsl_matrix_get (m, y, x);
5610 ? data_out (&(union value) { .f = f / scale}, NULL,
5611 &format, settings_get_fmt_settings ())
5612 : trimmed_string (f));
5613 ds_put_format (&line, " %s", s);
5616 output_log_nocopy (ds_steal_cstr (&line));
5619 string_array_destroy (rlabels);
5621 string_array_destroy (clabels);
5626 create_print_dimension (struct pivot_table *table, enum pivot_axis_type axis,
5627 const struct print_labels *print_labels, size_t n,
5628 const char *name, const char *prefix)
5630 struct string_array *labels = print_labels_get (print_labels, n, prefix,
5632 struct pivot_dimension *d = pivot_dimension_create (table, axis, name);
5633 for (size_t i = 0; i < n; i++)
5634 pivot_category_create_leaf (
5636 ? pivot_value_new_user_text (labels->strings[i], SIZE_MAX)
5637 : pivot_value_new_integer (i + 1)));
5639 d->hide_all_labels = true;
5640 string_array_destroy (labels);
5645 matrix_print_tables (const struct matrix_print *print, const gsl_matrix *m,
5646 struct fmt_spec format, int log_scale)
5648 struct pivot_table *table = pivot_table_create__ (
5649 pivot_value_new_user_text (print->title ? print->title : "", SIZE_MAX),
5652 create_print_dimension (table, PIVOT_AXIS_ROW, print->rlabels, m->size1,
5654 create_print_dimension (table, PIVOT_AXIS_COLUMN, print->clabels, m->size2,
5655 N_("Columns"), "col");
5657 struct pivot_footnote *footnote = NULL;
5660 char *s = xasprintf ("× 10**%d", log_scale);
5661 footnote = pivot_table_create_footnote (
5662 table, pivot_value_new_user_text_nocopy (s));
5665 double scale = pow (10.0, log_scale);
5666 bool numeric = fmt_is_numeric (format.type);
5667 for (size_t y = 0; y < m->size1; y++)
5668 for (size_t x = 0; x < m->size2; x++)
5670 double f = gsl_matrix_get (m, y, x);
5671 struct pivot_value *v;
5674 v = pivot_value_new_number (f / scale);
5675 v->numeric.format = format;
5678 v = pivot_value_new_user_text_nocopy (trimmed_string (f));
5680 pivot_value_add_footnote (v, footnote);
5681 pivot_table_put2 (table, y, x, v);
5684 pivot_table_submit (table);
5688 matrix_print_execute (const struct matrix_print *print)
5690 if (print->expression)
5692 gsl_matrix *m = matrix_expr_evaluate (print->expression);
5697 struct fmt_spec format = (print->use_default_format
5698 ? default_format (m, &log_scale)
5701 if (settings_get_mdisplay () == SETTINGS_MDISPLAY_TEXT)
5702 matrix_print_text (print, m, format, log_scale);
5704 matrix_print_tables (print, m, format, log_scale);
5706 gsl_matrix_free (m);
5710 matrix_print_space (print->space);
5712 output_log ("%s", print->title);
5719 matrix_do_if_clause_parse (struct matrix_state *s,
5720 struct matrix_do_if *ifc,
5721 bool parse_condition,
5722 size_t *allocated_clauses)
5724 if (ifc->n_clauses >= *allocated_clauses)
5725 ifc->clauses = x2nrealloc (ifc->clauses, allocated_clauses,
5726 sizeof *ifc->clauses);
5727 struct do_if_clause *c = &ifc->clauses[ifc->n_clauses++];
5728 *c = (struct do_if_clause) { .condition = NULL };
5730 if (parse_condition)
5732 c->condition = matrix_expr_parse (s);
5737 return matrix_commands_parse (s, &c->commands, "DO IF", "ELSE", "END IF");
5740 static struct matrix_command *
5741 matrix_do_if_parse (struct matrix_state *s)
5743 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5744 *cmd = (struct matrix_command) {
5746 .do_if = { .n_clauses = 0 }
5749 size_t allocated_clauses = 0;
5752 if (!matrix_do_if_clause_parse (s, &cmd->do_if, true, &allocated_clauses))
5755 while (lex_match_phrase (s->lexer, "ELSE IF"));
5757 if (lex_match_id (s->lexer, "ELSE")
5758 && !matrix_do_if_clause_parse (s, &cmd->do_if, false, &allocated_clauses))
5761 if (!lex_match_phrase (s->lexer, "END IF"))
5766 matrix_command_destroy (cmd);
5771 matrix_do_if_execute (struct matrix_do_if *cmd)
5773 for (size_t i = 0; i < cmd->n_clauses; i++)
5775 struct do_if_clause *c = &cmd->clauses[i];
5779 if (!matrix_expr_evaluate_scalar (c->condition,
5780 i ? "ELSE IF" : "DO IF",
5785 for (size_t j = 0; j < c->commands.n; j++)
5786 if (!matrix_command_execute (c->commands.commands[j]))
5795 static struct matrix_command *
5796 matrix_loop_parse (struct matrix_state *s)
5798 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5799 *cmd = (struct matrix_command) { .type = MCMD_LOOP, .loop = { .var = NULL } };
5801 struct matrix_loop *loop = &cmd->loop;
5802 if (lex_token (s->lexer) == T_ID && lex_next_token (s->lexer, 1) == T_EQUALS)
5804 struct substring name = lex_tokss (s->lexer);
5805 loop->var = matrix_var_lookup (s, name);
5807 loop->var = matrix_var_create (s, name);
5812 loop->start = matrix_expr_parse (s);
5813 if (!loop->start || !lex_force_match (s->lexer, T_TO))
5816 loop->end = matrix_expr_parse (s);
5820 if (lex_match (s->lexer, T_BY))
5822 loop->increment = matrix_expr_parse (s);
5823 if (!loop->increment)
5828 if (lex_match_id (s->lexer, "IF"))
5830 loop->top_condition = matrix_expr_parse (s);
5831 if (!loop->top_condition)
5835 bool was_in_loop = s->in_loop;
5837 bool ok = matrix_commands_parse (s, &loop->commands, "LOOP",
5839 s->in_loop = was_in_loop;
5843 if (!lex_match_phrase (s->lexer, "END LOOP"))
5846 if (lex_match_id (s->lexer, "IF"))
5848 loop->bottom_condition = matrix_expr_parse (s);
5849 if (!loop->bottom_condition)
5856 matrix_command_destroy (cmd);
5861 matrix_loop_execute (struct matrix_loop *cmd)
5865 long int increment = 1;
5868 if (!matrix_expr_evaluate_integer (cmd->start, "LOOP", &value)
5869 || !matrix_expr_evaluate_integer (cmd->end, "TO", &end)
5871 && !matrix_expr_evaluate_integer (cmd->increment, "BY",
5875 if (increment > 0 ? value > end
5876 : increment < 0 ? value < end
5881 int mxloops = settings_get_mxloops ();
5882 for (int i = 0; i < mxloops; i++)
5887 && (cmd->var->value->size1 != 1 || cmd->var->value->size2 != 1))
5889 gsl_matrix_free (cmd->var->value);
5890 cmd->var->value = NULL;
5892 if (!cmd->var->value)
5893 cmd->var->value = gsl_matrix_alloc (1, 1);
5895 gsl_matrix_set (cmd->var->value, 0, 0, value);
5898 if (cmd->top_condition)
5901 if (!matrix_expr_evaluate_scalar (cmd->top_condition, "LOOP IF",
5906 for (size_t j = 0; j < cmd->commands.n; j++)
5907 if (!matrix_command_execute (cmd->commands.commands[j]))
5910 if (cmd->bottom_condition)
5913 if (!matrix_expr_evaluate_scalar (cmd->bottom_condition,
5914 "END LOOP IF", &d) || d > 0)
5921 if (increment > 0 ? value > end : value < end)
5929 static struct matrix_command *
5930 matrix_break_parse (struct matrix_state *s)
5934 msg (SE, _("BREAK not inside LOOP."));
5938 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5939 *cmd = (struct matrix_command) { .type = MCMD_BREAK };
5945 static struct matrix_command *
5946 matrix_display_parse (struct matrix_state *s)
5948 while (lex_token (s->lexer) != T_ENDCMD)
5950 if (!lex_match_id (s->lexer, "DICTIONARY")
5951 && !lex_match_id (s->lexer, "STATUS"))
5953 lex_error_expecting (s->lexer, "DICTIONARY", "STATUS");
5958 struct matrix_command *cmd = xmalloc (sizeof *cmd);
5959 *cmd = (struct matrix_command) { .type = MCMD_DISPLAY, .display = { s } };
5964 compare_matrix_var_pointers (const void *a_, const void *b_)
5966 const struct matrix_var *const *ap = a_;
5967 const struct matrix_var *const *bp = b_;
5968 const struct matrix_var *a = *ap;
5969 const struct matrix_var *b = *bp;
5970 return strcmp (a->name, b->name);
5974 matrix_display_execute (struct matrix_display *cmd)
5976 const struct matrix_state *s = cmd->state;
5978 struct pivot_table *table = pivot_table_create (N_("Matrix Variables"));
5979 struct pivot_dimension *attr_dimension
5980 = pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Attribute"));
5981 pivot_category_create_group (attr_dimension->root, N_("Dimension"),
5982 N_("Rows"), N_("Columns"));
5983 pivot_category_create_leaves (attr_dimension->root, N_("Size (kB)"));
5985 const struct matrix_var **vars = xmalloc (hmap_count (&s->vars) * sizeof *vars);
5988 const struct matrix_var *var;
5989 HMAP_FOR_EACH (var, struct matrix_var, hmap_node, &s->vars)
5991 vars[n_vars++] = var;
5992 qsort (vars, n_vars, sizeof *vars, compare_matrix_var_pointers);
5994 struct pivot_dimension *rows = pivot_dimension_create (
5995 table, PIVOT_AXIS_ROW, N_("Variable"));
5996 for (size_t i = 0; i < n_vars; i++)
5998 const struct matrix_var *var = vars[i];
5999 pivot_category_create_leaf (
6000 rows->root, pivot_value_new_user_text (var->name, SIZE_MAX));
6002 size_t r = var->value->size1;
6003 size_t c = var->value->size2;
6004 double values[] = { r, c, r * c * sizeof (double) / 1024 };
6005 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
6006 pivot_table_put2 (table, j, i, pivot_value_new_integer (values[j]));
6009 pivot_table_submit (table);
6014 static struct matrix_command *
6015 matrix_release_parse (struct matrix_state *s)
6017 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6018 *cmd = (struct matrix_command) { .type = MCMD_RELEASE };
6020 size_t allocated_vars = 0;
6021 while (lex_token (s->lexer) == T_ID)
6023 struct matrix_var *var = matrix_var_lookup (s, lex_tokss (s->lexer));
6026 if (cmd->release.n_vars >= allocated_vars)
6027 cmd->release.vars = x2nrealloc (cmd->release.vars, &allocated_vars,
6028 sizeof *cmd->release.vars);
6029 cmd->release.vars[cmd->release.n_vars++] = var;
6032 lex_error (s->lexer, _("Variable name expected."));
6035 if (!lex_match (s->lexer, T_COMMA))
6043 matrix_release_execute (struct matrix_release *cmd)
6045 for (size_t i = 0; i < cmd->n_vars; i++)
6047 struct matrix_var *var = cmd->vars[i];
6048 gsl_matrix_free (var->value);
6055 static struct save_file *
6056 save_file_create (struct matrix_state *s, struct file_handle *fh,
6057 struct string_array *variables,
6058 struct matrix_expr *names,
6059 struct stringi_set *strings)
6061 for (size_t i = 0; i < s->n_save_files; i++)
6063 struct save_file *sf = s->save_files[i];
6064 if (fh_equal (sf->file, fh))
6068 string_array_destroy (variables);
6069 matrix_expr_destroy (names);
6070 stringi_set_destroy (strings);
6076 struct save_file *sf = xmalloc (sizeof *sf);
6077 *sf = (struct save_file) {
6079 .dataset = s->dataset,
6080 .variables = *variables,
6082 .strings = STRINGI_SET_INITIALIZER (sf->strings),
6085 stringi_set_swap (&sf->strings, strings);
6086 stringi_set_destroy (strings);
6088 s->save_files = xrealloc (s->save_files,
6089 (s->n_save_files + 1) * sizeof *s->save_files);
6090 s->save_files[s->n_save_files++] = sf;
6094 static struct casewriter *
6095 save_file_open (struct save_file *sf, gsl_matrix *m,
6096 const struct msg_location *save_location)
6098 if (sf->writer || sf->error)
6102 size_t n_variables = caseproto_get_n_widths (
6103 casewriter_get_proto (sf->writer));
6104 if (m->size2 != n_variables)
6106 const char *file_name = (sf->file == fh_inline_file () ? "*"
6107 : fh_get_name (sf->file));
6108 msg_at (SE, save_location,
6109 _("Cannot save %zu×%zu matrix to %s because the "
6110 "first SAVE to %s in this matrix program wrote a "
6111 "%zu-column matrix."),
6112 m->size1, m->size2, file_name, file_name, n_variables);
6113 msg_at (SE, sf->location,
6114 _("This is the location of the first SAVE to %s."),
6123 struct dictionary *dict = dict_create (get_default_encoding ());
6125 /* Fill 'names' with user-specified names if there were any, either from
6126 sf->variables or sf->names. */
6127 struct string_array names = { .n = 0 };
6128 if (sf->variables.n)
6129 string_array_clone (&names, &sf->variables);
6132 gsl_matrix *nm = matrix_expr_evaluate (sf->names);
6133 if (nm && is_vector (nm))
6135 gsl_vector nv = to_vector (nm);
6136 for (size_t i = 0; i < nv.size; i++)
6138 char *name = trimmed_string (gsl_vector_get (&nv, i));
6139 if (dict_id_is_valid (dict, name, true))
6140 string_array_append_nocopy (&names, name);
6145 gsl_matrix_free (nm);
6148 struct stringi_set strings;
6149 stringi_set_clone (&strings, &sf->strings);
6151 for (size_t i = 0; dict_get_n_vars (dict) < m->size2; i++)
6156 name = names.strings[i];
6159 snprintf (tmp_name, sizeof tmp_name, "COL%zu", i + 1);
6163 int width = stringi_set_delete (&strings, name) ? 8 : 0;
6164 struct variable *var = dict_create_var (dict, name, width);
6167 msg_at (ME, save_location,
6168 _("Duplicate variable name %s in SAVE statement."), name);
6173 if (!stringi_set_is_empty (&strings))
6175 size_t count = stringi_set_count (&strings);
6176 const char *example = stringi_set_node_get_string (
6177 stringi_set_first (&strings));
6180 msg_at (ME, save_location,
6181 _("The SAVE command STRINGS subcommand specifies an unknown "
6182 "variable %s."), example);
6184 msg_at (ME, save_location,
6185 ngettext ("The SAVE command STRINGS subcommand specifies %zu "
6186 "unknown variable: %s.",
6187 "The SAVE command STRINGS subcommand specifies %zu "
6188 "unknown variables, including %s.",
6193 stringi_set_destroy (&strings);
6194 string_array_destroy (&names);
6203 if (sf->file == fh_inline_file ())
6204 sf->writer = autopaging_writer_create (dict_get_proto (dict));
6206 sf->writer = any_writer_open (sf->file, dict);
6210 sf->location = msg_location_dup (save_location);
6221 save_file_destroy (struct save_file *sf)
6225 if (sf->file == fh_inline_file () && sf->writer && sf->dict)
6227 dataset_set_dict (sf->dataset, sf->dict);
6228 dataset_set_source (sf->dataset, casewriter_make_reader (sf->writer));
6232 casewriter_destroy (sf->writer);
6233 dict_unref (sf->dict);
6235 fh_unref (sf->file);
6236 string_array_destroy (&sf->variables);
6237 matrix_expr_destroy (sf->names);
6238 stringi_set_destroy (&sf->strings);
6239 msg_location_destroy (sf->location);
6244 static struct matrix_command *
6245 matrix_save_parse (struct matrix_state *s)
6247 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6248 *cmd = (struct matrix_command) {
6250 .save = { .expression = NULL },
6253 struct file_handle *fh = NULL;
6254 struct matrix_save *save = &cmd->save;
6256 struct string_array variables = STRING_ARRAY_INITIALIZER;
6257 struct matrix_expr *names = NULL;
6258 struct stringi_set strings = STRINGI_SET_INITIALIZER (strings);
6260 save->expression = matrix_parse_exp (s);
6261 if (!save->expression)
6264 while (lex_match (s->lexer, T_SLASH))
6266 if (lex_match_id (s->lexer, "OUTFILE"))
6268 lex_match (s->lexer, T_EQUALS);
6271 fh = (lex_match (s->lexer, T_ASTERISK)
6273 : fh_parse (s->lexer, FH_REF_FILE, s->session));
6277 else if (lex_match_id (s->lexer, "VARIABLES"))
6279 lex_match (s->lexer, T_EQUALS);
6283 struct dictionary *d = dict_create (get_default_encoding ());
6284 bool ok = parse_DATA_LIST_vars (s->lexer, d, &names, &n,
6285 PV_NO_SCRATCH | PV_NO_DUPLICATE);
6290 string_array_clear (&variables);
6291 variables = (struct string_array) {
6297 else if (lex_match_id (s->lexer, "NAMES"))
6299 lex_match (s->lexer, T_EQUALS);
6300 matrix_expr_destroy (names);
6301 names = matrix_parse_exp (s);
6305 else if (lex_match_id (s->lexer, "STRINGS"))
6307 lex_match (s->lexer, T_EQUALS);
6308 while (lex_token (s->lexer) == T_ID)
6310 stringi_set_insert (&strings, lex_tokcstr (s->lexer));
6312 if (!lex_match (s->lexer, T_COMMA))
6318 lex_error_expecting (s->lexer, "OUTFILE", "VARIABLES", "NAMES",
6326 if (s->prev_save_file)
6327 fh = fh_ref (s->prev_save_file);
6330 lex_sbc_missing ("OUTFILE");
6334 fh_unref (s->prev_save_file);
6335 s->prev_save_file = fh_ref (fh);
6337 if (variables.n && names)
6339 msg (SW, _("VARIABLES and NAMES both specified; ignoring NAMES."));
6340 matrix_expr_destroy (names);
6344 save->sf = save_file_create (s, fh, &variables, names, &strings);
6348 string_array_destroy (&variables);
6349 matrix_expr_destroy (names);
6350 stringi_set_destroy (&strings);
6352 matrix_command_destroy (cmd);
6357 matrix_save_execute (const struct matrix_command *cmd)
6359 const struct matrix_save *save = &cmd->save;
6361 gsl_matrix *m = matrix_expr_evaluate (save->expression);
6365 struct casewriter *writer = save_file_open (save->sf, m, cmd->location);
6368 gsl_matrix_free (m);
6372 const struct caseproto *proto = casewriter_get_proto (writer);
6373 for (size_t y = 0; y < m->size1; y++)
6375 struct ccase *c = case_create (proto);
6376 for (size_t x = 0; x < m->size2; x++)
6378 double d = gsl_matrix_get (m, y, x);
6379 int width = caseproto_get_width (proto, x);
6380 union value *value = case_data_rw_idx (c, x);
6384 memcpy (value->s, &d, width);
6386 casewriter_write (writer, c);
6388 gsl_matrix_free (m);
6393 static struct read_file *
6394 read_file_create (struct matrix_state *s, struct file_handle *fh)
6396 for (size_t i = 0; i < s->n_read_files; i++)
6398 struct read_file *rf = s->read_files[i];
6406 struct read_file *rf = xmalloc (sizeof *rf);
6407 *rf = (struct read_file) { .file = fh };
6409 s->read_files = xrealloc (s->read_files,
6410 (s->n_read_files + 1) * sizeof *s->read_files);
6411 s->read_files[s->n_read_files++] = rf;
6415 static struct dfm_reader *
6416 read_file_open (struct read_file *rf)
6419 rf->reader = dfm_open_reader (rf->file, NULL, rf->encoding);
6424 read_file_destroy (struct read_file *rf)
6428 fh_unref (rf->file);
6429 dfm_close_reader (rf->reader);
6430 free (rf->encoding);
6435 static struct matrix_command *
6436 matrix_read_parse (struct matrix_state *s)
6438 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6439 *cmd = (struct matrix_command) {
6441 .read = { .format = FMT_F },
6444 struct file_handle *fh = NULL;
6445 char *encoding = NULL;
6446 struct matrix_read *read = &cmd->read;
6447 read->dst = matrix_lvalue_parse (s);
6452 int repetitions = 0;
6453 int record_width = 0;
6454 bool seen_format = false;
6455 while (lex_match (s->lexer, T_SLASH))
6457 if (lex_match_id (s->lexer, "FILE"))
6459 lex_match (s->lexer, T_EQUALS);
6462 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
6466 else if (lex_match_id (s->lexer, "ENCODING"))
6468 lex_match (s->lexer, T_EQUALS);
6469 if (!lex_force_string (s->lexer))
6473 encoding = ss_xstrdup (lex_tokss (s->lexer));
6477 else if (lex_match_id (s->lexer, "FIELD"))
6479 lex_match (s->lexer, T_EQUALS);
6481 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
6483 read->c1 = lex_integer (s->lexer);
6485 if (!lex_force_match (s->lexer, T_TO)
6486 || !lex_force_int_range (s->lexer, "TO", read->c1, INT_MAX))
6488 read->c2 = lex_integer (s->lexer) + 1;
6491 record_width = read->c2 - read->c1;
6492 if (lex_match (s->lexer, T_BY))
6494 if (!lex_force_int_range (s->lexer, "BY", 1,
6495 read->c2 - read->c1))
6497 by = lex_integer (s->lexer);
6500 if (record_width % by)
6502 msg (SE, _("BY %d does not evenly divide record width %d."),
6510 else if (lex_match_id (s->lexer, "SIZE"))
6512 lex_match (s->lexer, T_EQUALS);
6513 matrix_expr_destroy (read->size);
6514 read->size = matrix_parse_exp (s);
6518 else if (lex_match_id (s->lexer, "MODE"))
6520 lex_match (s->lexer, T_EQUALS);
6521 if (lex_match_id (s->lexer, "RECTANGULAR"))
6522 read->symmetric = false;
6523 else if (lex_match_id (s->lexer, "SYMMETRIC"))
6524 read->symmetric = true;
6527 lex_error_expecting (s->lexer, "RECTANGULAR", "SYMMETRIC");
6531 else if (lex_match_id (s->lexer, "REREAD"))
6532 read->reread = true;
6533 else if (lex_match_id (s->lexer, "FORMAT"))
6537 lex_sbc_only_once ("FORMAT");
6542 lex_match (s->lexer, T_EQUALS);
6544 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
6547 const char *p = lex_tokcstr (s->lexer);
6548 if (c_isdigit (p[0]))
6550 repetitions = atoi (p);
6551 p += strspn (p, "0123456789");
6552 if (!fmt_from_name (p, &read->format))
6554 lex_error (s->lexer, _("Unknown format %s."), p);
6559 else if (fmt_from_name (p, &read->format))
6563 struct fmt_spec format;
6564 if (!parse_format_specifier (s->lexer, &format))
6566 read->format = format.type;
6572 lex_error_expecting (s->lexer, "FILE", "FIELD", "MODE",
6573 "REREAD", "FORMAT");
6580 lex_sbc_missing ("FIELD");
6584 if (!read->dst->n_indexes && !read->size)
6586 msg (SE, _("SIZE is required for reading data into a full matrix "
6587 "(as opposed to a submatrix)."));
6593 if (s->prev_read_file)
6594 fh = fh_ref (s->prev_read_file);
6597 lex_sbc_missing ("FILE");
6601 fh_unref (s->prev_read_file);
6602 s->prev_read_file = fh_ref (fh);
6604 read->rf = read_file_create (s, fh);
6608 free (read->rf->encoding);
6609 read->rf->encoding = encoding;
6613 /* Field width may be specified in multiple ways:
6616 2. The format on FORMAT.
6617 3. The repetition factor on FORMAT.
6619 (2) and (3) are mutually exclusive.
6621 If more than one of these is present, they must agree. If none of them is
6622 present, then free-field format is used.
6624 if (repetitions > record_width)
6626 msg (SE, _("%d repetitions cannot fit in record width %d."),
6627 repetitions, record_width);
6630 int w = (repetitions ? record_width / repetitions
6636 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
6637 "which implies field width %d, "
6638 "but BY specifies field width %d."),
6639 repetitions, record_width, w, by);
6641 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
6650 matrix_command_destroy (cmd);
6656 parse_error (const struct dfm_reader *reader, enum fmt_type format,
6657 struct substring data, size_t y, size_t x,
6658 int first_column, int last_column, char *error)
6660 int line_number = dfm_get_line_number (reader);
6661 struct msg_location location = {
6662 .file_name = intern_new (dfm_get_file_name (reader)),
6663 .start = { .line = line_number, .column = first_column },
6664 .end = { .line = line_number, .column = last_column },
6666 msg_at (DW, &location, _("Error reading \"%.*s\" as format %s "
6667 "for matrix row %zu, column %zu: %s"),
6668 (int) data.length, data.string, fmt_name (format),
6669 y + 1, x + 1, error);
6670 msg_location_uninit (&location);
6675 matrix_read_set_field (struct matrix_read *read, struct dfm_reader *reader,
6676 gsl_matrix *m, struct substring p, size_t y, size_t x,
6677 const char *line_start)
6679 const char *input_encoding = dfm_reader_get_encoding (reader);
6682 if (fmt_is_numeric (read->format))
6685 error = data_in (p, input_encoding, read->format,
6686 settings_get_fmt_settings (), &v, 0, NULL);
6687 if (!error && v.f == SYSMIS)
6688 error = xstrdup (_("Matrix data may not contain missing value."));
6693 uint8_t s[sizeof (double)];
6694 union value v = { .s = s };
6695 error = data_in (p, input_encoding, read->format,
6696 settings_get_fmt_settings (), &v, sizeof s, "UTF-8");
6697 memcpy (&f, s, sizeof f);
6702 int c1 = utf8_count_columns (line_start, p.string - line_start) + 1;
6703 int nc = ss_utf8_count_columns (p);
6704 int c2 = c1 + MAX (1, nc) - 1;
6705 parse_error (reader, read->format, p, y, x, c1, c2, error);
6709 gsl_matrix_set (m, y, x, f);
6710 if (read->symmetric && x != y)
6711 gsl_matrix_set (m, x, y, f);
6716 matrix_read_line (struct matrix_command *cmd, struct dfm_reader *reader,
6717 struct substring *line, const char **startp)
6719 struct matrix_read *read = &cmd->read;
6720 if (dfm_eof (reader))
6722 msg_at (SE, cmd->location,
6723 _("Unexpected end of file reading matrix data."));
6726 dfm_expand_tabs (reader);
6727 struct substring record = dfm_get_record (reader);
6728 /* XXX need to recode record into UTF-8 */
6729 *startp = record.string;
6730 *line = ss_utf8_columns (record, read->c1 - 1, read->c2 - read->c1);
6735 matrix_read (struct matrix_command *cmd, struct dfm_reader *reader,
6738 struct matrix_read *read = &cmd->read;
6739 for (size_t y = 0; y < m->size1; y++)
6741 size_t nx = read->symmetric ? y + 1 : m->size2;
6743 struct substring line = ss_empty ();
6744 const char *line_start = line.string;
6745 for (size_t x = 0; x < nx; x++)
6752 ss_ltrim (&line, ss_cstr (" ,"));
6753 if (!ss_is_empty (line))
6755 if (!matrix_read_line (cmd, reader, &line, &line_start))
6757 dfm_forward_record (reader);
6760 ss_get_bytes (&line, ss_cspan (line, ss_cstr (" ,")), &p);
6764 if (!matrix_read_line (cmd, reader, &line, &line_start))
6766 size_t fields_per_line = (read->c2 - read->c1) / read->w;
6767 int f = x % fields_per_line;
6768 if (f == fields_per_line - 1)
6769 dfm_forward_record (reader);
6771 p = ss_substr (line, read->w * f, read->w);
6774 matrix_read_set_field (read, reader, m, p, y, x, line_start);
6778 dfm_forward_record (reader);
6781 ss_ltrim (&line, ss_cstr (" ,"));
6782 if (!ss_is_empty (line))
6784 int line_number = dfm_get_line_number (reader);
6785 int c1 = utf8_count_columns (line_start,
6786 line.string - line_start) + 1;
6787 int c2 = c1 + ss_utf8_count_columns (line) - 1;
6788 struct msg_location location = {
6789 .file_name = intern_new (dfm_get_file_name (reader)),
6790 .start = { .line = line_number, .column = c1 },
6791 .end = { .line = line_number, .column = c2 },
6793 msg_at (DW, &location,
6794 _("Trailing garbage following data for matrix row %zu."),
6796 msg_location_uninit (&location);
6803 matrix_read_execute (struct matrix_command *cmd)
6805 struct matrix_read *read = &cmd->read;
6806 struct index_vector iv0, iv1;
6807 if (!matrix_lvalue_evaluate (read->dst, &iv0, &iv1))
6810 size_t size[2] = { SIZE_MAX, SIZE_MAX };
6813 gsl_matrix *m = matrix_expr_evaluate (read->size);
6819 msg_at (SE, matrix_expr_location (read->size),
6820 _("SIZE must evaluate to a scalar or a 2-element vector, "
6821 "not a %zu×%zu matrix."), m->size1, m->size2);
6822 gsl_matrix_free (m);
6823 index_vector_uninit (&iv0);
6824 index_vector_uninit (&iv1);
6828 gsl_vector v = to_vector (m);
6832 d[0] = gsl_vector_get (&v, 0);
6835 else if (v.size == 2)
6837 d[0] = gsl_vector_get (&v, 0);
6838 d[1] = gsl_vector_get (&v, 1);
6842 msg_at (SE, matrix_expr_location (read->size),
6843 _("SIZE must evaluate to a scalar or a 2-element vector, "
6844 "not a %zu×%zu matrix."),
6845 m->size1, m->size2),
6846 gsl_matrix_free (m);
6847 index_vector_uninit (&iv0);
6848 index_vector_uninit (&iv1);
6851 gsl_matrix_free (m);
6853 if (d[0] < 0 || d[0] > SIZE_MAX || d[1] < 0 || d[1] > SIZE_MAX)
6855 msg_at (SE, matrix_expr_location (read->size),
6856 _("Matrix dimensions %g×%g specified on SIZE "
6857 "are outside valid range."),
6859 index_vector_uninit (&iv0);
6860 index_vector_uninit (&iv1);
6868 if (read->dst->n_indexes)
6870 size_t submatrix_size[2];
6871 if (read->dst->n_indexes == 2)
6873 submatrix_size[0] = iv0.n;
6874 submatrix_size[1] = iv1.n;
6876 else if (read->dst->var->value->size1 == 1)
6878 submatrix_size[0] = 1;
6879 submatrix_size[1] = iv0.n;
6883 submatrix_size[0] = iv0.n;
6884 submatrix_size[1] = 1;
6889 if (size[0] != submatrix_size[0] || size[1] != submatrix_size[1])
6891 msg_at (SE, cmd->location,
6892 _("Dimensions specified on SIZE differ from dimensions "
6893 "of destination submatrix."));
6894 msg_at (SN, matrix_expr_location (read->size),
6895 _("SIZE specifies dimensions %zu×%zu."),
6897 msg_at (SN, read->dst->full_location,
6898 _("Destination submatrix has dimensions %zu×%zu."),
6899 submatrix_size[0], submatrix_size[1]);
6900 index_vector_uninit (&iv0);
6901 index_vector_uninit (&iv1);
6907 size[0] = submatrix_size[0];
6908 size[1] = submatrix_size[1];
6912 struct dfm_reader *reader = read_file_open (read->rf);
6914 dfm_reread_record (reader, 1);
6916 if (read->symmetric && size[0] != size[1])
6918 msg_at (SE, cmd->location,
6919 _("Cannot read non-square %zu×%zu matrix "
6920 "using READ with MODE=SYMMETRIC."),
6922 index_vector_uninit (&iv0);
6923 index_vector_uninit (&iv1);
6926 gsl_matrix *tmp = gsl_matrix_calloc (size[0], size[1]);
6927 matrix_read (cmd, reader, tmp);
6928 matrix_lvalue_assign (read->dst, &iv0, &iv1, tmp, cmd->location);
6933 static struct write_file *
6934 write_file_create (struct matrix_state *s, struct file_handle *fh)
6936 for (size_t i = 0; i < s->n_write_files; i++)
6938 struct write_file *wf = s->write_files[i];
6946 struct write_file *wf = xmalloc (sizeof *wf);
6947 *wf = (struct write_file) { .file = fh };
6949 s->write_files = xrealloc (s->write_files,
6950 (s->n_write_files + 1) * sizeof *s->write_files);
6951 s->write_files[s->n_write_files++] = wf;
6955 static struct dfm_writer *
6956 write_file_open (struct write_file *wf)
6959 wf->writer = dfm_open_writer (wf->file, wf->encoding);
6964 write_file_destroy (struct write_file *wf)
6970 dfm_put_record_utf8 (wf->writer, wf->held->s.ss.string,
6971 wf->held->s.ss.length);
6972 u8_line_destroy (wf->held);
6976 fh_unref (wf->file);
6977 dfm_close_writer (wf->writer);
6978 free (wf->encoding);
6983 static struct matrix_command *
6984 matrix_write_parse (struct matrix_state *s)
6986 struct matrix_command *cmd = xmalloc (sizeof *cmd);
6987 *cmd = (struct matrix_command) {
6991 struct file_handle *fh = NULL;
6992 char *encoding = NULL;
6993 struct matrix_write *write = &cmd->write;
6994 write->expression = matrix_parse_exp (s);
6995 if (!write->expression)
6999 int repetitions = 0;
7000 int record_width = 0;
7001 enum fmt_type format = FMT_F;
7002 bool has_format = false;
7003 while (lex_match (s->lexer, T_SLASH))
7005 if (lex_match_id (s->lexer, "OUTFILE"))
7007 lex_match (s->lexer, T_EQUALS);
7010 fh = fh_parse (s->lexer, FH_REF_FILE, NULL);
7014 else if (lex_match_id (s->lexer, "ENCODING"))
7016 lex_match (s->lexer, T_EQUALS);
7017 if (!lex_force_string (s->lexer))
7021 encoding = ss_xstrdup (lex_tokss (s->lexer));
7025 else if (lex_match_id (s->lexer, "FIELD"))
7027 lex_match (s->lexer, T_EQUALS);
7029 if (!lex_force_int_range (s->lexer, "FIELD", 1, INT_MAX))
7031 write->c1 = lex_integer (s->lexer);
7033 if (!lex_force_match (s->lexer, T_TO)
7034 || !lex_force_int_range (s->lexer, "TO", write->c1, INT_MAX))
7036 write->c2 = lex_integer (s->lexer) + 1;
7039 record_width = write->c2 - write->c1;
7040 if (lex_match (s->lexer, T_BY))
7042 if (!lex_force_int_range (s->lexer, "BY", 1,
7043 write->c2 - write->c1))
7045 by = lex_integer (s->lexer);
7048 if (record_width % by)
7050 msg (SE, _("BY %d does not evenly divide record width %d."),
7058 else if (lex_match_id (s->lexer, "MODE"))
7060 lex_match (s->lexer, T_EQUALS);
7061 if (lex_match_id (s->lexer, "RECTANGULAR"))
7062 write->triangular = false;
7063 else if (lex_match_id (s->lexer, "TRIANGULAR"))
7064 write->triangular = true;
7067 lex_error_expecting (s->lexer, "RECTANGULAR", "TRIANGULAR");
7071 else if (lex_match_id (s->lexer, "HOLD"))
7073 else if (lex_match_id (s->lexer, "FORMAT"))
7075 if (has_format || write->format)
7077 lex_sbc_only_once ("FORMAT");
7081 lex_match (s->lexer, T_EQUALS);
7083 if (lex_token (s->lexer) != T_STRING && !lex_force_id (s->lexer))
7086 const char *p = lex_tokcstr (s->lexer);
7087 if (c_isdigit (p[0]))
7089 repetitions = atoi (p);
7090 p += strspn (p, "0123456789");
7091 if (!fmt_from_name (p, &format))
7093 lex_error (s->lexer, _("Unknown format %s."), p);
7099 else if (fmt_from_name (p, &format))
7106 struct fmt_spec spec;
7107 if (!parse_format_specifier (s->lexer, &spec))
7109 write->format = xmemdup (&spec, sizeof spec);
7114 lex_error_expecting (s->lexer, "OUTFILE", "FIELD", "MODE",
7122 lex_sbc_missing ("FIELD");
7128 if (s->prev_write_file)
7129 fh = fh_ref (s->prev_write_file);
7132 lex_sbc_missing ("OUTFILE");
7136 fh_unref (s->prev_write_file);
7137 s->prev_write_file = fh_ref (fh);
7139 write->wf = write_file_create (s, fh);
7143 free (write->wf->encoding);
7144 write->wf->encoding = encoding;
7148 /* Field width may be specified in multiple ways:
7151 2. The format on FORMAT.
7152 3. The repetition factor on FORMAT.
7154 (2) and (3) are mutually exclusive.
7156 If more than one of these is present, they must agree. If none of them is
7157 present, then free-field format is used.
7159 if (repetitions > record_width)
7161 msg (SE, _("%d repetitions cannot fit in record width %d."),
7162 repetitions, record_width);
7165 int w = (repetitions ? record_width / repetitions
7166 : write->format ? write->format->w
7171 msg (SE, _("FORMAT specifies %d repetitions with record width %d, "
7172 "which implies field width %d, "
7173 "but BY specifies field width %d."),
7174 repetitions, record_width, w, by);
7176 msg (SE, _("FORMAT specifies field width %d but BY specifies %d."),
7180 if (w && !write->format)
7182 write->format = xmalloc (sizeof *write->format);
7183 *write->format = (struct fmt_spec) { .type = format, .w = w };
7185 if (!fmt_check_output (write->format))
7189 if (write->format && fmt_var_width (write->format) > sizeof (double))
7191 char s[FMT_STRING_LEN_MAX + 1];
7192 fmt_to_string (write->format, s);
7193 msg (SE, _("Format %s is too wide for %zu-byte matrix eleemnts."),
7194 s, sizeof (double));
7202 matrix_command_destroy (cmd);
7207 matrix_write_execute (struct matrix_write *write)
7209 gsl_matrix *m = matrix_expr_evaluate (write->expression);
7213 if (write->triangular && m->size1 != m->size2)
7215 msg_at (SE, matrix_expr_location (write->expression),
7216 _("WRITE with MODE=TRIANGULAR requires a square matrix but "
7217 "the matrix to be written has dimensions %zu×%zu."),
7218 m->size1, m->size2);
7219 gsl_matrix_free (m);
7223 struct dfm_writer *writer = write_file_open (write->wf);
7224 if (!writer || !m->size1)
7226 gsl_matrix_free (m);
7230 const struct fmt_settings *settings = settings_get_fmt_settings ();
7231 struct u8_line *line = write->wf->held;
7232 for (size_t y = 0; y < m->size1; y++)
7236 line = xmalloc (sizeof *line);
7237 u8_line_init (line);
7239 size_t nx = write->triangular ? y + 1 : m->size2;
7241 for (size_t x = 0; x < nx; x++)
7244 double f = gsl_matrix_get (m, y, x);
7248 if (fmt_is_numeric (write->format->type))
7251 v.s = (uint8_t *) &f;
7252 s = data_out (&v, NULL, write->format, settings);
7256 s = xmalloc (DBL_BUFSIZE_BOUND);
7257 if (c_dtoastr (s, DBL_BUFSIZE_BOUND, FTOASTR_UPPER_E, 0, f)
7258 >= DBL_BUFSIZE_BOUND)
7261 size_t len = strlen (s);
7262 int width = u8_width (CHAR_CAST (const uint8_t *, s), len, UTF8);
7263 if (width + x0 > write->c2)
7265 dfm_put_record_utf8 (writer, line->s.ss.string,
7267 u8_line_clear (line);
7270 u8_line_put (line, x0, x0 + width, s, len);
7273 x0 += write->format ? write->format->w : width + 1;
7276 if (y + 1 >= m->size1 && write->hold)
7278 dfm_put_record_utf8 (writer, line->s.ss.string, line->s.ss.length);
7279 u8_line_clear (line);
7283 u8_line_destroy (line);
7287 write->wf->held = line;
7289 gsl_matrix_free (m);
7294 static struct matrix_command *
7295 matrix_get_parse (struct matrix_state *s)
7297 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7298 *cmd = (struct matrix_command) {
7301 .dataset = s->dataset,
7302 .user = { .treatment = MGET_ERROR },
7303 .system = { .treatment = MGET_ERROR },
7307 struct matrix_get *get = &cmd->get;
7308 get->dst = matrix_lvalue_parse (s);
7312 while (lex_match (s->lexer, T_SLASH))
7314 if (lex_match_id (s->lexer, "FILE"))
7316 lex_match (s->lexer, T_EQUALS);
7318 fh_unref (get->file);
7319 if (lex_match (s->lexer, T_ASTERISK))
7323 get->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
7328 else if (lex_match_id (s->lexer, "ENCODING"))
7330 lex_match (s->lexer, T_EQUALS);
7331 if (!lex_force_string (s->lexer))
7334 free (get->encoding);
7335 get->encoding = ss_xstrdup (lex_tokss (s->lexer));
7339 else if (lex_match_id (s->lexer, "VARIABLES"))
7341 lex_match (s->lexer, T_EQUALS);
7345 lex_sbc_only_once ("VARIABLES");
7349 if (!var_syntax_parse (s->lexer, &get->vars, &get->n_vars))
7352 else if (lex_match_id (s->lexer, "NAMES"))
7354 lex_match (s->lexer, T_EQUALS);
7355 if (!lex_force_id (s->lexer))
7358 struct substring name = lex_tokss (s->lexer);
7359 get->names = matrix_var_lookup (s, name);
7361 get->names = matrix_var_create (s, name);
7364 else if (lex_match_id (s->lexer, "MISSING"))
7366 lex_match (s->lexer, T_EQUALS);
7367 if (lex_match_id (s->lexer, "ACCEPT"))
7368 get->user.treatment = MGET_ACCEPT;
7369 else if (lex_match_id (s->lexer, "OMIT"))
7370 get->user.treatment = MGET_OMIT;
7371 else if (lex_is_number (s->lexer))
7373 get->user.treatment = MGET_RECODE;
7374 get->user.substitute = lex_number (s->lexer);
7379 lex_error (s->lexer, NULL);
7383 else if (lex_match_id (s->lexer, "SYSMIS"))
7385 lex_match (s->lexer, T_EQUALS);
7386 if (lex_match_id (s->lexer, "OMIT"))
7387 get->system.treatment = MGET_OMIT;
7388 else if (lex_is_number (s->lexer))
7390 get->system.treatment = MGET_RECODE;
7391 get->system.substitute = lex_number (s->lexer);
7396 lex_error (s->lexer, NULL);
7402 lex_error_expecting (s->lexer, "FILE", "VARIABLES", "NAMES",
7403 "MISSING", "SYSMIS");
7408 if (get->user.treatment != MGET_ACCEPT)
7409 get->system.treatment = MGET_ERROR;
7414 matrix_command_destroy (cmd);
7419 matrix_get_execute__ (struct matrix_command *cmd, struct casereader *reader,
7420 const struct dictionary *dict)
7422 struct matrix_get *get = &cmd->get;
7423 struct variable **vars;
7428 if (!var_syntax_evaluate (get->vars, get->n_vars, dict,
7429 &vars, &n_vars, PV_NUMERIC))
7434 n_vars = dict_get_n_vars (dict);
7435 vars = xnmalloc (n_vars, sizeof *vars);
7436 for (size_t i = 0; i < n_vars; i++)
7438 struct variable *var = dict_get_var (dict, i);
7439 if (!var_is_numeric (var))
7441 msg_at (SE, cmd->location, _("Variable %s is not numeric."),
7442 var_get_name (var));
7452 gsl_matrix *names = gsl_matrix_alloc (n_vars, 1);
7453 for (size_t i = 0; i < n_vars; i++)
7455 char s[sizeof (double)];
7457 buf_copy_str_rpad (s, sizeof s, var_get_name (vars[i]), ' ');
7458 memcpy (&f, s, sizeof f);
7459 gsl_matrix_set (names, i, 0, f);
7462 gsl_matrix_free (get->names->value);
7463 get->names->value = names;
7467 gsl_matrix *m = gsl_matrix_alloc (4, n_vars);
7468 long long int casenum = 1;
7470 for (struct ccase *c = casereader_read (reader); c;
7471 c = casereader_read (reader), casenum++)
7473 if (n_rows >= m->size1)
7475 gsl_matrix *p = gsl_matrix_alloc (m->size1 * 2, n_vars);
7476 for (size_t y = 0; y < n_rows; y++)
7477 for (size_t x = 0; x < n_vars; x++)
7478 gsl_matrix_set (p, y, x, gsl_matrix_get (m, y, x));
7479 gsl_matrix_free (m);
7484 for (size_t x = 0; x < n_vars; x++)
7486 const struct variable *var = vars[x];
7487 double d = case_num (c, var);
7490 if (get->system.treatment == MGET_RECODE)
7491 d = get->system.substitute;
7492 else if (get->system.treatment == MGET_OMIT)
7496 msg_at (SE, cmd->location, _("Variable %s in case %lld "
7497 "is system-missing."),
7498 var_get_name (var), casenum);
7502 else if (var_is_num_missing (var, d, MV_USER))
7504 if (get->user.treatment == MGET_RECODE)
7505 d = get->user.substitute;
7506 else if (get->user.treatment == MGET_OMIT)
7508 else if (get->user.treatment != MGET_ACCEPT)
7510 msg_at (SE, cmd->location,
7511 _("Variable %s in case %lld has user-missing "
7513 var_get_name (var), casenum, d);
7517 gsl_matrix_set (m, n_rows, x, d);
7528 matrix_lvalue_evaluate_and_assign (get->dst, m, cmd->location);
7531 gsl_matrix_free (m);
7536 matrix_open_casereader (const struct matrix_command *cmd,
7537 const char *command_name, struct file_handle *file,
7538 const char *encoding, struct dataset *dataset,
7539 struct casereader **readerp, struct dictionary **dictp)
7543 *readerp = any_reader_open_and_decode (file, encoding, dictp, NULL);
7544 return *readerp != NULL;
7548 if (dict_get_n_vars (dataset_dict (dataset)) == 0)
7550 msg_at (ME, cmd->location,
7551 _("The %s command cannot read an empty active file."),
7555 *readerp = proc_open (dataset);
7556 *dictp = dict_ref (dataset_dict (dataset));
7562 matrix_close_casereader (struct file_handle *file, struct dataset *dataset,
7563 struct casereader *reader, struct dictionary *dict)
7566 casereader_destroy (reader);
7568 proc_commit (dataset);
7572 matrix_get_execute (struct matrix_command *cmd)
7574 struct matrix_get *get = &cmd->get;
7575 struct casereader *r;
7576 struct dictionary *d;
7577 if (matrix_open_casereader (cmd, "GET", get->file, get->encoding,
7578 get->dataset, &r, &d))
7580 matrix_get_execute__ (cmd, r, d);
7581 matrix_close_casereader (get->file, get->dataset, r, d);
7588 variables_changed (const char *keyword,
7589 const struct string_array *new,
7590 const struct string_array *old)
7596 msg (SE, _("%s may only be specified on MSAVE if it was specified "
7597 "on the first MSAVE within MATRIX."), keyword);
7600 else if (!string_array_equal_case (old, new))
7602 msg (SE, _("%s must specify the same variables each time within "
7603 "a given MATRIX."), keyword);
7611 msave_common_changed (const struct msave_common *old,
7612 const struct msave_common *new)
7614 if (new->outfile && !fh_equal (old->outfile, new->outfile))
7615 msg (SE, _("OUTFILE must name the same file on each MSAVE "
7616 "within a single MATRIX command."));
7617 else if (variables_changed ("VARIABLES", &new->variables, &old->variables)
7618 || variables_changed ("FNAMES", &new->fnames, &old->fnames)
7619 || variables_changed ("SNAMES", &new->snames, &old->snames))
7620 msg_at (SN, old->location,
7621 _("This is the location of the first MSAVE command."));
7629 msave_common_destroy (struct msave_common *common)
7633 msg_location_destroy (common->location);
7634 fh_unref (common->outfile);
7635 string_array_destroy (&common->variables);
7636 string_array_destroy (&common->fnames);
7637 string_array_destroy (&common->snames);
7639 for (size_t i = 0; i < common->n_factors; i++)
7640 matrix_expr_destroy (common->factors[i]);
7641 free (common->factors);
7643 for (size_t i = 0; i < common->n_splits; i++)
7644 matrix_expr_destroy (common->splits[i]);
7645 free (common->splits);
7647 dict_unref (common->dict);
7648 casewriter_destroy (common->writer);
7655 match_rowtype (struct lexer *lexer)
7657 static const char *rowtypes[] = {
7658 "COV", "CORR", "MEAN", "STDDEV", "N", "COUNT"
7660 size_t n_rowtypes = sizeof rowtypes / sizeof *rowtypes;
7662 for (size_t i = 0; i < n_rowtypes; i++)
7663 if (lex_match_id (lexer, rowtypes[i]))
7666 lex_error_expecting_array (lexer, rowtypes, n_rowtypes);
7671 parse_var_names (struct lexer *lexer, struct string_array *sa)
7673 lex_match (lexer, T_EQUALS);
7675 string_array_clear (sa);
7677 struct dictionary *dict = dict_create (get_default_encoding ());
7680 bool ok = parse_DATA_LIST_vars (lexer, dict, &names, &n_names,
7681 PV_NO_DUPLICATE | PV_NO_SCRATCH);
7686 for (size_t i = 0; i < n_names; i++)
7687 if (ss_equals_case (ss_cstr (names[i]), ss_cstr ("ROWTYPE_"))
7688 || ss_equals_case (ss_cstr (names[i]), ss_cstr ("VARNAME_")))
7690 msg (SE, _("Variable name %s is reserved."), names[i]);
7691 for (size_t j = 0; j < n_names; j++)
7697 string_array_clear (sa);
7698 sa->strings = names;
7699 sa->n = sa->allocated = n_names;
7704 static struct matrix_command *
7705 matrix_msave_parse (struct matrix_state *s)
7707 int start_ofs = lex_ofs (s->lexer);
7709 struct msave_common *common = xmalloc (sizeof *common);
7710 *common = (struct msave_common) { .outfile = NULL };
7712 struct matrix_command *cmd = xmalloc (sizeof *cmd);
7713 *cmd = (struct matrix_command) { .type = MCMD_MSAVE, .msave = { .expr = NULL } };
7715 struct matrix_expr *splits = NULL;
7716 struct matrix_expr *factors = NULL;
7718 struct matrix_msave *msave = &cmd->msave;
7719 msave->expr = matrix_parse_exp (s);
7723 while (lex_match (s->lexer, T_SLASH))
7725 if (lex_match_id (s->lexer, "TYPE"))
7727 lex_match (s->lexer, T_EQUALS);
7729 msave->rowtype = match_rowtype (s->lexer);
7730 if (!msave->rowtype)
7733 else if (lex_match_id (s->lexer, "OUTFILE"))
7735 lex_match (s->lexer, T_EQUALS);
7737 fh_unref (common->outfile);
7738 common->outfile = fh_parse (s->lexer, FH_REF_FILE, NULL);
7739 if (!common->outfile)
7742 else if (lex_match_id (s->lexer, "VARIABLES"))
7744 if (!parse_var_names (s->lexer, &common->variables))
7747 else if (lex_match_id (s->lexer, "FNAMES"))
7749 if (!parse_var_names (s->lexer, &common->fnames))
7752 else if (lex_match_id (s->lexer, "SNAMES"))
7754 if (!parse_var_names (s->lexer, &common->snames))
7757 else if (lex_match_id (s->lexer, "SPLIT"))
7759 lex_match (s->lexer, T_EQUALS);
7761 matrix_expr_destroy (splits);
7762 splits = matrix_parse_exp (s);
7766 else if (lex_match_id (s->lexer, "FACTOR"))
7768 lex_match (s->lexer, T_EQUALS);
7770 matrix_expr_destroy (factors);
7771 factors = matrix_parse_exp (s);
7777 lex_error_expecting (s->lexer, "TYPE", "OUTFILE", "VARIABLES",
7778 "FNAMES", "SNAMES", "SPLIT", "FACTOR");
7782 if (!msave->rowtype)
7784 lex_sbc_missing ("TYPE");
7788 if (!s->msave_common)
7790 if (common->fnames.n && !factors)
7792 msg (SE, _("FNAMES requires FACTOR."));
7795 if (common->snames.n && !splits)
7797 msg (SE, _("SNAMES requires SPLIT."));
7800 if (!common->outfile)
7802 lex_sbc_missing ("OUTFILE");
7805 common->location = lex_ofs_location (s->lexer, start_ofs,
7806 lex_ofs (s->lexer));
7807 msg_location_remove_columns (common->location);
7808 s->msave_common = common;
7812 if (msave_common_changed (s->msave_common, common))
7814 msave_common_destroy (common);
7816 msave->common = s->msave_common;
7818 struct msave_common *c = s->msave_common;
7821 if (c->n_factors >= c->allocated_factors)
7822 c->factors = x2nrealloc (c->factors, &c->allocated_factors,
7823 sizeof *c->factors);
7824 c->factors[c->n_factors++] = factors;
7826 if (c->n_factors > 0)
7827 msave->factors = c->factors[c->n_factors - 1];
7831 if (c->n_splits >= c->allocated_splits)
7832 c->splits = x2nrealloc (c->splits, &c->allocated_splits,
7834 c->splits[c->n_splits++] = splits;
7836 if (c->n_splits > 0)
7837 msave->splits = c->splits[c->n_splits - 1];
7842 matrix_expr_destroy (splits);
7843 matrix_expr_destroy (factors);
7844 msave_common_destroy (common);
7845 matrix_command_destroy (cmd);
7850 matrix_expr_evaluate_vector (const struct matrix_expr *e, const char *name)
7852 gsl_matrix *m = matrix_expr_evaluate (e);
7858 msg_at (SE, matrix_expr_location (e),
7859 _("%s expression must evaluate to vector, "
7860 "not a %zu×%zu matrix."),
7861 name, m->size1, m->size2);
7862 gsl_matrix_free (m);
7866 return matrix_to_vector (m);
7870 msave_add_vars (struct dictionary *d, const struct string_array *vars)
7872 for (size_t i = 0; i < vars->n; i++)
7873 if (!dict_create_var (d, vars->strings[i], 0))
7874 return vars->strings[i];
7878 static struct dictionary *
7879 msave_create_dict (const struct msave_common *common,
7880 const struct msg_location *location)
7882 struct dictionary *dict = dict_create (get_default_encoding ());
7884 const char *dup_split = msave_add_vars (dict, &common->snames);
7887 /* Should not be possible because the parser ensures that the names are
7892 dict_create_var_assert (dict, "ROWTYPE_", 8);
7894 const char *dup_factor = msave_add_vars (dict, &common->fnames);
7897 msg_at (SE, location, _("Duplicate or invalid FACTOR variable name %s."),
7902 dict_create_var_assert (dict, "VARNAME_", 8);
7904 const char *dup_var = msave_add_vars (dict, &common->variables);
7907 msg_at (SE, location, _("Duplicate or invalid variable name %s."),
7920 matrix_msave_execute (struct matrix_command *cmd)
7922 struct matrix_msave *msave = &cmd->msave;
7923 struct msave_common *common = msave->common;
7924 gsl_matrix *m = NULL;
7925 gsl_vector *factors = NULL;
7926 gsl_vector *splits = NULL;
7928 m = matrix_expr_evaluate (msave->expr);
7932 if (!common->variables.n)
7933 for (size_t i = 0; i < m->size2; i++)
7934 string_array_append_nocopy (&common->variables,
7935 xasprintf ("COL%zu", i + 1));
7936 else if (m->size2 != common->variables.n)
7938 msg_at (SE, matrix_expr_location (msave->expr),
7939 _("Matrix on MSAVE has %zu columns but there are %zu variables."),
7940 m->size2, common->variables.n);
7946 factors = matrix_expr_evaluate_vector (msave->factors, "FACTOR");
7950 if (!common->fnames.n)
7951 for (size_t i = 0; i < factors->size; i++)
7952 string_array_append_nocopy (&common->fnames,
7953 xasprintf ("FAC%zu", i + 1));
7954 else if (factors->size != common->fnames.n)
7956 msg_at (SE, matrix_expr_location (msave->factors),
7957 _("There are %zu factor variables, "
7958 "but %zu factor values were supplied."),
7959 common->fnames.n, factors->size);
7966 splits = matrix_expr_evaluate_vector (msave->splits, "SPLIT");
7970 if (!common->snames.n)
7971 for (size_t i = 0; i < splits->size; i++)
7972 string_array_append_nocopy (&common->snames,
7973 xasprintf ("SPL%zu", i + 1));
7974 else if (splits->size != common->snames.n)
7976 msg_at (SE, matrix_expr_location (msave->splits),
7977 _("There are %zu split variables, "
7978 "but %zu split values were supplied."),
7979 common->snames.n, splits->size);
7984 if (!common->writer)
7986 struct dictionary *dict = msave_create_dict (common, cmd->location);
7990 common->writer = any_writer_open (common->outfile, dict);
7991 if (!common->writer)
7997 common->dict = dict;
8000 bool matrix = (!strcmp (msave->rowtype, "COV")
8001 || !strcmp (msave->rowtype, "CORR"));
8002 for (size_t y = 0; y < m->size1; y++)
8004 struct ccase *c = case_create (dict_get_proto (common->dict));
8007 /* Split variables */
8009 for (size_t i = 0; i < splits->size; i++)
8010 case_data_rw_idx (c, idx++)->f = gsl_vector_get (splits, i);
8013 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8014 msave->rowtype, ' ');
8018 for (size_t i = 0; i < factors->size; i++)
8019 *case_num_rw_idx (c, idx++) = gsl_vector_get (factors, i);
8022 const char *varname_ = (matrix && y < common->variables.n
8023 ? common->variables.strings[y]
8025 buf_copy_str_rpad (CHAR_CAST (char *, case_data_rw_idx (c, idx++)->s), 8,
8028 /* Continuous variables. */
8029 for (size_t x = 0; x < m->size2; x++)
8030 case_data_rw_idx (c, idx++)->f = gsl_matrix_get (m, y, x);
8031 casewriter_write (common->writer, c);
8035 gsl_matrix_free (m);
8036 gsl_vector_free (factors);
8037 gsl_vector_free (splits);
8042 static struct matrix_command *
8043 matrix_mget_parse (struct matrix_state *s)
8045 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8046 *cmd = (struct matrix_command) {
8050 .rowtypes = STRINGI_SET_INITIALIZER (cmd->mget.rowtypes),
8054 struct matrix_mget *mget = &cmd->mget;
8056 lex_match (s->lexer, T_SLASH);
8057 while (lex_token (s->lexer) != T_ENDCMD)
8059 if (lex_match_id (s->lexer, "FILE"))
8061 lex_match (s->lexer, T_EQUALS);
8063 fh_unref (mget->file);
8064 mget->file = fh_parse (s->lexer, FH_REF_FILE, s->session);
8068 else if (lex_match_id (s->lexer, "ENCODING"))
8070 lex_match (s->lexer, T_EQUALS);
8071 if (!lex_force_string (s->lexer))
8074 free (mget->encoding);
8075 mget->encoding = ss_xstrdup (lex_tokss (s->lexer));
8079 else if (lex_match_id (s->lexer, "TYPE"))
8081 lex_match (s->lexer, T_EQUALS);
8082 while (lex_token (s->lexer) != T_SLASH
8083 && lex_token (s->lexer) != T_ENDCMD)
8085 const char *rowtype = match_rowtype (s->lexer);
8089 stringi_set_insert (&mget->rowtypes, rowtype);
8094 lex_error_expecting (s->lexer, "FILE", "TYPE");
8097 lex_match (s->lexer, T_SLASH);
8102 matrix_command_destroy (cmd);
8106 static const struct variable *
8107 get_a8_var (const struct msg_location *loc,
8108 const struct dictionary *d, const char *name)
8110 const struct variable *v = dict_lookup_var (d, name);
8113 msg_at (SE, loc, _("Matrix data file lacks %s variable."), name);
8116 if (var_get_width (v) != 8)
8118 msg_at (SE, loc, _("%s variable in matrix data file must be "
8119 "8-byte string, but it has width %d."),
8120 name, var_get_width (v));
8127 var_changed (const struct ccase *ca, const struct ccase *cb,
8128 const struct variable *var)
8131 ? !value_equal (case_data (ca, var), case_data (cb, var),
8132 var_get_width (var))
8137 vars_changed (const struct ccase *ca, const struct ccase *cb,
8138 const struct dictionary *d,
8139 size_t first_var, size_t n_vars)
8141 for (size_t i = 0; i < n_vars; i++)
8143 const struct variable *v = dict_get_var (d, first_var + i);
8144 if (var_changed (ca, cb, v))
8151 vars_all_missing (const struct ccase *c, const struct dictionary *d,
8152 size_t first_var, size_t n_vars)
8154 for (size_t i = 0; i < n_vars; i++)
8155 if (case_num (c, dict_get_var (d, first_var + i)) != SYSMIS)
8161 matrix_mget_commit_var (struct ccase **rows, size_t n_rows,
8162 const struct dictionary *d,
8163 const struct variable *rowtype_var,
8164 const struct stringi_set *accepted_rowtypes,
8165 struct matrix_state *s,
8166 size_t ss, size_t sn, size_t si,
8167 size_t fs, size_t fn, size_t fi,
8168 size_t cs, size_t cn,
8169 struct pivot_table *pt,
8170 struct pivot_dimension *var_dimension)
8175 /* Is this a matrix for pooled data, either where there are no factor
8176 variables or the factor variables are missing? */
8177 bool pooled = !fn || vars_all_missing (rows[0], d, fs, fn);
8179 struct substring rowtype = case_ss (rows[0], rowtype_var);
8180 ss_rtrim (&rowtype, ss_cstr (" "));
8181 if (!stringi_set_is_empty (accepted_rowtypes)
8182 && !stringi_set_contains_len (accepted_rowtypes,
8183 rowtype.string, rowtype.length))
8186 const char *prefix = (ss_equals_case (rowtype, ss_cstr ("COV")) ? "CV"
8187 : ss_equals_case (rowtype, ss_cstr ("CORR")) ? "CR"
8188 : ss_equals_case (rowtype, ss_cstr ("MEAN")) ? "MN"
8189 : ss_equals_case (rowtype, ss_cstr ("STDDEV")) ? "SD"
8190 : ss_equals_case (rowtype, ss_cstr ("N")) ? "NC"
8191 : ss_equals_case (rowtype, ss_cstr ("COUNT")) ? "CN"
8195 msg (SE, _("Matrix data file contains unknown ROWTYPE_ \"%.*s\"."),
8196 (int) rowtype.length, rowtype.string);
8200 struct string name = DS_EMPTY_INITIALIZER;
8201 ds_put_cstr (&name, prefix);
8203 ds_put_format (&name, "F%zu", fi);
8205 ds_put_format (&name, "S%zu", si);
8207 struct matrix_var *mv = matrix_var_lookup (s, ds_ss (&name));
8209 mv = matrix_var_create (s, ds_ss (&name));
8212 msg (SW, _("Matrix data file contains variable with existing name %s."),
8214 goto exit_free_name;
8217 gsl_matrix *m = gsl_matrix_alloc (n_rows, cn);
8218 size_t n_missing = 0;
8219 for (size_t y = 0; y < n_rows; y++)
8221 for (size_t x = 0; x < cn; x++)
8223 struct variable *var = dict_get_var (d, cs + x);
8224 double value = case_num (rows[y], var);
8225 if (var_is_num_missing (var, value, MV_ANY))
8230 gsl_matrix_set (m, y, x, value);
8234 int var_index = pivot_category_create_leaf (
8235 var_dimension->root, pivot_value_new_user_text (ds_cstr (&name), SIZE_MAX));
8236 double values[] = { n_rows, cn };
8237 for (size_t j = 0; j < sn; j++)
8239 struct variable *var = dict_get_var (d, ss + j);
8240 const union value *value = case_data (rows[0], var);
8241 pivot_table_put2 (pt, j, var_index,
8242 pivot_value_new_var_value (var, value));
8244 for (size_t j = 0; j < fn; j++)
8246 struct variable *var = dict_get_var (d, fs + j);
8247 const union value sysmis = { .f = SYSMIS };
8248 const union value *value = pooled ? &sysmis : case_data (rows[0], var);
8249 pivot_table_put2 (pt, j + sn, var_index,
8250 pivot_value_new_var_value (var, value));
8252 for (size_t j = 0; j < sizeof values / sizeof *values; j++)
8253 pivot_table_put2 (pt, j + sn + fn, var_index,
8254 pivot_value_new_integer (values[j]));
8257 msg (SE, ngettext ("Matrix data file variable %s contains a missing "
8258 "value, which was treated as zero.",
8259 "Matrix data file variable %s contains %zu missing "
8260 "values, which were treated as zero.", n_missing),
8261 ds_cstr (&name), n_missing);
8268 for (size_t y = 0; y < n_rows; y++)
8269 case_unref (rows[y]);
8273 matrix_mget_execute__ (struct matrix_command *cmd, struct casereader *r,
8274 const struct dictionary *d)
8276 struct matrix_mget *mget = &cmd->mget;
8277 const struct msg_location *loc = cmd->location;
8278 const struct variable *rowtype_ = get_a8_var (loc, d, "ROWTYPE_");
8279 const struct variable *varname_ = get_a8_var (loc, d, "VARNAME_");
8280 if (!rowtype_ || !varname_)
8283 if (var_get_dict_index (rowtype_) >= var_get_dict_index (varname_))
8286 _("ROWTYPE_ must precede VARNAME_ in matrix data file."));
8289 if (var_get_dict_index (varname_) + 1 >= dict_get_n_vars (d))
8291 msg_at (SE, loc, _("Matrix data file contains no continuous variables."));
8295 for (size_t i = 0; i < dict_get_n_vars (d); i++)
8297 const struct variable *v = dict_get_var (d, i);
8298 if (v != rowtype_ && v != varname_ && var_get_width (v) != 0)
8301 _("Matrix data file contains unexpected string variable %s."),
8307 /* SPLIT variables. */
8309 size_t sn = var_get_dict_index (rowtype_);
8310 struct ccase *sc = NULL;
8313 /* FACTOR variables. */
8314 size_t fs = var_get_dict_index (rowtype_) + 1;
8315 size_t fn = var_get_dict_index (varname_) - var_get_dict_index (rowtype_) - 1;
8316 struct ccase *fc = NULL;
8319 /* Continuous variables. */
8320 size_t cs = var_get_dict_index (varname_) + 1;
8321 size_t cn = dict_get_n_vars (d) - cs;
8322 struct ccase *cc = NULL;
8325 struct pivot_table *pt = pivot_table_create (
8326 N_("Matrix Variables Created by MGET"));
8327 struct pivot_dimension *attr_dimension = pivot_dimension_create (
8328 pt, PIVOT_AXIS_COLUMN, N_("Attribute"));
8329 struct pivot_dimension *var_dimension = pivot_dimension_create (
8330 pt, PIVOT_AXIS_ROW, N_("Variable"));
8333 struct pivot_category *splits = pivot_category_create_group (
8334 attr_dimension->root, N_("Split Values"));
8335 for (size_t i = 0; i < sn; i++)
8336 pivot_category_create_leaf (splits, pivot_value_new_variable (
8337 dict_get_var (d, ss + i)));
8341 struct pivot_category *factors = pivot_category_create_group (
8342 attr_dimension->root, N_("Factors"));
8343 for (size_t i = 0; i < fn; i++)
8344 pivot_category_create_leaf (factors, pivot_value_new_variable (
8345 dict_get_var (d, fs + i)));
8347 pivot_category_create_group (attr_dimension->root, N_("Dimensions"),
8348 N_("Rows"), N_("Columns"));
8351 struct ccase **rows = NULL;
8352 size_t allocated_rows = 0;
8356 while ((c = casereader_read (r)) != NULL)
8358 bool row_has_factors = fn && !vars_all_missing (c, d, fs, fn);
8368 = (sn && (!sc || vars_changed (sc, c, d, ss, sn)) ? SPLITS_CHANGED
8369 : fn && (!fc || vars_changed (fc, c, d, fs, fn)) ? FACTORS_CHANGED
8370 : !cc || var_changed (cc, c, rowtype_) ? ROWTYPE_CHANGED
8373 if (change != NOTHING_CHANGED)
8375 matrix_mget_commit_var (rows, n_rows, d,
8376 rowtype_, &mget->rowtypes,
8387 if (n_rows >= allocated_rows)
8388 rows = x2nrealloc (rows, &allocated_rows, sizeof *rows);
8391 if (change == SPLITS_CHANGED)
8397 /* Reset the factor number, if there are factors. */
8401 if (row_has_factors)
8407 else if (change == FACTORS_CHANGED)
8409 if (row_has_factors)
8415 matrix_mget_commit_var (rows, n_rows, d,
8416 rowtype_, &mget->rowtypes,
8428 if (var_dimension->n_leaves)
8429 pivot_table_submit (pt);
8431 pivot_table_unref (pt);
8435 matrix_mget_execute (struct matrix_command *cmd)
8437 struct matrix_mget *mget = &cmd->mget;
8438 struct casereader *r;
8439 struct dictionary *d;
8440 if (matrix_open_casereader (cmd, "MGET", mget->file, mget->encoding,
8441 mget->state->dataset, &r, &d))
8443 matrix_mget_execute__ (cmd, r, d);
8444 matrix_close_casereader (mget->file, mget->state->dataset, r, d);
8451 matrix_parse_dst_var (struct matrix_state *s, struct matrix_var **varp)
8453 if (!lex_force_id (s->lexer))
8456 *varp = matrix_var_lookup (s, lex_tokss (s->lexer));
8458 *varp = matrix_var_create (s, lex_tokss (s->lexer));
8463 static struct matrix_command *
8464 matrix_eigen_parse (struct matrix_state *s)
8466 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8467 *cmd = (struct matrix_command) {
8469 .eigen = { .expr = NULL }
8472 struct matrix_eigen *eigen = &cmd->eigen;
8473 if (!lex_force_match (s->lexer, T_LPAREN))
8475 eigen->expr = matrix_expr_parse (s);
8477 || !lex_force_match (s->lexer, T_COMMA)
8478 || !matrix_parse_dst_var (s, &eigen->evec)
8479 || !lex_force_match (s->lexer, T_COMMA)
8480 || !matrix_parse_dst_var (s, &eigen->eval)
8481 || !lex_force_match (s->lexer, T_RPAREN))
8487 matrix_command_destroy (cmd);
8492 matrix_eigen_execute (struct matrix_command *cmd)
8494 struct matrix_eigen *eigen = &cmd->eigen;
8495 gsl_matrix *A = matrix_expr_evaluate (eigen->expr);
8498 if (!is_symmetric (A))
8500 msg_at (SE, cmd->location, _("Argument of EIGEN must be symmetric."));
8501 gsl_matrix_free (A);
8505 gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc (A->size1);
8506 gsl_matrix *eval = gsl_matrix_alloc (A->size1, 1);
8507 gsl_vector v_eval = to_vector (eval);
8508 gsl_matrix *evec = gsl_matrix_alloc (A->size1, A->size2);
8509 gsl_eigen_symmv (A, &v_eval, evec, w);
8510 gsl_eigen_symmv_free (w);
8512 gsl_eigen_symmv_sort (&v_eval, evec, GSL_EIGEN_SORT_VAL_DESC);
8514 gsl_matrix_free (A);
8516 gsl_matrix_free (eigen->eval->value);
8517 eigen->eval->value = eval;
8519 gsl_matrix_free (eigen->evec->value);
8520 eigen->evec->value = evec;
8525 static struct matrix_command *
8526 matrix_setdiag_parse (struct matrix_state *s)
8528 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8529 *cmd = (struct matrix_command) {
8530 .type = MCMD_SETDIAG,
8531 .setdiag = { .dst = NULL }
8534 struct matrix_setdiag *setdiag = &cmd->setdiag;
8535 if (!lex_force_match (s->lexer, T_LPAREN) || !lex_force_id (s->lexer))
8538 setdiag->dst = matrix_var_lookup (s, lex_tokss (s->lexer));
8541 lex_error (s->lexer, _("Unknown variable %s."), lex_tokcstr (s->lexer));
8546 if (!lex_force_match (s->lexer, T_COMMA))
8549 setdiag->expr = matrix_expr_parse (s);
8553 if (!lex_force_match (s->lexer, T_RPAREN))
8559 matrix_command_destroy (cmd);
8564 matrix_setdiag_execute (struct matrix_command *cmd)
8566 struct matrix_setdiag *setdiag = &cmd->setdiag;
8567 gsl_matrix *dst = setdiag->dst->value;
8570 msg_at (SE, cmd->location,
8571 _("SETDIAG destination matrix %s is uninitialized."),
8572 setdiag->dst->name);
8576 gsl_matrix *src = matrix_expr_evaluate (setdiag->expr);
8580 size_t n = MIN (dst->size1, dst->size2);
8581 if (is_scalar (src))
8583 double d = to_scalar (src);
8584 for (size_t i = 0; i < n; i++)
8585 gsl_matrix_set (dst, i, i, d);
8587 else if (is_vector (src))
8589 gsl_vector v = to_vector (src);
8590 for (size_t i = 0; i < n && i < v.size; i++)
8591 gsl_matrix_set (dst, i, i, gsl_vector_get (&v, i));
8594 msg_at (SE, matrix_expr_location (setdiag->expr),
8595 _("SETDIAG argument 2 must be a scalar or a vector, "
8596 "not a %zu×%zu matrix."),
8597 src->size1, src->size2);
8598 gsl_matrix_free (src);
8603 static struct matrix_command *
8604 matrix_svd_parse (struct matrix_state *s)
8606 struct matrix_command *cmd = xmalloc (sizeof *cmd);
8607 *cmd = (struct matrix_command) {
8609 .svd = { .expr = NULL }
8612 struct matrix_svd *svd = &cmd->svd;
8613 if (!lex_force_match (s->lexer, T_LPAREN))
8615 svd->expr = matrix_expr_parse (s);
8617 || !lex_force_match (s->lexer, T_COMMA)
8618 || !matrix_parse_dst_var (s, &svd->u)
8619 || !lex_force_match (s->lexer, T_COMMA)
8620 || !matrix_parse_dst_var (s, &svd->s)
8621 || !lex_force_match (s->lexer, T_COMMA)
8622 || !matrix_parse_dst_var (s, &svd->v)
8623 || !lex_force_match (s->lexer, T_RPAREN))
8629 matrix_command_destroy (cmd);
8634 matrix_svd_execute (struct matrix_svd *svd)
8636 gsl_matrix *m = matrix_expr_evaluate (svd->expr);
8640 if (m->size1 >= m->size2)
8643 gsl_matrix *V = gsl_matrix_alloc (A->size2, A->size2);
8644 gsl_matrix *S = gsl_matrix_calloc (A->size2, A->size2);
8645 gsl_vector Sv = gsl_matrix_diagonal (S).vector;
8646 gsl_vector *work = gsl_vector_alloc (A->size2);
8647 gsl_linalg_SV_decomp (A, V, &Sv, work);
8648 gsl_vector_free (work);
8650 matrix_var_set (svd->u, A);
8651 matrix_var_set (svd->s, S);
8652 matrix_var_set (svd->v, V);
8656 gsl_matrix *At = gsl_matrix_alloc (m->size2, m->size1);
8657 gsl_matrix_transpose_memcpy (At, m);
8658 gsl_matrix_free (m);
8660 gsl_matrix *Vt = gsl_matrix_alloc (At->size2, At->size2);
8661 gsl_matrix *St = gsl_matrix_calloc (At->size2, At->size2);
8662 gsl_vector Stv = gsl_matrix_diagonal (St).vector;
8663 gsl_vector *work = gsl_vector_alloc (At->size2);
8664 gsl_linalg_SV_decomp (At, Vt, &Stv, work);
8665 gsl_vector_free (work);
8667 matrix_var_set (svd->v, At);
8668 matrix_var_set (svd->s, St);
8669 matrix_var_set (svd->u, Vt);
8673 /* The main MATRIX command logic. */
8676 matrix_command_execute (struct matrix_command *cmd)
8681 matrix_compute_execute (cmd);
8685 matrix_print_execute (&cmd->print);
8689 return matrix_do_if_execute (&cmd->do_if);
8692 matrix_loop_execute (&cmd->loop);
8699 matrix_display_execute (&cmd->display);
8703 matrix_release_execute (&cmd->release);
8707 matrix_save_execute (cmd);
8711 matrix_read_execute (cmd);
8715 matrix_write_execute (&cmd->write);
8719 matrix_get_execute (cmd);
8723 matrix_msave_execute (cmd);
8727 matrix_mget_execute (cmd);
8731 matrix_eigen_execute (cmd);
8735 matrix_setdiag_execute (cmd);
8739 matrix_svd_execute (&cmd->svd);
8747 matrix_command_destroy (struct matrix_command *cmd)
8752 msg_location_destroy (cmd->location);
8757 matrix_lvalue_destroy (cmd->compute.lvalue);
8758 matrix_expr_destroy (cmd->compute.rvalue);
8762 matrix_expr_destroy (cmd->print.expression);
8763 free (cmd->print.title);
8764 print_labels_destroy (cmd->print.rlabels);
8765 print_labels_destroy (cmd->print.clabels);
8769 for (size_t i = 0; i < cmd->do_if.n_clauses; i++)
8771 matrix_expr_destroy (cmd->do_if.clauses[i].condition);
8772 matrix_commands_uninit (&cmd->do_if.clauses[i].commands);
8774 free (cmd->do_if.clauses);
8778 matrix_expr_destroy (cmd->loop.start);
8779 matrix_expr_destroy (cmd->loop.end);
8780 matrix_expr_destroy (cmd->loop.increment);
8781 matrix_expr_destroy (cmd->loop.top_condition);
8782 matrix_expr_destroy (cmd->loop.bottom_condition);
8783 matrix_commands_uninit (&cmd->loop.commands);
8793 free (cmd->release.vars);
8797 matrix_expr_destroy (cmd->save.expression);
8801 matrix_lvalue_destroy (cmd->read.dst);
8802 matrix_expr_destroy (cmd->read.size);
8806 matrix_expr_destroy (cmd->write.expression);
8807 free (cmd->write.format);
8811 matrix_lvalue_destroy (cmd->get.dst);
8812 fh_unref (cmd->get.file);
8813 free (cmd->get.encoding);
8814 var_syntax_destroy (cmd->get.vars, cmd->get.n_vars);
8818 matrix_expr_destroy (cmd->msave.expr);
8822 fh_unref (cmd->mget.file);
8823 stringi_set_destroy (&cmd->mget.rowtypes);
8827 matrix_expr_destroy (cmd->eigen.expr);
8831 matrix_expr_destroy (cmd->setdiag.expr);
8835 matrix_expr_destroy (cmd->svd.expr);
8842 matrix_commands_parse (struct matrix_state *s, struct matrix_commands *c,
8843 const char *command_name,
8844 const char *stop1, const char *stop2)
8846 lex_end_of_command (s->lexer);
8847 lex_discard_rest_of_command (s->lexer);
8849 size_t allocated = 0;
8852 while (lex_token (s->lexer) == T_ENDCMD)
8855 if (lex_at_phrase (s->lexer, stop1)
8856 || (stop2 && lex_at_phrase (s->lexer, stop2)))
8859 if (lex_at_phrase (s->lexer, "END MATRIX"))
8861 msg (SE, _("Premature END MATRIX within %s."), command_name);
8865 struct matrix_command *cmd = matrix_command_parse (s);
8869 if (c->n >= allocated)
8870 c->commands = x2nrealloc (c->commands, &allocated, sizeof *c->commands);
8871 c->commands[c->n++] = cmd;
8876 matrix_commands_uninit (struct matrix_commands *cmds)
8878 for (size_t i = 0; i < cmds->n; i++)
8879 matrix_command_destroy (cmds->commands[i]);
8880 free (cmds->commands);
8883 struct matrix_command_name
8886 struct matrix_command *(*parse) (struct matrix_state *);
8889 static const struct matrix_command_name *
8890 matrix_command_name_parse (struct lexer *lexer)
8892 static const struct matrix_command_name commands[] = {
8893 { "COMPUTE", matrix_compute_parse },
8894 { "CALL EIGEN", matrix_eigen_parse },
8895 { "CALL SETDIAG", matrix_setdiag_parse },
8896 { "CALL SVD", matrix_svd_parse },
8897 { "PRINT", matrix_print_parse },
8898 { "DO IF", matrix_do_if_parse },
8899 { "LOOP", matrix_loop_parse },
8900 { "BREAK", matrix_break_parse },
8901 { "READ", matrix_read_parse },
8902 { "WRITE", matrix_write_parse },
8903 { "GET", matrix_get_parse },
8904 { "SAVE", matrix_save_parse },
8905 { "MGET", matrix_mget_parse },
8906 { "MSAVE", matrix_msave_parse },
8907 { "DISPLAY", matrix_display_parse },
8908 { "RELEASE", matrix_release_parse },
8910 static size_t n = sizeof commands / sizeof *commands;
8912 for (const struct matrix_command_name *c = commands; c < &commands[n]; c++)
8913 if (lex_match_phrase (lexer, c->name))
8918 static struct matrix_command *
8919 matrix_command_parse (struct matrix_state *s)
8921 int start_ofs = lex_ofs (s->lexer);
8922 size_t nesting_level = SIZE_MAX;
8924 struct matrix_command *c = NULL;
8925 const struct matrix_command_name *cmd = matrix_command_name_parse (s->lexer);
8927 lex_error (s->lexer, _("Unknown matrix command."));
8928 else if (!cmd->parse)
8929 lex_error (s->lexer, _("Matrix command %s is not yet implemented."),
8933 nesting_level = output_open_group (
8934 group_item_create_nocopy (utf8_to_title (cmd->name),
8935 utf8_to_title (cmd->name)));
8941 c->location = lex_ofs_location (s->lexer, start_ofs, lex_ofs (s->lexer));
8942 msg_location_remove_columns (c->location);
8943 lex_end_of_command (s->lexer);
8945 lex_discard_rest_of_command (s->lexer);
8946 if (nesting_level != SIZE_MAX)
8947 output_close_groups (nesting_level);
8953 cmd_matrix (struct lexer *lexer, struct dataset *ds)
8955 if (!lex_force_match (lexer, T_ENDCMD))
8958 struct matrix_state state = {
8960 .session = dataset_session (ds),
8962 .vars = HMAP_INITIALIZER (state.vars),
8967 while (lex_match (lexer, T_ENDCMD))
8969 if (lex_token (lexer) == T_STOP)
8971 msg (SE, _("Unexpected end of input expecting matrix command."));
8975 if (lex_match_phrase (lexer, "END MATRIX"))
8978 struct matrix_command *c = matrix_command_parse (&state);
8981 matrix_command_execute (c);
8982 matrix_command_destroy (c);
8986 struct matrix_var *var, *next;
8987 HMAP_FOR_EACH_SAFE (var, next, struct matrix_var, hmap_node, &state.vars)
8990 gsl_matrix_free (var->value);
8991 hmap_delete (&state.vars, &var->hmap_node);
8994 hmap_destroy (&state.vars);
8995 msave_common_destroy (state.msave_common);
8996 fh_unref (state.prev_read_file);
8997 for (size_t i = 0; i < state.n_read_files; i++)
8998 read_file_destroy (state.read_files[i]);
8999 free (state.read_files);
9000 fh_unref (state.prev_write_file);
9001 for (size_t i = 0; i < state.n_write_files; i++)
9002 write_file_destroy (state.write_files[i]);
9003 free (state.write_files);
9004 fh_unref (state.prev_save_file);
9005 for (size_t i = 0; i < state.n_save_files; i++)
9006 save_file_destroy (state.save_files[i]);
9007 free (state.save_files);