Add distance backend abstraction scaffold for issue #172 https://github.com/pysal/pointpats/issues/172

pointpats/distance.py

@@ -0,0 +1,48 @@
+"""
+Distance backends for point pattern analysis.
+
+This module defines an abstraction layer for distance computation
+used throughout pointpats. By default, Euclidean distances are used,
+but alternative backends (e.g., network-based distances) can be
+implemented by subclassing DistanceBackend.
+"""
+
+import abc
+from scipy.spatial.distance import cdist
+
+
+class DistanceBackend(abc.ABC):
+    """Abstract base class for distance computation."""
+
+    @abc.abstractmethod
+    def pairwise(self, points):
+        """
+        Compute pairwise distances between points.
+
+        Parameters
+        ----------
+        points : array-like, shape (n, d)
+
+        Returns
+        -------
+        distances : ndarray, shape (n, n)
+        """
+        raise NotImplementedError
+
+    def within(self, points, r):
+        """
+        Boolean matrix of distances <= r.
+
+        Parameters
+        ----------
+        points : array-like
+        r : float
+        """
+        return self.pairwise(points) <= r
+
+
+class EuclideanDistanceBackend(DistanceBackend):
+    """Default Euclidean distance backend."""
+
+    def pairwise(self, points):
+        return cdist(points, points)

pointpats/distance_statistics.py

@@ -449,7 +449,19 @@ def k(
                 f" is {distances.shape}, but required shape is ({upper_tri_n},) or ({n},{n})"
             )
     else:
-        upper_tri_distances = spatial.distance.pdist(coordinates, metric=metric)
+    # Use distance backend if available, otherwise fall back to scipy
+        try:
+            D = getattr(coordinates, "distance_backend", None)
+            if D is not None:
+                full_distances = D.pairwise(coordinates)
+                upper_tri_distances = full_distances[
+                    numpy.triu_indices_from(full_distances, k=1)
+                ]
+            else:
+                upper_tri_distances = spatial.distance.pdist(coordinates, metric=metric)
+        except AttributeError:
+            upper_tri_distances = spatial.distance.pdist(coordinates, metric=metric)
+
     n_pairs_less_than_d = (upper_tri_distances < support.reshape(-1, 1)).sum(axis=1)
     intensity = n / _area(hull)
     k_estimate = ((n_pairs_less_than_d * 2) / n) / intensity