Re: Monte Carlo API problem

[Szymon Jaroszewicz <sj@cs.umb.edu>]

[-- Attachment #1: Type: TEXT/PLAIN, Size: 837 bytes --]

Hi Brian,

Attached is the patch.  As there are many parameters for miser and vegas, 
I added a structure which holds all the parameters and added functions 
which copy the values in and out of the structure.  Otherwise we would 
need over 10 getter/setter functions  for each method.

Let me know if you like this approach, I'll update the manual then.

All the best,

Szymon

On Thu, 2009-07-09 at 21:06 +0100, Brian Gough wrote:
At Wed, 8 Jul 2009 06:54:17 -0400 (EDT),
> Szymon Jaroszewicz wrote:
> > Should the API be updated to allow setting at least some of the
> > parameters without messing with state fields directly?  At the
> > minimum I think there should be a way to get the chisq field for
> > Vegas which is used to test convergence.
>
> This makes sense, would you like to send me a patch for the functions
> you want.
>

[-- Attachment #2: Type: TEXT/PLAIN, Size: 8438 bytes --]

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 <stdlib.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_monte.h>
 
@@ -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, &params) ; params.alpha = 1.5 ; gsl_monte_miser_set_params(s, &params) ; 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, &params) ; params.alpha = 2 ; params.iterations = 3 ; gsl_monte_vegas_set_params(s, &params) ; 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[])
 {