Address code review: remove test coupling, fix log(0) warning, clean up mask evaluation

Agent-Logs-Url: https://github.com/FlorianPfaff/PyRecEst/sessions/3f0bd299-8203-4bff-91e8-b8e62be085a8

Co-authored-by: FlorianPfaff <6773539+FlorianPfaff@users.noreply.github.com>

Author: Copilot

pyrecest/distributions/hypersphere_subset/complex_watson_distribution.py

@@ -72,8 +72,10 @@ def log_norm(D, kappa):
 
         # Asymptotic formula for high kappa
         # log C ~ log(2) + D*log(pi) + (1-D)*log(kappa) + kappa
+        # log_c_high is evaluated for all kappa before masking; clip to avoid log(0) warning
         log_c_high = (
-            math.log(2) + D * math.log(math.pi) + (1 - D) * np.log(kappa + 1e-300) + kappa
+            math.log(2) + D * math.log(math.pi)
+            + (1 - D) * np.log(np.maximum(kappa, 1e-300)) + kappa
         )
 
         # Intermediate formula (Mardia1999 Eq. 3):
@@ -322,13 +324,14 @@ def _sample_diagonal_complex_bingham_magnitudes(Lambda, D):
     Lambda_pos = Lambda[: D - 1]  # first D-1 (positive) eigenvalues
 
     # Precompute for the truncated exponential inverse CDF
-    temp1 = np.where(Lambda_pos >= 0.03, -1.0 / np.where(Lambda_pos >= 0.03, Lambda_pos, 1.0), 0.0)
-    temp2 = np.where(Lambda_pos >= 0.03, 1.0 - np.exp(-Lambda_pos), 0.0)
+    large = Lambda_pos >= 0.03
+    safe_lambda = np.where(large, Lambda_pos, 1.0)
+    temp1 = np.where(large, -1.0 / safe_lambda, 0.0)
+    temp2 = np.where(large, 1.0 - np.exp(-Lambda_pos), 0.0)
 
     s = np.zeros(D)
     while True:
         U = np.random.rand(D - 1)
-        large = Lambda_pos >= 0.03
         if np.any(large):
             s[: D - 1][large] = temp1[large] * np.log(1.0 - U[large] * temp2[large])
         if np.any(~large):

pyrecest/tests/distributions/test_complex_watson_distribution.py

@@ -6,9 +6,6 @@
 from pyrecest.distributions.hypersphere_subset.complex_watson_distribution import (
     ComplexWatsonDistribution,
 )
-from pyrecest.distributions.hypersphere_subset.bayesian_complex_watson_mixture_model import (
-    _simplex_integral,
-)
 
 
 def _random_unit_vector(D, rng=None):
@@ -149,24 +146,6 @@ def test_fit_recovers_parameters(self):
         self.assertAlmostEqual(dist_hat.kappa, kappa, delta=2.0)
 
 
-class TestSimplexIntegral(unittest.TestCase):
-    def test_D1(self):
-        self.assertAlmostEqual(_simplex_integral(np.array([3.0])), np.exp(3.0))
-
-    def test_D2_known(self):
-        # int_0^1 exp(a*t + b*(1-t)) dt = (exp(a) - exp(b)) / (a - b)
-        a, b = 2.0, 1.0
-        expected = (np.exp(a) - np.exp(b)) / (a - b)
-        result = _simplex_integral(np.array([a, b]))
-        self.assertAlmostEqual(result, expected, places=8)
-
-    def test_D3_nonnegative(self):
-        result = _simplex_integral(np.array([2.0, 1.0, 0.0]))
-        self.assertGreater(result, 0.0)
-        # Known value from direct integration: exp(2)/2 - exp(1) + 1/2
-        expected = np.exp(2) / 2 - np.exp(1) + 0.5
-        self.assertAlmostEqual(result, expected, places=5)
-
 
 if __name__ == "__main__":
     unittest.main()