# MIT License
# Copyright (c) 2024 Niklas Hellmer
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
# This cell is copy-pasted from pyDowker
# by Niklas Hellmer and Jan Spalinski, but only the portion
# that's relevant to our project.
# See https://github.com/nihell/pyDowker
from gudhi import SimplexTree
import numpy as np
[docs]
class DowkerComplex:
"""
Class MNeighborComplex. Constructs Dowker's simplicial complex for a relation.
Filtrations can be added using filtered relations, or total weight, or combining both into a bifiltration.
"""
[docs]
def __init__(self, rel_matrix, max_filtration=float('inf')) -> None:
"""_summary_
Constructor for the filtered Dowker complex from the relations given by sublevels of the matrix rel_matrix.
The vertices in the complex will correspond to the rows of the matrix.
Parameters
----------
rel_matrix (Sequence[Sequence[float]]): distance matrix (full square or lower triangular).
max_filtration (float): specifies the maximal filtration value to be considered.
"""
self.rel_matrix = rel_matrix
self.st = None
def create_simplex_tree(self,
max_dimension,
filtration = 'None',
m=1,
level = 0,
max_filtration = np.inf):
"""
Creates a gudhi simplex tree storing a skeleton of the (filtered) simpicial complex.
Uses recursive algorithm with nummpy arrays, fast for small datasets but worse runtime and memory complexity.
Parameters
----------
max_dimension : int
Dimension of the skeleton to compute.
filtration : str, optional
valid values: "Sublevel", "TotalWeight", "None".
"Sublevel" takes the filtration of relations by sublevels of the matrix.
"Total weight" takes the filtration by sublevels of the negative total weight function.
By default 'None'
m : int, optional
restriction to this superlevel of total weight (this is only used if filtration!="TotalWeight"); m=1 corresponds to the whole Dowker complex, by default 1
level : int, optional
restriction to this sublevel of the matrix as relation (this is only used if filtration!="Sublevel"), by default 0
max_filtration : float, optional
cutoff for the filtration (only used if filtration="Sublevel"), by default np.inf
Returns
-------
gudhi.SimplexTree
The simplex tree storing the (filtered) simplicial complex
"""
self.st=SimplexTree()
LAMBDA = self.rel_matrix
num_points=len(LAMBDA)
if filtration == "Sublevel":
if LAMBDA.dtype != np.float64:
raise TypeError("Only float arrays are allowed with sublevel filtration")
def append_upper_cofaces(sigma, r, witness_values):
if r > max_filtration:
return
self.st.insert(sigma,r)
if len(sigma)<=max_dimension:
for j in range(np.max(sigma)+1,num_points):
tau = sigma+[j]
j_witness_values=LAMBDA[j,:]
common_witness_values = np.maximum(j_witness_values,witness_values)
new_r = np.partition(common_witness_values, m-1)[m-1]
append_upper_cofaces(tau, new_r, common_witness_values)
for k in range(num_points-1,-1,-1):
witness_values = LAMBDA[k,:]
r_new = np.partition(witness_values, m-1)[m-1]
append_upper_cofaces([k],r_new,witness_values)
return self.st
elif filtration == "TotalWeight":
if LAMBDA.dtype != np.bool_:
LAMBDA = LAMBDA <= level
def append_upper_cofaces(sigma, witnesses):
self.st.insert(sigma,-np.sum(witnesses))
if len(sigma)<=max_dimension:
for j in range(np.max(sigma)+1,num_points):
tau = sigma+[j]
j_witnesses=LAMBDA[j,:]
common_witnesses = np.logical_and(j_witnesses,witnesses)
if np.sum(common_witnesses>0):
append_upper_cofaces(tau, common_witnesses)
for k in range(num_points-1,-1,-1):
witnesses = LAMBDA[k,:]
append_upper_cofaces([k], witnesses)
return self.st
elif filtration == "None":
if LAMBDA.dtype != np.bool_:
LAMBDA = LAMBDA <= level
def append_upper_cofaces(sigma, witnesses):
if len(witnesses)<m:
return
self.st.insert(sigma)
if len(sigma)<=max_dimension:
for j in range(np.max(sigma)+1,num_points):
tau = sigma+[j]
j_witnesses=LAMBDA[j,:]
common_witnesses = np.logical_and(j_witnesses,witnesses)
if len(common_witnesses>0):
append_upper_cofaces(tau, common_witnesses)
for k in range(num_points-1,-1,-1):
witnesses = LAMBDA[k,:]
append_upper_cofaces([k], witnesses)
return self.st
else:
raise Exception("filtration parameter must be one of 'Sublevel', 'TotalWeight', 'None'")