multiple functions from scipy.special do not work with quantities

[namurphy]

I just tried to use SciPy's error function, and Astropy politely requested that I raise this issue here:

>>> from astropy import units
>>> from scipy.special import erf, erfc
>>> r1 = 5*u.m
>>> r2 = 6*u.m
>>> erf(r1/r2)
TypeError: Unknown ufunc erf.  Please raise issue on https://github.com/astropy/astropy
>>> erfc(r1/r2)
TypeError: Unknown ufunc erfc.  Please raise issue on https://github.com/astropy/astropy
>>> erf(r1.value/r2.value)  # workaround
0.76140717068356445

In scipy.special, ufuncs that I would prioritize include:

• [ ] Error function, etc
• [ ] Bessel functions
• [ ] Gamma and related functions
• [ ] Legendre functions
• [ ] zeta and zetac - Riemann zeta functions

Thanks!

[namurphy]

Note: this is a superset of issue https://github.com/astropy/astropy/issues/4092.

[mhvk]

6843 reminded me that it would be really nice to support scipy's special functions; it is not difficult in principle; see discussion in #4092 (which also gives a work-around).

[mhvk]

@kyleaoman - moving discussion to this open issue, and trying to be a bit more constructive: how about starting with a list of scipy.special functions that are dimensionless to dimensionless, and just add those? It could be something like the following at the end of quantity_helpers.py:

try:
    import scipy
except ImportError:
    pass
else:
    for special in ['erf', 'erfc', ...]:
        UFUNC_HELPERS[getattr(scipy.special, special)] = helper_dimensionless_to_dimensionless

but also needs some tests (sadly, I fear that very good requirement is the reason I didn't get around to this...)

[kyleaoman]

I had a go at starting to add support for scipy.special functions, I now have a branch in a state where I can pick away at it for the next couple of weeks in spare moments. This is my first foray into collaborative github work (though I use it for my own projects). I sense I should mark this as 'in progress' somehow?

[pllim]

@kyleaoman , you can always open a "WIP" pull request and add [ci skip] in your commits to skip the CI builds until you are ready.

[mhvk]

If you want me/others to have a look, the easiest would be to make it a pull request -- see http://docs.astropy.org/en/latest/development/workflow/development_workflow.html#pull-request; then we can add a label that it is in progress (do have a look at the documentation, and let us know if it is unclear). But if it is not in that state yet, you can also paste a link here.

... ah, I see @pllim was ahead of me ... was about to add something about testing.

[kyleaoman]

Sounds good, I think I can start a wip pull request shortly. I don't currently have [ci skip] on my commits, but each commit does pass the tests on my system. I suspect it will break on other systems since scipy.special seems to have shuffled around which functions are implemented as ufuncs over time (at least, the docs of the current version 1.0.0 only sort of reflect what's in the source...). Is it better to go back and add the tag to my commit messages in this scenario?

[pllim]

Is it better to go back and add the tag to my commit messages in this scenario?

Well, given that the CI traffic is not too bad right now, it does not really matter. If you want to do it anyway, only the last commit needs the "[ci skip]" text.

[mhvk]

6852 included a number, but far from all, scipy.special functions. With it merged, though, adding further functions is not much effort, so I changed the label. It would be particularly easy to just add more functions that take and give dimensionless numbers.

[pllim]

Is this still an issue, @mhvk or @namurphy ?

[bsipocz]

erf and erfc are working from the example, but I haven't checked the rest in the list.

[namurphy]

Many of the high priority functions from scipy.special do work now (thanks, @kyleaoman!). There are still a lot of special functions that are still incompatible with dimensionless Quantity instances, so this issue should stay open.

[bsipocz]

@namurphy - if you have a list of them it would be useful to update the original description above. Thanks!

[namurphy]

What I did to quickly check this was run through all of the non-private functions in scipy.special and see if they worked with 1 * u.dimensionless_unscaled as the sole argument. The list of functions that did not work was very long, but this list would have also included any functions that require two or more arguments and any functions for which 1 is an invalid input.

from astropy import units as u
from scipy import special

dimensionless_quantity = 1 * u.dimensionless_unscaled

for key in sorted(special.__dict__):
    try:
        if key[0] != '_':
            special.__dict__[key](dimensionless_quantity)
    except Exception:
        print(key)

[pllim]

Using astropy 3.2rc1 and numpy 1.16.3, attached is the output that I get. I am not convinced that we need to support all of them. Can you pick out those that are important to you?

specialsauce.txt

UPDATE: Ah, I see that you have a prioritized list above.

[namurphy]

Yeah, supporting all of them does seem like a bit much, though it is hard to know ahead of time which functions will be needed and which won't. I keep wondering if there is a way to resolve this more generally without having to take into account all of the different special cases, though that sounds difficult too. The main functions I remember needing are erf, erfc, and wofz, which all appear to work now.

[mhvk]

My feeling would be to just leave this open and see if someone has the energy to do it...

[luisval]

I'm interested, is it still open?

[kyleaoman]

@luisval yes it's still open, and afaik no one has worked on it since the last PR was merged, so go for it.

[mhvk]

Indeed, and with the infrastructure put in place by @kyleaoman it should be quite easy - be sure for every added function to also add corresponding tests (and think carefully about the required units, of course!)

[astrojuanlu]

Prospective contributors: take a look at the prioritized list on the first comment and then read https://github.com/astropy/astropy/blob/master/astropy/units/quantity_helper/scipy_special.py

[DiegoAsterio]

I am concerned my change won't close the issue but adds functionality to some error and gamma functions. I also realized there are some functions that are not being tested. On the other hand when I try to add zeta and zetac to @mhvk list I get an error, particularly that some intermediate function is not a universal function.

[mhvk]

@DiegoAsterio - indeed, your change won't close the issue, but it will make more functions work, so it is a good step!

In the end, I think this will work two ways - people finding functions they'd like to work and adding them, and someone going through the whole list carefully thinking which ones should work (and how).