[SciPy-Dev] Discussion on a new interface for multidimensional quadrature

Tue Aug 13 15:49:13 EDT 2013

Right, that's what I thought, too. The thing is that dbl- and tplquad DON'T do that. If you carefully at the calling sequence for tplquad, 
-----------
scipy.integrate.tplquad(func, a, b, gfun, hfun, qfun, rfun, args=(), epsabs=1.49e-08, epsrel=1.49e-08)[source]¶

Return the triple integral of func(z, y, x) from x=a..b, y=gfun(x)..hfun(x), and z=qfun(x,y)..rfun(x,y)
------------

The only way that the underlying quad routine can accept qfun(x,y) and rfun(x,y) as limits is if they are evaluated to be constants. That is, integration over z must be the innermost loop of integration, with fund(z,y,z), qfun(x,y), and rfun(x,y) receiving (x,y) passed in as arguments from the outer loops. 

I agree that the changes I just put in are kind of ridiculous, but (as far as I can tell) they exactly duplicate what is done in dbl- and tplquad.

If we really want to put in some kind of integration order, then that's also not too hard, though a bit more than just reversing the lists. It's just something that I don't see a real, compelling reason to do. 

N

On Aug 13, 2013, at 1:42 PM, Ralf Gommers <ralf.gommers at gmail.com> wrote:

> 
> 
> 
> On Tue, Aug 13, 2013 at 7:20 PM, Ralf Gommers <ralf.gommers at gmail.com> wrote:
> 
> 
> 
> On Tue, Aug 13, 2013 at 6:38 PM, Nathan Woods <charlesnwoods at gmail.com> wrote:
> Since I didn't hear from anyone, I went ahead and implemented the 'reversed' argument as requested. As expected, it was pretty simple. However, I really don't think it adds anything of value to the routine, and it does provide another avenue for error. My personal preference is to roll back the change and leave nquad as forward integration only. I'd appreciate anyone's thoughts on the matter.
> 
> N
> 
> On Aug 12, 2013, at 3:39 PM, Nathan Woods <charlesnwoods at gmail.com> wrote:
> 
>> In the process of implementing the optional argument, I noticed something that I had missed before.
>> 
>> tplquad(f(z,y,x),x0,x1,y0,y1,z0,z1) gives the same answer, implemented in the same way (nesting order) as 
>> 
>> nquad(f(z,y,x),[[z0,z1],[y0,y1],[x0,x1]]).
> 
> 
> This only helps if the integration limits are constants, so I think 'reversed' still helps.
> 
> Now that I looked at how you implemented this, indeed that doesn't help much. What I was thinking is that 'reversed' would reverse the order of integration, so first over y instead of first over x for a function f(x, y). I suspect that's also what Evgeni meant.
> 
> Ralf
> 
>> 
>> That means that the only change here is in the calling sequence. The actual nesting order of integration remains the same (first argument to last). This means that there would be no need to wrap an existing function to reverse arguments, for instance.
>> 
>> I'm not entirely sure what was meant by mapping to dblquad/tplquad, though. Does that mean something like having dblquad/tplquad become wrappers for nquad? It seems like that would be something more easily done in the wrappers, rather than complicating the underlying function. Thoughts?
>> 
>> Nate
>> 
>> 
>> On Aug 11, 2013, at 9:00 PM, Nathan Woods <charlesnwoods at gmail.com> wrote:
>> 
>>> I like the optional argument idea. I think it gives a good balance between a clean interface and legacy support. It's also not too difficult to implement. 
>>> 
>>> Nate
>>> 
>>> On Aug 10, 2013, at 8:48 AM, Ralf Gommers <ralf.gommers at gmail.com> wrote:
>>> 
>>>> 
>>>> 
>>>> 
>>>> On Fri, Aug 9, 2013 at 5:55 AM, Evgeni Burovski <evgeny.burovskiy at gmail.com> wrote:
>>>> How about adding an optional argument to nquad, specifying the order of integration limits? Defaulting to the 'forward' order. Something like
>>>> 
>>>> nquad(f(x,y), [a, b], [g(x), h(x)])
>>>> being the default, and
>>>> 
>>>> nquad(f(x,y), [a, b], [g(y), h(y)], integration_order='reversed')
>>>> interpreting the limits the way dblquad/tplquad do.
>>>> 
>>>> 
>>>> That makes sense to me, both to make it easier to map to dblquad/tplquad and to avoid writing an argument-reversing wrapper when you want to integrate an existing function.
>>>>  
>>>> Evgeni
>>>> 
>>>> On Aug 9, 2013 11:35 AM, "Nathan Woods" <charlesnwoods at gmail.com> wrote:
>>>> Clear documentation is always important, so I don't think there will be any problems there. Does anyone else have thoughts they'd like to share?
>>>> 
>>>> Nate
>>>> 
>>>> On Aug 5, 2013, at 11:35 PM, josef.pktd at gmail.com wrote:
>>>> 
>>>> > On Mon, Aug 5, 2013 at 4:24 PM, Ralf Gommers <ralf.gommers at gmail.com> wrote:
>>>> >>
>>>> >>
>>>> >>
>>>> >> On Fri, Aug 2, 2013 at 7:19 PM, Nathan Woods <charlesnwoods at gmail.com>
>>>> >> wrote:
>>>> >>>
>>>> >>> I'm copying the most relevant parts over from a discussion on GitHub,
>>>> >>> having to do with a new n-dimensional integration tool. The discussion is
>>>> >>> whether or not to make nquad's new interface consistent with dbl-and
>>>> >>> tplquad, even though their precedent is arguably more confusing to users.
>>>> >>>
>>>> >>> nquad is a recursively defined wrapper to scipy.integrate.quad that allows
>>>> >>> for n-dimensional integration with almost the entire list of arguments
>>>> >>> allowed by quad. It is therefore more general than either dbl- or tplquad,
>>>> >>> and encompasses their use cases, while allowing the user greater control
>>>> >>> over the integration. The interface is given: nquad( f(x0,x1,…xn,t0,…tm)
>>>> >>> ,[[x0_lo,x0_hi],…[xn_lo,xn_hi]],[t0,…tm],[opts0,…optsn]).
>>>> >>> The ranges of integration are given in the same order as the arguments in
>>>> >>> the function to be integrated. This ordering is consistent with both Matlab
>>>> >>> and Octave, and it is similar to what is encountered in other SciPy packages
>>>> >>> such as fft. It is also what makes sense to new users.
>>>> >>>
>>>> >>> dbl- and tplquad, however, reverse the order of arguments in f. From the
>>>> >>> documentation for tplquad:
>>>> >>>
>>>> >>> Return the triple integral of func(z, y, x) from x=a..b,
>>>> >>> y=gfun(x)..hfun(x), and z=qfun(x,y)..rfun(x,y)
>>>> >>>
>>>> >>> This is confusing, and contrary to common practice.
>>>> >>>
>>>> >>> It is proposed that nquad's calling sequence be left as is, and that the
>>>> >>> use of dbl- and tplquad be deprecated in some future version of SciPy, and
>>>> >>> that nquad be recommended for multidimensional integration in new code. dbl-
>>>> >>> and tplquad would of course be retained, though not recommended, in order to
>>>> >>> avoid breaking existing code.
>>>> >>>
>>>> >>> Wow, that came out all legalistic. Anyway, what do you all think about
>>>> >>> this change? I have a horse in the race (I wrote most of nquad), but I
>>>> >>> really think that maintaining the backwards style from dbl- and tplquad and
>>>> >>> enshrining it in new code is a bad idea. This is a perfect time to get away
>>>> >>> from it, if that's what we want to do.
>>>> >>
>>>> >>
>>>> >> My first instinct was that nquad should match what's already there in
>>>> >> scipy.integrate, but I've changed my mind. I always have to look at the
>>>> >> dblquad/tplquad documentation to get it right, so those functions don't have
>>>> >> a good API imho. Therefore just making nquad do the logical thing (i.e.
>>>> >> ``f(x, y ,z)``) makes sense to me.
>>>> >>
>>>> >> Let's not use the word "deprecate" for dblquad/tplquad, since they shouldn't
>>>> >> be removed. Instead just recommend `nquad` as the function to use in the
>>>> >> docs.
>>>> >>
>>>> >> Any other opinions?
>>>> >
>>>> > How do you know the integration order in the new interface?
>>>> >
>>>> > The current dblquad, tplquad looks to me like writing an integral
>>>> >
>>>> > integral_x integral_y integral_z  func3d(z, y,x)  dz dy dz
>>>> > where integration limits are x=a..b, y=gfun(x)..hfun(x), and
>>>> > z=qfun(x,y)..rfun(x,y)
>>>> >
>>>> > based on reading of the docstrings for tplquad. So this doesn't look
>>>> > very strange to me.
>>>> 
>>>> I guess that was the reason why dblquad/tplquad are implemented the way they are. The problem with that is that no one writes functions like f(z, y, x), everyone writes f(x, y, z). So it's plain unnatural.
>>>>  
>>>> >
>>>> > The suggested description
>>>> > ``
>>>> > The interface is given: nquad( f(x0,x1,…xn,t0,…tm)
>>>> > ,[[x0_lo,x0_hi],…[xn_lo,xn_hi]],[t0,…tm],[opts0,…optsn]).
>>>> > The ranges of integration are given in the same order as the arguments
>>>> > in the function to be integrated. 
>>>> > ```
>>>> >
>>>> > Does this mean we integrate over x0 first (innermost integral) or last
>>>> > (outermost integral) ?
>>>> 
>>>> x0 first, that should be made clearer in the docstring.
>>>> 
>>>> Ralf
>>>> 
>>>> _______________________________________________
>>>> 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/20130813/12903175/attachment.html>