#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Temporal segmentation
=====================

Recurrence and self-similarity
------------------------------
.. autosummary::
    :toctree: generated/

    cross_similarity
    recurrence_matrix
    recurrence_to_lag
    lag_to_recurrence
    timelag_filter
    path_enhance

Temporal clustering
-------------------
.. autosummary::
    :toctree: generated/

    agglomerative
    subsegment
"""

from decorator import decorator

import numpy as np
import scipy
import scipy.signal
import scipy.ndimage

import sklearn
import sklearn.cluster
import sklearn.feature_extraction
import sklearn.neighbors

from ._cache import cache
from . import util
from .filters import diagonal_filter
from .util.exceptions import ParameterError
from typing import Any, Callable, Optional, TypeVar, Union, overload
from typing_extensions import Literal
from ._typing import _WindowSpec, _FloatLike_co

__all__ = [
    "cross_similarity",
    "recurrence_matrix",
    "recurrence_to_lag",
    "lag_to_recurrence",
    "timelag_filter",
    "agglomerative",
    "subsegment",
    "path_enhance",
]


@overload
def cross_similarity(
    data: np.ndarray,
    data_ref: np.ndarray,
    *,
    k: Optional[int] = ...,
    metric: str = ...,
    sparse: Literal[False] = ...,
    mode: str = ...,
    bandwidth: Optional[Union[np.ndarray, _FloatLike_co, str]] = None,
    full: bool = False,
) -> np.ndarray:
    ...


@overload
def cross_similarity(
    data: np.ndarray,
    data_ref: np.ndarray,
    *,
    k: Optional[int] = ...,
    metric: str = ...,
    sparse: Literal[True] = ...,
    mode: str = ...,
    bandwidth: Optional[Union[np.ndarray, _FloatLike_co, str]] = None,
    full: bool = False,
) -> scipy.sparse.csc_matrix:
    ...


@cache(level=30)
def cross_similarity(
    data: np.ndarray,
    data_ref: np.ndarray,
    *,
    k: Optional[int] = None,
    metric: str = "euclidean",
    sparse: bool = False,
    mode: str = "connectivity",
    bandwidth: Optional[Union[np.ndarray, _FloatLike_co, str]] = None,
    full: bool = False,
) -> Union[np.ndarray, scipy.sparse.csc_matrix]:
    """Compute cross-similarity from one data sequence to a reference sequence.

    The output is a matrix ``xsim``, where ``xsim[i, j]`` is non-zero
    if ``data_ref[..., i]`` is a k-nearest neighbor of ``data[..., j]``.

    Parameters
    ----------
    data : np.ndarray [shape=(..., d, n)]
        A feature matrix for the comparison sequence.
        If the data has more than two dimensions (e.g., for multi-channel inputs),
        the leading dimensions are flattened prior to comparison.
        For example, a stereo input with shape `(2, d, n)` is
        automatically reshaped to `(2 * d, n)`.

    data_ref : np.ndarray [shape=(..., d, n_ref)]
        A feature matrix for the reference sequence
        If the data has more than two dimensions (e.g., for multi-channel inputs),
        the leading dimensions are flattened prior to comparison.
        For example, a stereo input with shape `(2, d, n_ref)` is
        automatically reshaped to `(2 * d, n_ref)`.

    k : int > 0 [scalar] or None
        the number of nearest-neighbors for each sample

        Default: ``k = 2 * ceil(sqrt(n_ref))``,
        or ``k = 2`` if ``n_ref <= 3``

    metric : str
        Distance metric to use for nearest-neighbor calculation.

        See `sklearn.neighbors.NearestNeighbors` for details.

    sparse : bool [scalar]
        if False, returns a dense type (ndarray)
        if True, returns a sparse type (scipy.sparse.csc_matrix)

    mode : str, {'connectivity', 'distance', 'affinity'}
        If 'connectivity', a binary connectivity matrix is produced.

        If 'distance', then a non-zero entry contains the distance between
        points.

        If 'affinity', then non-zero entries are mapped to
        ``exp( - distance(i, j) / bandwidth)`` where ``bandwidth`` is
        as specified below.

    bandwidth : None, float > 0, ndarray, or str
        str options include ``{'med_k_scalar', 'mean_k', 'gmean_k', 'mean_k_avg', 'gmean_k_avg', 'mean_k_avg_and_pair'}``

        If ndarray is supplied, use ndarray as bandwidth for each i,j pair.

        If using ``mode='affinity'``, this can be used to set the
        bandwidth on the affinity kernel.

        If no value is provided or ``None``, default to ``'med_k_scalar'``.

        If ``bandwidth='med_k_scalar'``, bandwidth is set automatically to the median
        distance to the k'th nearest neighbor of each ``data[:, i]``.

        If ``bandwidth='mean_k'``, bandwidth is estimated for each sample-pair (i, j) by taking the
        arithmetic mean between distances to the k-th nearest neighbor for sample i and sample j.

        If ``bandwidth='gmean_k'``, bandwidth is estimated for each sample-pair (i, j) by taking the
        geometric mean between distances to the k-th nearest neighbor for sample i and j [#z]_.

        If ``bandwidth='mean_k_avg'``, bandwidth is estimated for each sample-pair (i, j) by taking the
        arithmetic mean between the average distances to the first k-th nearest neighbors for
        sample i and sample j.
        This is similar to the approach in Wang et al. (2014) [#w]_ but does not include the distance
        between i and j.

        If ``bandwidth='gmean_k_avg'``, bandwidth is estimated for each sample-pair (i, j) by taking the
        geometric mean between the average distances to the first k-th nearest neighbors for
        sample i and sample j.

        If ``bandwidth='mean_k_avg_and_pair'``, bandwidth is estimated for each sample-pair (i, j) by
        taking the arithmetic mean between three terms: the average distances to the first
        k-th nearest neighbors for sample i and sample j respectively, as well as
        the distance between i and j.
        This is similar to the approach in Wang et al. (2014). [#w]_

        .. [#z] Zelnik-Manor, Lihi, and Pietro Perona. (2004).
            "Self-tuning spectral clustering." Advances in neural information processing systems 17.

        .. [#w] Wang, Bo, et al. (2014).
            "Similarity network fusion for aggregating data types on a genomic scale." Nat Methods 11, 333–337.
            https://doi.org/10.1038/nmeth.2810

    full : bool
        If using ``mode ='affinity'`` or ``mode='distance'``, this option can be used to compute
        the full affinity or distance matrix as opposed a sparse matrix with only none-zero terms
        for the first k-neighbors of each sample.
        This option has no effect when using ``mode='connectivity'``.

        When using ``mode='distance'``, setting ``full=True`` will ignore ``k`` and ``width``.
        When using ``mode='affinity'``, setting ``full=True`` will use ``k`` exclusively for
        bandwidth estimation, and ignore ``width``.

    Returns
    -------
    xsim : np.ndarray or scipy.sparse.csc_matrix, [shape=(n_ref, n)]
        Cross-similarity matrix

    See Also
    --------
    recurrence_matrix
    recurrence_to_lag
    librosa.feature.stack_memory
    sklearn.neighbors.NearestNeighbors
    scipy.spatial.distance.cdist

    Notes
    -----
    This function caches at level 30.

    Examples
    --------
    Find nearest neighbors in CQT space between two sequences

    >>> hop_length = 1024
    >>> y_ref, sr = librosa.load(librosa.ex('pistachio'))
    >>> y_comp, sr = librosa.load(librosa.ex('pistachio'), offset=10)
    >>> chroma_ref = librosa.feature.chroma_cqt(y=y_ref, sr=sr, hop_length=hop_length)
    >>> chroma_comp = librosa.feature.chroma_cqt(y=y_comp, sr=sr, hop_length=hop_length)
    >>> # Use time-delay embedding to get a cleaner recurrence matrix
    >>> x_ref = librosa.feature.stack_memory(chroma_ref, n_steps=10, delay=3)
    >>> x_comp = librosa.feature.stack_memory(chroma_comp, n_steps=10, delay=3)
    >>> xsim = librosa.segment.cross_similarity(x_comp, x_ref)

    Or fix the number of nearest neighbors to 5

    >>> xsim = librosa.segment.cross_similarity(x_comp, x_ref, k=5)

    Use cosine similarity instead of Euclidean distance

    >>> xsim = librosa.segment.cross_similarity(x_comp, x_ref, metric='cosine')

    Use an affinity matrix instead of binary connectivity

    >>> xsim_aff = librosa.segment.cross_similarity(x_comp, x_ref, metric='cosine', mode='affinity')

    Plot the feature and recurrence matrices

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(ncols=2, sharex=True, sharey=True)
    >>> imgsim = librosa.display.specshow(xsim, x_axis='s', y_axis='s',
    ...                          hop_length=hop_length, ax=ax[0])
    >>> ax[0].set(title='Binary cross-similarity (symmetric)')
    >>> imgaff = librosa.display.specshow(xsim_aff, x_axis='s', y_axis='s',
    ...                          cmap='magma_r', hop_length=hop_length, ax=ax[1])
    >>> ax[1].set(title='Cross-affinity')
    >>> ax[1].label_outer()
    >>> fig.colorbar(imgsim, ax=ax[0], orientation='horizontal', ticks=[0, 1])
    >>> fig.colorbar(imgaff, ax=ax[1], orientation='horizontal')
    """
    data_ref = np.atleast_2d(data_ref)
    data = np.atleast_2d(data)

    if not np.allclose(data_ref.shape[:-1], data.shape[:-1]):
        raise ParameterError(
            f"data_ref.shape={data_ref.shape} and data.shape={data.shape} do not match on leading dimension(s)"
        )

    # swap data axes so the feature axis is last
    data_ref = np.swapaxes(data_ref, -1, 0)
    n_ref = data_ref.shape[0]
    # Use F-ordering for reshape to preserve leading axis
    data_ref = data_ref.reshape((n_ref, -1), order="F")

    data = np.swapaxes(data, -1, 0)
    n = data.shape[0]
    data = data.reshape((n, -1), order="F")

    if mode not in ["connectivity", "distance", "affinity"]:
        raise ParameterError(
            (
                f"Invalid mode='{mode}'. Must be one of "
                "['connectivity', 'distance', 'affinity']"
            )
        )
    if k is None:
        k = min(n_ref, 2 * np.ceil(np.sqrt(n_ref)))

    k = int(k)

    # using k for bandwidth estimation also and decouple k for full mode
    bandwidth_k = k
    if full and (mode != "connectivity"):
        k = n

    # Build the neighbor search object
    # `auto` mode does not work with some choices of metric.  Rather than special-case
    # those here, we instead use a fall-back to brute force if auto fails.
    try:
        knn = sklearn.neighbors.NearestNeighbors(
            n_neighbors=min(n_ref, k), metric=metric, algorithm="auto"
        )
    except ValueError:
        knn = sklearn.neighbors.NearestNeighbors(
            n_neighbors=min(n_ref, k), metric=metric, algorithm="brute"
        )

    knn.fit(data_ref)

    # Get the knn graph
    if mode == "affinity":
        # sklearn's nearest neighbor doesn't support affinity,
        # so we use distance here and then do the conversion post-hoc
        kng_mode = "distance"
    else:
        kng_mode = mode

    xsim = knn.kneighbors_graph(X=data, mode=kng_mode).tolil()

    if not full:
        # Retain only the top-k links per point
        for i in range(n):
            # Get the links from point i
            links = xsim[i].nonzero()[1]

            # Order them ascending
            idx = links[np.argsort(xsim[i, links].toarray())][0]

            # Everything past the kth closest gets squashed
            xsim[i, idx[k:]] = 0

    # Convert a compressed sparse row (CSR) format
    xsim = xsim.tocsr()
    xsim.eliminate_zeros()

    if mode == "connectivity":
        xsim = xsim.astype(bool)
    elif mode == "affinity":
        aff_bandwidth = __affinity_bandwidth(xsim, bandwidth, bandwidth_k)
        xsim.data[:] = np.exp(xsim.data / (-1 * aff_bandwidth))

    # Transpose to n_ref by n
    xsim = xsim.T

    if not sparse:
        xsim = xsim.toarray()

    return xsim


@overload
def recurrence_matrix(
    data: np.ndarray,
    *,
    k: Optional[int] = ...,
    width: int = ...,
    metric: str = ...,
    sym: bool = ...,
    sparse: Literal[True] = ...,
    mode: str = ...,
    bandwidth: Optional[Union[np.ndarray, _FloatLike_co, str]] = ...,
    self: bool = ...,
    axis: int = ...,
    full: bool = False,
) -> scipy.sparse.csc_matrix:
    ...


@overload
def recurrence_matrix(
    data: np.ndarray,
    *,
    k: Optional[int] = ...,
    width: int = ...,
    metric: str = ...,
    sym: bool = ...,
    sparse: Literal[False] = ...,
    mode: str = ...,
    bandwidth: Optional[Union[np.ndarray, _FloatLike_co, str]] = ...,
    self: bool = ...,
    axis: int = ...,
    full: bool = False,
) -> np.ndarray:
    ...


@cache(level=30)
def recurrence_matrix(
    data: np.ndarray,
    *,
    k: Optional[int] = None,
    width: int = 1,
    metric: str = "euclidean",
    sym: bool = False,
    sparse: bool = False,
    mode: str = "connectivity",
    bandwidth: Optional[Union[np.ndarray, _FloatLike_co, str]] = None,
    self: bool = False,
    axis: int = -1,
    full: bool = False,
) -> Union[np.ndarray, scipy.sparse.csc_matrix]:
    """Compute a recurrence matrix from a data matrix.

    ``rec[i, j]`` is non-zero if ``data[..., i]`` is a k-nearest neighbor
    of ``data[..., j]`` and ``|i - j| >= width``

    The specific value of ``rec[i, j]`` can have several forms, governed
    by the ``mode`` parameter below:

        - Connectivity: ``rec[i, j] = 1 or 0`` indicates that frames ``i`` and ``j`` are repetitions

        - Affinity: ``rec[i, j] > 0`` measures how similar frames ``i`` and ``j`` are.  This is also
          known as a (sparse) self-similarity matrix.

        - Distance: ``rec[i, j] > 0`` measures how distant frames ``i`` and ``j`` are.  This is also
          known as a (sparse) self-distance matrix.

    The general term *recurrence matrix* can refer to any of the three forms above.

    Parameters
    ----------
    data : np.ndarray [shape=(..., d, n)]
        A feature matrix.
        If the data has more than two dimensions (e.g., for multi-channel inputs),
        the leading dimensions are flattened prior to comparison.
        For example, a stereo input with shape `(2, d, n)` is
        automatically reshaped to `(2 * d, n)`.

    k : int > 0 [scalar] or None
        the number of nearest-neighbors for each sample

        Default: ``k = 2 * ceil(sqrt(t - 2 * width + 1))``,
        or ``k = 2`` if ``t <= 2 * width + 1``

    width : int >= 1 [scalar]
        only link neighbors ``(data[..., i], data[..., j])``
        if ``|i - j| >= width``

        ``width`` cannot exceed the length of the data.

    metric : str
        Distance metric to use for nearest-neighbor calculation.

        See `sklearn.neighbors.NearestNeighbors` for details.

    sym : bool [scalar]
        set ``sym=True`` to only link mutual nearest-neighbors

    sparse : bool [scalar]
        if False, returns a dense type (ndarray)
        if True, returns a sparse type (scipy.sparse.csc_matrix)

    mode : str, {'connectivity', 'distance', 'affinity'}
        If 'connectivity', a binary connectivity matrix is produced.

        If 'distance', then a non-zero entry contains the distance between
        points.

        If 'affinity', then non-zero entries are mapped to
        ``exp( - distance(i, j) / bandwidth)`` where ``bandwidth`` is
        as specified below.

    bandwidth : None, float > 0, ndarray, or str
        str options include ``{'med_k_scalar', 'mean_k', 'gmean_k', 'mean_k_avg', 'gmean_k_avg', 'mean_k_avg_and_pair'}``

        If ndarray is supplied, use ndarray as bandwidth for each i,j pair.

        If using ``mode='affinity'``, the ``bandwidth`` option can be used to set the
        bandwidth on the affinity kernel.

        If no value is provided or ``None``, default to ``'med_k_scalar'``.

        If ``bandwidth='med_k_scalar'``, a scalar bandwidth is set to the median distance
        of the k-th nearest neighbor for all samples.

        If ``bandwidth='mean_k'``, bandwidth is estimated for each sample-pair (i, j) by taking the
        arithmetic mean between distances to the k-th nearest neighbor for sample i and sample j.

        If ``bandwidth='gmean_k'``, bandwidth is estimated for each sample-pair (i, j) by taking the
        geometric mean between distances to the k-th nearest neighbor for sample i and j [#z]_.

        If ``bandwidth='mean_k_avg'``, bandwidth is estimated for each sample-pair (i, j) by taking the
        arithmetic mean between the average distances to the first k-th nearest neighbors for
        sample i and sample j.
        This is similar to the approach in Wang et al. (2014) [#w]_ but does not include the distance
        between i and j.

        If ``bandwidth='gmean_k_avg'``, bandwidth is estimated for each sample-pair (i, j) by taking the
        geometric mean between the average distances to the first k-th nearest neighbors for
        sample i and sample j.

        If ``bandwidth='mean_k_avg_and_pair'``, bandwidth is estimated for each sample-pair (i, j) by
        taking the arithmetic mean between three terms: the average distances to the first
        k-th nearest neighbors for sample i and sample j respectively, as well as
        the distance between i and j.
        This is similar to the approach in Wang et al. (2014). [#w]_

        .. [#z] Zelnik-Manor, Lihi, and Pietro Perona. (2004).
            "Self-tuning spectral clustering." Advances in neural information processing systems 17.

        .. [#w] Wang, Bo, et al. (2014).
            "Similarity network fusion for aggregating data types on a genomic scale." Nat Methods 11, 333–337.
            https://doi.org/10.1038/nmeth.2810

    self : bool
        If ``True``, then the main diagonal is populated with self-links:
        0 if ``mode='distance'``, and 1 otherwise.

        If ``False``, the main diagonal is left empty.

    axis : int
        The axis along which to compute recurrence.
        By default, the last index (-1) is taken.

    full : bool
        If using ``mode ='affinity'`` or ``mode='distance'``, this option can be used to compute
        the full affinity or distance matrix as opposed a sparse matrix with only none-zero terms
        for the first k-neighbors of each sample.
        This option has no effect when using ``mode='connectivity'``.

        When using ``mode='distance'``, setting ``full=True`` will ignore ``k`` and ``width``.
        When using ``mode='affinity'``, setting ``full=True`` will use ``k`` exclusively for
        bandwidth estimation, and ignore ``width``.

    Returns
    -------
    rec : np.ndarray or scipy.sparse.csc_matrix, [shape=(t, t)]
        Recurrence matrix

    See Also
    --------
    sklearn.neighbors.NearestNeighbors
    scipy.spatial.distance.cdist
    librosa.feature.stack_memory
    recurrence_to_lag

    Notes
    -----
    This function caches at level 30.

    Examples
    --------
    Find nearest neighbors in CQT space

    >>> y, sr = librosa.load(librosa.ex('nutcracker'))
    >>> hop_length = 1024
    >>> chroma = librosa.feature.chroma_cqt(y=y, sr=sr, hop_length=hop_length)
    >>> # Use time-delay embedding to get a cleaner recurrence matrix
    >>> chroma_stack = librosa.feature.stack_memory(chroma, n_steps=10, delay=3)
    >>> R = librosa.segment.recurrence_matrix(chroma_stack)

    Or fix the number of nearest neighbors to 5

    >>> R = librosa.segment.recurrence_matrix(chroma_stack, k=5)

    Suppress neighbors within +- 7 frames

    >>> R = librosa.segment.recurrence_matrix(chroma_stack, width=7)

    Use cosine similarity instead of Euclidean distance

    >>> R = librosa.segment.recurrence_matrix(chroma_stack, metric='cosine')

    Require mutual nearest neighbors

    >>> R = librosa.segment.recurrence_matrix(chroma_stack, sym=True)

    Use an affinity matrix instead of binary connectivity

    >>> R_aff = librosa.segment.recurrence_matrix(chroma_stack, metric='cosine',
    ...                                           mode='affinity')

    Plot the feature and recurrence matrices

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(ncols=2, sharex=True, sharey=True)
    >>> imgsim = librosa.display.specshow(R, x_axis='s', y_axis='s',
    ...                          hop_length=hop_length, ax=ax[0])
    >>> ax[0].set(title='Binary recurrence (symmetric)')
    >>> imgaff = librosa.display.specshow(R_aff, x_axis='s', y_axis='s',
    ...                          hop_length=hop_length, cmap='magma_r', ax=ax[1])
    >>> ax[1].set(title='Affinity recurrence')
    >>> ax[1].label_outer()
    >>> fig.colorbar(imgsim, ax=ax[0], orientation='horizontal', ticks=[0, 1])
    >>> fig.colorbar(imgaff, ax=ax[1], orientation='horizontal')
    """
    data = np.atleast_2d(data)

    # Swap observations to the first dimension and flatten the rest
    data = np.swapaxes(data, axis, 0)
    t = data.shape[0]
    # Use F-ordering here to preserve leading axis layout
    data = data.reshape((t, -1), order="F")

    if width < 1 or width >= (t - 1) // 2:
        raise ParameterError(
            "width={} must be at least 1 and at most (data.shape[{}] - 1) // 2={}".format(
                width, axis, (t - 1) // 2
            )
        )

    if mode not in ["connectivity", "distance", "affinity"]:
        raise ParameterError(
            (
                f"Invalid mode='{mode}'. Must be one of "
                "['connectivity', 'distance', 'affinity']"
            )
        )
    if k is None:
        k = 2 * np.ceil(np.sqrt(t - 2 * width + 1))

    k = int(k)

    # using k for bandwidth estimation also and decouple k for full mode
    bandwidth_k = k
    if full and (mode != "connectivity"):
        k = t

    # Build the neighbor search object
    try:
        knn = sklearn.neighbors.NearestNeighbors(
            n_neighbors=min(t - 1, k + 2 * width), metric=metric, algorithm="auto"
        )
    except ValueError:
        knn = sklearn.neighbors.NearestNeighbors(
            n_neighbors=min(t - 1, k + 2 * width), metric=metric, algorithm="brute"
        )

    knn.fit(data)

    # Get the knn graph
    if mode == "affinity":
        kng_mode = "distance"
    else:
        kng_mode = mode

    rec = knn.kneighbors_graph(mode=kng_mode).tolil()

    if not full:
        # Remove connections within width
        for diag in range(-width + 1, width):
            rec.setdiag(0, diag)

        # Retain only the top-k links per point
        for i in range(t):
            # Get the links from point i
            links = rec[i].nonzero()[1]

            # Order them ascending
            idx = links[np.argsort(rec[i, links].toarray())][0]

            # Everything past the kth closest gets squashed
            rec[i, idx[k:]] = 0

    if self:
        if mode == "connectivity":
            rec.setdiag(1)
        elif mode == "affinity":
            # we need to keep the self-loop in here, but not mess up the
            # bandwidth estimation
            #
            # using negative distances here preserves the structure without changing
            # the statistics of the data
            rec.setdiag(-1)
    else:
        rec.setdiag(0)

    # symmetrize
    if sym:
        # Note: this operation produces a CSR (compressed sparse row) matrix!
        # This is why we have to do it after filling the diagonal in self-mode
        rec = rec.minimum(rec.T)

    rec = rec.tocsr()
    rec.eliminate_zeros()

    if mode == "connectivity":
        rec = rec.astype(bool)
    elif mode == "affinity":
        # Set all the negatives back to 0
        # Negatives are temporarily inserted above to preserve the sparsity structure
        # of the matrix without corrupting the bandwidth calculations
        rec.data[rec.data < 0] = 0.0
        aff_bandwidth = __affinity_bandwidth(rec, bandwidth, bandwidth_k)
        rec.data[:] = np.exp(rec.data / (-1 * aff_bandwidth))

    # Transpose to be column-major
    rec = rec.T

    if not sparse:
        rec = rec.toarray()

    return rec


_ArrayOrSparseMatrix = TypeVar(
    "_ArrayOrSparseMatrix", bound=Union[np.ndarray, scipy.sparse.spmatrix]
)


def recurrence_to_lag(
    rec: _ArrayOrSparseMatrix, *, pad: bool = True, axis: int = -1
) -> _ArrayOrSparseMatrix:
    """Convert a recurrence matrix into a lag matrix.

        ``lag[i, j] == rec[i+j, j]``

    This transformation turns diagonal structures in the recurrence matrix
    into horizontal structures in the lag matrix.
    These horizontal structures can be used to infer changes in the repetition
    structure of a piece, e.g., the beginning of a new section as done in [#]_.

    .. [#] Serra, J., Müller, M., Grosche, P., & Arcos, J. L. (2014).
           Unsupervised music structure annotation by time series structure
           features and segment similarity.
           IEEE Transactions on Multimedia, 16(5), 1229-1240.

    Parameters
    ----------
    rec : np.ndarray, or scipy.sparse.spmatrix [shape=(n, n)]
        A (binary) recurrence matrix, as returned by `recurrence_matrix`

    pad : bool
        If False, ``lag`` matrix is square, which is equivalent to
        assuming that the signal repeats itself indefinitely.

        If True, ``lag`` is padded with ``n`` zeros, which eliminates
        the assumption of repetition.

    axis : int
        The axis to keep as the ``time`` axis.
        The alternate axis will be converted to lag coordinates.

    Returns
    -------
    lag : np.ndarray
        The recurrence matrix in (lag, time) (if ``axis=1``)
        or (time, lag) (if ``axis=0``) coordinates

    Raises
    ------
    ParameterError : if ``rec`` is non-square

    See Also
    --------
    recurrence_matrix
    lag_to_recurrence
    util.shear

    Examples
    --------
    >>> y, sr = librosa.load(librosa.ex('nutcracker'))
    >>> hop_length = 1024
    >>> chroma = librosa.feature.chroma_cqt(y=y, sr=sr, hop_length=hop_length)
    >>> chroma_stack = librosa.feature.stack_memory(chroma, n_steps=10, delay=3)
    >>> recurrence = librosa.segment.recurrence_matrix(chroma_stack)
    >>> lag_pad = librosa.segment.recurrence_to_lag(recurrence, pad=True)
    >>> lag_nopad = librosa.segment.recurrence_to_lag(recurrence, pad=False)

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(nrows=2, sharex=True)
    >>> librosa.display.specshow(lag_pad, x_axis='time', y_axis='lag',
    ...                          hop_length=hop_length, ax=ax[0])
    >>> ax[0].set(title='Lag (zero-padded)')
    >>> ax[0].label_outer()
    >>> librosa.display.specshow(lag_nopad, x_axis='time', y_axis='lag',
    ...                          hop_length=hop_length, ax=ax[1])
    >>> ax[1].set(title='Lag (no padding)')
    """
    axis = np.abs(axis)

    if rec.ndim != 2 or rec.shape[0] != rec.shape[1]:
        raise ParameterError(f"non-square recurrence matrix shape: {rec.shape}")

    sparse = scipy.sparse.issparse(rec)

    if sparse:
        # suppress type check here, mypy doesn't know about issparse
        fmt = rec.format  # type: ignore

    t = rec.shape[axis]

    if pad:
        if sparse:
            padding = np.asarray([[1, 0]], dtype=rec.dtype).swapaxes(axis, 0)
            if axis == 0:
                rec_fmt = "csr"
            else:
                rec_fmt = "csc"
            rec = scipy.sparse.kron(padding, rec, format=rec_fmt)
        else:
            padding = np.array([(0, 0), (0, 0)])
            padding[(1 - axis), :] = [0, t]
            # Suppress type check, mypy doesn't know that rec is an ndarray here
            rec = np.pad(rec, padding, mode="constant")  # type: ignore

    lag: _ArrayOrSparseMatrix = util.shear(rec, factor=-1, axis=axis)

    if sparse:
        # Suppress type check, mypy doesn't know
        # that lag is sparse here
        lag = lag.asformat(fmt)  # type: ignore

    return lag


def lag_to_recurrence(
    lag: _ArrayOrSparseMatrix, *, axis: int = -1
) -> _ArrayOrSparseMatrix:
    """Convert a lag matrix into a recurrence matrix.

    Parameters
    ----------
    lag : np.ndarray or scipy.sparse.spmatrix
        A lag matrix, as produced by ``recurrence_to_lag``
    axis : int
        The axis corresponding to the time dimension.
        The alternate axis will be interpreted in lag coordinates.

    Returns
    -------
    rec : np.ndarray or scipy.sparse.spmatrix [shape=(n, n)]
        A recurrence matrix in (time, time) coordinates
        For sparse matrices, format will match that of ``lag``.

    Raises
    ------
    ParameterError : if ``lag`` does not have the correct shape

    See Also
    --------
    recurrence_to_lag

    Examples
    --------
    >>> y, sr = librosa.load(librosa.ex('nutcracker'))
    >>> hop_length = 1024
    >>> chroma = librosa.feature.chroma_cqt(y=y, sr=sr, hop_length=hop_length)
    >>> chroma_stack = librosa.feature.stack_memory(chroma, n_steps=10, delay=3)
    >>> recurrence = librosa.segment.recurrence_matrix(chroma_stack)
    >>> lag_pad = librosa.segment.recurrence_to_lag(recurrence, pad=True)
    >>> lag_nopad = librosa.segment.recurrence_to_lag(recurrence, pad=False)
    >>> rec_pad = librosa.segment.lag_to_recurrence(lag_pad)
    >>> rec_nopad = librosa.segment.lag_to_recurrence(lag_nopad)

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True)
    >>> librosa.display.specshow(lag_pad, x_axis='s', y_axis='lag',
    ...                          hop_length=hop_length, ax=ax[0, 0])
    >>> ax[0, 0].set(title='Lag (zero-padded)')
    >>> ax[0, 0].label_outer()
    >>> librosa.display.specshow(lag_nopad, x_axis='s', y_axis='time',
    ...                          hop_length=hop_length, ax=ax[0, 1])
    >>> ax[0, 1].set(title='Lag (no padding)')
    >>> ax[0, 1].label_outer()
    >>> librosa.display.specshow(rec_pad, x_axis='s', y_axis='time',
    ...                          hop_length=hop_length, ax=ax[1, 0])
    >>> ax[1, 0].set(title='Recurrence (with padding)')
    >>> librosa.display.specshow(rec_nopad, x_axis='s', y_axis='time',
    ...                          hop_length=hop_length, ax=ax[1, 1])
    >>> ax[1, 1].set(title='Recurrence (without padding)')
    >>> ax[1, 1].label_outer()
    """
    if axis not in [0, 1, -1]:
        raise ParameterError(f"Invalid target axis: {axis}")

    axis = np.abs(axis)

    if lag.ndim != 2 or (
        lag.shape[0] != lag.shape[1] and lag.shape[1 - axis] != 2 * lag.shape[axis]
    ):
        raise ParameterError(f"Invalid lag matrix shape: {lag.shape}")

    # Since lag must be 2-dimensional, abs(axis) = axis
    t = lag.shape[axis]

    rec = util.shear(lag, factor=+1, axis=axis)

    sub_slice = [slice(None)] * rec.ndim
    sub_slice[1 - axis] = slice(t)
    rec_slice: _ArrayOrSparseMatrix = rec[tuple(sub_slice)]
    return rec_slice


_F = TypeVar("_F", bound=Callable[..., Any])


def timelag_filter(function: _F, pad: bool = True, index: int = 0) -> _F:
    """Apply a filter in the time-lag domain.

    This is primarily useful for adapting image filters to operate on
    `recurrence_to_lag` output.

    Using `timelag_filter` is equivalent to the following sequence of
    operations:

    >>> data_tl = librosa.segment.recurrence_to_lag(data)
    >>> data_filtered_tl = function(data_tl)
    >>> data_filtered = librosa.segment.lag_to_recurrence(data_filtered_tl)

    Parameters
    ----------
    function : callable
        The filtering function to wrap, e.g., `scipy.ndimage.median_filter`
    pad : bool
        Whether to zero-pad the structure feature matrix
    index : int >= 0
        If ``function`` accepts input data as a positional argument, it should be
        indexed by ``index``

    Returns
    -------
    wrapped_function : callable
        A new filter function which applies in time-lag space rather than
        time-time space.

    Examples
    --------
    Apply a 31-bin median filter to the diagonal of a recurrence matrix.
    With default, parameters, this corresponds to a time window of about
    0.72 seconds.

    >>> y, sr = librosa.load(librosa.ex('nutcracker'), duration=30)
    >>> chroma = librosa.feature.chroma_cqt(y=y, sr=sr)
    >>> chroma_stack = librosa.feature.stack_memory(chroma, n_steps=3, delay=3)
    >>> rec = librosa.segment.recurrence_matrix(chroma_stack)
    >>> from scipy.ndimage import median_filter
    >>> diagonal_median = librosa.segment.timelag_filter(median_filter)
    >>> rec_filtered = diagonal_median(rec, size=(1, 31), mode='mirror')

    Or with affinity weights

    >>> rec_aff = librosa.segment.recurrence_matrix(chroma_stack, mode='affinity')
    >>> rec_aff_fil = diagonal_median(rec_aff, size=(1, 31), mode='mirror')

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(nrows=2, ncols=2, sharex=True, sharey=True)
    >>> librosa.display.specshow(rec, y_axis='s', x_axis='s', ax=ax[0, 0])
    >>> ax[0, 0].set(title='Raw recurrence matrix')
    >>> ax[0, 0].label_outer()
    >>> librosa.display.specshow(rec_filtered, y_axis='s', x_axis='s', ax=ax[0, 1])
    >>> ax[0, 1].set(title='Filtered recurrence matrix')
    >>> ax[0, 1].label_outer()
    >>> librosa.display.specshow(rec_aff, x_axis='s', y_axis='s',
    ...                          cmap='magma_r', ax=ax[1, 0])
    >>> ax[1, 0].set(title='Raw affinity matrix')
    >>> librosa.display.specshow(rec_aff_fil, x_axis='s', y_axis='s',
    ...                          cmap='magma_r', ax=ax[1, 1])
    >>> ax[1, 1].set(title='Filtered affinity matrix')
    >>> ax[1, 1].label_outer()
    """
    def __my_filter(wrapped_f, *args, **kwargs):
        """Wrap the filter with lag conversions"""
        # Map the input data into time-lag space
        args = list(args)

        args[index] = recurrence_to_lag(args[index], pad=pad)

        # Apply the filtering function
        result = wrapped_f(*args, **kwargs)

        # Map back into time-time and return
        return lag_to_recurrence(result)

    return decorator(__my_filter, function)  # type: ignore


@cache(level=30)
def subsegment(
    data: np.ndarray, frames: np.ndarray, *, n_segments: int = 4, axis: int = -1
) -> np.ndarray:
    """Sub-divide a segmentation by feature clustering.

    Given a set of frame boundaries (``frames``), and a data matrix (``data``),
    each successive interval defined by ``frames`` is partitioned into
    ``n_segments`` by constrained agglomerative clustering.

    .. note::
        If an interval spans fewer than ``n_segments`` frames, then each
        frame becomes a sub-segment.

    Parameters
    ----------
    data : np.ndarray
        Data matrix to use in clustering
    frames : np.ndarray [shape=(n_boundaries,)], dtype=int, non-negative]
        Array of beat or segment boundaries, as provided by
        `librosa.beat.beat_track`,
        `librosa.onset.onset_detect`,
        or `agglomerative`.
    n_segments : int > 0
        Maximum number of frames to sub-divide each interval.
    axis : int
        Axis along which to apply the segmentation.
        By default, the last index (-1) is taken.

    Returns
    -------
    boundaries : np.ndarray [shape=(n_subboundaries,)]
        List of sub-divided segment boundaries

    See Also
    --------
    agglomerative : Temporal segmentation
    librosa.onset.onset_detect : Onset detection
    librosa.beat.beat_track : Beat tracking

    Notes
    -----
    This function caches at level 30.

    Examples
    --------
    Load audio, detect beat frames, and subdivide in twos by CQT

    >>> y, sr = librosa.load(librosa.ex('choice'), duration=10)
    >>> tempo, beats = librosa.beat.beat_track(y=y, sr=sr, hop_length=512)
    >>> beat_times = librosa.frames_to_time(beats, sr=sr, hop_length=512)
    >>> cqt = np.abs(librosa.cqt(y, sr=sr, hop_length=512))
    >>> subseg = librosa.segment.subsegment(cqt, beats, n_segments=2)
    >>> subseg_t = librosa.frames_to_time(subseg, sr=sr, hop_length=512)

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots()
    >>> librosa.display.specshow(librosa.amplitude_to_db(cqt,
    ...                                                  ref=np.max),
    ...                          y_axis='cqt_hz', x_axis='time', ax=ax)
    >>> lims = ax.get_ylim()
    >>> ax.vlines(beat_times, lims[0], lims[1], color='lime', alpha=0.9,
    ...            linewidth=2, label='Beats')
    >>> ax.vlines(subseg_t, lims[0], lims[1], color='linen', linestyle='--',
    ...            linewidth=1.5, alpha=0.5, label='Sub-beats')
    >>> ax.legend()
    >>> ax.set(title='CQT + Beat and sub-beat markers')
    """
    frames = util.fix_frames(frames, x_min=0, x_max=data.shape[axis], pad=True)

    if n_segments < 1:
        raise ParameterError("n_segments must be a positive integer")

    boundaries = []
    idx_slices = [slice(None)] * data.ndim

    for seg_start, seg_end in zip(frames[:-1], frames[1:]):
        idx_slices[axis] = slice(seg_start, seg_end)
        boundaries.extend(
            seg_start
            + agglomerative(
                data[tuple(idx_slices)], min(seg_end - seg_start, n_segments), axis=axis
            )
        )

    return np.array(boundaries)


def agglomerative(
    data: np.ndarray,
    k: int,
    *,
    clusterer: Optional[sklearn.cluster.AgglomerativeClustering] = None,
    axis: int = -1,
) -> np.ndarray:
    """Bottom-up temporal segmentation.

    Use a temporally-constrained agglomerative clustering routine to partition
    ``data`` into ``k`` contiguous segments.

    Parameters
    ----------
    data : np.ndarray
        data to cluster
    k : int > 0 [scalar]
        number of segments to produce
    clusterer : sklearn.cluster.AgglomerativeClustering, optional
        An optional AgglomerativeClustering object.
        If `None`, a constrained Ward object is instantiated.
    axis : int
        axis along which to cluster.
        By default, the last axis (-1) is chosen.

    Returns
    -------
    boundaries : np.ndarray [shape=(k,)]
        left-boundaries (frame numbers) of detected segments. This
        will always include `0` as the first left-boundary.

    See Also
    --------
    sklearn.cluster.AgglomerativeClustering

    Examples
    --------
    Cluster by chroma similarity, break into 20 segments

    >>> y, sr = librosa.load(librosa.ex('nutcracker'), duration=15)
    >>> chroma = librosa.feature.chroma_cqt(y=y, sr=sr)
    >>> bounds = librosa.segment.agglomerative(chroma, 20)
    >>> bound_times = librosa.frames_to_time(bounds, sr=sr)
    >>> bound_times
    array([ 0.   ,  0.65 ,  1.091,  1.927,  2.438,  2.902,  3.924,
            4.783,  5.294,  5.712,  6.13 ,  7.314,  8.522,  8.916,
            9.66 , 10.844, 11.238, 12.028, 12.492, 14.095])

    Plot the segmentation over the chromagram

    >>> import matplotlib.pyplot as plt
    >>> import matplotlib.transforms as mpt
    >>> fig, ax = plt.subplots()
    >>> trans = mpt.blended_transform_factory(
    ...             ax.transData, ax.transAxes)
    >>> librosa.display.specshow(chroma, y_axis='chroma', x_axis='time', ax=ax)
    >>> ax.vlines(bound_times, 0, 1, color='linen', linestyle='--',
    ...           linewidth=2, alpha=0.9, label='Segment boundaries',
    ...           transform=trans)
    >>> ax.legend()
    >>> ax.set(title='Power spectrogram')
    """
    # Make sure we have at least two dimensions
    data = np.atleast_2d(data)

    # Swap data index to position 0
    data = np.swapaxes(data, axis, 0)

    # Flatten the features
    n = data.shape[0]
    data = data.reshape((n, -1), order="F")

    if clusterer is None:
        # Connect the temporal connectivity graph
        grid = sklearn.feature_extraction.image.grid_to_graph(n_x=n, n_y=1, n_z=1)

        # Instantiate the clustering object
        clusterer = sklearn.cluster.AgglomerativeClustering(
            n_clusters=k, connectivity=grid, memory=cache.memory
        )

    # Fit the model
    clusterer.fit(data)

    # Find the change points from the labels
    boundaries = [0]
    boundaries.extend(list(1 + np.nonzero(np.diff(clusterer.labels_))[0].astype(int)))
    return np.asarray(boundaries)


def path_enhance(
    R: np.ndarray,
    n: int,
    *,
    window: _WindowSpec = "hann",
    max_ratio: float = 2.0,
    min_ratio: Optional[float] = None,
    n_filters: int = 7,
    zero_mean: bool = False,
    clip: bool = True,
    **kwargs: Any,
) -> np.ndarray:
    """Multi-angle path enhancement for self- and cross-similarity matrices.

    This function convolves multiple diagonal smoothing filters with a self-similarity (or
    recurrence) matrix R, and aggregates the result by an element-wise maximum.

    Technically, the output is a matrix R_smooth such that::

        R_smooth[i, j] = max_theta (R * filter_theta)[i, j]

    where `*` denotes 2-dimensional convolution, and ``filter_theta`` is a smoothing filter at
    orientation theta.

    This is intended to provide coherent temporal smoothing of self-similarity matrices
    when there are changes in tempo.

    Smoothing filters are generated at evenly spaced orientations between min_ratio and
    max_ratio.

    This function is inspired by the multi-angle path enhancement of [#]_, but differs by
    modeling tempo differences in the space of similarity matrices rather than re-sampling
    the underlying features prior to generating the self-similarity matrix.

    .. [#] Müller, Meinard and Frank Kurth.
            "Enhancing similarity matrices for music audio analysis."
            2006 IEEE International Conference on Acoustics Speech and Signal Processing Proceedings.
            Vol. 5. IEEE, 2006.

    .. note:: if using recurrence_matrix to construct the input similarity matrix, be sure to include the main
              diagonal by setting ``self=True``.  Otherwise, the diagonal will be suppressed, and this is likely to
              produce discontinuities which will pollute the smoothing filter response.

    Parameters
    ----------
    R : np.ndarray
        The self- or cross-similarity matrix to be smoothed.
        Note: sparse inputs are not supported.

        If the recurrence matrix is multi-dimensional, e.g. `shape=(c, n, n)`,
        then enhancement is conducted independently for each leading channel.

    n : int > 0
        The length of the smoothing filter

    window : window specification
        The type of smoothing filter to use.  See `filters.get_window` for more information
        on window specification formats.

    max_ratio : float > 0
        The maximum tempo ratio to support

    min_ratio : float > 0
        The minimum tempo ratio to support.
        If not provided, it will default to ``1/max_ratio``

    n_filters : int >= 1
        The number of different smoothing filters to use, evenly spaced
        between ``min_ratio`` and ``max_ratio``.

        If ``min_ratio = 1/max_ratio`` (the default), using an odd number
        of filters will ensure that the main diagonal (ratio=1) is included.

    zero_mean : bool
        By default, the smoothing filters are non-negative and sum to one (i.e. are averaging
        filters).

        If ``zero_mean=True``, then the smoothing filters are made to sum to zero by subtracting
        a constant value from the non-diagonal coordinates of the filter.  This is primarily
        useful for suppressing blocks while enhancing diagonals.

    clip : bool
        If True, the smoothed similarity matrix will be thresholded at 0, and will not contain
        negative entries.

    **kwargs : additional keyword arguments
        Additional arguments to pass to `scipy.ndimage.convolve`

    Returns
    -------
    R_smooth : np.ndarray, shape=R.shape
        The smoothed self- or cross-similarity matrix

    See Also
    --------
    librosa.filters.diagonal_filter
    recurrence_matrix

    Examples
    --------
    Use a 51-frame diagonal smoothing filter to enhance paths in a recurrence matrix

    >>> y, sr = librosa.load(librosa.ex('nutcracker'))
    >>> hop_length = 2048
    >>> chroma = librosa.feature.chroma_cqt(y=y, sr=sr, hop_length=hop_length)
    >>> chroma_stack = librosa.feature.stack_memory(chroma, n_steps=10, delay=3)
    >>> rec = librosa.segment.recurrence_matrix(chroma_stack, mode='affinity', self=True)
    >>> rec_smooth = librosa.segment.path_enhance(rec, 51, window='hann', n_filters=7)

    Plot the recurrence matrix before and after smoothing

    >>> import matplotlib.pyplot as plt
    >>> fig, ax = plt.subplots(ncols=2, sharex=True, sharey=True)
    >>> img = librosa.display.specshow(rec, x_axis='s', y_axis='s',
    ...                          hop_length=hop_length, ax=ax[0])
    >>> ax[0].set(title='Unfiltered recurrence')
    >>> imgpe = librosa.display.specshow(rec_smooth, x_axis='s', y_axis='s',
    ...                          hop_length=hop_length, ax=ax[1])
    >>> ax[1].set(title='Multi-angle enhanced recurrence')
    >>> ax[1].label_outer()
    >>> fig.colorbar(img, ax=ax[0], orientation='horizontal')
    >>> fig.colorbar(imgpe, ax=ax[1], orientation='horizontal')
    """
    if min_ratio is None:
        min_ratio = 1.0 / max_ratio
    elif min_ratio > max_ratio:
        raise ParameterError(
            f"min_ratio={min_ratio} cannot exceed max_ratio={max_ratio}"
        )

    R_smooth = None
    for ratio in np.logspace(
        np.log2(min_ratio), np.log2(max_ratio), num=n_filters, base=2
    ):
        kernel = diagonal_filter(window, n, slope=ratio, zero_mean=zero_mean)

        # Expand leading dimensions to match R
        # This way, if R has shape, eg, [2, 3, n, n]
        # the expanded kernel will have shape [1, 1, m, m]

        # The following is valid for numpy >= 1.18
        # kernel = np.expand_dims(kernel, axis=list(np.arange(R.ndim - kernel.ndim)))

        # This is functionally equivalent, but works on numpy 1.17
        shape = [1] * R.ndim
        shape[-2:] = kernel.shape
        kernel = np.reshape(kernel, shape)

        if R_smooth is None:
            R_smooth = scipy.ndimage.convolve(R, kernel, **kwargs)
        else:
            # Compute the point-wise maximum in-place
            np.maximum(
                R_smooth, scipy.ndimage.convolve(R, kernel, **kwargs), out=R_smooth
            )

    if clip:
        # Clip the output in-place
        np.clip(R_smooth, 0, None, out=R_smooth)  # type: ignore

    return np.asanyarray(R_smooth)


def __affinity_bandwidth(
    rec: scipy.sparse.csr_matrix,
    bw_mode: Optional[Union[np.ndarray, _FloatLike_co, str]],
    k: int,
) -> Union[float, np.ndarray]:
    # rec should be a csr_matrix

    # the api allows users to specify a scalar bandwidth directly, besides the string based options.
    if isinstance(bw_mode, np.ndarray):
        bandwidth = bw_mode
        # check if bw is the right size
        if bandwidth.shape != rec.shape:
            raise ParameterError(
                f"Invalid matrix bandwidth shape: {bandwidth.shape}."
                f"Should be {rec.shape}."
            )
        if (bandwidth <= 0).any():
            raise ParameterError(
                "Invalid bandwidth. All entries must be strictly positive."
            )
        return np.array(bandwidth[rec.nonzero()])

    elif isinstance(bw_mode, (int, float)):
        scalar_bandwidth = float(bw_mode)
        if scalar_bandwidth <= 0:
            raise ParameterError(
                f"Invalid scalar bandwidth={scalar_bandwidth}. Must be strictly positive."
            )
        return scalar_bandwidth

    if bw_mode is None:
        bw_mode = "med_k_scalar"

    if bw_mode not in [
        "med_k_scalar",
        "mean_k",
        "gmean_k",
        "mean_k_avg",
        "gmean_k_avg",
        "mean_k_avg_and_pair",
    ]:
        raise ParameterError(
            f"Invalid bandwidth='{bw_mode}'. Must be either a positive scalar or one of "
            "['med_k_scalar', 'mean_k', 'gmean_k', 'mean_k_avg', 'gmean_k_avg', 'mean_k_avg_and_pair']"
        )

    # build a list of list that stores the distances to k nearest neighbors for all t points.
    t = rec.shape[0]
    knn_dists = []
    for i in range(t):
        # Get the links from point i
        links = rec[i].nonzero()[1]
        # catch empty dists lists in knn_dists
        if len(links) == 0:
            # Disconnected vertices are only a problem for point-wise bandwidth estimation
            if bw_mode not in ["med_k_scalar"]:
                raise ParameterError(f"The sample at time point {i} has no neighbors")
            else:
                # If we have no links, then there's no distance
                # shove a nan in here
                knn_dists.append(np.array([np.nan]))
        else:
            # Compute k nearest neighbors' distance and sort ascending
            knn_dist_row = np.sort(rec[i, links].toarray()[0])[:k]
            knn_dists.append(knn_dist_row)

    # take the last element of each list for the distance to kth neighbor
    dist_to_k = np.asarray([dists[-1] for dists in knn_dists])
    avg_dist_to_first_ks = np.asarray([np.mean(dists) for dists in knn_dists])

    if bw_mode == "med_k_scalar":
        if not np.any(np.isfinite(dist_to_k)):
            raise ParameterError("Cannot estimate bandwidth from an empty graph")
        return float(np.nanmedian(dist_to_k))

    if bw_mode in ["mean_k", "gmean_k"]:
        # building bandwidth components (sigma) using sparse matrix structures and indices
        sigma_i_data = np.empty_like(rec.data)
        sigma_j_data = np.empty_like(rec.data)
        for row in range(t):
            sigma_i_data[rec.indptr[row] : rec.indptr[row + 1]] = dist_to_k[row]
            col_idx = rec.indices[rec.indptr[row] : rec.indptr[row + 1]]
            sigma_j_data[rec.indptr[row] : rec.indptr[row + 1]] = dist_to_k[col_idx]

        if bw_mode == "mean_k":
            out = np.array((sigma_i_data + sigma_j_data) / 2)
        elif bw_mode == "gmean_k":
            out = np.array((sigma_i_data * sigma_j_data) ** 0.5)

    if bw_mode in ["mean_k_avg", "gmean_k_avg", "mean_k_avg_and_pair"]:
        # building bandwidth components (sigma) using sparse matrix structures and indices
        sigma_i_data = np.empty_like(rec.data)
        sigma_j_data = np.empty_like(rec.data)
        for row in range(t):
            sigma_i_data[rec.indptr[row] : rec.indptr[row + 1]] = avg_dist_to_first_ks[
                row
            ]
            col_idx = rec.indices[rec.indptr[row] : rec.indptr[row + 1]]
            sigma_j_data[rec.indptr[row] : rec.indptr[row + 1]] = avg_dist_to_first_ks[
                col_idx
            ]

        if bw_mode == "mean_k_avg":
            out = np.array((sigma_i_data + sigma_j_data) / 2)
        elif bw_mode == "gmean_k_avg":
            out = np.array((sigma_i_data * sigma_j_data) ** 0.5)
        elif bw_mode == "mean_k_avg_and_pair":
            out = np.array((sigma_i_data + sigma_j_data + rec.data) / 3)

    return out