[SciPy-Dev] Request for extension to scipy.integrate
Nathan Woods
charlesnwoods at gmail.com
Sun May 5 08:42:39 EDT 2013
Ralf,
Thank you for the feedback and suggestions. I'll get started on a few of
them right away. There are a few questions you brought up that I'll try and
answer, while I think about ways to improve them.
First, the debug arguments and the global variable abserr. The debug
arguments are actually required by the code, but they are not intended to
be used by the user. They carry around information that is required by the
nested integration, and that would otherwise have to be included as class
variables. They are not intended to be used by the user at all, really. I
could certainly get rid of them if I wrote this routine as a class instead
of a def, but I've heard that there are performance reasons not to use
classes in a nested, iterative, routine like this one. Thoughts?
For options specifications, it sounds like you're advocating a dict if
options containing keys like "all", 1, 2, ... to allow the user to specify
options only at a particular level of integration, or at them all. Did I
understand you right? Or were you just suggesting that a single opts
specification for everything should be a viable input?
As for the full-output option, there would need to be some way to reduce
the potentially thousands of outputs in a total integration to something
manageable. I'll certainly look at it again, though.
Thanks again for your feedback.
N
On May 5, 2013, at 3:08, Ralf Gommers <ralf.gommers at gmail.com> wrote:
Hi Nathan,
On Fri, Apr 26, 2013 at 2:53 AM, Nathan Woods <charlesnwoods at gmail.com>wrote:
> 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.
>
Including an n-dim integration function sounds good to me. I've tried your
implementation with a few simple functions and it seems to give the right
results.
Some comments on the API and code:
- the debug keyword arguments should be removed
- function name: maybe `nquad` or `ndimquad`? "mul" can mean multiple or
multiply
- don't use a global for abserr
- why must `opts` be of the same length as `ranges`? Repeating it in case
it has only one element would make sense. For example, this should work:
``opts={'full_output', True})``, now I have to write
``opts=[{'full_output', True}, {}, {}, {}])``.
- why is full_output disabled?
Next steps would be to add good unit tests, adapt the documentation to
Numpy format (see
https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt), add
an example or two to the docstring, and submit a PR for this.
Cheers,
Ralf
>
> 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,{},{},{}])
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>
>
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev at scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20130505/a8b6371e/attachment.html>
More information about the SciPy-Dev
mailing list