Monte Carlo API problem

Monte Carlo API problem

[Szymon Jaroszewicz]

Hi,

I was wrapping monte carlo routines in Python in ctypesGSL and
hit an API problem esp. in Vegas.

The gsl_monte_vegas_state structure has several public fields
which can be changed to modify the algorithm.  However to make
these fields available in Python through ctypes I would need to
wrap several private members of the structure - quite unelegant.

The situation is somewhat better in Miser where all public fields
are at the beginning of the structure and private fields can
simply be skipped.

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.

-- Szymon

Re: Monte Carlo API problem

[Brian Gough]

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.

Re: Monte Carlo API problem

[Brian Gough]

At Mon, 20 Jul 2009 10:10:41 -0400 (EDT),
Szymon Jaroszewicz wrote:
> 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.

Thanks for the patch.  That looks fine. Sorry for the delay in replying.

Minor point, using 
gsl_monte_vegas_params_get
instead of
gsl_monte_vegas_get_params
would be closer to the other _get,_set functions in the library.

It is nice that it only adds the the API and is backwards compatible,
so I can put it in the 1.13 release which I am currently preparing.
I've fixed all the outstanding bugs I wanted to clear so it is almost
ready to go.

-- 
Brian Gough

Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/

Re: Monte Carlo API problem

[Szymon Jaroszewicz]

[-- 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[])
 {