[SciPy-Dev] Request for extension to scipy.integrate
Nathan Woods
charlesnwoods@gmail....
Thu Apr 25 19:53:24 CDT 2013
SciPy's multidimensional integration capabilities are somewhat limited, as mentioned previously: https://github.com/scipy/scipy/issues/2098. Although dblquad and tplquad attempt to address this problem, they do not allow any detailed control over the underlying quad algorithm and there is no option for 4-dimensional integration at all.
One obvious instance where this functionality would be desired is in the evaluation of volume integrals in 4-dimensional space-time. Another instance is the integration of discontinuous functions, where access to the "points" option of quad is critical. As it stands, SciPy does not provide any good functionality for either of these situations.
I've written a recursive wrapper function for quad that resolves both of these problems, allowing n-dimensional integration with complete specification of options for quad. I've included a test problem that demonstrates the resolution of both of the above problems (just execute the script). I've done my best to document everything, and I've verified the result of the test problem against Mathematica.
I would like this code (or equivalent functionality) to be included in the integrate package. I'm open to any feedback or suggestions, particularly with the interface.
Nathan Woods
"""
Recursive integration using SciPy's quad function.
Contains the following:
class Error - module wrapper for Exception. For future expansion only.
function mul_quad - Evaluate multiple integrals using recursive calls to quad.
"""
from scipy.integrate import quad
import numpy
class Error(Exception):
pass
def mul_quad(func,ranges,args=(),opts=(),_depth=0,_int_vars=()):
"""
Evaluate multiple integrals through recursive calls to scipy.integrate.quad.
mul_quad takes care of the programming required to perform nested
integration using the 1-dimensional integration routines provided
by SciPy's adaptive quadrature function. It extends the capabilities
of dblquad and tplquad by allowing for more levels of integration,
and allowing the user to specify nearly the full range of options
allowed by quad, for each level of integration. Users are cautioned
that nested Gaussian integration of this kind is computationally
intensive, and may be unsuitable for many nested integrals.
Usage: mul_quad(func,ranges,args=(),opts=(),_depth=0,_int_vars=())
Inputs:
func - callable, acceptable to SciPy's quad, returning a number.
Should accept a float, followed by the contents of _int_vars and
args, e.g. if x is a float, args=(1.4,string), and _int_vars =
(1,1,3), then func should be of the form
func(x,1,1,3,1.4,string).
ranges - sequence describing the ranges of integration. Integrals
are performed in order, so ranges[0] corresponds to the first
argument of func, ranges[1] to the second, and so on. Each
element of ranges may be either a constant sequence of length 2
or else a function that returns such a sequence. If a function,
then it will be called with all of the integration arguments
available to that point. e.g. for func = f(x0,x1,x2,x3), the
range of integration for x0 may be defined as either a constant
such as (0,1) or as a function range0(x1,x2,x3). The functional
range of integration for x1 will be range1(x2,x3), x2 will be
range2(x3), and so on.
args - optional sequence of arguments. Contains only arguments of
func beyond those over which the integration is being performed.
opts - optional sequence of options for SciPy's quad. Each element
of opts may be specified as either a dictionary or as a function
that returns a dictionary similarly to ranges. opts must either
be left empty (), or it must be the same length as ranges.
Options are passed in the same order as ranges, so opts[0]
corresponds to integration over the first argument of func, and
so on. The full_output option from quad is not available, due to
the difficulty of consolidating the large number of additional
outputs. For reference, the default options from quad are:
- epsabs = 1.49e-08
- epsrel = 1.49e-08
- limit = 50
- points = None
- weight = None
- wvar = None
- wopts = None
(As of Apr 2013)
_depth - used to determine level of integration. Should be omitted
by the user, except for debugging purposes.
_int_vars - contains values of integration variables in inner
integration loops. Should not be used manually except for
debugging.
Returns:
out - value of multiple integral in the specified range.
abserr - estimate of the absolute error in the result. The
maximum value of abserr among all the SciPy quad evaluations.
"""
global abserr
if _depth == 0:
abserr = None
if not (len(opts) in [0,len(ranges)]):
raise Error('opts must be given for all integration levels or none!')
total_args = _int_vars+args
# Select the range and opts for the given depth of integration.
ind = -_depth-1
if callable(ranges[ind]):
current_range = ranges[ind](*total_args)
else:
current_range = ranges[ind]
if len(opts) != 0:
if callable(opts[ind]):
current_opts = opts[ind](*total_args)
else:
current_opts = opts[ind]
else:
current_opts = {}
try:
if current_opts["full_output"] != 0:
raise Error('full_output option is disabled!')
except(KeyError):
pass
# Check to make sure that all points lie within the specified range.
try:
for point in current_opts["points"]:
if point < current_range[0] or point > current_range[1]:
current_opts["points"].remove(point)
except(KeyError):
pass
if current_range is ranges[0]:# Check to see if down to 1-D integral.
func_new = func
else:
# Define a new integrand.
def func_new(*_int_vars):
return mul_quad(func,ranges,args=args,opts=opts,
_depth=_depth+1,_int_vars=_int_vars)
out = quad(func_new,*current_range,args=_int_vars+args,**current_opts)
# Track the maximum value of abserr from quad
if abserr is None:
abserr = out[1]
if out[1] > abserr:
abserr = out[1]
if _depth == 0:
return out[0],abserr
else:
return out[0]
if __name__=='__main__':
func = lambda x0,x1,x2,x3 : x0**2+x1*x2-x3**3+numpy.sin(x0)+(
1 if (x0-.2*x3-.5-.25*x1>0) else 0)
points=[[lambda (x1,x2,x3) : .2*x3+.5+.25*x1],[],[],[]]
def opts0(*args,**kwargs):
return {'points':[.2*args[2]+.5+.25*args[0]]}
print mul_quad(func,[[0,1],[-1,1],[.13,.8],[-.15,1]],
opts=[opts0,{},{},{}])
-------------- next part --------------
A non-text attachment was scrubbed...
Name: myquad.py
Type: text/x-python-script
Size: 5736 bytes
Desc: not available
Url : http://mail.scipy.org/pipermail/scipy-dev/attachments/20130425/4bf9bee8/attachment-0001.bin
More information about the SciPy-Dev
mailing list