===================================================================
RCS file: /home/cvs/OpenXM/src/ox_gsl/call_gsl.c,v
retrieving revision 1.6
retrieving revision 1.7
diff -u -p -r1.6 -r1.7
--- OpenXM/src/ox_gsl/call_gsl.c	2018/06/06 10:46:00	1.6
+++ OpenXM/src/ox_gsl/call_gsl.c	2018/06/07 01:53:33	1.7
@@ -1,4 +1,4 @@
-/* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.5 2018/06/06 07:40:32 takayama Exp $ 
+/* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.6 2018/06/06 10:46:00 takayama Exp $ 
 */
 //#include <gsl/gsl_types.h>
 //#include <gsl/gsl_sys.h>
@@ -43,7 +43,8 @@ void  call_gsl_sf_lngamma_complex_e() {
   push(ans);
 }
 
-cmo *Func_x=NULL;
+cmo *Func_x=NULL;  // function to evalute.
+cmo *Vlist=NULL;   // list of variables (not implemented yet).
 double f_x(double x,void *params) {
   double d;
   if (Debug) ox_printf("f_x\n");
@@ -100,25 +101,27 @@ void  call_gsl_integration_qags() {
 }  
 
 double f_n(double x[],size_t dim,void *params) {
+  int debug_f_n=0;
   double d;
   int i;
   char *xx[30]={"x0","x1","x2","x3","x4","x5","x6","x7","x8","x9",
                 "x10","x11","x12","x13","x14","x15","x16","x17","x18","x19",
                 "x20","x21","x22","x23","x24","x25","x26","x27","x28","x29"
                };
+//  debug_f_n = Debug;
   if (dim > 30) GSL_ERROR("f_n supports functions with args <= 30",GSL_ETOL);
   init_dic();
-  if (Debug) ox_printf("f_n, dim=%d\n",dim);
+  if (debug_f_n) ox_printf("f_n, dim=%d\n",dim);
   for (i=0; i<dim; i++) {
-    if (Debug) ox_printf("x%d=%lg, ",i,x[i]);
+    if (debug_f_n) ox_printf("x%d=%lg, ",i,x[i]);
     register_entry(xx[i],x[i]);
   }
-  if (Debug) { ox_printf("\n"); fflush(NULL); }
+  if (debug_f_n) { ox_printf("\n"); fflush(NULL); }
   if (eval_cmo(Func_x,&d)==0) GSL_ERROR("eval_cmo fails in f_n",GSL_ETOL);
-  if (Debug) ox_printf("\nf_x(...) -> d=%lg\n",d);
+  if (debug_f_n) ox_printf("\nf_x(...) -> d=%lg\n",d);
   return(d);
 }
-void  call_gsl_monte_plain_integrate() {
+void  call_gsl_monte_plain_miser_vegas_integrate(int type) {
   int argc;
   int dim;
   int dim2;
@@ -132,7 +135,10 @@ void  call_gsl_monte_plain_integrate() {
   const gsl_rng_type *T;
   gsl_rng *r;
   //  size_t calls = 500000;
-  size_t calls = 500; // for test
+  size_t calls = 100000; // for test  1000
+  gsl_monte_plain_state *s_p;
+  gsl_monte_miser_state *s_m;
+  gsl_monte_vegas_state *s_v;
   gsl_rng_env_setup ();
   T = gsl_rng_default;
   r = gsl_rng_alloc (T);
@@ -140,24 +146,45 @@ void  call_gsl_monte_plain_integrate() {
   //  gsl_set_error_handler_off();
   gsl_set_error_handler((gsl_error_handler_t *)myhandler);
   argc = get_i(); // number of args
-  if (argc != 3) {
+  if (argc < 3) {
     pops(argc);
-    push(make_error2("The argc must be 3 for gsl_integration_qags.",NULL,0,-1));
+    push(make_error2("The argc must be >=3 for gsl_integration_qags.",NULL,0,-1));
     return;
   }
   Func_x = pop();
   xl = get_double_list(&dim);
   xu = get_double_list(&dim2);
+  if (argc >= 4) {
+    calls = (int) get_double();
+  }
+  if (argc >= 5) {
+    Vlist = pop();  // list of variables. Not implemented yet.
+  }
   if (dim != dim2) {
     push(make_error2("gsl_monte_plain: dim of interval differs.",NULL,0,-1));
     return;
   }
   G.dim=dim;
-  gsl_monte_plain_state *s = gsl_monte_plain_alloc (dim);
-  gsl_monte_plain_integrate (&G, xl, xu, dim, calls, r, s, 
-                             &res, &err);
-  gsl_monte_plain_free (s);
-
+  switch(type) {
+  case 0: // plain
+    s_p = gsl_monte_plain_alloc (dim);
+    gsl_monte_plain_integrate (&G, xl, xu, dim, calls, r, s_p, 
+                               &res, &err);
+    gsl_monte_plain_free (s_p);
+    break;
+  case 1: // miser
+    s_m = gsl_monte_miser_alloc (dim);
+    gsl_monte_miser_integrate (&G, xl, xu, dim, calls, r, s_m, 
+                               &res, &err);
+    gsl_monte_miser_free (s_m);
+    break;
+  default: // vegas;
+    s_v = gsl_monte_vegas_alloc (dim);
+    gsl_monte_vegas_integrate (&G, xl, xu, dim, calls, r, s_v, 
+                               &res, &err);
+    gsl_monte_vegas_free (s_v);
+    break;
+  }
   if (Debug) ox_printf("result = %lg\n",res);
   cr[0] = (cmo *)new_cmo_double(res);
   cr[1] = (cmo *)new_cmo_double(err);