From 9d5cfbc144ae8fb9d83ae59b8e37b069c462fd5d Mon Sep 17 00:00:00 2001 From: Umberto Lupo Date: Wed, 3 Feb 2021 10:36:41 +0100 Subject: [PATCH 01/15] First draft of lower-star persistence and support for extended persistence --- doc/modules/homology.rst | 1 + doc/modules/plotting.rst | 1 + gtda/homology/__init__.py | 3 +- gtda/homology/_utils.py | 28 ++-- gtda/homology/simplicial.py | 179 +++++++++++++++++++++++++- gtda/plotting/__init__.py | 3 +- gtda/plotting/persistence_diagrams.py | 156 +++++++++++++++++++--- 7 files changed, 343 insertions(+), 28 deletions(-) diff --git a/doc/modules/homology.rst b/doc/modules/homology.rst index 36c38b41f..b65a9096c 100644 --- a/doc/modules/homology.rst +++ b/doc/modules/homology.rst @@ -18,6 +18,7 @@ Undirected simplicial homology homology.SparseRipsPersistence homology.WeakAlphaPersistence homology.EuclideanCechPersistence + homology.LowerStarFlagPersistence Directed simplicial homology ---------------------------- diff --git a/doc/modules/plotting.rst b/doc/modules/plotting.rst index 2f8154057..045730b0d 100644 --- a/doc/modules/plotting.rst +++ b/doc/modules/plotting.rst @@ -14,5 +14,6 @@ plotting.plot_point_cloud plotting.plot_heatmap plotting.plot_diagram + plotting.plot_extended_diagram plotting.plot_betti_curves plotting.plot_betti_surfaces diff --git a/gtda/homology/__init__.py b/gtda/homology/__init__.py index c491ebe45..6de39435c 100644 --- a/gtda/homology/__init__.py +++ b/gtda/homology/__init__.py @@ -4,7 +4,7 @@ from .simplicial import VietorisRipsPersistence, WeightedRipsPersistence, \ SparseRipsPersistence, WeakAlphaPersistence, EuclideanCechPersistence, \ - FlagserPersistence + FlagserPersistence, LowerStarFlagPersistence from .cubical import CubicalPersistence __all__ = [ @@ -14,5 +14,6 @@ 'WeakAlphaPersistence', 'EuclideanCechPersistence', 'FlagserPersistence', + 'LowerStarFlagPersistence', 'CubicalPersistence', ] diff --git a/gtda/homology/_utils.py b/gtda/homology/_utils.py index 007d3ef58..f3259f1e4 100644 --- a/gtda/homology/_utils.py +++ b/gtda/homology/_utils.py @@ -21,6 +21,7 @@ def replace_infinity_values(subdiagram): for dim in homology_dimensions} Xt = [{dim: replace_infinity_values(diagram[dim][slices[dim]]) for dim in homology_dimensions} for diagram in Xt] + feature_vect_len = 3 elif format == "gudhi": # Input is list of list of [dim, (birth, death)] # In H0, remove one infinite bar placed at the beginning by GUDHI only # if `reduce` is True @@ -31,23 +32,28 @@ def replace_infinity_values(subdiagram): if pers_info[0] == dim]).reshape(-1, 2)[slices[dim]] ) for dim in homology_dimensions} for diagram in Xt] + feature_vect_len = 3 + elif format == "extended": # Input is list of list of subdiagrams + Xt = [{dim: diagram[dim] + for dim in homology_dimensions} for diagram in Xt] + feature_vect_len = 4 else: - raise ValueError( - f"Unknown input format {format} for collection of diagrams." - ) + raise ValueError(f"Unknown input format {format} for collection of " + f"diagrams.") - # Conversion to array of triples with padding triples + # Conversion to array of triples/quadruples with padding start_idx_per_dim = np.cumsum( - [0] + [np.max([len(diagram[dim]) for diagram in Xt] + [1]) - for dim in homology_dimensions] - ) + [0] + [np.max([len(diagram[dim]) for diagram in Xt] + [1]) + for dim in homology_dimensions] + ) min_values = [min([np.min(diagram[dim][:, 0]) if diagram[dim].size else np.inf for diagram in Xt]) for dim in homology_dimensions] min_values = [min_value if min_value != np.inf else 0 for min_value in min_values] n_features = start_idx_per_dim[-1] - Xt_padded = np.empty((len(Xt), n_features, 3), dtype=float) + Xt_padded = np.empty((len(Xt), n_features, feature_vect_len), dtype=float) + Xt_padded[:, :, 3:] = 1. # Only applies to extended persistence for i, dim in enumerate(homology_dimensions): start_idx, end_idx = start_idx_per_dim[i:i + 2] @@ -58,8 +64,10 @@ def replace_infinity_values(subdiagram): subdiagram = diagram[dim] end_idx_nontrivial = start_idx + len(subdiagram) # Populate nontrivial part of the subdiagram - Xt_padded[j, start_idx:end_idx_nontrivial, :2] = subdiagram - # Insert padding triples + Xt_padded[j, start_idx:end_idx_nontrivial, :2] = subdiagram[:, :2] + Xt_padded[j, start_idx:end_idx_nontrivial, 3:] = \ + subdiagram[:, 2:] # Only applies to extended persistence + # Insert padding triples/quadruples Xt_padded[j, end_idx_nontrivial:end_idx, :2] = [padding_value] * 2 return Xt_padded diff --git a/gtda/homology/simplicial.py b/gtda/homology/simplicial.py index 18702913d..a4b69d9bb 100644 --- a/gtda/homology/simplicial.py +++ b/gtda/homology/simplicial.py @@ -16,7 +16,7 @@ from ._utils import _postprocess_diagrams from ..base import PlotterMixin from ..externals.python import ripser, SparseRipsComplex, CechComplex -from ..plotting import plot_diagram +from ..plotting import plot_diagram, plot_extended_diagram from ..utils._docs import adapt_fit_transform_docs from ..utils.intervals import Interval from ..utils.validation import validate_params, check_point_clouds @@ -1772,3 +1772,180 @@ def plot(Xt, sample=0, homology_dimensions=None, plotly_params=None): Xt[sample], homology_dimensions=homology_dimensions, plotly_params=plotly_params ) + + +# @adapt_fit_transform_docs +class LowerStarFlagPersistence(BaseEstimator, TransformerMixin, PlotterMixin): + """TODO""" + _hyperparameters = { + "homology_dimensions": { + "type": (list, tuple), + "of": {"type": int, "in": Interval(0, np.inf, closed="left")} + }, + "coeff": {"type": int, "in": Interval(2, np.inf, closed="left")}, + "extended": {"type": bool}, + "infinity_values": {"type": (Real, type(None))}, + "reduced_homology": {"type": bool}, + "collapse_edges": {"type": bool} + } + + def __init__(self, homology_dimensions=(0, 1), coeff=2, extended=True, + infinity_values=np.inf, reduced_homology=False, + collapse_edges=False, n_jobs=None): + self.homology_dimensions = homology_dimensions + self.coeff = coeff + self.extended = extended + self.infinity_values = infinity_values + self.reduced_homology = reduced_homology + self.collapse_edges = collapse_edges + self.n_jobs = n_jobs + + def _lower_star_extended_diagram(self, X): + n_points = max(X.shape) + n_points_cone = n_points + 1 + n_diag = min(X.shape) + X = X.tocoo() + + data_diag = np.zeros(n_points_cone, dtype=X.dtype) + data_diag[:n_diag] = X.diagonal() + max_value = data_diag[:n_diag].max() + min_value = data_diag[:n_diag].min() + data_diag[-1] = min_value - 1 + + off_diag = X.row != X.col + row_off_diag, col_off_diag = X.row[off_diag], X.col[off_diag] + + row = np.concatenate([row_off_diag, np.arange(n_points)]) + col = np.concatenate([col_off_diag, np.full(n_points, n_points)]) + data = np.concatenate([ + np.maximum(data_diag[row_off_diag], data_diag[col_off_diag]), + 2 * max_value + 1 - data_diag[:n_points] + ]) + + X = coo_matrix((data, (row, col)), + shape=(n_points_cone, n_points_cone)) + X.setdiag(data_diag) + + Xdgms = ripser( + X, metric="precomputed", maxdim=self._max_homology_dimension, + coeff=self.coeff, collapse_edges=self.collapse_edges + )["dgms"] + + for i in range(len(Xdgms)): + mask_down_sweep = Xdgms[i] > max_value + sgn = 2 * np.logical_not( + np.logical_xor.reduce(mask_down_sweep, axis=1, + keepdims=True)).astype(int) - 1 + Xdgms[i][mask_down_sweep] = \ + 2 * max_value + 1 - Xdgms[i][mask_down_sweep] + if not i: + Xdgms[i] = np.hstack([Xdgms[i][:-1, :], sgn[:-1, :]]) + else: + Xdgms[i] = np.hstack([Xdgms[i], sgn]) + + return Xdgms + + def _lower_star_diagram(self, X): + n_points = max(X.shape) + n_diag = min(X.shape) + X = X.tocoo() + + data_diag = np.zeros(n_points, dtype=X.dtype) + data_diag[:n_diag] = X.diagonal() + + off_diag = X.row != X.col + row_off_diag, col_off_diag = X.row[off_diag], X.col[off_diag] + + row = np.concatenate([row_off_diag, np.arange(n_points)]) + col = np.concatenate([col_off_diag, np.arange(n_points)]) + data = np.concatenate([ + np.maximum(data_diag[row_off_diag], data_diag[col_off_diag]), + data_diag] + ) + + X = coo_matrix((data, (row, col))) + + Xdgms = ripser( + X, metric="precomputed", maxdim=self._max_homology_dimension, + coeff=self.coeff, collapse_edges=self.collapse_edges + )["dgms"] + + return Xdgms + + def fit(self, X, y=None): + """TODO""" + check_point_clouds(X, accept_sparse=True, distance_matrices=True) + validate_params( + self.get_params(), self._hyperparameters, exclude=["n_jobs"]) + + self._homology_dimensions = sorted(self.homology_dimensions) + self._max_homology_dimension = self._homology_dimensions[-1] + + if self.extended: + self._diagram_computer = self._lower_star_extended_diagram + self._format = "extended" + else: + self._diagram_computer = self._lower_star_diagram + self._format = "ripser" + + self._is_fitted = True + + return self + + def transform(self, X, y=None): + """TODO""" + check_is_fitted(self, "_is_fitted") + X = check_point_clouds(X, accept_sparse=True, distance_matrices=True) + + Xt = Parallel(n_jobs=self.n_jobs)( + delayed(self._diagram_computer)(x) for x in X) + + format = "extended" if self.extended else "ripser" + Xt = _postprocess_diagrams( + Xt, format, self._homology_dimensions, self.infinity_values, + self.reduced_homology + ) + return Xt + + @staticmethod + def plot(Xt, sample=0, homology_dimensions=None, plotly_params=None): + """Plot a sample from a collection of persistence diagrams, with + homology in multiple dimensions. + + Parameters + ---------- + Xt : ndarray of shape (n_samples, n_features, 3) + Collection of persistence diagrams, such as returned by + :meth:`transform`. + + sample : int, optional, default: ``0`` + Index of the sample in `Xt` to be plotted. + + homology_dimensions : list, tuple or None, optional, default: ``None`` + Which homology dimensions to include in the plot. ``None`` means + plotting all dimensions present in ``Xt[sample]``. + + plotly_params : dict or None, optional, default: ``None`` + Custom parameters to configure the plotly figure. Allowed keys are + ``"traces"`` and ``"layout"``, and the corresponding values should + be dictionaries containing keyword arguments as would be fed to the + :meth:`update_traces` and :meth:`update_layout` methods of + :class:`plotly.graph_objects.Figure`. + + Returns + ------- + fig : :class:`plotly.graph_objects.Figure` object + Plotly figure. + + """ + Xt_sample = Xt[sample] + if Xt_sample.shape[1] == 4: + return plot_extended_diagram( + Xt_sample, homology_dimensions=homology_dimensions, + plotly_params=plotly_params + ) + + return plot_diagram( + Xt_sample, homology_dimensions=homology_dimensions, + plotly_params=plotly_params + ) diff --git a/gtda/plotting/__init__.py b/gtda/plotting/__init__.py index c98e17011..71ffee549 100644 --- a/gtda/plotting/__init__.py +++ b/gtda/plotting/__init__.py @@ -2,13 +2,14 @@ giotto-tda transformers.""" from .point_clouds import plot_point_cloud -from .persistence_diagrams import plot_diagram +from .persistence_diagrams import plot_diagram, plot_extended_diagram from .diagram_representations import plot_betti_curves, plot_betti_surfaces from .images import plot_heatmap __all__ = [ 'plot_point_cloud', 'plot_diagram', + 'plot_extended_diagram', 'plot_heatmap', 'plot_betti_curves', 'plot_betti_surfaces' diff --git a/gtda/plotting/persistence_diagrams.py b/gtda/plotting/persistence_diagrams.py index a01fd4afd..c860b87d0 100644 --- a/gtda/plotting/persistence_diagrams.py +++ b/gtda/plotting/persistence_diagrams.py @@ -2,7 +2,10 @@ # License: GNU AGPLv3 import numpy as np -import plotly.graph_objs as gobj +import plotly.graph_objs as go +from plotly.subplots import make_subplots + +from plotly.colors import DEFAULT_PLOTLY_COLORS as default_colors def plot_diagram(diagram, homology_dimensions=None, plotly_params=None): @@ -36,12 +39,12 @@ def plot_diagram(diagram, homology_dimensions=None, plotly_params=None): if homology_dimensions is None: homology_dimensions = np.unique(diagram[:, 2]) - diagram = diagram[diagram[:, 0] != diagram[:, 1]] - diagram_no_dims = diagram[:, :2] - posinfinite_mask = np.isposinf(diagram_no_dims) - neginfinite_mask = np.isneginf(diagram_no_dims) - max_val = np.max(np.where(posinfinite_mask, -np.inf, diagram_no_dims)) - min_val = np.min(np.where(neginfinite_mask, np.inf, diagram_no_dims)) + diag = diagram[diagram[:, 0] != diagram[:, 1]] + diag_no_dims = diag[:, :2] + posinfinite_mask = np.isposinf(diag_no_dims) + neginfinite_mask = np.isneginf(diag_no_dims) + max_val = np.max(np.where(posinfinite_mask, -np.inf, diag_no_dims)) + min_val = np.min(np.where(neginfinite_mask, np.inf, diag_no_dims)) parameter_range = max_val - min_val extra_space_factor = 0.02 has_posinfinite_death = np.any(posinfinite_mask[:, 1]) @@ -52,8 +55,8 @@ def plot_diagram(diagram, homology_dimensions=None, plotly_params=None): min_val_display = min_val - extra_space max_val_display = max_val + extra_space - fig = gobj.Figure() - fig.add_trace(gobj.Scatter( + fig = go.Figure() + fig.add_trace(go.Scatter( x=[min_val_display, max_val_display], y=[min_val_display, max_val_display], mode="lines", @@ -64,9 +67,9 @@ def plot_diagram(diagram, homology_dimensions=None, plotly_params=None): for dim in homology_dimensions: name = f"H{int(dim)}" if dim != np.inf else "Any homology dimension" - subdiagram = diagram[diagram[:, 2] == dim] + subdiag = diag[diag[:, 2] == dim] unique, inverse, counts = np.unique( - subdiagram, axis=0, return_inverse=True, return_counts=True + subdiag, axis=0, return_inverse=True, return_counts=True ) hovertext = [ f"{tuple(unique[unique_row_index][:2])}" + @@ -76,11 +79,11 @@ def plot_diagram(diagram, homology_dimensions=None, plotly_params=None): ) for unique_row_index in inverse ] - y = subdiagram[:, 1] + y = subdiag[:, 1] if has_posinfinite_death: y[np.isposinf(y)] = posinfinity_val - fig.add_trace(gobj.Scatter( - x=subdiagram[:, 0], y=y, mode="markers", + fig.add_trace(go.Scatter( + x=subdiag[:, 0], y=y, mode="markers", hoverinfo="text", hovertext=hovertext, name=name )) @@ -122,7 +125,7 @@ def plot_diagram(diagram, homology_dimensions=None, plotly_params=None): # Add a horizontal dashed line for points with infinite death if has_posinfinite_death: - fig.add_trace(gobj.Scatter( + fig.add_trace(go.Scatter( x=[min_val_display, max_val_display], y=[posinfinity_val, posinfinity_val], mode="lines", @@ -138,3 +141,126 @@ def plot_diagram(diagram, homology_dimensions=None, plotly_params=None): fig.update_layout(plotly_params.get("layout", None)) return fig + + +def plot_extended_diagram(diagram, homology_dimensions=None, + plotly_params=None): + """Plot a single extended persistence diagram. + + Parameters + ---------- + diagram : ndarray of shape (n_points, 4) + The persistence diagram to plot, where the third dimension along axis 1 + contains homology dimensions, the fourth is either 1 or -1 depending on + whether the feature was born and died during the same phase (upwards or + downwards sweep) and the first two contain (birth, death) pairs to be + used as coordinates in the two-dimensional plot. + + homology_dimensions : list of int or None, optional, default: ``None`` + Homology dimensions which will appear on the plot. If ``None``, all + homology dimensions which appear in `diagram` will be plotted. + + plotly_params : dict or None, optional, default: ``None`` + Custom parameters to configure the plotly figure. Allowed keys are + ``"traces"`` and ``"layout"``, and the corresponding values should be + dictionaries containing keyword arguments as would be fed to the + :meth:`update_traces` and :meth:`update_layout` methods of + :class:`plotly.graph_objects.Figure`. + + Returns + ------- + fig : :class:`plotly.graph_objects.Figure` object + Figure representing the persistence diagram. + + """ + # TODO: increase the marker size + if homology_dimensions is None: + homology_dimensions = np.unique(diagram[:, 2]) + + diag = diagram[diagram[:, 0] != diagram[:, 1]] + diag_no_dims = diag[:, :2] + posinfinite_mask = np.isposinf(diag_no_dims) + neginfinite_mask = np.isneginf(diag_no_dims) + max_val = np.max(np.where(posinfinite_mask, -np.inf, diag_no_dims)) + min_val = np.min(np.where(neginfinite_mask, np.inf, diag_no_dims)) + parameter_range = max_val - min_val + extra_space_factor = 0.02 + extra_space = extra_space_factor * parameter_range + min_val_display = min_val - extra_space + max_val_display = max_val + extra_space + + subplot_titles = ["Same sweep", "Different sweep"] + fig = make_subplots(rows=2, cols=1, + subplot_titles=subplot_titles) + + for row in [1, 2]: + fig.add_trace(go.Scatter( + x=[min_val_display, max_val_display], + y=[min_val_display, max_val_display], + mode="lines", + line={"dash": "dash", "width": 1, "color": "black"}, + showlegend=False, + hoverinfo="none" + ), row=row, col=1) + + for i, sweep in enumerate([1., -1.]): + diag_sweep = diag[diag[:, 3] == sweep] + for j, dim in enumerate(homology_dimensions): + name = f"H{int(dim)}" \ + if dim != np.inf else "Any homology dimension" + subdiag = diag_sweep[diag_sweep[:, 2] == dim] + unique, inverse, counts = np.unique( + subdiag, axis=0, return_inverse=True, return_counts=True + ) + hovertext = [ + f"{tuple(unique[unique_row_index][:2])}" + + ( + f", multiplicity: {counts[unique_row_index]}" + if counts[unique_row_index] > 1 else "" + ) + for unique_row_index in inverse + ] + fig.add_trace(go.Scatter( + x=subdiag[:, 0], y=subdiag[:, 1], mode="markers", + marker={"color": default_colors[j]}, hoverinfo="text", + hovertext=hovertext, name=name + ), row=i + 1, col=1) + + fig.update_layout(width=500, + height=1000, + plot_bgcolor="white") + fig.update_xaxes(title="Birth", + side="bottom", + type="linear", + range=[min_val_display, max_val_display], + autorange=False, + ticks="outside", + showline=True, + zeroline=True, + linewidth=1, + linecolor="black", + mirror=False, + showexponent="all", + exponentformat="e") + fig.update_yaxes(title="Death", + side="left", + type="linear", + range=[min_val_display, max_val_display], + autorange=False, + scaleanchor="x", + scaleratio=1, + ticks="outside", + showline=True, + zeroline=True, + linewidth=1, + linecolor="black", + mirror=False, + showexponent="all", + exponentformat="e") + + # Update traces and layout according to user input + if plotly_params: + fig.update_traces(plotly_params.get("traces", None)) + fig.update_layout(plotly_params.get("layout", None)) + + return fig From 1374f861ac1479f4b001530a699fb5c2601f0787 Mon Sep 17 00:00:00 2001 From: Umberto Lupo <46537483+ulupo@users.noreply.github.com> Date: Fri, 5 Feb 2021 12:09:41 +0100 Subject: [PATCH 02/15] Fix error following @wreise's review --- gtda/homology/simplicial.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/gtda/homology/simplicial.py b/gtda/homology/simplicial.py index a4b69d9bb..977d3fe92 100644 --- a/gtda/homology/simplicial.py +++ b/gtda/homology/simplicial.py @@ -1900,9 +1900,8 @@ def transform(self, X, y=None): Xt = Parallel(n_jobs=self.n_jobs)( delayed(self._diagram_computer)(x) for x in X) - format = "extended" if self.extended else "ripser" Xt = _postprocess_diagrams( - Xt, format, self._homology_dimensions, self.infinity_values, + Xt, self._format, self._homology_dimensions, self.infinity_values, self.reduced_homology ) return Xt From d749f4e9d373502c2b6e01ea327e7396268c34e8 Mon Sep 17 00:00:00 2001 From: wreise Date: Mon, 8 Feb 2021 20:40:54 +0100 Subject: [PATCH 03/15] Add a particular test for extended persistence --- gtda/homology/tests/test_simplicial.py | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py index a2c7d9017..43d87ff73 100644 --- a/gtda/homology/tests/test_simplicial.py +++ b/gtda/homology/tests/test_simplicial.py @@ -5,14 +5,14 @@ import plotly.io as pio import pytest from numpy.testing import assert_almost_equal -from scipy.sparse import csr_matrix +from scipy.sparse import csr_matrix, coo_matrix from scipy.spatial.distance import pdist, squareform from scipy.spatial.qhull import QhullError from sklearn.exceptions import NotFittedError from gtda.homology import VietorisRipsPersistence, WeightedRipsPersistence, \ SparseRipsPersistence, WeakAlphaPersistence, EuclideanCechPersistence, \ - FlagserPersistence + FlagserPersistence, LowerStarFlagPersistence pio.renderers.default = 'plotly_mimetype' @@ -466,3 +466,23 @@ def test_fp_fit_transform_plot(X, hom_dims): FlagserPersistence(directed=False).fit_transform_plot( X_dist, sample=0, homology_dimensions=hom_dims ) + + +X_lsp_cp = coo_matrix((np.array([1, 2, -1, 3, -2, 0.5, + 1, 1, 1, 1, 1, 1]), + (np.array([0, 1, 2, 3, 4, 5, + 0, 1, 2, 3, 4, 0]), + np.array([0, 1, 2, 3, 4, 5, + 1, 2, 3, 4, 5, 5]))) + ) +diag_lsp_cp = np.array([[-2, 3, 0, -1], + [-1, 2, 0, 1], + [2, -1, 1, 1], + [3, -2, 1, 1]]) + + +def test_lsp_fit_transform(): + lp = LowerStarFlagPersistence() + result = lp.fit_transform([X_lsp_cp])[0] + assert_almost_equal(np.sort(result, axis=1), + np.sort(diag_lsp_cp, axis=1)) From 755979c6cd1ff685015b72c051d69c45100a2307 Mon Sep 17 00:00:00 2001 From: wreise Date: Mon, 8 Feb 2021 21:11:38 +0100 Subject: [PATCH 04/15] Correct sign in test --- gtda/homology/tests/test_simplicial.py | 36 ++++++++++++++------------ 1 file changed, 19 insertions(+), 17 deletions(-) diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py index 43d87ff73..840c37869 100644 --- a/gtda/homology/tests/test_simplicial.py +++ b/gtda/homology/tests/test_simplicial.py @@ -113,7 +113,7 @@ def test_vrp_low_infinity_values(X, metric): def test_vrp_fit_transform_plot(X, metric, hom_dims): VietorisRipsPersistence(metric=metric).fit_transform_plot( X, sample=0, homology_dimensions=hom_dims - ) + ) def test_wrp_params(): @@ -215,7 +215,7 @@ def test_wrp_infinity_error(): def test_wrp_fit_transform_plot(X, metric, hom_dims): WeightedRipsPersistence( metric=metric, weight_params={'n_neighbors': 1} - ).fit_transform_plot(X, sample=0, homology_dimensions=hom_dims) + ).fit_transform_plot(X, sample=0, homology_dimensions=hom_dims) def test_srp_params(): @@ -338,7 +338,7 @@ def test_wap_low_infinity_values(X): def test_wap_fit_transform_plot(X, hom_dims): WeakAlphaPersistence().fit_transform_plot( X, sample=0, homology_dimensions=hom_dims - ) + ) def test_cp_params(): @@ -380,7 +380,7 @@ def test_cp_transform(X): def test_cp_fit_transform_plot(X, hom_dims): EuclideanCechPersistence().fit_transform_plot( X, sample=0, homology_dimensions=hom_dims - ) + ) def test_fp_params(): @@ -410,7 +410,7 @@ def test_fp_not_fitted(): [0., 0.34546959, 0.], [0., 0.40065941, 0.], [0., 0.43094373, 0.], - [0.5117411, 0.51976681, 1.]]]) + [0.5117411, 0.51976681, 1.]]]) @pytest.mark.parametrize('X', @@ -465,24 +465,26 @@ def test_fp_transform_high_hom_dim(delta): def test_fp_fit_transform_plot(X, hom_dims): FlagserPersistence(directed=False).fit_transform_plot( X_dist, sample=0, homology_dimensions=hom_dims - ) + ) X_lsp_cp = coo_matrix((np.array([1, 2, -1, 3, -2, 0.5, - 1, 1, 1, 1, 1, 1]), - (np.array([0, 1, 2, 3, 4, 5, - 0, 1, 2, 3, 4, 0]), - np.array([0, 1, 2, 3, 4, 5, - 1, 2, 3, 4, 5, 5]))) - ) + 1, 1, 1, 1, 1, 1]), + (np.array([0, 1, 2, 3, 4, 5, + 0, 1, 2, 3, 4, 0]), + np.array([0, 1, 2, 3, 4, 5, + 1, 2, 3, 4, 5, 5]))) + ) + diag_lsp_cp = np.array([[-2, 3, 0, -1], - [-1, 2, 0, 1], - [2, -1, 1, 1], - [3, -2, 1, 1]]) + [-1, 2, 0, 1], + [2, -1, 1, 1], + [3, -2, 1, -1]], dtype=float) def test_lsp_fit_transform(): lp = LowerStarFlagPersistence() result = lp.fit_transform([X_lsp_cp])[0] - assert_almost_equal(np.sort(result, axis=1), - np.sort(diag_lsp_cp, axis=1)) + assert_almost_equal(np.sort(result, axis=0), + np.sort(diag_lsp_cp, axis=0)) + From f6e7427088ca0a474c4fec1d98f0807e03d1c5bd Mon Sep 17 00:00:00 2001 From: wreise Date: Tue, 9 Feb 2021 10:11:13 +0100 Subject: [PATCH 05/15] Propose a (failing) duality test --- gtda/homology/tests/test_simplicial.py | 48 +++++++++++++++++++++++++- 1 file changed, 47 insertions(+), 1 deletion(-) diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py index 840c37869..6cf845d49 100644 --- a/gtda/homology/tests/test_simplicial.py +++ b/gtda/homology/tests/test_simplicial.py @@ -4,6 +4,9 @@ import numpy as np import plotly.io as pio import pytest +from hypothesis import given +from hypothesis.extra.numpy import arrays, array_shapes +from hypothesis.strategies import composite, integers, floats, booleans, sets, tuples, lists from numpy.testing import assert_almost_equal from scipy.sparse import csr_matrix, coo_matrix from scipy.spatial.distance import pdist, squareform @@ -483,8 +486,51 @@ def test_fp_fit_transform_plot(X, hom_dims): def test_lsp_fit_transform(): - lp = LowerStarFlagPersistence() + lp = LowerStarFlagPersistence(extended=True) result = lp.fit_transform([X_lsp_cp])[0] assert_almost_equal(np.sort(result, axis=0), np.sort(diag_lsp_cp, axis=0)) + +@composite +def get_lsp_matrix(draw): + """Generate a 1d array of floats, of a given shape. If the shape is not + given, generate a shape of at least (4,).""" + n_points = draw(integers(3, 10)) + diag = draw(arrays(dtype=np.float, + elements=floats(allow_nan=False, + allow_infinity=False, + min_value=1, + max_value=1e2), + shape=(n_points,), unique=True)) + n_edges = draw(integers(2, int(n_points*(n_points-1)/2))) + list_vertices = lists(integers(min_value=0, max_value=n_points), + min_size=n_edges, max_size=n_edges) + edges = draw(lists(list_vertices, min_size=2, max_size=2, + unique_by=tuple(lambda x: x[k] for k in range(n_edges)))) + + edges = np.array(edges) + X = coo_matrix((np.concatenate([diag, np.ones(n_edges)]), + (np.concatenate([np.arange(n_points), edges[0]]), + np.concatenate([np.arange(n_points), edges[1]])))) + return X.toarray() + + +@given(X=get_lsp_matrix()) +def test_lsp_transform_duality(X): + X = coo_matrix(X) + hom_dims = (0, 1) + lp = LowerStarFlagPersistence(homology_dimensions=hom_dims, extended=True) + result = lp.fit_transform([X])[0] + result = result[np.where(np.logical_not( + np.logical_and(np.isclose(result[:, 0] - result[:, 1], 0), result[:, 3] == 1.)) + )] + same_sweep = result[np.where(result[:, 3] == 1)] + max_dim = max(hom_dims) + for dim in hom_dims: + dual_dim = max_dim - dim + primary, dual = [same_sweep[np.where(same_sweep[:, 2] == d)] + for d in [dim, dual_dim]] + assert_almost_equal(np.sort(primary[:, [0, 1]], axis=0), + np.sort(dual[:, [1, 0]], axis=0)) + From dab9c953cf82c787c543d4b8e8228e1337585964 Mon Sep 17 00:00:00 2001 From: Umberto Lupo Date: Fri, 12 Feb 2021 11:36:33 +0100 Subject: [PATCH 06/15] Rename variable after @wreise's suggestion --- gtda/homology/_utils.py | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/gtda/homology/_utils.py b/gtda/homology/_utils.py index f3259f1e4..e6ebdbf07 100644 --- a/gtda/homology/_utils.py +++ b/gtda/homology/_utils.py @@ -4,9 +4,8 @@ import numpy as np -def _postprocess_diagrams( - Xt, format, homology_dimensions, infinity_values, reduced - ): +def _postprocess_diagrams(Xt, format, homology_dimensions, infinity_values, + reduced): # NOTE: `homology_dimensions` must be sorted in ascending order def replace_infinity_values(subdiagram): np.nan_to_num(subdiagram, posinf=infinity_values, copy=False) @@ -21,7 +20,7 @@ def replace_infinity_values(subdiagram): for dim in homology_dimensions} Xt = [{dim: replace_infinity_values(diagram[dim][slices[dim]]) for dim in homology_dimensions} for diagram in Xt] - feature_vect_len = 3 + n_diagram_cols = 3 elif format == "gudhi": # Input is list of list of [dim, (birth, death)] # In H0, remove one infinite bar placed at the beginning by GUDHI only # if `reduce` is True @@ -32,11 +31,11 @@ def replace_infinity_values(subdiagram): if pers_info[0] == dim]).reshape(-1, 2)[slices[dim]] ) for dim in homology_dimensions} for diagram in Xt] - feature_vect_len = 3 + n_diagram_cols = 3 elif format == "extended": # Input is list of list of subdiagrams Xt = [{dim: diagram[dim] for dim in homology_dimensions} for diagram in Xt] - feature_vect_len = 4 + n_diagram_cols = 4 else: raise ValueError(f"Unknown input format {format} for collection of " f"diagrams.") @@ -52,7 +51,7 @@ def replace_infinity_values(subdiagram): min_values = [min_value if min_value != np.inf else 0 for min_value in min_values] n_features = start_idx_per_dim[-1] - Xt_padded = np.empty((len(Xt), n_features, feature_vect_len), dtype=float) + Xt_padded = np.empty((len(Xt), n_features, n_diagram_cols), dtype=float) Xt_padded[:, :, 3:] = 1. # Only applies to extended persistence for i, dim in enumerate(homology_dimensions): From 1837c6f0d8912271773b487869415d3f96dd69c4 Mon Sep 17 00:00:00 2001 From: Umberto Lupo Date: Fri, 12 Feb 2021 11:37:08 +0100 Subject: [PATCH 07/15] Convert to float32 in _lower_star_extended_diagram to avoid numerical issues --- gtda/homology/simplicial.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/gtda/homology/simplicial.py b/gtda/homology/simplicial.py index 977d3fe92..71f996713 100644 --- a/gtda/homology/simplicial.py +++ b/gtda/homology/simplicial.py @@ -302,8 +302,8 @@ def transform(self, X, y=None): X = check_point_clouds(X, accept_sparse=True, distance_matrices=self._is_precomputed) - Xt = Parallel(n_jobs=self.n_jobs)( - delayed(self._ripser_diagram)(x) for x in X) + Xt = Parallel(n_jobs=self.n_jobs)(delayed(self._ripser_diagram)(x) + for x in X) Xt = _postprocess_diagrams( Xt, "ripser", self._homology_dimensions, self.infinity_values_, @@ -691,8 +691,8 @@ def transform(self, X, y=None): X = check_point_clouds(X, accept_sparse=True, force_all_finite=True, distance_matrices=self._is_precomputed) - Xt = Parallel(n_jobs=self.n_jobs)( - delayed(self._ripser_diagram)(x) for x in X) + Xt = Parallel(n_jobs=self.n_jobs)(delayed(self._ripser_diagram)(x) + for x in X) Xt = _postprocess_diagrams( Xt, "ripser", self._homology_dimensions, self.infinity_values_, @@ -962,8 +962,8 @@ def transform(self, X, y=None): X = check_point_clouds(X, accept_sparse=True, distance_matrices=self._is_precomputed) - Xt = Parallel(n_jobs=self.n_jobs)( - delayed(self._gudhi_diagram)(x) for x in X) + Xt = Parallel(n_jobs=self.n_jobs)(delayed(self._gudhi_diagram)(x) + for x in X) Xt = _postprocess_diagrams( Xt, "gudhi", self._homology_dimensions, self.infinity_values_, @@ -1213,8 +1213,8 @@ def transform(self, X, y=None): check_is_fitted(self) X = check_point_clouds(X) - Xt = Parallel(n_jobs=self.n_jobs)( - delayed(self._weak_alpha_diagram)(x) for x in X) + Xt = Parallel(n_jobs=self.n_jobs)(delayed(self._weak_alpha_diagram)(x) + for x in X) Xt = _postprocess_diagrams( Xt, "ripser", self._homology_dimensions, self.infinity_values_, @@ -1728,8 +1728,8 @@ def transform(self, X, y=None): check_is_fitted(self) X = check_point_clouds(X, accept_sparse=True, distance_matrices=True) - Xt = Parallel(n_jobs=self.n_jobs)( - delayed(self._flagser_diagram)(x) for x in X) + Xt = Parallel(n_jobs=self.n_jobs)(delayed(self._flagser_diagram)(x) + for x in X) Xt = _postprocess_diagrams( Xt, "flagser", self._homology_dimensions, self.infinity_values_, @@ -1804,7 +1804,7 @@ def _lower_star_extended_diagram(self, X): n_points = max(X.shape) n_points_cone = n_points + 1 n_diag = min(X.shape) - X = X.tocoo() + X = X.tocoo().astype(np.float32) # ripser needs single precision data_diag = np.zeros(n_points_cone, dtype=X.dtype) data_diag[:n_diag] = X.diagonal() @@ -1897,8 +1897,8 @@ def transform(self, X, y=None): check_is_fitted(self, "_is_fitted") X = check_point_clouds(X, accept_sparse=True, distance_matrices=True) - Xt = Parallel(n_jobs=self.n_jobs)( - delayed(self._diagram_computer)(x) for x in X) + Xt = Parallel(n_jobs=self.n_jobs)(delayed(self._diagram_computer)(x) + for x in X) Xt = _postprocess_diagrams( Xt, self._format, self._homology_dimensions, self.infinity_values, From 10c4926a60e748a8246bc975a69a6830abd40396 Mon Sep 17 00:00:00 2001 From: Umberto Lupo Date: Fri, 12 Feb 2021 12:09:53 +0100 Subject: [PATCH 08/15] Plot zero-persistence points when from different sweeps --- gtda/plotting/persistence_diagrams.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/gtda/plotting/persistence_diagrams.py b/gtda/plotting/persistence_diagrams.py index c860b87d0..7de80fbbb 100644 --- a/gtda/plotting/persistence_diagrams.py +++ b/gtda/plotting/persistence_diagrams.py @@ -177,7 +177,10 @@ def plot_extended_diagram(diagram, homology_dimensions=None, if homology_dimensions is None: homology_dimensions = np.unique(diagram[:, 2]) - diag = diagram[diagram[:, 0] != diagram[:, 1]] + nontrivial_mask = np.logical_not( + np.logical_and(diagram[:, 0] == diagram[:, 1], diagram[:, 3] == 1.) + ) + diag = diagram[nontrivial_mask] diag_no_dims = diag[:, :2] posinfinite_mask = np.isposinf(diag_no_dims) neginfinite_mask = np.isneginf(diag_no_dims) From a5672a67b7d2eb6c9398f877363c224875acce91 Mon Sep 17 00:00:00 2001 From: wreise Date: Sun, 14 Feb 2021 17:18:48 +0100 Subject: [PATCH 09/15] Simplify tests --- gtda/homology/simplicial.py | 6 +-- gtda/homology/tests/test_simplicial.py | 57 ++++++++++++++++++++------ 2 files changed, 47 insertions(+), 16 deletions(-) diff --git a/gtda/homology/simplicial.py b/gtda/homology/simplicial.py index 71f996713..a80e8dee9 100644 --- a/gtda/homology/simplicial.py +++ b/gtda/homology/simplicial.py @@ -1809,7 +1809,7 @@ def _lower_star_extended_diagram(self, X): data_diag = np.zeros(n_points_cone, dtype=X.dtype) data_diag[:n_diag] = X.diagonal() max_value = data_diag[:n_diag].max() - min_value = data_diag[:n_diag].min() + min_value = data_diag[:n_points].min() data_diag[-1] = min_value - 1 off_diag = X.row != X.col @@ -1819,7 +1819,7 @@ def _lower_star_extended_diagram(self, X): col = np.concatenate([col_off_diag, np.full(n_points, n_points)]) data = np.concatenate([ np.maximum(data_diag[row_off_diag], data_diag[col_off_diag]), - 2 * max_value + 1 - data_diag[:n_points] + np.sign(max_value)*2 * max_value + 1 - data_diag[:n_points] ]) X = coo_matrix((data, (row, col)), @@ -1837,7 +1837,7 @@ def _lower_star_extended_diagram(self, X): np.logical_xor.reduce(mask_down_sweep, axis=1, keepdims=True)).astype(int) - 1 Xdgms[i][mask_down_sweep] = \ - 2 * max_value + 1 - Xdgms[i][mask_down_sweep] + np.sign(max_value)*2 * max_value + 1 - Xdgms[i][mask_down_sweep] if not i: Xdgms[i] = np.hstack([Xdgms[i][:-1, :], sgn[:-1, :]]) else: diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py index 6cf845d49..ad9839467 100644 --- a/gtda/homology/tests/test_simplicial.py +++ b/gtda/homology/tests/test_simplicial.py @@ -6,7 +6,7 @@ import pytest from hypothesis import given from hypothesis.extra.numpy import arrays, array_shapes -from hypothesis.strategies import composite, integers, floats, booleans, sets, tuples, lists +from hypothesis.strategies import composite, integers, floats, booleans, sets, tuples, lists, permutations from numpy.testing import assert_almost_equal from scipy.sparse import csr_matrix, coo_matrix from scipy.spatial.distance import pdist, squareform @@ -497,11 +497,8 @@ def get_lsp_matrix(draw): """Generate a 1d array of floats, of a given shape. If the shape is not given, generate a shape of at least (4,).""" n_points = draw(integers(3, 10)) - diag = draw(arrays(dtype=np.float, - elements=floats(allow_nan=False, - allow_infinity=False, - min_value=1, - max_value=1e2), + diag = draw(arrays(dtype=np.float32, + elements=integers(min_value=1, max_value=int(1e2)), shape=(n_points,), unique=True)) n_edges = draw(integers(2, int(n_points*(n_points-1)/2))) list_vertices = lists(integers(min_value=0, max_value=n_points), @@ -516,21 +513,55 @@ def get_lsp_matrix(draw): return X.toarray() -@given(X=get_lsp_matrix()) +@composite +def get_circle_matrix(draw): + n_points = draw(integers(3, 10)) + diag = draw(arrays(dtype=np.float32, + elements=integers(min_value=1, max_value=int(1e2)), + shape=(n_points,), unique=True)) + n_edges = n_points + rows = np.arange(n_points) + cols = draw(permutations(rows)) + + X = coo_matrix((np.concatenate([diag, np.ones(n_edges)]), + (np.concatenate([np.arange(n_points), rows]), + np.concatenate([np.arange(n_points), cols])))) + return X.toarray() + + +@given(X=get_circle_matrix()) +def test_lsp_transform_symmetry(X): + X = coo_matrix(X) + hom_dims = (0, 1) + lp = LowerStarFlagPersistence(homology_dimensions=hom_dims, extended=True) + result, result_m = lp.fit_transform([X])[0], lp.fit_transform([-X])[0] + result, result_m = [r[np.where(np.logical_not( + np.logical_and(np.isclose(r[:, 0] - r[:, 1], 0), r[:, 3] == 1.)))] + for r in [result, result_m]] + same_sweep, same_sweep_m = [r[np.where(r[:, 3] == 1)] + for r in [result, result_m]] + max_dim = max(hom_dims) + for dim in range(max_dim): + dual_dim = max_dim - dim - 1 + primary, dual = [s[np.where(s[:, 2] == d)] + for s, d in zip([same_sweep, same_sweep_m], [dim, dual_dim])] + assert_almost_equal(np.sort(primary[:, [0, 1]], axis=0), + np.sort(- dual[:, [1, 0]], axis=0), decimal=3) + + +@given(X=get_circle_matrix()) def test_lsp_transform_duality(X): X = coo_matrix(X) hom_dims = (0, 1) lp = LowerStarFlagPersistence(homology_dimensions=hom_dims, extended=True) result = lp.fit_transform([X])[0] - result = result[np.where(np.logical_not( - np.logical_and(np.isclose(result[:, 0] - result[:, 1], 0), result[:, 3] == 1.)) - )] + result = result[np.where(np.logical_not(np.logical_and( + np.isclose(result[:, 0] - result[:, 1], 0), result[:, 3] == 1.)))] same_sweep = result[np.where(result[:, 3] == 1)] max_dim = max(hom_dims) - for dim in hom_dims: + for dim in range(max_dim): dual_dim = max_dim - dim primary, dual = [same_sweep[np.where(same_sweep[:, 2] == d)] for d in [dim, dual_dim]] assert_almost_equal(np.sort(primary[:, [0, 1]], axis=0), - np.sort(dual[:, [1, 0]], axis=0)) - + np.sort(dual[:, [1, 0]], axis=0), decimal=3) From 7a3a517b31ab058913e31eacdc0e77bf10fa0cfb Mon Sep 17 00:00:00 2001 From: wreise Date: Sun, 14 Feb 2021 17:20:44 +0100 Subject: [PATCH 10/15] Negative values are possible --- gtda/homology/tests/test_simplicial.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py index ad9839467..f0892db88 100644 --- a/gtda/homology/tests/test_simplicial.py +++ b/gtda/homology/tests/test_simplicial.py @@ -519,11 +519,12 @@ def get_circle_matrix(draw): diag = draw(arrays(dtype=np.float32, elements=integers(min_value=1, max_value=int(1e2)), shape=(n_points,), unique=True)) + sign = 2*int(draw(booleans())) - 1 n_edges = n_points rows = np.arange(n_points) cols = draw(permutations(rows)) - X = coo_matrix((np.concatenate([diag, np.ones(n_edges)]), + X = coo_matrix((np.concatenate([sign*diag, np.ones(n_edges)]), (np.concatenate([np.arange(n_points), rows]), np.concatenate([np.arange(n_points), cols])))) return X.toarray() From 672d23b39cd7e640e5509030d6ac9329e620f63c Mon Sep 17 00:00:00 2001 From: wreise Date: Mon, 15 Feb 2021 08:34:19 +0100 Subject: [PATCH 11/15] Apply @ulupo's comments --- gtda/homology/simplicial.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/gtda/homology/simplicial.py b/gtda/homology/simplicial.py index a80e8dee9..7f3e78e14 100644 --- a/gtda/homology/simplicial.py +++ b/gtda/homology/simplicial.py @@ -1808,7 +1808,7 @@ def _lower_star_extended_diagram(self, X): data_diag = np.zeros(n_points_cone, dtype=X.dtype) data_diag[:n_diag] = X.diagonal() - max_value = data_diag[:n_diag].max() + max_value = data_diag[:n_points].max() min_value = data_diag[:n_points].min() data_diag[-1] = min_value - 1 @@ -1819,7 +1819,8 @@ def _lower_star_extended_diagram(self, X): col = np.concatenate([col_off_diag, np.full(n_points, n_points)]) data = np.concatenate([ np.maximum(data_diag[row_off_diag], data_diag[col_off_diag]), - np.sign(max_value)*2 * max_value + 1 - data_diag[:n_points] + min_value + 2 * (max_value - min_value) + 1 + - (data_diag[:n_points]-min_value) ]) X = coo_matrix((data, (row, col)), @@ -1837,7 +1838,8 @@ def _lower_star_extended_diagram(self, X): np.logical_xor.reduce(mask_down_sweep, axis=1, keepdims=True)).astype(int) - 1 Xdgms[i][mask_down_sweep] = \ - np.sign(max_value)*2 * max_value + 1 - Xdgms[i][mask_down_sweep] + min_value + 2 * (max_value - min_value) + 1\ + - (Xdgms[i][mask_down_sweep]-min_value) if not i: Xdgms[i] = np.hstack([Xdgms[i][:-1, :], sgn[:-1, :]]) else: From ed605699fd41ab356164e7dae6d8393b2b5fb077 Mon Sep 17 00:00:00 2001 From: wreise Date: Mon, 15 Feb 2021 08:34:32 +0100 Subject: [PATCH 12/15] Linting --- gtda/homology/tests/test_simplicial.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py index f0892db88..514c2aaaa 100644 --- a/gtda/homology/tests/test_simplicial.py +++ b/gtda/homology/tests/test_simplicial.py @@ -5,8 +5,9 @@ import plotly.io as pio import pytest from hypothesis import given -from hypothesis.extra.numpy import arrays, array_shapes -from hypothesis.strategies import composite, integers, floats, booleans, sets, tuples, lists, permutations +from hypothesis.extra.numpy import arrays +from hypothesis.strategies import composite, integers, booleans,\ + lists, permutations from numpy.testing import assert_almost_equal from scipy.sparse import csr_matrix, coo_matrix from scipy.spatial.distance import pdist, squareform @@ -504,7 +505,8 @@ def get_lsp_matrix(draw): list_vertices = lists(integers(min_value=0, max_value=n_points), min_size=n_edges, max_size=n_edges) edges = draw(lists(list_vertices, min_size=2, max_size=2, - unique_by=tuple(lambda x: x[k] for k in range(n_edges)))) + unique_by=tuple(lambda x: x[k] + for k in range(n_edges)))) edges = np.array(edges) X = coo_matrix((np.concatenate([diag, np.ones(n_edges)]), @@ -515,7 +517,7 @@ def get_lsp_matrix(draw): @composite def get_circle_matrix(draw): - n_points = draw(integers(3, 10)) + n_points = draw(integers(3, 50)) diag = draw(arrays(dtype=np.float32, elements=integers(min_value=1, max_value=int(1e2)), shape=(n_points,), unique=True)) @@ -545,7 +547,8 @@ def test_lsp_transform_symmetry(X): for dim in range(max_dim): dual_dim = max_dim - dim - 1 primary, dual = [s[np.where(s[:, 2] == d)] - for s, d in zip([same_sweep, same_sweep_m], [dim, dual_dim])] + for s, d in zip([same_sweep, same_sweep_m], + [dim, dual_dim])] assert_almost_equal(np.sort(primary[:, [0, 1]], axis=0), np.sort(- dual[:, [1, 0]], axis=0), decimal=3) From 6d7d77d0e8f0bd08e3a34e0e8ae9bcc6828bec02 Mon Sep 17 00:00:00 2001 From: wreise Date: Mon, 15 Feb 2021 21:24:19 +0100 Subject: [PATCH 13/15] Add equivariance test --- gtda/homology/tests/test_simplicial.py | 32 +++++++++++++++++++------- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py index 514c2aaaa..a028a764f 100644 --- a/gtda/homology/tests/test_simplicial.py +++ b/gtda/homology/tests/test_simplicial.py @@ -6,8 +6,8 @@ import pytest from hypothesis import given from hypothesis.extra.numpy import arrays -from hypothesis.strategies import composite, integers, booleans,\ - lists, permutations +from hypothesis.strategies import composite, integers, booleans, \ + lists, permutations, floats from numpy.testing import assert_almost_equal from scipy.sparse import csr_matrix, coo_matrix from scipy.spatial.distance import pdist, squareform @@ -495,8 +495,7 @@ def test_lsp_fit_transform(): @composite def get_lsp_matrix(draw): - """Generate a 1d array of floats, of a given shape. If the shape is not - given, generate a shape of at least (4,).""" + """Generate a matrix""" n_points = draw(integers(3, 10)) diag = draw(arrays(dtype=np.float32, elements=integers(min_value=1, max_value=int(1e2)), @@ -509,14 +508,31 @@ def get_lsp_matrix(draw): for k in range(n_edges)))) edges = np.array(edges) + c = draw(floats(min_value=-10, max_value=10, allow_nan=False, + allow_infinity=False, width=32)) X = coo_matrix((np.concatenate([diag, np.ones(n_edges)]), (np.concatenate([np.arange(n_points), edges[0]]), np.concatenate([np.arange(n_points), edges[1]])))) - return X.toarray() + X2 = coo_matrix((np.concatenate([diag + c, np.ones(n_edges)]), + (np.concatenate([np.arange(n_points), edges[0]]), + np.concatenate([np.arange(n_points), edges[1]])))) + return X.toarray(), X2.toarray(), c + + +@given(Xs=get_lsp_matrix()) +def test_lsp_transform_equivariance(Xs): + X, X2, c = Xs + X, X2 = coo_matrix(X), coo_matrix(X2) + hom_dims = (0, 1) + lp = LowerStarFlagPersistence(homology_dimensions=hom_dims, extended=True) + result, result_m = lp.fit_transform([X, X2]) + result_m[:, 0:2] -= c + assert_almost_equal(np.sort(result, axis=0), + np.sort(result_m, axis=0), decimal=3) @composite -def get_circle_matrix(draw): +def get_sparse_matrix(draw): n_points = draw(integers(3, 50)) diag = draw(arrays(dtype=np.float32, elements=integers(min_value=1, max_value=int(1e2)), @@ -532,7 +548,7 @@ def get_circle_matrix(draw): return X.toarray() -@given(X=get_circle_matrix()) +@given(X=get_sparse_matrix()) def test_lsp_transform_symmetry(X): X = coo_matrix(X) hom_dims = (0, 1) @@ -553,7 +569,7 @@ def test_lsp_transform_symmetry(X): np.sort(- dual[:, [1, 0]], axis=0), decimal=3) -@given(X=get_circle_matrix()) +@given(X=get_sparse_matrix()) def test_lsp_transform_duality(X): X = coo_matrix(X) hom_dims = (0, 1) From b5308a76b3d6aed65cb286b426cc894e25d917c0 Mon Sep 17 00:00:00 2001 From: wreise Date: Tue, 16 Feb 2021 08:53:24 +0100 Subject: [PATCH 14/15] Make sure that there are no implicit zeros --- gtda/homology/tests/test_simplicial.py | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py index a028a764f..681525613 100644 --- a/gtda/homology/tests/test_simplicial.py +++ b/gtda/homology/tests/test_simplicial.py @@ -501,7 +501,7 @@ def get_lsp_matrix(draw): elements=integers(min_value=1, max_value=int(1e2)), shape=(n_points,), unique=True)) n_edges = draw(integers(2, int(n_points*(n_points-1)/2))) - list_vertices = lists(integers(min_value=0, max_value=n_points), + list_vertices = lists(integers(min_value=0, max_value=n_points-1), min_size=n_edges, max_size=n_edges) edges = draw(lists(list_vertices, min_size=2, max_size=2, unique_by=tuple(lambda x: x[k] @@ -519,13 +519,21 @@ def get_lsp_matrix(draw): return X.toarray(), X2.toarray(), c +def filter_true_zero_persistence(r): + """Filter out the rows of a persistence diagram, where the points are + truly of 0 persistence.""" + return r[np.where(np.logical_not( + np.logical_and(np.isclose(r[:, 0] - r[:, 1], 0), r[:, 3] == 1.)))] + + @given(Xs=get_lsp_matrix()) def test_lsp_transform_equivariance(Xs): X, X2, c = Xs X, X2 = coo_matrix(X), coo_matrix(X2) hom_dims = (0, 1) lp = LowerStarFlagPersistence(homology_dimensions=hom_dims, extended=True) - result, result_m = lp.fit_transform([X, X2]) + result, result_m = [filter_true_zero_persistence(r) + for r in lp.fit_transform([X, X2])] result_m[:, 0:2] -= c assert_almost_equal(np.sort(result, axis=0), np.sort(result_m, axis=0), decimal=3) @@ -554,8 +562,7 @@ def test_lsp_transform_symmetry(X): hom_dims = (0, 1) lp = LowerStarFlagPersistence(homology_dimensions=hom_dims, extended=True) result, result_m = lp.fit_transform([X])[0], lp.fit_transform([-X])[0] - result, result_m = [r[np.where(np.logical_not( - np.logical_and(np.isclose(r[:, 0] - r[:, 1], 0), r[:, 3] == 1.)))] + result, result_m = [filter_true_zero_persistence(r) for r in [result, result_m]] same_sweep, same_sweep_m = [r[np.where(r[:, 3] == 1)] for r in [result, result_m]] @@ -575,8 +582,7 @@ def test_lsp_transform_duality(X): hom_dims = (0, 1) lp = LowerStarFlagPersistence(homology_dimensions=hom_dims, extended=True) result = lp.fit_transform([X])[0] - result = result[np.where(np.logical_not(np.logical_and( - np.isclose(result[:, 0] - result[:, 1], 0), result[:, 3] == 1.)))] + result = filter_true_zero_persistence(result) same_sweep = result[np.where(result[:, 3] == 1)] max_dim = max(hom_dims) for dim in range(max_dim): From bcc81b6e31dad10d4bdc01874ae8d024f6c26b30 Mon Sep 17 00:00:00 2001 From: wreise Date: Tue, 16 Feb 2021 18:04:50 +0100 Subject: [PATCH 15/15] Add test for lowerStarPersistence without extended homology --- gtda/homology/tests/test_simplicial.py | 27 +++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/gtda/homology/tests/test_simplicial.py b/gtda/homology/tests/test_simplicial.py index 681525613..6c4e20632 100644 --- a/gtda/homology/tests/test_simplicial.py +++ b/gtda/homology/tests/test_simplicial.py @@ -519,21 +519,30 @@ def get_lsp_matrix(draw): return X.toarray(), X2.toarray(), c -def filter_true_zero_persistence(r): +def filter_true_zero_persistence(r, is_extended=True): """Filter out the rows of a persistence diagram, where the points are truly of 0 persistence.""" - return r[np.where(np.logical_not( - np.logical_and(np.isclose(r[:, 0] - r[:, 1], 0), r[:, 3] == 1.)))] - - -@given(Xs=get_lsp_matrix()) -def test_lsp_transform_equivariance(Xs): + is_close = np.isclose(r[:, 0] - r[:, 1], 0) + if is_extended: + is_close_and_same_sweep = np.logical_and(np.isclose( + r[:, 0] - r[:, 1], 0), r[:, 3] == 1.) + mask = is_close_and_same_sweep + else: + mask = is_close + return r[np.where(np.logical_not(mask))] + + +@given(Xs=get_lsp_matrix(), is_extended=booleans()) +def test_lsp_transform_equivariance(Xs, is_extended): X, X2, c = Xs X, X2 = coo_matrix(X), coo_matrix(X2) hom_dims = (0, 1) - lp = LowerStarFlagPersistence(homology_dimensions=hom_dims, extended=True) - result, result_m = [filter_true_zero_persistence(r) + lp = LowerStarFlagPersistence(homology_dimensions=hom_dims, + extended=is_extended) + result, result_m = [filter_true_zero_persistence(r, + is_extended=is_extended) for r in lp.fit_transform([X, X2])] + result_m[:, 0:2] -= c assert_almost_equal(np.sort(result, axis=0), np.sort(result_m, axis=0), decimal=3)