}
\f
#define MATRIX_FUNCTIONS \
- F(ABS, "ABS", m_m_e, NULL) \
+ F(ABS, "ABS", 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(ARSIN, "ARSIN", m_e, "a[-1,1]") \
+ F(ARTAN, "ARTAN", 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(COS, "COS", 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(EXP, "EXP", 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(LG10, "LG10", m_e, "a>0") \
+ F(LN, "LN", m_e, "a>0") \
F(MAGIC, "MAGIC", m_d, NULL) \
F(MAKE, "MAKE", m_ddd, NULL) \
F(MDIAG, "MDIAG", m_v, 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(RND, "RND", 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(SIN, "SIN", m_e, NULL) \
F(SOLVE, "SOLVE", m_mm, NULL) \
F(SQRT, "SQRT", m_m, NULL) \
F(SSCP, "SSCP", m_m, 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(TRUNC, "TRUNC", 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(PDF_BETA, "PDF.BETA", m_edd, "a[0,1] b>0 c>0") \
+ F(CDF_BETA, "CDF.BETA", m_edd, "a[0,1] b>0 c>0") \
+ F(IDF_BETA, "IDF.BETA", m_edd, "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(NCDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
+ F(NPDF_BETA, "NCDF.BETA", m_eddd, "a>=0 b>0 c>0 d>0") \
+ F(CDF_BVNOR, "CDF.BVNOR", m_eed, "c[-1,1]") \
+ F(PDF_BVNOR, "PDF.BVNOR", m_eed, "c[-1,1]") \
+ F(CDF_CAUCHY, "CDF.CAUCHY", m_edd, "c>0") \
+ F(IDF_CAUCHY, "IDF.CAUCHY", m_edd, "a(0,1) c>0") \
+ F(PDF_CAUCHY, "PDF.CAUCHY", m_edd, "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(CDF_CHISQ, "CDF.CHISQ", m_ed, "a>=0 b>0") \
+ F(CHICDF, "CHICDF", m_ed, "a>=0 b>0") \
+ F(IDF_CHISQ, "IDF.CHISQ", m_ed, "a[0,1) b>0") \
+ F(PDF_CHISQ, "PDF.CHISQ", m_ed, "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(SIG_CHISQ, "SIG.CHISQ", m_ed, "a>=0 b>0") \
+ F(CDF_EXP, "CDF.EXP", m_ed, "a>=0 b>=0") \
+ F(IDF_EXP, "IDF.EXP", m_ed, "a[0,1) b>0") \
+ F(PDF_EXP, "PDF.EXP", m_ed, "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(PDF_XPOWER, "PDF.XPOWER", m_edd, "b>0 c>=0") \
F(RV_XPOWER, "RV.XPOWER", d_dd, "a>0 c>=0") \
- F(CDF_F, "CDF.F", m_mdd_e, "a>=0 b>0 c>0") \
- F(IDF_F, "IDF.F", m_mdd_e, "a[0,1) b>0 c>0") \
- F(PDF_F, "PDF.F", m_mdd_e, "a>=0 b>0 c>0") \
+ F(CDF_F, "CDF.F", m_edd, "a>=0 b>0 c>0") \
+ F(FCDF, "FCDF", m_edd, "a>=0 b>0 c>0") \
+ F(IDF_F, "IDF.F", m_edd, "a[0,1) b>0 c>0") \
+ F(PDF_F, "PDF.F", m_edd, "a>=0 b>0 c>0") \
F(RV_F, "RV.F", d_dd, "a>0 b>0") \
- F(SIG_F, "SIG.F", m_mdd_e, "a>=0 b>0 c>0")
+ F(SIG_F, "SIG.F", m_edd, "a>=0 b>0 c>0") \
+ F(CDF_GAMMA, "CDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
+ F(IDF_GAMMA, "IDF.GAMMA", m_edd, "a[0,1] b>0 c>0") \
+ F(PDF_GAMMA, "PDF.GAMMA", m_edd, "a>=0 b>0 c>0") \
+ F(RV_GAMMA, "RV.GAMMA", d_dd, "a>0 b>0") \
+ F(PDF_LANDAU, "PDF.LANDAU", m_e, NULL) \
+ F(RV_LANDAU, "RV.LANDAU", d_none, NULL) \
+ F(CDF_LAPLACE, "CDF.LAPLACE", m_edd, "c>0") \
+ F(IDF_LAPLACE, "IDF.LAPLACE", m_edd, "a(0,1) c>0") \
+ F(PDF_LAPLACE, "PDF.LAPLACE", m_edd, "c>0") \
+ F(RV_LAPLACE, "RV.LAPLACE", d_dd, "b>0") \
+ F(RV_LEVY, "RV.LEVY", d_dd, "b(0,2]") \
+ F(RV_LVSKEW, "RV.LVSKEW", d_ddd, "b(0,2] c[-1,1]") \
+ F(CDF_LOGISTIC, "CDF.LOGISTIC", m_edd, "c>0") \
+ F(IDF_LOGISTIC, "IDF.LOGISTIC", m_edd, "a(0,1) c>0") \
+ F(PDF_LOGISTIC, "PDF.LOGISTIC", m_edd, "c>0") \
+ F(RV_LOGISTIC, "RV.LOGISTIC", d_dd, "b>0") \
+ F(CDF_LNORMAL, "CDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
+ F(IDF_LNORMAL, "IDF.LNORMAL", m_edd, "a[0,1) b>0 c>0") \
+ F(PDF_LNORMAL, "PDF.LNORMAL", m_edd, "a>=0 b>0 c>0") \
+ F(RV_LNORMAL, "RV.LNORMAL", d_dd, "a>0 b>0") \
+ F(CDF_NORMAL, "CDF.NORMAL", m_edd, "c>0") \
+ F(IDF_NORMAL, "IDF.NORMAL", m_edd, "a(0,1) c>0") \
+ F(PDF_NORMAL, "PDF.NORMAL", m_edd, "c>0") \
+ F(RV_NORMAL, "RV.NORMAL", d_dd, "b>0") \
+ F(CDFNORM, "CDFNORM", m_e, NULL) \
+ F(PROBIT, "PROBIT", m_e, "a(0,1)") \
+ F(NORMAL, "NORMAL", m_e, "a>0") \
+ F(PDF_NTAIL, "PDF.NTAIL", m_edd, "b>0 c>0") \
+ F(RV_NTAIL, "RV.NTAIL", d_dd, "a>0 b>0") \
+ F(CDF_PARETO, "CDF.PARETO", m_edd, "a>=b b>0 c>0") \
+ F(IDF_PARETO, "IDF.PARETO", m_edd, "a[0,1) b>0 c>0") \
+ F(PDF_PARETO, "PDF.PARETO", m_edd, "a>=b b>0 c>0") \
+ F(RV_PARETO, "RV.PARETO", d_dd, "a>0 b>0") \
+ F(CDF_RAYLEIGH, "CDF.RAYLEIGH", m_ed, "b>0") \
+ F(IDF_RAYLEIGH, "IDF.RAYLEIGH", m_ed, "a[0,1] b>0") \
+ F(PDF_RAYLEIGH, "PDF.RAYLEIGH", m_ed, "b>0") \
+ F(RV_RAYLEIGH, "RV.RAYLEIGH", d_d, "a>0") \
+ F(PDF_RTAIL, "PDF.RTAIL", m_edd, NULL) \
+ F(RV_RTAIL, "RV.RTAIL", d_dd, NULL) \
+ F(CDF_T, "CDF.T", m_ed, "b>0") \
+ F(TCDF, "TCDF", m_ed, "b>0") \
+ F(IDF_T, "IDF.T", m_ed, "a(0,1) b>0") \
+ F(PDF_T, "PDF.T", m_ed, "b>0") \
+ F(RV_T, "RV.T", d_d, "a>0") \
+ F(CDF_T1G, "CDF.T1G", m_edd, NULL) \
+ F(IDF_T1G, "IDF.T1G", m_edd, "a(0,1)") \
+ F(PDF_T1G, "PDF.T1G", m_edd, NULL) \
+ F(RV_T1G, "RV.T1G", d_dd, NULL) \
+ F(CDF_T2G, "CDF.T2G", m_edd, NULL) \
+ F(IDF_T2G, "IDF.T2G", m_edd, "a(0,1)") \
+ F(PDF_T2G, "PDF.T2G", m_edd, NULL) \
+ F(RV_T2G, "RV.T2G", d_dd, NULL) \
+ F(CDF_UNIFORM, "CDF.UNIFORM", m_edd, "a<=c b<=c") \
+ F(IDF_UNIFORM, "IDF.UNIFORM", m_edd, "a[0,1] b<=c") \
+ F(PDF_UNIFORM, "PDF.UNIFORM", m_edd, "a<=c b<=c") \
+ F(RV_UNIFORM, "RV.UNIFORM", d_dd, "a<=b") \
+ F(CDF_WEIBULL, "CDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
+ F(IDF_WEIBULL, "IDF.WEIBULL", m_edd, "a[0,1) b>0 c>0") \
+ F(PDF_WEIBULL, "PDF.WEIBULL", m_edd, "a>=0 b>0 c>0") \
+ F(RV_WEIBULL, "RV.WEIBULL", d_dd, "a>0 b>0") \
+ F(CDF_BERNOULLI, "CDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
+ F(PDF_BERNOULLI, "PDF.BERNOULLI", m_ed, "ai[0,1] b[0,1]") \
+ F(RV_BERNOULLI, "RV.BERNOULLI", d_d, "a[0,1]") \
+ F(CDF_BINOM, "CDF.BINOM", m_edd, "bi>0 c[0,1]") \
+ F(PDF_BINOM, "PDF.BINOM", m_edd, "ai>=0<=b bi>0 c[0,1]") \
+ F(RV_BINOM, "RV.BINOM", d_dd, "ai>0 b[0,1]") \
+ F(CDF_GEOM, "CDF.GEOM", m_ed, "ai>=1 b[0,1]") \
+ F(PDF_GEOM, "PDF.GEOM", m_ed, "ai>=1 b[0,1]") \
+ F(RV_GEOM, "RV.GEOM", d_d, "a[0,1]") \
+ F(CDF_HYPER, "CDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
+ F(PDF_HYPER, "PDF.HYPER", m_eddd, "ai>=0<=d bi>0 ci>0<=b di>0<=b") \
+ F(RV_HYPER, "RV.HYPER", d_ddd, "ai>0 bi>0<=a ci>0<=a") \
+ F(PDF_LOG, "PDF.LOG", m_ed, "a>=1 b(0,1]") \
+ F(RV_LOG, "RV.LOG", d_d, "a(0,1]") \
+ F(CDF_NEGBIN, "CDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
+ F(PDF_NEGBIN, "PDF.NEGBIN", m_edd, "a>=1 bi c(0,1]") \
+ F(RV_NEGBIN, "RV.NEGBIN", d_dd, "ai b(0,1]") \
+ F(CDF_POISSON, "CDF.POISSON", m_ed, "ai>=0 b>0") \
+ F(PDF_POISSON, "PDF.POISSON", m_ed, "ai>=0 b>0") \
+ F(RV_POISSON, "RV.POISSON", d_d, "a>0")
struct matrix_function_properties
{
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_none_MIN_ARGS = 0, d_none_MAX_ARGS = 0 };
enum { d_d_MIN_ARGS = 1, d_d_MAX_ARGS = 1 };
enum { d_dd_MIN_ARGS = 2, d_dd_MAX_ARGS = 2 };
+enum { d_ddd_MIN_ARGS = 3, d_ddd_MAX_ARGS = 3 };
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_e_MIN_ARGS = 1, 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_ed_MIN_ARGS = 2, m_ed_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_edd_MIN_ARGS = 3, m_edd_MAX_ARGS = 3 };
+enum { m_eddd_MIN_ARGS = 4, m_eddd_MAX_ARGS = 4 };
+enum { m_eed_MIN_ARGS = 3, m_eed_MAX_ARGS = 3 };
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_none (void);
typedef double matrix_proto_d_d (double);
typedef double matrix_proto_d_dd (double, double);
+typedef double matrix_proto_d_dd (double, double);
+typedef double matrix_proto_d_ddd (double, 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 double matrix_proto_m_e (double);
typedef gsl_matrix *matrix_proto_m_md (gsl_matrix *, double);
-typedef double matrix_proto_m_md_e (double, double);
+typedef double matrix_proto_m_ed (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 double matrix_proto_m_edd (double, double, double);
+typedef double matrix_proto_m_eddd (double, double, double, double);
+typedef double matrix_proto_m_eed (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 *);
return npdf_beta (x, a, b, lambda);
}
+static double
+matrix_eval_CDF_BVNOR (double x0, double x1, double r)
+{
+ return cdf_bvnor (x0, x1, r);
+}
+
static double
matrix_eval_PDF_BVNOR (double x0, double x1, double r)
{
return gsl_cdf_chisq_P (x, df);
}
+static double
+matrix_eval_CHICDF (double x, double df)
+{
+ return matrix_eval_CDF_CHISQ (x, df);
+}
+
static double
matrix_eval_IDF_CHISQ (double P, double df)
{
return gsl_cdf_fdist_P (x, df1, df2);
}
+static double
+matrix_eval_FCDF (double x, double df1, double df2)
+{
+ return matrix_eval_CDF_F (x, df1, df2);
+}
+
static double
matrix_eval_IDF_F (double P, double df1, double df2)
{
return gsl_cdf_fdist_Q (x, df1, df2);
}
+static double
+matrix_eval_CDF_GAMMA (double x, double a, double b)
+{
+ return gsl_cdf_gamma_P (x, a, 1. / b);
+}
+
+static double
+matrix_eval_IDF_GAMMA (double P, double a, double b)
+{
+ return gsl_cdf_gamma_Pinv (P, a, 1. / b);
+}
+
+static double
+matrix_eval_PDF_GAMMA (double x, double a, double b)
+{
+ return gsl_ran_gamma_pdf (x, a, 1. / b);
+}
+
+static double
+matrix_eval_RV_GAMMA (double a, double b)
+{
+ return gsl_ran_gamma (get_rng (), a, 1. / b);
+}
+
+static double
+matrix_eval_PDF_LANDAU (double x)
+{
+ return gsl_ran_landau_pdf (x);
+}
+
+static double
+matrix_eval_RV_LANDAU (void)
+{
+ return gsl_ran_landau (get_rng ());
+}
+
+static double
+matrix_eval_CDF_LAPLACE (double x, double a, double b)
+{
+ return gsl_cdf_laplace_P ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_IDF_LAPLACE (double P, double a, double b)
+{
+ return a + b * gsl_cdf_laplace_Pinv (P, 1);
+}
+
+static double
+matrix_eval_PDF_LAPLACE (double x, double a, double b)
+{
+ return gsl_ran_laplace_pdf ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_RV_LAPLACE (double a, double b)
+{
+ return a + b * gsl_ran_laplace (get_rng (), 1);
+}
+
+static double
+matrix_eval_RV_LEVY (double c, double alpha)
+{
+ return gsl_ran_levy (get_rng (), c, alpha);
+}
+
+static double
+matrix_eval_RV_LVSKEW (double c, double alpha, double beta)
+{
+ return gsl_ran_levy_skew (get_rng (), c, alpha, beta);
+}
+
+static double
+matrix_eval_CDF_LOGISTIC (double x, double a, double b)
+{
+ return gsl_cdf_logistic_P ((x - a) / b, 1);
+}
+
+static double
+matrix_eval_IDF_LOGISTIC (double P, double a, double b)
+{
+ return a + b * gsl_cdf_logistic_Pinv (P, 1);
+}
+
+static double
+matrix_eval_PDF_LOGISTIC (double x, double a, double b)
+{
+ return gsl_ran_logistic_pdf ((x - a) / b, 1) / b;
+}
+
+static double
+matrix_eval_RV_LOGISTIC (double a, double b)
+{
+ return a + b * gsl_ran_logistic (get_rng (), 1);
+}
+
+static double
+matrix_eval_CDF_LNORMAL (double x, double m, double s)
+{
+ return gsl_cdf_lognormal_P (x, log (m), s);
+}
+
+static double
+matrix_eval_IDF_LNORMAL (double P, double m, double s)
+{
+ return gsl_cdf_lognormal_Pinv (P, log (m), s);;
+}
+
+static double
+matrix_eval_PDF_LNORMAL (double x, double m, double s)
+{
+ return gsl_ran_lognormal_pdf (x, log (m), s);
+}
+
+static double
+matrix_eval_RV_LNORMAL (double m, double s)
+{
+ return gsl_ran_lognormal (get_rng (), log (m), s);
+}
+
+static double
+matrix_eval_CDF_NORMAL (double x, double u, double s)
+{
+ return gsl_cdf_gaussian_P (x - u, s);
+}
+
+static double
+matrix_eval_IDF_NORMAL (double P, double u, double s)
+{
+ return u + gsl_cdf_gaussian_Pinv (P, s);
+}
+
+static double
+matrix_eval_PDF_NORMAL (double x, double u, double s)
+{
+ return gsl_ran_gaussian_pdf ((x - u) / s, 1) / s;
+}
+
+static double
+matrix_eval_RV_NORMAL (double u, double s)
+{
+ return u + gsl_ran_gaussian (get_rng (), s);
+}
+
+static double
+matrix_eval_CDFNORM (double x)
+{
+ return gsl_cdf_ugaussian_P (x);
+}
+
+static double
+matrix_eval_PROBIT (double P)
+{
+ return gsl_cdf_ugaussian_Pinv (P);
+}
+
+static double
+matrix_eval_NORMAL (double s)
+{
+ return gsl_ran_gaussian (get_rng (), s);
+}
+
+static double
+matrix_eval_PDF_NTAIL (double x, double a, double sigma)
+{
+ return gsl_ran_gaussian_tail_pdf (x, a, sigma);;
+}
+
+static double
+matrix_eval_RV_NTAIL (double a, double sigma)
+{
+ return gsl_ran_gaussian_tail (get_rng (), a, sigma);
+}
+
+static double
+matrix_eval_CDF_PARETO (double x, double a, double b)
+{
+ return gsl_cdf_pareto_P (x, b, a);
+}
+
+static double
+matrix_eval_IDF_PARETO (double P, double a, double b)
+{
+ return gsl_cdf_pareto_Pinv (P, b, a);
+}
+
+static double
+matrix_eval_PDF_PARETO (double x, double a, double b)
+{
+ return gsl_ran_pareto_pdf (x, b, a);
+}
+
+static double
+matrix_eval_RV_PARETO (double a, double b)
+{
+ return gsl_ran_pareto (get_rng (), b, a);
+}
+
+static double
+matrix_eval_CDF_RAYLEIGH (double x, double sigma)
+{
+ return gsl_cdf_rayleigh_P (x, sigma);
+}
+
+static double
+matrix_eval_IDF_RAYLEIGH (double P, double sigma)
+{
+ return gsl_cdf_rayleigh_Pinv (P, sigma);
+}
+
+static double
+matrix_eval_PDF_RAYLEIGH (double x, double sigma)
+{
+ return gsl_ran_rayleigh_pdf (x, sigma);
+}
+
+static double
+matrix_eval_RV_RAYLEIGH (double sigma)
+{
+ return gsl_ran_rayleigh (get_rng (), sigma);
+}
+
+static double
+matrix_eval_PDF_RTAIL (double x, double a, double sigma)
+{
+ return gsl_ran_rayleigh_tail_pdf (x, a, sigma);
+}
+
+static double
+matrix_eval_RV_RTAIL (double a, double sigma)
+{
+ return gsl_ran_rayleigh_tail (get_rng (), a, sigma);
+}
+
+static double
+matrix_eval_CDF_T (double x, double df)
+{
+ return gsl_cdf_tdist_P (x, df);
+}
+
+static double
+matrix_eval_TCDF (double x, double df)
+{
+ return matrix_eval_CDF_T (x, df);
+}
+
+static double
+matrix_eval_IDF_T (double P, double df)
+{
+ return gsl_cdf_tdist_Pinv (P, df);
+}
+
+static double
+matrix_eval_PDF_T (double x, double df)
+{
+ return gsl_ran_tdist_pdf (x, df);
+}
+
+static double
+matrix_eval_RV_T (double df)
+{
+ return gsl_ran_tdist (get_rng (), df);
+}
+
+static double
+matrix_eval_CDF_T1G (double x, double a, double b)
+{
+ return gsl_cdf_gumbel1_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_T1G (double P, double a, double b)
+{
+ return gsl_cdf_gumbel1_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_T1G (double x, double a, double b)
+{
+ return gsl_ran_gumbel1_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_T1G (double a, double b)
+{
+ return gsl_ran_gumbel1 (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_T2G (double x, double a, double b)
+{
+ return gsl_cdf_gumbel1_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_T2G (double P, double a, double b)
+{
+ return gsl_cdf_gumbel1_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_T2G (double x, double a, double b)
+{
+ return gsl_ran_gumbel1_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_T2G (double a, double b)
+{
+ return gsl_ran_gumbel1 (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_UNIFORM (double x, double a, double b)
+{
+ return gsl_cdf_flat_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_UNIFORM (double P, double a, double b)
+{
+ return gsl_cdf_flat_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_UNIFORM (double x, double a, double b)
+{
+ return gsl_ran_flat_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_UNIFORM (double a, double b)
+{
+ return gsl_ran_flat (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_WEIBULL (double x, double a, double b)
+{
+ return gsl_cdf_weibull_P (x, a, b);
+}
+
+static double
+matrix_eval_IDF_WEIBULL (double P, double a, double b)
+{
+ return gsl_cdf_weibull_Pinv (P, a, b);
+}
+
+static double
+matrix_eval_PDF_WEIBULL (double x, double a, double b)
+{
+ return gsl_ran_weibull_pdf (x, a, b);
+}
+
+static double
+matrix_eval_RV_WEIBULL (double a, double b)
+{
+ return gsl_ran_weibull (get_rng (), a, b);
+}
+
+static double
+matrix_eval_CDF_BERNOULLI (double k, double p)
+{
+ return k ? 1 : 1 - p;
+}
+
+static double
+matrix_eval_PDF_BERNOULLI (double k, double p)
+{
+ return gsl_ran_bernoulli_pdf (k, p);
+}
+
+static double
+matrix_eval_RV_BERNOULLI (double p)
+{
+ return gsl_ran_bernoulli (get_rng (), p);
+}
+
+static double
+matrix_eval_CDF_BINOM (double k, double n, double p)
+{
+ return gsl_cdf_binomial_P (k, p, n);
+}
+
+static double
+matrix_eval_PDF_BINOM (double k, double n, double p)
+{
+ return gsl_ran_binomial_pdf (k, p, n);
+}
+
+static double
+matrix_eval_RV_BINOM (double n, double p)
+{
+ return gsl_ran_binomial (get_rng (), p, n);
+}
+
+static double
+matrix_eval_CDF_GEOM (double k, double p)
+{
+ return gsl_cdf_geometric_P (k, p);
+}
+
+static double
+matrix_eval_PDF_GEOM (double k, double p)
+{
+ return gsl_ran_geometric_pdf (k, p);
+}
+
+static double
+matrix_eval_RV_GEOM (double p)
+{
+ return gsl_ran_geometric (get_rng (), p);
+}
+
+static double
+matrix_eval_CDF_HYPER (double k, double a, double b, double c)
+{
+ return gsl_cdf_hypergeometric_P (k, c, a - c, b);
+}
+
+static double
+matrix_eval_PDF_HYPER (double k, double a, double b, double c)
+{
+ return gsl_ran_hypergeometric_pdf (k, c, a - c, b);
+}
+
+static double
+matrix_eval_RV_HYPER (double a, double b, double c)
+{
+ return gsl_ran_hypergeometric (get_rng (), c, a - c, b);
+}
+
+static double
+matrix_eval_PDF_LOG (double k, double p)
+{
+ return gsl_ran_logarithmic_pdf (k, p);
+}
+
+static double
+matrix_eval_RV_LOG (double p)
+{
+ return gsl_ran_logarithmic (get_rng (), p);
+}
+
+static double
+matrix_eval_CDF_NEGBIN (double k, double n, double p)
+{
+ return gsl_cdf_negative_binomial_P (k, p, n);
+}
+
+static double
+matrix_eval_PDF_NEGBIN (double k, double n, double p)
+{
+ return gsl_ran_negative_binomial_pdf (k, p, n);
+}
+
+static double
+matrix_eval_RV_NEGBIN (double n, double p)
+{
+ return gsl_ran_negative_binomial (get_rng (), p, n);
+}
+
+static double
+matrix_eval_CDF_POISSON (double k, double mu)
+{
+ return gsl_cdf_poisson_P (k, mu);
+}
+
+static double
+matrix_eval_PDF_POISSON (double k, double mu)
+{
+ return gsl_ran_poisson_pdf (k, mu);
+}
+
+static double
+matrix_eval_RV_POISSON (double mu)
+{
+ return gsl_ran_poisson (get_rng (), mu);
+}
+
struct matrix_function
{
const char *name;
{
char *tail;
long retval = strtol (*constraintsp, &tail, 10);
+ assert (tail > *constraintsp);
*constraintsp = tail;
return retval;
}
+static void
+argument_invalid (const struct matrix_function_properties *props,
+ size_t arg_index, gsl_matrix *a, size_t y, size_t x,
+ struct string *out)
+{
+ if (is_scalar (a))
+ ds_put_format (out, _("Argument %zu to matrix function %s "
+ "has invalid value %g."),
+ arg_index, props->name, gsl_matrix_get (a, y, x));
+ else
+ ds_put_format (out, _("Row %zu, column %zu of argument %zu "
+ "to matrix function %s has "
+ "invalid value %g."),
+ y + 1, x + 1, arg_index, props->name,
+ gsl_matrix_get (a, y, x));
+}
+
+static void
+argument_inequality_error (const struct matrix_function_properties *props,
+ size_t a_index, gsl_matrix *a, size_t y, size_t x,
+ size_t b_index, double b,
+ bool greater, bool equal)
+{
+ struct string s = DS_EMPTY_INITIALIZER;
+
+ argument_invalid (props, a_index, a, y, x, &s);
+ ds_put_cstr (&s, " ");
+ if (greater && equal)
+ ds_put_format (&s, _("This argument must be greater than or "
+ "equal to argument %zu, but that has value %g."),
+ b_index, b);
+ else if (greater && !equal)
+ ds_put_format (&s, _("This argument must be greater than argument %zu, "
+ "but that has value %g."),
+ b_index, b);
+ else if (equal)
+ ds_put_format (&s, _("This argument must be less than or "
+ "equal to argument %zu, but that has value %g."),
+ b_index, b);
+ else
+ ds_put_format (&s, _("This argument must be less than argument %zu, "
+ "but that has value %g."),
+ b_index, b);
+ msg (ME, ds_cstr (&s));
+ ds_destroy (&s);
+}
+
+static void
+argument_value_error (const struct matrix_function_properties *props,
+ size_t a_index, gsl_matrix *a, size_t y, size_t x,
+ double b,
+ bool greater, bool equal)
+{
+ struct string s = DS_EMPTY_INITIALIZER;
+ argument_invalid (props, a_index, a, y, x, &s);
+ ds_put_cstr (&s, " ");
+ if (greater && equal)
+ ds_put_format (&s, _("This argument must be greater than or equal to %g."),
+ b);
+ else if (greater && !equal)
+ ds_put_format (&s, _("This argument must be greater than %g."), b);
+ else if (equal)
+ ds_put_format (&s, _("This argument must be less than or equal to %g."), b);
+ else
+ ds_put_format (&s, _("This argument must be less than %g."), b);
+ msg (ME, ds_cstr (&s));
+ ds_destroy (&s);
+}
+
static bool
check_constraints (const struct matrix_function_properties *props,
gsl_matrix *args[], size_t n_args)
}
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;
- }
+ if (*constraints >= 'a' && *constraints <= 'd')
+ {
+ size_t a_index = arg_index;
+ size_t b_index = *constraints - 'a';
+ assert (a_index < n_args);
+ assert (b_index < n_args);
+
+ /* We only support one of the two arguments being non-scalar.
+ It's easier to support only the first one being non-scalar, so
+ flip things around if it's the other way. */
+ if (!is_scalar (args[b_index]))
+ {
+ assert (is_scalar (args[a_index]));
+ size_t tmp_index = a_index;
+ a_index = b_index;
+ b_index = tmp_index;
+
+ greater = !greater;
+ }
+
+ double b = to_scalar (args[b_index]);
+ MATRIX_FOR_ALL_ELEMENTS (a, y, x, args[a_index])
+ if (greater
+ ? (equal ? !(*a >= b) : !(*a > b))
+ : (equal ? !(*a <= b) : !(*a < b)))
+ {
+ argument_inequality_error (props,
+ a_index + 1, args[a_index], y, x,
+ b_index + 1, b,
+ greater, equal);
+ return false;
+ }
+ }
+ else
+ {
+ int comparand = parse_constraint_value (&constraints);
+
+ MATRIX_FOR_ALL_ELEMENTS (d, y, x, args[arg_index])
+ if (greater
+ ? (equal ? !(*d >= comparand) : !(*d > comparand))
+ : (equal ? !(*d <= comparand) : !(*d < comparand)))
+ {
+ argument_value_error (props,
+ arg_index + 1, args[arg_index], y, x,
+ comparand,
+ greater, equal);
+ return false;
+ }
+ }
}
else
{
assert (*constraints == ' ');
constraints++;
+ arg_index = SIZE_MAX;
}
}
return true;
}
+static gsl_matrix *
+matrix_expr_evaluate_d_none (
+ const struct matrix_function_properties *props UNUSED,
+ gsl_matrix *subs[] UNUSED, size_t n_subs,
+ matrix_proto_d_none *f)
+{
+ assert (n_subs == 0);
+
+ gsl_matrix *m = gsl_matrix_alloc (1, 1);
+ gsl_matrix_set (m, 0, 0, f ());
+ return m;
+}
+
static gsl_matrix *
matrix_expr_evaluate_d_d (const struct matrix_function_properties *props,
- gsl_matrix *subs[], size_t n_subs,
- matrix_proto_d_d *f)
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_d_d *f)
{
assert (n_subs == 1);
return m;
}
+static gsl_matrix *
+matrix_expr_evaluate_d_ddd (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_d_ddd *f)
+{
+ assert (n_subs == 3);
+
+ double d[3];
+ 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], d[2]));
+ return m;
+}
+
static gsl_matrix *
matrix_expr_evaluate_m_d (const struct matrix_function_properties *props,
gsl_matrix *subs[], size_t n_subs,
}
static gsl_matrix *
-matrix_expr_evaluate_m_m_e (const struct matrix_function_properties *props,
+matrix_expr_evaluate_m_e (const struct matrix_function_properties *props,
gsl_matrix *subs[], size_t n_subs,
- matrix_proto_m_m_e *f)
+ matrix_proto_m_e *f)
{
assert (n_subs == 1);
}
static gsl_matrix *
-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)
+matrix_expr_evaluate_m_ed (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_ed *f)
{
assert (n_subs == 2);
if (!check_scalar_arg (props->name, subs, 1))
}
static gsl_matrix *
-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)
+matrix_expr_evaluate_m_edd (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_edd *f)
{
assert (n_subs == 3);
if (!check_scalar_arg (props->name, subs, 1) || !check_scalar_arg (props->name, subs, 2))
}
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)
+matrix_expr_evaluate_m_eddd (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_eddd *f)
{
assert (n_subs == 4);
for (size_t i = 1; i < 4; i++)
return subs[0];
}
+static gsl_matrix *
+matrix_expr_evaluate_m_eed (const struct matrix_function_properties *props,
+ gsl_matrix *subs[], size_t n_subs,
+ matrix_proto_m_eed *f)
+{
+ assert (n_subs == 3);
+ if (!check_scalar_arg (props->name, subs, 2))
+ return NULL;
+
+ if (!is_scalar (subs[0]) && !is_scalar (subs[1])
+ && (subs[0]->size1 != subs[1]->size1 || subs[0]->size2 != subs[1]->size2))
+ {
+ msg (ME, _("Arguments 1 and 2 to %s have dimensions %zu×%zu and "
+ "%zu×%zu, but %s requires these arguments either to have "
+ "the same dimensions or for one of them to be a scalar."),
+ props->name,
+ subs[0]->size1, subs[0]->size2,
+ subs[1]->size1, subs[1]->size2,
+ props->name);
+ return NULL;
+ }
+
+ if (!check_constraints (props, subs, n_subs))
+ return NULL;
+
+ double c = to_scalar (subs[2]);
+
+ if (is_scalar (subs[0]))
+ {
+ double a = to_scalar (subs[0]);
+ MATRIX_FOR_ALL_ELEMENTS (b, y, x, subs[1])
+ *b = f (a, *b, c);
+ return subs[1];
+ }
+ else
+ {
+ double b = to_scalar (subs[1]);
+ MATRIX_FOR_ALL_ELEMENTS (a, y, x, subs[0])
+ *a = f (*a, b, c);
+ 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,
projects / pspp / commitdiff