Ask Your Question
0

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

asked 2011-01-19 09:45:40 -0500

Shashank gravatar image

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 flag offensive close merge delete

2 answers

Sort by ยป oldest newest most voted
1

answered 2011-01-19 12:55:21 -0500

DSM gravatar image

updated 2011-01-19 12:58:09 -0500

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

edit flag offensive delete link more
0

answered 2011-01-19 12:35:39 -0500

Shashank gravatar image

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
edit flag offensive delete link more

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

Stats

Asked: 2011-01-19 09:45:40 -0500

Seen: 378 times

Last updated: Jan 19 '11