ENH: special: extend Spheroid Wave Functions domain to negative input

[danchendrickson]

Under SciPy.Special the spheroid wave functions are returning 'nan' for negative values of x. Should return a value.

Values to be returned are out of: https://dlmf.nist.gov/30.7#ii https://dlmf.nist.gov/30.7#iii

Reproducing code example:

import scipy.special as sp
sp.pro_rad1(0,2,.25,-0.5)[0]

Error message:

No error message is created, just the wrong returned value

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  ...

Scipy/Numpy/Python version information:

1.5.0 1.18.5 sys.version_info(major=3, minor=8, micro=3, releaselevel='final', serial=0)

[-]

[tupui]

Thanks for reporting. On latest this does not raise but return nan. Actually looking at the C (prolate_radial1_nocv_wrap), one must use x>1 otherwise it will just return nan. This seems to be incorrect as suggested by the docstrings. Changing the check, I have the following results:

import scipy.special as sp
sp.pro_rad1(0,2,.25,-0.5)
 (-0.0003420668035636386, -0.0041723915799936956)

Is this correct? If yes I will fix it.

cc @person142

[person142]

Something of an aside, but I’ll note that I have ambitions to replace the implementations of the spheriodal wave functions in the next month or so-the current implementations are limited in many ways. So they might not be the best place to sink effort into.

[tupui]

I am glad to read that you plan this @person142. As we are branching off very soon for 1.7, if the fix I proposed is correct, I would just do it to correct 1.7.

[danchendrickson]

I am not sure what the correct values are. The authoritative source on what they should be is the NIST Handbook of Mathematical functions: https://www.cambridge.org/catalogue/catalogue.asp?isbn=9780521192255

The 3 pages of tables of values for the functions are attached here.

[tupui]

I am not sure what the correct values are. The authoritative source on what they should be is the NIST Handbook of Mathematical functions: https://www.cambridge.org/catalogue/catalogue.asp?isbn=9780521192255

Thanks for that! Just checking before commenting further, @rgommers can we use this to validate results?

[rgommers]

Sure, it's just data - you can always use raw numbers, that's not copyrightable.

[tupui]

Sure, it's just data - you can always use raw numbers, that's not copyrightable.

Thanks, I was not sure about that.

In this case I will do a small PR to mitigate the issue. Mitigate as the values I get a close but not equal (2-3 correct digits only). Still better than just outputting NaNs IMO.

[tupui]

Ok here are my finding according to this ressource:

• obl_rad1: looks ok
• obl_rad2: seems only correct for x=0
• pro_rad1 and pro_rad2: according to the table, x can be >1 although the documentation says no. Without constrains seems ok.

So I can propose to remove the check on x for pro_rad and use abs(x)>=1.0 for obl_rad (so only valid for abs(x)<1). Does it make sense @person142? Or we just wait for your re-write? I would not add tests as you would re-write it and, most important, I have no clue about this.

[person142]

Sorry I can’t be of much help before the 1.7 branching-I’m on the road and mobile-only at the moment.

[tupui]

Sure no problem. Well in that case as I really don't know what I am doing I better not do anything and release something even more wrong 😅

@danchendrickson feel free to open a PR if you can find someone to check this. Everything is in scipy/special/specfun_wrappers.c and the tests scipy/special/test_basic.py. Then for example for pro_rad1 the test on x is L829.

[dschmitz89]

Can this be closed now after #21441 @steppi ?

[steppi]

Can this be closed now after #21441 @steppi ?

Thanks for pointing this out. Not yet. The issue here is that we only allow evaluating this functions when the arguments give valid points in prolate or oblate spheroidal coordinates respectively, but they can be analytically continued to extend their domains, and @danchendrickson seems to have a use case for that. I think it should be recategorized as an enhancement request though.

[dschmitz89]

Someone should train a LLM on @ilayn s F77 translations, so that the world finally gets something better than f2c.