diff --git a/monte/demo.c b/monte/demo.c index ed41151..6b58e6b 100644 --- a/monte/demo.c +++ b/monte/demo.c @@ -81,9 +81,9 @@ main () { gsl_monte_vegas_integrate (&G, xl, xu, 3, calls, r, s, &res, &err); printf ("result = % .6f sigma = % .6f chisq/dof = %.1f\n", - res, err, s->chisq); + res, err, gsl_monte_vegas_chisq (s)); } - while (fabs (s->chisq - 1.0) > 0.5); + while (fabs (gsl_monte_vegas_chisq (s) - 1.0) > 0.5); display_results ("vegas final", res, err); diff --git a/monte/gsl_monte_miser.h b/monte/gsl_monte_miser.h index b053063..76dcb3f 100644 --- a/monte/gsl_monte_miser.h +++ b/monte/gsl_monte_miser.h @@ -77,6 +77,17 @@ int gsl_monte_miser_init(gsl_monte_miser_state* state); void gsl_monte_miser_free(gsl_monte_miser_state* state); +typedef struct { + double estimate_frac; + size_t min_calls; + size_t min_calls_per_bisection; + double alpha; + double dither; +} gsl_monte_miser_params; + +void gsl_monte_miser_get_params(const gsl_monte_miser_state* state, gsl_monte_miser_params* params); + +void gsl_monte_miser_set_params(gsl_monte_miser_state* state, const gsl_monte_miser_params* params); __END_DECLS diff --git a/monte/gsl_monte_vegas.h b/monte/gsl_monte_vegas.h index 8ecef46..35b8090 100644 --- a/monte/gsl_monte_vegas.h +++ b/monte/gsl_monte_vegas.h @@ -22,6 +22,7 @@ #ifndef __GSL_MONTE_VEGAS_H__ #define __GSL_MONTE_VEGAS_H__ +#include #include #include @@ -99,6 +100,21 @@ int gsl_monte_vegas_init(gsl_monte_vegas_state* state); void gsl_monte_vegas_free (gsl_monte_vegas_state* state); +double gsl_monte_vegas_chisq (const gsl_monte_vegas_state* state); + +typedef struct { + double alpha; + size_t iterations; + int stage; + int mode; + int verbose; + FILE * ostream; +} gsl_monte_vegas_params; + +void gsl_monte_vegas_get_params(const gsl_monte_vegas_state* state, gsl_monte_vegas_params* params); + +void gsl_monte_vegas_set_params(gsl_monte_vegas_state* state, const gsl_monte_vegas_params* params); + __END_DECLS #endif /* __GSL_MONTE_VEGAS_H__ */ diff --git a/monte/miser.c b/monte/miser.c index f7b29c3..c468a64 100644 --- a/monte/miser.c +++ b/monte/miser.c @@ -567,6 +567,27 @@ gsl_monte_miser_free (gsl_monte_miser_state * s) free (s); } +void +gsl_monte_miser_get_params(const gsl_monte_miser_state* s, gsl_monte_miser_params* p) +{ + p->estimate_frac = s->estimate_frac; + p->min_calls = s->min_calls; + p->min_calls_per_bisection = s->min_calls_per_bisection; + p->alpha = s->alpha; + p->dither = s->dither; +} + +void +gsl_monte_miser_set_params(gsl_monte_miser_state* s, const gsl_monte_miser_params* p) +{ + s->estimate_frac = p->estimate_frac; + s->min_calls = p->min_calls; + s->min_calls_per_bisection = p->min_calls_per_bisection; + s->alpha = p->alpha; + s->dither = p->dither; +} + + static int estimate_corrmc (gsl_monte_function * f, const double xl[], const double xu[], diff --git a/monte/test.c b/monte/test.c index c1dd0a7..72d2810 100644 --- a/monte/test.c +++ b/monte/test.c @@ -282,6 +282,26 @@ main (void) #undef MONTE_SPEEDUP #endif +#ifdef MISER +#define NAME "miser(params)" +#define MONTE_STATE gsl_monte_miser_state +#define MONTE_ALLOC gsl_monte_miser_alloc +#define MONTE_PARAMS gsl_monte_miser_params +#define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_miser_get_params(s, ¶ms) ; params.alpha = 1.5 ; gsl_monte_miser_set_params(s, ¶ms) ; gsl_monte_miser_integrate(f,xl,xu,dim,calls,r,s,res,err); } +#define MONTE_FREE gsl_monte_miser_free +#define MONTE_SPEEDUP 2 +#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 5.0 * expected, NAME ", %s, abserr[%d] (obs %g vs plain %g)", I->description, i, err, expected) +#include "test_main.c" +#undef NAME +#undef MONTE_STATE +#undef MONTE_ALLOC +#undef MONTE_PARAMS +#undef MONTE_INTEGRATE +#undef MONTE_FREE +#undef MONTE_ERROR_TEST +#undef MONTE_SPEEDUP +#endif + #ifdef VEGAS #define NAME "vegas" #define MONTE_STATE gsl_monte_vegas_state @@ -289,7 +309,7 @@ main (void) #define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err) ; } #define MONTE_FREE gsl_monte_vegas_free #define MONTE_SPEEDUP 3 -#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 3.0 * (expected == 0 ? 1.0/(I->calls/MONTE_SPEEDUP) : expected), NAME ", %s, abserr[%d] (obs %g vs exp %g)", I->description, i, err, expected) ; gsl_test(s->chisq < 0, NAME " returns valid chisq (%g)", s->chisq) +#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 3.0 * (expected == 0 ? 1.0/(I->calls/MONTE_SPEEDUP) : expected), NAME ", %s, abserr[%d] (obs %g vs exp %g)", I->description, i, err, expected) ; gsl_test(gsl_monte_vegas_chisq(s) < 0, NAME " returns valid chisq (%g)", gsl_monte_vegas_chisq(s)) #include "test_main.c" #undef NAME #undef MONTE_STATE @@ -308,11 +328,32 @@ main (void) #define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err) ; gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err); } #define MONTE_FREE gsl_monte_vegas_free #define MONTE_SPEEDUP 3 -#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 3.0 * (expected == 0 ? 1.0/(I->calls/MONTE_SPEEDUP) : expected), NAME ", %s, abserr[%d] (obs %g vs exp %g)", I->description, i, err, expected); gsl_test(s->chisq < 0, NAME " returns valid chisq (%g)", s->chisq) +#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 3.0 * (expected == 0 ? 1.0/(I->calls/MONTE_SPEEDUP) : expected), NAME ", %s, abserr[%d] (obs %g vs exp %g)", I->description, i, err, expected); gsl_test(gsl_monte_vegas_chisq(s) < 0, NAME " returns valid chisq (%g)", gsl_monte_vegas_chisq(s)) +#include "test_main.c" +#undef NAME +#undef MONTE_STATE +#undef MONTE_ALLOC +#undef MONTE_INTEGRATE +#undef MONTE_FREE +#undef MONTE_ERROR_TEST +#undef MONTE_SPEEDUP +#endif + + +#ifdef VEGAS +#define NAME "vegas(params)" +#define MONTE_STATE gsl_monte_vegas_state +#define MONTE_ALLOC gsl_monte_vegas_alloc +#define MONTE_PARAMS gsl_monte_vegas_params +#define MONTE_INTEGRATE(f,xl,xu,dim,calls,r,s,res,err) { gsl_monte_vegas_get_params(s, ¶ms) ; params.alpha = 2 ; params.iterations = 3 ; gsl_monte_vegas_set_params(s, ¶ms) ; gsl_monte_vegas_integrate(f,xl,xu,dim,calls,r,s,res,err); } +#define MONTE_FREE gsl_monte_vegas_free +#define MONTE_SPEEDUP 3 +#define MONTE_ERROR_TEST(err,expected) gsl_test(err > 3.0 * (expected == 0 ? 1.0/(I->calls/MONTE_SPEEDUP) : expected), NAME ", %s, abserr[%d] (obs %g vs exp %g)", I->description, i, err, expected); gsl_test(gsl_monte_vegas_chisq(s) < 0, NAME " returns valid chisq (%g)", gsl_monte_vegas_chisq(s)) #include "test_main.c" #undef NAME #undef MONTE_STATE #undef MONTE_ALLOC +#undef MONTE_PARAMS #undef MONTE_INTEGRATE #undef MONTE_FREE #undef MONTE_ERROR_TEST diff --git a/monte/test_main.c b/monte/test_main.c index a797e8a..7cd71d1 100644 --- a/monte/test_main.c +++ b/monte/test_main.c @@ -15,7 +15,10 @@ for (I = problems ; I->f != 0; I++) for (i = 0; i < TRIALS ; i++) { MONTE_STATE *s = MONTE_ALLOC (I->dim); - +#ifdef MONTE_PARAMS + MONTE_PARAMS params; +#endif + I->f->dim = I->dim; MONTE_INTEGRATE (I->f, I->xl, I->xu, diff --git a/monte/vegas.c b/monte/vegas.c index e8c4d38..87d91fa 100644 --- a/monte/vegas.c +++ b/monte/vegas.c @@ -520,6 +520,32 @@ gsl_monte_vegas_free (gsl_monte_vegas_state * s) free (s); } +double +gsl_monte_vegas_chisq (const gsl_monte_vegas_state* s) +{ + return s->chisq; +} + +void gsl_monte_vegas_get_params(const gsl_monte_vegas_state* s, gsl_monte_vegas_params* p) +{ + p->alpha = s->alpha; + p->iterations = s->iterations; + p->stage = s->stage; + p->mode = s->mode; + p->verbose = s->verbose; + p->ostream = s->ostream; +} + +void gsl_monte_vegas_set_params(gsl_monte_vegas_state* s, const gsl_monte_vegas_params* p) +{ + s->alpha = p->alpha; + s->iterations = p->iterations; + s->stage = p->stage; + s->mode = p->mode; + s->verbose = p->verbose; + s->ostream = p->ostream; +} + static void init_box_coord (gsl_monte_vegas_state * s, coord box[]) {