Ias15 implentation

[johnwez1]

Implemented the ias15 algorithm in c with python wrappers, and added to relevant tests.
Some notes:

• Didnt implement the “begins to oscillate”
• Used "global" calculations.
• Didn't implement user warning for 12 iterations
• My implementation treats each step as an independant integration, with dynamic stepping in between. This part seems a bit silly but I thought at least I would give a working version at first. Otherwise, I would need to create an interpolation function to match the other integrators.
• Slight concern that it might get caught in infinite loop, in the while time remaining > 0 loop.

[codecov]

Codecov Report

Attention: Patch coverage is 98.12030% with 5 lines in your changes missing coverage. Please review.

Project coverage is 99.89%. Comparing base (7c5ad82) to head (4176e74).
Report is 34 commits behind head on main.

Files with missing lines	Patch %	Lines
galpy/util/wez_ias15.c	97.95%	5 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #643      +/-   ##
==========================================
- Coverage   99.91%   99.89%   -0.02%     
==========================================
  Files         200      202       +2     
  Lines       29269    29677     +408     
  Branches      564      607      +43     
==========================================
+ Hits        29243    29646     +403     
- Misses         26       31       +5     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[jobovy]

Thanks for this PR! I'll need a bit of time to work through this and familiarize myself better with IAS15 to be able to review this. But in the meantime, it would be good to fix some of the failing tests. It seems like some of the orbit integration tests fail for the new integrator (these are run on different OSs and for different Python versions, so lead to a bunch of failed jobs here). Some comments:

• One issue with the tests is that there's a bug in some of the error messages that leads to an error in the test (e.g., https://github.com/jobovy/galpy/actions/runs/9051119773/job/24867266058?pr=643#step:19:382). I could push a fix and you can rebase, but it might be easiest if you just fix this in this PR so we can be sure that it fully fixes the issue. The main issue is that pot should be used instead of p in the error message that fails
• The error message is only displayed because energy is not conserved well enough. So fixing the error message should allow you to see how big of an issue this is. To iterate on the tests, you might want to run them locally before pushing to GitHub. You can do this with

  pytest -v tests/test_orbit.py -k test_energy_jacobi_conservation -xs
  

  for the main orbit integration tests and similar for the other tests that fail:

  pytest -v tests/test_orbit.py -k test_integrate_negative_time -xs
  

  and

  pytest -v tests/test_orbit.py -k test_integrate_backwards -xs
  

[jobovy]

Also, to avoid the pre-commit commits here (which you should pull in before making any further changes), install pre-commit locally:

pip install pre-commit
pre-commit install

This will run the pre-commit hooks whenever you commit locally, automatically fix issues (like indentation and code style), and tell you when you need to fix something manually.

[jobovy]

Thanks for fixing the error message bug! Looks like there's some flakiness in the integrator, because the energy conservation tests pass for some Python/OS versions but fail for others (at least that's what it seems like from the test runs). It seems like the integrator sometimes just returns NaNs. Also the backwards integration tests seem to always be failing, maybe the integrators doesn't immediately work for backwards integration?

[johnwez1]

I made the PR after every test in test_orbit was passing on local, so yeah I guess something dependant on OS/Python version. I've recreated some issues with a fresh python 3.12 so hopefully that will lead to a fix. Anyway, leave it with me! I don't have a whole lot of free time at the moment to work on it but hopefully it shouldn't take too long.

[github-actions]

This PR has been automatically marked as stale because it has not had recent activity. Please close the PR if it is no longer relevant. If no further activity occurs, the PR may be closed. Thank you for your contributions.

[hannorein]

@dangcpham just pointed out this PR to me. A couple of comments:

• The pow(*,1/7) function call is most certainly OS/machine dependent and could explain the difference you see in local/CI runs. Have a look at the REBOUND source code for a simple and fast machine independent implementation.
• I understand why treating each timestep separately without taking data from previous timesteps into account is desirable. However, in that case the performance will be quite quite bad as you'll need many more predictor corrector loops.
• You might want to switch to the new time stepping criteria that @dangcpham came up with: https://arxiv.org/abs/2404.07160. It should perform better in all cases.

[jobovy]

Thanks for the helpful comments @hannorein!

I'll pin this PR so it doesn't get automatically closed, in case @johnwez1 or someone else wants to finish it soon.

[jobovy]

Thanks for continuing to push on this @johnwez1 ! The tests are not failing now, but they seem to be stuck and not finishing, so I think I'm going to cancel them to not waste too much compute. Do they run locally in a short amount of time? The problem seems to be consistent across Python versions and OS.

Also, I've noticed that sometimes the workflows start running automatically and sometimes I have to manually approve them. If they don't start running soon after you push, please ping me, because I might not be aware that I have to approve them.

[jobovy]

Looks like the orbit integration tests are still hanging, so I think I'll cancel them again to not waste compute. Do they finish locally for you in a short time? I assume it's that they're going into an infinite loop, is this possible in the code?

[johnwez1]

All tests are passing locally in reasonable time. An infinite loop may be occurring, so I guess I will have to investigate on a different machine.

[jobovy]

All tests are passing locally in reasonable time. An infinite loop may be occurring, so I guess I will have to investigate on a different machine.

Okay, that's strange! In the CI, it seems to happen across Python versions and OS (happening for Linux, Mac, and Windows), so it seems like a fundamental issue in the code.

It's very frustrating debugging issues that only show up in the CI, but not locally and I unfortunately don't have any great suggestions for how to do it.

One option is to use nektos/act to run workflows locally, which should work to debug the Linux builds and then maybe fixes the issue overall (https://nektosact.com/introduction.html); but personally I haven't had much luck with this in the past and it's somewhat complicated to set up.

Another option would be to enable Actions for your galpy fork (probably here, but maybe in your settings instead) and then branch to another local branch to experiment. That way the workflows here don't all run every time you push the code (and you don't have to wait for me to approve them) and you can have more fine-grained control over which workflows run. If you then figure out the issue, you can merge your new local branch back into your ias15 branch and the PR here would be updated. Personally, I would put a bunch of debug print statements in the code to see where it hangs (note that you then probably have to add -s to the pytest command in build.yml to see this output, because otherwise it gets captured and not shown by pytest.

[jobovy]

Hi @johnwez1,

Thanks for these latest changes and apologies for my slow response, things have been very busy at the end of the term. Looks like all the tests are passing now 🎉 ! I also briefly tested the code locally and things seem to work. On a simple test problem, the ias15_c integrator is about 6 times slower than dop853_c, but it is more precise by a few orders of magnitudes.

There are just a few things left to do to fully finish this PR:

• There are five new lines that aren't covered by the tests, all in galpy/util/wez_ias15.c. Three lines are related to being able to interrupt the code; to cover these you should add ias15_c to the integrators in the test_orbit_c_sigint_full test in test_orbit.py and in the other ...sigint... tests as well in that file. The other two lines that are not covered are in the double seventhroot function, specifically the lines in the while(a<1e-7 && isnormal(a)){ loop. I don't know whether you know a test case that would hit those lines, that would be ideal to add. If not, we should probably just be okay with slightly increasing the number of uncovered lines.
• Could you add the ias15_c integrator to the docstrings of the relevant methods in Orbits.py? You can just search for all mentions of dop853_c and add to those docstrings
• It would also be good to add to the documentation under doc/source/orbit.rst to mention the new integrator. The most minimal addition is to mention it in the "Fast orbit integration and available integrators" section. You might want to add a note on when using ias15_c would be appropriate; I imagine it's mostly for cases where very high precision is required.
• Please add an entry to the HISTORY.txt file.

I think that should be it and hopefully this will be pretty straightforward. Thanks for your contributions so far!