[SciPy-Dev] SciPy-Dev Digest, Vol 195, Issue 3

Fri Jan 10 21:52:31 EST 2020

> On Jan 4, 2020, at 12:00 PM, scipy-dev-request at python.org wrote:
> 
> Send SciPy-Dev mailing list submissions to
>    scipy-dev at python.org
> 
> To subscribe or unsubscribe via the World Wide Web, visit
>    https://mail.python.org/mailman/listinfo/scipy-dev
> or, via email, send a message with subject or body 'help' to
>    scipy-dev-request at python.org
> 
> You can reach the person managing the list at
>    scipy-dev-owner at python.org
> 
> When replying, please edit your Subject line so it is more specific
> than "Re: Contents of SciPy-Dev digest..."
> 
> 
> Today's Topics:
> 
>   1. Wright's Generalized Bessel Function (Christian Lorentzen)
>   2. Re: Symmetric/Hermitian Eigenvalue Problem API and LAPACK
>      wrappers change proposal (Ilhan Polat)
>   3. Re: Wright's Generalized Bessel Function (Joshua Wilson)
> 
> 
> ----------------------------------------------------------------------
> 
> Message: 1
> Date: Sat, 4 Jan 2020 12:17:19 +0100
> From: Christian Lorentzen <lorentzen.ch at gmail.com>
> To: scipy-dev at python.org
> Subject: [SciPy-Dev] Wright's Generalized Bessel Function
> Message-ID: <7d671672-cbe7-7517-6228-f1137e307377 at googlemail.com>
> Content-Type: text/plain; charset="utf-8"; Format="flowed"
> 
> Dear Scipy Developers
> 
> I started the PR [1] in order to make Wright's generalized Bessel 
> function a member of scipy.special (see PR for more infos). This is a 
> necessary ingredient for the implementation of Tweedie distributions [2] 
> in scipy.stats and an interesting function in its own right.
> 
> AFAIK there is no public reference implementation. As it is an entire 
> function in z, it is possible to use its defining series expansion:
> - Every term needs a call to 1/Gamma(...).
> - It is not very clear when to stop the expansion.
> - There are asymptotic expansions, too. I got very unsatisfactory 
> results from them.
> 
> My questions:
> 1) Do you want to include this function in scipy.special?
> 2) API: Optional arguments `tol` and `max_iter` seem to fit a dynamic 
> convergence check,
> |lambertw also has `tol`.
> |3) Implementation: Is it ok to start with a pure python implementation 
> and later on change to, e.g. Cython? Or is a Cython version mandatory 
> for merging a PR?
> 4) How to test this function? There are very little special cases (one 
> with scipy.special.iv)?
> 
> References:
> [1] https://github.com/scipy/scipy/pull/11313
> [2] https://github.com/scipy/scipy/issues/11291
> 
> 
> Kind regards,
> Christian
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20200104/b6ba89d2/attachment-0001.html>
> 
> ------------------------------
> 
> Message: 2
> Date: Sat, 4 Jan 2020 14:12:15 +0100
> From: Ilhan Polat <ilhanpolat at gmail.com>
> To: SciPy Developers List <scipy-dev at python.org>
> Subject: Re: [SciPy-Dev] Symmetric/Hermitian Eigenvalue Problem API
>    and LAPACK wrappers change proposal
> Message-ID:
>    <CAEBuzr-+THd6U52vjerA8-4ivxNywJ2XFHsTCEhO6OOzy81PkA at mail.gmail.com>
> Content-Type: text/plain; charset="utf-8"
> 
> Dear all,
> 
> Regarding this, PR is seemingly stabilized. Please review it if you have
> time https://github.com/scipy/scipy/pull/11304 as it is kind-of long due to
> the new exposed functionality.
> 
> Best,
> ilhan
> 
> 
> 
>> On Mon, Dec 23, 2019 at 12:00 AM Ilhan Polat <ilhanpolat at gmail.com> wrote:
>> 
>> Dear all,
>> 
>> I've finally managed to create some time to put together the current
>> status of the eigh, eigvalsh functions and possible ways to improve and
>> maintain them.
>> 
>> The issue is here
>> 
>> https://github.com/scipy/scipy/issues/6502#issuecomment-241437506
>> 
>> I would really appreciate any feedback for this. It's in WIP format and
>> I'll add more information as I go along.
>> 
>> Happy holidays to all!
>> 
>> Best,
>> ilhan
>> 
> -------------- next part --------------
> An HTML attachment was scrubbed...
> URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20200104/dd32dba6/attachment-0001.html>
> 
> ------------------------------
> 
> Message: 3
> Date: Sat, 4 Jan 2020 08:58:08 -0800
> From: Joshua Wilson <josh.craig.wilson at gmail.com>
> To: SciPy Developers List <scipy-dev at python.org>
> Subject: Re: [SciPy-Dev] Wright's Generalized Bessel Function
> Message-ID:
>    <CAKFGQGzTyAB9aYB3WEjXm1m2Gv6txBx826DW5VCWWZCuScZRkw at mail.gmail.com>
> Content-Type: text/plain; charset="UTF-8"
> 
> Some thoughts here:
> 
>> it is possible to use its defining series expansion
> 
> It looks like that series is going to suffer from cancellation for (at
> least) z along the negative real axis, so the asymptotic series will
> be needed in some fashion. In general, this function is reasonably
> complex. You will probably need to use multiple methods across various
> regions.

There are a couple of papers on numerical implementations of wrights Bessel. They don’t appear to be much work in the area nor applications beyond fractional differential equations. So I’d second @person142 on the private method in special.

> 
>> 2) API: Optional arguments `tol` and `max_iter` seem to fit a dynamic convergence check, lambertw also has `tol`.
> 
> Including tol in lambertw was a mistake IMO, and I would not like to
> see it here. Generally we strive to have all functions be as correct
> as possible instead of fiddling with tols.
> 
>> 3) Implementation: Is it ok to start with a pure python implementation
> 
> As mentioned in the PR, this is not ok IMO, for some reasons outlined there.
> 
>> 4) How to test this function? There are very little special cases (one with scipy.special.iv)?
> 
> One approach is to do what Boost does in situations like this-they
> make a very naive implementation using arbitrary precision (use mpmath
> for this) and then compare that implementation to the optimized one.

I second the mpmath approach. Another approach is to find some papers that have some example data and use that data as test cases and make sure code reproduces results. 
Upshot if this is the lit review process you might find an algorithmic approach that works better in a region.

>> 1) Do you want to include this function in scipy.special?
> 
> It's in the DLMF, so I'd say it's reasonable to include. *But* I think
> you are underestimating how hard it is going to be to create a good
> implementation of this function (and the bar for special is high).
> Depending on your level of experience with implementing special
> functions, it could potentially be  *a lot* of work. I don't know
> about the details of the Tweedie distributions, but you might want to
> think about whether you need a generic implementation of this
> function. If you only need it for some regions e.g. we can just create
> a private implementation for those specific regions, potentially
> reducing the amount of work involved.

For the tweedie the key parameter is p which is supported on (1,nifty) and (-nifty, 0). From the brief lit review I did on google scholar publicly available docs, numerically coverage of the whole real support of p is probably going to be a lot of work. The most common applications of tweedie I’ve used in the past were with 1<=p<=2.
P=1 so it might make sense to work towards supporting those argument values first. Then expanding coverage to other values of p.

> 
>> On Sat, Jan 4, 2020 at 3:18 AM Christian Lorentzen
>> <lorentzen.ch at gmail.com> wrote:
>> 
>> Dear Scipy Developers
>> 
>> I started the PR [1] in order to make Wright's generalized Bessel function a member of scipy.special (see PR for more infos). This is a necessary ingredient for the implementation of Tweedie distributions [2] in scipy.stats and an interesting function in its own right.
>> 
>> AFAIK there is no public reference implementation. As it is an entire function in z, it is possible to use its defining series expansion:
>> - Every term needs a call to 1/Gamma(...).
>> - It is not very clear when to stop the expansion.
>> - There are asymptotic expansions, too. I got very unsatisfactory results from them.
>> 
>> My questions:
>> 1) Do you want to include this function in scipy.special?
>> 2) API: Optional arguments `tol` and `max_iter` seem to fit a dynamic convergence check,
>> lambertw also has `tol`.
>> 3) Implementation: Is it ok to start with a pure python implementation and later on change to, e.g. Cython? Or is a Cython version mandatory for merging a PR?
>> 4) How to test this function? There are very little special cases (one with scipy.special.iv)?
>> 
>> References:
>> [1] https://github.com/scipy/scipy/pull/11313
>> [2] https://github.com/scipy/scipy/issues/11291
>> 
>> 
>> Kind regards,
>> Christian
>> _______________________________________________
>> SciPy-Dev mailing list
>> SciPy-Dev at python.org
>> https://mail.python.org/mailman/listinfo/scipy-dev
> 
> 
> ------------------------------
> 
> Subject: Digest Footer
> 
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at python.org
> https://mail.python.org/mailman/listinfo/scipy-dev
> 
> 
> ------------------------------
> 
> End of SciPy-Dev Digest, Vol 195, Issue 3
> *****************************************