Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions libpysal/graph/_kernel.py
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ def _parabolic(distances, bandwidth):

def _gaussian(distances, bandwidth):
u = distances / bandwidth
return numpy.exp(-((u / 2) ** 2)) / (numpy.sqrt(2 * numpy.pi))
return numpy.exp(-(u**2)/2) / (numpy.sqrt(2 * numpy.pi))


def _bisquare(distances, bandwidth):
Expand Down Expand Up @@ -140,9 +140,9 @@ def _kernel(
coordinates, ids=ids, valid_geometry_types=_VALID_GEOMETRY_TYPES
)
else:
assert coordinates.shape[0] == coordinates.shape[1], (
"coordinates should represent a distance matrix if metric='precomputed'"
)
assert (
coordinates.shape[0] == coordinates.shape[1]
), "coordinates should represent a distance matrix if metric='precomputed'"
if ids is None:
ids = numpy.arange(coordinates.shape[0])

Expand Down
40 changes: 20 additions & 20 deletions libpysal/graph/tests/test_builders.py
Original file line number Diff line number Diff line change
Expand Up @@ -237,26 +237,26 @@ def test_kernel_precompute(self):
g = graph.Graph.build_kernel(distmat, metric="precomputed")
expected = np.array(
[
0.126,
0.266,
0.174,
0.071,
0.126,
0.329,
0.311,
0.291,
0.266,
0.329,
0.31,
0.205,
0.174,
0.311,
0.31,
0.339,
0.071,
0.291,
0.205,
0.339,
0.04,
0.177,
0.076,
0.013,
0.04,
0.271,
0.242,
0.212,
0.177,
0.271,
0.241,
0.105,
0.076,
0.242,
0.241,
0.288,
0.013,
0.212,
0.105,
0.288,
]
)

Expand Down
4 changes: 2 additions & 2 deletions libpysal/graph/tests/test_kernel.py
Original file line number Diff line number Diff line change
Expand Up @@ -210,8 +210,8 @@ def test_kernels(kernel, grocs):
assert weight.mean() == pytest.approx(0.10312196315841769)
assert weight.max() == pytest.approx(0.749881829575671)
elif kernel == "gaussian":
assert weight.mean() == pytest.approx(0.19932294761630429)
assert weight.max() == pytest.approx(0.3989265663183409)
assert weight.mean() == pytest.approx(0.13787969156713978)
assert weight.max() == pytest.approx(0.39891085285421685)
elif kernel == "bisquare":
assert weight.mean() == pytest.approx(0.09084085210598618)
assert weight.max() == pytest.approx(0.9372045972129259)
Expand Down
33 changes: 17 additions & 16 deletions libpysal/graph/tests/test_triangulation.py
Original file line number Diff line number Diff line change
Expand Up @@ -257,24 +257,25 @@ def test_ids(ids, stores_unique):

def test_kernel():
_, _, weight = _delaunay(cau_coords, kernel="gaussian")

expected = np.array(
[
0.231415,
0.307681,
0.395484,
0.231415,
0.237447,
0.057774,
0.307681,
0.237447,
0.123151,
0.319563,
0.057774,
0.123151,
0.035525,
0.395484,
0.319563,
0.035525,
0.134237,
0.237297,
0.392056,
0.134237,
0.141327,
0.008367,
0.237297,
0.141327,
0.038016,
0.255978,
0.008367,
0.038016,
0.003163,
0.392056,
0.255978,
0.003163,
]
)

Expand Down
265 changes: 265 additions & 0 deletions libpysal/kernels.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,265 @@
"""
kernels.py

This module defines a collection of common kernel functions used for
distance-based weighting in spatial analysis, nonparametric regression,
and density estimation.

Each kernel function takes as input an array of distances and a bandwidth
parameter and returns an array of weights according to the shape of the kernel.

A general ``kernel()`` dispatcher is provided to apply a named kernel or a
user-supplied callable.

Available kernels:
- ``triangular``
- ``parabolic`` (Epanechnikov)
- ``gaussian``
- ``bisquare`` (quartic)
- ``cosine``
- ``exponential``
- ``boxcar`` (uniform)
- ``identity`` (raw distances)

Mathematical Formulation
------------------------

All kernels are evaluated as:

.. math::

K(z), \\quad \\text{where} \\ z = \\frac{d}{h}

- :math:`d` is the distance between points.
- :math:`h` is the kernel bandwidth.
- For :math:`z > 1`, all kernels return :math:`K(z) = 0`.

"""

import numpy


def _trim(d, bandwidth):
"""
Normalize and clip distances to the range [0, 1].

Parameters
----------
d : ndarray
Array of distances.
bandwidth : float
Bandwidth parameter.

Returns
-------
ndarray
Clipped and normalized distances.
"""
return numpy.clip(numpy.abs(d) / bandwidth, 0, 1)


def _triangular(distances, bandwidth):
"""
Triangular kernel function.

Parameters
----------
distances : ndarray
Array of distances.
bandwidth : float
Kernel bandwidth.

Returns
-------
ndarray
Triangular kernel weights.
"""
return 1 - _trim(distances, bandwidth)


def _parabolic(distances, bandwidth):
"""
Parabolic (Epanechnikov) kernel function.

Parameters
----------
distances : ndarray
Array of distances.
bandwidth : float
Kernel bandwidth.

Returns
-------
ndarray
Parabolic kernel weights.
"""
z = _trim(distances, bandwidth)
return 0.75 * (1 - z**2)


def _gaussian(distances, bandwidth):
"""
Gaussian kernel function (truncated at z=1).

Parameters
----------
distances : ndarray
Array of distances.
bandwidth : float
Kernel bandwidth.

Returns
-------
ndarray
Gaussian kernel weights.
"""
z = distances / bandwidth
exponent_term = -0.5 * (z**2)
c = 1 / (bandwidth * numpy.sqrt(2 * numpy.pi))
k = c * numpy.exp(exponent_term)
return numpy.where(z <= 1, k, 0)


def _bisquare(distances, bandwidth):
"""
Bisquare (quartic) kernel function.

Parameters
----------
distances : ndarray
Array of distances.
bandwidth : float
Kernel bandwidth.

Returns
-------
ndarray
Bisquare kernel weights.
"""
z = numpy.clip(distances / bandwidth, 0, 1)
return (15 / 16) * (1 - z**2) ** 2


def _cosine(distances, bandwidth):
"""
Cosine kernel function.

Parameters
----------
distances : ndarray
Array of distances.
bandwidth : float
Kernel bandwidth.

Returns
-------
ndarray
Cosine kernel weights.
"""
z = numpy.clip(distances / bandwidth, 0, 1)
return (numpy.pi / 4) * numpy.cos(numpy.pi / 2 * z)


def _exponential(distances, bandwidth):
"""
Exponential kernel function, truncated at z=1.

Parameters
----------
distances : ndarray
Array of distances.
bandwidth : float
Kernel bandwidth.

Returns
-------
ndarray
Exponential kernel weights.
"""
z = distances / bandwidth
k = numpy.exp(-z)
return numpy.where(z <= 1, k, 0)


def _boxcar(distances, bandwidth):
"""
Boxcar (uniform) kernel function.

Parameters
----------
distances : ndarray
Array of distances.
bandwidth : float
Kernel bandwidth.

Returns
-------
ndarray
Binary weights: 1 if distance < bandwidth, else 0.
"""
return (distances < bandwidth).astype(int)


def _identity(distances, _):
"""
Identity kernel (no weighting, returns raw distances).

Parameters
----------
distances : ndarray
Array of distances.
_ : float
Unused bandwidth parameter.

Returns
-------
ndarray
The raw input 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(distances, bandwidth, kernel="gaussian"):
"""
Evaluate a kernel function over a distance array.

Parameters
----------
distances : ndarray
Array of distances.
bandwidth : float
Kernel bandwidth.
kernel : str or callable, optional
The kernel function to use. If a string, must be one of the predefined
kernel names: 'triangular', 'parabolic', 'gaussian', 'bisquare',
'cosine', 'boxcar', 'discrete', 'exponential', 'identity'.
If callable, it should have the signature `(distances, bandwidth)`.
If None, the 'identity' kernel is used.

Returns
-------
ndarray
Kernel weights.
"""
if callable(kernel):
return kernel(distances, bandwidth)
elif kernel is None:
func = _kernel_functions[None]
else:
func = _kernel_functions[kernel.lower()]

return func(distances, bandwidth)
Loading