OpenXM/src/ox_gsl/call_gsl.c - diff

Return to call_gsl.c CVS log

Up to [local] / OpenXM / src / ox_gsl

Diff for /OpenXM/src/ox_gsl/call_gsl.c between version 1.4 and 1.7

-version 1.4, 2018/04/18 02:20:51
+version 1.7, 2018/06/07 01:53:33
 Line 1
 Line 1
 Line 1
- /* $OpenXM: OpenXM/src/ox_gsl/call_gsl.c,v 1.3 2018/04/17 00:56:38 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>
 Line 7
 Line 7
 Line 7
  #include <gsl/gsl_errno.h>
  #include <gsl/gsl_sf_gamma.h>
  #include <gsl/gsl_integration.h>
+ #include <gsl/gsl_monte.h>
+ #include <gsl/gsl_monte_plain.h>
+ #include <gsl/gsl_monte_miser.h>
+ #include <gsl/gsl_monte_vegas.h>
  #include "ox_gsl.h"
  extern int Debug;
  // local prototype declarations
-Line 39  void  call_gsl_sf_lngamma_complex_e() {
+Line 43  void  call_gsl_sf_lngamma_complex_e() {
 Line 39  void  call_gsl_sf_lngamma_complex_e() {
 Line 43  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");
-Line 94  void  call_gsl_integration_qags() {
+Line 99  void  call_gsl_integration_qags() {
 Line 94  void  call_gsl_integration_qags() {
 Line 99  void  call_gsl_integration_qags() {
    ans = (cmo *)new_cmo_list_array((void *)r,3);
    push(ans);
  }
+ 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_f_n) ox_printf("f_n, dim=%d\n",dim);
+   for (i=0; i<dim; i++) {
+     if (debug_f_n) ox_printf("x%d=%lg, ",i,x[i]);
+     register_entry(xx[i],x[i]);
+   }
+   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_f_n) ox_printf("\nf_x(...) -> d=%lg\n",d);
+   return(d);
+ }
+ void  call_gsl_monte_plain_miser_vegas_integrate(int type) {
+   int argc;
+   int dim;
+   int dim2;
+   cmo *cr[2];
+   cmo *ans;
+   double *xl;
+   double *xu;
+   gsl_monte_function G = {&f_n,0,0};
+   double res, err;
+   const gsl_rng_type *T;
+   gsl_rng *r;
+   //  size_t calls = 500000;
+   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);
+   //  gsl_set_error_handler_off();
+   gsl_set_error_handler((gsl_error_handler_t *)myhandler);
+   argc = get_i(); // number of args
+   if (argc < 3) {
+     pops(argc);
+     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;
+   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);
+   ans = (cmo *)new_cmo_list_array((void *)cr,2);
+   push(ans);
+ }

FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>