-
-
Notifications
You must be signed in to change notification settings - Fork 28
ENH: Updated and refactored Mathieu function implementation (part 1). #144
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Open
steppi
wants to merge
105
commits into
scipy:main
Choose a base branch
from
steppi:brorson-mathieu-part1
base: main
Could not load branches
Branch not found: {{ refName }}
Loading
Could not load tags
Nothing to show
Loading
Are you sure you want to change the base?
Some commits from the old base branch may be removed from the timeline,
and old review comments may become outdated.
Open
Changes from all commits
Commits
Show all changes
105 commits
Select commit
Hold shift + click to select a range
0adc5a6
Update mathieu.h with improved docstrings and comments.
brorson 9f62cf0
Add new directory of Mathieu fcn impls.
brorson 204e09d
Make sure we return nan when error is encountered.
brorson b09c280
Remove non-standard d in double constant literals.
brorson 41c1f28
Merge branch 'scipy:main' into main
brorson 93a41df
Remove license file and place code under whatever license xsf uses.
brorson b8911c1
Merge branch 'main' of https://github.com/brorson/xsf_mathieu
brorson 5cbf45c
I removed my Makefile since it's not used in the Scipy CI tool. I moved
brorson cdf598b
Run clang-format on test_mathieu.cpp so it passes CI.
brorson cec5ca1
Merge branch 'scipy:main' into main
brorson 9c850da
Replaced malloc/free with std::vector, updated some comments, added
brorson 72a5146
Merge branch 'main' of https://github.com/brorson/xsf_mathieu
brorson b25ba42
Merge branch 'scipy:main' into main
brorson e32d63c
Remove mathieu.h since it doesn't belong in this dir.
brorson 7614c4e
run clang-format on all .h files.
brorson f96e616
Change dsyev -> dsyevd to speed up calc.
brorson 5ddc975
Change dsyevd to dstevd since dstevd is for tridiagonal matrices,
brorson 67314ef
Revert back to dsyevd. Use correct calling sig this time.
brorson c5e098b
Clean up files using clang-format before checking in.
brorson 7873a34
Merge branch 'scipy:main' into main
brorson 78ce632
I used a more recent version of clang-format. Hopefully this will
brorson f3074d5
Modify meson.build so it includes Lapack as a dependency.
brorson 536959f
Don't use _Float128 if not defined on this platform.
brorson eac37fe
Make _Float128 directive obey clang-format.
brorson 1d0ec45
Fix some constant doubles -- don't use d suffix since Mac doesn't lik…
brorson 7eb1c86
Run clang-format on mathieu_fcns to satisfy Scipy ci thing.
brorson d200bd4
Merge branch 'scipy:main' into main
brorson eb44cb3
Remove d suffix from constant doubles in test fcn to appease the CI t…
brorson 8a401b3
Fix meson build problem where standalone build of xsf doesn't find blas
brorson 678e532
Updates to the build files again so I can build both stand-alone xsf
brorson 2c0710f
Updates to mathieu.h and make_matrix.h based on reviewer comments.
brorson ce4baf2
Change impl to use dstevd instead of dsyevd. The old dsyevd did an
brorson cb4b6d6
This creates the tridiagonal recursion matrix suitable for dstevd.
brorson f08887b
Udpate the tests to match the new API for the dstevd impl.
brorson caab938
I ran clang-format on the .h files since I forgot to do this last time.
brorson 3c17ce4
Big changes including fixing row/col confusion in the eigenvector
brorson 36f71c9
Run clang-format on *.h files.
brorson 80d930d
Old s=0 code removed from mathieu_fcns. Also added s-max aproach to
brorson d1da899
Fix bug in se_2m+2 by changing coeff generation -- ask for matrix
brorson d6b2c3f
Add get_partial_sum_N to mathieu_fcns. The goal is to make the matrix
brorson bfd8ba5
Must update the mathieu test using clang-format.
brorson 78c3ef9
Fix bugs found by Ho-Chul Shin.
brorson b9df532
Merge branch 'main' into Brorson-Mathieu-impl
brorson a8cacdd
Merge branch 'main' into Brorson-Mathieu-impl
steppi f656268
TST: delete soon to be unnecessary mathieu tests
steppi 5b98901
MAINT: remove duplicate file
steppi a36e379
MAINT: Update make_matrix_tridiagonal to use mdspan
steppi 49ce54a
Move test_mathieu.cpp out of tests
steppi 9e5ec75
Get rid of include of now removed header
steppi 38650ea
TST: Add test for mathieu tridiagonal
steppi 666904a
Update docstring style comments
steppi a0cce79
keep old mathieu.h header and create new one so scipy still builds
steppi abfe40f
Use mathieu_new.h in debug program
steppi 4f6ac4c
Make type of q generic
steppi 08cc326
Merge branch 'main' into Brorson-Mathieu-impl-test-updates
steppi 319ccae
Add as_mdspan helper
steppi 6616001
Start reorganizing mathieu code
steppi adb3795
Use templates to reduce duplication in mathieu rr matrix
steppi 3dc49c7
Some fixes for mathieu stuff
steppi f0bc516
Add comment explaining what get_partial_sum_N does
steppi 0ceb45e
Add fourier series summing for mathieu + simplifications
steppi de819ac
Adding missing semicolons
steppi 75fff8a
Remove now unnecessary using statement
steppi d1df6fc
Fix some small mistakes
steppi bdbbe5c
Remove unused variable
steppi 765369c
Fix sign bug
steppi af53522
Add missing type
steppi 553d7ae
Compute fourier sum in degrees internally
steppi 2546fb8
Fix namespace for sinpi and cospi
steppi ede3b68
Add more comments on new mathieu implementations
steppi e60356a
Add AngleUnitPolicy
steppi ef6ac22
Add docstring like comment for sum_fourier_series
steppi 51fcdb8
Remove pre-(refactored + isolated) mathieu implementations
steppi 15e0022
Rename mathieu headers
steppi 005737b
Add comment explaining mathieu_legacy.h
steppi e82b4b9
Fix typo and minor error
steppi a447f5b
Mark wrap_PyUFunc_getfperr inline to work around ODR
steppi 2b8ddbb
Revert making wrap_PUFunc_getfperr inline
steppi cd6a277
Drop as_mdspan from numpy.h
steppi ffb73ce
Run clang format on numpy.h
steppi 117d9e5
Update mathieu tests
steppi 0ce474e
Add first sum_fourier_series test for mathieu
steppi ec7fb0c
Add another sum_fourier_series test
steppi 2d8cbd8
Add more sum_fourier_series tests
steppi fcc8b53
Add more fourier series mathieu tests
steppi eed6bf4
Add a header comment about mathieu tests in xsf
steppi 2829f04
Clean up some unnecessary things from #99
steppi f40b80b
Remove solve-group default for CuPy
steppi 203e303
Update pixi.lock
steppi 1d24a1b
Fix ufunc bug
steppi f4f416b
Run clang format
steppi eac5808
Merge branch 'ufunc-loop-bug' into brorson-mathieu-part1
steppi 9b8b852
Merge branch 'main' into brorson-mathieu-part1
steppi 9fdf040
Merge branch 'main' into brorson-mathieu-part1
steppi fa6756c
Merge branch 'main' into brorson-mathieu-part1
steppi 95e1000
Apply suggestions from code review
steppi 6ff3839
make capitalization consistent throughout
steppi b35edb8
use |q| in calculation of num Fourier terms
steppi 7faa939
fix normalization condition for sem and add comments
steppi 27f87a9
update comments on sign for mathieu
steppi a7e21f7
Apply suggestions from code review
steppi 64b52cb
fix more typos
steppi 652d508
rung clang format
steppi b51b24b
Apply suggestions from code review
steppi 788adeb
Merge branch 'main' into brorson-mathieu-part1
steppi File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
Large diffs are not rendered by default.
Oops, something went wrong.
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
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,230 @@ | ||
| /* Less accurate (but faster) legacy versions of Mathieu functions. */ | ||
|
|
||
| #pragma once | ||
|
|
||
| #include "error.h" | ||
| #include "specfun/specfun.h" | ||
|
|
||
| namespace xsf { | ||
|
|
||
| template <typename T> | ||
| T sem_cva(T m, T q); | ||
|
|
||
| template <typename T> | ||
| void sem(T m, T q, T x, T &csf, T &csd); | ||
|
|
||
| /* Mathieu functions */ | ||
| /* Characteristic values */ | ||
| template <typename T> | ||
| T cem_cva(T m, T q) { | ||
| int int_m, kd = 1; | ||
|
|
||
| if ((m < 0) || (m != floor(m))) { | ||
| set_error("mathieu_a", SF_ERROR_DOMAIN, NULL); | ||
| return std::numeric_limits<T>::quiet_NaN(); | ||
| } | ||
| int_m = (int)m; | ||
| if (q < 0) { | ||
| /* https://dlmf.nist.gov/28.2#E26 */ | ||
| if (int_m % 2 == 0) { | ||
| return cem_cva(m, -q); | ||
| } else { | ||
| return sem_cva(m, -q); | ||
| } | ||
| } | ||
|
|
||
| if (int_m % 2) { | ||
| kd = 2; | ||
| } | ||
| return specfun::cva2(kd, int_m, q); | ||
| } | ||
|
|
||
| template <typename T> | ||
| T sem_cva(T m, T q) { | ||
| int int_m, kd = 4; | ||
|
|
||
| if ((m <= 0) || (m != floor(m))) { | ||
| set_error("mathieu_b", SF_ERROR_DOMAIN, NULL); | ||
| return std::numeric_limits<T>::quiet_NaN(); | ||
| } | ||
| int_m = (int)m; | ||
| if (q < 0) { | ||
| /* https://dlmf.nist.gov/28.2#E26 */ | ||
| if (int_m % 2 == 0) { | ||
| return sem_cva(m, -q); | ||
| } else { | ||
| return cem_cva(m, -q); | ||
| } | ||
| } | ||
| if (int_m % 2) { | ||
| kd = 3; | ||
| } | ||
| return specfun::cva2(kd, int_m, q); | ||
| } | ||
|
|
||
| /* Mathieu functions */ | ||
| template <typename T> | ||
| void cem(T m, T q, T x, T &csf, T &csd) { | ||
| int int_m, kf = 1, sgn; | ||
| T f = 0.0, d = 0.0; | ||
| if ((m < 0) || (m != floor(m))) { | ||
| csf = std::numeric_limits<T>::quiet_NaN(); | ||
| csd = std::numeric_limits<T>::quiet_NaN(); | ||
| set_error("mathieu_cem", SF_ERROR_DOMAIN, NULL); | ||
| } else { | ||
| int_m = (int)m; | ||
| if (q < 0) { | ||
| /* https://dlmf.nist.gov/28.2#E34 */ | ||
| if (int_m % 2 == 0) { | ||
| sgn = ((int_m / 2) % 2 == 0) ? 1 : -1; | ||
| cem(m, -q, 90 - x, f, d); | ||
| csf = sgn * f; | ||
| csd = -sgn * d; | ||
|
|
||
| } else { | ||
| sgn = ((int_m / 2) % 2 == 0) ? 1 : -1; | ||
| sem(m, -q, 90 - x, f, d); | ||
| csf = sgn * f; | ||
| csd = -sgn * d; | ||
| } | ||
| } else { | ||
| using specfun::Status; | ||
| Status status = specfun::mtu0(kf, int_m, q, x, &csf, &csd); | ||
| if (status != Status::OK) { | ||
| csf = std::numeric_limits<T>::quiet_NaN(); | ||
| csd = std::numeric_limits<T>::quiet_NaN(); | ||
| sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; | ||
| set_error("mathieu_cem", sf_error, NULL); | ||
| } | ||
| } | ||
| } | ||
| } | ||
|
|
||
| template <typename T> | ||
| void sem(T m, T q, T x, T &csf, T &csd) { | ||
| int int_m, kf = 2, sgn; | ||
| T f = 0.0, d = 0.0; | ||
| if ((m < 0) || (m != floor(m))) { | ||
| csf = std::numeric_limits<T>::quiet_NaN(); | ||
| csd = std::numeric_limits<T>::quiet_NaN(); | ||
| set_error("mathieu_sem", SF_ERROR_DOMAIN, NULL); | ||
| } else { | ||
| int_m = (int)m; | ||
| if (int_m == 0) { | ||
| csf = 0; | ||
| csd = 0; | ||
| } else if (q < 0) { | ||
| /* https://dlmf.nist.gov/28.2#E34 */ | ||
| if (int_m % 2 == 0) { | ||
| sgn = ((int_m / 2) % 2 == 0) ? -1 : 1; | ||
| sem(m, -q, 90 - x, f, d); | ||
| csf = sgn * f; | ||
| csd = -sgn * d; | ||
| } else { | ||
| sgn = ((int_m / 2) % 2 == 0) ? 1 : -1; | ||
| cem(m, -q, 90 - x, f, d); | ||
| csf = sgn * f; | ||
| csd = -sgn * d; | ||
| } | ||
| } else { | ||
| using specfun::Status; | ||
| Status status = specfun::mtu0(kf, int_m, q, x, &csf, &csd); | ||
| if (status != Status::OK) { | ||
| csf = std::numeric_limits<T>::quiet_NaN(); | ||
| csd = std::numeric_limits<T>::quiet_NaN(); | ||
| sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; | ||
| set_error("mathieu_sem", sf_error, NULL); | ||
| } | ||
| } | ||
| } | ||
| } | ||
|
|
||
| template <typename T> | ||
| void mcm1(T m, T q, T x, T &f1r, T &d1r) { | ||
| int int_m, kf = 1, kc = 1; | ||
| T f2r = 0.0, d2r = 0.0; | ||
|
|
||
| if ((m < 0) || (m != floor(m)) || (q < 0)) { | ||
| f1r = std::numeric_limits<T>::quiet_NaN(); | ||
| d1r = std::numeric_limits<T>::quiet_NaN(); | ||
| set_error("mathieu_modcem1", SF_ERROR_DOMAIN, NULL); | ||
| } else { | ||
| using specfun::Status; | ||
| int_m = (int)m; | ||
| Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r); | ||
| if (status != Status::OK) { | ||
| f1r = std::numeric_limits<T>::quiet_NaN(); | ||
| d1r = std::numeric_limits<T>::quiet_NaN(); | ||
| sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; | ||
| set_error("mathieu_modcem1", sf_error, NULL); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| template <typename T> | ||
| void msm1(T m, T q, T x, T &f1r, T &d1r) { | ||
| int int_m, kf = 2, kc = 1; | ||
| T f2r = 0.0, d2r = 0.0; | ||
|
|
||
| if ((m < 1) || (m != floor(m)) || (q < 0)) { | ||
| f1r = std::numeric_limits<T>::quiet_NaN(); | ||
| d1r = std::numeric_limits<T>::quiet_NaN(); | ||
| set_error("mathieu_modsem1", SF_ERROR_DOMAIN, NULL); | ||
| } else { | ||
| using specfun::Status; | ||
| int_m = (int)m; | ||
| Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r); | ||
| if (status != Status::OK) { | ||
| f1r = std::numeric_limits<T>::quiet_NaN(); | ||
| d1r = std::numeric_limits<T>::quiet_NaN(); | ||
| sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; | ||
| set_error("mathieu_modsem1", sf_error, NULL); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| template <typename T> | ||
| void mcm2(T m, T q, T x, T &f2r, T &d2r) { | ||
| int int_m, kf = 1, kc = 2; | ||
| T f1r = 0.0, d1r = 0.0; | ||
|
|
||
| if ((m < 0) || (m != floor(m)) || (q < 0)) { | ||
| f2r = std::numeric_limits<T>::quiet_NaN(); | ||
| d2r = std::numeric_limits<T>::quiet_NaN(); | ||
| set_error("mathieu_modcem2", SF_ERROR_DOMAIN, NULL); | ||
| } else { | ||
| using specfun::Status; | ||
| int_m = (int)m; | ||
| Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r); | ||
| if (status != Status::OK) { | ||
| f2r = std::numeric_limits<T>::quiet_NaN(); | ||
| d2r = std::numeric_limits<T>::quiet_NaN(); | ||
| sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; | ||
| set_error("mathieu_modcem2", sf_error, NULL); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| template <typename T> | ||
| void msm2(T m, T q, T x, T &f2r, T &d2r) { | ||
| int int_m, kf = 2, kc = 2; | ||
| T f1r = 0.0, d1r = 0.0; | ||
|
|
||
| if ((m < 1) || (m != floor(m)) || (q < 0)) { | ||
| f2r = std::numeric_limits<T>::quiet_NaN(); | ||
| d2r = std::numeric_limits<T>::quiet_NaN(); | ||
| set_error("mathieu_modsem2", SF_ERROR_DOMAIN, NULL); | ||
| } else { | ||
| using specfun::Status; | ||
| int_m = (int)m; | ||
| Status status = specfun::mtu12(kf, kc, int_m, q, x, &f1r, &d1r, &f2r, &d2r); | ||
| if (status != Status::OK) { | ||
| f2r = std::numeric_limits<T>::quiet_NaN(); | ||
| d2r = std::numeric_limits<T>::quiet_NaN(); | ||
| sf_error_t sf_error = status == Status::NoMemory ? SF_ERROR_MEMORY : SF_ERROR_OTHER; | ||
| set_error("mathieu_modsem2", sf_error, NULL); | ||
| } | ||
| } | ||
| } | ||
|
|
||
| } // namespace xsf | ||
This file was deleted.
Oops, something went wrong.
This file was deleted.
Oops, something went wrong.
This file was deleted.
Oops, something went wrong.
This file was deleted.
Oops, something went wrong.
Oops, something went wrong.
Oops, something went wrong.
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.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
what is the plan for this file? what about backward-compatibility?
should we not keep old tests for smoke coverage?
Uh oh!
There was an error while loading. Please reload this page.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This file remains for now because it's still being used for the modified (radial) mathieu functions, and as a large order fallback for
mathieu_cv. Once the modified mathieu functions are refactored like this, and we figure out something else for the fallback, then we can get rid of the file. We can't really keep the old tests because they rely on computing the functions from end to end andxsfdoesn't offer that capability. The full implementations for CPU live in SciPy and rely on kernels from here that the CuPy implementations will also use. All of those old test cases (using degrees as angular unit) remain in SciPy.