BUG: special.hyp2f1 now allows neg int c value #8110

FormerPhysicist · 2017-10-31T00:03:15Z

This PR fixes several issues associated with the complex z version of hyp2f1 and relating to negative integer values of c:

Previously, if c=0 and z is complex, hyp2f1 would attempt to calculate a value and would return nan+nan*j. This PR fixes it so that (inf + 0j) is returned. This is similar in behavior to what happens when c=0 and z is real (in which case hyp2f1 returns inf).
Now if z is complex and c=-m, this PR fixes hyp2f1 so that it returns the correctly computed value if and only if either a or b equals a negative integer -n where n<m. When c=-m and neither a nor b are negative integers with smaller magnitude than c, then hyp2f1 still returns (inf +0j) as expected.
The code in specfun_wrappers.c duplicated one of the safety checks that is already present in specfun.f (i.e., the check if z=1 and c-a-b>=0). I removed this code so that HYGFZ is now solely responsible for validating the values of a,b,c.

The original issue for this PR was: #7340

Closes #7340
Closes #8119

person142 · 2017-10-31T22:05:31Z

scipy/special/specfun/specfun.f

        X=DBLE(Z)
        Y=DIMAG(Z)
        EPS=1.0D-15
        ISFER=0
-        L0=C.EQ.INT(C).AND.C.LT.0.0D0
+        L0=C.EQ.INT(C).AND.C.LE.0.0D0


Hm... this problem was here already, but checking C.EQ.INT(C) is bad practice because INT(C) overflows before C does. The floor check in the C wrapper is better.

Would it be better to use a bigger sized integer in the fortran code? Something like IDINT? Or does the floor function itself have better overflow protections?

My only concern here is with duplicating code in both the wrapper file and the fortran file which seems like it would make the code base hard to follow and maintain.

But to be clear I can also re-produce my checks in the c wrapper code if that's what needs to be done...

Or does the floor function itself have better overflow protections?

The idea is that in C floor returns a double, so overflow is not an issue. The closest Fortran functions CEILING and FLOOR return integers. I don't know of a way to round in Fortran that returns REAL(8).

So I think this change is fine; the code was already doing the checks like that anyway.

Just for clarification, you think my change is fine. I do not need to revert the safety checks that I removed from specfun_wrappers.c

Just for clarification, you think my change is fine. I do not need to revert the safety checks that I removed from specfun_wrappers.c

Yes.

person142 · 2017-10-31T22:08:07Z

Could you add test cases? They should be here:

https://github.com/scipy/scipy/blob/master/scipy/special/tests/test_mpmath.py#L105

FormerPhysicist · 2017-10-31T23:50:01Z

I actually do have quite a few test cases inside of a txt file that I was going to add to the folder:
https://github.com/scipy/scipy/tree/master/scipy/special/tests/data/local

These test cases would be loaded into local.npz and then the test script test_data.py would run them.

But I haven't added them to this PR, because the file also includes a lot of test cases for my upcoming PR that fixes #8054 .

But I can either add a smaller txt file with just the c=-n test cases as described above, or I can add test cases to the file you linked. I'm fine either way.

person142 · 2017-11-01T22:39:16Z

But I can either add a smaller txt file with just the c=-n test cases as described above, or I can add test cases to the file you linked. I'm fine either way.

I think we try to only add to test_data.py when we can't reproduce results on-the-fly with Mpmath (if for example Mpmath takes too long or doesn't have a particular function). So knowing that, use your best judgement as to where it should go. But do add a test specifically for this PR.

If you do end up adding to test_data.py and use something like Mpmath to generate the data, add the script to special/_precompute.

FormerPhysicist · 2017-11-03T19:00:44Z

I've added test cases to test_mpmath.py.

I've also fixed the complex version of the bug fixed in PR #7981 as well as the new issue #8119

person142 · 2017-11-06T14:17:51Z

scipy/special/tests/test_mpmath.py

+
+    # Test both the complex and real versions of hyp2f1
+    for z in (0.5+0.5j,0.5):
+        curr_pts=[p+(z,) for p in pts]


Note that the Pep 8 checker is failing because of the lack of whitespace around the =. Remember to pip install pep8 and run pep8 scipy/special to check style.

In general, remember to include spaces between function arguments and around operators like +.

@person142 Okay cool! I did not know about Pep 8...

person142 · 2017-11-06T14:19:04Z

scipy/special/tests/test_mpmath.py

+        (1.23, -4, -4),
+        (2.12, 0, -3),
+        (-2, 0.89, -4)
+        ]


Less indentation on the ].

person142 · 2017-11-06T14:20:38Z

scipy/special/tests/test_mpmath.py

+        scipy_values=[sc.hyp2f1(*p) for p in curr_pts]
+        scipy_values=np.array(scipy_values, dtype=curr_type)
+
+        assert_allclose(scipy_values,mp_values,rtol=1.e-10)


Look at the uses of FuncData elsewhere in the module; probably best to use the same pattern here. (You will get better formatted error reporting, which is always nice.)

I tried doing that, but the problem was that FuncData requires a single array of a single type in order to evaluate the function. In the case of hyp2f1, I would need to pass in an array that contains both real numbers (a,b,c) and complex numbers (z,mpmath's value). (Note that I can't just pass in a single array of complex numbers since this array can't be passed to hyp2f1 since it doesn't match hyp2f1's signature of real,real,real,complex).

Now I can pass a single array of real numbers (basically the columns=a,b,c,z_real,z_imag,expected_real,expected_imag) and pass in a wrapper function to FuncData that converts this array into something that can be passed to hyp2f1. This is what I did originally did with my text file of 4000 or so test cases. But I was just afraid that people might think that was overkill for just checking 12 or so values....

So let me know what you prefer. I already have code that uses FuncData; it's just a matter of cut and paste...

I vote FuncData, whenever assert_allclose fails you have to go hack up the test to figure out what actually went wrong. The other option (probably more common in the tests) is to use the columns a, b, c, z, expected where you make a, b, and c complex and write a wrapper that does hyp2f1(a.real, b.real, c.real, z).

FormerPhysicist · 2017-11-06T21:01:08Z

I've fixed the PEP 8 style issues.

@person142 I can also use FuncData to run the test cases if that's what needs to be done; see my reply to your previous comment.

FormerPhysicist · 2017-11-06T21:52:14Z

@person142 Fixed tests so that they use FuncData

FormerPhysicist · 2017-11-08T20:18:09Z

Fixed a small bug that was failing one of the tests in the 'continuous-integration/travis-ci/pr...' check...

ghost · 2017-11-08T21:08:35Z

Failure due to #8136.

pv · 2018-04-08T18:29:43Z

scipy/special/specfun/specfun.f

@@ -6228,6 +6230,14 @@ SUBROUTINE HYGFZ(A,B,C,Z,ZHF,ISFER)
           CALL GAMMA2(1.0D0+A/2.0D0-B,G2)
           CALL GAMMA2(0.5D0+0.5D0*A,G3)
           ZHF=G0*G1/(G2*G3)
+        ELSE IF (L7.AND.((A.EQ.C).OR.(B.EQ.C))) THEN
+           IF (A.EQ.C) A=BB


That modifies the input argument --- probably should set AA=BB and use AA instead below.

Also, is the intent of if statements that a=c=negint or b=c=negint?
If yes, that does not seem to be guaranteed by the above.

…-c-fix

FormerPhysicist · 2018-08-19T23:50:58Z

Also, with respect to the comment about changing the values of A and B, note that the original code already changes the values of A and B as required. Before exiting, A and B are returned to their original values (as stored in AA and BB). My code changes follow this convention.

tylerjereddy · 2020-05-19T23:25:35Z

@FormerPhysicist @pv @person142 No activity here in almost two years--I will probably bump the milestone soon given branching in ~1 week unless you tell me otherwise.

tylerjereddy · 2020-11-11T22:29:22Z

@FormerPhysicist @person142 Have to bump the milestone again here soon?

h-vetinari · 2020-12-30T17:40:35Z

scipy/special/specfun/specfun.f

+        L7=(L3.AND.(DABS(A).LE.DABS(C)))
+     &       .OR.(L4.AND.(DABS(B).LE.DABS(C)))


@FormerPhysicist, you write [my emphasis]:

Now if z is complex and c=-m, this PR fixes hyp2f1 so that it returns the correctly computed value if and only if either a or b equals a negative integer -n where n<m.

Shouldn't it then be

Suggested change

L7=(L3.AND.(DABS(A).LE.DABS(C)))

& .OR.(L4.AND.(DABS(B).LE.DABS(C)))

L7=(L3.AND.(DABS(A).LT.DABS(C)))

& .OR.(L4.AND.(DABS(B).LT.DABS(C)))

@h-vetinari No, I think it should be LE (less than or equal).

When a=c (or b=c) and both are negative integers, then we need to let the method fall through and hit this section of code.

That part of the code implements formula 15.4.2 from Abramowitz and Stegun.

Otherwise, when a (or b) is a negative integer and abs(a) < abs(c), then the hypergeometric function reduces to a polynomial that is analytic everywhere. And again in this case we need to let the method fall through and run the code here

Note that the line you've highlighted is just for detecting infinities. I believe I got those conditions right, but feel free to review them again.

Thanks for the answer. I get that this was for the infinity-detection - what confused me was the relation to the statement in the OP that I quoted, which might be better reflected as (because the n=m part already exists):

Now if z is complex and c=-m, this PR fixes hyp2f1 so that it returns the correctly computed value ~~if and only~~ also if either a or b equals a negative integer -n where n<m.

tylerjereddy · 2021-05-13T02:08:40Z

Bumping to 1.8.0 for the usual reasons--if activity does suddenly pick up in the next few weeks feel free to ping me to reconsider the milestone.

tylerjereddy · 2021-11-30T02:32:59Z

Inactive for the last 6 months so I'll bump the milestone.

Also, @steppi is this PR still relevant with your related work in C/Cython?

steppi · 2021-11-30T03:41:57Z

Inactive for the last 6 months so I'll bump the milestone.

Also, @steppi is this PR still relevant with your related work in C/Cython?

This has been superseded by #14454.

BUG: special.hyp2f1 now allows neg int c value; see issue scipy#7340

b8a4ee7

person142 reviewed Oct 31, 2017

View reviewed changes

person142 added defect A clear bug or issue that prevents SciPy from being installed or used as expected scipy.special labels Nov 1, 2017

FormerPhysicist mentioned this pull request Nov 3, 2017

Bug: special.hyp2f1, wrong value returned when a=c=-m, where m is pos integer #8119

Open

FormerPhysicist added 2 commits November 3, 2017 11:55

BUG: special.hyp2f1, fixed issue scipy#8119

9af8ba1

TST: special.hyp2f1, added test cases for c=neg integer

b7d3076

FormerPhysicist changed the title ~~BUG: special.hyp2f1 now allows neg int c value; see issue #7340~~ BUG: special.hyp2f1 now allows neg int c value; closes #7340, #8119 Nov 3, 2017

person142 reviewed Nov 6, 2017

View reviewed changes

scipy/special/tests/test_mpmath.py Outdated

(1.23, -4, -4),

(2.12, 0, -3),

(-2, 0.89, -4)

]

Copy link

Member

person142 Nov 6, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Less indentation on the ].

person142 reviewed Nov 6, 2017

View reviewed changes

BUG: fixing some white space issues

6ec88ca

TST: use FuncData to do unit testing

adcde78

ilayn changed the title ~~BUG: special.hyp2f1 now allows neg int c value; closes #7340, #8119~~ BUG: special.hyp2f1 now allows neg int c value Nov 8, 2017

BUG: fixing small bug in code

3f468fb

BUG: no longer return inf when c-a-b<0 but hyp2f1 is a polynomial

41aec8a

pv reviewed Apr 8, 2018

View reviewed changes

FormerPhysicist added 2 commits August 18, 2018 15:14

BUG: fixing small bug in hygfz

4d2466f

Merge branch 'master' of https://github.com/scipy/scipy into negative…

5837dae

…-c-fix

tylerjereddy added this to the 1.5.0 milestone Jan 13, 2020

tylerjereddy modified the milestones: 1.5.0, 1.6.0 May 25, 2020

tylerjereddy modified the milestones: 1.6.0, 1.7.0 Nov 22, 2020

h-vetinari reviewed Dec 30, 2020

View reviewed changes

tylerjereddy modified the milestones: 1.7.0, 1.8.0 May 13, 2021

steppi mentioned this pull request Jun 19, 2021

MAINT: Begin translation of hyp2f1 for complex numbers into Cython. #14256

Closed

tylerjereddy modified the milestones: 1.8.0, 1.9.0 Nov 30, 2021

person142 closed this Nov 30, 2021

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

BUG: special.hyp2f1 now allows neg int c value #8110

BUG: special.hyp2f1 now allows neg int c value #8110

FormerPhysicist commented Oct 31, 2017 •

edited by ilayn

Loading

person142 Oct 31, 2017

FormerPhysicist Oct 31, 2017

FormerPhysicist Nov 1, 2017

person142 Nov 1, 2017

FormerPhysicist Nov 3, 2017

person142 Nov 3, 2017

person142 commented Oct 31, 2017

FormerPhysicist commented Oct 31, 2017

person142 commented Nov 1, 2017

FormerPhysicist commented Nov 3, 2017

person142 Nov 6, 2017

person142 Nov 6, 2017

FormerPhysicist Nov 6, 2017

person142 Nov 6, 2017

person142 Nov 6, 2017

FormerPhysicist Nov 6, 2017

person142 Nov 6, 2017

FormerPhysicist commented Nov 6, 2017

FormerPhysicist commented Nov 6, 2017

FormerPhysicist commented Nov 8, 2017

ghost commented Nov 8, 2017

pv Apr 8, 2018

FormerPhysicist commented Aug 19, 2018

tylerjereddy commented May 19, 2020

tylerjereddy commented Nov 11, 2020

h-vetinari Dec 30, 2020

FormerPhysicist Dec 31, 2020 •

edited

Loading

h-vetinari Dec 31, 2020

tylerjereddy commented May 13, 2021

tylerjereddy commented Nov 30, 2021

steppi commented Nov 30, 2021

		L7=(L3.AND.(DABS(A).LE.DABS(C)))
		& .OR.(L4.AND.(DABS(B).LE.DABS(C)))

BUG: special.hyp2f1 now allows neg int c value #8110

BUG: special.hyp2f1 now allows neg int c value #8110

Conversation

FormerPhysicist commented Oct 31, 2017 • edited by ilayn Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

person142 commented Oct 31, 2017

FormerPhysicist commented Oct 31, 2017

person142 commented Nov 1, 2017

FormerPhysicist commented Nov 3, 2017

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

FormerPhysicist commented Nov 6, 2017

FormerPhysicist commented Nov 6, 2017

FormerPhysicist commented Nov 8, 2017

ghost commented Nov 8, 2017

Choose a reason for hiding this comment

FormerPhysicist commented Aug 19, 2018

tylerjereddy commented May 19, 2020

tylerjereddy commented Nov 11, 2020

Choose a reason for hiding this comment

FormerPhysicist Dec 31, 2020 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tylerjereddy commented May 13, 2021

tylerjereddy commented Nov 30, 2021

steppi commented Nov 30, 2021

FormerPhysicist commented Oct 31, 2017 •

edited by ilayn

Loading

FormerPhysicist Dec 31, 2020 •

edited

Loading