[SciPy-Dev] Hankel transforms, again

[Ralf Gommers]

On Mon, Mar 17, 2014 at 11:48 AM, Tom Grydeland <tom.grydeland at gmail.com>wrote:

> On Mon, Mar 17, 2014 at 12:12 AM, Ralf Gommers <ralf.gommers at gmail.com>
> wrote:
> > Hi Tom. The commit looks fine. I guggest you send this as a pull
> request, so
> > it's easier to review.
>
> Will do.
>
> > A few comments already:
> > - the API looks slightly awkward, instantiating the class is basically a
> > do-nothing operation. You'd normally do this with a plain function that
> has
> > a ``method='anderson'`` keyword.
>
> Agreed -- will change
>
> (I have had more stuff in the class previously, such as keeping track
> of scaled function arguments for various arguments 'B' of the
> transform, but the resulting interface was sufficiently hard to
> explain that I decided to cut it down to its minimum, I didn't realise
> there was nothing left :D )
>
> > - the hankel0, hankel1 and hankel01 methods look unnecessary.
>
> Here, I disagree.  They are merely convenience functions, but it is
> much more obvious to say "hfunc = hankel0(func, bvals)" than "hfunc =
> hankel(func, bvals, order=[0])[0]", and yet I wanted the actual worker
> routine to be factored out, since the function evaluations can be
> reused for both transforms.  I could agree that hankel01 is redundant,
> if the default order argument to 'hankel' becomes '[0, 1]'.
>

It still doesn't look right to me. Either you have three functions, or one
function which takes an `order` keyword. Not both. Also the keywords
`method, ybase, wt0, wt1` are all doing nothing right now. hankel() is
simple enough that there's little value in providing those for the few
souls that take the effort to track down alternative coefficients.

After looking at this in some more detail, it looks to me like `hankel` is
too generic a name - that's what I would expect as the name for `f(x, nu)`
instead of `(f(x, 0), f(x, 1))`.

 > - The file names of all the Python files you add in scipy/signal/ should
> > start with an underscore, so it's clear that they are private.
>
> Okay, will change.  How do I expose the callable functions then?
>

- $ git mv hankel.py _hankel.py
- in signal/__init__.py:   from _hankel import hankel0, ....


> > - The docstrings could use an example and should be formatted according
> to
> > https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt
>
> Okay, will change.  Will the unit tests be suitable examples?
>

I would actually define a normal function instead of a lambda, one which
has a closed-form solution like f(x) = x**(-0.5) * exp(-x), then plot that
solution together with the numerical evaluation of it.


> > - the ``try: iter(B)`` blocks can be written simply as ``B =
> > np.asarray(B)``. The way it's now, B as a list will raise an error later
> on.
>
> Okay, will change.
>
> > - The for-loop over all values of B looks quite inefficient.
>
> It is.
>
> It is possible to arrange your sequence of 'B' values such that you
> can reuse a large portion of function evaluations for each additional
> transformed point, but that requires the caller to know in some detail
> the properties of the grid being used, and I have not implemented
> this.
>
> Various tricks are imaginable to overcome this, e.g. creating a grid
> covering all the desired B values and using some sort of interpolation
> on the results.  It is also possible, when you want to transform
> several functions that are related in some way, to keep intermediate
> evaluations to perform all transforms at once.  The interfaces quickly
> become muddled, however.


That indeed doesn't sound good.


>  For an initial inclusion, I'd like to make
> sure the interface is simple and usable, and the results predictable.
> An interface for better efficiency can come as a result of usage and
> experience.
>

Hmm, aiming for first time right would be better.


> There are also certain ranges of inputs where the transform is better
> performed using contour integration techniques.  I have not
> implemented that either.
>

There's actually a quite large body of literature on numerical evaluation
of the Hankel transform (see
http://www.sciencedirect.com/science/article/pii/0898122193900816). Some
digging turned up for example an improved algorithm from Anderson himself (
http://library.seg.org/doi/abs/10.1190/1.1442650) as well as
http://dl.acm.org/citation.cfm?id=317284 which the author's website says is
public domain (I think, see
http://homepages.tu-darmstadt.de/~wieder/hankel/hankel.html).

Some motivation for why the original Anderson algorithm is still suitable
to add to Scipy would be useful.

Ralf


> Thank you, and best regards,
>
> --
> Tom Grydeland
>   <Tom.Grydeland@(gmail.com)>
> _______________________________________________
> 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/20140321/c24e5faf/attachment.html>