ENH: Add decay and taper arguments to normalize kernel in distance-based weights (#791)

* Allow for normalize=True for gaussian kernel as pdf, False K(0)=1

* bool not binary for type

* Add tests for normalize kw in Kernel

* numpy testing

* Correct reference for kernel defs

* prototyping new _kernels.py

* triangular (#29)

Co-authored-by: knaaptime <knaaptime@Mac.attlocal.net>

* move primative kernels out to libpysal._kernels.py

* renaming kernels.py

* populate self-weight with K(0)

* modify imports

* roll back graph kernels

* Add tests and illustrative nb

* Correct Gaussian kernel function

* Update tests for gaussian correction

* ruff format rules and remove test nbs

* ruff tests

* Relocate kernels from graph to libpysal for reuse and centralization

* move taper/decay logic out of graph module (#30)

* replay levi's change to not remove self-weight (#32)

Co-authored-by: knaaptime <knaaptime@magneto.local>

* gaussian constant term

* update tests for trim

* import proper kernel.py in graph

* Add exclude_loops=True for kernels on graph

* Fix gaussian kernel.

* The coplanar tests assume we don't taper. In main this could be a bug?

* clean up doc strings

* update tests for correct Gaussian kernel

* update tests for gaussian kernel

* Do not decay boxcar kernel

* Update reference for kernels

* cast boxcar kernel to float

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* formatting and import sorting

* skip tests if pyproj not available

* Spelling

* np.in1d is deprecated, using np.isin instead

* Don't assign to unused variable in pytest

* Docs for kernels

* ruff formatting

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* ruff take 2

* more narrative on kernels

* Pass decay and trim kw args into graph kernels

* add decay to graph/_kernel.py

* If k is specified we need the bandwith set to max distance

* Handle taper appropriately

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* ruff

* Propagate decay and taper kwargs

* ruff

* ruff

* use lps kernel

* Propagate decay and taper to triangulation builders

* update build_triangulation

* pass decay and taper in triangulation functions

---------

Co-authored-by: eli knaap <ek@knaaptime.com>
Co-authored-by: knaaptime <knaaptime@Mac.attlocal.net>
Co-authored-by: Levi John Wolf <levi.john.wolf@gmail.com>
Co-authored-by: knaaptime <knaaptime@magneto.local>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: James Gaboardi <jgaboardi@gmail.com>

Author: 7 people

docs/api.rst

@@ -255,6 +255,17 @@ io
    libpysal.io.fileio.FileIO
 
 
+kernels
+-------
+
+
+.. autosummary::
+   :toctree: generated/
+
+   libpysal.kernels.kernel
+
+   
+
 examples
 --------
 

docs/index.rst

@@ -3,9 +3,9 @@
 libpysal: Python Spatial Analysis Library Core
 ==============================================
 
-.. image:: https://github.com/pysal/libpysal/workflows/.github/workflows/unittests.yml/badge.svg
-   :target: https://github.com/pysal/libpysal/actions?query=workflow%3A.github%2Fworkflows%2Funittests.yml
-
+.. image:: https://github.com/pysal/libpysal/actions/workflows/.github/workflows/unittests.yml/badge.svg
+   :target: https://github.com/pysal/libpysal/actions/workflows/unittests.yml
+      
 .. image:: https://img.shields.io/badge/Discord-join%20chat-7289da?style=flat&logo=discord&logoColor=cccccc
    :target: https://discord.gg/BxFTEPFFZn
 

docs/migration.rst

@@ -18,7 +18,7 @@ Members common to W and Graph
 
 
 +-----------------------------------------------------------------------------------------+------------------+
-| Member                                                                                  |   Typee          |
+| Member                                                                                  |   Type           |
 +=========================================================================================+==================+
 | `asymmetry <generated/libpysal.graph.Graph.html#libpysal.graph.Graph.asymmetry>`_       |  builtins.method |
 +-----------------------------------------------------------------------------------------+------------------+

docs/user-guide/intro.rst

@@ -7,9 +7,10 @@ This user guide covers essential features of libpysal, mostly in the form of int
 - work with built-in example spatial data sets
 - how to use spatial weights
 - how to use spatial graphs
+- how to use kernels
 
 Notebooks cover just a small selection of functions as an illustration of
-principles. For a full overview of momepy capabilities, head to the `API <../api.rst>`_.
+principles. For a full overview of libpysal's capabilities, head to the `API <../api.rst>`_.
 
 
 .. toctree::
@@ -18,3 +19,4 @@ principles. For a full overview of momepy capabilities, head to the `API <../api
    Datasets <data/examples>
    Weights <weights/intro>
    Graphs  <graph/intro>
+   Kernels <kernels>

docs/user-guide/kernels.ipynb

@@ -2,6 +2,9 @@
 import pandas
 from scipy import optimize, sparse, spatial, stats
 
+from libpysal.kernels import _kernel_functions
+from libpysal.kernels import kernel as _lps_kernel
+
 from ._utils import (
     CoplanarError,
     _build_coplanarity_lookup,
@@ -22,59 +25,6 @@
 _VALID_GEOMETRY_TYPES = ["Point"]
 
 
-def _triangular(distances, bandwidth):
-    u = numpy.clip(distances / bandwidth, 0, 1)
-    return 1 - u
-
-
-def _parabolic(distances, bandwidth):
-    u = numpy.clip(distances / bandwidth, 0, 1)
-    return 0.75 * (1 - u**2)
-
-
-def _gaussian(distances, bandwidth):
-    u = distances / bandwidth
-    return numpy.exp(-((u / 2) ** 2)) / (numpy.sqrt(2 * numpy.pi))
-
-
-def _bisquare(distances, bandwidth):
-    u = numpy.clip(distances / bandwidth, 0, 1)
-    return (15 / 16) * (1 - u**2) ** 2
-
-
-def _cosine(distances, bandwidth):
-    u = numpy.clip(distances / bandwidth, 0, 1)
-    return (numpy.pi / 4) * numpy.cos(numpy.pi / 2 * u)
-
-
-def _exponential(distances, bandwidth):
-    u = distances / bandwidth
-    return numpy.exp(-u)
-
-
-def _boxcar(distances, bandwidth):
-    r = (distances < bandwidth).astype(int)
-    return r
-
-
-def _identity(distances, _):
-    return distances
-
-
-_kernel_functions = {
-    "triangular": _triangular,
-    "parabolic": _parabolic,
-    "gaussian": _gaussian,
-    "bisquare": _bisquare,
-    "cosine": _cosine,
-    "boxcar": _boxcar,
-    "discrete": _boxcar,
-    "exponential": _exponential,
-    "identity": _identity,
-    None: _identity,
-}
-
-
 def _kernel(
     coordinates,
     bandwidth=None,
@@ -84,8 +34,10 @@ def _kernel(
     ids=None,
     p=2,
     taper=True,
+    decay=False,
     coplanar="raise",
     resolve_isolates=True,
+    exclude_self_weights=True,
 ):
     """
     Compute a kernel function over a distance matrix.
@@ -132,8 +84,18 @@ def _kernel(
         parameter for minkowski metric, ignored if metric != "minkowski".
     taper : bool (default: True)
         remove links with a weight equal to zero
+    decay : bool (default: False)
+        whether to calculate the kernel using the decay formulation.
+        In the decay form, a kernel measures the distance decay in
+        similarity between observations. It varies from from maximal
+        similarity (1) at a distance of zero to minimal similarity (0
+        or negative) at some very large (possibly infinite) distance.
+        Otherwise, kernel functions are treated as proper
+        volume-preserving probability distributions.
     resolve_isolates : bool
         Try to resolve isolates. Can be disabled if we are dealing with cliques later.
+    exclude_self_weights : bool (default: True)
+        Remove self-weights
     """
     if metric != "precomputed":
         coordinates, ids, _ = _validate_geometry_input(
@@ -188,16 +150,17 @@ def _kernel(
             data = sq.flatten()
             i = numpy.tile(numpy.arange(sq.shape[0]), sq.shape[0])
             j = numpy.repeat(numpy.arange(sq.shape[0]), sq.shape[0])
-            # remove diagonal
-            data = numpy.delete(data, numpy.arange(0, data.size, sq.shape[0] + 1))
-            i = numpy.delete(i, numpy.arange(0, i.size, sq.shape[0] + 1))
-            j = numpy.delete(j, numpy.arange(0, j.size, sq.shape[0] + 1))
-            # construct sparse
+
+            if exclude_self_weights:
+                data = numpy.delete(data, numpy.arange(0, data.size, sq.shape[0] + 1))
+                i = numpy.delete(i, numpy.arange(0, i.size, sq.shape[0] + 1))
+                j = numpy.delete(j, numpy.arange(0, j.size, sq.shape[0] + 1))
+
             d = sparse.csc_array((data, (i, j)))
         else:
             d = sparse.csc_array(coordinates)
     if bandwidth is None:
-        bandwidth = numpy.percentile(d.data, 25)
+        bandwidth = numpy.percentile(d.data, 25) if k is None else d.data.max()
     elif bandwidth == "auto":
         if (kernel == "identity") or (kernel is None):
             bandwidth = numpy.nan  # ignored by identity
@@ -206,11 +169,9 @@ def _kernel(
     if callable(kernel):
         d.data = kernel(d.data, bandwidth)
     else:
-        d.data = _kernel_functions[kernel](d.data, bandwidth)
-
+        d.data = _lps_kernel(d.data, bandwidth, kernel=kernel, taper=taper, decay=decay)
     if taper:
         d.eliminate_zeros()
-
     return _sparse_to_arrays(d, ids=ids, resolve_isolates=resolve_isolates)
 
 

libpysal/graph/_kernel.py

@@ -1,6 +1,6 @@
 import numpy as np
 
-from ._kernel import _kernel_functions
+from ._kernel import _lps_kernel
 from ._utils import _induce_cliques, _validate_geometry_input
 
 
@@ -62,7 +62,9 @@ def pdna_to_adj(origins, network, node_ids, threshold):
     return adj
 
 
-def build_travel_graph(df, network, threshold, mapping_distance, kernel=None):
+def build_travel_graph(
+    df, network, threshold, mapping_distance, kernel=None, taper=True, decay=False
+):
     """Compute the shortest path between gdf centroids via a pandana.Network
     and return an adjacency list with weight=cost. Note unlike distance_band,
     :math:`G_{ij}` and :math:`G_{ji}` are often different because travel networks
@@ -87,6 +89,16 @@ def build_travel_graph(df, network, threshold, mapping_distance, kernel=None):
         libpysal.graph.Graph.build_kernel for more information on kernel
         transformation options. Default is None, in which case the Graph weight
         is pure distance between focal and neighbor
+    taper : bool (default: True)
+        remove links with a weight equal to zero
+    decay : bool (default: False)
+        whether to calculate the kernel using the decay formulation.
+        In the decay form, a kernel measures the distance decay in
+        similarity between observations. It varies from from maximal
+        similarity (1) at a distance of zero to minimal similarity (0
+        or negative) at some very large (possibly infinite) distance.
+        Otherwise, kernel functions are treated as proper
+        volume-preserving probability distributions.
 
     Returns
     -------
@@ -125,10 +137,16 @@ def build_travel_graph(df, network, threshold, mapping_distance, kernel=None):
     )
 
     if callable(kernel):
-        adj_cliques["weight"] = kernel(adj_cliques["weight"], threshold)
+        adj_cliques["weight"] = kernel(
+            adj_cliques["weight"], threshold, decay=decay, taper=taper
+        )
     else:
-        adj_cliques["weight"] = _kernel_functions[kernel](
-            adj_cliques["weight"], threshold
+        adj_cliques["weight"] = _lps_kernel(
+            adj_cliques["weight"].values,
+            kernel=kernel,
+            bandwidth=threshold,
+            decay=decay,
+            taper=taper,
         )
 
     return adj_cliques

libpysal/graph/_network.py

@@ -50,6 +50,8 @@ def tri_with_validation(
         kernel=None,
         bandwidth=None,
         seed=None,
+        decay=False,
+        taper=False,
         **kwargs,
     ):
         if coplanar not in ["raise", "jitter", "clique"]:
@@ -85,8 +87,9 @@ def tri_with_validation(
                 metric="precomputed",
                 kernel=kernel,
                 bandwidth=bandwidth,
-                taper=False,
+                taper=taper,
                 resolve_isolates=False,  # no isolates in triangulation
+                decay=decay,
             )
         # create adjacency
         adjtable = pandas.DataFrame.from_dict(