ToO mode: ~40x faster Target runs (ephemeris grid, dedup, serial EAS below threshold)

src/nuspacesim/compute.py

@@ -187,12 +187,6 @@ def add_meta(self, name: str, value: Any, comment: str):
 
     logv(f"Running NuSpaceSim with Energy Spectrum ({config.simulation.spectrum})")
 
-    # The single EAS optical __call__ needs a process-based dask cluster whose
-    # ~2s spawn would otherwise be paid serially at that stage. Start it now, in
-    # the background, so the spawn overlaps the geometry/spectra/tau/decay
-    # stages; it's handed to eas() warm and torn down right after.
-    optical_cluster = BackgroundCluster() if config.detector.optical.enable else None
-
     logv("Computing [green] Geometries.[/]")
     beta_tr, thetaArr, pathLenArr, *_ = geom(
         config.simulation.thrown_events, store=sw, plot=to_plot
@@ -208,10 +202,19 @@ def add_meta(self, name: str, value: Any, comment: str):
         console.log(
             "\t[red] WARNING: No valid events thrown! Exiting early! Check geometry![/]"
         )
-        if optical_cluster is not None:
-            optical_cluster.close()
         return sim
 
+    # The single EAS optical __call__ can fan out over a process-based dask
+    # cluster, but its ~1s spawn/teardown only pays for itself on large valid
+    # batches. Decide now that the valid count is known; when a cluster is
+    # warranted, start it in the background so the spawn overlaps the
+    # spectra/tau/decay stages, then hand it to eas() warm and tear it down
+    # right after.
+    use_cluster = config.detector.optical.enable and (
+        beta_tr.size >= config.simulation.eas_parallel_threshold
+    )
+    optical_cluster = BackgroundCluster() if use_cluster else None
+
     init_lat, init_long = geom.find_lat_long_along_traj(np.zeros_like(beta_tr))
     sw(
         ("init_lat", "init_lon"),
@@ -243,14 +246,16 @@ def add_meta(self, name: str, value: Any, comment: str):
             init_lat,
             init_long,
             cloudf=cloud,
-            client=optical_cluster.client(),
+            client=optical_cluster.client() if use_cluster else None,
+            serial=not use_cluster,
             store=sw,
             plot=to_plot,
         )
 
         # Single consumer is done; release the warm cluster immediately.
-        optical_cluster.close()
-        optical_cluster = None
+        if use_cluster:
+            optical_cluster.close()
+            optical_cluster = None
 
         logv("Computing [green] Optical Monte Carlo Integral.[/]")
         mcint, mcintgeo, passEV, mcunc = geom.mcintegral(

src/nuspacesim/config.py

@@ -293,13 +293,21 @@ class TargetOfOpportunity(BaseModel):
         """Date of the event and format"""
         source_obst: float = 86400  # 24.0 * 60.0 * 60.0
         """Observation time (s). Default = 1 day"""
+        ephemeris_step: float = 60.0
+        """Grid spacing (s) for source/sun/moon sky positions, which are
+        evaluated exactly on the grid and cubic-interpolated to thrown times.
+        0 evaluates every thrown time exactly. Default = 60 s"""
 
     ################################################################################
 
     mode: Literal["Diffuse", "Target"] = "Diffuse"
     """ The Simulation Mode """
     thrown_events: int = 1000
     """ Number of thrown event trajectories. """
+    eas_parallel_threshold: int = 150_000
+    """Valid showers at or above which the EAS optical stage runs on a
+    process-based dask cluster; below it the stage runs in-process, avoiding
+    the cluster's ~1 s spawn/teardown. 0 always uses the cluster."""
     max_cherenkov_angle: Radians = np.radians(3)
     """ Maximum Cherenkov Angle (Radians). """
     max_azimuth_angle: Radians = np.radians(360)

src/nuspacesim/simulation/eas_optical/cphotang.py

@@ -863,6 +863,7 @@ def __call__(
         n_slant_sub=8,
         n_energy_low=3,
         n_energy_high=8,
+        serial=False,
     ):
         """
         Iterate over the list of events and return the result as pair of
@@ -900,6 +901,11 @@ def __call__(
         Gauss-Legendre quadrature node counts forwarded to :meth:`run` (defaults
         match ``run()``). The pipeline wires these from
         ``config.simulation.cherenkov_quadrature``.
+
+        ``serial=True`` runs the whole batch in this process through the same
+        kernel the workers use and ignores ``client``; the pipeline picks it
+        below ``config.simulation.eas_parallel_threshold`` showers, where a
+        cluster's spawn and teardown would outweigh the work.
         """
 
         if (
@@ -940,13 +946,18 @@ def chunk_worker(b, a, e, lat, lon):
             d_rows = d_batch.T if per_wavelength else d_batch[None, :]
             return np.concatenate([d_rows, c_batch[None, :]], axis=0)
 
-        results = map_showers_distributed(
-            chunk_worker,
-            (betaE, alt, Eshow100PeV, init_lat, init_long),
-            n_rows=n_den + 1,
-            chunks=chunks,
-            client=client,
-        )
+        if serial:
+            results = chunk_worker(
+                *(np.asarray(x) for x in (betaE, alt, Eshow100PeV, init_lat, init_long))
+            )
+        else:
+            results = map_showers_distributed(
+                chunk_worker,
+                (betaE, alt, Eshow100PeV, init_lat, init_long),
+                n_rows=n_den + 1,
+                chunks=chunks,
+                client=client,
+            )
 
         # Unpack (n_rows, N): density rows then the Cang row. Collapsed ->
         # (N,); per-wavelength -> (N, n_wl) (transpose back the n_den rows).

src/nuspacesim/simulation/eas_optical/eas.py

@@ -92,14 +92,16 @@ def __call__(
         *args,
         cloudf=None,
         client=None,
+        serial=False,
         **kwargs,
     ):
         """
         Electromagnetic Air Shower operation.
 
             ``client`` is an optional pre-built distributed client forwarded to
             :meth:`CphotAng.__call__`; see :class:`BackgroundCluster`. When ``None``
-            CphotAng spins up its own LocalCluster.
+            CphotAng spins up its own LocalCluster. ``serial=True`` skips the
+            cluster entirely and runs CphotAng in this process.
         """
 
         # Mask out-of-bounds events. Do not pass to CphotAng. Instead use
@@ -121,6 +123,7 @@ def __call__(
             init_long[mask],
             cloudf,
             client=client,
+            serial=serial,
             n_nodes=quad.n_nodes,
             n_slant_sub=quad.n_slant_sub,
             n_energy_low=quad.n_energy_low,

src/nuspacesim/simulation/geometry/too.py

@@ -35,6 +35,7 @@
 import astropy.time
 import astropy.units as u
 import numpy as np
+from scipy.interpolate import CubicSpline
 
 
 class ToOEvent:
@@ -56,6 +57,7 @@ def __init__(self, config):
         self.sourceDATE = self.config.simulation.target.source_date
         self.sourceDateFormat = self.config.simulation.target.source_date_format
         self.sourceOBSTime = self.config.simulation.target.source_obst
+        self.ephemeris_step = self.config.simulation.target.ephemeris_step
 
         self.eventtime = astropy.time.Time(
             self.sourceDATE, format=self.sourceDateFormat, scale="utc"
@@ -74,44 +76,118 @@ def __init__(self, config):
             height=self.detalt * 1000 * u.m,
         )
 
+    def detframe(self, time):
+        return astropy.coordinates.AltAz(obstime=time, location=self.detcords)
+
+    def _ephemeris_grid(self, time):
+        """Coarse time grid spanning ``time`` for interpolated sky positions.
+
+        Returns ``(grid_time, grid_x, x)`` with ``x`` the query abscissae in
+        seconds, or ``None`` when exact evaluation is requested or is no more
+        work than the grid itself.
+        """
+        step = self.ephemeris_step
+        if step <= 0 or time.isscalar:
+            return None
+        x = time.utc.unix
+        x0, x1 = x.min(), x.max()
+        # Pad one step each side so every query is interior to the spline.
+        n_grid = int(np.ceil((x1 - x0) / step)) + 3
+        if n_grid >= x.size:
+            return None
+        grid_x = (x0 - step) + step * np.arange(n_grid)
+        grid_time = astropy.time.Time(grid_x, format="unix", scale="utc")
+        return grid_time, grid_x, x
+
     def localcoords(self, time):
-        detframe = astropy.coordinates.AltAz(obstime=time, location=self.detcords)
-        return self.eventcoords.transform_to(detframe)
+        grid = self._ephemeris_grid(time)
+        if grid is None:
+            return self.eventcoords.transform_to(self.detframe(time))
+        grid_time, grid_x, x = grid
+        exact = self.eventcoords.transform_to(self.detframe(grid_time))
+        # Interpolate the unit vector, not (alt, az): azimuth wraps at 2pi.
+        cos_alt = np.cos(exact.alt.rad)
+        vec = np.stack(
+            [
+                cos_alt * np.cos(exact.az.rad),
+                cos_alt * np.sin(exact.az.rad),
+                np.sin(exact.alt.rad),
+            ]
+        )
+        vx, vy, vz = CubicSpline(grid_x, vec, axis=1)(x)
+        alt = np.arctan2(vz, np.hypot(vx, vy))
+        az = np.arctan2(vy, vx) % (2.0 * np.pi)
+        return astropy.coordinates.AltAz(
+            alt=alt * u.rad, az=az * u.rad, obstime=time, location=self.detcords
+        )
 
     def get_sun(self, time):
         sun_coord = astropy.coordinates.get_body("sun", time)
-        detframe = astropy.coordinates.AltAz(obstime=time, location=self.detcords)
-        return sun_coord.transform_to(detframe)
+        return sun_coord.transform_to(self.detframe(time))
 
     def get_moon(self, time):
         moon_coord = astropy.coordinates.get_body("moon", time)
-        detframe = astropy.coordinates.AltAz(obstime=time, location=self.detcords)
-        return moon_coord.transform_to(detframe)
+        return moon_coord.transform_to(self.detframe(time))
 
     @staticmethod
-    def moon_phase_angle(time: astropy.time.Time) -> float:
+    def phase_angle_from_bodies(sun, moon):
         """
-        Returns the moon phase angle in rad
+        Moon phase angle in rad from geocentric sun and moon coordinates
         0 -> full moon
         pi -> new moon
         """
-        sun = astropy.coordinates.get_body("sun", time)
-        moon = astropy.coordinates.get_body("moon", time)
         elongation = sun.separation(moon)
         return np.arctan2(
             sun.distance * np.sin(elongation),
             moon.distance - sun.distance * np.cos(elongation),
         )
 
+    @classmethod
+    def moon_phase_angle(cls, time: astropy.time.Time) -> float:
+        """
+        Returns the moon phase angle in rad
+        0 -> full moon
+        pi -> new moon
+        """
+        sun = astropy.coordinates.get_body("sun", time)
+        moon = astropy.coordinates.get_body("moon", time)
+        return cls.phase_angle_from_bodies(sun, moon)
+
+    def _sun_moon_state(self, time):
+        """Exact (sun altitude, moon altitude, moon phase angle) in rad."""
+        # One ephemeris lookup per body and one AltAz frame serve the altitude
+        # cuts and the phase angle; the lookups dominate ToO-mode runtime.
+        sun = astropy.coordinates.get_body("sun", time)
+        moon = astropy.coordinates.get_body("moon", time)
+        detframe = self.detframe(time)
+        return np.stack(
+            [
+                sun.transform_to(detframe).alt.rad,
+                moon.transform_to(detframe).alt.rad,
+                self.phase_angle_from_bodies(sun, moon).value,
+            ]
+        )
+
+    def sun_moon_state(self, time):
+        """(sun altitude, moon altitude, moon phase angle) in rad at ``time``.
+
+        Evaluated on the ephemeris grid and cubic-interpolated when that is
+        cheaper than evaluating every requested time.
+        """
+        grid = self._ephemeris_grid(time)
+        if grid is None:
+            return self._sun_moon_state(time)
+        grid_time, grid_x, x = grid
+        return CubicSpline(grid_x, self._sun_moon_state(grid_time), axis=1)(x)
+
     def sun_moon_cut(self, time: astropy.time.Time) -> bool:
         """
         Function to calculate the time during which sun and moon allow observation
         True -> observation possible
         False -> no observation posible
         """
-        sun_alt = self.get_sun(time).alt.rad < self.sun_alt_cut
-        moon_alt = self.get_moon(time).alt.rad < self.moon_alt_cut
-        moon_phase = self.moon_phase_angle(time).value > self.MoonMinPhaseAngleCut
-        moon_cut = np.logical_or(moon_phase, moon_alt)
-
-        return np.logical_and(sun_alt, moon_cut)
+        sun_alt, moon_alt, moon_phase = self.sun_moon_state(time)
+        moon_cut = np.logical_or(
+            moon_phase > self.MoonMinPhaseAngleCut, moon_alt < self.moon_alt_cut
+        )
+        return np.logical_and(sun_alt < self.sun_alt_cut, moon_cut)

test/core/test_config.py

@@ -309,6 +309,7 @@ def test_default_simulation():
     assert a.model_dump() == {
         "mode": "Diffuse",
         "thrown_events": 1000,
+        "eas_parallel_threshold": 150000,
         "max_cherenkov_angle": "3.0000000000000004 deg",
         "max_azimuth_angle": "360.0 deg",
         "angle_from_limb": "7.0 deg",
@@ -336,6 +337,7 @@ def test_default_simulation():
             "source_date": "2022-06-02T01:00:00",
             "source_date_format": "isot",
             "source_obst": 86400,
+            "ephemeris_step": 60.0,
         },
     }
 

test/simulation/geometry/test_too.py

@@ -1,3 +1,4 @@
+import astropy.units as u
 import numpy as np
 import pytest
 from astropy.time import Time
@@ -100,3 +101,37 @@ def test_sun_moon_cuts(too_event):
     darkmoon_date = "2022-11-29"
     dark = Time(darkmoon_date + "T22:00:00", format="isot", scale="utc")
     assert too_event.sun_moon_cut(dark)
+
+
+def test_ephemeris_grid_matches_exact(nss_config_event):
+    """Grid-interpolated sky positions must reproduce exact evaluation."""
+    import copy
+
+    exact_conf = copy.deepcopy(nss_config_event)
+    exact_conf.simulation.target.ephemeris_step = 0.0
+    exact = too.ToOEvent(exact_conf)
+    interp = too.ToOEvent(nss_config_event)
+    assert interp.ephemeris_step == 60.0
+
+    # Random times over one observation window: far more queries than grid points.
+    rng = np.random.default_rng(7)
+    times = exact.eventtime + rng.random(5000) * exact.sourceOBSTime * u.s
+
+    a, b = exact.localcoords(times), interp.localcoords(times)
+    sep = np.arccos(
+        np.clip(
+            np.sin(a.alt.rad) * np.sin(b.alt.rad)
+            + np.cos(a.alt.rad) * np.cos(b.alt.rad) * np.cos(a.az.rad - b.az.rad),
+            -1,
+            1,
+        )
+    )
+    assert sep.max() < 1e-6  # rad; ~0.2 arcsec
+
+    sa, sb = exact.sun_moon_state(times), interp.sun_moon_state(times)
+    assert np.abs(sa - sb).max() < 1e-6
+    assert np.array_equal(exact.sun_moon_cut(times), interp.sun_moon_cut(times))
+
+    # Scalar and few-point queries fall back to exact evaluation.
+    assert interp._ephemeris_grid(times[0]) is None
+    assert interp._ephemeris_grid(times[:3]) is None