Fix #333 by filtering jags out of discrete search

Our discrete search was returning several tightly clustered points if
the jagged edge of a curve happened to move a function momentarily back
across a boundary before proceeding forward, so let’s filter out small
series of solutions that cluster close together.

This also introduces the very first unit tests of the search functions
themselves, which lets us catch edge cases without needing to test edge
cases separately for each almanac routine.

skyfield/searchlib.py

@@ -65,11 +65,17 @@ def _find_discrete(ts, jd, f, epsilon, num):
         # below epsilon at around the same time; so for efficiency we
         # only test the first pair.
         if ends[0] - starts[0] <= epsilon:
+            y = y.take(indices + 1)
+            # Keep only the first of several zero crossings that might
+            # possibly be separated by less than epsilon.
+            mask = concatenate(((True,), diff(ends) > 3.0 * epsilon))
+            ends = ends[mask]
+            y = y[mask]
             break
 
         jd = o(starts, start_mask).flatten() + o(ends, end_mask).flatten()
 
-    return ts.tt_jd(ends), y.take(indices + 1)
+    return ts.tt_jd(ends), y
 
 def _find_maxima(start_time, end_time, f, epsilon, num):
     ts = start_time.ts

skyfield/tests/test_almanac_searches.py

@@ -1,3 +1,66 @@
 """Low-level tests of the almanac search routines."""
 
+from numpy import where, sin
+
+from skyfield.api import load
+from skyfield.constants import tau
 from skyfield.searchlib import find_discrete, _find_maxima as find_maxima
+
+bump = 1e-5
+epsilon = 1e-10
+
+def make_t():
+    ts = load.timescale(builtin=True)
+    t0 = ts.tt_jd(0)
+    t1 = ts.tt_jd(1)
+    return t0, t1
+
+def make_f(offset):
+    """Make a sine wave from 0 to 1, with `t` offset by `offset` days."""
+    def f(t):
+        return sin((t.tt + offset) * tau) >= 0.0
+    f.rough_period = 1.0
+    return f
+
+def is_close(value, expected):
+    return (abs(value - expected) < epsilon).all()
+
+def test_find_discrete_near_left_edge():
+    t0, t1 = make_t()
+    f = make_f(-bump)  # cross zero barely past t0
+    t, y = find_discrete(t0, t1, f, epsilon)
+    assert is_close(t.tt, (bump, 0.5 + bump))
+    assert list(y) == [1, 0]
+
+def test_find_discrete_near_right_edge():
+    t0, t1 = make_t()
+    f = make_f(bump)  # cross zero almost at the end of the range
+    t, y = find_discrete(t0, t1, f, epsilon)
+    assert is_close(t.tt, (0.5 - bump, 1.0 - bump))
+    assert list(y) == [0, 1]
+
+def test_find_discrete_with_a_barely_detectable_jag_right_at_zero():
+    t0, t1 = make_t()
+    def f(t):
+        n = t.tt
+        n = where(n < 0.5, n + 3.1 * epsilon, n - 3.1 * epsilon)
+        return sin(n * tau) >= 0.0
+    f.rough_period = 1.0
+    t, y = find_discrete(t0, t1, f, epsilon)
+    assert is_close(t.tt, (0.5 - 3.1 * epsilon, 0.5, 0.5 + 3.1 * epsilon))
+    assert list(y) == [0, 1, 0]
+
+def test_find_discrete_with_a_sub_epsilon_jag_right_at_zero():
+    t0, t1 = make_t()
+    def f(t):
+        n = t.tt
+        n = where(n < 0.5, n + epsilon * 0.99, n - epsilon * 0.99)
+        return sin(n * tau) >= 0.0
+    f.rough_period = 1.0
+
+    # We hard-code num=12, just in case the default ever changes to
+    # another value that might not trigger the symptom.
+    t, y = find_discrete(t0, t1, f, epsilon, 12)
+
+    assert is_close(t.tt, (0.5,))
+    assert list(y) == [0]