Fix: sepal numba compilation options changes the results significantly

src/squidpy/gr/_sepal.py

@@ -202,30 +202,28 @@ def _score_helper(
 @njit(fastmath=True)
 def _diffusion(
     conc: NDArrayA,
-    laplacian: Callable[[NDArrayA, NDArrayA, NDArrayA], float],
+    laplacian: Callable[[NDArrayA, NDArrayA], float],
     n_iter: int,
     sat: NDArrayA,
     sat_idx: NDArrayA,
     unsat: NDArrayA,
     unsat_idx: NDArrayA,
     dt: float = 0.001,
-    D: float = 1.0,
     thresh: float = 1e-8,
 ) -> float:
     """Simulate diffusion process on a regular graph."""
     sat_shape, conc_shape = sat.shape[0], conc.shape[0]
     entropy_arr = np.zeros(n_iter)
     prev_ent = 1.0
     nhood = np.zeros(sat_shape)
-    weights = np.ones(sat_shape)
 
     for i in range(n_iter):
         for j in range(sat_shape):
             nhood[j] = np.sum(conc[sat_idx[j]])
-        d2 = laplacian(conc[sat], nhood, weights)
+        d2 = laplacian(conc[sat], nhood)
 
         dcdt = np.zeros(conc_shape)
-        dcdt[sat] = D * d2
+        dcdt[sat] = d2
         conc[sat] += dcdt[sat] * dt
         conc[unsat] += dcdt[unsat_idx] * dt
         # set values below zero to 0
@@ -246,16 +244,13 @@ def _diffusion(
 def _laplacian_rect(
     centers: NDArrayA,
     nbrs: NDArrayA,
-    h: float,
 ) -> NDArrayA:
     """
     Five point stencil approximation on rectilinear grid.
 
     See `Wikipedia <https://en.wikipedia.org/wiki/Five-point_stencil>`_ for more information.
     """
     d2f: NDArrayA = nbrs - 4 * centers
-    d2f = d2f / h**2
-
     return d2f
 
 
@@ -264,7 +259,6 @@ def _laplacian_rect(
 def _laplacian_hex(
     centers: NDArrayA,
     nbrs: NDArrayA,
-    h: float,
 ) -> NDArrayA:
     """
     Seven point stencil approximation on hexagonal grid.
@@ -275,10 +269,7 @@ def _laplacian_hex(
     Curtis D. Benster, L.V. Kantorovich, V.I. Krylov,
     ISBN-13: 978-0486821603.
     """
-    d2f: NDArrayA = nbrs - 6 * centers
-    d2f = d2f / h**2
-    d2f = (d2f * 2) / 3
-
+    d2f: NDArrayA = (2.0 * nbrs - 12.0 * centers) / 3.0
     return d2f
 
 
@@ -287,11 +278,18 @@ def _laplacian_hex(
 def _entropy(
     xx: NDArrayA,
 ) -> float:
-    """Get entropy of an array."""
+    """Compute Shannon entropy of an array of probability values (in nats)."""
     xnz = xx[xx > 0]
     xs: np.float64 = np.sum(xnz)
+    eps = np.finfo(np.float64).eps  # ~2.22e-16
+    if xs < eps:
+        # 0 because
+        # xn represents probabilities
+        # and p(x)=0 is taken as 0 entropy
+        # see https://stats.stackexchange.com/a/433096
+        return 0.0
     xn = xnz / xs
-    xl = np.log(xn)
+    xl = np.log(np.maximum(xn, eps))
     return float((-xl * xn).sum())