[SciPy-Dev] Request for extension to scipy.integrate
Ralf Gommers
ralf.gommers at gmail.com
Sun May 5 18:12:46 EDT 2013
On Sun, May 5, 2013 at 2:42 PM, Nathan Woods <charlesnwoods at gmail.com>wrote:
> 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?
>
A class would be preferable to globals and keyword arguments not intended
for the user. I don't see a performance issue with that here.
> 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?
>
The latter, in addition to what you currently have. If I want to pass the
same keyword to each integration steps, then I shouldn't be forced to
create N times the same dict.
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.
>
Maybe it's unavoidable that it becomes a mess - dblquad/tplquad also don't
have a full_output kw. Not sure if there's any point in just returning the
same thing as quad for the outer integration step maybe.
Ralf
> 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
>
>
> _______________________________________________
> 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/20130506/1b0d417f/attachment.html>
More information about the SciPy-Dev
mailing list