diff --git a/ext/gsl_native/extconf.rb b/ext/gsl_native/extconf.rb index a966a45a..4533b29c 100644 --- a/ext/gsl_native/extconf.rb +++ b/ext/gsl_native/extconf.rb @@ -115,4 +115,12 @@ def gsl_gem_config(target, dir = 'ext') gsl_have_header('tamuanova', 'tamu_anova/tamu_anova.h') end +have_struct_member('gsl_multifit_fdfsolver', 'J', 'gsl/gsl_multifit_nlin.h') +have_func('gsl_sf_mathieu_a_e', 'gsl/gsl_sf_mathieu.h'); +have_func('gsl_sf_mathieu_b_e', 'gsl/gsl_sf_mathieu.h'); +have_func('gsl_sf_mathieu_ce_e', 'gsl/gsl_sf_mathieu.h'); +have_func('gsl_sf_mathieu_se_e', 'gsl/gsl_sf_mathieu.h'); +have_func('gsl_sf_mathieu_Mc_e', 'gsl/gsl_sf_mathieu.h'); +have_func('gsl_sf_mathieu_Ms_e', 'gsl/gsl_sf_mathieu.h'); + create_makefile('gsl_native') diff --git a/ext/gsl_native/histogram.c b/ext/gsl_native/histogram.c index 9324bf8e..84fb2ae5 100644 --- a/ext/gsl_native/histogram.c +++ b/ext/gsl_native/histogram.c @@ -1169,6 +1169,9 @@ static VALUE rb_gsl_histogram_fit_gaussian(int argc, VALUE *argv, VALUE obj) size_t n, dof; /* # of data points */ size_t p = 3; /* # of fitting parameters */ gsl_multifit_function_fdf f; +#ifndef GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix *J = NULL; +#endif gsl_matrix *covar = NULL; gsl_vector *x = NULL; double sigma, mean, height, errs, errm, errh, chi2; @@ -1197,6 +1200,9 @@ static VALUE rb_gsl_histogram_fit_gaussian(int argc, VALUE *argv, VALUE obj) hh.binend = binend; n = binend - binstart + 1; +#ifndef GSL_MULTIFIT_FDFSOLVER_J + J = gsl_matrix_alloc(n, p); +#endif covar = gsl_matrix_alloc(p, p); f.f = Gaussian_f; @@ -1219,7 +1225,12 @@ static VALUE rb_gsl_histogram_fit_gaussian(int argc, VALUE *argv, VALUE obj) sigma = sqrt(gsl_vector_get(s->x, 0)); mean = gsl_vector_get(s->x, 1); height = gsl_vector_get(s->x, 2)*sigma*sqrt(2*M_PI); +#ifdef GSL_MULTIFIT_FDFSOLVER_J gsl_multifit_covar(s->J, 0.0, covar); +#else + gsl_multifit_fdfsolver_jac(s, J); + gsl_multifit_covar(J, 0.0, covar); +#endif chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f)); /* not reduced chi-square */ dof = n - p; errs = sqrt(chi2/dof*gsl_matrix_get(covar, 0, 0))/sigma/2; @@ -1228,6 +1239,9 @@ static VALUE rb_gsl_histogram_fit_gaussian(int argc, VALUE *argv, VALUE obj) gsl_multifit_fdfsolver_free(s); gsl_vector_free(x); +#ifndef GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix_free(J); +#endif gsl_matrix_free(covar); return rb_ary_new3(8, rb_float_new(sigma), rb_float_new(mean), rb_float_new(height), rb_float_new(errs), @@ -1305,6 +1319,9 @@ static VALUE rb_gsl_histogram_fit_rayleigh(int argc, VALUE *argv, VALUE obj) size_t n, dof; /* # of data points */ size_t p = 2; /* # of fitting parameters */ gsl_multifit_function_fdf f; +#ifndef GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix *J = NULL; +#endif gsl_matrix *covar = NULL; gsl_vector *x = NULL; double sigma, height, errs, errh, chi2; @@ -1332,6 +1349,9 @@ static VALUE rb_gsl_histogram_fit_rayleigh(int argc, VALUE *argv, VALUE obj) hh.binend = binend; n = binend - binstart + 1; +#ifndef GSL_MULTIFIT_FDFSOLVER_J + J = gsl_matrix_alloc(n, p); +#endif covar = gsl_matrix_alloc(p, p); f.f = Rayleigh_f; @@ -1353,7 +1373,12 @@ static VALUE rb_gsl_histogram_fit_rayleigh(int argc, VALUE *argv, VALUE obj) } while (status == GSL_CONTINUE); sigma = sqrt(gsl_vector_get(s->x, 0)); height = gsl_vector_get(s->x, 1)*sigma*sigma; +#ifdef GSL_MULTIFIT_FDFSOLVER_J gsl_multifit_covar(s->J, 0.0, covar); +#else + gsl_multifit_fdfsolver_jac(s, J); + gsl_multifit_covar(J, 0.0, covar); +#endif chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f)); /* not reduced chi-square */ dof = n - p; errs = sqrt(chi2/dof*gsl_matrix_get(covar, 0, 0))/sigma/2; @@ -1361,6 +1386,9 @@ static VALUE rb_gsl_histogram_fit_rayleigh(int argc, VALUE *argv, VALUE obj) gsl_multifit_fdfsolver_free(s); gsl_vector_free(x); +#ifndef GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix_free(J); +#endif gsl_matrix_free(covar); return rb_ary_new3(6, rb_float_new(sigma), rb_float_new(height), rb_float_new(errs), @@ -1441,6 +1469,9 @@ static VALUE rb_gsl_histogram_fit_xexponential(int argc, VALUE *argv, VALUE obj) size_t n, dof; /* # of data points */ size_t p = 2; /* # of fitting parameters */ gsl_multifit_function_fdf f; +#ifndef GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix *J = NULL; +#endif gsl_matrix *covar = NULL; gsl_vector *x = NULL; double b, height, errs, errh, chi2; @@ -1468,6 +1499,9 @@ static VALUE rb_gsl_histogram_fit_xexponential(int argc, VALUE *argv, VALUE obj) hh.binend = binend; n = binend - binstart + 1; +#ifndef GSL_MULTIFIT_FDFSOLVER_J + J = gsl_matrix_alloc(n, p); +#endif covar = gsl_matrix_alloc(p, p); f.f = xExponential_f; @@ -1489,7 +1523,12 @@ static VALUE rb_gsl_histogram_fit_xexponential(int argc, VALUE *argv, VALUE obj) } while (status == GSL_CONTINUE); b = gsl_vector_get(s->x, 0); height = gsl_vector_get(s->x, 1); +#ifdef GSL_MULTIFIT_FDFSOLVER_J gsl_multifit_covar(s->J, 0.0, covar); +#else + gsl_multifit_fdfsolver_jac(s, J); + gsl_multifit_covar(J, 0.0, covar); +#endif chi2 = gsl_pow_2(gsl_blas_dnrm2(s->f)); /* not reduced chi-square */ dof = n - p; errs = sqrt(chi2/dof*gsl_matrix_get(covar, 0, 0)); @@ -1497,6 +1536,9 @@ static VALUE rb_gsl_histogram_fit_xexponential(int argc, VALUE *argv, VALUE obj) gsl_multifit_fdfsolver_free(s); gsl_vector_free(x); +#ifndef GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix_free(J); +#endif gsl_matrix_free(covar); return rb_ary_new3(6, rb_float_new(b), rb_float_new(height), rb_float_new(errs), diff --git a/ext/gsl_native/multifit.c b/ext/gsl_native/multifit.c index 55e4f026..251035d7 100644 --- a/ext/gsl_native/multifit.c +++ b/ext/gsl_native/multifit.c @@ -324,6 +324,9 @@ static VALUE rb_gsl_multifit_fdfsolver_test_gradient(int argc, VALUE *argv, VALU { gsl_multifit_fdfsolver *solver = NULL; gsl_vector *g = NULL; +#ifndef HAVE_GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix *J = NULL; +#endif int status; double epsabs; Data_Get_Struct(obj, gsl_multifit_fdfsolver, solver); @@ -331,7 +334,14 @@ static VALUE rb_gsl_multifit_fdfsolver_test_gradient(int argc, VALUE *argv, VALU case 1: Need_Float(argv[0]); g = gsl_vector_alloc(solver->x->size); +#ifdef HAVE_GSL_MULTIFIT_FDFSOLVER_J gsl_multifit_gradient(solver->J, solver->f, g); +#else + J = gsl_matrix_alloc(solver->fdf->n, solver->fdf->p); + gsl_multifit_fdfsolver_jac(solver, J); + gsl_multifit_gradient(J, solver->f, g); + gsl_matrix_free(J); +#endif epsabs = NUM2DBL(argv[0]); status = gsl_multifit_test_gradient(g, epsabs); gsl_vector_free(g); @@ -352,16 +362,37 @@ static VALUE rb_gsl_multifit_fdfsolver_test_gradient(int argc, VALUE *argv, VALU static VALUE rb_gsl_multifit_fdfsolver_gradient(int argc, VALUE *argv, VALUE obj) { gsl_multifit_fdfsolver *solver = NULL; +#ifndef HAVE_GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix *J = NULL; +#endif gsl_vector *g = NULL; // local variable "status" declared and set, but never used //int status; Data_Get_Struct(obj, gsl_multifit_fdfsolver, solver); +#ifndef HAVE_GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix_alloc(solver->fdf->n, solver->fdf->p); + gsl_multifit_fdfsolver_jac(solver, J); +#endif if (argc == 1) { +#ifndef HAVE_GSL_MULTIFIT_FDFSOLVER_J + int retval = 0; +#endif Data_Get_Vector(argv[0], g); +#ifdef HAVE_GSL_MULTIFIT_FDFSOLVER_J return INT2FIX(gsl_multifit_gradient(solver->J, solver->f, g)); +#else + retval = gsl_multifit_gradient(J, solver->f, g); + gsl_matrix_free(J); + return INT2FIX(retval); +#endif } else { g = gsl_vector_alloc(solver->x->size); +#ifdef HAVE_GSL_MULTIFIT_FDFSOLVER_J /*status =*/ gsl_multifit_gradient(solver->J, solver->f, g); +#else + /*status =*/ gsl_multifit_gradient(J, solver->f, g); + gsl_matrix_free(J); +#endif return Data_Wrap_Struct(cgsl_vector, 0, gsl_vector_free, g); } } @@ -369,6 +400,9 @@ static VALUE rb_gsl_multifit_fdfsolver_gradient(int argc, VALUE *argv, VALUE obj static VALUE rb_gsl_multifit_fdfsolver_covar(int argc, VALUE *argv, VALUE obj) { gsl_multifit_fdfsolver *solver = NULL; +#ifndef HAVE_GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix *J = NULL; +#endif gsl_matrix *covar = NULL; double epsrel; // local variable "status" declared and set, but never used @@ -378,14 +412,32 @@ static VALUE rb_gsl_multifit_fdfsolver_covar(int argc, VALUE *argv, VALUE obj) Data_Get_Struct(obj, gsl_multifit_fdfsolver, solver); epsrel = NUM2DBL(argv[0]); switch (argc) { +#ifndef HAVE_GSL_MULTIFIT_FDFSOLVER_J + int retval = 0; +#endif case 1: covar = gsl_matrix_alloc(solver->x->size, solver->x->size); +#ifdef HAVE_GSL_MULTIFIT_FDFSOLVER_J /*status =*/ gsl_multifit_covar(solver->J, epsrel, covar); +#else + J = gsl_matrix_alloc(solver->fdf->n, solver->fdf->p); + gsl_multifit_fdfsolver_jac(solver, J); + /*status =*/ gsl_multifit_covar(J, epsrel, covar); + gsl_matrix_free(J); +#endif return Data_Wrap_Struct(cgsl_matrix, 0, gsl_matrix_free, covar); break; case 2: Data_Get_Matrix(argv[1], covar); +#ifdef HAVE_GSL_MULTIFIT_FDFSOLVER_J return INT2FIX(gsl_multifit_covar(solver->J, epsrel, covar)); +#else + J = gsl_matrix_alloc(solver->fdf->n, solver->fdf->p); + gsl_multifit_fdfsolver_jac(solver, J); + retval = gsl_multifit_covar(J, epsrel, covar); + gsl_matrix_free(J); + return INT2FIX(retval); +#endif break; default: rb_raise(rb_eArgError, "wrong number of arguments"); @@ -417,8 +469,17 @@ static VALUE rb_gsl_multifit_fdfsolver_f(VALUE obj) static VALUE rb_gsl_multifit_fdfsolver_J(VALUE obj) { gsl_multifit_fdfsolver *solver = NULL; +#ifndef HAVE_GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix *J = NULL; +#endif Data_Get_Struct(obj, gsl_multifit_fdfsolver, solver); +#ifdef HAVE_GSL_MULTIFIT_FDFSOLVER_J return Data_Wrap_Struct(cgsl_matrix_view_ro, 0, NULL, solver->J); +#else + J = gsl_matrix_alloc(solver->fdf->n, solver->fdf->p); + gsl_multifit_fdfsolver_jac(solver, J); + return Data_Wrap_Struct(cgsl_matrix_view_ro, 0, gsl_matrix_free, J); +#endif } /* singleton */ @@ -1626,6 +1687,9 @@ static VALUE rb_gsl_multifit_fit(int argc, VALUE *argv, VALUE klass) size_t n, dof; /* # of data points */ size_t p; /* # of fitting parameters */ gsl_multifit_function_fdf f; +#ifndef HAVE_GSL_MULTIFIT_FDFSOLVER_J + gsl_matrix *J = NULL; +#endif gsl_matrix *covar = NULL; gsl_vector *v = NULL; gsl_vector *x, *y, *w = NULL; @@ -1699,7 +1763,14 @@ static VALUE rb_gsl_multifit_fit(int argc, VALUE *argv, VALUE klass) covar = gsl_matrix_alloc(p, p); chi2 = gsl_pow_2(gsl_blas_dnrm2(solver->f)); /* not reduced chi-square */ dof = n - p; +#ifdef HAVE_GSL_MULTIFIT_FDFSOLVER_J gsl_multifit_covar(solver->J, 0.0, covar); +#else + J = gsl_matrix_alloc(n, p); + gsl_multifit_fdfsolver_jac(solver, J); + gsl_multifit_covar(J, 0.0, covar); + gsl_matrix_free(J); +#endif for (i = 0; i < p; i++) gsl_vector_set(verr, i, sqrt(chi2/dof*gsl_matrix_get(covar, i, i))); gsl_matrix_free(covar); diff --git a/ext/gsl_native/sf_ellint.c b/ext/gsl_native/sf_ellint.c index fd2e8baf..66169d21 100644 --- a/ext/gsl_native/sf_ellint.c +++ b/ext/gsl_native/sf_ellint.c @@ -82,17 +82,25 @@ static VALUE rb_gsl_sf_ellint_P_e(VALUE obj, VALUE phi, VALUE k, static VALUE rb_gsl_sf_ellint_D(int argc, VALUE *argv, VALUE obj) { if (argc == 3) +#if GSL_MAJOR_VERSION > 1 + return rb_gsl_sf_eval_double2_m(gsl_sf_ellint_D, argv[0], argv[1], +#else return rb_gsl_sf_eval_double3_m(gsl_sf_ellint_D, argv[0], argv[1], argv[2], +#endif INT2FIX(GSL_PREC_DOUBLE)); else +#if GSL_MAJOR_VERSION > 1 + return rb_gsl_sf_eval_double2_m(gsl_sf_ellint_D, argv[0], argv[1], +#else return rb_gsl_sf_eval_double3_m(gsl_sf_ellint_D, argv[0], argv[1], argv[2], +#endif argv[3]); } static VALUE rb_gsl_sf_ellint_D_e(VALUE obj, VALUE phi, VALUE k, - VALUE n, VALUE m) + VALUE m) { - return rb_gsl_sf_eval_e_double3_m(gsl_sf_ellint_D_e, phi, k, n, m); + return rb_gsl_sf_eval_e_double2_m(gsl_sf_ellint_D_e, phi, k, m); } static VALUE rb_gsl_sf_ellint_RC(int argc, VALUE *argv, VALUE obj) diff --git a/ext/gsl_native/sf_mathieu.c b/ext/gsl_native/sf_mathieu.c index 8af91b7a..b8d33e68 100644 --- a/ext/gsl_native/sf_mathieu.c +++ b/ext/gsl_native/sf_mathieu.c @@ -133,11 +133,19 @@ static VALUE sf_mathieu_eval_e_int2_double2(VALUE n1, VALUE n2, VALUE qq, VALUE /**********/ static VALUE rb_gsl_sf_mathieu_a_e(VALUE module, VALUE order, VALUE qq) { +#ifdef HAVE_GSL_SF_MATHIEU_A_E + return rb_gsl_sf_eval_e_int_double(gsl_sf_mathieu_a_e, order, qq); +#else return rb_gsl_sf_eval_e_int_double(gsl_sf_mathieu_a, order, qq); +#endif } static VALUE rb_gsl_sf_mathieu_a(VALUE module, VALUE order, VALUE qq) { +#ifdef HAVE_GSL_SF_MATHIEU_A_E + return sf_mathieu_eval(order, qq, gsl_sf_mathieu_a_e); +#else return sf_mathieu_eval(order, qq, gsl_sf_mathieu_a); +#endif } static VALUE rb_gsl_sf_mathieu_a_array(VALUE module, int argc, VALUE *argv) { @@ -145,11 +153,19 @@ static VALUE rb_gsl_sf_mathieu_a_array(VALUE module, int argc, VALUE *argv) } static VALUE rb_gsl_sf_mathieu_b_e(VALUE module, VALUE order, VALUE qq) { +#ifdef HAVE_GSL_SF_MATHIEU_B_E + return rb_gsl_sf_eval_e_int_double(gsl_sf_mathieu_b_e, order, qq); +#else return rb_gsl_sf_eval_e_int_double(gsl_sf_mathieu_b, order, qq); +#endif } static VALUE rb_gsl_sf_mathieu_b(VALUE module, VALUE order, VALUE qq) { +#ifdef HAVE_GSL_SF_MATHIEU_B_E + return sf_mathieu_eval(order, qq, gsl_sf_mathieu_b_e); +#else return sf_mathieu_eval(order, qq, gsl_sf_mathieu_b); +#endif } static VALUE rb_gsl_sf_mathieu_b_array(VALUE module, int argc, VALUE *argv) { @@ -157,11 +173,19 @@ static VALUE rb_gsl_sf_mathieu_b_array(VALUE module, int argc, VALUE *argv) } static VALUE rb_gsl_sf_mathieu_ce_e(VALUE module, VALUE order, VALUE qq, VALUE zz) { +#ifdef HAVE_GSL_SF_MATHIEU_CE_E + return sf_mathieu_eval_e_int_double2(order, qq, zz, gsl_sf_mathieu_ce_e); +#else return sf_mathieu_eval_e_int_double2(order, qq, zz, gsl_sf_mathieu_ce); +#endif } static VALUE rb_gsl_sf_mathieu_ce(VALUE module, VALUE order, VALUE qq, VALUE zz) { +#ifdef HAVE_GSL_SF_MATHIEU_CE_E + return sf_mathieu_eval_int_double2(order, qq, zz, gsl_sf_mathieu_ce_e); +#else return sf_mathieu_eval_int_double2(order, qq, zz, gsl_sf_mathieu_ce); +#endif } static VALUE rb_gsl_sf_mathieu_ce_array(VALUE module, int argc, VALUE *argv) { @@ -169,11 +193,19 @@ static VALUE rb_gsl_sf_mathieu_ce_array(VALUE module, int argc, VALUE *argv) } static VALUE rb_gsl_sf_mathieu_se_e(VALUE module, VALUE order, VALUE qq, VALUE zz) { +#ifdef HAVE_GSL_SF_MATHIEU_SE_E + return sf_mathieu_eval_e_int_double2(order, qq, zz, gsl_sf_mathieu_se_e); +#else return sf_mathieu_eval_e_int_double2(order, qq, zz, gsl_sf_mathieu_se); +#endif } static VALUE rb_gsl_sf_mathieu_se(VALUE module, VALUE order, VALUE qq, VALUE zz) { +#ifdef HAVE_GSL_SF_MATHIEU_SE_E + return sf_mathieu_eval_int_double2(order, qq, zz, gsl_sf_mathieu_se_e); +#else return sf_mathieu_eval_int_double2(order, qq, zz, gsl_sf_mathieu_se); +#endif } static VALUE rb_gsl_sf_mathieu_se_array(VALUE module, int argc, VALUE *argv) { @@ -183,11 +215,19 @@ static VALUE rb_gsl_sf_mathieu_se_array(VALUE module, int argc, VALUE *argv) /*****/ static VALUE rb_gsl_sf_mathieu_Mc_e(VALUE module, VALUE n1, VALUE n2, VALUE q, VALUE x) { +#ifdef HAVE_GSL_SF_MATHIEU_MC_E + return sf_mathieu_eval_e_int2_double2(n1, n2, q, x, gsl_sf_mathieu_Mc_e); +#else return sf_mathieu_eval_e_int2_double2(n1, n2, q, x, gsl_sf_mathieu_Mc); +#endif } static VALUE rb_gsl_sf_mathieu_Mc(VALUE module, VALUE n1, VALUE n2, VALUE q, VALUE x) { +#ifdef HAVE_GSL_SF_MATHIEU_MC_E + return sf_mathieu_eval2(n1, n2, q, x, gsl_sf_mathieu_Mc_e); +#else return sf_mathieu_eval2(n1, n2, q, x, gsl_sf_mathieu_Mc); +#endif } static VALUE rb_gsl_sf_mathieu_Mc_array(VALUE module, int argc, VALUE *argv) { @@ -195,11 +235,19 @@ static VALUE rb_gsl_sf_mathieu_Mc_array(VALUE module, int argc, VALUE *argv) } static VALUE rb_gsl_sf_mathieu_Ms_e(VALUE module, VALUE n1, VALUE n2, VALUE q, VALUE x) { +#ifdef HAVE_GSL_SF_MATHIEU_MS_E + return sf_mathieu_eval_e_int2_double2(n1, n2, q, x, gsl_sf_mathieu_Ms_e); +#else return sf_mathieu_eval_e_int2_double2(n1, n2, q, x, gsl_sf_mathieu_Ms); +#endif } static VALUE rb_gsl_sf_mathieu_Ms(VALUE module, VALUE n1, VALUE n2, VALUE q, VALUE x) { +#ifdef HAVE_GSL_SF_MATHIEU_MS_E + return sf_mathieu_eval2(n1, n2, q, x, gsl_sf_mathieu_Ms_e); +#else return sf_mathieu_eval2(n1, n2, q, x, gsl_sf_mathieu_Ms); +#endif } static VALUE rb_gsl_sf_mathieu_Ms_array(VALUE module, int argc, VALUE *argv) { diff --git a/ext/gsl_native/tamu_anova.c b/ext/gsl_native/tamu_anova.c index 123ef20d..5f9c9f88 100644 --- a/ext/gsl_native/tamu_anova.c +++ b/ext/gsl_native/tamu_anova.c @@ -32,7 +32,7 @@ VALUE rb_tamu_anova_alloc(int argc, VALUE *argv, VALUE klass) VALUE rb_tamu_anova_printtable(VALUE *vTable) { struct tamu_anova_table *table; - Data_Get_Struct(vTable, struct tamu_anova_table, table); + Data_Get_Struct(*vTable, struct tamu_anova_table, table); tamu_anova_printtable(*table); return Qtrue; }