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

WIP:trying the new solvers #368

Merged
merged 7 commits into from May 16, 2018

Conversation

2 participants
@nikita-astronaut
Copy link
Contributor

nikita-astronaut commented May 6, 2018

Hello! This is the PR to show how changing the API would affect poliastro.

  1. The new API does not support callbacks as the previous one. Before it was the function called at every step, now it requires some event (something becomes zero). So for now I removed callbacks stuff.

  2. In test_propagate_to_date_has_proper_epoch I removed cowell, because maintaining the rtol=1e-10 makes run of this test fantastically long.

  3. In test_sample.py I had to mark two cowell tests xfail, because somehow the required precision 1e-7 is not achieved. Perhaps I'll manage to fix that at some point.

  4. In the new API, the jac parameter is only required by some solvers, and RK45 is not one of them, so it is now used anywhere. But this is a useful function for later.

  5. I could not manage to fix the precision problem by changing the integration method.

To my opinion, if we need to refactor and use the new API, we could say that the Cowell method is used only for perturbation stuff, and mean_motion is used for analytical gravitational orbits.

@nikita-astronaut nikita-astronaut changed the title trying the new solvers WIP:trying the new solvers May 6, 2018

@nikita-astronaut

This comment has been minimized.

Copy link
Contributor Author

nikita-astronaut commented May 6, 2018

Hello! I rolled back a bit. Seems like using method='Radau' resolves all the issues. I returned tests back and removed some xfails. Seems a lot better now.

@Juanlu001
Copy link
Member

Juanlu001 left a comment

I added some questions here, although I would prefer to see the tests results as well (and for that we need to fix the mypy error that appeared)

@@ -143,28 +143,6 @@ def test_apply_zero_maneuver_returns_equal_state():
ss.v.to(u.km / u.s).value)


def test_cowell_propagation_callback():

This comment has been minimized.

@Juanlu001

Juanlu001 May 6, 2018

Member

Great, this remained as an old way to store the results.

@@ -240,7 +218,6 @@ def test_propagate_to_date_has_proper_epoch():
@pytest.mark.filterwarnings('ignore::UserWarning')
@pytest.mark.parametrize('method', [
mean_motion,
pytest.param(cowell, marks=pytest.mark.xfail),

This comment has been minimized.

@Juanlu001

Juanlu001 May 6, 2018

Member

Why are you removing this case? Is it because of #367?

This comment has been minimized.

@nikita-astronaut

nikita-astronaut May 7, 2018

Author Contributor

No, it just runs very-very long with Cowell and new API. I will investigate. 100 years seem to be too much

This comment has been minimized.

@nikita-astronaut

nikita-astronaut May 7, 2018

Author Contributor

With Cowell and rtol=1e-10, this test runs for more than 5 minutes. But reducing rtol in Cowell will cause precision problems.


from scipy.integrate import ode
from scipy.integrate import ode, solve_ivp

This comment has been minimized.

@Juanlu001

Juanlu001 May 6, 2018

Member

If scipy.integrate.ode is not used, let's just remove it from the import

if rr.successful():
r, v = rr.y[:3], rr.y[3:]
f_with_ad = functools.partial(func_twobody, k=k, ad=ad, ad_kwargs=ad_kwargs)
jac_with_k = functools.partial(gravitational_jac, k=k)

This comment has been minimized.

@Juanlu001

Juanlu001 May 6, 2018

Member

I see that this is not being used. Convergence problems? The tests never ran because of this new MyPy error https://travis-ci.org/poliastro/poliastro/jobs/375477928#L854

This comment has been minimized.

@nikita-astronaut

nikita-astronaut May 7, 2018

Author Contributor

No, the reason is that this Jacobian is gravitational only. For other petrurbations I will add later, when everything works here.

So I need to a check whether ad is provided, and if no, use gravitational_jac.

r, v = rr.y[:3], rr.y[3:]
f_with_ad = functools.partial(func_twobody, k=k, ad=ad, ad_kwargs=ad_kwargs)
jac_with_k = functools.partial(gravitational_jac, k=k)
rr = solve_ivp(f_with_ad, (0, tof), u0, t_eval=[tof], rtol=rtol, method='Radau')

This comment has been minimized.

@Juanlu001

Juanlu001 May 6, 2018

Member

If we are not bringing back "dop853" (which would require some extra coding, since it's not available in the new API), is there a particular reason to use a stiff solver like 'Radau'? Does RK45 work in this case?

This comment has been minimized.

@nikita-astronaut

nikita-astronaut May 7, 2018

Author Contributor

RK45 doesn't get enough rtol in a couple of tests run in test_sample.py. But I don't know why using Radau helps, it looks like magic here


if rr.successful():
r, v = rr.y[:3], rr.y[3:]
f_with_ad = functools.partial(func_twobody, k=k, ad=ad, ad_kwargs=ad_kwargs)

This comment has been minimized.

@Juanlu001

Juanlu001 May 6, 2018

Member

Great to see the extra *args going away :) (see discussion that followed scipy/scipy#6326 (comment))

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented May 6, 2018

Sorry, I made the MyPy change myself :) Now you will have to git stash && git fetch origin && git merge origin/RK45vsDORPRI835 && git stash pop before you go on, otherwise you will have ugly error messages on push.

@codecov

This comment has been minimized.

Copy link

codecov bot commented May 6, 2018

Codecov Report

Merging #368 into master will increase coverage by 0.29%.
The diff coverage is 87.34%.

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #368      +/-   ##
==========================================
+ Coverage   82.45%   82.75%   +0.29%     
==========================================
  Files          31       33       +2     
  Lines        1573     1722     +149     
  Branches      124      137      +13     
==========================================
+ Hits         1297     1425     +128     
- Misses        245      258      +13     
- Partials       31       39       +8
Impacted Files Coverage Δ
src/poliastro/integrator_params.py 100% <100%> (ø)
src/poliastro/twobody/propagation.py 85.56% <56.25%> (-7.77%) ⬇️
src/poliastro/integrators.py 90.37% <90.37%> (ø)

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 135274b...f52d88a. Read the comment docs.

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented May 7, 2018

Also, it would be nice to have a precision comparison with the old solver. For example, using this snippet:

import numpy as np
from astropy import units as u

from poliastro.twobody import Orbit
from poliastro.bodies import Earth
from poliastro.twobody.propagation import cowell

np.set_printoptions(precision=20)

r0 = [-2384.46, 5729.01, 3050.46] * u.km
v0 = [-7.36138, -2.98997, 1.64354] * u.km / u.s

initial = Orbit.from_vectors(Earth, r0, v0)

def accel(t0, state, k):
    v_vec = state[3:]
    norm_v = (v_vec * v_vec).sum() ** .5
    return 1e-5 * v_vec / norm_v

final = initial.propagate(3 * u.day, method=cowell, ad=accel)
print(final.r)
print(final.v)

I get:

>>> print(final.r)
[ 13179.39566663877121754922 -13026.25123408228319021873
  -9852.66213692844394245185] km
>>> print(final.v)
[ 2.78170542314378943516  3.21596786944631274352  0.16327165546278937791] km / s
@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented May 7, 2018

Regarding the poor performance of Radau that we discussed privately, I found an interesting plot:

screenshot-2018-5-7 untitled - 3817ff090acce3164049f11182e827fc8bdb pdf

(from https://ieeexplore.ieee.org/document/5641330/)

So we are not the only ones to notice a slowdown when using an implicit scheme for a non-stiff problem.

@Juanlu001
Copy link
Member

Juanlu001 left a comment

I made some comments, but anyway I am very happy with the result. We can add numba at a later stage. Well done!

norm_v = (v_vec * v_vec).sum() ** .5
return 1e-6 * v_vec / norm_v

final = initial.propagate(3 * u.day, method=cowell, ad=accel)

This comment has been minimized.

@Juanlu001

Juanlu001 May 16, 2018

Member

This test has no assertions?

This comment has been minimized.

@nikita-astronaut

nikita-astronaut May 16, 2018

Author Contributor

yes, the idea of the test is that the integrator simply converges and does not fail) I have no data to compare to in this case. But perhaps I could make up an analytical solution to compare. Necessary?

assert_quantity_allclose(final.v, v_expected)


def test_integrator_converges_with_small_perturbations():

This comment has been minimized.

@Juanlu001

Juanlu001 May 16, 2018

Member

Could we try with 0 * v_vec / norm_v instead?

@@ -65,3 +65,37 @@ def test_atmospheric_drag():
r, v = cowell(orbit, tof, ad=atmospheric_drag, R=R, C_D=C_D, A=A, m=m, H0=H0, rho0=rho0)

assert_quantity_allclose(norm(r) - norm(r0), dr_expected, rtol=1e-2)


def test_new_integrator_works_with_small_perturbations():

This comment has been minimized.

@Juanlu001

Juanlu001 May 16, 2018

Member

No need to call it "new_integrator", just "cowell"


class DOP835DenseOutput(DenseOutput):
def __init__(self, t_old, t_new, y_old, interpolation):
super(DOP835DenseOutput, self).__init__(t_old, t_new)

This comment has been minimized.

@Juanlu001

Juanlu001 May 16, 2018

Member

In Python 3 you can use super().__init__(t_old, t_new), please change

self.rtol, self.atol = validate_tol(rtol, atol, self.n)
self.min_step_change = min_step_change
self.max_step_change = max_step_change
self.order = 7

This comment has been minimized.

@Juanlu001
@@ -0,0 +1,185 @@
from scipy.integrate import OdeSolver, DenseOutput

This comment has been minimized.

@Juanlu001

Juanlu001 May 16, 2018

Member

This whole module is just impressive. Good job!

This comment has been minimized.

@nikita-astronaut

nikita-astronaut May 16, 2018

Author Contributor

thanks) that was a great thing to me to recall how all these finite-difference suffisticated methods work

@@ -0,0 +1,100 @@
import numpy as np

A = []

This comment has been minimized.

@Juanlu001

Juanlu001 May 16, 2018

Member

Could we do

A = [
    np.array([5.26001519587677318785587544488e-2]),
    np.array([1.97250569845378994544595329183e-2, 5.91751709536136983633785987549e-2]),
    ...
]

instead? Same for the rest of the lists.

Also, is there a reason they have to be arrays? Can we just do

A = [
    [5.26001519587677318785587544488e-2],
    [1.97250569845378994544595329183e-2, 5.91751709536136983633785987549e-2],
    [2.95875854768068491816892993775e-2, 0, 8.87627564304205475450678981324e-2]
    ...
]

?

This comment has been minimized.

@nikita-astronaut

nikita-astronaut May 16, 2018

Author Contributor

The first question is <>, I will do that. No, list of lists won't do because I then use np.dot and .T routines

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented May 16, 2018

Also, before merging this:

  • Rebase the branch on current master: git fetch upstream && git rebase upstream/master && git push -f origin RK45vsDORPRI835
  • Write a brief wiki page explaining the options that we considered and their problems, so we don't forget

And then I'll hit the button. Let's see if we can accelerate this before the next release, but we will likely need numba/numba#2952 fixed for that.

nikita-astronaut and others added some commits May 6, 2018

Update for MyPy 0.600 without changing the code
See http://mypy-lang.blogspot.com.es/2018/05/mypy-0600-released.html

Complying with MyPy more strictly should be a separate effort,
so we are "temporarily" setting the new flag as explained in the
release notes.
norm_v = (v_vec * v_vec).sum() ** .5
return 0.0 * v_vec / norm_v

final = initial.propagate(3 * u.day, method=cowell, ad=accel)

This comment has been minimized.

@Juanlu001

Juanlu001 May 16, 2018

Member

What about checking that final is approximately equal to initial?

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented May 16, 2018

Aaaaaaaaaaaaaaand I'm merging! 🎉 Thanks @nikita-astronaut!

@Juanlu001 Juanlu001 merged commit aefb5b6 into poliastro:master May 16, 2018

5 checks passed

codeclimate Approved by Juan Luis Cano Rodríguez.
Details
codecov/patch 87.34% of diff hit (target 82.45%)
Details
codecov/project 82.75% (+0.29%) compared to 135274b
Details
continuous-integration/appveyor/pr AppVeyor build succeeded
Details
continuous-integration/travis-ci/pr The Travis CI build passed
Details

@wafflebot wafflebot bot removed the 2 - In Progress label May 16, 2018

@nikita-astronaut nikita-astronaut deleted the nikita-astronaut:RK45vsDORPRI835 branch May 16, 2018

@Juanlu001

This comment has been minimized.

Copy link
Member

Juanlu001 commented Sep 19, 2018

Opened discussion to contribute this code directly to SciPy: scipy/scipy#9290

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment