Revision history [back]

Although some gsl functions exist in Sage, it doesn't appear that the monte module has been wrapped.

So I downloaded pygsl from http://sourceforge.net/projects/pygsl/ and installed it using "sage -python setup.py build --gsl-prefix=/opt/local/; sage -python setup.py install". Note I had to specify the prefix expliclty (which is a little weird, because gsl-config was on my path, but whatever.)

After that I could use the monte functions (modified from the pygsl example):

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

def g(xx, params):
    v = 3*sum(xx)
    return v

G = monte.gsl_monte_function(g, None, 3)

xl = np.array([0,0,0])
xu = np.array([1,1,1],dtype=float)

grand = pygsl.rng.mt19937_1999()
calls = int(10**6)

var("x y z")
exact_int = integrate(integrate(integrate(3*(x+y+z), x),y),z)
exact = exact_int.subs(dict(zip([x,y,z], map(float,xu)))) - exact_int.subs(dict(zip([x,y,z], map(float, xl))))

def display_results (title, result, error):
    print "%s ==================" %  title
    print "result = % .6f" % result
    print "sigma  = % .6f" % error
    print "exact  = % .6f" % exact
    t = (result - exact, np.abs(result - exact) / error)
    print "error  = % .6f = %.1g sigma" % t

def run_vegas():
    s = monte.vegas(int(3))
    s.init()

    res, err = s.integrate(G, xl, xu, int(10000), grand)
    display_results ("vegas warm-up", res, err)

    print "converging..."
    while 1:
        res, err = s.integrate(G, xl, xu, int(calls/1000), grand)
        chisq = s.get_chisq()
        print "result = % .6f sigma = % .6f chisq/dof = %.1f" % (res, err, chisq);
        display_results("vegas:", res, err)
        if (abs(chisq - 1.0) <= 0.5):    
            break

run_vegas()

produced

sage: load "vegasexample.sage"
vegas warm-up ==================
result =  4.500225
sigma  =  0.000317
exact  =  4.500000
error  =  0.000225 = 0.7 sigma
converging...
result =  4.497538 sigma =  0.004253 chisq/dof = 1.0
vegas: ==================
result =  4.497538
sigma  =  0.004253
exact  =  4.500000
error  = -0.002462 = 0.6 sigma

You should note that I had to put in some "int()" calls there, which is because the preprocessor was turning integer literals into sage Integers, which was causing GSL problems.

While the above works, 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:

## 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)

Although some gsl functions exist in Sage, it doesn't appear that the monte module has been wrapped.

So I downloaded pygsl from http://sourceforge.net/projects/pygsl/ and installed it using "sage -python setup.py build --gsl-prefix=/opt/local/; sage -python setup.py install". Note I had to specify the prefix expliclty (which is a little weird, because gsl-config was on my path, but whatever.)

After that I could use the monte functions (modified from the pygsl example):

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

def g(xx, params):
    v = 3*sum(xx)
    return v

G = monte.gsl_monte_function(g, None, 3)

xl = np.array([0,0,0])
xu = np.array([1,1,1],dtype=float)

grand = pygsl.rng.mt19937_1999()
calls = int(10**6)

var("x y z")
exact_int = integrate(integrate(integrate(3*(x+y+z), x),y),z)
exact = exact_int.subs(dict(zip([x,y,z], map(float,xu)))) - exact_int.subs(dict(zip([x,y,z], map(float, xl))))

def display_results (title, result, error):
    print "%s ==================" %  title
    print "result = % .6f" % result
    print "sigma  = % .6f" % error
    print "exact  = % .6f" % exact
    t = (result - exact, np.abs(result - exact) / error)
    print "error  = % .6f = %.1g sigma" % t

def run_vegas():
    s = monte.vegas(int(3))
    s.init()

    res, err = s.integrate(G, xl, xu, int(10000), grand)
    display_results ("vegas warm-up", res, err)

    print "converging..."
    while 1:
        res, err = s.integrate(G, xl, xu, int(calls/1000), grand)
        chisq = s.get_chisq()
        print "result = % .6f sigma = % .6f chisq/dof = %.1f" % (res, err, chisq);
        display_results("vegas:", res, err)
        if (abs(chisq - 1.0) <= 0.5):    
            break

run_vegas()

produced

sage: load "vegasexample.sage"
vegas warm-up ==================
result =  4.500225
sigma  =  0.000317
exact  =  4.500000
error  =  0.000225 = 0.7 sigma
converging...
result =  4.497538 sigma =  0.004253 chisq/dof = 1.0
vegas: ==================
result =  4.497538
sigma  =  0.004253
exact  =  4.500000
error  = -0.002462 = 0.6 sigma

You should note that I had to put in some "int()" calls there, which is because the preprocessor was turning integer literals into sage Integers, which was causing GSL problems.

While the above works, if 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: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)

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)