[SciPy-Dev] Request for extension to scipy.integrate
Nathan Woods
charlesnwoods at gmail.com
Mon May 20 21:02:26 EDT 2013
A new PR is up for an n-dimensional integration wrapper for scipy.integrate.quad, along with unit tests and supporting documentation. Feel free to take a look and let me know what you think.
N
On May 8, 2013, at 11:19 AM, Nathan Woods <charlesnwoods at gmail.com> wrote:
> Thanks again for the feedback from before. I've implemented the changes suggested, and reformatted the documentation as requested. I admit that I'm not exactly sure what you mean by unit tests, other than a few comparisons with another numerical integrator such as Mathematica. Nor do I really know how to implement something like that. Also, what's a PR? Pull request, maybe?
>
> N
>
> On May 5, 2013, at 4:12 PM, Ralf Gommers <ralf.gommers at gmail.com> wrote:
>
>>
>>
>>
>> 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
>>
>>
>> _______________________________________________
>> 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/20130520/073f47ee/attachment.html>
More information about the SciPy-Dev
mailing list