[SciPy-dev] double exponential integrator

[Thouis (Ray) Jones]

Renamed the function to quad_de.
Changed its signature to be closer to quadrature().
Improved reporting when the func() returns nonfinite values.
Added vectorization (via the vectorize1 functions in quadrature.py.
Filed a trac ticket: #895.

The generation of the weights and abscissas is not particularly
costly, but I believe they do fail in double precision arithmetic.  If
it's really desirable to have them created dynamically (to allow more
levels of integration, for instance), it would probably require
writing a Taylor expansion or similar.  I've left them in a separate
file for now, due to their reliance on mpmath.

In the future, this code could accept arbitrary limits (0, 1, or 2
infinite limits) and use whichever abscissa and weight functions are
appropriate, with C code to generate them on demand to arbitrary
levels.  I would like to gauge interest in this code in general before
adding more functionality.  The interface would be backwards
compatible, with possibly one new keyword argument to handle the two
different versions of integrals from 0 to infinity (standard versus
exponential falloff).

Ray Jones


On Fri, Mar 20, 2009 at 05:46, Pauli Virtanen <pav at iki.fi> wrote:
> Hi,
>
> Thu, 19 Mar 2009 13:55:58 -0400, Thouis (Ray) Jones wrote:
>> I'm soliciting feedback on an implementation of integration using the
>> double exponential transform.  I've tentatively placed it in
>> scipy.integrate as de_integrate.
>
> +1 for including this in Scipy 0.8.0.
>
> Some quick comments (I'll try to find time for better comments later):
>
> - Vectorization: you are using a list comprehension in
>  doubleexp.py:contribution_at_level and elsewhere
>
>  These statements could be vectorized -- I believe you can also require
>  that the integrand function `f` can evaluate many points at the same
>  time and return an array. Could be an useful speedup.
>
> - Is generating the _abscissas_raw and _weights_raw costly?
>
>  I see that you use `mpmath` to prepare these. Does the generation
>  fail in double precision?
>
>  If not, it might be better to generate them when the integration
>  function is first used. (Also, it might be nice to put also the
>  generator functions in the same doubleexp.py; these are not long files.)
>
> - Function name `de_integrate`: Perhaps it should be `quad_de` or
>  something similar, since it's a quadrature, and the "basic" quadrature
>  in scipy.integrate is called `quad`.
>
> - I'm not sure about if a `full_output` switch is good API.
>  Several Scipy functions do currently use something like that,
>  but perhaps it would be better not to introduce more...
>
>> It's available at this url and branch:
>> http://broad.mit.edu/~thouis/scipy.git DEintegrator
>>
>> (assuming I've set up my git repository correctly, which is quite
>> possibly not the case.)
>
> The repository works OK.
>
> --
> Pauli Virtanen
>
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>