X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Fmatrix.c;h=18fcfdd8197de1a373d263fe2aa86a35c3285867;hb=8d33786b32ea7f1d709c33d5a7c2c442fbeaafff;hp=a33e7693cd55dfb9558345ba6b33ddb1e603b1d7;hpb=d5c2af4608820668efdf4d40da5f5e597cebc05c;p=pspp diff --git a/src/language/stats/matrix.c b/src/language/stats/matrix.c index a33e7693cd..18fcfdd819 100644 --- a/src/language/stats/matrix.c +++ b/src/language/stats/matrix.c @@ -17,6 +17,7 @@ #include #include +#include #include #include #include @@ -53,12 +54,14 @@ #include "libpspp/string-array.h" #include "libpspp/stringi-set.h" #include "libpspp/u8-line.h" +#include "math/distributions.h" #include "math/random.h" #include "output/driver.h" #include "output/output-item.h" #include "output/pivot-table.h" #include "gl/c-ctype.h" +#include "gl/c-strcase.h" #include "gl/ftoastr.h" #include "gl/minmax.h" #include "gl/xsize.h" @@ -177,87 +180,187 @@ matrix_var_set (struct matrix_var *var, gsl_matrix *value) var->value = value; } -#define MATRIX_FUNCTIONS \ - F(ABS, "ABS", m_m) \ - F(ALL, "ALL", d_m) \ - F(ANY, "ANY", d_m) \ - F(ARSIN, "ARSIN", m_m) \ - F(ARTAN, "ARTAN", m_m) \ - F(BLOCK, "BLOCK", m_any) \ - F(CHOL, "CHOL", m_m) \ - F(CMIN, "CMIN", m_m) \ - F(CMAX, "CMAX", m_m) \ - F(COS, "COS", m_m) \ - F(CSSQ, "CSSQ", m_m) \ - F(CSUM, "CSUM", m_m) \ - F(DESIGN, "DESIGN", m_m) \ - F(DET, "DET", d_m) \ - F(DIAG, "DIAG", m_m) \ - F(EVAL, "EVAL", m_m) \ - F(EXP, "EXP", m_m) \ - F(GINV, "GINV", m_m) \ - F(GRADE, "GRADE", m_m) \ - F(GSCH, "GSCH", m_m) \ - F(IDENT, "IDENT", IDENT) \ - F(INV, "INV", m_m) \ - F(KRONEKER, "KRONEKER", m_mm) \ - F(LG10, "LG10", m_m) \ - F(LN, "LN", m_m) \ - F(MAGIC, "MAGIC", m_d) \ - F(MAKE, "MAKE", m_ddd) \ - F(MDIAG, "MDIAG", m_v) \ - F(MMAX, "MMAX", d_m) \ - F(MMIN, "MMIN", d_m) \ - F(MOD, "MOD", m_md) \ - F(MSSQ, "MSSQ", d_m) \ - F(MSUM, "MSUM", d_m) \ - F(NCOL, "NCOL", d_m) \ - F(NROW, "NROW", d_m) \ - F(RANK, "RANK", d_m) \ - F(RESHAPE, "RESHAPE", m_mdd) \ - F(RMAX, "RMAX", m_m) \ - F(RMIN, "RMIN", m_m) \ - F(RND, "RND", m_m) \ - F(RNKORDER, "RNKORDER", m_m) \ - F(RSSQ, "RSSQ", m_m) \ - F(RSUM, "RSUM", m_m) \ - F(SIN, "SIN", m_m) \ - F(SOLVE, "SOLVE", m_mm) \ - F(SQRT, "SQRT", m_m) \ - F(SSCP, "SSCP", m_m) \ - F(SVAL, "SVAL", m_m) \ - F(SWEEP, "SWEEP", m_md) \ - F(T, "T", m_m) \ - F(TRACE, "TRACE", d_m) \ - F(TRANSPOS, "TRANSPOS", m_m) \ - F(TRUNC, "TRUNC", m_m) \ - F(UNIFORM, "UNIFORM", m_dd) +#define MATRIX_FUNCTIONS \ + F(ABS, "ABS", m_m_e, NULL) \ + F(ALL, "ALL", d_m, NULL) \ + F(ANY, "ANY", d_m, NULL) \ + F(ARSIN, "ARSIN", m_m_e, "a[-1,1]") \ + F(ARTAN, "ARTAN", m_m_e, NULL) \ + F(BLOCK, "BLOCK", m_any, NULL) \ + F(CHOL, "CHOL", m_m, NULL) \ + F(CMIN, "CMIN", m_m, NULL) \ + F(CMAX, "CMAX", m_m, NULL) \ + F(COS, "COS", m_m_e, NULL) \ + F(CSSQ, "CSSQ", m_m, NULL) \ + F(CSUM, "CSUM", m_m, NULL) \ + F(DESIGN, "DESIGN", m_m, NULL) \ + F(DET, "DET", d_m, NULL) \ + F(DIAG, "DIAG", m_m, NULL) \ + F(EVAL, "EVAL", m_m, NULL) \ + F(EXP, "EXP", m_m_e, NULL) \ + F(GINV, "GINV", m_m, NULL) \ + F(GRADE, "GRADE", m_m, NULL) \ + F(GSCH, "GSCH", m_m, NULL) \ + F(IDENT, "IDENT", IDENT, NULL) \ + F(INV, "INV", m_m, NULL) \ + F(KRONEKER, "KRONEKER", m_mm, NULL) \ + F(LG10, "LG10", m_m_e, "a>0") \ + F(LN, "LN", m_m_e, "a>0") \ + F(MAGIC, "MAGIC", m_d, NULL) \ + F(MAKE, "MAKE", m_ddd, NULL) \ + F(MDIAG, "MDIAG", m_v, NULL) \ + F(MMAX, "MMAX", d_m, NULL) \ + F(MMIN, "MMIN", d_m, NULL) \ + F(MOD, "MOD", m_md, NULL) \ + F(MSSQ, "MSSQ", d_m, NULL) \ + F(MSUM, "MSUM", d_m, NULL) \ + F(NCOL, "NCOL", d_m, NULL) \ + F(NROW, "NROW", d_m, NULL) \ + F(RANK, "RANK", d_m, NULL) \ + F(RESHAPE, "RESHAPE", m_mdd, NULL) \ + F(RMAX, "RMAX", m_m, NULL) \ + F(RMIN, "RMIN", m_m, NULL) \ + F(RND, "RND", m_m_e, NULL) \ + F(RNKORDER, "RNKORDER", m_m, NULL) \ + F(RSSQ, "RSSQ", m_m, NULL) \ + F(RSUM, "RSUM", m_m, NULL) \ + F(SIN, "SIN", m_m_e, NULL) \ + F(SOLVE, "SOLVE", m_mm, NULL) \ + F(SQRT, "SQRT", m_m, NULL) \ + F(SSCP, "SSCP", m_m, NULL) \ + F(SVAL, "SVAL", m_m, NULL) \ + F(SWEEP, "SWEEP", m_md, NULL) \ + F(T, "T", m_m, NULL) \ + F(TRACE, "TRACE", d_m, NULL) \ + F(TRANSPOS, "TRANSPOS", m_m, NULL) \ + F(TRUNC, "TRUNC", m_m_e, NULL) \ + F(UNIFORM, "UNIFORM", m_dd, NULL) \ + F(PDF_BETA, "PDF.BETA", m_mdd_e, "a[0,1] b>0 c>0") \ + F(CDF_BETA, "CDF.BETA", m_mdd_e, "a[0,1] b>0 c>0") \ + F(IDF_BETA, "IDF.BETA", m_mdd_e, "a[0,1] b>0 c>0") \ + F(RV_BETA, "RV.BETA", d_dd, "a>0 b>0") \ + F(NCDF_BETA, "NCDF.BETA", m_mddd_e, "a>=0 b>0 c>0 d>0") \ + F(NPDF_BETA, "NCDF.BETA", m_mddd_e, "a>=0 b>0 c>0 d>0") \ + /* XXX CDF.BVNOR */ \ + F(PDF_BVNOR, "PDF.BVNOR", m_mdd_e, "c[-1,1]") \ + F(CDF_CAUCHY, "CDF.CAUCHY", m_mdd_e, "c>0") \ + F(IDF_CAUCHY, "IDF.CAUCHY", m_mdd_e, "a(0,1) c>0") \ + F(PDF_CAUCHY, "PDF.CAUCHY", m_mdd_e, "c>0") \ + F(RV_CAUCHY, "RV.CAUCHY", d_dd, "b>0") \ + F(CDF_CHISQ, "CDF.CHISQ", m_md_e, "a>=0 b>0") \ + F(IDF_CHISQ, "IDF.CHISQ", m_md_e, "a[0,1) b>0") \ + F(PDF_CHISQ, "PDF.CHISQ", m_md_e, "a>=0 b>0") \ + F(RV_CHISQ, "RV.CHISQ", d_d, "a>0") \ + F(SIG_CHISQ, "SIG.CHISQ", m_md_e, "a>=0 b>0") \ + F(CDF_EXP, "CDF.EXP", m_md_e, "a>=0 b>=0") \ + F(IDF_EXP, "IDF.EXP", m_md_e, "a[0,1) b>0") \ + F(PDF_EXP, "PDF.EXP", m_md_e, "a>=0 b>0") \ + F(RV_EXP, "RV.EXP", d_d, "a>0") \ + F(PDF_XPOWER, "PDF.XPOWER", m_mdd_e, "b>0 c>=0") \ + F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") + +struct matrix_function_properties + { + const char *name; + const char *constraints; + }; +/* +() +(a) +(a > 0) +(a >= 0) +(a > 0 && a < 1) +(a > 0 && a <= 1) +(a >= 0 && a <= 1) +(a > 0 && a < 1, b > 0) +(a >= 0 && a < 1, b > 0) +(a >= 0 && a <= 1, b > 0) +(a == 0 || a == 1, b >= 0 && b <= 1) +(a >= 0 && a < 1, b > 0, c > 0) +(a >= 0 && a <= 1, b > 0, c > 0) +(a >= 0 && a < 1, b >= 1, c >= 1) +(a >= 0 && a <= 1, b, c) +(a > 0 && a < 1, b, c > 0) +(a >= 0 && a <= 1, b <= c, c) +(a > 0 && a == floor (a), +(a >= 0 && a == floor (a) && a <= b, +(a >= 0 && a == floor (a) && a <= d, +(a >= 0 && a == floor (a), b > 0) +(a > 0 && a == floor (a), b >= 0 && b <= 1) +(a > 0, b > 0) +(a > 0, b >= 0) +(a >= 0, b > 0) +(a >= 0, b > 0, c) +(a > 0, b > 0, c > 0) +(a >= 0, b > 0, c > 0) +(a >= 0, b > 0, c > 0, d > 0) +(a >= 0, b > 0, c > 0, d >= 0) +(a > 0, b >= 1, c >= 1) +(a >= 1 && a == floor (a), b >= 0 && b <= 1) +(a >= 1, b > 0 && b <= 1) +(a >= 1, b == floor (b), c > 0 && c <= 1) +(a, b) +(a, b > 0) +(a, b > 0 && b <= 2) +(a, b > 0 && b <= 2, c >= -1 && c <= 1) +(a, b > 0 && b == floor (c), c >= 0 && c <= 1) +(a, b > 0, c) +(a, b > 0, c > 0) +(a, b > 0, c >= 0) +(a <= b, b) +(a >= b, b > 0, c > 0) +(a, b, c) +(a, b, c > 0) +(a, b, c >= -1 && c <= 1) +(a <= c, b <= a, c) +(a == floor (a), b > 0 && b <= 1) +b > 0 && b == floor (b), +b > 0 && b == floor (b) && b <= a, +c >= 0 && c <= 1) +c > 0 && c == floor (c) && c <= a) +c > 0 && c == floor (c) && c <= b, +d > 0 && d == floor (d) && d <= b) +*/ + +enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 }; +enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 }; enum { m_d_MIN_ARGS = 1, m_d_MAX_ARGS = 1 }; enum { m_dd_MIN_ARGS = 2, m_dd_MAX_ARGS = 2 }; enum { m_ddd_MIN_ARGS = 3, m_ddd_MAX_ARGS = 3 }; enum { m_m_MIN_ARGS = 1, m_m_MAX_ARGS = 1 }; +enum { m_m_e_MIN_ARGS = 1, m_m_e_MAX_ARGS = 1 }; enum { m_md_MIN_ARGS = 2, m_md_MAX_ARGS = 2 }; +enum { m_md_e_MIN_ARGS = 2, m_md_e_MAX_ARGS = 2 }; enum { m_mdd_MIN_ARGS = 3, m_mdd_MAX_ARGS = 3 }; +enum { m_mdd_e_MIN_ARGS = 3, m_mdd_e_MAX_ARGS = 3 }; +enum { m_mddd_e_MIN_ARGS = 4, m_mddd_e_MAX_ARGS = 4 }; enum { m_mm_MIN_ARGS = 2, m_mm_MAX_ARGS = 2 }; enum { m_v_MIN_ARGS = 1, m_v_MAX_ARGS = 1 }; enum { d_m_MIN_ARGS = 1, d_m_MAX_ARGS = 1 }; enum { m_any_MIN_ARGS = 1, m_any_MAX_ARGS = INT_MAX }; enum { IDENT_MIN_ARGS = 1, IDENT_MAX_ARGS = 2 }; +typedef double matrix_proto_d_d (double); +typedef double matrix_proto_d_dd (double, double); typedef gsl_matrix *matrix_proto_m_d (double); typedef gsl_matrix *matrix_proto_m_dd (double, double); typedef gsl_matrix *matrix_proto_m_ddd (double, double, double); typedef gsl_matrix *matrix_proto_m_m (gsl_matrix *); +typedef double matrix_proto_m_m_e (double); typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double); +typedef double matrix_proto_m_md_e (double, double); typedef gsl_matrix *matrix_proto_m_mdd (gsl_matrix *, double, double); +typedef double matrix_proto_m_mdd_e (double, double, double); +typedef double matrix_proto_m_mddd_e (double, double, double, double); typedef gsl_matrix *matrix_proto_m_mm (gsl_matrix *, gsl_matrix *); typedef gsl_matrix *matrix_proto_m_v (gsl_vector *); typedef double matrix_proto_d_m (gsl_matrix *); typedef gsl_matrix *matrix_proto_m_any (gsl_matrix *[], size_t n); typedef gsl_matrix *matrix_proto_IDENT (double, double); -#define F(ENUM, STRING, PROTO) static matrix_proto_##PROTO matrix_eval_##ENUM; +#define F(ENUM, STRING, PROTO, CONSTRAINTS) \ + static matrix_proto_##PROTO matrix_eval_##ENUM; MATRIX_FUNCTIONS #undef F @@ -268,7 +371,7 @@ struct matrix_expr enum matrix_op { /* Functions. */ -#define F(ENUM, STRING, PROTO) MOP_F_##ENUM, +#define F(ENUM, STRING, PROTO, CONSTRAINTS) MOP_F_##ENUM, MATRIX_FUNCTIONS #undef F @@ -343,7 +446,7 @@ matrix_expr_destroy (struct matrix_expr *e) switch (e->op) { -#define F(ENUM, STRING, PROTO) case MOP_F_##ENUM: +#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM: MATRIX_FUNCTIONS #undef F case MOP_NEGATE: @@ -502,12 +605,10 @@ to_vector (gsl_matrix *m) } -static gsl_matrix * -matrix_eval_ABS (gsl_matrix *m) +static double +matrix_eval_ABS (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = fabs (*d); - return m; + return fabs (d); } static double @@ -528,20 +629,16 @@ matrix_eval_ANY (gsl_matrix *m) return 0.0; } -static gsl_matrix * -matrix_eval_ARSIN (gsl_matrix *m) +static double +matrix_eval_ARSIN (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = asin (*d); - return m; + return asin (d); } -static gsl_matrix * -matrix_eval_ARTAN (gsl_matrix *m) +static double +matrix_eval_ARTAN (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = atan (*d); - return m; + return atan (d); } static gsl_matrix * @@ -623,12 +720,10 @@ matrix_eval_CMIN (gsl_matrix *m) return matrix_eval_col_extremum (m, true); } -static gsl_matrix * -matrix_eval_COS (gsl_matrix *m) +static double +matrix_eval_COS (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = cos (*d); - return m; + return cos (d); } static gsl_matrix * @@ -788,12 +883,10 @@ matrix_eval_EVAL (gsl_matrix *m) return eval; } -static gsl_matrix * -matrix_eval_EXP (gsl_matrix *m) +static double +matrix_eval_EXP (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = exp (*d); - return m; + return exp (d); } /* From https://gist.github.com/turingbirds/5e99656e08dbe1324c99, where it was @@ -1039,20 +1132,16 @@ matrix_eval_KRONEKER (gsl_matrix *a, gsl_matrix *b) return k; } -static gsl_matrix * -matrix_eval_LG10 (gsl_matrix *m) +static double +matrix_eval_LG10 (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = log10 (*d); - return m; + return log10 (d); } -static gsl_matrix * -matrix_eval_LN (gsl_matrix *m) +static double +matrix_eval_LN (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = log (*d); - return m; + return log (d); } static void @@ -1378,12 +1467,10 @@ matrix_eval_RMIN (gsl_matrix *m) return matrix_eval_row_extremum (m, true); } -static gsl_matrix * -matrix_eval_RND (gsl_matrix *m) +static double +matrix_eval_RND (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = rint (*d); - return m; + return rint (d); } struct rank @@ -1464,12 +1551,10 @@ matrix_eval_RSUM (gsl_matrix *m) return matrix_eval_row_sum (m, false); } -static gsl_matrix * -matrix_eval_SIN (gsl_matrix *m) +static double +matrix_eval_SIN (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = sin (*d); - return m; + return sin (d); } static gsl_matrix * @@ -1627,12 +1712,10 @@ matrix_eval_TRANSPOS (gsl_matrix *m) } } -static gsl_matrix * -matrix_eval_TRUNC (gsl_matrix *m) +static double +matrix_eval_TRUNC (double d) { - MATRIX_FOR_ALL_ELEMENTS (d, y, x, m) - *d = trunc (*d); - return m; + return trunc (d); } static gsl_matrix * @@ -1657,6 +1740,138 @@ matrix_eval_UNIFORM (double r_, double c_) return m; } +static double +matrix_eval_PDF_BETA (double x, double a, double b) +{ + return gsl_ran_beta_pdf (x, a, b); +} + +static double +matrix_eval_CDF_BETA (double x, double a, double b) +{ + return gsl_cdf_beta_P (x, a, b); +} + +static double +matrix_eval_IDF_BETA (double P, double a, double b) +{ + return gsl_cdf_beta_Pinv (P, a, b); +} + +static double +matrix_eval_RV_BETA (double a, double b) +{ + return gsl_ran_beta (get_rng (), a, b); +} + +static double +matrix_eval_NCDF_BETA (double x, double a, double b, double lambda) +{ + return ncdf_beta (x, a, b, lambda); +} + +static double +matrix_eval_NPDF_BETA (double x, double a, double b, double lambda) +{ + return npdf_beta (x, a, b, lambda); +} + +static double +matrix_eval_PDF_BVNOR (double x0, double x1, double r) +{ + return gsl_ran_bivariate_gaussian_pdf (x0, x1, 1, 1, r); +} + +static double +matrix_eval_CDF_CAUCHY (double x, double a, double b) +{ + return gsl_cdf_cauchy_P ((x - a) / b, 1); +} + +static double +matrix_eval_IDF_CAUCHY (double P, double a, double b) +{ + return a + b * gsl_cdf_cauchy_Pinv (P, 1); +} + +static double +matrix_eval_PDF_CAUCHY (double x, double a, double b) +{ + return gsl_ran_cauchy_pdf ((x - a) / b, 1) / b; +} + +static double +matrix_eval_RV_CAUCHY (double a, double b) +{ + return a + b * gsl_ran_cauchy (get_rng (), 1); +} + +static double +matrix_eval_CDF_CHISQ (double x, double df) +{ + return gsl_cdf_chisq_P (x, df); +} + +static double +matrix_eval_IDF_CHISQ (double P, double df) +{ + return gsl_cdf_chisq_Pinv (P, df); +} + +static double +matrix_eval_PDF_CHISQ (double x, double df) +{ + return gsl_ran_chisq_pdf (x, df); +} + +static double +matrix_eval_RV_CHISQ (double df) +{ + return gsl_ran_chisq (get_rng (), df); +} + +static double +matrix_eval_SIG_CHISQ (double x, double df) +{ + return gsl_cdf_chisq_Q (x, df); +} + +static double +matrix_eval_CDF_EXP (double x, double a) +{ + return gsl_cdf_exponential_P (x, 1. / a); +} + +static double +matrix_eval_IDF_EXP (double P, double a) +{ + return gsl_cdf_exponential_Pinv (P, 1. / a); +} + +static double +matrix_eval_PDF_EXP (double x, double a) +{ + return gsl_ran_exponential_pdf (x, 1. / a); +} + +static double +matrix_eval_RV_EXP (double a) +{ + return gsl_ran_exponential (get_rng (), 1. / a); +} + +static double +matrix_eval_PDF_XPOWER (double x, double a, double b) +{ + return gsl_ran_exppow_pdf (x, a, b); +} + +static double +matrix_eval_RV_XPOWER (double a, double b) +{ + return gsl_ran_exppow (get_rng (), a, b); +} + struct matrix_function { const char *name; @@ -1666,13 +1881,58 @@ struct matrix_function static struct matrix_expr *matrix_parse_expr (struct matrix_state *); +static bool +word_matches (const char **test, const char **name) +{ + size_t test_len = strcspn (*test, "."); + size_t name_len = strcspn (*name, "."); + if (test_len == name_len) + { + if (buf_compare_case (*test, *name, test_len)) + return false; + } + else if (test_len < 3 || test_len > name_len) + return false; + else + { + if (buf_compare_case (*test, *name, test_len)) + return false; + } + + *test += test_len; + *name += name_len; + if (**test != **name) + return false; + + if (**test == '.') + { + (*test)++; + (*name)++; + } + return true; +} + +/* Returns 0 if TOKEN and FUNC do not match, + 1 if TOKEN is an acceptable abbreviation for FUNC, + 2 if TOKEN equals FUNC. */ +static int +compare_function_names (const char *token_, const char *func_) +{ + const char *token = token_; + const char *func = func_; + while (*token || *func) + if (!word_matches (&token, &func)) + return 0; + return !c_strcasecmp (token_, func_) ? 2 : 1; +} + static const struct matrix_function * matrix_parse_function_name (const char *token) { static const struct matrix_function functions[] = { -#define F(ENUM, STRING, PROTO) \ - { #ENUM, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS }, +#define F(ENUM, STRING, PROTO, CONSTRAINTS) \ + { STRING, MOP_F_##ENUM, PROTO##_MIN_ARGS, PROTO##_MAX_ARGS }, MATRIX_FUNCTIONS #undef F }; @@ -1680,7 +1940,7 @@ matrix_parse_function_name (const char *token) for (size_t i = 0; i < N_FUNCTIONS; i++) { - if (lex_id_match_n (ss_cstr (functions[i].name), ss_cstr (token), 3)) + if (compare_function_names (token, functions[i].name) > 0) return &functions[i]; } return NULL; @@ -2217,7 +2477,7 @@ matrix_op_eval (enum matrix_op op, double a, double b) case MOP_OR: return (a > 0) || (b > 0); case MOP_XOR: return (a > 0) != (b > 0); -#define F(ENUM, STRING, PROTO) case MOP_F_##ENUM: +#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM: MATRIX_FUNCTIONS #undef F case MOP_NEGATE: @@ -2262,7 +2522,7 @@ matrix_op_name (enum matrix_op op) case MOP_OR: return "OR"; case MOP_XOR: return "XOR"; -#define F(ENUM, STRING, PROTO) case MOP_F_##ENUM: +#define F(ENUM, STRING, PROTO, CONSTRAINTS) case MOP_F_##ENUM: MATRIX_FUNCTIONS #undef F case MOP_NEGATE: @@ -2720,9 +2980,9 @@ matrix_expr_evaluate_mat_index (gsl_matrix *sm, gsl_matrix *im0, return dm; } -#define F(ENUM, STRING, PROTO) \ - static gsl_matrix *matrix_expr_evaluate_##PROTO ( \ - const char *name, gsl_matrix *[], size_t, \ +#define F(ENUM, STRING, PROTO, CONSTRAINTS) \ + static gsl_matrix *matrix_expr_evaluate_##PROTO ( \ + const struct matrix_function_properties *, gsl_matrix *[], size_t, \ matrix_proto_##PROTO *); MATRIX_FUNCTIONS #undef F @@ -2765,41 +3025,194 @@ to_scalar_args (const char *name, gsl_matrix *subs[], size_t n_subs, double d[]) return true; } +static int +parse_constraint_value (const char **constraintsp) +{ + char *tail; + long retval = strtol (*constraintsp, &tail, 10); + *constraintsp = tail; + return retval; +} + +static bool +check_constraints (const struct matrix_function_properties *props, + gsl_matrix *args[], size_t n_args) +{ + const char *constraints = props->constraints; + if (!constraints) + return true; + + size_t arg_index = SIZE_MAX; + while (*constraints) + { + if (*constraints >= 'a' && *constraints <= 'd') + { + arg_index = *constraints++ - 'a'; + assert (arg_index < n_args); + } + else if (*constraints == '[' || *constraints == '(') + { + assert (arg_index < n_args); + bool open_lower = *constraints++ == '('; + int minimum = parse_constraint_value (&constraints); + assert (*constraints == ','); + constraints++; + int maximum = parse_constraint_value (&constraints); + assert (*constraints == ']' || *constraints == ')'); + bool open_upper = *constraints++ == ')'; + + MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index]) + if ((open_lower ? *d <= minimum : *d < minimum) + || (open_upper ? *d >= maximum : *d > maximum)) + { + if (!is_scalar (args[arg_index])) + msg (ME, _("Row %zu, column %zu of argument %zu to matrix " + "function %s has value %g, which is outside " + "the valid range %c%d,%d%c."), + y + 1, x + 1, arg_index + 1, props->name, *d, + open_lower ? '(' : '[', + minimum, maximum, + open_upper ? ')' : ']'); + else + msg (ME, _("Argument %zu to matrix function %s has value %g, " + "which is outside the valid range %c%d,%d%c."), + arg_index + 1, props->name, *d, + open_lower ? '(' : '[', + minimum, maximum, + open_upper ? ')' : ']'); + return false; + } + } + else if (*constraints == '>' || *constraints == '<') + { + bool greater = *constraints++ == '>'; + bool equal; + if (*constraints == '=') + { + equal = true; + constraints++; + } + else + equal = false; + int comparand = parse_constraint_value (&constraints); + + assert (arg_index < n_args); + MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index]) + if (greater + ? (equal ? !(*d >= comparand) : !(*d > comparand)) + : (equal ? !(*d <= comparand) : !(*d < comparand))) + { + struct string s = DS_EMPTY_INITIALIZER; + if (!is_scalar (args[arg_index])) + ds_put_format (&s, _("Row %zu, column %zu of argument %zu " + "to matrix function %s has " + "invalid value %g."), + y + 1, x + 1, arg_index + 1, props->name, *d); + else + ds_put_format (&s, _("Argument %zu to matrix function %s " + "has invalid value %g."), + arg_index + 1, props->name, *d); + + ds_put_cstr (&s, " "); + if (greater && equal) + ds_put_format (&s, _("This argument must be greater than or " + "equal to %d."), comparand); + else if (greater && !equal) + ds_put_format (&s, _("This argument must be greater than %d."), + comparand); + else if (equal) + ds_put_format (&s, _("This argument must be less than or " + "equal to %d."), comparand); + else + ds_put_format (&s, _("This argument must be less than %d."), + comparand); + msg (ME, ds_cstr (&s)); + ds_destroy (&s); + return false; + } + } + else + { + assert (*constraints == ' '); + constraints++; + } + } + return true; +} + static gsl_matrix * -matrix_expr_evaluate_m_d (const char *name, +matrix_expr_evaluate_d_d (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_d_d *f) +{ + assert (n_subs == 1); + + double d; + if (!to_scalar_args (props->name, subs, n_subs, &d)) + return NULL; + + if (!check_constraints (props, subs, n_subs)) + return NULL; + + gsl_matrix *m = gsl_matrix_alloc (1, 1); + gsl_matrix_set (m, 0, 0, f (d)); + return m; +} + +static gsl_matrix * +matrix_expr_evaluate_d_dd (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_d_dd *f) +{ + assert (n_subs == 2); + + double d[2]; + if (!to_scalar_args (props->name, subs, n_subs, d)) + return NULL; + + if (!check_constraints (props, subs, n_subs)) + return NULL; + + gsl_matrix *m = gsl_matrix_alloc (1, 1); + gsl_matrix_set (m, 0, 0, f (d[0], d[1])); + return m; +} + +static gsl_matrix * +matrix_expr_evaluate_m_d (const struct matrix_function_properties *props, gsl_matrix *subs[], size_t n_subs, matrix_proto_m_d *f) { assert (n_subs == 1); double d; - return to_scalar_args (name, subs, n_subs, &d) ? f(d) : NULL; + return to_scalar_args (props->name, subs, n_subs, &d) ? f (d) : NULL; } static gsl_matrix * -matrix_expr_evaluate_m_dd (const char *name, +matrix_expr_evaluate_m_dd (const struct matrix_function_properties *props, gsl_matrix *subs[], size_t n_subs, matrix_proto_m_dd *f) { assert (n_subs == 2); double d[2]; - return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1]) : NULL; + return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1]) : NULL; } static gsl_matrix * -matrix_expr_evaluate_m_ddd (const char *name, +matrix_expr_evaluate_m_ddd (const struct matrix_function_properties *props, gsl_matrix *subs[], size_t n_subs, matrix_proto_m_ddd *f) { assert (n_subs == 3); double d[3]; - return to_scalar_args (name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL; + return to_scalar_args (props->name, subs, n_subs, d) ? f(d[0], d[1], d[2]) : NULL; } static gsl_matrix * -matrix_expr_evaluate_m_m (const char *name UNUSED, +matrix_expr_evaluate_m_m (const struct matrix_function_properties *props UNUSED, gsl_matrix *subs[], size_t n_subs, matrix_proto_m_m *f) { @@ -2808,29 +3221,102 @@ matrix_expr_evaluate_m_m (const char *name UNUSED, } static gsl_matrix * -matrix_expr_evaluate_m_md (const char *name UNUSED, +matrix_expr_evaluate_m_m_e (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_m_e *f) +{ + assert (n_subs == 1); + + if (!check_constraints (props, subs, n_subs)) + return NULL; + + MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0]) + *a = f (*a); + return subs[0]; +} + +static gsl_matrix * +matrix_expr_evaluate_m_md (const struct matrix_function_properties *props UNUSED, gsl_matrix *subs[], size_t n_subs, matrix_proto_m_md *f) { assert (n_subs == 2); - if (!check_scalar_arg (name, subs, 1)) + if (!check_scalar_arg (props->name, subs, 1)) return NULL; return f (subs[0], to_scalar (subs[1])); } static gsl_matrix * -matrix_expr_evaluate_m_mdd (const char *name UNUSED, +matrix_expr_evaluate_m_md_e (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_md_e *f) +{ + assert (n_subs == 2); + if (!check_scalar_arg (props->name, subs, 1)) + return NULL; + + if (!check_constraints (props, subs, n_subs)) + return NULL; + + double b = to_scalar (subs[1]); + MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0]) + *a = f (*a, b); + return subs[0]; +} + +static gsl_matrix * +matrix_expr_evaluate_m_mdd (const struct matrix_function_properties *props UNUSED, gsl_matrix *subs[], size_t n_subs, matrix_proto_m_mdd *f) { assert (n_subs == 3); - if (!check_scalar_arg (name, subs, 1) || !check_scalar_arg (name, subs, 2)) + if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2)) return NULL; return f (subs[0], to_scalar (subs[1]), to_scalar (subs[2])); } static gsl_matrix * -matrix_expr_evaluate_m_mm (const char *name UNUSED, +matrix_expr_evaluate_m_mdd_e (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_mdd_e *f) +{ + assert (n_subs == 3); + if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2)) + return NULL; + + if (!check_constraints (props, subs, n_subs)) + return NULL; + + double b = to_scalar (subs[1]); + double c = to_scalar (subs[2]); + MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0]) + *a = f (*a, b, c); + return subs[0]; +} + +static gsl_matrix * +matrix_expr_evaluate_m_mddd_e (const struct matrix_function_properties *props, + gsl_matrix *subs[], size_t n_subs, + matrix_proto_m_mddd_e *f) +{ + assert (n_subs == 4); + for (size_t i = 1; i < 4; i++) + if (!check_scalar_arg (props->name, subs, i)) + return NULL; + + if (!check_constraints (props, subs, n_subs)) + return NULL; + + double b = to_scalar (subs[1]); + double c = to_scalar (subs[2]); + double d = to_scalar (subs[3]); + MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0]) + *a = f (*a, b, c, d); + return subs[0]; +} + +static gsl_matrix * +matrix_expr_evaluate_m_mm (const struct matrix_function_properties *props UNUSED, gsl_matrix *subs[], size_t n_subs, matrix_proto_m_mm *f) { @@ -2839,19 +3325,19 @@ matrix_expr_evaluate_m_mm (const char *name UNUSED, } static gsl_matrix * -matrix_expr_evaluate_m_v (const char *name, +matrix_expr_evaluate_m_v (const struct matrix_function_properties *props, gsl_matrix *subs[], size_t n_subs, matrix_proto_m_v *f) { assert (n_subs == 1); - if (!check_vector_arg (name, subs, 0)) + if (!check_vector_arg (props->name, subs, 0)) return NULL; gsl_vector v = to_vector (subs[0]); return f (&v); } static gsl_matrix * -matrix_expr_evaluate_d_m (const char *name UNUSED, +matrix_expr_evaluate_d_m (const struct matrix_function_properties *props UNUSED, gsl_matrix *subs[], size_t n_subs, matrix_proto_d_m *f) { @@ -2863,7 +3349,7 @@ matrix_expr_evaluate_d_m (const char *name UNUSED, } static gsl_matrix * -matrix_expr_evaluate_m_any (const char *name UNUSED, +matrix_expr_evaluate_m_any (const struct matrix_function_properties *props UNUSED, gsl_matrix *subs[], size_t n_subs, matrix_proto_m_any *f) { @@ -2871,14 +3357,14 @@ matrix_expr_evaluate_m_any (const char *name UNUSED, } static gsl_matrix * -matrix_expr_evaluate_IDENT (const char *name, +matrix_expr_evaluate_IDENT (const struct matrix_function_properties *props, gsl_matrix *subs[], size_t n_subs, matrix_proto_IDENT *f) { assert (n_subs <= 2); double d[2]; - if (!to_scalar_args (name, subs, n_subs, d)) + if (!to_scalar_args (props->name, subs, n_subs, d)) return NULL; return f (d[0], d[n_subs - 1]); @@ -2937,10 +3423,16 @@ matrix_expr_evaluate (const struct matrix_expr *e) gsl_matrix *result = NULL; switch (e->op) { -#define F(ENUM, STRING, PROTO) \ +#define F(ENUM, STRING, PROTO, CONSTRAINTS) \ case MOP_F_##ENUM: \ - result = matrix_expr_evaluate_##PROTO (STRING, subs, e->n_subs, \ - matrix_eval_##ENUM); \ + { \ + static const struct matrix_function_properties props = { \ + .name = STRING, \ + .constraints = CONSTRAINTS, \ + }; \ + result = matrix_expr_evaluate_##PROTO (&props, subs, e->n_subs, \ + matrix_eval_##ENUM); \ + } \ break; MATRIX_FUNCTIONS #undef F