ENH: add cyl_hankel_1_all#181
Open
JaRoSchm wants to merge 2 commits into
Open
Conversation
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Reference issue
Towards #50, #99 (comment), scipy/scipy#19079
(Probably) requires #93, #132 (which are both ready from my pov)
What does this implement/fix?
This adds the function
cyl_hankel_1_allwhich surfaces AMOS' functionality to compute sequential orders of Bessel functions in one go. This gives huge speed-ups (depending on the number of orders needed) as outlined in #50. The reason is that Bessel functions are recursively computed in AMOS and this functionality simply saves the intermediate results which are otherwise recomputed. The interface is inspired by thelegendre_*_allfunctions contained in xsf/scipy.Additional information
I plan to submit similar PRs for the other Bessel functions but first wanted to agree on the details to avoid duplicate work. The (probably unusual) choice of starting with the Hankel function of the first kind was simply motivated by my research where I use these the most.
The thing I failed to come up with is the correct way to create a float overload of the function. The reason is the
cyparameter which expects anmdspanof the correct length. However, I would expect this to be a float which can not be handed to thedoublecode. The original AMOS code contains float versions of all the functions. But currently I do not think, that porting them to xsf is worth the effort. So I am happy for help here (or the decision to leave out the float overload for now). Another option would be to handle this in the Python wrapper for scipy.The test code compares the vectorized version to the scalar ones for many different parameters. One important observation is, that AMOS returns a non-zero
ierrif the computation failed for the given order. In this case,cyl_hankel_1returns nan. However, for the vectorized version a non-zeroierris returned if the computation fails for any order. Because of this, all entries of the array are set to nan as AMOS does not return the order where the computation failed exactly (which probably makes sense as it needs correct starting values for the recursion). This required special handling in the test code.Additionally, I want to emphasize that I am not really familiar with C++. So everything you can find here is an interpolation of my Python and C knowledge, lots of web searches, and imitating other xsf code.
Concerning the inclusion in scipy I have a VERY rough outline of the code here. But I did not try to finish this, especially as scipy with the current xsf main branch does not seem to compile. Currently I would plan to include the
diff_nfunctionality of the existing*_allfunctions in the Python wrapper and not in xsf as outlined in the linked branch (if we decide to include it at all for now).AI Generation Disclosure
The test code was originally generated using a self-hosted (by my university) Kimi K-2.6 from a Python version of the basic idea which I created. Then I modified it, e.g., to include the nan-handling.