Updated Mathieu function implementation#99
Conversation
main.c to test_mathieu.cpp in the tests dir. I ran clang-format on all .h files.
functions to replace duplicated code.
and that's what I have.
fix the failed builds in CI.
and lapack which are requred for Mathieu fcns. Also fix the same build problem when using pixi. Also rerun clang-format on all files.
as well as the entirety of scipy.
eigendecomposition on the full NxN recursion matrix, the new dstevd is for tridiagonal, symmetric matrices, so it runs significantly faster.
matrices (this didn't cause a bug because the problem was mostly in variable naming), updates to the choice of c, removal of the c=0 branch for the modified fcns, and others.
setting c. Minor comments added to make_matrix_tridiagonal.h.
col (m-2)/2 in the oe coeff fcn.
size scale with q somehow for better accuracy at high q. Also run clang-format on everything.
|
I figured I'd pop in here again and show some test results. All tests are on GitHub under https://github.com/brorson/ScipyMathieuTesting You can use them to check my impl if desired. Meanwhile here is the eigenvalue spectrum of the existing SciPy Mathieu implementation and my new one:
Here are the Wronskian errors RMS(W(Mx1, Mx2) - 2/pi) for Mc and Ms computed by the existing SciPy impl: And the new impl: |
|
Fixes #47 |
rgommers
left a comment
rgommers
left a comment
There was a problem hiding this comment.
Thanks for continuing to work on this @brorson! The algorithmic part of this looks pretty solid to me, and the accuracy results you posted on the issue look like a significant improvement over what we have today. @steppi may have more thoughts, but also looking at the review of the previous PR, this seems to be close to final.
There are quite a few issues with the C++ code that would need resolving, and I believe @steppi offered some help. I did a high-level review - my impression is that the size of the diff in this PR can be at least halved with some refactoring. Tests, public API, and LAPACK usage all need work to make this mergeable - I added high-level comments on those inline. None of that is super difficult, and some of the refactoring to reduce code duplication could even be done quickly with an LLM (it's straightforward enough that they'll probably get it first-time-right). It'd probably be better to leave any refactoring until after all tests pass though, to make it easy to verify that the refactoring didn't change anything. Hence, I'd suggest this sequence of steps:
- Rewrite the added tests to fit the
xsftest suite pattern, and ensure that test pass/fail works as expected. - Update or remove the existing tests as needed to make them pass - probably loosening tolerances or changing reference values for a bunch of them.
- At this point, the whole test suite should pass, at least with OpenBLAS on Linux (it won't build on all platforms or with all LAPACK libraries)
- Then, do all the refactors for code duplication, dead code, and controlling the public API better.
- At the end, deal with LAPACK linkage and build system integration (I may be able to help with this).
@steppi could you have a quick look, does that seem like the right plan?
| printf("===========================================\n"); | ||
| printf("At end of run, pass = %d, fail = %d\n", pass, fail); | ||
|
|
||
| return 0; |
There was a problem hiding this comment.
Returning 0 unconditionally looks like the tests will always appear to pass. Taking a step back: these tests look more like a local debug helper - they should be using TEST_CASE, CAPTURE and REQUIRE like the other test files, and no printf.
It looks like these tests need a complete overhaul; I think @steppi had offered some help with that: #60 (comment)
| cpp = meson.get_compiler('cpp') | ||
|
|
||
| # ======================================== | ||
| # BLAS/LAPACK Dependencies -- needed for Mathieu fcns. |
There was a problem hiding this comment.
Adding a LAPACK dependency to all of xsf won't work, that's way too invasive - and it also actually won't work in its current form, since this is missing handling of symbol mangling which is very much nontrivial to add. Instead, there are two ways of dealing with this that I can think of:
- An adaptor pattern. Use the
dstevdsignature, but with anxsf-specific name. Then make the consumer responsible for adding a header file to their build that defines that name to thedstevdin the LAPACK library they use in their build. - Avoid LAPACK completely, and use a handwritten algorithm to obtain the eigenvalues + eigenvectors for order
m(one more numerical routine, but seems feasible) and might end up faster.
There was a problem hiding this comment.
What I'm thinking here is thatxsf doesn't always need to contain standalone implementations of functions.
Along the lines of what I'd suggested here, scipy/scipy#24871 (comment), I think it can just supply numerical kernels necessary for implementing the functions that can be glued together on the SciPy/CuPy sides. I'd like to write the handwritten algorithm in a way that could be used with what I'm proposing here cupy/cupy#9839, but for now I think the LAPACK calls and full implementation can stay in SciPy.
| namespace mathieu { | ||
|
|
||
| //------------------------------------------------------ | ||
| int mathieu_coeffs_ee(int N, double q, int m, double *AA) { |
There was a problem hiding this comment.
There is a lot of code duplication in this file, should be possible to factor that out, which will also make the code and in particular the actual differences between these four functions easier to understand.
| return retcode; | ||
| } | ||
|
|
||
| // The peak Fourier coeff tracks m. Therefore, |
There was a problem hiding this comment.
There is a lot of code duplication for this Fourier summation. The size of this file could be cut in half it looks like by factoring that out.
|
|
||
| // Sort ww vector from lowest to highest | ||
| // quickSort(ww, 0, N-1); | ||
| // print_matrix(ww, N, 1); |
There was a problem hiding this comment.
Aside from this commented out code needing to be removed: this whole code block is essentially duplicated from the one higher up, can be factored out.
|
|
||
| // Min and max macros for scalars. | ||
| #define MIN(a, b) (((a) < (b)) ? (a) : (b)) | ||
| #define MAX(a, b) (((a) > (b)) ? (a) : (b)) |
There was a problem hiding this comment.
Can be replaced with std::min and std::max
| } | ||
| } // for (k=(N-1) ... | ||
|
|
||
| // I should do a sort before doing the sum |
There was a problem hiding this comment.
For this kind of comment, it'd be good to remove the ones that are no longer relevant, and for the ones that are relevant but can be deferred for future work, change them to // TODO: ... with a bit more context on why it would help to address the TODO.
| #include "mathieu/mathieu_coeffs.h" | ||
| #include "mathieu/mathieu_eigs.h" | ||
| #include "mathieu/mathieu_fcns.h" | ||
| #include "mathieu/matrix_utils.h" |
There was a problem hiding this comment.
There's a lot that leaks into the public API here. Using namespace detail would be a solution that is done in many places in xsf.
| namespace mathieu { | ||
|
|
||
| //================================================================== | ||
| double besselj(int k, double z) { |
There was a problem hiding this comment.
These wrappers should probably be removed, xsf::cyl_bessel_j & co can be called directly instead.
| #include <vector> | ||
|
|
||
| // In case of gcc use Float128 for more accuracy. | ||
| #ifndef _Float128 |
There was a problem hiding this comment.
This should use HAVE_FLOAT128, which is written to a config header as part of this PR exactly to enable the usage in this file. This ifndef as written will always be true, since _Float128 isn't a macro.
That's an impressive amount of work, and valuable to have such a paper as a reference. |
|
I looked through everything closely. Algorithmically I think things look solid. I think your plan is the right general idea @rgommers, except that I'd like to skip dealing with LAPACK in I will take things from here and work towards getting this ready to merge. I will break the work up into a series of small PRs against this branch. |
|
@brorson. I now plan to submit two PRs against |
|
@steppi, thanks for tagging me and nice to see progress on this front and that it's also helpful for others! Just be aware that you might stumble across cases, where these "vectorized order" variants give different results than the "single order" ones. These are probably due to bugs from the Fortran to C translation (#92, #93, #138 and at least one more I already have found locally). If you find such cases, please feel free to ping me. I think I managed to find a quite efficient way for locating these errors (by comparing to an automatic Fortran-C translation in https://github.com/JaRoSchm/pyamos/). Actually, I am in the process of slowly working through my own "test suite" for locating such cases (https://github.com/JaRoSchm/pyamos/blob/main/compare_vectorized_xsf.py). After getting this to pass I also planned to getting back to #50. |
5 participants
I am back. I withdrew my last pull request for my updated Mathieu function impl since I decided to polish the implementation further. Now I am back after making improvements to the impl and fixing some bugs. I received some help from Ho-Chul Shin, whose paper on computing the Mathieu functions appeared in the SIAM Journal in Dec 2025.
Besides improving my implementation I wrote a paper on the details of my implementation. I realize it's difficult to accept a big blob of code implementing a difficult special function without understanding what's in the code. I attempt to explain my implementation in this paper:
https://github.com/brorson/ScipyMathieuPaper
I maintain both Latex and .pdf copies on that GitHub page.
As part of my impl I subjected my functions to some tests to validate their correctness. Results of testing are presented in that paper.
Please reach out to me with any questions -- I am eager to help get this into the SciPy/xsf trunk.
Stuart Brorson
Northeatern University
Boston, Mass, USA