API Reference

Auto-generated from the docstrings in the source. See Concepts & coming from MATLAB for the MATLAB function-name mapping.

Core pipeline

tmapper.tknndigraph

Port of tmapper_tools/tknndigraph.m from the MATLAB toolbox.

tknndigraph(X_or_D, k, tidx, *, reciprocal=True, time_exclude_space=True, time_exclude_range=1, max_neighbor_dist=np.inf, max_neighbor_dist_prct=100.0)

Construct a directed graph based on k-nearest neighbors that also includes temporal neighbors (t -> t+1 links).

Port of MATLAB's tknndigraph.m. Node i of the returned graph corresponds to row i of X_or_D / element i of tidx (0-indexed, unlike the 1-indexed MATLAB original).

Parameters:

Name	Type	Description	Default
X_or_D	(array_like, shape(N, d) or (N, N))	Either raw coordinates (N points in d-dim space) or a precomputed (N, N) distance matrix. Auto-detected: treated as a distance matrix only if square and symmetric; otherwise treated as raw coordinates and converted to a Euclidean distance matrix.	required
k	int	Max number of spatial neighbors per point. Must be a positive integer strictly less than N.	required
tidx	array_like of int, shape (N,)	Time index per point. Two points x, y are temporal neighbors iff tidx[y] == tidx[x] + 1.	required
reciprocal	bool	Whether to require spatial k-NN links to be mutual (both directions) to be kept.	True
time_exclude_space	bool	Whether temporal neighbors are barred from also counting as spatial neighbors.	True
time_exclude_range	int	How many following time points are excluded from spatial-kNN consideration.	1
max_neighbor_dist	float	Absolute distance cutoff for spatial neighbors.	inf
max_neighbor_dist_prct	float	Percentile distance cutoff for spatial neighbors. The stricter (smaller) of this and max_neighbor_dist is applied.	100

Returns:

Name	Type	Description
g	DiGraph	Unweighted directed graph, one node per input point.
par	dict	Parameters actually used, including the resolved max_neighbor_dist.

Source code in src/tmapper/tknndigraph.py

def tknndigraph(
    X_or_D,
    k,
    tidx,
    *,
    reciprocal=True,
    time_exclude_space=True,
    time_exclude_range=1,
    max_neighbor_dist=np.inf,
    max_neighbor_dist_prct=100.0,
):
    """Construct a directed graph based on k-nearest neighbors that also
    includes temporal neighbors (t -> t+1 links).

    Port of MATLAB's ``tknndigraph.m``. Node i of the returned graph
    corresponds to row i of ``X_or_D`` / element i of ``tidx`` (0-indexed,
    unlike the 1-indexed MATLAB original).

    Parameters
    ----------
    X_or_D : array_like, shape (N, d) or (N, N)
        Either raw coordinates (N points in d-dim space) or a precomputed
        (N, N) distance matrix. Auto-detected: treated as a distance
        matrix only if square *and* symmetric; otherwise treated as raw
        coordinates and converted to a Euclidean distance matrix.
    k : int
        Max number of spatial neighbors per point. Must be a positive
        integer strictly less than N.
    tidx : array_like of int, shape (N,)
        Time index per point. Two points x, y are temporal neighbors iff
        ``tidx[y] == tidx[x] + 1``.
    reciprocal : bool, default True
        Whether to require spatial k-NN links to be mutual (both
        directions) to be kept.
    time_exclude_space : bool, default True
        Whether temporal neighbors are barred from also counting as
        spatial neighbors.
    time_exclude_range : int, default 1
        How many following time points are excluded from spatial-kNN
        consideration.
    max_neighbor_dist : float, default inf
        Absolute distance cutoff for spatial neighbors.
    max_neighbor_dist_prct : float, default 100
        Percentile distance cutoff for spatial neighbors. The stricter
        (smaller) of this and ``max_neighbor_dist`` is applied.

    Returns
    -------
    g : networkx.DiGraph
        Unweighted directed graph, one node per input point.
    par : dict
        Parameters actually used, including the resolved
        ``max_neighbor_dist``.
    """
    X_or_D = np.asarray(X_or_D, dtype=float)
    if X_or_D.ndim != 2:
        raise ValueError("X_or_D must be a 2D array.")
    if not (np.isscalar(k) and float(k) == round(float(k)) and k >= 1):
        raise ValueError("k must be a positive integer scalar.")
    k = int(k)

    tidx = np.asarray(tidx, dtype=float).ravel()
    if tidx.shape[0] != X_or_D.shape[0]:
        raise ValueError(
            "tidx must have the same number of elements as rows of X_or_D."
        )

    par = {
        "reciprocal": reciprocal,
        "time_exclude_space": time_exclude_space,
        "time_exclude_range": time_exclude_range,
        "max_neighbor_dist": max_neighbor_dist,
        "max_neighbor_dist_prct": max_neighbor_dist_prct,
        "k": k,
    }

    # -- check input and obtain distance matrix D
    nr, nc = X_or_D.shape
    if nr != nc or not np.allclose(X_or_D, X_or_D.T):
        D = cdist(X_or_D, X_or_D)
    else:
        D = X_or_D.copy()
    Nn = D.shape[0]

    if k >= Nn:
        raise ValueError(f"k must be smaller than the number of points ({Nn}).")

    np.fill_diagonal(D, np.inf)  # exclude self-loops

    # -- find indices for temporal links D[i(t), i(t+1)]
    # t_wafter[i] is True iff point i has a point immediately after it in time
    # (mirrors MATLAB's circshift(tidx,-1,1) - 1 == tidx, including its
    # wraparound behavior at the last index, which is virtually always False
    # for non-cyclic tidx and is harmless here for the same reason it is in
    # the original).
    tidx_next = np.roll(tidx, -1)
    t_wafter = tidx_next - 1 == tidx

    # t_after_idx1: unconditional immediate temporal edges i -> i+1
    t_after_idx1 = np.zeros((Nn, Nn), dtype=bool)
    for i in range(Nn):
        if t_wafter[i]:
            t_after_idx1[i, (i + 1) % Nn] = True

    # t_after_idx: temporal-exclusion mask, extended up to time_exclude_range
    # steps, restricted to strictly-forward (i < j) positions
    t_after_idx = np.zeros((Nn, Nn), dtype=bool)
    for n in range(1, time_exclude_range + 1):
        for i in range(Nn):
            if t_wafter[i]:
                t_after_idx[i, (i + n) % Nn] = True
    upper = np.triu(np.ones((Nn, Nn), dtype=bool), k=1)
    t_after_idx = t_after_idx & upper

    if time_exclude_space:
        D[t_after_idx] = np.inf

    # -- compute adjacency matrix (k nearest neighbors per row)
    order = np.argsort(D, axis=1, kind="stable")
    D_sorted = np.take_along_axis(D, order, axis=1)
    A = np.zeros((Nn, Nn), dtype=bool)
    rows = np.repeat(np.arange(Nn), k)
    cols = order[:, :k].ravel()
    A[rows, cols] = True

    # -- check for duplicate points: any other point tied with the k-th
    # nearest neighbor's distance is also included
    dmax = D_sorted[:, k - 1][:, None]
    A |= D <= dmax

    # -- get distance threshold (the stricter of the absolute and percentile caps).
    # NOTE: matches MATLAB's prctile(D(:), Prct) exactly by computing the
    # percentile over the *full* D, Inf entries included (the masked
    # diagonal/temporal-exclusion entries) -- not just the finite values.
    prct_dist = _percentile_with_inf(D, max_neighbor_dist_prct)
    resolved_max_dist = min(prct_dist, max_neighbor_dist)
    par["max_neighbor_dist"] = resolved_max_dist

    # -- remove neighbors that exceed max distance
    A[D > resolved_max_dist] = False

    # -- exclude or retain temporal links as spatial links
    if time_exclude_space:
        A_space = A & ~t_after_idx
    else:
        A_space = A.copy()

    # -- enforce symmetry of spatial links
    if reciprocal:
        A_space = A_space & A_space.T
    else:
        A_space = A_space | A_space.T

    # -- (re-)incorporate temporal links
    A_final = t_after_idx1 | A_space

    g = nx.from_numpy_array(A_final, create_using=nx.DiGraph)
    return g, par

tmapper.filtergraph

Port of tmapper_tools/filtergraph.m from the MATLAB toolbox.

filtergraph(g, d, *, reciprocal=True)

Contract a graph g into a simplified graph (a Mapper-style shape graph) by merging nodes that are within geodesic distance d of one another.

Port of MATLAB's filtergraph.m. This is the Mapper-style contraction step: nodes of g within geodesic distance d of each other are connected in an intermediate graph, and its connected components become the nodes of the simplified graph.

Parameters:

Name	Type	Description	Default
g	Graph or DiGraph	Graph to simplify.	required
d	float	Threshold under which original nodes are collapsed together. Must be a positive real number.	required
reciprocal	bool	Whether to require both the path length from x to y and from y to x to be under d (True), or either direction (False).	True

Returns:

Name	Type	Description
g_simp	DiGraph	The simplified graph -- the attractor transition network.
members	list of list	members[n] is the list of original node labels (from g) belonging to new node n. If g's nodes are the default 0..N-1 integers straight from :func:tknndigraph, these are positions into your original data, not real time labels.
nodesize	(ndarray, shape(n_new))	Number of original-graph members in each new node.
D_simp	(ndarray, shape(n_new, n_new))	Shortest "distance" (from the original graph) between each pair of new nodes' member sets.

Source code in src/tmapper/filtergraph.py

def filtergraph(g, d, *, reciprocal=True):
    """Contract a graph ``g`` into a simplified graph (a Mapper-style shape
    graph) by merging nodes that are within geodesic distance ``d`` of one
    another.

    Port of MATLAB's ``filtergraph.m``. This is the Mapper-style contraction
    step: nodes of ``g`` within geodesic distance ``d`` of each other are
    connected in an intermediate graph, and its connected components become
    the nodes of the simplified graph.

    Parameters
    ----------
    g : networkx.Graph or networkx.DiGraph
        Graph to simplify.
    d : float
        Threshold under which original nodes are collapsed together. Must
        be a positive real number.
    reciprocal : bool, default True
        Whether to require both the path length from x to y *and* from y
        to x to be under ``d`` (True), or either direction (False).

    Returns
    -------
    g_simp : networkx.DiGraph
        The simplified graph -- the attractor transition network.
    members : list of list
        ``members[n]`` is the list of original node labels (from ``g``)
        belonging to new node n. If ``g``'s nodes are the default 0..N-1
        integers straight from :func:`tknndigraph`, these are *positions*
        into your original data, not real time labels.
    nodesize : numpy.ndarray, shape (n_new,)
        Number of original-graph members in each new node.
    D_simp : numpy.ndarray, shape (n_new, n_new)
        Shortest "distance" (from the original graph) between each pair
        of new nodes' member sets.
    """
    if not isinstance(g, (nx.Graph, nx.DiGraph)):
        raise ValueError("g must be a networkx Graph or DiGraph.")
    if not (np.isscalar(d) and np.isreal(d) and d > 0):
        raise ValueError("d must be a positive real scalar.")

    nodelist = list(g.nodes())
    Nn = len(nodelist)
    idx_of = {node: i for i, node in enumerate(nodelist)}

    A = nx.to_numpy_array(g, nodelist=nodelist)  # weighted adjacency (weight 1 if unweighted)
    D = all_pairs_distance(g, nodelist)  # geodesic distance, inf if unreachable

    # -- connectivity within a distance threshold
    if reciprocal:
        A_ = (D < d) & (D.T < d)
    else:
        A_ = (D < d) | (D.T < d)
    np.fill_diagonal(A_, False)  # remove self-loops

    # -- create graph out of nodes within said threshold, find connected components
    g_ = nx.from_numpy_array(A_)
    components = list(nx.connected_components(g_))
    # sort components by their smallest member index, for deterministic/stable ordering
    components = sorted(components, key=min)

    idx_newnodes = np.empty(Nn, dtype=int)
    for new_idx, comp in enumerate(components):
        for i in comp:
            idx_newnodes[i] = new_idx
    n_new = len(components)

    # -- group original nodes by new-node label via a stable sort, so each
    # group is a contiguous run. This lets the block min/mean below use
    # ufunc.reduceat instead of an O(n_new^2) Python loop over Nn-sized
    # boolean masks (the previous approach was O(n_new^2 * Nn) and became
    # the dominant cost at realistic graph sizes).
    order = np.argsort(idx_newnodes, kind="stable")
    group_start = np.searchsorted(idx_newnodes[order], np.arange(n_new))
    group_bounds = np.append(group_start, Nn)

    members = [
        [nodelist[i] for i in order[group_bounds[n]:group_bounds[n + 1]]]
        for n in range(n_new)
    ]
    nodesize = np.array([len(m) for m in members])

    D_sorted = D[np.ix_(order, order)]
    A_sorted = A[np.ix_(order, order)]

    # -- define distance between new nodes: shortest path between member sets
    D_row = np.minimum.reduceat(D_sorted, group_start, axis=0)
    D_simp = np.minimum.reduceat(D_row, group_start, axis=1)

    # -- construct simplified graph: A_simp(n,m) = average connectivity between blocks
    A_row_sum = np.add.reduceat(A_sorted, group_start, axis=0)
    A_col_sum = np.add.reduceat(A_row_sum, group_start, axis=1)
    group_sizes = np.diff(group_bounds)
    A_simp = A_col_sum / np.outer(group_sizes, group_sizes)

    g_simp = nx.from_numpy_array(A_simp, create_using=nx.DiGraph)
    g_simp.remove_edges_from(nx.selfloop_edges(g_simp))  # OmitSelfLoops

    return g_simp, members, nodesize, D_simp

tmapper.find_node_label(members, x_label, *, labelmethod='mode')

Aggregate a per-time-point label to a per-node label.

Port of MATLAB's findnodelabel.m.

Parameters:

Name	Type	Description	Default
members	sequence of sequence of int	members[n] gives the positional indices belonging to node n.	required
x_label	array_like	Label/value for each time point (indexed the same way as the indices inside members).	required
labelmethod	(mode, mean, median, none)	How to aggregate each node's members' x_label values. A callable is applied as labelmethod(x_label[members[n]]) per node and must return a scalar. 'none' gives every node the same label (0).	'mode'

Returns:

Type	Description
(ndarray, shape(len(members)))	The aggregated label for each node.

Source code in src/tmapper/labeling.py

def find_node_label(members, x_label, *, labelmethod="mode"):
    """Aggregate a per-time-point label to a per-node label.

    Port of MATLAB's ``findnodelabel.m``.

    Parameters
    ----------
    members : sequence of sequence of int
        ``members[n]`` gives the positional indices belonging to node n.
    x_label : array_like
        Label/value for each time point (indexed the same way as the
        indices inside ``members``).
    labelmethod : {'mode', 'mean', 'median', 'none'} or callable, default 'mode'
        How to aggregate each node's members' x_label values. A callable
        is applied as ``labelmethod(x_label[members[n]])`` per node and
        must return a scalar. 'none' gives every node the same label (0).

    Returns
    -------
    numpy.ndarray, shape (len(members),)
        The aggregated label for each node.
    """
    x_label = np.asarray(x_label)

    if callable(labelmethod):
        return np.array([labelmethod(x_label[list(m)]) for m in members])

    if labelmethod == "mode":
        # scipy.stats.mode returns the smallest value among ties, matching
        # MATLAB's mode().
        return np.array([stats.mode(x_label[list(m)], keepdims=False).mode for m in members])
    elif labelmethod == "mean":
        return np.array([np.mean(x_label[list(m)]) for m in members])
    elif labelmethod == "median":
        return np.array([np.median(x_label[list(m)]) for m in members])
    elif labelmethod == "none":
        return np.zeros(len(members))
    else:
        raise ValueError(f"Unknown labelmethod: {labelmethod!r}")

tmapper.tcm_distance

Port of tmapper_tools/TCMdistance.m from the MATLAB toolbox.

tcm_distance(g, nodet, weighted=False)

Compute a temporal connectivity matrix: for every pair of original time points, the shortest path length between the nodes (of a simplified graph g) that contain them.

Port of MATLAB's TCMdistance.m. Sizes its output by the full range covered by nodet (max - min + 1), not just the count of distinct points referenced -- so time points genuinely not covered by any node in nodet correctly stay NaN rather than being silently zero-filled.

Parameters:

Name	Type	Description	Default
g	Graph or DiGraph	The simplified graph.	required
nodet	sequence of sequence of int	nodet[i] gives the original time-point indices belonging to node i of g (in the same node order as g.nodes()).	required
weighted	bool	Whether to use g's edge weights. If False, every edge is treated as weight 1 (hop count).	False

Returns:

Type	Description
(ndarray, shape(Nt, Nt))	The temporal connectivity matrix, where Nt = max(nodet) - min(nodet) + 1. Entries for time points never covered by any node stay NaN.

Source code in src/tmapper/tcm_distance.py

def tcm_distance(g, nodet, weighted=False):
    """Compute a temporal connectivity matrix: for every pair of original
    time points, the shortest path length between the nodes (of a
    simplified graph ``g``) that contain them.

    Port of MATLAB's ``TCMdistance.m``. Sizes its output by the full range
    covered by ``nodet`` (``max - min + 1``), not just the count of
    distinct points referenced -- so time points genuinely not covered by
    any node in ``nodet`` correctly stay ``NaN`` rather than being silently
    zero-filled.

    Parameters
    ----------
    g : networkx.Graph or networkx.DiGraph
        The simplified graph.
    nodet : sequence of sequence of int
        ``nodet[i]`` gives the original time-point indices belonging to
        node i of ``g`` (in the same node order as ``g.nodes()``).
    weighted : bool, default False
        Whether to use ``g``'s edge weights. If False, every edge is
        treated as weight 1 (hop count).

    Returns
    -------
    numpy.ndarray, shape (Nt, Nt)
        The temporal connectivity matrix, where Nt = max(nodet) - min(nodet) + 1.
        Entries for time points never covered by any node stay NaN.
    """
    nodelist = list(g.nodes())
    all_t = np.concatenate([np.asarray(list(nt), dtype=int) for nt in nodet]) if nodet else np.array([], dtype=int)
    if all_t.size == 0:
        return np.zeros((0, 0))
    t_0 = all_t.min()
    Nt = all_t.max() - t_0 + 1

    distmat = all_pairs_distance(g, nodelist, weight="weight" if weighted else None)

    tcm = np.full((Nt, Nt), np.nan)

    n_nodes = len(nodelist)
    idx_lists = [np.asarray(list(nt), dtype=int) - t_0 for nt in nodet]

    for i in range(n_nodes):
        ii = idx_lists[i]
        tcm[np.ix_(ii, ii)] = 0

    for i in range(n_nodes):
        for j in range(i + 1, n_nodes):
            ii, jj = idx_lists[i], idx_lists[j]
            block_ij = tcm[np.ix_(ii, jj)]
            tcm[np.ix_(ii, jj)] = np.fmin(block_ij, distmat[i, j])
            block_ji = tcm[np.ix_(jj, ii)]
            tcm[np.ix_(jj, ii)] = np.fmin(block_ji, distmat[j, i])

    return tcm

tmapper.plot_tmgraph(g, x_label=None, nodemembers=None, *, ax=None, nodesizerange=(1, 10), nodesizemode='log', colorlabel='x_label', cmap='jet', labelmethod='mode', nodeclim=None, nodescatter=False, center_expand=4.0)

Plot a temporal mapper graph (without a recurrence plot).

Port of MATLAB's plottmgraph.m.

Parameters:

Name	Type	Description	Default
g	Graph or DiGraph	The graph to plot.	required
x_label	array_like	A label for each member of each node, indexed the same way as the indices inside nodemembers. Defaults to a constant array of ones sized to the number of unique members across all nodes.	None
nodemembers	sequence of sequence of int	nodemembers[n] gives the indices belonging to node n. Defaults to one singleton member per node (0..N-1).	None
ax	Axes	Axes to draw into. A new figure/axes is created if not given.	None
nodesizerange	(float, float)	(min, max) marker size.	(1, 10)
nodesizemode	(log, rank, original)	How to transform node sizes before rescaling to nodesizerange.	'log'
colorlabel	str	Colorbar label.	'x_label'
cmap	str or Colormap	Colormap for node coloring.	'jet'
labelmethod	(mode, mean, median, none)	See :func:tmapper.labeling.find_node_label.	'mode'
nodeclim	(float, float)	Color axis limits. Defaults to (min(x_label), max(x_label)).	None
nodescatter	bool	Whether to overlay a scatter plot on top of the graph nodes.	False
center_expand	float	Log-radial de-clumping strength applied to the layout after normalizing it to the unit circle (0 disables it). Force-directed layouts tend to pack most nodes into a dense central clump with a few outliers stretching the frame; this spreads the dense core outward without disturbing the outer structure. See :func:_expand_center.	4.0

Returns:

Name	Type	Description
ax	Axes
cbar	Colorbar
node_collection	PathCollection	The drawn graph nodes.
scatter_collection	PathCollection or None	The scatter overlay, if nodescatter is True.

Source code in src/tmapper/plotting.py

def plot_tmgraph(
    g,
    x_label=None,
    nodemembers=None,
    *,
    ax=None,
    nodesizerange=(1, 10),
    nodesizemode="log",
    colorlabel="x_label",
    cmap="jet",
    labelmethod="mode",
    nodeclim=None,
    nodescatter=False,
    center_expand=4.0,
):
    """Plot a temporal mapper graph (without a recurrence plot).

    Port of MATLAB's ``plottmgraph.m``.

    Parameters
    ----------
    g : networkx.Graph or networkx.DiGraph
        The graph to plot.
    x_label : array_like, optional
        A label for each member of each node, indexed the same way as the
        indices inside ``nodemembers``. Defaults to a constant array of
        ones sized to the number of unique members across all nodes.
    nodemembers : sequence of sequence of int, optional
        ``nodemembers[n]`` gives the indices belonging to node n. Defaults
        to one singleton member per node (0..N-1).
    ax : matplotlib.axes.Axes, optional
        Axes to draw into. A new figure/axes is created if not given.
    nodesizerange : (float, float), default (1, 10)
        (min, max) marker size.
    nodesizemode : {'log', 'rank', 'original'}, default 'log'
        How to transform node sizes before rescaling to ``nodesizerange``.
    colorlabel : str, default 'x_label'
        Colorbar label.
    cmap : str or matplotlib.colors.Colormap, default 'jet'
        Colormap for node coloring.
    labelmethod : {'mode', 'mean', 'median', 'none'} or callable, default 'mode'
        See :func:`tmapper.labeling.find_node_label`.
    nodeclim : (float, float), optional
        Color axis limits. Defaults to (min(x_label), max(x_label)).
    nodescatter : bool, default False
        Whether to overlay a scatter plot on top of the graph nodes.
    center_expand : float, default 4.0
        Log-radial de-clumping strength applied to the layout after
        normalizing it to the unit circle (0 disables it). Force-directed
        layouts tend to pack most nodes into a dense central clump with a
        few outliers stretching the frame; this spreads the dense core
        outward without disturbing the outer structure. See
        :func:`_expand_center`.

    Returns
    -------
    ax : matplotlib.axes.Axes
    cbar : matplotlib.colorbar.Colorbar
    node_collection : matplotlib.collections.PathCollection
        The drawn graph nodes.
    scatter_collection : matplotlib.collections.PathCollection or None
        The scatter overlay, if ``nodescatter`` is True.
    """
    import matplotlib.pyplot as plt

    nodelist = list(g.nodes())
    n_nodes = len(nodelist)

    if nodemembers is None:
        nodemembers = [[i] for i in range(n_nodes)]

    if x_label is None:
        all_members = sorted({m for members in nodemembers for m in members})
        x_label = np.ones(len(all_members))
    x_label = np.asarray(x_label, dtype=float)

    if nodeclim is None:
        nodeclim = (float(np.min(x_label)), float(np.max(x_label)))

    # -- define node size and labels/colors
    nodesize = _prepare_node_size(nodemembers, nodesizerange, nodesizemode)
    nodelabel = find_node_label(nodemembers, x_label, labelmethod=labelmethod)

    # -- plotting
    if ax is None:
        _, ax = plt.subplots()

    pos = _compute_layout(g, nodelist, center_expand)

    xs = np.array([pos[n][0] for n in nodelist])
    ys = np.array([pos[n][1] for n in nodelist])

    # Edges visible but behind (zorder=1); nodes outlined and drawn in
    # front (zorder=3), smallest last, so small nodes aren't buried
    # under large ones and stay crisp regardless of edge clutter.
    edge_collection = nx.draw_networkx_edges(
        g, pos, ax=ax, alpha=0.6, edge_color="#888888", width=0.8,
        arrowsize=6, nodelist=nodelist,
    )
    for obj in (edge_collection if isinstance(edge_collection, list) else [edge_collection]):
        if obj is not None:
            obj.set_zorder(1)

    order = np.argsort(nodesize, kind="stable")[::-1]  # largest first, smallest drawn last (on top)
    node_collection = ax.scatter(
        xs[order], ys[order], s=(nodesize[order] ** 2), c=nodelabel[order], cmap=cmap,
        vmin=nodeclim[0] if nodeclim[1] != nodeclim[0] else None,
        vmax=nodeclim[1] if nodeclim[1] != nodeclim[0] else None,
        edgecolors="black", linewidths=0.5, zorder=3,
    )
    ax.set_aspect("equal")
    ax.axis("off")

    sm = plt.cm.ScalarMappable(cmap=cmap)
    sm.set_array(nodelabel)
    if nodeclim[1] != nodeclim[0]:
        sm.set_clim(nodeclim[0], nodeclim[1])
    cbar = plt.colorbar(sm, ax=ax)
    cbar.set_label(colorlabel)

    scatter_collection = None
    if nodescatter:
        order = np.argsort(nodesize, kind="stable")
        scatter_collection = ax.scatter(
            xs[order], ys[order], s=(nodesize[order] ** 2),
            c=nodelabel[order], cmap=cmap,
            vmin=nodeclim[0] if nodeclim[1] != nodeclim[0] else None,
            vmax=nodeclim[1] if nodeclim[1] != nodeclim[0] else None,
            edgecolors="k", linewidths=0.2,
        )

    return ax, cbar, node_collection, scatter_collection

tmapper.plot_tmgraph_tcm(g, x_label, t, nodemembers, **kwargs)

Plot a temporal mapper graph alongside its geodesic recurrence plot.

Port of MATLAB's plotgraphtcm.m. Calls :func:plot_tmgraph internally for the left subplot (all **kwargs are passed through), so it accepts the same styling parameters.

Parameters:

Name	Type	Description	Default
g	Graph or DiGraph	The (simplified) graph to plot.	required
x_label	array_like	A label for each time point in the time series.	required
t	array_like	Actual time (or a time index) associated with each time point, used as the recurrence plot's axis labels.	required
nodemembers	sequence of sequence of int	nodemembers[n] gives the original time-point indices belonging to node n.	required
**kwargs		Passed through to :func:plot_tmgraph.	{}

Returns:

Name	Type	Description
	ax1, ax2 : matplotlib.axes.Axes	The network plot and recurrence plot axes.
	cbar, cbar2 : matplotlib.colorbar.Colorbar
	(node_collection, scatter_collection)	As returned by :func:plot_tmgraph.
D_geo	ndarray	The recurrence plot matrix.

Source code in src/tmapper/plotting.py

def plot_tmgraph_tcm(g, x_label, t, nodemembers, **kwargs):
    """Plot a temporal mapper graph alongside its geodesic recurrence plot.

    Port of MATLAB's ``plotgraphtcm.m``. Calls :func:`plot_tmgraph`
    internally for the left subplot (all ``**kwargs`` are passed through),
    so it accepts the same styling parameters.

    Parameters
    ----------
    g : networkx.Graph or networkx.DiGraph
        The (simplified) graph to plot.
    x_label : array_like
        A label for each time point in the time series.
    t : array_like
        Actual time (or a time index) associated with each time point,
        used as the recurrence plot's axis labels.
    nodemembers : sequence of sequence of int
        ``nodemembers[n]`` gives the original time-point indices
        belonging to node n.
    **kwargs
        Passed through to :func:`plot_tmgraph`.

    Returns
    -------
    ax1, ax2 : matplotlib.axes.Axes
        The network plot and recurrence plot axes.
    cbar, cbar2 : matplotlib.colorbar.Colorbar
    node_collection, scatter_collection
        As returned by :func:`plot_tmgraph`.
    D_geo : numpy.ndarray
        The recurrence plot matrix.
    """
    import matplotlib.pyplot as plt
    import matplotlib.dates as mdates

    nodesize = np.array([len(m) for m in nodemembers])
    bsinglemember = np.all(nodesize == 1)

    # constrained_layout resolves spacing between subplots/colorbars
    # automatically, avoiding overlap between cbar's label and ax2's ylabel.
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.5), constrained_layout=True)

    ax1, cbar, node_collection, scatter_collection = plot_tmgraph(
        g, x_label, nodemembers, ax=ax1, **kwargs
    )

    if bsinglemember:
        D_geo = all_pairs_distance(g, list(g.nodes()), weight=None)
    else:
        D_geo = tcm_distance(g, nodemembers)

    t = np.asarray(t)
    is_datetime = np.issubdtype(t.dtype, np.datetime64)
    t_num = mdates.date2num(t) if is_datetime else t.astype(float)

    im = ax2.imshow(
        D_geo, cmap="hot",
        extent=[t_num.min(), t_num.max(), t_num.max(), t_num.min()],
        aspect="equal", interpolation="nearest",  # matches MATLAB's imagesc (flat, unsmoothed cells)
    )
    cbar2 = plt.colorbar(im, ax=ax2)
    cbar2.set_label("path length")

    if is_datetime:
        for axis in (ax2.xaxis, ax2.yaxis):
            locator = mdates.AutoDateLocator()
            axis.set_major_locator(locator)
            axis.set_major_formatter(mdates.ConciseDateFormatter(locator))
        ax2.set_xlabel("time")
        ax2.set_ylabel("time")
    else:
        ax2.set_xlabel("time (s)")
        ax2.set_ylabel("time (s)")
    ax2.set_title("geodesic recurrence plot")

    return ax1, ax2, cbar, cbar2, node_collection, scatter_collection, D_geo

tmapper.plot_tmgraph_interactive(g, x_label=None, nodemembers=None, *, nodesizerange=(6, 40), nodesizemode='log', colorlabel='x_label', cmap='jet', labelmethod='mode', nodeclim=None, center_expand=4.0, title='', output_path=None)

Plot a temporal mapper graph as a draggable/zoomable/hoverable standalone HTML page (via pyvis/vis.js), instead of a static image.

Uses the same layout, node-size, and coloring logic as :func:plot_tmgraph (so the two should look like the same network), with physics disabled -- nodes stay exactly where the layout put them, but can still be dragged, and hovering shows each node's size and color-label value.

Parameters:

Name	Type	Description	Default
g	Graph or DiGraph	The graph to plot.	required
x_label	array_like	See :func:plot_tmgraph.	None
nodemembers	sequence of sequence of int	See :func:plot_tmgraph.	None
nodesizerange	(float, float)	(min, max) node radius in pixels (pyvis/vis.js sizes nodes by radius, unlike matplotlib's area-based marker size -- so this default is not directly comparable to :func:plot_tmgraph's).	(6, 40)
nodesizemode	(log, rank, original)	See :func:plot_tmgraph.	'log'
colorlabel	str	Legend label (rendered as a small embedded colorbar image).	'x_label'
cmap	str or Colormap	Colormap for node coloring -- anything matplotlib accepts, including a discrete :class:~matplotlib.colors.ListedColormap for categorical labels.	'jet'
labelmethod	(mode, mean, median, none)	See :func:tmapper.labeling.find_node_label.	'mode'
nodeclim	(float, float)	Color axis limits. Defaults to (min(x_label), max(x_label)).	None
center_expand	float	See :func:plot_tmgraph.	4.0
title	str	Page heading (e.g. summarize the construction parameters here).	''
output_path	str	If given, also write the standalone HTML page to this path. Default None -- no file is written; use the returned html string instead (e.g. pass it directly to streamlit.components.v1.html(html, height=900) rather than writing to disk and reading it back).	None

Returns:

Name	Type	Description
net	Network	The populated network.
html	str	The standalone HTML page as a string (already includes the embedded colorbar legend) -- identical to what gets written to output_path when one is given.

Source code in src/tmapper/plotting.py

def plot_tmgraph_interactive(
    g,
    x_label=None,
    nodemembers=None,
    *,
    nodesizerange=(6, 40),
    nodesizemode="log",
    colorlabel="x_label",
    cmap="jet",
    labelmethod="mode",
    nodeclim=None,
    center_expand=4.0,
    title="",
    output_path=None,
):
    """Plot a temporal mapper graph as a draggable/zoomable/hoverable
    standalone HTML page (via pyvis/vis.js), instead of a static image.

    Uses the same layout, node-size, and coloring logic as
    :func:`plot_tmgraph` (so the two should look like the same network),
    with physics disabled -- nodes stay exactly where the layout put
    them, but can still be dragged, and hovering shows each node's size
    and color-label value.

    Parameters
    ----------
    g : networkx.Graph or networkx.DiGraph
        The graph to plot.
    x_label : array_like, optional
        See :func:`plot_tmgraph`.
    nodemembers : sequence of sequence of int, optional
        See :func:`plot_tmgraph`.
    nodesizerange : (float, float), default (6, 40)
        (min, max) node radius in pixels (pyvis/vis.js sizes nodes by
        radius, unlike matplotlib's area-based marker size -- so this
        default is not directly comparable to :func:`plot_tmgraph`'s).
    nodesizemode : {'log', 'rank', 'original'}, default 'log'
        See :func:`plot_tmgraph`.
    colorlabel : str, default 'x_label'
        Legend label (rendered as a small embedded colorbar image).
    cmap : str or matplotlib.colors.Colormap, default 'jet'
        Colormap for node coloring -- anything matplotlib accepts,
        including a discrete :class:`~matplotlib.colors.ListedColormap`
        for categorical labels.
    labelmethod : {'mode', 'mean', 'median', 'none'} or callable, default 'mode'
        See :func:`tmapper.labeling.find_node_label`.
    nodeclim : (float, float), optional
        Color axis limits. Defaults to (min(x_label), max(x_label)).
    center_expand : float, default 4.0
        See :func:`plot_tmgraph`.
    title : str, default ''
        Page heading (e.g. summarize the construction parameters here).
    output_path : str, optional
        If given, also write the standalone HTML page to this path.
        Default None -- no file is written; use the returned ``html``
        string instead (e.g. pass it directly to
        ``streamlit.components.v1.html(html, height=900)`` rather than
        writing to disk and reading it back).

    Returns
    -------
    net : pyvis.network.Network
        The populated network.
    html : str
        The standalone HTML page as a string (already includes the
        embedded colorbar legend) -- identical to what gets written to
        ``output_path`` when one is given.
    """
    import base64
    from io import BytesIO

    import matplotlib.pyplot as plt
    import matplotlib.colors as mcolors
    from pyvis.network import Network

    nodelist = list(g.nodes())
    n_nodes = len(nodelist)

    if nodemembers is None:
        nodemembers = [[i] for i in range(n_nodes)]

    if x_label is None:
        all_members = sorted({m for members in nodemembers for m in members})
        x_label = np.ones(len(all_members))
    x_label = np.asarray(x_label, dtype=float)

    if nodeclim is None:
        nodeclim = (float(np.min(x_label)), float(np.max(x_label)))

    nodesize = _prepare_node_size(nodemembers, nodesizerange, nodesizemode)
    nodelabel = find_node_label(nodemembers, x_label, labelmethod=labelmethod)
    pos = _compute_layout(g, nodelist, center_expand)

    # -- per-node hex colors via matplotlib's colormap machinery, so this
    # accepts the exact same `cmap`/`nodeclim` (continuous or discrete)
    # as plot_tmgraph.
    cmap_obj = plt.get_cmap(cmap) if isinstance(cmap, str) else cmap
    vmin, vmax = nodeclim if nodeclim[1] != nodeclim[0] else (nodeclim[0] - 0.5, nodeclim[0] + 0.5)
    norm = mcolors.Normalize(vmin=vmin, vmax=vmax)
    node_colors_hex = [mcolors.to_hex(cmap_obj(norm(v))) for v in nodelabel]

    net = Network(
        height="900px", width="100%", directed=g.is_directed(), bgcolor="#ffffff",
        cdn_resources="in_line", notebook=False, heading=title,
    )
    net.set_options('{"physics": {"enabled": false}}')

    SCALE = 800  # pyvis coordinates are pixels; our layout is normalized to roughly [-1.x, 1.x]
    order = np.argsort(nodesize, kind="stable")[::-1]  # largest first, smallest added last (drawn on top)
    for i in order:
        i = int(i)
        n = nodelist[i]
        x, y = pos[n]
        net.add_node(
            i,
            label="",
            title=f"node {n} | members={int(len(nodemembers[i]))} | {colorlabel}={nodelabel[i]:.3g}",
            x=float(x * SCALE), y=float(-y * SCALE),  # flip y: vis.js y-axis points down
            size=float(nodesize[i]), color=node_colors_hex[i], physics=False,
        )
    idx_of = {n: i for i, n in enumerate(nodelist)}
    edge_kwargs = {"arrows": "to"} if g.is_directed() else {}
    for u, v in g.edges():
        net.add_edge(idx_of[u], idx_of[v], color="#bbbbbb", width=1, **edge_kwargs)

    html = net.generate_html(notebook=False)

    # pyvis's own template renders "<h1>{{heading}}</h1>" twice
    # unconditionally -- drop the second copy.
    if title:
        heading_block = f"<h1>{title}</h1>"
        first = html.find(heading_block)
        second = html.find(heading_block, first + 1)
        if second != -1:
            html = html[:second] + html[second + len(heading_block):]

    # -- embed a small matplotlib colorbar as the legend, reusing the
    # exact cmap/norm above so it always matches the node colors.
    fig_cb, ax_cb = plt.subplots(figsize=(4, 0.5))
    cbar = fig_cb.colorbar(
        plt.cm.ScalarMappable(cmap=cmap_obj, norm=norm), cax=ax_cb, orientation="horizontal"
    )
    cbar.set_label(colorlabel)
    buf = BytesIO()
    fig_cb.savefig(buf, format="png", dpi=130, bbox_inches="tight", transparent=True)
    plt.close(fig_cb)
    legend_b64 = base64.b64encode(buf.getvalue()).decode("ascii")
    legend_html = (
        f'<div style="text-align:center;">'
        f'<img src="data:image/png;base64,{legend_b64}"></div>'
    )
    anchor = "</h1>" if title else "<body>"
    html = html.replace(anchor, anchor + "\n" + legend_html, 1)

    if output_path is not None:
        with open(output_path, "w", encoding="utf-8") as f:
            f.write(html)

    return net, html

Standalone graph builders

tmapper.knngraph

Port of tmapper_tools/knngraph.m from the MATLAB toolbox.

knngraph(X_or_D, k, *, reciprocal=True)

Construct an undirected graph based on k-nearest neighbors: each point is a node, linked to its k nearest neighbors.

Port of MATLAB's knngraph.m. A standalone graph builder with no temporal component (unlike :func:tknndigraph).

Parameters:

Name	Type	Description	Default
X_or_D	(array_like, shape(N, d) or (N, N))	Either raw coordinates (N points in d-dim space) or a precomputed (N, N) distance matrix. Auto-detected: treated as a distance matrix only if square and symmetric.	required
k	int	Number of nearest neighbors. Must be a positive integer strictly less than N.	required
reciprocal	bool	Whether to require a k-NN link to be mutual (both directions) to be kept.	True

Returns:

Type	Description
Graph	Unweighted undirected graph, one node per input point.

Source code in src/tmapper/knngraph.py

def knngraph(X_or_D, k, *, reciprocal=True):
    """Construct an undirected graph based on k-nearest neighbors: each
    point is a node, linked to its k nearest neighbors.

    Port of MATLAB's ``knngraph.m``. A standalone graph builder with no
    temporal component (unlike :func:`tknndigraph`).

    Parameters
    ----------
    X_or_D : array_like, shape (N, d) or (N, N)
        Either raw coordinates (N points in d-dim space) or a precomputed
        (N, N) distance matrix. Auto-detected: treated as a distance
        matrix only if square *and* symmetric.
    k : int
        Number of nearest neighbors. Must be a positive integer strictly
        less than N.
    reciprocal : bool, default True
        Whether to require a k-NN link to be mutual (both directions) to
        be kept.

    Returns
    -------
    networkx.Graph
        Unweighted undirected graph, one node per input point.
    """
    X_or_D = np.asarray(X_or_D, dtype=float)
    if X_or_D.ndim != 2:
        raise ValueError("X_or_D must be a 2D array.")
    if not (np.isscalar(k) and float(k) == round(float(k)) and k >= 1):
        raise ValueError("k must be a positive integer scalar.")
    k = int(k)

    nr, nc = X_or_D.shape
    if nr != nc or not np.allclose(X_or_D, X_or_D.T):
        D = cdist(X_or_D, X_or_D)
    else:
        D = X_or_D.copy()
    Nn = D.shape[0]

    if k >= Nn:
        raise ValueError(f"k must be smaller than the number of points ({Nn}).")

    order = np.argsort(D, axis=1, kind="stable")
    A = np.zeros((Nn, Nn), dtype=bool)
    rows = np.repeat(np.arange(Nn), k)
    cols = order[:, 1:1 + k].ravel()  # skip column 0: self (distance 0)
    A[rows, cols] = True

    if reciprocal:
        A = A & A.T
    else:
        A = A | A.T

    return nx.from_numpy_array(A)

tmapper.cknngraph

Port of tmapper_tools/cknngraph.m from the MATLAB toolbox.

cknngraph(X_or_D, k, delta, *, average=False)

Construct a graph based on continuous k-nearest neighbors (Berry & Sauer, 2019).

Port of MATLAB's cknngraph.m. Connects points x, y iff::

D(x, y) < delta * sqrt(D(x, x_k) * D(y, y_k))

where D(~, ~) is the distance and ~_k is the k-th nearest neighbor of ~.

Parameters:

Name	Type	Description	Default
X_or_D	(array_like, shape(N, d) or (N, N))	Either raw coordinates (N points in d-dim space) or a precomputed (N, N) distance matrix. Auto-detected: treated as a distance matrix only if square and symmetric.	required
k	int	Used to normalize distance by local point density. Must be a positive integer strictly less than N.	required
delta	float	Threshold for linking two nodes (must be a positive real number).	required
average	bool	Whether x_k is the average distance over the first k nearest neighbors (True) or just the k-th nearest neighbor (False).	False

Returns:

Type	Description
Graph	Unweighted undirected graph, one node per input point.

Source code in src/tmapper/cknngraph.py

def cknngraph(X_or_D, k, delta, *, average=False):
    """Construct a graph based on continuous k-nearest neighbors (Berry &
    Sauer, 2019).

    Port of MATLAB's ``cknngraph.m``. Connects points x, y iff::

        D(x, y) < delta * sqrt(D(x, x_k) * D(y, y_k))

    where ``D(~, ~)`` is the distance and ``~_k`` is the k-th nearest
    neighbor of ``~``.

    Parameters
    ----------
    X_or_D : array_like, shape (N, d) or (N, N)
        Either raw coordinates (N points in d-dim space) or a precomputed
        (N, N) distance matrix. Auto-detected: treated as a distance
        matrix only if square *and* symmetric.
    k : int
        Used to normalize distance by local point density. Must be a
        positive integer strictly less than N.
    delta : float
        Threshold for linking two nodes (must be a positive real
        number).
    average : bool, default False
        Whether ``x_k`` is the average distance over the first k nearest
        neighbors (True) or just the k-th nearest neighbor (False).

    Returns
    -------
    networkx.Graph
        Unweighted undirected graph, one node per input point.
    """
    X_or_D = np.asarray(X_or_D, dtype=float)
    if X_or_D.ndim != 2:
        raise ValueError("X_or_D must be a 2D array.")
    if not (np.isscalar(k) and float(k) == round(float(k)) and k >= 1):
        raise ValueError("k must be a positive integer scalar.")
    k = int(k)
    if not (np.isscalar(delta) and np.isreal(delta) and delta > 0):
        raise ValueError("delta must be a positive real scalar.")

    nr, nc = X_or_D.shape
    if nr != nc or not np.allclose(X_or_D, X_or_D.T):
        D = cdist(X_or_D, X_or_D)
    else:
        D = X_or_D.copy()
    Nn = D.shape[0]

    if k >= Nn:
        raise ValueError(f"k must be smaller than the number of points ({Nn}).")

    # -- find nearest neighbors and compute distance
    D_sorted = np.sort(D, axis=1)
    if average:
        Dk = D_sorted[:, 1:1 + k].mean(axis=1)  # col 0 is distance to self
    else:
        Dk = D_sorted[:, k]

    # -- normalize D by local density
    D_norm = D / np.sqrt(np.outer(Dk, Dk))

    # -- construct graph from adjacency matrix
    A = D_norm < delta
    np.fill_diagonal(A, False)

    return nx.from_numpy_array(A)

Graph & data utilities

tmapper.node_size(nodemembers)

Number of members in each node. Port of nodesize.m.

Source code in src/tmapper/graph_utils.py

def node_size(nodemembers):
    """Number of members in each node. Port of ``nodesize.m``."""
    return np.array([len(m) for m in nodemembers])

tmapper.node_measure(nodemembers)

Number of members in each node, normalized to sum to 1. Port of nodemeasure.m.

Source code in src/tmapper/graph_utils.py

def node_measure(nodemembers):
    """Number of members in each node, normalized to sum to 1. Port of
    ``nodemeasure.m``."""
    n = node_size(nodemembers).astype(float)
    return n / n.sum()

tmapper.normalize_geodesic(geod, nsize=None, *, exclude_diag=False)

Normalize a geodesic distance matrix by node measure.

Port of normgeo.m. Note: unlike the MATLAB original, an omitted nsize defaults to a uniform vector of ones (equal-size nodes), which is clearly the intended behavior from the docstring -- the MATLAB version's ones(N_nodes) (a matrix, not ones(N_nodes,1)) appears to be an unexercised bug in that rarely-hit default-value path.

Parameters:

Name	Type	Description	Default
geod	(array_like, shape(N, N))	Geodesic distance matrix.	required
nsize	(array_like, shape(N))	Size of each node. Defaults to a uniform vector of ones.	None
exclude_diag	bool	Whether to exclude the diagonal when computing the normalizing factor.	False

Returns:

Name	Type	Description
geod_n	(ndarray, shape(N, N))	Normalized geodesic distance matrix.
nm	(ndarray, shape(N))	Node measure (probability), sums to 1.

Source code in src/tmapper/graph_utils.py

def normalize_geodesic(geod, nsize=None, *, exclude_diag=False):
    """Normalize a geodesic distance matrix by node measure.

    Port of ``normgeo.m``. Note: unlike the MATLAB original, an omitted
    ``nsize`` defaults to a uniform *vector* of ones (equal-size nodes),
    which is clearly the intended behavior from the docstring -- the
    MATLAB version's ``ones(N_nodes)`` (a matrix, not ``ones(N_nodes,1)``)
    appears to be an unexercised bug in that rarely-hit default-value path.

    Parameters
    ----------
    geod : array_like, shape (N, N)
        Geodesic distance matrix.
    nsize : array_like, shape (N,), optional
        Size of each node. Defaults to a uniform vector of ones.
    exclude_diag : bool, default False
        Whether to exclude the diagonal when computing the normalizing
        factor.

    Returns
    -------
    geod_n : numpy.ndarray, shape (N, N)
        Normalized geodesic distance matrix.
    nm : numpy.ndarray, shape (N,)
        Node measure (probability), sums to 1.
    """
    geod = np.array(geod, dtype=float, copy=True)
    n_nodes = geod.shape[0]

    if nsize is None:
        nsize = np.ones(n_nodes)
    nsize = np.asarray(nsize, dtype=float)

    # -- handle inf
    if np.any(np.isinf(geod)):
        finite = geod[np.isfinite(geod)]
        geod[np.isinf(geod)] = finite.max() if finite.size else 0.0

    # -- weight nodes and geodesics
    nm = nsize / nsize.sum()
    geod_n = np.outer(nm, nm) * geod ** 2

    # -- normalizing factor: the sum of weighted geodesics
    if exclude_diag:
        mask = ~np.eye(n_nodes, dtype=bool)
        normfactor = np.sqrt(geod_n[mask].sum())
    else:
        normfactor = np.sqrt(geod_n.sum())

    # -- normalize
    if normfactor != 0:
        geod_n = geod / normfactor

    return geod_n, nm

tmapper.normalize_tcm(tcm, *, normtype='max', infreplace='max')

Normalize a temporal connectivity matrix by its maximal finite value (or its norm). Port of normtcm.m.

Parameters:

Name	Type	Description	Default
tcm	array_like	Temporal connectivity matrix.	required
normtype	(max, norm)	How to compute the normalizing factor: the matrix's maximum finite value, or its (Frobenius/vector) norm.	'max'
infreplace	(max, nan)	How to handle inf entries before normalizing: replace with the matrix's maximum finite value, or with NaN.	'max'

Returns:

Type	Description
ndarray

Source code in src/tmapper/graph_utils.py

def normalize_tcm(tcm, *, normtype="max", infreplace="max"):
    """Normalize a temporal connectivity matrix by its maximal finite
    value (or its norm). Port of ``normtcm.m``.

    Parameters
    ----------
    tcm : array_like
        Temporal connectivity matrix.
    normtype : {'max', 'norm'}, default 'max'
        How to compute the normalizing factor: the matrix's maximum
        finite value, or its (Frobenius/vector) norm.
    infreplace : {'max', 'nan'}, default 'max'
        How to handle inf entries before normalizing: replace with the
        matrix's maximum finite value, or with NaN.

    Returns
    -------
    numpy.ndarray
    """
    tcm = np.array(tcm, dtype=float, copy=True)

    if infreplace == "max":
        finite = tcm[np.isfinite(tcm)]
        if finite.size:
            tcm[np.isinf(tcm)] = finite.max()
    elif infreplace == "nan":
        tcm[np.isinf(tcm)] = np.nan
    else:
        raise ValueError(f"Unknown infreplace: {infreplace!r}")

    if normtype == "max":
        normfactor = np.nanmax(tcm)
    elif normtype == "norm":
        normfactor = np.linalg.norm(tcm.ravel())
    else:
        raise ValueError(f"Unknown normtype: {normtype!r}")

    if normfactor != 0:
        return tcm / normfactor
    return tcm

tmapper.members_to_tidx(members, tidx)

Translate positional-index members to real tidx values.

Port of members2tidx.m. Use this when tidx is non-contiguous or offset (e.g. real timestamps with gaps): :func:filtergraph's members output gives positional indices into your original data, not tidx values.

Parameters:

Name	Type	Description	Default
members	sequence of sequence of int	Positional indices into the original data, e.g. from :func:filtergraph.	required
tidx	array_like	The same tidx array originally passed to :func:tknndigraph.	required

Returns:

Type	Description
list of numpy.ndarray

Source code in src/tmapper/graph_utils.py

def members_to_tidx(members, tidx):
    """Translate positional-index members to real ``tidx`` values.

    Port of ``members2tidx.m``. Use this when ``tidx`` is non-contiguous
    or offset (e.g. real timestamps with gaps): :func:`filtergraph`'s
    ``members`` output gives positional indices into your original data,
    not ``tidx`` values.

    Parameters
    ----------
    members : sequence of sequence of int
        Positional indices into the original data, e.g. from
        :func:`filtergraph`.
    tidx : array_like
        The same ``tidx`` array originally passed to :func:`tknndigraph`.

    Returns
    -------
    list of numpy.ndarray
    """
    tidx = np.asarray(tidx)
    return [tidx[np.asarray(list(m), dtype=int)] for m in members]

tmapper.subgraph_from_members(g_simp, members, include_members)

Extract the subgraph of a transition network restricted to a subset of original time points. Port of subgraphFromMembers.m.

Parameters:

Name	Type	Description	Default
g_simp	Graph or DiGraph	The simplified transition network.	required
members	sequence of sequence of int	members[n] gives the original time-point indices belonging to node n (in the same node order as g_simp.nodes()).	required
include_members	iterable of int	The subset of original time points to keep.	required

Returns:

Name	Type	Description
g_sub	Graph or DiGraph	The induced subgraph containing only nodes with at least one member in include_members.
members_sub	list of list	Members of each node of g_sub, restricted to include_members.
sub_orig_nodeidx	list of int	Positional indices (into the original g_simp.nodes()) of the nodes kept in g_sub.

Source code in src/tmapper/graph_utils.py

def subgraph_from_members(g_simp, members, include_members):
    """Extract the subgraph of a transition network restricted to a
    subset of original time points. Port of ``subgraphFromMembers.m``.

    Parameters
    ----------
    g_simp : networkx.Graph or networkx.DiGraph
        The simplified transition network.
    members : sequence of sequence of int
        ``members[n]`` gives the original time-point indices belonging
        to node n (in the same node order as ``g_simp.nodes()``).
    include_members : iterable of int
        The subset of original time points to keep.

    Returns
    -------
    g_sub : networkx.Graph or networkx.DiGraph
        The induced subgraph containing only nodes with at least one
        member in ``include_members``.
    members_sub : list of list
        Members of each node of ``g_sub``, restricted to
        ``include_members``.
    sub_orig_nodeidx : list of int
        Positional indices (into the original ``g_simp.nodes()``) of the
        nodes kept in ``g_sub``.
    """
    include_set = set(include_members)
    members_sub_all = [sorted(set(m) & include_set) for m in members]
    sub_orig_nodeidx = [i for i, m in enumerate(members_sub_all) if len(m) > 0]

    nodelist = list(g_simp.nodes())
    sub_nodes = [nodelist[i] for i in sub_orig_nodeidx]
    g_sub = g_simp.subgraph(sub_nodes).copy()
    members_sub = [members_sub_all[i] for i in sub_orig_nodeidx]

    return g_sub, members_sub, sub_orig_nodeidx

tmapper.sym_dyn_to_digraph(sym_dyn)

Construct a directed graph from a single time series of symbolic dynamics. Port of symDyn2digraph.m.

Parameters:

Name	Type	Description	Default
sym_dyn	array_like of int	A vector of N labels (purely nominal -- numeric differences between labels are irrelevant), one per time point.	required

Returns:

Name	Type	Description
dg	DiGraph	One node per unique state in sym_dyn; an edge wherever a transition between two states is observed.
dwelltime	ndarray	Number of time points belonging to each state (in the same order as dg.nodes()).
nodemembers	list of numpy.ndarray	nodemembers[n] gives the time-point indices where the system was in dg's n-th state.

Source code in src/tmapper/graph_utils.py

def sym_dyn_to_digraph(sym_dyn):
    """Construct a directed graph from a single time series of symbolic
    dynamics. Port of ``symDyn2digraph.m``.

    Parameters
    ----------
    sym_dyn : array_like of int
        A vector of N labels (purely nominal -- numeric differences
        between labels are irrelevant), one per time point.

    Returns
    -------
    dg : networkx.DiGraph
        One node per unique state in ``sym_dyn``; an edge wherever a
        transition between two states is observed.
    dwelltime : numpy.ndarray
        Number of time points belonging to each state (in the same
        order as ``dg.nodes()``).
    nodemembers : list of numpy.ndarray
        ``nodemembers[n]`` gives the time-point indices where the
        system was in ``dg``'s n-th state.
    """
    sym_dyn = np.asarray(sym_dyn)
    tidx = np.arange(len(sym_dyn))

    state_names, state_idx = np.unique(sym_dyn, return_inverse=True)
    dwelltime = np.bincount(state_idx, minlength=len(state_names))

    dg = nx.DiGraph()
    dg.add_nodes_from(state_names.tolist())

    if len(sym_dyn) > 1:
        changed = np.diff(sym_dyn) != 0
        tidx_trans = tidx[:-1][changed]
        state_s = sym_dyn[tidx_trans]
        state_t = sym_dyn[tidx_trans + 1]
        trans = np.unique(np.stack([state_s, state_t], axis=1), axis=0) if tidx_trans.size else np.empty((0, 2))
        dg.add_edges_from((s, t) for s, t in trans.tolist())

    nodemembers = [tidx[state_idx == n] for n in range(len(state_names))]

    return dg, dwelltime, nodemembers

tmapper.digraph_to_graph(dg)

Convert a directed graph to an undirected graph, averaging the weights of the two edges between each pair of nodes. Port of digraph2graph.m.

Parameters:

Name	Type	Description	Default
dg	DiGraph	The directed graph to symmetrize.	required

Returns:

Type	Description
Graph

Source code in src/tmapper/graph_utils.py

def digraph_to_graph(dg):
    """Convert a directed graph to an undirected graph, averaging the
    weights of the two edges between each pair of nodes. Port of
    ``digraph2graph.m``.

    Parameters
    ----------
    dg : networkx.DiGraph
        The directed graph to symmetrize.

    Returns
    -------
    networkx.Graph
    """
    nodelist = list(dg.nodes())
    A = nx.to_numpy_array(dg, nodelist=nodelist, weight="weight")
    A_sym = (A + A.T) / 2
    g = nx.from_numpy_array(A_sym)
    return nx.relabel_nodes(g, {i: nodelist[i] for i in range(len(nodelist))})

tmapper.find_blocks(indicator)

Find the start, end, and size of contiguous True/1 runs ("blocks") in a 0/1 indicator array. Port of findtaskn.m.

Parameters:

Name	Type	Description	Default
indicator	array_like	A 1D array of 0s/1s (or booleans).	required

Returns:

Name	Type	Description
block_start	ndarray	0-indexed position of the first sample of each block.
block_end	ndarray	0-indexed position of the last sample of each block.
block_size	ndarray	Number of samples in each block.

Source code in src/tmapper/graph_utils.py

def find_blocks(indicator):
    """Find the start, end, and size of contiguous True/1 runs
    ("blocks") in a 0/1 indicator array. Port of ``findtaskn.m``.

    Parameters
    ----------
    indicator : array_like
        A 1D array of 0s/1s (or booleans).

    Returns
    -------
    block_start : numpy.ndarray
        0-indexed position of the first sample of each block.
    block_end : numpy.ndarray
        0-indexed position of the last sample of each block.
    block_size : numpy.ndarray
        Number of samples in each block.
    """
    ind = np.asarray(indicator).astype(int).ravel()

    padded_start = np.concatenate([[0], ind])
    block_start = np.flatnonzero(np.diff(padded_start) == 1)

    padded_end = np.concatenate([ind, [0]])
    block_end = np.flatnonzero(np.diff(padded_end) == -1)

    if len(block_start) != len(block_end):
        raise ValueError("Block structure incomplete.")

    block_size = block_end - block_start + 1
    return block_start, block_end, block_size

Cycle / path analysis toolkit

tmapper.cycle_count

Port of tmapper_tools/CycleCount.m from the MATLAB toolbox.

CycleCount.m is itself a third-party algorithm (not original to the tmapper toolbox), reproduced here under its original license:

Authors: P.-L. Giscard, N. Kriege, R. Wilson, Septembre 2016

Copyright (c) 2017, Pierre-Louis Giscard
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

    * Redistributions of source code must retain the above
      copyright notice, this list of conditions and the following
      disclaimer.
    * Redistributions in binary form must reproduce the above
      copyright notice, this list of conditions and the following
      disclaimer in the documentation and/or other materials
      provided with the distribution.
    * Neither the name of the University of York nor the names of
      its contributors may be used to endorse or promote products
      derived from this software without specific prior written
      permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.

cycle_count(A, L0)

Count all simple cycles of length up to L0 (inclusive) on a graph whose adjacency matrix is A, via the combinatorial-sieve algorithm of Giscard, Kriege & Wilson (2016).

Parameters:

Name	Type	Description	Default
A	(array_like, shape(N, N))	Adjacency matrix of the graph (directed or undirected, weighted or unweighted).	required
L0	int	Maximum length of the simple cycles to count.	required

Returns:

Type	Description
(ndarray, shape(L0))	primes[i] is the number of simple cycles of length i + 1, for i in 0..L0-1.

Source code in src/tmapper/cycle_count.py

def cycle_count(A, L0):
    """Count all simple cycles of length up to ``L0`` (inclusive) on a
    graph whose adjacency matrix is ``A``, via the combinatorial-sieve
    algorithm of Giscard, Kriege & Wilson (2016).

    Parameters
    ----------
    A : array_like, shape (N, N)
        Adjacency matrix of the graph (directed or undirected, weighted
        or unweighted).
    L0 : int
        Maximum length of the simple cycles to count.

    Returns
    -------
    numpy.ndarray, shape (L0,)
        ``primes[i]`` is the number of simple cycles of length ``i + 1``,
        for ``i`` in ``0..L0-1``.
    """
    A = np.asarray(A, dtype=float).copy()
    L0 = int(L0)

    primes = np.zeros(L0, dtype=complex)
    primes[0] += np.trace(A)  # self-loops (length-1 cycles)
    np.fill_diagonal(A, 0)
    A = _clean_matrix(A)

    directed = not np.array_equal(A, A.T)
    if directed:
        A_undir = ((A != 0) | (A.T != 0)).astype(float)
    else:
        A_undir = (A != 0).astype(float)

    size = A.shape[0]
    if L0 > size:
        L0 = size

    allowed_vert = np.ones(size, dtype=bool)
    for i in range(size - 1):
        allowed_vert[i] = False
        neighbourhood = np.zeros(size)
        neighbourhood[i] = 1
        neighbourhood = neighbourhood + A_undir[i, :]
        primes = _recursive_subgraphs(A, A_undir, L0, [i], allowed_vert, primes, neighbourhood)

    return primes.real

tmapper.cycle_count2p

Port of tmapper_tools/CycleCount2p.m and reorgCycles.m from the MATLAB toolbox.

cycle_count2p(A, *, simple=True)

Estimate the number of cycles of different lengths by finding the smallest cycle passing through every pair of vertices.

Port of MATLAB's CycleCount2p.m. Only counts unique cycles, and by default only unique simple cycles (no repeated nodes other than the start/end).

Parameters:

Name	Type	Description	Default
A	(array_like, shape(N, N))	Adjacency matrix for a simple directed graph.	required
simple	bool	Whether to only count simple cycles.	True

Returns:

Name	Type	Description
cyc_count	ndarray	Number of cycles per unique length (shortest to longest).
cyc_len	ndarray	The unique cycle lengths themselves (shortest to longest).
cyc_path	list of numpy.ndarray	cyc_path[m] is an (Nm, Lm) array of the Nm cycles of length cyc_len[m].
all_cycles	list of list	Every cycle's path, flattened into a single list.

Source code in src/tmapper/cycle_count2p.py

def cycle_count2p(A, *, simple=True):
    """Estimate the number of cycles of different lengths by finding the
    smallest cycle passing through every pair of vertices.

    Port of MATLAB's ``CycleCount2p.m``. Only counts unique cycles, and
    by default only unique simple cycles (no repeated nodes other than
    the start/end).

    Parameters
    ----------
    A : array_like, shape (N, N)
        Adjacency matrix for a simple directed graph.
    simple : bool, default True
        Whether to only count simple cycles.

    Returns
    -------
    cyc_count : numpy.ndarray
        Number of cycles per unique length (shortest to longest).
    cyc_len : numpy.ndarray
        The unique cycle lengths themselves (shortest to longest).
    cyc_path : list of numpy.ndarray
        ``cyc_path[m]`` is an (Nm, Lm) array of the Nm cycles of length
        ``cyc_len[m]``.
    all_cycles : list of list
        Every cycle's path, flattened into a single list.
    """
    g = nx.from_numpy_array(np.asarray(A), create_using=nx.DiGraph)
    n = g.number_of_nodes()

    all_closed = []
    for s in range(n):
        for t in range(s + 1, n):
            if not (nx.has_path(g, s, t) and nx.has_path(g, t, s)):
                continue
            spf = nx.shortest_path(g, s, t)  # forward path
            spb = nx.shortest_path(g, t, s)  # backward path
            cycle = spf + spb[1:-1]
            if not simple or len(set(cycle)) == len(cycle):
                root_idx = cycle.index(min(cycle))
                all_closed.append(cycle[root_idx:] + cycle[:root_idx])  # start at smallest node

    cyc_length = np.array([len(c) for c in all_closed], dtype=int)
    cyc_len = np.unique(cyc_length) if cyc_length.size else np.array([], dtype=int)

    cyc_path = []
    cyc_count = []
    for L in cyc_len:
        rows = sorted({tuple(c) for c, l in zip(all_closed, cyc_length) if l == L})
        cyc_path.append(np.array(rows, dtype=int))
        cyc_count.append(len(rows))
    cyc_count = np.array(cyc_count, dtype=int)

    all_cycles = [list(row) for block in cyc_path for row in block]

    return cyc_count, cyc_len, cyc_path, all_cycles

reorg_cycles(cyc_path)

Unpack the cyc_path output of :func:cycle_count2p into a flat list where each element is a single cycle's path.

Port of MATLAB's reorgCycles.m.

Source code in src/tmapper/cycle_count2p.py

def reorg_cycles(cyc_path):
    """Unpack the ``cyc_path`` output of :func:`cycle_count2p` into a flat
    list where each element is a single cycle's path.

    Port of MATLAB's ``reorgCycles.m``.
    """
    return [list(row) for block in cyc_path for row in np.asarray(block)]

tmapper.reorg_cycles(cyc_path)

Unpack the cyc_path output of :func:cycle_count2p into a flat list where each element is a single cycle's path.

Port of MATLAB's reorgCycles.m.

Source code in src/tmapper/cycle_count2p.py

def reorg_cycles(cyc_path):
    """Unpack the ``cyc_path`` output of :func:`cycle_count2p` into a flat
    list where each element is a single cycle's path.

    Port of MATLAB's ``reorgCycles.m``.
    """
    return [list(row) for block in cyc_path for row in np.asarray(block)]

tmapper.cycle_path_overlap(c, *, cycle=True, overlap_type='edge', grpvar=None)

Calculate the overlap between cycles or paths.

Port of MATLAB's CyclePathOverlap.m.

Parameters:

Name	Type	Description	Default
c	sequence of sequence of int	c[n] is a path of node indices.	required
cycle	bool	Whether the paths in c are cycles (adds a wraparound edge from the last node back to the first when overlap_type is 'edge').	True
overlap_type	(edge, node)	Whether to calculate overlap between edges or between nodes.	'edge'
grpvar	array_like of int	A grouping variable for the paths. By default, each path is its own group.	None

Returns:

Name	Type	Description
CO	(ndarray, shape(P, P))	Percent overlap between every pair of groups (1 on the diagonal).
grpnames	ndarray	The unique group labels, in the same order as CO's rows.

Source code in src/tmapper/cycle_overlap.py

def cycle_path_overlap(c, *, cycle=True, overlap_type="edge", grpvar=None):
    """Calculate the overlap between cycles or paths.

    Port of MATLAB's ``CyclePathOverlap.m``.

    Parameters
    ----------
    c : sequence of sequence of int
        ``c[n]`` is a path of node indices.
    cycle : bool, default True
        Whether the paths in ``c`` are cycles (adds a wraparound edge
        from the last node back to the first when ``overlap_type``
        is 'edge').
    overlap_type : {'edge', 'node'}, default 'edge'
        Whether to calculate overlap between edges or between nodes.
    grpvar : array_like of int, optional
        A grouping variable for the paths. By default, each path is its
        own group.

    Returns
    -------
    CO : numpy.ndarray, shape (P, P)
        Percent overlap between every pair of groups (1 on the
        diagonal).
    grpnames : numpy.ndarray
        The unique group labels, in the same order as ``CO``'s rows.
    """
    if overlap_type == "edge":
        def as_items(x):
            x = list(x)
            edges = list(zip(x, x[1:] + [x[0]]))
            if not cycle:
                edges = edges[:-1]  # drop the wraparound edge
            return edges
    elif overlap_type == "node":
        def as_items(x):
            return list(x)
    else:
        raise ValueError(f"Unknown type: {overlap_type!r}")

    c_repr = [as_items(x) for x in c]

    Nc = len(c_repr)
    if grpvar is None:
        grpvar = np.arange(Nc)
    grpvar = np.asarray(grpvar)
    grpnames = np.unique(grpvar)
    Ngrp = len(grpnames)

    groups = []
    for name in grpnames:
        items = set()
        for idx in np.flatnonzero(grpvar == name):
            items |= set(c_repr[idx])
        groups.append(items)

    CO = np.full((Ngrp, Ngrp), np.nan)
    for ii in range(Ngrp):
        for jj in range(ii):
            inter = len(groups[ii] & groups[jj])
            union = len(groups[ii] | groups[jj])
            CO[ii, jj] = inter / union
            CO[jj, ii] = CO[ii, jj]
    np.fill_diagonal(CO, 1)

    return CO, grpnames

tmapper.cycle_cluster

Port of tmapper_tools/CycleCluster.m from the MATLAB toolbox.

cycle_cluster(allcycles, thres, *, plotmat=True, plotmds=False, plothist=True, reordermat=True, ax=None, return_ax=False)

Cluster cycles based on the fraction of overlap between them (shared nodes / union of nodes).

Port of MATLAB's CycleCluster.m.

Parameters:

Name	Type	Description	Default
allcycles	sequence of sequence of int	allcycles[n] is the path of one cycle.	required
thres	float	Cut-off threshold (0-1) for single-linkage clustering: if overlap > thres, two cycles belong to the same cluster.	required
plotmat	bool	Whether to plot the overlap matrix.	True
plotmds	bool	Whether to plot a 2D classical-MDS projection of the cycles.	False
plothist	bool	Whether to plot the histogram of linkage distances.	True
reordermat	bool	Whether to reorder the overlap matrix by cluster assignment.	True
ax	Axes	Axes for the overlap-matrix plot (only used if plotmat). A new figure is created if not given.	None
return_ax	bool	If True, also return the axes used for the overlap-matrix plot (None if plotmat is False), so a caller can add further annotations to it (e.g. :func:cycle_path_decomp's cluster-block outlines).	False

Returns:

Name	Type	Description
	ndarray	1-indexed cluster label for each cycle (matching MATLAB's cluster() convention).
ax	Axes or None	Only returned if return_ax is True.

Source code in src/tmapper/cycle_cluster.py

def cycle_cluster(allcycles, thres, *, plotmat=True, plotmds=False, plothist=True,
                   reordermat=True, ax=None, return_ax=False):
    """Cluster cycles based on the fraction of overlap between them
    (shared nodes / union of nodes).

    Port of MATLAB's ``CycleCluster.m``.

    Parameters
    ----------
    allcycles : sequence of sequence of int
        ``allcycles[n]`` is the path of one cycle.
    thres : float
        Cut-off threshold (0-1) for single-linkage clustering: if
        overlap > thres, two cycles belong to the same cluster.
    plotmat : bool, default True
        Whether to plot the overlap matrix.
    plotmds : bool, default False
        Whether to plot a 2D classical-MDS projection of the cycles.
    plothist : bool, default True
        Whether to plot the histogram of linkage distances.
    reordermat : bool, default True
        Whether to reorder the overlap matrix by cluster assignment.
    ax : matplotlib.axes.Axes, optional
        Axes for the overlap-matrix plot (only used if ``plotmat``).
        A new figure is created if not given.
    return_ax : bool, default False
        If True, also return the axes used for the overlap-matrix plot
        (None if ``plotmat`` is False), so a caller can add further
        annotations to it (e.g. :func:`cycle_path_decomp`'s cluster-block
        outlines).

    Returns
    -------
    numpy.ndarray
        1-indexed cluster label for each cycle (matching MATLAB's
        ``cluster()`` convention).
    ax : matplotlib.axes.Axes or None
        Only returned if ``return_ax`` is True.
    """
    Nc = len(allcycles)

    if Nc <= 1:
        cluster_idx = np.ones(Nc, dtype=int)
        return (cluster_idx, None) if return_ax else cluster_idx

    prct_overlap = np.zeros((Nc, Nc))
    sets = [set(c) for c in allcycles]
    for ii in range(Nc):
        for jj in range(Nc):
            prct_overlap[ii, jj] = len(sets[ii] & sets[jj]) / len(sets[ii] | sets[jj])

    if plotmds:
        import matplotlib.pyplot as plt
        Y = _classical_mds(1 - prct_overlap, 2)
        plt.figure()
        plt.scatter(Y[:, 0], Y[:, 1])

    Z = linkage(squareform(1 - prct_overlap, checks=False), method="complete")

    if plothist:
        import matplotlib.pyplot as plt
        plt.figure()
        plt.hist(1 - Z[:, 2], bins=np.linspace(0, 1, 21))
        ylim = plt.ylim()
        plt.plot([thres, thres], ylim)
        plt.xlabel("overlap")
        plt.ylabel("count")

    cluster_idx = fcluster(Z, t=1 - thres, criterion="distance")

    if plotmat:
        import matplotlib.pyplot as plt
        if ax is None:
            _, ax = plt.subplots()
        if reordermat:
            order = np.argsort(cluster_idx, kind="stable")
            im = ax.imshow(prct_overlap[np.ix_(order, order)], origin="lower", vmin=0, vmax=1)
            ax.set_xlabel("cycle/path index (reordered)")
            ax.set_ylabel("cycle/path index (reordered)")
        else:
            im = ax.imshow(prct_overlap, origin="lower", vmin=0, vmax=1)
            ax.set_xlabel("cycle/path index")
            ax.set_ylabel("cycle/path index")
        cbar = plt.colorbar(im, ax=ax)
        cbar.set_label("fraction of overlap")
        ax.set_aspect("equal")
    else:
        ax = None

    return (cluster_idx, ax) if return_ax else cluster_idx

tmapper.cycle_cluster_conn

Port of tmapper_tools/CycleClusterConn.m from the MATLAB toolbox.

cycle_cluster_conn(dg, allcycles, cluster_idx)

Connectivity between M clusters of N cycles.

Port of MATLAB's CycleClusterConn.m.

Parameters:

Name	Type	Description	Default
dg	DiGraph	The directed graph the cycles/paths live on.	required
allcycles	sequence of sequence of int	allcycles[n] is the path of one cycle.	required
cluster_idx	array_like of int	Cluster label of each cycle (e.g. from :func:cycle_cluster).	required

Returns:

Name	Type	Description
cluster_conn	list of list of set	cluster_conn[i][j] (== cluster_conn[j][i]) contains nodes on the boundary between loop-clusters i and j.
cluster_conn_dir	list of list of set	cluster_conn_dir[i][j] contains nodes in the boundary of cluster j that receive a link from at least one node of cluster i (excluding shared nodes).
clusters_nodes	list of list	All nodes in each cluster (first-seen order).
clusters_boundary	list of set	Boundary nodes of each cluster.
clusters_interior	list of list	Nodes of each cluster that are not on its boundary.
clusters_crtpts	list of list	Critical points (nodes with more than one source or target) in each cluster.
clusters_intcrtpts	list of list	Critical points in each cluster that are not on its boundary.

Source code in src/tmapper/cycle_cluster_conn.py

def cycle_cluster_conn(dg, allcycles, cluster_idx):
    """Connectivity between M clusters of N cycles.

    Port of MATLAB's ``CycleClusterConn.m``.

    Parameters
    ----------
    dg : networkx.DiGraph
        The directed graph the cycles/paths live on.
    allcycles : sequence of sequence of int
        ``allcycles[n]`` is the path of one cycle.
    cluster_idx : array_like of int
        Cluster label of each cycle (e.g. from :func:`cycle_cluster`).

    Returns
    -------
    cluster_conn : list of list of set
        ``cluster_conn[i][j]`` (== ``cluster_conn[j][i]``) contains
        nodes on the boundary between loop-clusters i and j.
    cluster_conn_dir : list of list of set
        ``cluster_conn_dir[i][j]`` contains nodes in the boundary of
        cluster j that receive a link from at least one node of cluster
        i (excluding shared nodes).
    clusters_nodes : list of list
        All nodes in each cluster (first-seen order).
    clusters_boundary : list of set
        Boundary nodes of each cluster.
    clusters_interior : list of list
        Nodes of each cluster that are not on its boundary.
    clusters_crtpts : list of list
        Critical points (nodes with more than one source or target) in
        each cluster.
    clusters_intcrtpts : list of list
        Critical points in each cluster that are not on its boundary.
    """
    cluster_idx = np.asarray(cluster_idx)
    clusters = np.unique(cluster_idx)
    n_clusters = len(clusters)

    out_degree = dict(dg.out_degree())
    in_degree = dict(dg.in_degree())
    crtpts = {n for n in dg.nodes() if out_degree[n] > 1 or in_degree[n] > 1}

    clusters_nodes = [None] * n_clusters
    clusters_crtpts = [None] * n_clusters
    clusters_intcrtpts = [None] * n_clusters
    clusters_inneighbors = [None] * n_clusters
    clusters_outneighbors = [None] * n_clusters
    clusters_boundary = [None] * n_clusters
    clusters_interior = [None] * n_clusters

    for ii, c in enumerate(clusters):
        cycles_ii = [allcycles[k] for k in range(len(allcycles)) if cluster_idx[k] == c]

        seen, seen_set = [], set()
        for cyc in cycles_ii:
            for node in cyc:
                if node not in seen_set:
                    seen_set.add(node)
                    seen.append(node)
        clusters_nodes[ii] = seen

        crtpts_ii = [n for n in seen if n in crtpts]
        clusters_crtpts[ii] = crtpts_ii

        innbg, outnbg = set(), set()
        for cp in crtpts_ii:
            innbg |= set(dg.predecessors(cp))
            outnbg |= set(dg.successors(cp))
        innbg -= seen_set
        outnbg -= seen_set

        bd = set()
        for n in innbg:
            bd |= set(dg.successors(n))
        for n in outnbg:
            bd |= set(dg.predecessors(n))
        bd &= seen_set

        clusters_inneighbors[ii] = innbg
        clusters_outneighbors[ii] = outnbg
        clusters_boundary[ii] = bd
        clusters_interior[ii] = [n for n in seen if n not in bd]
        clusters_intcrtpts[ii] = [n for n in crtpts_ii if n not in bd]

    cluster_conn = [[None] * n_clusters for _ in range(n_clusters)]
    cluster_conn_dir = [[None] * n_clusters for _ in range(n_clusters)]
    cluster_conn_in = [[None] * n_clusters for _ in range(n_clusters)]
    cluster_conn_out = [[None] * n_clusters for _ in range(n_clusters)]

    for ii in range(n_clusters):
        for jj in range(n_clusters):
            cluster_conn[ii][jj] = clusters_boundary[ii] & set(clusters_nodes[jj])
            cluster_conn_in[ii][jj] = clusters_inneighbors[jj] & set(clusters_nodes[ii])
            cluster_conn_out[ii][jj] = clusters_outneighbors[ii] & set(clusters_nodes[jj])

    for ii in range(n_clusters):
        for jj in range(n_clusters):
            bd = set()
            for nn in cluster_conn[ii][jj]:
                preds = set(dg.predecessors(nn))
                if preds & cluster_conn_in[ii][jj]:
                    bd.add(nn)
            cluster_conn_dir[ii][jj] = bd & cluster_conn[ii][jj]

    return (cluster_conn, cluster_conn_dir, clusters_nodes, clusters_boundary,
            clusters_interior, clusters_crtpts, clusters_intcrtpts)

tmapper.cycle_cutter

Port of tmapper_tools/CycleCutter.m from the MATLAB toolbox.

cycle_cutter(cyc, node_name)

Cut a cycle into multiple paths at given nodes.

Port of MATLAB's CycleCutter.m.

Parameters:

Name	Type	Description	Default
cyc	sequence	A single cycle's path (node names, usually integers).	required
node_name	sequence or scalar	Names of nodes that are the cutting points of the cycle. Names not contained in the cycle are ignored.	required

Returns:

Type	Description
list of list	Each element is a path from one cutting point to the next. Wraps around: if there are M cutting points, this list has M elements, and the last one wraps from the last cutting point, through the end of cyc, back through the start, up to (and including) the first cutting point.

Source code in src/tmapper/cycle_cutter.py

def cycle_cutter(cyc, node_name):
    """Cut a cycle into multiple paths at given nodes.

    Port of MATLAB's ``CycleCutter.m``.

    Parameters
    ----------
    cyc : sequence
        A single cycle's path (node names, usually integers).
    node_name : sequence or scalar
        Names of nodes that are the cutting points of the cycle. Names
        not contained in the cycle are ignored.

    Returns
    -------
    list of list
        Each element is a path from one cutting point to the next.
        Wraps around: if there are ``M`` cutting points, this list has
        ``M`` elements, and the last one wraps from the last cutting
        point, through the end of ``cyc``, back through the start, up to
        (and including) the first cutting point.
    """
    cyc = list(cyc)

    if np.isscalar(node_name):
        node_set = {node_name}
    else:
        node_set = set(node_name)

    node_idx = [i for i, v in enumerate(cyc) if v in node_set]
    if not node_idx:
        return [cyc]

    n_path = len(node_idx)
    nodepath = [None] * n_path

    for n in range(n_path - 1):
        nodepath[n] = cyc[node_idx[n]: node_idx[n + 1] + 1]

    nodepath[n_path - 1] = cyc[node_idx[-1]:] + cyc[:node_idx[0] + 1]

    return nodepath

tmapper.cycles_to_paths

Port of tmapper_tools/Cycles2Paths.m from the MATLAB toolbox.

cycles_to_paths(allcycles, cutpts)

Cut a set of cycles to obtain a set of unique paths connecting a set of cutting points (nodes).

Port of MATLAB's Cycles2Paths.m.

Parameters:

Name	Type	Description	Default
allcycles	sequence of sequence of int	A list of cycle paths.	required
cutpts	sequence of int	Node indices at which each cycle should be cut (cutting points).	required

Returns:

Type	Description
list of list	Unique paths linking the cutting points, sorted by (length, then lexicographically within each length).

Source code in src/tmapper/cycles_to_paths.py

def cycles_to_paths(allcycles, cutpts):
    """Cut a set of cycles to obtain a set of unique paths connecting a
    set of cutting points (nodes).

    Port of MATLAB's ``Cycles2Paths.m``.

    Parameters
    ----------
    allcycles : sequence of sequence of int
        A list of cycle paths.
    cutpts : sequence of int
        Node indices at which each cycle should be cut (cutting points).

    Returns
    -------
    list of list
        Unique paths linking the cutting points, sorted by (length, then
        lexicographically within each length).
    """
    if not allcycles:
        return []

    allpath = []
    for x in allcycles:
        allpath.extend(cycle_cutter(x, cutpts))

    pathlen = [len(p) for p in allpath]
    upathlen = sorted(set(pathlen))

    allupath = []
    for L in upathlen:
        rows = sorted({tuple(p) for p, l in zip(allpath, pathlen) if l == L})
        allupath.extend(list(r) for r in rows)

    return allupath

tmapper.cycle_path_decomp

Port of tmapper_tools/CyclePathDecomp.m from the MATLAB toolbox.

The MATLAB original's diagonal cluster-block overlay (addDiagBlock.m) is not ported literally -- it's tightly coupled to MATLAB axes internals (ax.Children(end).CData) with no Python equivalent. The same visual annotation is reproduced directly here with matplotlib patches, reusing :func:tmapper.graph_utils.find_blocks.

cycle_path_decomp(dg, *, clusterthres=0.5, plotmat=True, plotmds=False, plothist=False, reordermat=True)

Cycle-path decomposition of a directed graph.

Port of MATLAB's CyclePathDecomp.m. Cycles on the graph are first computed, then clustered such that cycles with clusterthres amount of overlap belong to the same cluster. Boundaries between clusters of cycles (nodes where one cycle enters into another) are used to cut cycles into paths, which are then clustered themselves.

Parameters:

Name	Type	Description	Default
dg	DiGraph	The graph to decompose.	required
clusterthres	float	Overlap threshold for cycle/path clustering.	0.5
plotmat	bool	Whether to plot the overlap matrices (with cluster-block outlines) for diagnostics.	True
plotmds	bool	Whether to plot a 2D classical-MDS projection of cycles/paths.	False
plothist	bool	Whether to plot linkage-distance histograms.	False
reordermat	bool	Whether to reorder the overlap matrices by cluster assignment.	True

Returns:

Name	Type	Description
allbd	list	All boundary points between cycle clusters (used to cut cycles into paths).
pclusters_nodes	list of list	Nodes in each path cluster.
pclusters_interior	list of list	Interior (non-boundary) nodes of each path cluster.
pclusters_boundary	list of set	Boundary nodes of each path cluster.
pcluster_conn_dir	list of list of set	Directed boundary connectivity between path clusters.
pclusters_intcrtpts	list of list	Interior critical points of each path cluster.
allupath	list of list	All unique decomposed paths.
Tp	ndarray	Cluster assignment of each path in allupath.

Source code in src/tmapper/cycle_path_decomp.py

def cycle_path_decomp(dg, *, clusterthres=0.5, plotmat=True, plotmds=False,
                       plothist=False, reordermat=True):
    """Cycle-path decomposition of a directed graph.

    Port of MATLAB's ``CyclePathDecomp.m``. Cycles on the graph are
    first computed, then clustered such that cycles with ``clusterthres``
    amount of overlap belong to the same cluster. Boundaries between
    clusters of cycles (nodes where one cycle enters into another) are
    used to cut cycles into paths, which are then clustered themselves.

    Parameters
    ----------
    dg : networkx.DiGraph
        The graph to decompose.
    clusterthres : float, default 0.5
        Overlap threshold for cycle/path clustering.
    plotmat : bool, default True
        Whether to plot the overlap matrices (with cluster-block
        outlines) for diagnostics.
    plotmds : bool, default False
        Whether to plot a 2D classical-MDS projection of cycles/paths.
    plothist : bool, default False
        Whether to plot linkage-distance histograms.
    reordermat : bool, default True
        Whether to reorder the overlap matrices by cluster assignment.

    Returns
    -------
    allbd : list
        All boundary points between cycle clusters (used to cut cycles
        into paths).
    pclusters_nodes : list of list
        Nodes in each path cluster.
    pclusters_interior : list of list
        Interior (non-boundary) nodes of each path cluster.
    pclusters_boundary : list of set
        Boundary nodes of each path cluster.
    pcluster_conn_dir : list of list of set
        Directed boundary connectivity between path clusters.
    pclusters_intcrtpts : list of list
        Interior critical points of each path cluster.
    allupath : list of list
        All unique decomposed paths.
    Tp : numpy.ndarray
        Cluster assignment of each path in ``allupath``.
    """
    A = nx.to_numpy_array(dg, nodelist=list(dg.nodes()), weight=None)
    _, _, _, allcycles = cycle_count2p(A)

    T, ax1 = cycle_cluster(allcycles, clusterthres, plotmat=plotmat, plotmds=plotmds,
                           plothist=plothist, reordermat=reordermat, return_ax=True)
    _, _, _, clusters_boundary, _, _, _ = cycle_cluster_conn(dg, allcycles, T)

    if plotmat and ax1 is not None:
        _add_diag_block(ax1, T, reordermat=reordermat)
    print(f"# cycle clusters = {int(T.max()) if len(T) else 0}")

    allbd = sorted(set().union(*clusters_boundary)) if clusters_boundary else []
    allupath = cycles_to_paths(allcycles, allbd)

    Tp, ax2 = cycle_cluster(allupath, clusterthres, plotmat=plotmat, plotmds=plotmds,
                            plothist=plothist, reordermat=reordermat, return_ax=True)
    (_, pcluster_conn_dir, pclusters_nodes, pclusters_boundary,
     pclusters_interior, _, pclusters_intcrtpts) = cycle_cluster_conn(dg, allupath, Tp)

    if plotmat and ax2 is not None:
        _add_diag_block(ax2, Tp, reordermat=reordermat)
    print(f"# path clusters = {int(Tp.max()) if len(Tp) else 0}")

    return (allbd, pclusters_nodes, pclusters_interior, pclusters_boundary,
            pcluster_conn_dir, pclusters_intcrtpts, allupath, Tp)

tmapper.path_traffic

Port of tmapper_tools/pathtraffic.m from the MATLAB toolbox.

path_traffic(allpath, nodesize)

Compute traffic statistics (based on node size) along each path.

Port of MATLAB's pathtraffic.m. Note: matches MATLAB's std() convention of dividing by N-1 (sample standard deviation), not numpy's default N.

Parameters:

Name	Type	Description	Default
allpath	sequence of sequence of int	Each element is a path of node indices.	required
nodesize	array_like	Size of each node, indexed the same way as the indices in allpath.	required

Returns:

Type	Description
traf_mean, traf_med, traf_min, traf_max, traf_std : numpy.ndarray	One value per path in allpath.

Source code in src/tmapper/path_traffic.py

def path_traffic(allpath, nodesize):
    """Compute traffic statistics (based on node size) along each path.

    Port of MATLAB's ``pathtraffic.m``. Note: matches MATLAB's ``std()``
    convention of dividing by ``N-1`` (sample standard deviation), not
    numpy's default ``N``.

    Parameters
    ----------
    allpath : sequence of sequence of int
        Each element is a path of node indices.
    nodesize : array_like
        Size of each node, indexed the same way as the indices in
        ``allpath``.

    Returns
    -------
    traf_mean, traf_med, traf_min, traf_max, traf_std : numpy.ndarray
        One value per path in ``allpath``.
    """
    nodesize = np.asarray(nodesize, dtype=float)
    path_nodesize = [nodesize[np.asarray(p, dtype=int)] for p in allpath]

    traf_mean = np.array([p.mean() for p in path_nodesize])
    traf_med = np.array([np.median(p) for p in path_nodesize])
    traf_min = np.array([p.min() for p in path_nodesize])
    traf_max = np.array([p.max() for p in path_nodesize])
    traf_std = np.array([p.std(ddof=1) if len(p) > 1 else 0.0 for p in path_nodesize])

    return traf_mean, traf_med, traf_min, traf_max, traf_std

tmapper.qasym(A, C)

Modularity measure of an asymmetric, weighted network.

Port of MATLAB's Qasym.m. Adapted from the canonical modularity definition (Fortunato 2010, Physics Reports) for asymmetric networks with weighted edges.

Parameters:

Name	Type	Description	Default
A	(array_like, shape(N, N))	Adjacency matrix, which may or may not be symmetric.	required
C	(array_like, shape(N))	Community assignment of each node.	required

Returns:

Type	Description
float	Modularity Q. Returns 0 for a zero-edge network (neither modular nor non-modular).

Source code in src/tmapper/modularity.py

def qasym(A, C):
    """Modularity measure of an asymmetric, weighted network.

    Port of MATLAB's ``Qasym.m``. Adapted from the canonical modularity
    definition (Fortunato 2010, Physics Reports) for asymmetric networks
    with weighted edges.

    Parameters
    ----------
    A : array_like, shape (N, N)
        Adjacency matrix, which may or may not be symmetric.
    C : array_like, shape (N,)
        Community assignment of each node.

    Returns
    -------
    float
        Modularity Q. Returns 0 for a zero-edge network (neither modular
        nor non-modular).
    """
    A = np.asarray(A, dtype=float)
    N_edges = A.sum()  # "2m" in other notations

    if N_edges == 0:
        return 0.0

    k_source = A.sum(axis=1, keepdims=True)  # source degree
    k_sink = A.sum(axis=0, keepdims=True)  # sink degree
    P_ij = k_source @ k_sink / N_edges  # null model

    C = np.asarray(C).ravel()
    same_community = C[:, None] == C[None, :]

    return float(np.sum((A - P_ij) * same_community) / N_edges)

tmapper.cal_mod(W, m0)

Modularity score of a (typically symmetric) network given a community assignment.

Port of MATLAB's calMod.m.

Parameters:

Name	Type	Description	Default
W	(array_like, shape(N, N))	Graph adjacency matrix.	required
m0	(array_like, shape(N))	Categorical node labels / community assignment.	required

Returns:

Type	Description
float	Modularity score. Returns 0 for a zero-edge network (neither modular nor non-modular).

Source code in src/tmapper/modularity.py

def cal_mod(W, m0):
    """Modularity score of a (typically symmetric) network given a
    community assignment.

    Port of MATLAB's ``calMod.m``.

    Parameters
    ----------
    W : array_like, shape (N, N)
        Graph adjacency matrix.
    m0 : array_like, shape (N,)
        Categorical node labels / community assignment.

    Returns
    -------
    float
        Modularity score. Returns 0 for a zero-edge network (neither
        modular nor non-modular).
    """
    W = np.asarray(W, dtype=float)
    s = W.sum()  # sum of degrees of nodes (assuming symmetric)

    if s == 0:
        return 0.0

    gamma = 1  # scaling parameter
    B = (W - gamma * (W.sum(axis=1, keepdims=True) @ W.sum(axis=0, keepdims=True)) / s) / s
    B = (B + B.T) / 2  # symmetrize

    m0 = np.asarray(m0).ravel()
    same_community = m0[:, None] == m0[None, :]

    return float(B[same_community].sum())