Added checks for negative RRR counts

following from the multipole decomposition

Author: misharash

RascalC/correction_function_3pcf.py

@@ -31,6 +31,19 @@ def compute_inv_phi_aperiodic_3pcf(n: int, m: int, n_multipoles: int, r_bins: np
     for ell in range(n_multipoles):
         # (NB: we've absorbed a factor of delta_mu into RRR_true here)
         leg_triple[:, :, ell] += (2.*ell+1.) * np.sum(legendre(ell)(mu_cen)[None, None, :] * RRR_true, axis=-1)
+    
+    # as a precaution, check for negative counts, which should be problematic
+    n_mu_check = 2001
+    triple_counts_check = np.zeros([n, n, n_mu_check])
+    mu_values_check = np.linspace(-1, 1, n_mu_check)
+    for ell in range(n_multipoles):
+        triple_counts_check += leg_triple[:, :, ell][:, :, None] * legendre(ell)(mu_values_check)[None, None, :]
+    problem_indices = np.argwhere(triple_counts_check <= 0)
+    if len(problem_indices) > 0:
+        rbins, mu_counts = np.unique(problem_indices[:, :2], axis=0, return_counts=True)
+        for ((rbin1, rbin2), mu_count) in zip(rbins, mu_counts):
+            if rbin1 > rbin2: continue # this case can be skipped by symmetry
+            print_function(f"WARNING: counts are not positive for radial bin pair {rbin1}, {rbin2} for {mu_count} mu values of {n_mu_check} checked")
 
     vol_r = 4 * np.pi / 3 * (r_bins[:, 1] **3 - r_bins[:, 0] ** 3)
 
@@ -155,7 +168,20 @@ def compute_3pcf_correction_function_from_encore(randoms_pos: np.ndarray[float],
     # fill the edge/corner diagonal elements which are more tricky
     leg_triple[0, 0] = (2 * (2 * leg_triple[1, 0] - leg_triple[2, 0]) + (2 * leg_triple[1, 1] - leg_triple[2, 2])) / 3
     leg_triple[-1, -1] = (2 * (2 * leg_triple[-2, -1] - leg_triple[-3, -1]) + (2 * leg_triple[-2, -2] - leg_triple[-3, -3])) / 3
-    # hopefully this leads to no negative bin counts, worth checking later
+    
+    # check for negative counts, which should be problematic
+    n_mu_check = 2001
+    triple_counts_check = np.zeros([n, n, n_mu_check])
+    mu_values_check = np.linspace(-1, 1, n_mu_check)
+    for ell in range(n_multipoles):
+        triple_counts_check += leg_triple[:, :, ell][:, :, None] * legendre(ell)(mu_values_check)[None, None, :]
+    problem_indices = np.argwhere(triple_counts_check <= 0)
+    if len(problem_indices) > 0:
+        rbins, mu_counts = np.unique(problem_indices[:, :2], axis=0, return_counts=True)
+        for ((rbin1, rbin2), mu_count) in zip(rbins, mu_counts):
+            if rbin1 > rbin2: continue # this case can be skipped by symmetry
+            print_function(("INFO" if rbin1 == rbin2 else "WARNING") + f": counts are not positive for radial bin pair {rbin1}, {rbin2} for {mu_count} mu values of {n_mu_check} checked")
+            # the problem for same-bin pairs is less critical (and seems more likely), for different-bin pairs it is more critical
 
     vol_r = 4 * np.pi / 3 * (r_bins[:, 1] ** 3 - r_bins[:, 0] ** 3) # volume of radial/separation bins as 1D array