# Is is possible to call the GSL vegas monte carlo routine from sage?

I am trying to perform a monte-carlo integration of a function defined in sage. It is possible to call the vegas monte-carlo routine in GSL library from sage? Or is there a vegas routine in some other package like numpy accessible from sage?

edit retag close merge delete

Sort by ยป oldest newest most voted

If I were planning on using this regularly I'd probably write a wrapper function to turn a sage expression into a pygsl-compatible function automatically (N.B. I installed pygsl to make life a little easier):


## warning: untested!

import numpy as np
import pygsl.monte as monte
import pygsl.rng

def sage_expr_wrap(some_expr,varorder=None):
if varorder == None: varorder=sorted(some_expr.variables())
def numpy_fn(xx, params):
return np.float(some_expr.subs(dict(zip(varorder, xx))))
return monte.gsl_monte_function(numpy_fn, None, len(varorder))

def sage_vegas(some_expr, xl, xu, ncalls):
xl = np.array(xl,dtype=float)
xu = np.array(xu,dtype=float)
ncalls = int(ncalls)
grand = pygsl.rng.mt19937_1999()
G = sage_expr_wrap(some_expr) # warning: defaults to alphabetical order
s = monte.vegas(len(some_expr.variables()))
s.init()
res, err = s.integrate(G, xl, xu, ncalls, grand)
return res, err

sage: var("x y z")
(x, y, z)
sage: print sage_vegas(3*(x+y+z),[0,0,0],[1,1,1],100)
(4.5695212037401234, 0.0432533513425419)
sage: print sage_vegas(3*(x+y+z),[0,0,0],[1,1,1],1000)
(4.5050689799289199, 0.0051246768259351639)
sage: print sage_vegas(3*(x+y+z),[0,0,0],[1,1,1],10000)
(4.5002245924643995, 0.00031658038710181448)


but if speed's a concern obviously the cythonic version will win. :^)

more

I am sorry for asking a question that is already present in the examples file. There is an example explaining the implementation of Vegas monte-carlo in sage in the folder

\$SAGE_ROOT/examples/gsl/examples/monte.pyx

more