Add Hypothesis property-based tests for core tHMM invariants (#1010)

* Add Hypothesis property-based tests for core tHMM invariants

Covers transition-matrix stochasticity, CellVar.divide, the stationary
distribution, weighted Bernoulli/Gamma emission estimators, LineageTree
generation, and Viterbi/log-score consistency. One test is a pinned,
strict xfail documenting a real bug it found: StateDistribution.logpdf
can raise FloatingPointError (NaN from +inf/-inf) when an unobserved
(-1) Bernoulli sentinel coincides with an uncensored, exactly-zero
Gamma observation under a shape parameter < 1.

* Fix FloatingPointError crash in StateDistribution.logpdf

Exclude sentinel/negative Bernoulli and Gamma observations from the
Gamma and Bernoulli log-likelihood terms up front, instead of only
zeroing them out afterwards. Previously an unobserved (-1) Bernoulli
value on a row with an uncensored, exactly-zero Gamma observation and
a shape parameter < 1 (where the Gamma density diverges to +inf at
x=0) could produce a (+inf) + (-inf) = NaN, which crashes under
np.seterr(all="raise").

Un-xfails test_gamma_logpdf_never_produces_nan now that the bug it
documented is fixed.

Author: aarmey

lineage/states/StateDistributionGamma.py

@@ -62,17 +62,27 @@ def logpdf(self, x: np.ndarray) -> np.ndarray:
         """
         ll = np.zeros(x.shape[0])
 
+        # Negative Bernoulli/Gamma observations are sentinels for "unobserved"
+        # and must not contribute to the log-likelihood. Excluding them up
+        # front (rather than only zeroing them out afterwards) avoids ever
+        # combining a diverging Gamma density (e.g. +inf at x=0 for shape < 1)
+        # with a -inf sentinel Bernoulli log-likelihood into a NaN.
+        valid = (x[:, 0] >= 0) & (x[:, 1] >= 0)
+
         # Update uncensored Gamma
-        ll[x[:, 2] == 1] += sp.gamma.logpdf(x[x[:, 2] == 1, 1], a=self.params[1], scale=self.params[2])
+        uncensored = valid & (x[:, 2] == 1)
+        ll[uncensored] += sp.gamma.logpdf(x[uncensored, 1], a=self.params[1], scale=self.params[2])
 
         # Update censored Gamma
-        ll[x[:, 2] == 0] += sp.gamma.logsf(x[x[:, 2] == 0, 1], a=self.params[1], scale=self.params[2])
+        censored = valid & (x[:, 2] == 0)
+        ll[censored] += sp.gamma.logsf(x[censored, 1], a=self.params[1], scale=self.params[2])
 
         # Remove dead cells
         ll[x[:, 0] == 0] = 0.0
 
         # Update for observed Bernoulli
-        ll[np.isfinite(x[:, 0])] += sp.bernoulli.logpmf(x[np.isfinite(x[:, 0]), 0], self.params[0])
+        observed_bern = valid & np.isfinite(x[:, 0])
+        ll[observed_bern] += sp.bernoulli.logpmf(x[observed_bern, 0], self.params[0])
 
         # Log likelihood of negative values should be zero
         ll[x[:, 1] < 0] = 0.0