Vectorize StateDistributionCR.logpdf against scipy.special

logpdf() runs once per state on every E step over the whole lineage
array, but still went through the frozen rv_continuous machinery
(self.div_clock.logpdf/logsf, self.death_clock.logpdf/logsf) for that -
the same class of overhead already fixed in division_probability(),
just not applied here. Replaced with direct vectorized formulas against
scipy.special (gammaln/gammaincc), matching the old output to float
precision (verified numerically) while cutting test_works_through_crossval
from ~14.7s to ~6.3s in isolation.

The div_clock/death_clock properties are left in place since
compare_emissions.py still uses them outside this hot path.

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

Author: aarmey

lineage/states/StateDistributionCR.py

@@ -39,6 +39,16 @@
 TIME_FLOOR = 1e-10
 
 
+def _gamma_logpdf(t: np.ndarray, a: float, scale: float) -> np.ndarray:
+    """log Gamma(a, scale) pdf, vectorized directly against scipy.special."""
+    return (a - 1.0) * np.log(t) - t / scale - a * np.log(scale) - gammaln(a)
+
+
+def _gamma_logsf(t: np.ndarray, a: float, scale: float) -> np.ndarray:
+    """log Gamma(a, scale) survival function, vectorized directly against scipy.special."""
+    return np.log(gammaincc(a, t / scale))
+
+
 def event_masks(x: np.ndarray) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
     """Split cells into the three competing-risks cases.
 
@@ -242,19 +252,26 @@ def dof(self) -> int:
         return 3 if self.fixed_death_shape else 4
 
     def logpdf(self, x: np.ndarray) -> np.ndarray:
-        """Competing-risks log likelihood of each cell's phase observation."""
+        """Competing-risks log likelihood of each cell's phase observation.
+
+        Called once per state on every E step, over the whole lineage array, so it is
+        vectorized directly against ``scipy.special`` rather than going through the
+        frozen ``rv_continuous`` machinery (``self.div_clock.logpdf`` etc.), which pays
+        generic validation/broadcasting overhead on every call regardless of array size.
+        """
         divided, died, censored = event_masks(x)
         timed = divided | died | censored
         t = np.clip(x[:, 1], TIME_FLOOR, None)
 
-        div, death = self.div_clock, self.death_clock
+        a1, s1 = self.params[1], self.params[2]
+        a2, s2 = self.params[3], self.params[4]
 
         ll = np.zeros(x.shape[0])
         # Every timed cell survived both clocks up to t; the one that fired then
         # swaps its survival term for a density.
-        ll[timed] += div.logsf(t[timed]) + death.logsf(t[timed])
-        ll[divided] += div.logpdf(t[divided]) - div.logsf(t[divided])
-        ll[died] += death.logpdf(t[died]) - death.logsf(t[died])
+        ll[timed] += _gamma_logsf(t[timed], a1, s1) + _gamma_logsf(t[timed], a2, s2)
+        ll[divided] += _gamma_logpdf(t[divided], a1, s1) - _gamma_logsf(t[divided], a1, s1)
+        ll[died] += _gamma_logpdf(t[died], a2, s2) - _gamma_logsf(t[died], a2, s2)
 
         assert not np.any(np.isnan(ll))
         return ll