examples/histogram/xexp.rb examples/poly/eval_derivs.rb tests/monte/miser.rb tests/monte/vegas.rb @@ -1,3 +1,10 @@ +Sat Aug 8 22:36:24 BOT 2009 + * Added GSL::Monte::Miser#params_get, GSL::Monte::Vegas#params_get + * Added GSL::BSpline::greville_abscissa(i) + * Added GSL::Poly#eval_derivs(x) + * Added a new minimization algorithm "quad_golden" + * Added a new multimin algorithm "nmsimplex2rand" + Tue Jun 9 10:42:30 JST 2009 * Fixed bug in linalg.c (bug [#25842]) ChangeLog @@ -1 +1,2 @@ -1.11.2 +1.12.90 + VERSION @@ -98,19 +98,33 @@ static VALUE rb_gsl_bspline_eval(int argc, VALUE *argv, VALUE obj) return vB; } +#ifdef GSL_1_13_LATER +static VALUE rb_gsl_bspline_greville_abscissa(VALUE obj, VALUE i) +{ + gsl_bspline_workspace *w; + Data_Get_Struct(obj, gsl_bspline_workspace, w); + return rb_float_new(gsl_bspline_greville_abscissa(i, w)); +} +#endif + void Init_bspline(VALUE module) { - cBSWS = rb_define_class_under(module, "BSpline", cGSL_Object); - - rb_define_singleton_method(cBSWS, "alloc", rb_gsl_bspline_alloc, 2); - - rb_define_method(cBSWS, "ncoeffs", rb_gsl_bspline_ncoeffs, 0); - rb_define_method(cBSWS, "order", rb_gsl_bspline_order, 0); - rb_define_method(cBSWS, "nbreak", rb_gsl_bspline_nbreak, 0); - rb_define_method(cBSWS, "breakpoint", rb_gsl_bspline_breakpoint, 1); - rb_define_method(cBSWS, "knots", rb_gsl_bspline_knots, 1); - rb_define_method(cBSWS, "knots_uniform", rb_gsl_bspline_knots_uniform, -1); - rb_define_singleton_method(cBSWS, "knots_uniform", rb_gsl_bspline_knots_uniform, -1); - rb_define_method(cBSWS, "eval", rb_gsl_bspline_eval, -1); + cBSWS = rb_define_class_under(module, "BSpline", cGSL_Object); + + rb_define_singleton_method(cBSWS, "alloc", rb_gsl_bspline_alloc, 2); + + rb_define_method(cBSWS, "ncoeffs", rb_gsl_bspline_ncoeffs, 0); + rb_define_method(cBSWS, "order", rb_gsl_bspline_order, 0); + rb_define_method(cBSWS, "nbreak", rb_gsl_bspline_nbreak, 0); + rb_define_method(cBSWS, "breakpoint", rb_gsl_bspline_breakpoint, 1); + rb_define_method(cBSWS, "knots", rb_gsl_bspline_knots, 1); + rb_define_method(cBSWS, "knots_uniform", rb_gsl_bspline_knots_uniform, -1); + rb_define_singleton_method(cBSWS, "knots_uniform", rb_gsl_bspline_knots_uniform, -1); + rb_define_method(cBSWS, "eval", rb_gsl_bspline_eval, -1); + +#ifdef GSL_1_13_LATER + rb_define_method(cBSWS, "greville_abscissa", rb_gsl_bspline_greville_abscissa, 1); +#endif + } #endif ext/bspline.c @@ -119,6 +119,10 @@ def check_version(configfile) if ver >= "1.11" configfile.printf("#ifndef GSL_1_11_LATER\n#define GSL_1_11_LATER\n#endif\n") end + + if ver >= "1.12.90" + configfile.printf("#ifndef GSL_1_13_LATER\n#define GSL_1_13_LATER\n#endif\n") + end end end ext/extconf.rb @@ -18,6 +18,9 @@ double rb_gsl_function_f(double x, void *p); enum { GSL_MIN_FMINIMIZER_GOLDENSECTION, GSL_MIN_FMINIMIZER_BRENT, +#ifdef GSL_1_13_LATER + GSL_MIN_FMINIMIZER_QUAD_GOLDEN, +#endif }; static const gsl_min_fminimizer_type* rb_gsl_min_fminimizer_type_get(VALUE t); @@ -31,8 +34,12 @@ static const gsl_min_fminimizer_type* rb_gsl_min_fminimizer_type_get(VALUE t) return gsl_min_fminimizer_goldensection; else if (str_tail_grep(name, "brent") == 0) return gsl_min_fminimizer_brent; +#ifdef GSL_1_13_LATER + else if (str_tail_grep(name, "quad_golden") == 0) + return gsl_min_fminimizer_quad_golden; +#endif else - rb_raise(rb_eTypeError, "unknown type %s (goldensection or brent expected)", + rb_raise(rb_eTypeError, "unknown type %s (goldensection, brent or quad_golden expected)", name); break; case T_FIXNUM: @@ -43,8 +50,13 @@ static const gsl_min_fminimizer_type* rb_gsl_min_fminimizer_type_get(VALUE t) case GSL_MIN_FMINIMIZER_BRENT: return gsl_min_fminimizer_brent; break; +#ifdef GSL_1_13_LATER + case GSL_MIN_FMINIMIZER_QUAD_GOLDEN: + return gsl_min_fminimizer_quad_golden; + break; +#endif default: - rb_raise(rb_eTypeError, "unknown type (GOLDENSECION or BRENT expected)"); + rb_raise(rb_eTypeError, "unknown type (GOLDENSECION or BRENT or QUAD_GOLDEN expected)"); break; } break; @@ -210,6 +222,10 @@ void Init_gsl_min(VALUE module) INT2FIX(GSL_MIN_FMINIMIZER_BRENT)); rb_define_const(cgsl_fminimizer, "Brent", INT2FIX(GSL_MIN_FMINIMIZER_BRENT)); +#ifdef GSL_1_13_LATER + rb_define_const(cgsl_fminimizer, "QUAD_GOLDEN", + INT2FIX(GSL_MIN_FMINIMIZER_QUAD_GOLDEN)); +#endif rb_define_singleton_method(cgsl_fminimizer, "new", rb_gsl_min_fminimizer_new, 1); rb_define_singleton_method(cgsl_fminimizer, "alloc", rb_gsl_min_fminimizer_new, 1); ext/min.c @@ -24,7 +24,9 @@ static VALUE cgsl_monte_plain; static VALUE cgsl_monte_miser; static VALUE cgsl_monte_vegas; static VALUE cgsl_monte_function; - +#ifdef GSL_1_13_LATER +static VALUE cgsl_monte_miser_params, cgsl_monte_vegas_params; +#endif EXTERN VALUE cgsl_vector; enum { @@ -659,10 +661,193 @@ static VALUE rb_gsl_monte_vegas_state(VALUE obj) INT2FIX(s->mode), INT2FIX(s->verbose)); } +#ifdef GSL_1_13_LATER +static VALUE rb_gsl_monte_miser_params_get(VALUE obj) +{ + gsl_monte_miser_state *s = NULL; + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_state, s); + p = (gsl_monte_miser_params *) malloc(sizeof(gsl_monte_miser_params)); + gsl_monte_miser_params_get(s, p); + return Data_Wrap_Struct(cgsl_monte_miser_params, 0, free, p); +} +static VALUE rb_gsl_monte_miser_params_set(VALUE obj, VALUE params) +{ + gsl_monte_miser_state *s = NULL; + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_state, s); + Data_Get_Struct(params, gsl_monte_miser_params, p); + gsl_monte_miser_params_set(s, p); + return Qtrue; +} +static VALUE rb_gsl_monte_miser_params_get_estimate_frac(VALUE obj) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + return rb_float_new(p->estimate_frac); +} +static VALUE rb_gsl_monte_miser_params_set_estimate_frac(VALUE obj, VALUE val) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + p->estimate_frac = NUM2DBL(val); + return val; +} +static VALUE rb_gsl_monte_miser_params_get_min_calls(VALUE obj) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + return INT2FIX(p->min_calls); +} +static VALUE rb_gsl_monte_miser_params_set_min_calls(VALUE obj, VALUE val) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + p->min_calls = (size_t) FIX2INT(val); + return val; +} +static VALUE rb_gsl_monte_miser_params_get_min_calls_per_bisection(VALUE obj) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + return INT2FIX(p->min_calls_per_bisection); +} +static VALUE rb_gsl_monte_miser_params_set_min_calls_per_bisection(VALUE obj, VALUE val) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + p->min_calls_per_bisection = (size_t) FIX2INT(val); + return val; +} +static VALUE rb_gsl_monte_miser_params_get_alpha(VALUE obj) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + return rb_float_new(p->alpha); +} +static VALUE rb_gsl_monte_miser_params_set_alpha(VALUE obj, VALUE val) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + p->alpha = NUM2DBL(val); + return val; +} +static VALUE rb_gsl_monte_miser_params_get_dither(VALUE obj) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + return rb_float_new(p->dither); +} +static VALUE rb_gsl_monte_miser_params_set_dither(VALUE obj, VALUE val) +{ + gsl_monte_miser_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_miser_params, p); + p->dither = NUM2DBL(val); + return val; +} + +static VALUE rb_gsl_monte_vegas_params_get(VALUE obj) +{ + gsl_monte_vegas_state *s = NULL; + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_state, s); + p = (gsl_monte_vegas_params *) malloc(sizeof(gsl_monte_vegas_params)); + gsl_monte_vegas_params_get(s, p); + return Data_Wrap_Struct(cgsl_monte_vegas_params, 0, free, p); +} +static VALUE rb_gsl_monte_vegas_params_set(VALUE obj, VALUE params) +{ + gsl_monte_vegas_state *s = NULL; + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_state, s); + Data_Get_Struct(params, gsl_monte_vegas_params, p); + gsl_monte_vegas_params_set(s, p); + return Qtrue; +} +static VALUE rb_gsl_monte_vegas_params_get_alpha(VALUE obj) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + return rb_float_new(p->alpha); +} +static VALUE rb_gsl_monte_vegas_params_set_alpha(VALUE obj, VALUE val) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + p->alpha = NUM2DBL(val); + return val; +} +static VALUE rb_gsl_monte_vegas_params_get_iterations(VALUE obj) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + return INT2FIX(p->iterations); +} +static VALUE rb_gsl_monte_vegas_params_set_iterations(VALUE obj, VALUE val) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + p->iterations = (size_t) FIX2INT(val); + return val; +} +static VALUE rb_gsl_monte_vegas_params_get_stage(VALUE obj) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + return INT2FIX(p->stage); +} +static VALUE rb_gsl_monte_vegas_params_set_stage(VALUE obj, VALUE val) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + p->stage = FIX2INT(val); + return val; +} +static VALUE rb_gsl_monte_vegas_params_get_mode(VALUE obj) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + return INT2FIX(p->mode); +} +static VALUE rb_gsl_monte_vegas_params_set_mode(VALUE obj, VALUE val) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + p->mode = FIX2INT(val); + return val; +} +static VALUE rb_gsl_monte_vegas_params_get_verbose(VALUE obj) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + return INT2FIX(p->verbose); +} +static VALUE rb_gsl_monte_vegas_params_set_verbose(VALUE obj, VALUE val) +{ + gsl_monte_vegas_params *p = NULL; + Data_Get_Struct(obj, gsl_monte_vegas_params, p); + p->verbose = FIX2INT(val); + return val; +} +static VALUE rb_gsl_monte_vegas_runval(VALUE obj) +{ + gsl_monte_vegas_state *s = NULL; + double res, sig; + VALUE ary; + Data_Get_Struct(obj, gsl_monte_vegas_state, s); + gsl_monte_vegas_runval(s, &res, &sig); + ary = rb_ary_new2(2); + rb_ary_store(ary, 0, rb_float_new(res)); + rb_ary_store(ary, 1, rb_float_new(sig)); + return ary; +} +#endif + void Init_gsl_monte(VALUE module) { VALUE mgsl_monte; + mgsl_monte = rb_define_module_under(module, "Monte"); rb_define_const(mgsl_monte, "PLAIN", INT2FIX(GSL_MONTE_PLAIN_STATE)); @@ -753,6 +938,54 @@ void Init_gsl_monte(VALUE module) rb_gsl_monte_integrate, -1); rb_define_method(cgsl_monte_vegas, "integrate", rb_gsl_monte_vegas_integrate, -1); + +#ifdef GSL_1_13_LATER + cgsl_monte_miser_params = rb_define_class_under(cgsl_monte_miser, "Params", cGSL_Object); + cgsl_monte_vegas_params = rb_define_class_under(cgsl_monte_vegas, "Params", cGSL_Object); + + rb_define_method(cgsl_monte_miser, "params_get", rb_gsl_monte_miser_params_get, 0); + rb_define_method(cgsl_monte_miser, "params_set", rb_gsl_monte_miser_params_set, 1); + rb_define_method(cgsl_monte_miser_params, "estimate_frac", rb_gsl_monte_miser_params_get_estimate_frac, 0); + rb_define_method(cgsl_monte_miser_params, "set_estimate_frac", rb_gsl_monte_miser_params_set_estimate_frac, 1); + rb_define_alias(cgsl_monte_miser_params, "estimate_frac=", "set_estimate_frac"); + rb_define_method(cgsl_monte_miser_params, "min_calls", rb_gsl_monte_miser_params_get_min_calls, 0); + rb_define_method(cgsl_monte_miser_params, "set_min_calls", rb_gsl_monte_miser_params_set_min_calls, 1); + rb_define_alias(cgsl_monte_miser_params, "min_calls=", "set_min_calls"); + rb_define_method(cgsl_monte_miser_params, "min_calls_per_bisection", rb_gsl_monte_miser_params_get_min_calls_per_bisection, 0); + rb_define_method(cgsl_monte_miser_params, "set_min_calls_per_bisection", rb_gsl_monte_miser_params_set_min_calls_per_bisection, 1); + rb_define_alias(cgsl_monte_miser_params, "min_calls_per_bisection=", "set_min_calls_per_bisection"); + rb_define_method(cgsl_monte_miser_params, "alpha", rb_gsl_monte_miser_params_get_alpha, 0); + rb_define_method(cgsl_monte_miser_params, "set_alpha", rb_gsl_monte_miser_params_set_alpha, 1); + rb_define_alias(cgsl_monte_miser_params, "alpha=", "set_alpha"); + rb_define_method(cgsl_monte_miser_params, "dither", rb_gsl_monte_miser_params_get_dither, 0); + rb_define_method(cgsl_monte_miser_params, "set_dither", rb_gsl_monte_miser_params_set_dither, 1); + rb_define_alias(cgsl_monte_miser_params, "dither=", "set_dither"); + + rb_define_method(cgsl_monte_vegas, "params_get", rb_gsl_monte_vegas_params_get, 0); + rb_define_method(cgsl_monte_vegas, "params_set", rb_gsl_monte_vegas_params_set, 1); + rb_define_method(cgsl_monte_vegas_params, "alpha", rb_gsl_monte_vegas_params_get_alpha, 0); + rb_define_method(cgsl_monte_vegas_params, "set_alpha", rb_gsl_monte_vegas_params_set_alpha, 1); + rb_define_alias(cgsl_monte_vegas_params, "alpha=", "set_alpha"); + rb_define_method(cgsl_monte_vegas_params, "iterations", rb_gsl_monte_vegas_params_get_iterations, 0); + rb_define_method(cgsl_monte_vegas_params, "set_iterations", rb_gsl_monte_vegas_params_set_iterations, 1); + rb_define_alias(cgsl_monte_vegas_params, "iterations=", "set_iterations"); + rb_define_method(cgsl_monte_vegas_params, "stage", rb_gsl_monte_vegas_params_get_stage, 0); + rb_define_method(cgsl_monte_vegas_params, "set_stage", rb_gsl_monte_vegas_params_set_stage, 1); + rb_define_alias(cgsl_monte_vegas_params, "stage=", "set_stage"); + rb_define_method(cgsl_monte_vegas_params, "mode", rb_gsl_monte_vegas_params_get_mode, 0); + rb_define_method(cgsl_monte_vegas_params, "set_mode", rb_gsl_monte_vegas_params_set_mode, 1); + rb_define_alias(cgsl_monte_vegas_params, "mode=", "set_mode"); + rb_define_method(cgsl_monte_vegas_params, "verbose", rb_gsl_monte_vegas_params_get_verbose, 0); + rb_define_method(cgsl_monte_vegas_params, "set_verbose", rb_gsl_monte_vegas_params_set_verbose, 1); + rb_define_alias(cgsl_monte_vegas_params, "verbose=", "set_verbose"); + + rb_define_method(cgsl_monte_vegas, "runval", rb_gsl_monte_vegas_runval, 0); + + rb_define_const(cgsl_monte_vegas, "MODE_IMPORTANCE", INT2FIX(GSL_VEGAS_MODE_IMPORTANCE)); + rb_define_const(cgsl_monte_vegas, "MODE_IMPORTANCE_ONLY", INT2FIX(GSL_VEGAS_MODE_IMPORTANCE_ONLY)); + rb_define_const(cgsl_monte_vegas, "MODE_STRATIFIED", INT2FIX(GSL_VEGAS_MODE_STRATIFIED)); +#endif + } #ifdef CHECK_MONTE_FUNCTION #undef CHECK_MONTE_FUNCTION ext/monte.c @@ -38,7 +38,10 @@ enum { GSL_FDFMINIMIZER_STEEPEST_DESCENT, GSL_FMINIMIZER_NMSIMPLEX, #ifdef GSL_1_9_LATER - GSL_FDFMINIMIZER_VECTOR_BFGS2, + GSL_FDFMINIMIZER_VECTOR_BFGS2, +#endif +#ifdef GSL_1_13_LATER + GSL_FMINIMIZER_NMSIMPLEX2RAND, #endif }; @@ -446,6 +449,10 @@ static void define_const(VALUE klass1, VALUE klass2) rb_define_const(klass1, "VECTOR_BFGS2", INT2FIX(GSL_FDFMINIMIZER_VECTOR_BFGS2)); #endif +#ifdef GSL_1_13_LATER + rb_define_const(klass2, + "NMSIMPLEX2RAND", INT2FIX(GSL_FMINIMIZER_NMSIMPLEX2RAND)); +#endif } static const gsl_multimin_fdfminimizer_type* get_fdfminimizer_type(VALUE t) @@ -606,13 +613,21 @@ static const gsl_multimin_fminimizer_type* get_fminimizer_type(VALUE t) strcpy(name, STR2CSTR(t)); if (str_tail_grep(name, "nmsimplex") == 0) return gsl_multimin_fminimizer_nmsimplex; +#ifdef GSL_1_13_LATER + if (str_tail_grep(name, "nmsimplex2rand") == 0) + return gsl_multimin_fminimizer_nmsimplex2rand; +#endif else - rb_raise(rb_eTypeError, "unknown type %s (only nmsimplex supported)", name); + rb_raise(rb_eTypeError, "unknown type %s (nmsimplex and nmsimplex2rand supported)", name); break; case T_FIXNUM: switch (FIX2INT(t)) { case GSL_FMINIMIZER_NMSIMPLEX: return gsl_multimin_fminimizer_nmsimplex; break; +#ifdef GSL_1_13_LATER + case GSL_FMINIMIZER_NMSIMPLEX2RAND: + return gsl_multimin_fminimizer_nmsimplex2rand; break; +#endif default: rb_raise(rb_eTypeError, "%d: unknown type (not supported)", FIX2INT(t)); break; ext/multimin.c @@ -1624,6 +1624,86 @@ static VALUE rb_gsl_poly_wfit(int argc, VALUE *argv, VALUE obj) } #endif +#ifdef BASE_DOUBLE +#ifdef GSL_1_13_LATER +static VALUE rb_gsl_poly_eval_derivs_singleton(int argc, VALUE *argv, VALUE klass) +{ + VALUE ary; + gsl_vector *v = NULL, *v2 = NULL; + size_t i, lenc, lenres; +#ifdef HAVE_NARRAY_H + struct NARRAY *na; + double *ptr1, *ptr2; + int shape[1]; +#endif + + if (argc < 2) rb_raise(rb_eArgError, "Wrong number of arguments (%d for >= 2)"); + if (rb_obj_is_kind_of(argv[0], rb_cArray)) { + v = gsl_vector_alloc(RARRAY_LEN(argv[0])); + lenc = v->size; + for (i = 0; i < lenc; i++) { + gsl_vector_set(v, i, NUM2DBL(rb_ary_entry(argv[0], i))); + } + if (argc == 2) lenres = lenc + 1; + else lenres = FIX2INT(argv[2]); + v2 = gsl_vector_alloc(lenres); + gsl_poly_eval_derivs(v->data, lenc, NUM2DBL(argv[1]), v2->data, lenres); + ary = rb_ary_new2(lenres); + for (i = 0; i < lenres; i++) { + rb_ary_store(ary, i, rb_float_new(gsl_vector_get(v2, i))); + } + gsl_vector_free(v2); + gsl_vector_free(v); + return ary; + } + if (rb_obj_is_kind_of(argv[0], cgsl_vector)) { + Data_Get_Struct(argv[0], gsl_vector, v); + lenc = v->size; + if (argc == 2) lenres = lenc + 1; + else lenres = FIX2INT(argv[2]); + v2 = gsl_vector_alloc(lenres); + gsl_poly_eval_derivs(v->data, lenc, NUM2DBL(argv[1]), v2->data, lenres); + return Data_Wrap_Struct(cgsl_poly, 0, gsl_vector_free, v2); + } +#ifdef HAVE_NARRAY_H + if (NA_IsNArray(argv[0])) { + GetNArray(argv[0], na); + ptr1 = (double*) na->ptr; + lenc = na->total; + if (argc == 2) lenres = lenc + 1; + else lenres = FIX2INT(argv[2]); + shape[0] = lenres; + ary = na_make_object(NA_DFLOAT, na->rank, shape, CLASS_OF(argv[0])); + ptr2 = NA_PTR_TYPE(ary,double*); + gsl_poly_eval_derivs(ptr1, lenc, NUM2DBL(argv[1]), ptr2, lenres); + return ary; + } +#endif + return Qnil; // Never comes here +} +static VALUE rb_gsl_poly_eval_derivs(int argc, VALUE *argv, VALUE obj) +{ + gsl_vector *v, *v2; + size_t lenc, lenres; + Data_Get_Struct(obj, gsl_vector, v); + lenc = v->size; + switch (argc) { + case 1: + lenres = lenc + 1; + break; + case 2: + lenres = FIX2INT(argv[1]); + break; + default: + rb_raise(rb_eArgError, "Wrong number of arguments (%d for > 1)"); + } + v2 = gsl_vector_alloc(lenres); + gsl_poly_eval_derivs(v->data, lenc, NUM2DBL(argv[0]), v2->data, lenres); + return Data_Wrap_Struct(cgsl_poly, 0, gsl_vector_free, v2); +} +#endif +#endif + void FUNCTION(Init_gsl_poly,init)(VALUE module) { #ifdef BASE_DOUBLE @@ -1775,6 +1855,12 @@ void FUNCTION(Init_gsl_poly,init)(VALUE module) FUNCTION(rb_gsl_poly,fit), -1); rb_define_singleton_method(GSL_TYPE(cgsl_poly), "wfit", FUNCTION(rb_gsl_poly,wfit), -1); + +#ifdef GSL_1_13_LATER + rb_define_singleton_method(cgsl_poly, "eval_derivs", rb_gsl_poly_eval_derivs_singleton, -1); + rb_define_method(cgsl_vector, "eval_derivs", rb_gsl_poly_eval_derivs, -1); +#endif + #endif } ext/poly_source.c @@ -55,17 +55,21 @@ The function to minimize is given as an instance of the <a href="function.html"> <dt><a name="label-3" id="label-3"><code>GSL::Min::FMinimizer.alloc(<var>t</var>)</code></a></dt> <dd> <p>These method create an instance of the <code>GSL::Min::FMinimizer</code> class of -type <var>t</var>. The type <var>t</var> is given by a Ruby constant,</p> +type <var>t</var>. The type <var>t</var> is given by a String,</p> +<ul> +<li>"goldensection"</li> +<li>"brent"</li> +<li>"quad_golden"</li> +</ul> +<p>or by a Ruby constant,</p> <ul> <li>GSL::Min::FMinimizer::GOLDENSECTION</li> <li>GSL::Min::FMinimizer::BRENT</li> +<li>GSL::Min::FMinimizer::QUAD_GOLDEN (GSL-1.13)</li> </ul> -<p>ex1)</p> -<pre>include GSL -s1 = Min::FMinimizer.alloc(Min::FMinimizer::GOLDENSECTION)</pre> -<p>ex2)</p> +<p>ex)</p> <pre>include GSL::Min -s2 = FMinimizer.alloc(FMinimizer::BRENT)</pre></dd> +s = FMinimizer.alloc(FMinimizer::BRENT)</pre></dd> <dt><a name="label-4" id="label-4"><code>GSL::Min::FMinimizer#set(<var>f</var>, <var>xmin</var>, <var>xlow</var>, <var>xup</var>)</code></a></dt> <dd> <p>This method sets, or resets, an existing minimizer <var>self</var> to use html/min.html @@ -94,18 +94,90 @@ finishes.</p> <dt><a name="label-40" id="label-40"><code>GSL::Monte::Vegas#result</code></a></dt> <dt><a name="label-41" id="label-41"><code>GSL::Monte::Vegas#sigma</code></a></dt> <dt><a name="label-42" id="label-42"><code>GSL::Monte::Vegas#chisq</code></a></dt> -<dt><a name="label-43" id="label-43"><code>GSL::Monte::Vegas#iterations</code></a></dt> -<dt><a name="label-44" id="label-44"><code>GSL::Monte::Vegas#iterations=</code></a></dt> -<dt><a name="label-45" id="label-45"><code>GSL::Monte::Vegas#alpha</code></a></dt> -<dt><a name="label-46" id="label-46"><code>GSL::Monte::Vegas#alpha=</code></a></dt> -<dt><a name="label-47" id="label-47"><code>GSL::Monte::Vegas#stage</code></a></dt> -<dt><a name="label-48" id="label-48"><code>GSL::Monte::Vegas#stage=</code></a></dt> -<dt><a name="label-49" id="label-49"><code>GSL::Monte::Vegas#mode</code></a></dt> -<dt><a name="label-50" id="label-50"><code>GSL::Monte::Vegas#mode=</code></a></dt> -<dt><a name="label-51" id="label-51"><code>GSL::Monte::Vegas#verbose</code></a></dt> -<dt><a name="label-52" id="label-52"><code>GSL::Monte::Vegas#verbose=</code></a></dt> +<dd> +Returns the chi-squared per degree of freedom for the weighted estimate of the integral. The returned value should be close to 1. A value which differs significantly from 1 indicates that the values from different iterations are inconsistent. In this case the weighted error will be under-estimated, and further iterations of the algorithm are needed to obtain reliable results.</dd> +<dt><a name="label-43" id="label-43"><code>GSL::Monte::Vegas#runval</code></a></dt> +<dd> +Returns the raw (unaveraged) values of the integral and its error <var>[result, sigma]</var> from the most recent iteration of the algorithm.</dd> +<dt><a name="label-44" id="label-44"><code>GSL::Monte::Vegas#iterations</code></a></dt> +<dt><a name="label-45" id="label-45"><code>GSL::Monte::Vegas#iterations=</code></a></dt> +<dt><a name="label-46" id="label-46"><code>GSL::Monte::Vegas#alpha</code></a></dt> +<dt><a name="label-47" id="label-47"><code>GSL::Monte::Vegas#alpha=</code></a></dt> +<dt><a name="label-48" id="label-48"><code>GSL::Monte::Vegas#stage</code></a></dt> +<dt><a name="label-49" id="label-49"><code>GSL::Monte::Vegas#stage=</code></a></dt> +<dt><a name="label-50" id="label-50"><code>GSL::Monte::Vegas#mode</code></a></dt> +<dt><a name="label-51" id="label-51"><code>GSL::Monte::Vegas#mode=</code></a></dt> +<dt><a name="label-52" id="label-52"><code>GSL::Monte::Vegas#verbose</code></a></dt> +<dt><a name="label-53" id="label-53"><code>GSL::Monte::Vegas#verbose=</code></a></dt> +</dl> +<h2><a name="5" href="#5">5</a> <a name="label-54">Miser Parameters (GSL-1.13 or later)</a></h2> + +<dl> +<dt><a name="label-55" id="label-55"><code>GSL::Monte::Miser#params_get</code></a></dt> +<dd> +Returns the parameters of the integrator state as an instance of <code>GSL::Monte::Miser::Params</code> class.</dd> +<dt><a name="label-56" id="label-56"><code>GSL::Monte::Miser#params_set(<var>params</var>)</code></a></dt> +<dd> +Sets the integrator parameters based on values provided in an object of the <code>GSL::Monte::Miser::Params</code> class <var>params</var>.</dd> +</dl> +<h3><a name="5.1" href="#5.1">5.1</a> <a name="label-57">Accessors of <code>GSL::Monte::Miser::Params</code></a></h3> + +<dl> +<dt><a name="label-58" id="label-58"><code>GSL::Monte::Miser::Params#estimate_frac</code></a></dt> +<dt><a name="label-59" id="label-59"><code>GSL::Monte::Miser::Params#estimate_frac=</code></a></dt> +<dd> +The fraction of the currently available number of function calls which are allocated to estimating the variance at each recursive step. The default value is 0.1.</dd> +<dt><a name="label-60" id="label-60"><code>GSL::Monte::Miser::Params#min_calls</code></a></dt> +<dt><a name="label-61" id="label-61"><code>GSL::Monte::Miser::Params#min_calls=</code></a></dt> +<dd> +The minimum number of function calls required for each estimate of the variance. If the number of function calls allocated to the estimate using <var>estimate_frac</var> falls below <var>min_calls</var> then <var>min_calls</var> are used instead. This ensures that each estimate maintains a reasonable level of accuracy. The default value of min_calls is 16 * dim.</dd> +<dt><a name="label-62" id="label-62"><code>GSL::Monte::Miser::Params#min_calls_per_bisection</code></a></dt> +<dt><a name="label-63" id="label-63"><code>GSL::Monte::Miser::Params#min_calls_per_bisection=</code></a></dt> +<dd> +The minimum number of function calls required to proceed with a bisection step. When a recursive step has fewer calls available than <var>min_calls_per_bisection</var> it performs a plain Monte Carlo estimate of the current sub-region and terminates its branch of the recursion. The default value of this parameter is 32 * min_calls.</dd> +<dt><a name="label-64" id="label-64"><code>GSL::Monte::Miser::Params#alpha</code></a></dt> +<dt><a name="label-65" id="label-65"><code>GSL::Monte::Miser::Params#alpha=</code></a></dt> +<dd> +This parameter controls how the estimated variances for the two sub-regions of a bisection are combined when allocating points. With recursive sampling the overall variance should scale better than 1/N, since the values from the sub-regions will be obtained using a procedure which explicitly minimizes their variance. To accommodate this behavior the MISER algorithm allows the total variance to depend on a scaling parameter <var>alpha</var>, The authors of the original paper describing MISER recommend the value <var>alpha</var> = 2 as a good choice, obtained from numerical experiments, and this is used as the default value in this implementation.</dd> +<dt><a name="label-66" id="label-66"><code>GSL::Monte::Miser::Params#dither</code></a></dt> +<dt><a name="label-67" id="label-67"><code>GSL::Monte::Miser::Params#dither=</code></a></dt> +</dl> +<p>This parameter introduces a random fractional variation of size dither into each bisection, which can be used to break the symmetry of integrands which are concentrated near the exact center of the hypercubic integration region. The default value of dither is zero, so no variation is introduced. If needed, a typical value of dither is 0.1.</p> +<h2><a name="6" href="#6">6</a> <a name="label-68">Vegas Parameters (GSL-1.13 or later)</a></h2> + +<dl> +<dt><a name="label-69" id="label-69"><code>GSL::Monte::Vegas#params_get</code></a></dt> +<dd> +Returns the parameters of the integrator state as an instance of <code>GSL::Monte::Vegas::Params</code> class.</dd> +<dt><a name="label-70" id="label-70"><code>GSL::Monte::Vegas#params_set(<var>params</var>)</code></a></dt> +<dd> +Sets the integrator parameters based on values provided in an object of the <code>GSL::Monte::Vegas::Params</code> class <var>params</var>.</dd> +</dl> +<h3><a name="6.1" href="#6.1">6.1</a> <a name="label-71">Accessors of <code>GSL::Monte::Vegas::Params</code></a></h3> + +<dl> +<dt><a name="label-72" id="label-72"><code>GSL::Monte::Vegas::Params#alpha</code></a></dt> +<dt><a name="label-73" id="label-73"><code>GSL::Monte::Vegas::Params#alpha=</code></a></dt> +<dd> +Controls the stiffness of the rebinning algorithm. It is typically set between one and two. A value of zero prevents rebinning of the grid. The default value is 1.5.</dd> +<dt><a name="label-74" id="label-74"><code>GSL::Monte::Vegas::Params#iterations</code></a></dt> +<dt><a name="label-75" id="label-75"><code>GSL::Monte::Vegas::Params#iterations=</code></a></dt> +<dd> +The number of iterations to perform for each call to the routine. The default value is 5 iterations.</dd> +<dt><a name="label-76" id="label-76"><code>GSL::Monte::Vegas::Params#stage</code></a></dt> +<dt><a name="label-77" id="label-77"><code>GSL::Monte::Vegas::Params#stage=</code></a></dt> +<dd> +Setting this determines the stage of the calculation. Normally, stage = 0 which begins with a new uniform grid and empty weighted average. Calling vegas with stage = 1 retains the grid from the previous run but discards the weighted average, so that one can "tune" the grid using a relatively small number of points and then do a large run with stage = 1 on the optimized grid. Setting stage = 2 keeps the grid and the weighted average from the previous run, but may increase (or decrease) the number of histogram bins in the grid depending on the number of calls available. Choosing stage = 3 enters at the main loop, so that nothing is changed, and is equivalent to performing additional iterations in a previous call.</dd> +<dt><a name="label-78" id="label-78"><code>GSL::Monte::Vegas::Params#mode</code></a></dt> +<dt><a name="label-79" id="label-79"><code>GSL::Monte::Vegas::Params#mode=</code></a></dt> +<dd> +The possible choices are <code>GSL::VEGAS::MODE_IMPORTANCE</code>, <code>GSL::VEGAS::MODE_STRATIFIED</code>, <code>GSL::VEGAS::MODE_IMPORTANCE_ONLY</code>. This determines whether VEGAS will use importance sampling or stratified sampling, or whether it can pick on its own. In low dimensions VEGAS uses strict stratified sampling (more precisely, stratified sampling is chosen if there are fewer than 2 bins per box).</dd> +<dt><a name="label-80" id="label-80"><code>GSL::Monte::Vegas::Params#verbose</code></a></dt> +<dt><a name="label-81" id="label-81"><code>GSL::Monte::Vegas::Params#verbose=</code></a></dt> +<dd> +Set the level of information printed by VEGAS. All information is written to the stream ostream. The default setting of verbose is -1, which turns off all output. A verbose value of 0 prints summary information about the weighted average and final result, while a value of 1 also displays the grid coordinates. A value of 2 prints information from the rebinning procedure for each iteration.</dd> </dl> -<h2><a name="5" href="#5">5</a> <a name="label-53">Example</a></h2> +<h2><a name="7" href="#7">7</a> <a name="label-82">Example</a></h2> <pre>#!/usr/bin/env ruby require("gsl") html/monte.html @@ -81,7 +81,7 @@ the function in question. </p> <dt><a name="label-5" id="label-5"><code>GSL::MultiMin::FMinimizer.alloc(<var>type</var>, <var>n</var>)</code></a></dt> <dd> These method create a minimizer of type <var>type</var> for an <var>n</var>-dimension function. -The type is given by a name of the algorithm, or by a Ruby constant. +The type is given by a string, or by a Ruby constant. <ul> <li><code>GSL::MultiMin::FdfMinimizer::CONJUGATE_FR</code> or <code>"conjugate_fr"</code></li> <li><code>GSL::MultiMin::FdfMinimizer::CONJUGATE_PR</code> or <code>"conjugate_pr"</code></li> @@ -89,6 +89,7 @@ The type is given by a name of the algorithm, or by a Ruby constant. <li><code>GSL::MultiMin::FdfMinimizer::VECTOR_BFGS2</code> or <code>"vector_bfgs2"</code> (GSL-1.9 or later)</li> <li><code>GSL::MultiMin::FdfMinimizer::STEEPEST_DESCENT</code> or <code>"steepest_descent"</code></li> <li><code>GSL::MultiMin::FMinimizer::NMSIMPLEX</code> or <code>"nmsimplex"</code></li> +<li><code>GSL::MultiMin::FMinimizer::NMSIMPLEX2RAND</code> or <code>"nmsimplex2rand"</code> (GSL-1.13)</li> <li><p>ex:</p> <pre>include GSL::MultiMin m1 = FdfMinimizer.alloc(FdfMinimizer::CONJUGATE_FR, 2) html/multimin.html @@ -59,14 +59,45 @@ irb(main):006:0> GSL::Poly.eval([1, 2, 3], GSL::Vector[1, 2, 3]) irb(main):007:0> GSL::Poly.eval([1, 2, 3], NArray[1.0, 2, 3]) => NArray.float(3): [ 6.0, 17.0, 34.0 ]</pre></dd> +<dt><a name="label-3" id="label-3"><code>GSL::Poly.eval_derivs(<var>c</var>, <var>x</var>)</code></a></dt> +<dt><a name="label-4" id="label-4"><code>GSL::Poly.eval_derivs(<var>c</var>, <var>x</var>, <var>lenres</var>)</code></a></dt> +<dd> +(GSL-1.13) Evaluate and return a polynomial and its derivatives. The output contains the values of d^k P/d x^k for the specified value of x starting with k = 0. The input polynomial <var>c</var> can be an <code>Array</code>, <code>GSL::Poly</code> or an <code>NArray</code>. If <var>lenres</var> is not given, <var>lenres = LENGTH(c) + 1</var> is used, therefore the last element of the output is 0.</dd> +<dt><a name="label-5" id="label-5"><code>GSL::Poly#eval_derivs(<var>x</var>)</code></a></dt> +<dt><a name="label-6" id="label-6"><code>GSL::Poly#eval_derivs(<var>x</var>, <var>lenres</var>)</code></a></dt> +<dd> +<p>(GSL-1.13) Evaluate and return a polynomial and its derivatives. The output contains the values of d^k P/d x^k for the specified value of x starting with k = 0. If <var>lenres</var> is not given, <var>lenres = LENGTH(self) + 1</var> is used, therefore the last element of the output is 0.</p> +<p>Ex.)</p> +<pre>>> ary = [1, 2, 3] +=> [1, 2, 3] +>> GSL::Poly.eval_derivs(ary, 1) +=> [6.0, 8.0, 6.0, 0.0] +>> na = NArray[1.0, 2, 3] +=> NArray.float(3): +[ 1.0, 2.0, 3.0 ] +>> GSL::Poly.eval_derivs(na, 1) +=> NArray.float(4): +[ 6.0, 8.0, 6.0, 0.0 ] +>> poly = GSL::Poly[1.0, 2, 3] +=> GSL::Poly +[ 1.000e+00 2.000e+00 3.000e+00 ] +>> GSL::Poly.eval_derivs(poly, 1) +=> GSL::Poly +[ 6.000e+00 8.000e+00 6.000e+00 0.000e+00 ] +>> poly.eval_derivs(1) +=> GSL::Poly +[ 6.000e+00 8.000e+00 6.000e+00 0.000e+00 ] +>> poly.eval_derivs(1, 3) +=> GSL::Poly +[ 6.000e+00 8.000e+00 6.000e+00 ]</pre></dd> </dl> -<h2><a name="2" href="#2">2</a> <a name="label-3">Solving polynomial equations</a></h2> +<h2><a name="2" href="#2">2</a> <a name="label-7">Solving polynomial equations</a></h2> -<h3><a name="2.1" href="#2.1">2.1</a> <a name="label-4">Quadratic Equations</a></h3> +<h3><a name="2.1" href="#2.1">2.1</a> <a name="label-8">Quadratic Equations</a></h3> <dl> -<dt><a name="label-5" id="label-5"><code>GSL::Poly::solve_quadratic(<var>a</var>, <var>b</var>, <var>c</var>)</code></a></dt> -<dt><a name="label-6" id="label-6"><code>GSL::Poly::solve_quadratic([<var>a</var>, <var>b</var>, <var>c</var>])</code></a></dt> +<dt><a name="label-9" id="label-9"><code>GSL::Poly::solve_quadratic(<var>a</var>, <var>b</var>, <var>c</var>)</code></a></dt> +<dt><a name="label-10" id="label-10"><code>GSL::Poly::solve_quadratic([<var>a</var>, <var>b</var>, <var>c</var>])</code></a></dt> <dd> <p>Find the real roots of the quadratic equation,</p> <pre>a x^2 + b x + c = 0</pre> @@ -78,8 +109,8 @@ or a <code>GSL::Vector</code> object. The roots are returned as a <code>GSL::Vec => GSL::Vector: [ 1.000e+00 2.000e+00 ]</pre></li> </ul></dd> -<dt><a name="label-7" id="label-7"><code>GSL::Poly::complex_solve_quadratic(<var>a</var>, <var>b</var>, <var>c</var>)</code></a></dt> -<dt><a name="label-8" id="label-8"><code>GSL::Poly::complex_solve_quadratic([<var>a</var>, <var>b</var>, <var>c</var>])</code></a></dt> +<dt><a name="label-11" id="label-11"><code>GSL::Poly::complex_solve_quadratic(<var>a</var>, <var>b</var>, <var>c</var>)</code></a></dt> +<dt><a name="label-12" id="label-12"><code>GSL::Poly::complex_solve_quadratic([<var>a</var>, <var>b</var>, <var>c</var>])</code></a></dt> <dd> <p>Find the complex roots of the quadratic equation,</p> <pre>a z^2 + b z + z = 0</pre> @@ -98,23 +129,23 @@ irb(main):003:0> GSL::Poly::complex_solve_quadratic(1, -3, 2).real <--- R [ 1.000e+00 2.000e+00 ]</pre></li> </ul></dd> </dl> -<h3><a name="2.2" href="#2.2">2.2</a> <a name="label-9">Cubic Equations</a></h3> +<h3><a name="2.2" href="#2.2">2.2</a> <a name="label-13">Cubic Equations</a></h3> <dl> -<dt><a name="label-10" id="label-10"><code>GSL::Poly::solve_cubic(<var>same</var> <var>as</var> <var>solve_quadratic</var>)</code></a></dt> +<dt><a name="label-14" id="label-14"><code>GSL::Poly::solve_cubic(<var>same</var> <var>as</var> <var>solve_quadratic</var>)</code></a></dt> <dd> <p>This method finds the real roots of the cubic equation,</p> <pre>x^3 + a x^2 + b x + c = 0</pre></dd> -<dt><a name="label-11" id="label-11"><code>GSL::Poly::complex_solve_cubic(<var>same</var> <var>as</var> <var>solve_cubic</var>)</code></a></dt> +<dt><a name="label-15" id="label-15"><code>GSL::Poly::complex_solve_cubic(<var>same</var> <var>as</var> <var>solve_cubic</var>)</code></a></dt> <dd> <p>This method finds the complex roots of the cubic equation,</p> <pre>z^3 + a z^2 + b z + c = 0</pre></dd> </dl> -<h3><a name="2.3" href="#2.3">2.3</a> <a name="label-12">General Polynomial Equations</a></h3> +<h3><a name="2.3" href="#2.3">2.3</a> <a name="label-16">General Polynomial Equations</a></h3> <dl> -<dt><a name="label-13" id="label-13"><code>GSL::Poly::complex_solve(<var>c0</var>, <var>c1</var>, <var>c2</var>,,, )</code></a></dt> -<dt><a name="label-14" id="label-14"><code>GSL::Poly::solve(<var>c0</var>, <var>c1</var>, <var>c2</var>,,, )</code></a></dt> +<dt><a name="label-17" id="label-17"><code>GSL::Poly::complex_solve(<var>c0</var>, <var>c1</var>, <var>c2</var>,,, )</code></a></dt> +<dt><a name="label-18" id="label-18"><code>GSL::Poly::solve(<var>c0</var>, <var>c1</var>, <var>c2</var>,,, )</code></a></dt> <dd> Find the complex roots of the polynomial equation. Note that the coefficients are given by "ascending" order. @@ -125,14 +156,14 @@ the coefficients are given by "ascending" order. => #<GSL::Vector::Complex:0x75e614></pre></li> </ul></dd> </dl> -<h2><a name="3" href="#3">3</a> <a name="label-15">GSL::Poly Class</a></h2> +<h2><a name="3" href="#3">3</a> <a name="label-19">GSL::Poly Class</a></h2> <p>This class expresses polynomials of arbitrary orders.</p> -<h3><a name="3.1" href="#3.1">3.1</a> <a name="label-16">Constructors</a></h3> +<h3><a name="3.1" href="#3.1">3.1</a> <a name="label-20">Constructors</a></h3> <dl> -<dt><a name="label-17" id="label-17"><code>GSL::Poly.alloc(<var>c0</var>, <var>c1</var>, <var>c2</var>, ....)</code></a></dt> -<dt><a name="label-18" id="label-18"><code>GSL::Poly[c0, <var>c1</var>, <var>c2</var>, ....]</code></a></dt> +<dt><a name="label-21" id="label-21"><code>GSL::Poly.alloc(<var>c0</var>, <var>c1</var>, <var>c2</var>, ....)</code></a></dt> +<dt><a name="label-22" id="label-22"><code>GSL::Poly[c0, <var>c1</var>, <var>c2</var>, ....]</code></a></dt> <dd> <p>This creates an instance of the <code>GSL::Poly</code> class, which represents a polynomial</p> @@ -143,17 +174,17 @@ which represents a polynomial</p> <pre>poly = GSL::Poly.alloc([2, -3, 1])</pre></li> </ul></dd> </dl> -<h3><a name="3.2" href="#3.2">3.2</a> <a name="label-19">Instance Methods</a></h3> +<h3><a name="3.2" href="#3.2">3.2</a> <a name="label-23">Instance Methods</a></h3> <dl> -<dt><a name="label-20" id="label-20"><code>GSL::Poly#eval(<var>x</var>)</code></a></dt> -<dt><a name="label-21" id="label-21"><code>GSL::Poly#at(<var>x</var>)</code></a></dt> +<dt><a name="label-24" id="label-24"><code>GSL::Poly#eval(<var>x</var>)</code></a></dt> +<dt><a name="label-25" id="label-25"><code>GSL::Poly#at(<var>x</var>)</code></a></dt> <dd> Evaluates the polynomial c[0] + c[1] x + c[2] x^2 + ... + c[len-1] x^{len-1} using Horner's method for stability. The argument <var>x</var> is a <code>Numeric</code>, <code>GSL::Vector, Matrix</code> or an <code>Array</code>.</dd> -<dt><a name="label-22" id="label-22"><code>GSL::Poly#solve_quadratic</code></a></dt> +<dt><a name="label-26" id="label-26"><code>GSL::Poly#solve_quadratic</code></a></dt> <dd> Solve the quadratic equation. <ul> @@ -165,12 +196,12 @@ irb(main):005:0> a.solve_quadratic => GSL::Vector: [ 1.000e+00 2.000e+00 ]</pre></li> </ul></dd> -<dt><a name="label-23" id="label-23"><code>GSL::Poly#solve_cubic</code></a></dt> +<dt><a name="label-27" id="label-27"><code>GSL::Poly#solve_cubic</code></a></dt> <dd> Solve the cubic equation.</dd> -<dt><a name="label-24" id="label-24"><code>GSL::Poly#complex_solve</code></a></dt> -<dt><a name="label-25" id="label-25"><code>GSL::Poly#solve</code></a></dt> -<dt><a name="label-26" id="label-26"><code>GSL::Poly#roots</code></a></dt> +<dt><a name="label-28" id="label-28"><code>GSL::Poly#complex_solve</code></a></dt> +<dt><a name="label-29" id="label-29"><code>GSL::Poly#solve</code></a></dt> +<dt><a name="label-30" id="label-30"><code>GSL::Poly#roots</code></a></dt> <dd> <p>These methods find the complex roots of the quadratic equation,</p> <pre>c0 + c1 z + c2 z^2 + .... = 0</pre> @@ -184,11 +215,11 @@ irb(main):010:0> a.solve => #<GSL::Vector::Complex:0x35db28></pre></li> </ul></dd> </dl> -<h2><a name="4" href="#4">4</a> <a name="label-27">Polynomial fitting</a></h2> +<h2><a name="4" href="#4">4</a> <a name="label-31">Polynomial fitting</a></h2> <dl> -<dt><a name="label-28" id="label-28"><code>GSL::Poly.fit(<var>x</var>, <var>y</var>, <var>order</var>)</code></a></dt> -<dt><a name="label-29" id="label-29"><code>GSL::Poly.wfit(<var>x</var>, <var>w</var>, <var>y</var>, <var>order</var>)</code></a></dt> +<dt><a name="label-32" id="label-32"><code>GSL::Poly.fit(<var>x</var>, <var>y</var>, <var>order</var>)</code></a></dt> +<dt><a name="label-33" id="label-33"><code>GSL::Poly.wfit(<var>x</var>, <var>w</var>, <var>y</var>, <var>order</var>)</code></a></dt> <dd> <p>Finds the coefficient of a polynomial of order <var>order</var> that fits the vector data (<var>x, y</var>) in a least-square sense. @@ -206,29 +237,29 @@ coef, cov, chisq, status = Poly.fit(x, y, 3) x2 = GSL::Vector.linspace(1, 5, 20) graph([x, y], [x2, coef.eval(x2)], "-C -g 3 -S 4")</pre></dd> </dl> -<h2><a name="5" href="#5">5</a> <a name="label-30">Divided-difference representations</a></h2> +<h2><a name="5" href="#5">5</a> <a name="label-34">Divided-difference representations</a></h2> <dl> -<dt><a name="label-31" id="label-31"><code>GSL::Poly::dd_init(<var>xa</var>, <var>ya</var>)</code></a></dt> +<dt><a name="label-35" id="label-35"><code>GSL::Poly::dd_init(<var>xa</var>, <var>ya</var>)</code></a></dt> <dd> This method computes a divided-difference representation of the interpolating polynomial for the points <var>(xa, ya)</var>.</dd> -<dt><a name="label-32" id="label-32"><code>GSL::Poly::DividedDifference#eval(<var>x</var>)</code></a></dt> +<dt><a name="label-36" id="label-36"><code>GSL::Poly::DividedDifference#eval(<var>x</var>)</code></a></dt> <dd> This method evaluates the polynomial stored in divided-difference form <var>self</var> at the point <var>x</var>.</dd> -<dt><a name="label-33" id="label-33"><code>GSL::Poly::DividedDifference#taylor(<var>xp</var>)</code></a></dt> +<dt><a name="label-37" id="label-37"><code>GSL::Poly::DividedDifference#taylor(<var>xp</var>)</code></a></dt> <dd> This method converts the divided-difference representation of a polynomial to a Taylor expansion. On output the Taylor coefficients of the polynomial expanded about the point <var>xp</var> are returned.</dd> </dl> -<h2><a name="6" href="#6">6</a> <a name="label-34">Extensions</a></h2> +<h2><a name="6" href="#6">6</a> <a name="label-38">Extensions</a></h2> -<h3><a name="6.1" href="#6.1">6.1</a> <a name="label-35">Special Polynomials</a></h3> +<h3><a name="6.1" href="#6.1">6.1</a> <a name="label-39">Special Polynomials</a></h3> <dl> -<dt><a name="label-36" id="label-36"><code>GSL::Poly.hermite(<var>n</var>)</code></a></dt> +<dt><a name="label-40" id="label-40"><code>GSL::Poly.hermite(<var>n</var>)</code></a></dt> <dd> <p>This returns coefficients of the <var>n</var>-th order Hermite polynomial, <var>H(x; n)</var>. For order of <var>n</var> >= 3, this method uses the recurrence relation</p> @@ -245,25 +276,25 @@ irb(main):015:0> GSL::Poly.hermite(7) => GSL::Poly::Int: [ 0 -1680 0 3360 0 -1344 0 128 ]</pre></li> </ul></dd> -<dt><a name="label-37" id="label-37"><code>GSL::Poly.cheb(<var>n</var>)</code></a></dt> -<dt><a name="label-38" id="label-38"><code>GSL::Poly.chebyshev(<var>n</var>)</code></a></dt> +<dt><a name="label-41" id="label-41"><code>GSL::Poly.cheb(<var>n</var>)</code></a></dt> +<dt><a name="label-42" id="label-42"><code>GSL::Poly.chebyshev(<var>n</var>)</code></a></dt> <dd> <p>Return the coefficients of the <var>n</var>-th order Chebyshev polynomial, <var>T(x; n</var>. For order of <var>n</var> >= 3, this method uses the recurrence relation</p> <pre>T(x; n+1) = 2 x T(x; n) - T(x; n-1)</pre></dd> -<dt><a name="label-39" id="label-39"><code>GSL::Poly.cheb_II(<var>n</var>)</code></a></dt> -<dt><a name="label-40" id="label-40"><code>GSL::Poly.chebyshev_II(<var>n</var>)</code></a></dt> +<dt><a name="label-43" id="label-43"><code>GSL::Poly.cheb_II(<var>n</var>)</code></a></dt> +<dt><a name="label-44" id="label-44"><code>GSL::Poly.chebyshev_II(<var>n</var>)</code></a></dt> <dd> <p>Return the coefficients of the <var>n</var>-th order Chebyshev polynomial of type II, <var>U(x; n</var>.</p> <pre>U(x; n+1) = 2 x U(x; n) - U(x; n-1)</pre></dd> -<dt><a name="label-41" id="label-41"><code>GSL::Poly.bell(<var>n</var>)</code></a></dt> +<dt><a name="label-45" id="label-45"><code>GSL::Poly.bell(<var>n</var>)</code></a></dt> <dd> Bell polynomial</dd> -<dt><a name="label-42" id="label-42"><code>GSL::Poly.bessel(<var>n</var>)</code></a></dt> +<dt><a name="label-46" id="label-46"><code>GSL::Poly.bessel(<var>n</var>)</code></a></dt> <dd> Bessel polynomial</dd> -<dt><a name="label-43" id="label-43"><code>GSL::Poly.laguerre(<var>n</var>)</code></a></dt> +<dt><a name="label-47" id="label-47"><code>GSL::Poly.laguerre(<var>n</var>)</code></a></dt> <dd> <p>Retunrs the coefficients of the <var>n</var>-th order Laguerre polynomial multiplied by n!.</p> @@ -286,15 +317,15 @@ irb(main):006:0> GSL::Poly.laguerre(4) => GSL::Poly::Int: [ 24 -96 72 -16 1 ] <--- (x^4 - 16x^3 + 72x^2 - 96x + 24)/4!</pre></dd> </dl> -<h3><a name="6.2" href="#6.2">6.2</a> <a name="label-44">Polynomial Operations</a></h3> +<h3><a name="6.2" href="#6.2">6.2</a> <a name="label-48">Polynomial Operations</a></h3> <dl> -<dt><a name="label-45" id="label-45"><code>GSL::Poly#conv</code></a></dt> -<dt><a name="label-46" id="label-46"><code>GSL::Poly#deconv</code></a></dt> -<dt><a name="label-47" id="label-47"><code>GSL::Poly#reduce</code></a></dt> -<dt><a name="label-48" id="label-48"><code>GSL::Poly#deriv</code></a></dt> -<dt><a name="label-49" id="label-49"><code>GSL::Poly#integ</code></a></dt> -<dt><a name="label-50" id="label-50"><code>GSL::Poly#compan</code></a></dt> +<dt><a name="label-49" id="label-49"><code>GSL::Poly#conv</code></a></dt> +<dt><a name="label-50" id="label-50"><code>GSL::Poly#deconv</code></a></dt> +<dt><a name="label-51" id="label-51"><code>GSL::Poly#reduce</code></a></dt> +<dt><a name="label-52" id="label-52"><code>GSL::Poly#deriv</code></a></dt> +<dt><a name="label-53" id="label-53"><code>GSL::Poly#integ</code></a></dt> +<dt><a name="label-54" id="label-54"><code>GSL::Poly#compan</code></a></dt> </dl> <p><a href="complex.html">prev</a> <a href="sf.html">next</a></p> html/poly.html @@ -14,7 +14,7 @@ <h2><a name="1" href="#1">1</a> <a name="label-1">Front Matter</a></h2> <p>This document describes the modules, classes and the methods of Ruby/GSL. -This includes cut&paste from the <a href="http://www.gnu.org/software/gsl/manual/html_node/">GSL Reference</a>, and documents of the extention libraries (section 41).</p> +This includes cut&paste from the <a href="http://www.gnu.org/software/gsl/manual/html_node/">GSL Reference</a>, and documents of the extention libraries.</p> <h3><a name="1.1" href="#1.1">1.1</a> <a name="label-2">Copyright of the GSL Reference</a></h3> <pre>Copyright 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 The GSL Team. @@ -25,6 +25,11 @@ and "Free Software Needs Free Documentation", the Front-Cover text being "A GNU and with the Back-Cover Text being (a) (see below). A copy of the license is included in the section entitled "GNU Free Documentation License". (a) The Back-Cover Text is: "You have the freedom to copy and modify this GNU Manual."</pre> +<h3><a name="1.2" href="#1.2">1.2</a> <a name="label-3">Copyright of this reference</a></h3> + +<pre>2009 Yoshiki Tsunesada +Permission is granted to copy, distribute and/or modify this document under the terms of +the GNU Free Documentation License.</pre> <h2><a name="2" href="#2">2</a> <a name="label-0">Ruby/GSL Reference</a></h2> <ol> html/ref.html @@ -47,15 +47,18 @@ The function to minimize is given as an instance of the ((<GSL::Function|URL:fun == GSL::Min::FMinimizer class --- GSL::Min::FMinimizer.alloc(t) These method create an instance of the (({GSL::Min::FMinimizer})) class of - type ((|t|)). The type ((|t|)) is given by a Ruby constant, + type ((|t|)). The type ((|t|)) is given by a String, + * "goldensection" + * "brent" + * "quad_golden" + or by a Ruby constant, * GSL::Min::FMinimizer::GOLDENSECTION * GSL::Min::FMinimizer::BRENT - ex1) - include GSL - s1 = Min::FMinimizer.alloc(Min::FMinimizer::GOLDENSECTION) - ex2) + * GSL::Min::FMinimizer::QUAD_GOLDEN (GSL-1.13) + + ex) include GSL::Min - s2 = FMinimizer.alloc(FMinimizer::BRENT) + s = FMinimizer.alloc(FMinimizer::BRENT) --- GSL::Min::FMinimizer#set(f, xmin, xlow, xup) This method sets, or resets, an existing minimizer ((|self|)) to use rd/min.rd @@ -70,6 +70,9 @@ The function to be integrated has its own datatype, the (({GSL::Monte::Function} --- GSL::Monte::Vegas#result --- GSL::Monte::Vegas#sigma --- GSL::Monte::Vegas#chisq + Returns the chi-squared per degree of freedom for the weighted estimate of the integral. The returned value should be close to 1. A value which differs significantly from 1 indicates that the values from different iterations are inconsistent. In this case the weighted error will be under-estimated, and further iterations of the algorithm are needed to obtain reliable results. +--- GSL::Monte::Vegas#runval + Returns the raw (unaveraged) values of the integral and its error ((|[result, sigma]|)) from the most recent iteration of the algorithm. --- GSL::Monte::Vegas#iterations --- GSL::Monte::Vegas#iterations= --- GSL::Monte::Vegas#alpha @@ -81,6 +84,51 @@ The function to be integrated has its own datatype, the (({GSL::Monte::Function} --- GSL::Monte::Vegas#verbose --- GSL::Monte::Vegas#verbose= +== Miser Parameters (GSL-1.13 or later) +--- GSL::Monte::Miser#params_get + Returns the parameters of the integrator state as an instance of (({GSL::Monte::Miser::Params})) class. +--- GSL::Monte::Miser#params_set(params) + Sets the integrator parameters based on values provided in an object of the (({GSL::Monte::Miser::Params})) class ((|params|)). +=== Accessors of (({GSL::Monte::Miser::Params})) +--- GSL::Monte::Miser::Params#estimate_frac +--- GSL::Monte::Miser::Params#estimate_frac= + The fraction of the currently available number of function calls which are allocated to estimating the variance at each recursive step. The default value is 0.1. +--- GSL::Monte::Miser::Params#min_calls +--- GSL::Monte::Miser::Params#min_calls= + The minimum number of function calls required for each estimate of the variance. If the number of function calls allocated to the estimate using ((|estimate_frac|)) falls below ((|min_calls|)) then ((|min_calls|)) are used instead. This ensures that each estimate maintains a reasonable level of accuracy. The default value of min_calls is 16 * dim. +--- GSL::Monte::Miser::Params#min_calls_per_bisection +--- GSL::Monte::Miser::Params#min_calls_per_bisection= + The minimum number of function calls required to proceed with a bisection step. When a recursive step has fewer calls available than ((|min_calls_per_bisection|)) it performs a plain Monte Carlo estimate of the current sub-region and terminates its branch of the recursion. The default value of this parameter is 32 * min_calls. +--- GSL::Monte::Miser::Params#alpha +--- GSL::Monte::Miser::Params#alpha= + This parameter controls how the estimated variances for the two sub-regions of a bisection are combined when allocating points. With recursive sampling the overall variance should scale better than 1/N, since the values from the sub-regions will be obtained using a procedure which explicitly minimizes their variance. To accommodate this behavior the MISER algorithm allows the total variance to depend on a scaling parameter ((|alpha|)), The authors of the original paper describing MISER recommend the value ((|alpha|)) = 2 as a good choice, obtained from numerical experiments, and this is used as the default value in this implementation. +--- GSL::Monte::Miser::Params#dither +--- GSL::Monte::Miser::Params#dither= +This parameter introduces a random fractional variation of size dither into each bisection, which can be used to break the symmetry of integrands which are concentrated near the exact center of the hypercubic integration region. The default value of dither is zero, so no variation is introduced. If needed, a typical value of dither is 0.1. + +== Vegas Parameters (GSL-1.13 or later) +--- GSL::Monte::Vegas#params_get + Returns the parameters of the integrator state as an instance of (({GSL::Monte::Vegas::Params})) class. +--- GSL::Monte::Vegas#params_set(params) + Sets the integrator parameters based on values provided in an object of the (({GSL::Monte::Vegas::Params})) class ((|params|)). + +=== Accessors of (({GSL::Monte::Vegas::Params})) +--- GSL::Monte::Vegas::Params#alpha +--- GSL::Monte::Vegas::Params#alpha= + Controls the stiffness of the rebinning algorithm. It is typically set between one and two. A value of zero prevents rebinning of the grid. The default value is 1.5. +--- GSL::Monte::Vegas::Params#iterations +--- GSL::Monte::Vegas::Params#iterations= + The number of iterations to perform for each call to the routine. The default value is 5 iterations. +--- GSL::Monte::Vegas::Params#stage +--- GSL::Monte::Vegas::Params#stage= + Setting this determines the stage of the calculation. Normally, stage = 0 which begins with a new uniform grid and empty weighted average. Calling vegas with stage = 1 retains the grid from the previous run but discards the weighted average, so that one can "tune" the grid using a relatively small number of points and then do a large run with stage = 1 on the optimized grid. Setting stage = 2 keeps the grid and the weighted average from the previous run, but may increase (or decrease) the number of histogram bins in the grid depending on the number of calls available. Choosing stage = 3 enters at the main loop, so that nothing is changed, and is equivalent to performing additional iterations in a previous call. +--- GSL::Monte::Vegas::Params#mode +--- GSL::Monte::Vegas::Params#mode= + The possible choices are (({GSL::VEGAS::MODE_IMPORTANCE})), (({GSL::VEGAS::MODE_STRATIFIED})), (({GSL::VEGAS::MODE_IMPORTANCE_ONLY})). This determines whether VEGAS will use importance sampling or stratified sampling, or whether it can pick on its own. In low dimensions VEGAS uses strict stratified sampling (more precisely, stratified sampling is chosen if there are fewer than 2 bins per box). +--- GSL::Monte::Vegas::Params#verbose +--- GSL::Monte::Vegas::Params#verbose= + Set the level of information printed by VEGAS. All information is written to the stream ostream. The default setting of verbose is -1, which turns off all output. A verbose value of 0 prints summary information about the weighted average and final result, while a value of 1 also displays the grid coordinates. A value of 2 prints information from the rebinning procedure for each iteration. + == Example #!/usr/bin/env ruby rd/monte.rd @@ -71,7 +71,7 @@ the function in question. --- GSL::MultiMin::FdfMinimizer.alloc(type, n) --- GSL::MultiMin::FMinimizer.alloc(type, n) These method create a minimizer of type ((|type|)) for an ((|n|))-dimension function. - The type is given by a name of the algorithm, or by a Ruby constant. + The type is given by a string, or by a Ruby constant. * (({GSL::MultiMin::FdfMinimizer::CONJUGATE_FR})) or (({"conjugate_fr"})) * (({GSL::MultiMin::FdfMinimizer::CONJUGATE_PR})) or (({"conjugate_pr"})) @@ -79,6 +79,7 @@ the function in question. * (({GSL::MultiMin::FdfMinimizer::VECTOR_BFGS2})) or (({"vector_bfgs2"})) (GSL-1.9 or later) * (({GSL::MultiMin::FdfMinimizer::STEEPEST_DESCENT})) or (({"steepest_descent"})) * (({GSL::MultiMin::FMinimizer::NMSIMPLEX})) or (({"nmsimplex"})) + * (({GSL::MultiMin::FMinimizer::NMSIMPLEX2RAND})) or (({"nmsimplex2rand"})) (GSL-1.13) * ex: include GSL::MultiMin rd/multimin.rd @@ -42,6 +42,38 @@ Contents: => NArray.float(3): [ 6.0, 17.0, 34.0 ] +--- GSL::Poly.eval_derivs(c, x) +--- GSL::Poly.eval_derivs(c, x, lenres) + (GSL-1.13) Evaluate and return a polynomial and its derivatives. The output contains the values of d^k P/d x^k for the specified value of x starting with k = 0. The input polynomial ((|c|)) can be an (({Array})), (({GSL::Poly})) or an (({NArray})). If ((|lenres|)) is not given, ((|lenres = LENGTH(c) + 1|)) is used, therefore the last element of the output is 0. + +--- GSL::Poly#eval_derivs(x) +--- GSL::Poly#eval_derivs(x, lenres) + (GSL-1.13) Evaluate and return a polynomial and its derivatives. The output contains the values of d^k P/d x^k for the specified value of x starting with k = 0. If ((|lenres|)) is not given, ((|lenres = LENGTH(self) + 1|)) is used, therefore the last element of the output is 0. + + Ex.) + >> ary = [1, 2, 3] + => [1, 2, 3] + >> GSL::Poly.eval_derivs(ary, 1) + => [6.0, 8.0, 6.0, 0.0] + >> na = NArray[1.0, 2, 3] + => NArray.float(3): + [ 1.0, 2.0, 3.0 ] + >> GSL::Poly.eval_derivs(na, 1) + => NArray.float(4): + [ 6.0, 8.0, 6.0, 0.0 ] + >> poly = GSL::Poly[1.0, 2, 3] + => GSL::Poly + [ 1.000e+00 2.000e+00 3.000e+00 ] + >> GSL::Poly.eval_derivs(poly, 1) + => GSL::Poly + [ 6.000e+00 8.000e+00 6.000e+00 0.000e+00 ] + >> poly.eval_derivs(1) + => GSL::Poly + [ 6.000e+00 8.000e+00 6.000e+00 0.000e+00 ] + >> poly.eval_derivs(1, 3) + => GSL::Poly + [ 6.000e+00 8.000e+00 6.000e+00 ] + == Solving polynomial equations === Quadratic Equations --- GSL::Poly::solve_quadratic(a, b, c) rd/poly.rd @@ -85,7 +85,7 @@ Func3 = GSL::Function.alloc { |x| } Func4 = GSL::Function.alloc { |x| x - 30.0/(1.0 + 1e5*pow(x - 0.8, 2.0)) } -types = ["goldensection", "brent"] +types = ["goldensection", "brent", "quad_golden"] types.each do |t| test_f(t, "cos(x) [0 (3) 6]", F_cos, 0.0, 3.0, 6.0, M_PI) test_f(t, "x^4 - 1 [-3 (-1) 17]", Func1, -3.0, -1.0, 17.0, 0.0); tests/min.rb @@ -36,6 +36,17 @@ def test_f(desc, f, initpt) status |= s.fval.abs > 1e-5 ? 1 : 0 GSL::Test::test(status, "#{s.name}, on #{desc}: #{iter} iterations, f(x)=#{s.fval}") + + s = GSL::MultiMin::FMinimizer.alloc("nmsimplex2rand", f.n) + s.set(f, x, step_size) + iter = 0 + begin + status = s.iterate + status = GSL::MultiMin.test_size(s.size, 1e-3) + end while iter < 5000 and status == GSL::CONTINUE + + status |= s.fval.abs > 1e-5 ? 1 : 0 + GSL::Test::test(status, "#{s.name}, on #{desc}: #{iter} iterations, f(x)=#{s.fval}") end def roth_initpt tests/multimin.rb f056bdce4e7e1add9586b7aa22c13bb398a4968a ytsunesada ytsunesada@6e764f74-f39f-46f8-8c54-8307d62afe8d http://github.com/romanbsd/rb-gsl/commit/fcf16906687c03a8e31b376d583b6a65879d941d fcf16906687c03a8e31b376d583b6a65879d941d

2009-08-08T21:21:11-07:00

Added gsl-1.12.90 features git-svn-id: http://rb-gsl.rubyforge.org/svn/trunk/rb-gsl@96 6e764f74-f39f-46f8-8c54-8307d62afe8d 3873080fa84071b3848a664a8e8d4d6549ca3767 ytsunesada ytsunesada@6e764f74-f39f-46f8-8c54-8307d62afe8d