ENH:MAINT:special:Rewrite amos F77 code #19587
Conversation
07223f2 to
4817951
Compare
|
I managed to make some progress in a quite unexpected direction. Because the code is so unorthodox, doing real arithmetic with complex ops, jumping from one for loop into the middle of another and feeling no shame and so on, Cythonization didn't cut it so I switched to C as another goto infested lang. Anyways, I quickly realized why I stopped coding in these langs because it's been roughly 30 years and still nobody managed to agree on what a complex number type should be. So probably at some point In fact quoting Amos himself ;
If only he knew... I unrolled as many GOTOs as I could. I'll clean the cython code when I'm done 😢 If you have objections to style or strategy please tell me as I hardly know what I am doing. |
|
165e716 finished the rewrite. I am pretty confident about I'll test the the remaining Bessel ones and start dealing with the wrappers. |
|
Tested on Win10 and Linux Mint. Only one test is barely failing rel err I would like to know what we should do to define missing macros in different platforms like Happy new year all 🎉 🎆 |
|
Status update;
If the first one is resolved and second item is acceptable, then this is ready to be reviewed. |
I'm just guessing but is numpy/numpy#24085 related? |
|
Ah yes, probably. Let me use the |
|
Ok seems like we are game. All feedback welcome. |
|
Impressive effort @ilayn! This isn't quite reviewable line by line, so I guess we should look at test coverage, build warnings, etc. It maybe also be useful to run clang-format over the new code now. |
There was a problem hiding this comment.
It was quite a feat to untangle and translate all of this, the Fortran is not at all easy to understand. Bessel functions can be tricky to compute and despite its tortured code style, AMOS is actually pretty solid, so I wanted to give this a very close look.
What I decided to do is pick one of the most difficult functions and study the code in depth, comparing to the original Fortran and Donald Amos' paper. I chose amos_binu which computes the modified Bessel function of the first kind in the right half plane. If this part had been translated perfectly, I would have had more confidence in everything, since it is among the most difficult parts to get right.
I found an incorrect condition which causes degraded performance in some cases. Just the kind of little mistake I could have made too if I was the one who translated this. It makes me think I'll need to look over everything with a similar level of care. I don't quite have bandwidth for that at the moment. It would probably be several days of work.
@ilayn, I don't know what your bandwidth looks like, but you might want to trace through all of the control flow again to check if you made any other little mistakes like that. Also, I found a comment in the wrong place, it'll be helpful if you make sure all of those are where they belong.
I found it helpful to look through Amos' paper I linked to above. Despite all of the spaghetti, most of what's happening here is that different methods are being used in different regions of
|
Ha, it's not less feat to go through this so thanks @steppi Funny that you picked So we are looking at this flow scipy/scipy/special/amos/zbinu.f Lines 34 to 47 in 87f8428 So first jumps, then falls through then jumps further and if not then falls through. The reason I have those extra conditions are these if else sillyness. If it falls through we finish off with 120 or 130. But if not further things happen. So the idea was whether it is coming from the 50 jump or coming from 40 with false IF condition. Fortunately we don't have another occurence of this so I am pretty confident about the remaining parts since I tested quite a number of them. I'll check what the difference is between your version and what I wrote |
|
OK apparently I fixed this in my local work but did not pass it over to the PR This binu is such an annoying piece of code. |
Yeah, it's the trickiest one in my opinion. If that had been perfect, I think I would have just gave everything else a spot check and merged. I'm feeling better about things again after your explanations. I should be able to find time to finish reviewing within the coming week. |
|
No worries, that's really good catch. Let's hope I didn't make more labor intensive mistakes. |
|
@steppi apologies for getting back to this again but I think you found a bug in the fortran code itself. Because when I check it with the fortran code itself this is what I get; this is the fortran code annotated with lots of print statements and gives 1e-13 difference This is what you get from the wrong SciPy result above. The reason I am convinced that this is a Fortran bug is because my code also goes into the same branches and not into UOIK as you mentioned it should go immediately. They do go indeed but return the same result. I guess copied the Fortran code a bit too good. |
|
Hmm I am more confused now because then my fixed version is diverged from the Fortran code itself. And if I add your fix then we have many tests failing. |
|
OK the current SciPy which uses Fortran also does not pass 1e-15 rtol level too C:\Users\ilhan> ipython
Python 3.10.6 (tags/v3.10.6:9c7b4bd, Aug 1 2022, 21:53:49) [MSC v.1932 64 bit (AMD64)]
Type 'copyright', 'credits' or 'license' for more information
IPython 8.19.0 -- An enhanced Interactive Python. Type '?' for help.
In [1]: import scipy
In [2]: scipy.__version__
Out[2]: '1.12.0rc1'
In [3]: import scipy.special as ssp
In [4]: ssp.iv(99.2, 15+15j)
Out[4]: (-2.0645908843310788e-55-1.1127773731040388e-55j)
In [11]: assert_allclose(ssp.iv(99.2, 15+15j), mpmath,rtol=1-15)
[..snip..]
AssertionError:
Not equal to tolerance rtol=-14, atol=0
Mismatched elements: 1 / 1 (100%)
Max absolute difference: 1.28411573e-68
Max relative difference: 5.47508521e-14
x: array(-2.064591e-55-1.112777e-55j)
y: array(-2.064591e-55-1.112777e-55j)So I think I'll just take a break for today and tomorrow I'll do another round. |
* MAINT: special: remove goto in acai * MAINT: special: remove some gotos from bknu * MAINT: simplify if statements * MAINT: special: remove a few gotos unk1 * MAINT: remove duplicate code in acai * STY:special:amos:Basic fixups for do loops --------- Co-authored-by: Ilhan Polat <ilhanpolat@gmail.com>
MAINT:special:amos: Add const specifier for array data STY:special:amos: Cleanup and organize data and comments MAINT:special:amos: Silence compiler for fabs(integer) use
|
I am not sure where I confused myself but I used my working code and swapped places following your lead and indeed tests pass. And moreover, I reverted the modified test since that also passes too. Let's see how it goes. Overall, I think this indeed was a very good catch. And the test suite is quite comprehensive it seems. And just for the ease of testing, you can choose to test only the public functions as they go through multiple hoops and different domains should catch any remaining discrepencies. |
There was a problem hiding this comment.
Thanks @ilayn. I've looked over everything else with a level of attention somewhere between spot check and exhaustive study while referring to the original paper. The rest looks good to me, and none of it has such twisted logic in the original as binu. It's possible I missed something still, but I feel confident now, especially with the level of testing we have.
|
Thanks @steppi much appreciated. |
|
@steppi do you think we should update the docstrings as currently they state along the lines of: |
Sorry. I didn't find any sign of ilayn#7 mentioned in the comment chain above. I would have waited to merge this if I'd known you had this concern. Honestly, I thought the docstrings were fine as is. I was thinking of the AMOS routines in the abstract, not tied to any particular implementation. Now that you mention it though, letting users know these are translations might save people from some confusion around our functions giving slightly different results from AMOS in some cases. For the future. I expect these functions to diverge further and further from AMOS over time as we make improvements. If so, I think the docstrings should be updated again to say our implementations are based on those of Amos, citing D.E. Amos' Computation of Bessel functions of a complex argument, which discusses the math behind the implementations. |
No problem at all. There is just too much to look at with these Fortran stuff so we are always forgetting stuff. But indeed I would also using the abstract terminology of citing Amos' paper. The interested readers would discover where we got the implementations at the top of the license comments anyways. |
Reference issue
Part of #18566
Relevant issues will be added
What does this implement/fix?
The rewrite of the AMOS library.
Additional information
Not as big as cdflib, should be tested meticulously for the edge cases.
Total gotos left: ~40 (out of too many)