(X, maxdim=1, thresh=np.inf, coeff=2, metric="euclidean",
           metric_params={}, weights=None, weight_params=None,
           collapse_edges=False, n_perm=None)
return ret
Compute persistence diagrams for X data array using Ripser [1]_.

    If X is not a distance matrix, it will be converted to a distance matrix
    using the chosen metric.

    Parameters
    ----------
    X : ndarray of shape (n_samples, n_features)
        A numpy array of either data or distance matrix. Can also be a sparse
        distance matrix of type scipy.sparse

    maxdim : int, optional, default: ``1``
        Maximum homology dimension computed. Will compute all dimensions lower
        than and equal to this value. For 1, H_0 and H_1 will be computed.

    thresh : float, optional, default: ``numpy.inf``
        Maximum distances considered when constructing filtration. If
        ``numpy.inf``, compute the entire filtration.

    coeff : int prime, optional, default: ``2``
        Compute homology with coefficients in the prime field Z/pZ for p=coeff.

    metric : string or callable, optional, default: ``'euclidean'``
        The metric to use when calculating distance between instances in a
        feature array. If set to ``'precomputed'``, input data is interpreted
        as a distance matrix or of adjacency matrices of a weighted undirected
        graph. If a string, it must be one of the options allowed by
        :func:`scipy.spatial.distance.pdist` for its metric parameter, or a
        or a metric listed in
        :obj:`sklearn.pairwise.PAIRWISE_DISTANCE_FUNCTIONS`, including
        ``'euclidean'``, ``'manhattan'`` or ``'cosine'``. If a callable, it
        should take pairs of vectors (1D arrays) as input and, for each two
        vectors in a pair, it should return a scalar indicating the
        distance/dissimilarity between them.

    metric_params : dict, optional, default: ``{}``
        Additional parameters to be passed to the distance function.

    weights : ``"DTM"``, ndarray or None, optional, default: ``None``
        If not ``None``, the persistence of a weighted Vietoris-Rips filtration
        is computed as described in [3]_, and this parameter determines the
        vertex weights in the modified adjacency matrix. ``"DTM"`` denotes the
        empirical distance-to-measure function defined, following [3]_, by

        .. math:: w(x) = 2\\left\\(\\frac{1}{n+1} \\sum_{k=1}^n
           \\mathrm{dist}(x, x_k)^r \\right)^{1/r}.

        Here, :math:`\\mathrm{dist}` is the distance metric used, :math:`x_k`
        is the :math:`k`-th :math:`\\mathrm{dist}`-nearest neighbour of
        :math:`x` (:math:`x` is not considered a neighbour of itself),
        :math:`n` is the number of nearest neighbors to include, and :math:`r`
        is a parameter (see `weight_params`). If an ndarray is passed, it is
        interpreted as a user-defined list of vertex weights for the modified
        adjacency matrix. In either case, the edge weights
        :math:`\\{w_{ij}\\}_{i, j}` for the modified adjacency matrix are
        computed from the original distances and the new vertex weights
        :math:`\\{w_i\\}_i` as follows:

        .. math:: w_{ij} = \\begin{cases} \\max\\{ w_i, w_j \\}
           &\\text{if } 2\\mathrm{dist}_{ij} \\leq
           |w_i^p - w_j^p|^{\\frac{1}{p}} \\
           t &\\text{otherwise} \\end{cases}

        where :math:`t` is the only positive root of

        .. math:: 2 \\mathrm{dist}_{ij} = (t^p - w_i^p)^\\frac{1}{p} +
           (t^p - w_j^p)^\\frac{1}{p}

        and :math:`p` is a parameter specified in `metric_params`.

    weight_params : dict or None, optional, default: ``None``
        Parameters to be used in the case of weighted filtrations, see
        `weights`. In this case, the key ``"p"`` determines the power to be
        used in computing edge weights from vertex weights. It can be one of
        ``1``, ``2`` or ``np.inf`` and defaults to ``1``. If `weights` is
        ``"DTM"``, the additional keys ``"r"`` (default: ``2``) and
        ``"n_neighbors"`` (default: ``3``) are available (see `weights`,
        where the latter corresponds to :math:`n`).


    collapse_edges : bool, optional, default: ``False``
        Whether to use the edge collapse algorithm as described in [2]_ prior
        to calling ``ripser``.

    n_perm : int or None, optional, default: ``None``
        The number of points to subsample in a "greedy permutation", or a
        furthest point sampling of the points. These points will be used in
        lieu of the full point cloud for a faster computation, at the expense
        of some accuracy, which can be bounded as a maximum bottleneck distance
        to all diagrams on the original point set.

    Returns
    -------
    A dictionary holding all of the results of the computation
    {
        'dgms': list (size maxdim) of ndarray (n_pairs, 2)
            A list of persistence diagrams, one for each dimension less
            than maxdim. Each diagram is an ndarray of size (n_pairs, 2)
            with the first column representing the birth time and the
            second column representing the death time of each pair.
        'num_edges': int
            The number of edges added during the computation
        'dperm2all': None or ndarray (n_perm, n_samples)
            ``None`` if n_perm is ``None``. Otherwise, the distance from all
            points in the permutation to all points in the dataset.
        'idx_perm': ndarray(n_perm) if n_perm > 0
            Index into the original point cloud of the points used
            as a subsample in the greedy permutation
        'r_cover': float
            Covering radius of the subsampled points.
            If n_perm <= 0, then the full point cloud was used and this is 0
    }

    Notes
    -----
    `Ripser <https://github.com/Ripser/ripser>`_ is used as a C++ backend
    for computing Vietoris–Rips persistent homology. Python bindings were
    modified for performance from the `ripser.py
    <https://github.com/scikit-tda/ripser.py>`_ package.

    `GUDHI <https://github.com/GUDHI/gudhi-devel>`_ is used as a C++ backend
    for the edge collapse algorithm described in [2]_.

    References
    ----------
    .. [1] U. Bauer, "Ripser: efficient computation of Vietoris–Rips
           persistence barcodes", 2019; `arXiv:1908.02518
           <https://arxiv.org/abs/1908.02518>`_.

    .. [2] J.-D. Boissonnat and S. Pritam, "Edge Collapse and Persistence of
           Flag Complexes"; in *36th International Symposium on Computational
           Geometry (SoCG 2020)*, pp. 19:1–19:15, Schloss
           Dagstuhl-Leibniz–Zentrum für Informatik, 2020;
           `DOI: 10.4230/LIPIcs.SoCG.2020.19
           <https://doi.org/10.4230/LIPIcs.SoCG.2020.19>`_.

    .. [3] H. Anai et al, "DTM-Based Filtrations"; in *Topological Data
           Analysis* (Abel Symposia, vol 15), Springer, 2020;
           `DOI: 10.1007/978-3-030-43408-3_2
           <https://doi.org/10.1007/978-3-030-43408-3_2>`_.
Compute persistence diagrams for X data array using Ripser [1]_.
Compute persistence diagrams for X data array using Ripser [ 1 ] _ .
def ripser(X, maxdim=1, thresh=np.inf, coeff=2, metric="euclidean",
           metric_params={}, weights=None, weight_params=None,
           collapse_edges=False, n_perm=None):
    """Compute persistence diagrams for X data array using Ripser [1]_.

    If X is not a distance matrix, it will be converted to a distance matrix
    using the chosen metric.

    Parameters
    ----------
    X : ndarray of shape (n_samples, n_features)
        A numpy array of either data or distance matrix. Can also be a sparse
        distance matrix of type scipy.sparse

    maxdim : int, optional, default: ``1``
        Maximum homology dimension computed. Will compute all dimensions lower
        than and equal to this value. For 1, H_0 and H_1 will be computed.

    thresh : float, optional, default: ``numpy.inf``
        Maximum distances considered when constructing filtration. If
        ``numpy.inf``, compute the entire filtration.

    coeff : int prime, optional, default: ``2``
        Compute homology with coefficients in the prime field Z/pZ for p=coeff.

    metric : string or callable, optional, default: ``'euclidean'``
        The metric to use when calculating distance between instances in a
        feature array. If set to ``'precomputed'``, input data is interpreted
        as a distance matrix or of adjacency matrices of a weighted undirected
        graph. If a string, it must be one of the options allowed by
        :func:`scipy.spatial.distance.pdist` for its metric parameter, or a
        or a metric listed in
        :obj:`sklearn.pairwise.PAIRWISE_DISTANCE_FUNCTIONS`, including
        ``'euclidean'``, ``'manhattan'`` or ``'cosine'``. If a callable, it
        should take pairs of vectors (1D arrays) as input and, for each two
        vectors in a pair, it should return a scalar indicating the
        distance/dissimilarity between them.

    metric_params : dict, optional, default: ``{}``
        Additional parameters to be passed to the distance function.

    weights : ``"DTM"``, ndarray or None, optional, default: ``None``
        If not ``None``, the persistence of a weighted Vietoris-Rips filtration
        is computed as described in [3]_, and this parameter determines the
        vertex weights in the modified adjacency matrix. ``"DTM"`` denotes the
        empirical distance-to-measure function defined, following [3]_, by

        .. math:: w(x) = 2\\left\\(\\frac{1}{n+1} \\sum_{k=1}^n
           \\mathrm{dist}(x, x_k)^r \\right)^{1/r}.

        Here, :math:`\\mathrm{dist}` is the distance metric used, :math:`x_k`
        is the :math:`k`-th :math:`\\mathrm{dist}`-nearest neighbour of
        :math:`x` (:math:`x` is not considered a neighbour of itself),
        :math:`n` is the number of nearest neighbors to include, and :math:`r`
        is a parameter (see `weight_params`). If an ndarray is passed, it is
        interpreted as a user-defined list of vertex weights for the modified
        adjacency matrix. In either case, the edge weights
        :math:`\\{w_{ij}\\}_{i, j}` for the modified adjacency matrix are
        computed from the original distances and the new vertex weights
        :math:`\\{w_i\\}_i` as follows:

        .. math:: w_{ij} = \\begin{cases} \\max\\{ w_i, w_j \\}
           &\\text{if } 2\\mathrm{dist}_{ij} \\leq
           |w_i^p - w_j^p|^{\\frac{1}{p}} \\
           t &\\text{otherwise} \\end{cases}

        where :math:`t` is the only positive root of

        .. math:: 2 \\mathrm{dist}_{ij} = (t^p - w_i^p)^\\frac{1}{p} +
           (t^p - w_j^p)^\\frac{1}{p}

        and :math:`p` is a parameter specified in `metric_params`.

    weight_params : dict or None, optional, default: ``None``
        Parameters to be used in the case of weighted filtrations, see
        `weights`. In this case, the key ``"p"`` determines the power to be
        used in computing edge weights from vertex weights. It can be one of
        ``1``, ``2`` or ``np.inf`` and defaults to ``1``. If `weights` is
        ``"DTM"``, the additional keys ``"r"`` (default: ``2``) and
        ``"n_neighbors"`` (default: ``3``) are available (see `weights`,
        where the latter corresponds to :math:`n`).


    collapse_edges : bool, optional, default: ``False``
        Whether to use the edge collapse algorithm as described in [2]_ prior
        to calling ``ripser``.

    n_perm : int or None, optional, default: ``None``
        The number of points to subsample in a "greedy permutation", or a
        furthest point sampling of the points. These points will be used in
        lieu of the full point cloud for a faster computation, at the expense
        of some accuracy, which can be bounded as a maximum bottleneck distance
        to all diagrams on the original point set.

    Returns
    -------
    A dictionary holding all of the results of the computation
    {
        'dgms': list (size maxdim) of ndarray (n_pairs, 2)
            A list of persistence diagrams, one for each dimension less
            than maxdim. Each diagram is an ndarray of size (n_pairs, 2)
            with the first column representing the birth time and the
            second column representing the death time of each pair.
        'num_edges': int
            The number of edges added during the computation
        'dperm2all': None or ndarray (n_perm, n_samples)
            ``None`` if n_perm is ``None``. Otherwise, the distance from all
            points in the permutation to all points in the dataset.
        'idx_perm': ndarray(n_perm) if n_perm > 0
            Index into the original point cloud of the points used
            as a subsample in the greedy permutation
        'r_cover': float
            Covering radius of the subsampled points.
            If n_perm <= 0, then the full point cloud was used and this is 0
    }

    Notes
    -----
    `Ripser <https://github.com/Ripser/ripser>`_ is used as a C++ backend
    for computing Vietoris–Rips persistent homology. Python bindings were
    modified for performance from the `ripser.py
    <https://github.com/scikit-tda/ripser.py>`_ package.

    `GUDHI <https://github.com/GUDHI/gudhi-devel>`_ is used as a C++ backend
    for the edge collapse algorithm described in [2]_.

    References
    ----------
    .. [1] U. Bauer, "Ripser: efficient computation of Vietoris–Rips
           persistence barcodes", 2019; `arXiv:1908.02518
           <https://arxiv.org/abs/1908.02518>`_.

    .. [2] J.-D. Boissonnat and S. Pritam, "Edge Collapse and Persistence of
           Flag Complexes"; in *36th International Symposium on Computational
           Geometry (SoCG 2020)*, pp. 19:1–19:15, Schloss
           Dagstuhl-Leibniz–Zentrum für Informatik, 2020;
           `DOI: 10.4230/LIPIcs.SoCG.2020.19
           <https://doi.org/10.4230/LIPIcs.SoCG.2020.19>`_.

    .. [3] H. Anai et al, "DTM-Based Filtrations"; in *Topological Data
           Analysis* (Abel Symposia, vol 15), Springer, 2020;
           `DOI: 10.1007/978-3-030-43408-3_2
           <https://doi.org/10.1007/978-3-030-43408-3_2>`_.

    """
    if n_perm and issparse(X):
        raise Exception("Greedy permutation is not supported for sparse "
                        "distance matrices")
    if n_perm and n_perm > X.shape[0]:
        raise Exception("Number of points in greedy permutation is greater "
                        "than number of points in the point cloud")
    if n_perm and n_perm < 0:
        raise Exception("There should be a strictly positive number of points "
                        "in the greedy permutation")

    idx_perm = np.arange(X.shape[0])
    r_cover = 0.0
    if n_perm:
        idx_perm, lambdas, dperm2all = \
            get_greedy_perm(X, n_perm=n_perm, metric=metric)
        r_cover = lambdas[-1]
        dm = dperm2all[:, idx_perm]
    else:
        if metric == 'precomputed':
            dm = X
        else:
            dm = pairwise_distances(X, metric=metric, **metric_params)
        dperm2all = None

    n_points = max(dm.shape)

    use_sparse_computer = True
    if issparse(dm):
        row, col, data = _resolve_symmetry_conflicts(dm.tocoo())  # Upper diag

        if weights is not None:
            if (dm < 0).nnz:
                raise ValueError("Distance matrix has negative entries. "
                                 "Weighted Rips filtration unavailable.")

            weight_params = {} if weight_params is None else weight_params
            weights_p = weight_params.get("p", 1)

            # Restrict to off-diagonal entries for weights computation since
            # diagonal ones are given by `weights`. Explicitly set the diagonal
            # to 0 -- this is also important for DTM since otherwise
            # kneighbors_graph with include_self=False skips the first true
            # neighbor.
            off_diag = row != col
            row, col, data = (np.hstack([row[off_diag], np.arange(n_points)]),
                              np.hstack([col[off_diag], np.arange(n_points)]),
                              np.hstack([data[off_diag], np.zeros(n_points)]))

            if isinstance(weights, str) and (weights == "DTM"):
                n_neighbors = weight_params.get("n_neighbors", 3)
                weights_r = weight_params.get("r", 2)

                # CSR matrix must be symmetric for kneighbors_graph to give
                # correct results
                dm = csr_matrix((np.hstack([data, data[:-n_points]]),
                                 (np.hstack([row, col[:-n_points]]),
                                  np.hstack([col, row[:-n_points]]))))
                weights = _compute_dtm_weights(dm, n_neighbors, weights_r)
            else:
                weights = _check_weights(weights, n_points)

            data = _weight_filtration_sparse(row, col, data, weights,
                                             weights_p)

        if collapse_edges:
            row, col, data = _collapse_coo(row, col, data, thresh)

    else:
        if weights is not None:
            if (dm < 0).any():
                raise ValueError("Distance matrix has negative entries. "
                                 "Weighted Rips filtration unavailable.")

            weight_params = {} if weight_params is None else weight_params
            weights_p = weight_params.get("p", 1)

            if isinstance(weights, str) and (weights == "DTM"):
                n_neighbors = weight_params.get("n_neighbors", 3)
                weights_r = weight_params.get("r", 2)

                if not np.array_equal(dm, dm.T):
                    dm = np.triu(dm, k=1)
                    dm += dm.T

                weights = _compute_dtm_weights(dm, n_neighbors, weights_r)
            else:
                weights = _check_weights(weights, n_points)

            dm = _weight_filtration_dense(dm, weights, weights_p)
            np.fill_diagonal(dm, weights)

        if (dm.diagonal() != 0).any():
            # Convert to sparse format, because currently that's the only
            # one handling nonzero births
            (row, col) = np.triu_indices_from(dm)
            data = dm[(row, col)]
            if collapse_edges:
                row, col, data = _collapse_coo(row, col, data, thresh)
        elif collapse_edges:
            row, col, data = gtda_collapser.\
                flag_complex_collapse_edges_dense(dm, thresh)
        else:
            use_sparse_computer = False

    if use_sparse_computer:
        res = DRFDMSparse(np.asarray(row, dtype=np.int32, order="C"),
                          np.asarray(col, dtype=np.int32, order="C"),
                          np.asarray(data, dtype=np.float32, order="C"),
                          n_points,
                          maxdim,
                          thresh,
                          coeff)
    else:
        # Only consider strict upper diagonal
        DParam = squareform(dm, checks=False).astype(np.float32)
        # Run garbage collector to free up memory taken by `dm`
        del dm
        gc.collect()
        res = DRFDM(DParam, maxdim, thresh, coeff)

    # Unwrap persistence diagrams
    dgms = res.births_and_deaths_by_dim
    for dim in range(len(dgms)):
        N = int(len(dgms[dim]) / 2)
        dgms[dim] = np.reshape(np.array(dgms[dim]), [N, 2])

    ret = {"dgms": dgms,
           "num_edges": res.num_edges,
           "dperm2all": dperm2all,
           "idx_perm": idx_perm,
           "r_cover": r_cover}

    return ret
def ripser ( X , maxdim = 1 , thresh = np . inf , coeff = 2 , metric = "euclidean" , metric_params = { } , weights = None , weight_params = None , collapse_edges = False , n_perm = None ) : if n_perm and issparse ( X ) : raise Exception ( "Greedy permutation is not supported for sparse " "distance matrices" ) if n_perm and n_perm > X . shape [ 0 ] : raise Exception ( "Number of points in greedy permutation is greater " "than number of points in the point cloud" ) if n_perm and n_perm < 0 : raise Exception ( "There should be a strictly positive number of points " "in the greedy permutation" ) idx_perm = np . arange ( X . shape [ 0 ] ) r_cover = 0.0 if n_perm : idx_perm , lambdas , dperm2all = get_greedy_perm ( X , n_perm = n_perm , metric = metric ) r_cover = lambdas [ - 1 ] dm = dperm2all [ : , idx_perm ] else : if metric == 'precomputed' : dm = X else : dm = pairwise_distances ( X , metric = metric , * * metric_params ) dperm2all = None n_points = max ( dm . shape ) use_sparse_computer = True if issparse ( dm ) : row , col , data = _resolve_symmetry_conflicts ( dm . tocoo ( ) ) # Upper diag if weights is not None : if ( dm < 0 ) . nnz : raise ValueError ( "Distance matrix has negative entries. " "Weighted Rips filtration unavailable." ) weight_params = { } if weight_params is None else weight_params weights_p = weight_params . get ( "p" , 1 ) # Restrict to off-diagonal entries for weights computation since # diagonal ones are given by `weights`. Explicitly set the diagonal # to 0 -- this is also important for DTM since otherwise # kneighbors_graph with include_self=False skips the first true # neighbor. off_diag = row != col row , col , data = ( np . hstack ( [ row [ off_diag ] , np . arange ( n_points ) ] ) , np . hstack ( [ col [ off_diag ] , np . arange ( n_points ) ] ) , np . hstack ( [ data [ off_diag ] , np . zeros ( n_points ) ] ) ) if isinstance ( weights , str ) and ( weights == "DTM" ) : n_neighbors = weight_params . get ( "n_neighbors" , 3 ) weights_r = weight_params . get ( "r" , 2 ) # CSR matrix must be symmetric for kneighbors_graph to give # correct results dm = csr_matrix ( ( np . hstack ( [ data , data [ : - n_points ] ] ) , ( np . hstack ( [ row , col [ : - n_points ] ] ) , np . hstack ( [ col , row [ : - n_points ] ] ) ) ) ) weights = _compute_dtm_weights ( dm , n_neighbors , weights_r ) else : weights = _check_weights ( weights , n_points ) data = _weight_filtration_sparse ( row , col , data , weights , weights_p ) if collapse_edges : row , col , data = _collapse_coo ( row , col , data , thresh ) else : if weights is not None : if ( dm < 0 ) . any ( ) : raise ValueError ( "Distance matrix has negative entries. " "Weighted Rips filtration unavailable." ) weight_params = { } if weight_params is None else weight_params weights_p = weight_params . get ( "p" , 1 ) if isinstance ( weights , str ) and ( weights == "DTM" ) : n_neighbors = weight_params . get ( "n_neighbors" , 3 ) weights_r = weight_params . get ( "r" , 2 ) if not np . array_equal ( dm , dm . T ) : dm = np . triu ( dm , k = 1 ) dm += dm . T weights = _compute_dtm_weights ( dm , n_neighbors , weights_r ) else : weights = _check_weights ( weights , n_points ) dm = _weight_filtration_dense ( dm , weights , weights_p ) np . fill_diagonal ( dm , weights ) if ( dm . diagonal ( ) != 0 ) . any ( ) : # Convert to sparse format, because currently that's the only # one handling nonzero births ( row , col ) = np . triu_indices_from ( dm ) data = dm [ ( row , col ) ] if collapse_edges : row , col , data = _collapse_coo ( row , col , data , thresh ) elif collapse_edges : row , col , data = gtda_collapser . flag_complex_collapse_edges_dense ( dm , thresh ) else : use_sparse_computer = False if use_sparse_computer : res = DRFDMSparse ( np . asarray ( row , dtype = np . int32 , order = "C" ) , np . asarray ( col , dtype = np . int32 , order = "C" ) , np . asarray ( data , dtype = np . float32 , order = "C" ) , n_points , maxdim , thresh , coeff ) else : # Only consider strict upper diagonal DParam = squareform ( dm , checks = False ) . astype ( np . float32 ) # Run garbage collector to free up memory taken by `dm` del dm gc . collect ( ) res = DRFDM ( DParam , maxdim , thresh , coeff ) # Unwrap persistence diagrams dgms = res . births_and_deaths_by_dim for dim in range ( len ( dgms ) ) : N = int ( len ( dgms [ dim ] ) / 2 ) dgms [ dim ] = np . reshape ( np . array ( dgms [ dim ] ) , [ N , 2 ] ) ret = { "dgms" : dgms , "num_edges" : res . num_edges , "dperm2all" : dperm2all , "idx_perm" : idx_perm , "r_cover" : r_cover } return ret