Source code for cechmate.solver

import itertools
import time

import numpy as np
import phat


[docs]def phat_diagrams(simplices, show_inf=False, verbose=True): """ Compute the persistence diagram for :code:`simplices` using Phat. Parameters ----------- simplices: A list of lists of simplices and their distances the kth element is itself a list of tuples ([idx1, ..., idxk], dist) where [idx1, ..., idxk] is a list of vertices involved in the simplex and "dist" is the distance at which the simplex is added show_inf: Boolean Determines whether or not to return points that never die. Returns -------- dgms: list of diagrams the persistence diagram for Hk """ ## Convert simplices representation to sparse pivot column # -- sort by birth time, if tie, use order of simplex ordered_simplices = sorted(simplices, key=lambda x: (x[1], len(x[0]))) columns = _simplices_to_sparse_pivot_column(ordered_simplices, verbose) ## Setup boundary matrix and reduce if verbose: print("Computing persistence pairs...") tic = time.time() boundary_matrix = phat.boundary_matrix( columns=columns, representation=phat.representations.sparse_pivot_column ) pairs = boundary_matrix.compute_persistence_pairs() pairs.sort() if verbose: print( "Finished computing persistence pairs (Elapsed Time %.3g)" % (time.time() - tic) ) ## Setup persistence diagrams by reading off distances dgms = _process_distances(pairs, ordered_simplices) ## Add all unpaired simplices as infinite points if show_inf: dgms = _add_unpaired(dgms, pairs, simplices) ## Convert to arrays: dgms = [np.array(dgm) for dgm in dgms.values()] return dgms
def _simplices_to_sparse_pivot_column(ordered_simplices, verbose=False): """ """ idx = 0 columns = [] idxs2order = {} if verbose: print("Constructing boundary matrix...") tic = time.time() for idxs, dist in ordered_simplices: k = len(idxs) idxs = sorted(idxs) idxs2order[tuple(idxs)] = idx idxs = np.array(idxs) if len(idxs) == 1: columns.append((0, [])) else: # Get all faces with k-1 vertices collist = [] for fidxs in itertools.combinations(range(k), k - 1): fidxs = np.array(list(fidxs)) fidxs = tuple(idxs[fidxs]) if not fidxs in idxs2order: raise Exception( "Error: Not a proper filtration: %s added before %s" % (idxs, fidxs) ) return None collist.append(idxs2order[fidxs]) collist = sorted(collist) columns.append((k - 1, collist)) idx += 1 if verbose: print( "Finished constructing boundary matrix (Elapsed Time %.3g)" % (time.time() - tic) ) return columns def _process_distances(pairs, ordered_simplices): """ Setup persistence diagrams by reading off distances """ dgms = {} posneg = np.zeros(len(ordered_simplices)) for [bi, di] in pairs: bidxs, bd = ordered_simplices[bi] didxs, dd = ordered_simplices[di] assert posneg[bi] == 0 and posneg[di] == 0 posneg[bi], posneg[di] = 1, -1 assert dd >= bd # assert len(bidxs) == len(didxs) - 1 p = len(bidxs) - 1 # Don't add zero persistence pairs if bd != dd: dgms.setdefault(p, []).append([bd, dd]) return dgms def _add_unpaired(dgms, pairs, simplices): posneg = np.zeros(len(simplices)) for [bi, di] in pairs: assert posneg[bi] == 0 assert posneg[di] == 0 posneg[bi] = 1 posneg[di] = -1 for i in range(len(posneg)): if posneg[i] == 0: (idxs, dist) = simplices[i] p = len(idxs) - 1 if not p in dgms: dgms[p] = [] dgms[p].append([dist, np.inf]) return dgms