Skip to content

Fix for special.hyp0f1, closes tickets 1750/1752. #345

Closed
wants to merge 5 commits into from

3 participants

@rgommers
SciPy member

Also add a better test for fdtri (function still needs work).

@pv pv and 1 other commented on an outdated diff Oct 31, 2012
scipy/special/basic.py
z = asarray(z)
- if issubdtype(z.dtype, complexfloating):
- arg = 2*sqrt(abs(z))
- num = where(z>=0, iv(v-1,arg), jv(v-1,arg))
- den = abs(z)**((v-1.0)/2)
- else:
- num = iv(v-1,2*sqrt(z))
- den = z**((v-1.0)/2.0)
- num *= gamma(v)
- return where(z==0,1.0,num/ asarray(den))
+ arg = 2 * sqrt(abs(z))
+ num = where(z >= 0, iv(a - 1, arg), jv(a - 1, arg))
+ den = abs(z)**((a - 1.0) / 2)
+ num *= gamma(a)
+ den[z == 0] = 1 # Avoid RuntimeWarning in next line
@pv
SciPy member
pv added a note Oct 31, 2012

Maybe better than using the where clause: num[z == 0] = 0, den[z == 0] = 1?

@rgommers
SciPy member
rgommers added a note Oct 31, 2012

Should be num[z == 0] = 1, but yes that's better than where.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv pv commented on the diff Oct 31, 2012
scipy/special/tests/test_basic.py
@@ -1143,7 +1146,15 @@ def test_h2vp(self):
assert_almost_equal(h2,h2real,8)
def test_hyp0f1(self):
- pass
+ # float input, expected values match mpmath
+ x = special.hyp0f1(3.0, [-1.5, -1, 0, 1, 1.5])
+ expected = np.array([0.58493659229143, 0.70566805723127, 1.0,
+ 1.37789689539747, 1.60373685288480])
+ assert_allclose(x, expected, rtol=1e-12)
+
+ # complex input
+ x = special.hyp0f1(3.0, np.array([-1.5, -1, 0, 1, 1.5]) + 0.j)
+ assert_allclose(x, expected.astype(np.complex), rtol=1e-12)
@pv
SciPy member
pv added a note Oct 31, 2012

A test testing that broadcasting in the function works would be good. That this works correctly is not obvious from the function implementation.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@pv
SciPy member
pv commented Oct 31, 2012

Looks good. Although z >= 0 probably should be taken as z.real >= 0.

@pv
SciPy member
pv commented Oct 31, 2012

Ideally, the infrastructure of scipy.special should be such that this function can be easily converted to a ufunc implemented in Cython. This probably means adding support for fused types and a .pxd definition of functions in Cephes.

@rgommers
SciPy member

Just for making all functions look the same? Sounds like a good principle, however in this particular case the Cython version won't be all that much faster, and with fused types you still get the exploding binary size issue.

@pv
SciPy member
pv commented Oct 31, 2012

I would be surprised if fused types in themselves lead to significantly more binary size growth than manually writing separate versions of the functions for double and complex double. Would have to test to see this...

@rgommers
SciPy member

That was what I concluded from the discussion of the ndimage.label PR (#261). I could be wrong though.

@rgommers
SciPy member

Indexing num/den with z == 0 doesn't work for scalar input, because asarray turns them into 0-D arrays. It's a minor backwards compatibility issue, but I think returning 0-D arrays is undesired anyway. Scalar input should either return scalar output or a 1-D array of size 1 imho. Do you agree?

If not, it's probably best to leave where in and filter the warning.

@pv
SciPy member
pv commented Oct 31, 2012

Sounds OK.

@rgommers
SciPy member

Hmm, broadcasting really doesn't work for this function. Example illustrating the problem with z == 0 indexing:

In [52]: x1 = np.arange(6).reshape((2,3))

In [53]: x2 = np.arange(3)+10

In [54]: x1 + x2
Out[54]: 
array([[10, 12, 14],
       [13, 15, 17]])

In [55]: (x1+x2)[x2 == 10]  # want [10, 13] here
Out[55]: array([[10, 12, 14]])

polygamma seems to have the same issue. Most straightforward way to make broadcasting work is with np.broadcast, but it will be slow. Other option is to not allow this input (require a.shape == b.shape if a and b are non-scalar).

Did I miss an easier way to make this work?

@pv
SciPy member
pv commented Nov 1, 2012

np.broadcast_arrays?

@rgommers
SciPy member
rgommers commented Nov 1, 2012

Duh, hadn't looked at that closely enough. Thanks.

@rgommers
SciPy member
rgommers commented Nov 1, 2012

Should be all fixed now. Polygamma function too.

@ewmoore ewmoore and 1 other commented on an outdated diff Nov 1, 2012
scipy/special/basic.py
@@ -446,20 +447,40 @@ def fresnel_zeros(nt):
raise ValueError("Argument must be positive scalar integer.")
return specfun.fcszo(2,nt), specfun.fcszo(1,nt)
-def hyp0f1(v,z):
- """Confluent hypergeometric limit function 0F1.
- Limit as q->infinity of 1F1(q;a;z/q)
- """
- z = asarray(z)
- if issubdtype(z.dtype, complexfloating):
- arg = 2*sqrt(abs(z))
- num = where(z>=0, iv(v-1,arg), jv(v-1,arg))
- den = abs(z)**((v-1.0)/2)
- else:
- num = iv(v-1,2*sqrt(z))
- den = z**((v-1.0)/2.0)
- num *= gamma(v)
- return where(z==0,1.0,num/ asarray(den))
+def hyp0f1(a, z):
@ewmoore
SciPy member
ewmoore added a note Nov 1, 2012

Do we consider this v -> a change breaking since it means code like hyp0f1(v=0.3, z=0.3) no longer works?

@rgommers
SciPy member
rgommers added a note Nov 1, 2012

Good point, I should have changed the docstring instead. I'll change it back.

@rgommers
SciPy member
rgommers added a note Nov 1, 2012

done

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
@rgommers rgommers MAINT: change parameter name "a" back to "v" in special.hyp0f1.
This could have broken code like ``special.hyp0f1(v=1, x=2)``.
270e8cd
@pv
SciPy member
pv commented Nov 4, 2012

Merged in b05c922

@pv pv closed this Nov 4, 2012
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.