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

edit

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. :^)