Get root secant converging to wrong solution

[AndreaMalgarini]

Problem description and minimal code This issue is located in the secant function within get_next_passes. It is not properly an error in the sense that the solution of the zero finding function get_root_secant given the initial bounds (which are guessed) is not able to find the zero but instead diverges leading to an unfeasible solution. What should be done is either change method, going for a more robust (such as Scipy's brentq as @rubendibattista suggested) or write checks for this case so to have an adaptive initial guess.

This particular issue can be reproduced with the satellite LEMUR-2-BROWNCOW [43125] with the TLEs reported here:

1 43125U 18004Q   18251.42128650 +.00001666 +00000-0 +73564-4 0  9991
2 43125 097.5269 314.3317 0010735 157.6344 202.5362 15.23132245036381
1 43125U 18004Q   18251.74976761 +.00001622 +00000-0 +71703-4 0  9995
2 43125 097.5269 314.6610 0010766 156.5384 203.6346 15.23133161036434
1 43125U 18004Q   18252.14394450  .00001595  00000-0  70556-4 0  9993
2 43125  97.5268 315.0561 0010805 155.2806 204.8950 15.23134406 36495
1 43125U 18004Q   18252.86660073 +.00001552 +00000-0 +68737-4 0  9990
2 43125 097.5268 315.7805 0010867 152.9961 207.1843 15.23136929036609

The coordinates for which passages are trying to be found are :

latitude=51.953319
longitude=-8.174163
altitude=0.05

The day for which the problem happens is 2018-09-08 with a time of propagation set to 72 hours and an horizon set to 5 degrees. The output is the following:

Exception: ('Satellite crased at time %s', numpy.datetime64('2108-01-24T02:35:08.798417'))

What can be noticed is the date, which is 2108. What happens is that the method get_root_secant fails to arrive at a feasible zero, hence, it diverges and goes in a search region, in terms of time, where the satellite is crashed.

The fact this was not noticed can be explained by its rarity. I discovered it performing a huge number of simulations; it happened once for a list of 30 satellites in 180 days propagation scenario. So it has a very low probability to happen. Still it should be better to fix it.

Possible solutions One of the proposed solution, a sort of adaptive secant method, can be formulated like this:

        def get_root_secant(fun, start, end, tol=0.01):
            """Secant method."""
            x_0 = end
            x_1 = start
            fx_0 = fun(end)
            fx_1 = fun(start)
            it = 0
            counter = 0
            it_max = 10
            min_toll = 10

            if abs(fx_0) < abs(fx_1):
                fx_0, fx_1 = fx_1, fx_0
                x_0, x_1 = x_1, x_0
            while abs(x_0 - x_1) > tol:
                x_n = x_1 - fx_1 * ((x_1 - x_0) / (fx_1 - fx_0))
                x_0, x_1 = x_1, x_n
                try:
                    fx_0, fx_1 = fx_1, fun(x_n)
                    it += 1
                except ValueError:
                    it == it_max

                if counter == 9:
                    raise ValueError('Secant method failed to converge after '
                                     '{} attempts'.format(counter+1))

                if it >= it_max:
                    it = 0
                    counter +=1
                    x_0, x_1 =  start + np.random.rand(), start - np.random.rand()
                    fx_0, fx_1 = fun(x_0), fun(x_1)
                    if abs(fx_0) < abs(fx_1):
                        fx_0, fx_1 = fx_1, fx_0
                        x_0, x_1 = x_1, x_0

            if x_1-end > min_toll or start- x_1 > min_toll:
                raise ValueError('Secant method converged to a value outside '
                                 'the initial bounds')

            return x_1

It is just an idea, quite raw, but I tested on different scenarios and it works.

A second possibility is to employ brentq like this:

        def get_root(fun, start, end, tol=0.01):
            """Root finding scheme"""
            x_0 = end
            x_1 = start
            fx_0 = fun(end)
            fx_1 = fun(start)
            if abs(fx_0) < abs(fx_1):
                fx_0, fx_1 = fx_1, fx_0
                x_0, x_1 = x_1, x_0

            x_n = optimize.brentq(fun, x_0, x_1)
            return x_n

This solution is quite cleaner and definitely more robust On the benchmark note it showed to be slightly slower than the previous.

Versions of Python 2.7.15rc1, package at hand and relevant dependencies pyorbital v1.3.1.

[mraspaud]

Thank you very much for reporting this ! The thorough investigation and proposed solution is very much appreciated! I believe this is related to what was reported in #22 too. I will have a closer look a bit later as I am on business trip until the end of the week.

If you have the time, a pull request would be great :) I don't have any particular preference for the method used, but if we can avoid depending in scipy, that would be nice. Not that scipy is a problem it itself, it's just adding a dependency to the package that troubles me.

[rdbisme]

This is the reference code for the brentq method in Scipy.

We're gonna try to see if we can easily implement it exploiting numpy vectorization.

If performances are not satisfying, @mraspaud what is you opinion about Cython wrapping? I think that at that moment would be easier to add scipy as dep.

PS: I work w/ @AndreaMalgarini :)

[djhoese]

In my opinion, I'd rather depend on scipy than replicate the functionality of the brentq function inside it. That is, if there are no significant modifications to the algorithm.

[kevincroissant]

@AndreaMalgarini and @rubendibattista thanks for your work on this. I've been encountering issues that I am pretty sure are related to this for months, but they've been rare enough and intermittent enough that they never bothered me enough to look for the root cause until now. If you would like me to test anything, let me know and I will put it through it's paces as well.

[rdbisme]

@kevincroissant Thanks for the kind words. @AndreaMalgarini already tested in several cases but having more tests is clearly appreciated.