Issue #136: Speed up root finding

docs/sitecustomize.py

@@ -42,4 +42,3 @@
 
 # Mock some other imports that might cause issues
 sys.modules["Basilisk"].__path__ = "not/a/real/path"
-sys.modules["chebpy"] = MagicMock()

src/bsk_rl/_finish_install.py

@@ -1,7 +1,4 @@
 import io
-import os
-import subprocess
-import sys
 import zipfile
 from pathlib import Path
 
@@ -11,28 +8,6 @@
 
 
 def pck_install():
-    if os.uname().sysname == "Darwin" and os.uname().machine == "arm64":
-        subprocess.check_call(
-            [
-                sys.executable,
-                "-m",
-                "pip",
-                "install",
-                "--no-deps",
-                "git+https://github.com/chebpy/chebpy.git",
-            ]
-        )
-    else:
-        subprocess.check_call(
-            [
-                sys.executable,
-                "-m",
-                "pip",
-                "install",
-                "git+https://github.com/chebpy/chebpy.git",
-            ]
-        )
-
     r = requests.get(
         "https://simplemaps.com/static/data/world-cities/basic/simplemaps_worldcities_basicv1.76.zip"
     )

src/bsk_rl/envs/general_satellite_tasking/scenario/satellites.py

@@ -14,10 +14,10 @@
         Simulator,
     )
 
-import chebpy
 import numpy as np
 from Basilisk.utilities import macros
 from gymnasium import spaces
+from scipy.optimize import minimize_scalar, root_scalar
 
 from bsk_rl.envs.general_satellite_tasking.scenario.data import (
     DataStore,
@@ -294,9 +294,8 @@ class AccessSatellite(Satellite):
     def __init__(
         self,
         *args,
-        generation_duration: float = 60 * 95,
+        generation_duration: float = 600.0,
         initial_generation_duration: Optional[float] = None,
-        access_dist_threshold: float = 4e6,
         **kwargs,
     ) -> None:
         """Construct an AccessSatellite.
@@ -307,15 +306,12 @@ def __init__(
                 `time_limit` unless the simulation is infinite. [s]
             initial_generation_duration: Duration to initially calculate imaging windows
                 [s]
-            access_dist_threshold: Distance bound [m] for evaluating imaging windows
-                more exactly. 4e6 will capture >10 elevation windows for a 500 km orbit.
             args: Passed through to Satellite constructor
             kwargs: Passed through to Satellite constructor
         """
         super().__init__(*args, **kwargs)
         self.generation_duration = generation_duration
         self.initial_generation_duration = initial_generation_duration
-        self.access_dist_threshold = access_dist_threshold
 
     def reset_pre_sim(self) -> None:
         """Reset satellite window calculations and lists."""
@@ -388,9 +384,15 @@ def calculate_additional_windows(self, duration: float) -> None:
         times = r_BP_P_interp.x[window_calc_span]
         positions = r_BP_P_interp.y[window_calc_span]
 
+        r_max = np.max(np.linalg.norm(positions, axis=-1))
+        access_dist_thresh_multiplier = 1.1
         for location in self.locations_for_access_checking:
+            alt_est = r_max - np.linalg.norm(location["location"])
+            access_dist_threshold = (
+                access_dist_thresh_multiplier * alt_est / np.sin(location["min_elev"])
+            )
             candidate_windows = self._find_candidate_windows(
-                location["location"], times, positions, self.access_dist_threshold
+                location["location"], times, positions, access_dist_threshold
             )
 
             for candidate_window in candidate_windows:
@@ -419,6 +421,7 @@ def _find_elevation_roots(
         location: np.ndarray,
         min_elev: float,
         window: tuple[float, float],
+        min_duration: float = 0.1,
     ):
         """Find times where the elevation is equal to the minimum elevation.
 
@@ -427,15 +430,27 @@ def _find_elevation_roots(
         """
 
         def root_fn(t):
-            return elevation(position_interp(t), location) - min_elev
+            return -(elevation(position_interp(t), location) - min_elev)
 
-        settings = chebpy.UserPreferences()
-        with settings:
-            settings.eps = 1e-6
-            settings.sortroots = True
-            roots = chebpy.chebfun(root_fn, window).roots()
+        elev_0, elev_1 = root_fn(window[0]), root_fn(window[1])
 
-        return roots
+        if elev_0 < 0 and elev_1 < 0:
+            logging.warning(
+                "initial_generation_duration is shorter than the maximum window length; some windows may be neglected."
+            )
+            return []
+        elif elev_0 < 0 or elev_1 < 0:
+            return [root_scalar(root_fn, bracket=window).root]
+        else:
+            res = minimize_scalar(root_fn, bracket=window, tol=1e-4)
+            if res.fun < 0:
+                window_mid = res.x
+                r_open = root_scalar(root_fn, bracket=(window[0], window_mid)).root
+                r_close = root_scalar(root_fn, bracket=(window_mid, window[1])).root
+                if r_close - r_open > min_duration:
+                    return [r_open, r_close]
+
+        return []
 
     @staticmethod
     def _find_candidate_windows(