6566946791

Committed 18 Oct 2023 09:47PM UTC coverage: 82.528% (-2.2%) from 84.734%

Build # 6566946791

Build Type

push

github

Committed by

peekxc

Commit Message

Working on the abc for filtration types. I think I have a good idea of what a good wrapper type should implement. Started the Filtration abc, which going to be similar to a Sequence type mixed with a Set type, but have some additional constraints

Run Details

58 of 58 new or added lines in 5 files covered. (100.0%)

666 of 807 relevant lines covered (82.53%)

2.48 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.24

/src/splex/geometry.py

import numpy as np
from scipy.spatial.distance import pdist, squareform

from .generics import *
from .complexes import * 
from .filtrations import *
from .predicates import *
from combin import rank_to_comb, comb_to_rank, inverse_choose
from dataclassy import dataclass

def as_pairwise_dist(x: ArrayLike) -> ArrayLike:
  if is_point_cloud(x):
    pd = pdist(x)
  elif is_dist_like(x):
    pd = np.tril(x) if is_distance_matrix(x) else x
  else: 
    raise ValueError("Unknown input shape 'x' ")
  return pd

def enclosing_radius(x: ArrayLike) -> float:
  """Returns the smallest 'r' such that the Rips complex on the union of balls of radius 'r' is contractible to a point. """
  if is_point_cloud(x):
    d = squareform(pdist(x))
    return 0.5*np.min(np.amax(d, axis = 0))
  elif is_dist_like(x):
    assert is_distance_matrix(x) or is_pairwise_distances(x), "Must be valid distances"
    d = x if is_distance_matrix(x) else squareform(x)
    return 0.5*np.min(np.amax(d, axis = 0))
  else:
    raise ValueError("Unknown input type")

def rips_complex(x: ArrayLike, radius: float = None, p: int = 1) -> FiltrationLike:
  """Constructs the Vietoris-Rips complex from _x_ by unioning balls of diameter at most 2 * _radius_ 
  
  Parameters: 
    x: point cloud, pairwise distance vector, or distance matrix
    radius: scale parameter for the Rips complex. 
    p: highest dimension of simplices to consider in the expansion. 
  
  Returns: 
    rips complex, returned as a simplex tree
  """ 
  from simplextree import SimplexTree
  pd = as_pairwise_dist(x)
  n = inverse_choose(len(pd), 2)
  radius = enclosing_radius(squareform(pd)) if radius is None else float(radius)
  ind = np.flatnonzero(pd <= 2*radius)
  edges = rank_to_comb(ind, n=n, k=2, order="lex")
  st = SimplexTree([[i] for i in range(n)])
  st.insert(edges)
  st.expand(p)
  return st

def flag_weight(x: ArrayLike, vertex_weights: Optional[ArrayLike] = None) -> Callable:
  """Filter function factory method for constructing flag/clique filter functions. 

  Parameters: 
    x: point cloud, vector of pairwise weights, or square matrix. 
    vertex_weights: optional weights to use for vertices. Defaults to None, which sets vertex weights to 0.

  Returns: 
    callable which takes as input a simplex or set of simplices and returns their clique weights. 
  """
  pd = as_pairwise_dist(x)
  n = inverse_choose(len(pd), 2)
  vertex_weights = np.zeros(n) if vertex_weights is None else vertex_weights
  assert len(vertex_weights) == n, "Invalid vertex weights"
  @dataclass(frozen=True, slots=True, init=False, repr=False, eq=False)
  class _clique_weight:
    n: int = 0
    vertex_weights: np.ndarray = np.empty(0, dtype=np.float32)
    edge_weights: np.ndarray = np.empty(0, dtype=np.float32)
    def __init__(self, v: np.ndarray, pd: np.ndarray, n: int) -> None:
      object.__setattr__(self, 'n', n)
      object.__setattr__(self, 'vertex_weights', v)
      object.__setattr__(self, 'edge_weights', pd)
    def __call__(self, s: Union[SimplexConvertible, ArrayLike]) -> Union[float, np.ndarray]:
      if hasattr(s, "__array__") and is_complex_like(s):
        ## Handles numpy matrices of simplices OR array_convertible containers, so long as they are complex-like
        s = np.asarray(s)
        if s.ndim == 1 or (1 in s.shape):
          # print("hello")
          return np.ravel(self.vertex_weights[s])
        else: 
          if s.shape[1] == 2: 
            ind = comb_to_rank(s, n=self.n, order='lex')
            return self.edge_weights[ind]
          else:
            fw = np.zeros(s.shape[0])
            for i,j in combinations(range(s.shape[1]), 2):
              np.maximum(fw, self.edge_weights[comb_to_rank(s[:,[i,j]], n=self.n, order='lex')], out=fw)
            return fw
      elif is_simplex_like(s):
        if len(s) == 1: 
          return self.vertex_weights[s] 
        else: 
          ind = np.array(comb_to_rank(combinations(s, 2), n=self.n, order='lex'), dtype=np.uint64)
          # ind = np.fromiter((rank_lex(e, n=self.n) for e in combinations(s, 2)), dtype=np.uint32)
          return np.max(self.edge_weights[ind])
      else:
        # assert is_complex_like(s), "Input simplices must be complex like" # this is unneeded since not not be sized or repeatable
        rank_boundary = lambda f: np.array(comb_to_rank(combinations(f, 2), n=self.n, order='lex'), dtype=np.uint64)
        return np.array([np.max(self.edge_weights[rank_boundary(f)]) if dim(f) >= 1 else np.take(self.vertex_weights[f],0) for f in s], dtype=float)
  C = _clique_weight(vertex_weights, pd, n)
  return C

## TODO: revamp to include support for arbitrary simplicial complexes via index tracking with hirola 
def lower_star_weight(x: ArrayLike) -> Callable:
  """Constructs a simplex-parameterized _Callable_ that evaluates its lower star value based on _x_. 

  Vertex labels are assumed to be 0-indexed for now. 
  
  If simplex-like, use 0-indexed vertex labels to index vertex values directly. 
  
  Otherwise assumes a 2d array of simplex labels is given and vectorizes the computation.
  """
  @dataclass(frozen=True, slots=True, init=False, repr=False, eq=False)
  class LS:
    vertex_weights: np.ndarray = np.empty(0, dtype=np.float32)

    def __init__(self, v: np.ndarray) -> None:
      object.__setattr__(self, 'vertex_weights', v)
    
    def __call__(self, S: Union[SimplexConvertible, ArrayLike]) -> Union[float, np.ndarray]:
      if is_simplex_like(S):
        return np.max(self.vertex_weights[Simplex(S)]) # Simplices can be used for indexing!
      elif hasattr(S, "__array__") and is_complex_like(S):
        S = np.asarray(S)
        return np.max(self.vertex_weights[S], axis=-1) ## vectorized form
      else:
        assert isinstance(S, Iterable), "simplex iterable must be supplied"
        return np.array([np.max(self.vertex_weights[s]) for s in map(Simplex, S)])
  return LS(x)

def rips_filtration(x: ArrayLike, radius: float = None, p: int = 1, **kwargs) -> FiltrationLike:
  """Constructs a _p_-dimensional rips filtration from _x_ by unioning balls of diameter at most 2 * _radius_
  """
  from simplextree import SimplexTree
  pd = as_pairwise_dist(x)
  radius = enclosing_radius(pd) if radius is None else float(radius)
  ind = np.flatnonzero(pd <= 2*radius)
  n = inverse_choose(len(pd), 2)
  st = SimplexTree([[i] for i in range(n)])
  st.insert(rank_to_comb(ind, n=n, k=2, order="lex"))
  st.expand(p)
  f = flag_weight(pd)
  G = ((f([s]), s) for s in st)
  K = filtration(G, **kwargs)
  return K

def delaunay_complex(x: ArrayLike):
  from scipy.spatial import Delaunay
  dt = Delaunay(x)
  T = dt.simplices
  V = np.fromiter(range(x.shape[0]), dtype=np.int32)
  S = simplicial_complex(chain(V, T))
  return(S)



1	import numpy as np	3✔
2	from scipy.spatial.distance import pdist, squareform	3✔
3
4	from .generics import *	3✔
5	from .complexes import *	3✔
6	from .filtrations import *	3✔
7	from .predicates import *	3✔
8	from combin import rank_to_comb, comb_to_rank, inverse_choose	3✔
9	from dataclassy import dataclass	3✔
10
11	def as_pairwise_dist(x: ArrayLike) -> ArrayLike:	3✔
12	if is_point_cloud(x):	3✔
13	pd = pdist(x)	3✔
14	elif is_dist_like(x):	3✔
15	pd = np.tril(x) if is_distance_matrix(x) else x	3✔
16	else:
17	raise ValueError("Unknown input shape 'x' ")	×
18	return pd	3✔
19
20	def enclosing_radius(x: ArrayLike) -> float:	3✔
21	"""Returns the smallest 'r' such that the Rips complex on the union of balls of radius 'r' is contractible to a point. """
22	if is_point_cloud(x):	3✔
23	d = squareform(pdist(x))	3✔
24	return 0.5*np.min(np.amax(d, axis = 0))	3✔
25	elif is_dist_like(x):	3✔
26	assert is_distance_matrix(x) or is_pairwise_distances(x), "Must be valid distances"	3✔
27	d = x if is_distance_matrix(x) else squareform(x)	3✔
28	return 0.5*np.min(np.amax(d, axis = 0))	3✔
29	else:
30	raise ValueError("Unknown input type")	×
31
32	def rips_complex(x: ArrayLike, radius: float = None, p: int = 1) -> FiltrationLike:	3✔
33	"""Constructs the Vietoris-Rips complex from _x_ by unioning balls of diameter at most 2 * _radius_
34
35	Parameters:
36	x: point cloud, pairwise distance vector, or distance matrix
37	radius: scale parameter for the Rips complex.
38	p: highest dimension of simplices to consider in the expansion.
39
40	Returns:
41	rips complex, returned as a simplex tree
42	"""
43	from simplextree import SimplexTree	3✔
44	pd = as_pairwise_dist(x)	3✔
45	n = inverse_choose(len(pd), 2)	3✔
46	radius = enclosing_radius(squareform(pd)) if radius is None else float(radius)	3✔
47	ind = np.flatnonzero(pd <= 2*radius)	3✔
48	edges = rank_to_comb(ind, n=n, k=2, order="lex")	3✔
49	st = SimplexTree([[i] for i in range(n)])	3✔
50	st.insert(edges)	3✔
51	st.expand(p)	3✔
52	return st	3✔
53
54	def flag_weight(x: ArrayLike, vertex_weights: Optional[ArrayLike] = None) -> Callable:	3✔
55	"""Filter function factory method for constructing flag/clique filter functions.
56
57	Parameters:
58	x: point cloud, vector of pairwise weights, or square matrix.
59	vertex_weights: optional weights to use for vertices. Defaults to None, which sets vertex weights to 0.
60
61	Returns:
62	callable which takes as input a simplex or set of simplices and returns their clique weights.
63	"""
64	pd = as_pairwise_dist(x)	3✔
65	n = inverse_choose(len(pd), 2)	3✔
66	vertex_weights = np.zeros(n) if vertex_weights is None else vertex_weights	3✔
67	assert len(vertex_weights) == n, "Invalid vertex weights"	3✔
68	@dataclass(frozen=True, slots=True, init=False, repr=False, eq=False)	3✔
69	class _clique_weight:	3✔
70	n: int = 0	3✔
71	vertex_weights: np.ndarray = np.empty(0, dtype=np.float32)	3✔
72	edge_weights: np.ndarray = np.empty(0, dtype=np.float32)	3✔
73	def __init__(self, v: np.ndarray, pd: np.ndarray, n: int) -> None:	3✔
74	object.__setattr__(self, 'n', n)	3✔
75	object.__setattr__(self, 'vertex_weights', v)	3✔
76	object.__setattr__(self, 'edge_weights', pd)	3✔
77	def __call__(self, s: Union[SimplexConvertible, ArrayLike]) -> Union[float, np.ndarray]:	3✔
78	if hasattr(s, "__array__") and is_complex_like(s):	3✔
79	## Handles numpy matrices of simplices OR array_convertible containers, so long as they are complex-like
80	s = np.asarray(s)	3✔
81	if s.ndim == 1 or (1 in s.shape):	3✔
82	# print("hello")
83	return np.ravel(self.vertex_weights[s])	×
84	else:
85	if s.shape[1] == 2:	3✔
86	ind = comb_to_rank(s, n=self.n, order='lex')	3✔
87	return self.edge_weights[ind]	3✔
88	else:
89	fw = np.zeros(s.shape[0])	3✔
90	for i,j in combinations(range(s.shape[1]), 2):	3✔
91	np.maximum(fw, self.edge_weights[comb_to_rank(s[:,[i,j]], n=self.n, order='lex')], out=fw)	3✔
92	return fw	3✔
93	elif is_simplex_like(s):	3✔
94	if len(s) == 1:	3✔
95	return self.vertex_weights[s]	3✔
96	else:
97	ind = np.array(comb_to_rank(combinations(s, 2), n=self.n, order='lex'), dtype=np.uint64)	3✔
98	# ind = np.fromiter((rank_lex(e, n=self.n) for e in combinations(s, 2)), dtype=np.uint32)
99	return np.max(self.edge_weights[ind])	3✔
100	else:
101	# assert is_complex_like(s), "Input simplices must be complex like" # this is unneeded since not not be sized or repeatable
102	rank_boundary = lambda f: np.array(comb_to_rank(combinations(f, 2), n=self.n, order='lex'), dtype=np.uint64)	3✔
103	return np.array([np.max(self.edge_weights[rank_boundary(f)]) if dim(f) >= 1 else np.take(self.vertex_weights[f],0) for f in s], dtype=float)	3✔
104	C = _clique_weight(vertex_weights, pd, n)	3✔
105	return C	3✔
106
107	## TODO: revamp to include support for arbitrary simplicial complexes via index tracking with hirola
108	def lower_star_weight(x: ArrayLike) -> Callable:	3✔
109	"""Constructs a simplex-parameterized _Callable_ that evaluates its lower star value based on _x_.
110
111	Vertex labels are assumed to be 0-indexed for now.
112
113	If simplex-like, use 0-indexed vertex labels to index vertex values directly.
114
115	Otherwise assumes a 2d array of simplex labels is given and vectorizes the computation.
116	"""
117	@dataclass(frozen=True, slots=True, init=False, repr=False, eq=False)	3✔
118	class LS:	3✔
119	vertex_weights: np.ndarray = np.empty(0, dtype=np.float32)	3✔
120
121	def __init__(self, v: np.ndarray) -> None:	3✔
122	object.__setattr__(self, 'vertex_weights', v)	3✔
123
124	def __call__(self, S: Union[SimplexConvertible, ArrayLike]) -> Union[float, np.ndarray]:	3✔
125	if is_simplex_like(S):	3✔
126	return np.max(self.vertex_weights[Simplex(S)]) # Simplices can be used for indexing!	3✔
127	elif hasattr(S, "__array__") and is_complex_like(S):	3✔
128	S = np.asarray(S)	3✔
129	return np.max(self.vertex_weights[S], axis=-1) ## vectorized form	3✔
130	else:
131	assert isinstance(S, Iterable), "simplex iterable must be supplied"	×
132	return np.array([np.max(self.vertex_weights[s]) for s in map(Simplex, S)])	×
133	return LS(x)	3✔
134
135	def rips_filtration(x: ArrayLike, radius: float = None, p: int = 1, **kwargs) -> FiltrationLike:	3✔
136	"""Constructs a _p_-dimensional rips filtration from _x_ by unioning balls of diameter at most 2 * _radius_
137	"""
138	from simplextree import SimplexTree	3✔
139	pd = as_pairwise_dist(x)	3✔
140	radius = enclosing_radius(pd) if radius is None else float(radius)	3✔
141	ind = np.flatnonzero(pd <= 2*radius)	3✔
142	n = inverse_choose(len(pd), 2)	3✔
143	st = SimplexTree([[i] for i in range(n)])	3✔
144	st.insert(rank_to_comb(ind, n=n, k=2, order="lex"))	3✔
145	st.expand(p)	3✔
146	f = flag_weight(pd)	3✔
147	G = ((f([s]), s) for s in st)	3✔
148	K = filtration(G, **kwargs)	3✔
149	return K	3✔
150
151	def delaunay_complex(x: ArrayLike):	3✔
152	from scipy.spatial import Delaunay	3✔
153	dt = Delaunay(x)	3✔
154	T = dt.simplices	3✔
155	V = np.fromiter(range(x.shape[0]), dtype=np.int32)	3✔
156	S = simplicial_complex(chain(V, T))	3✔
157	return(S)	3✔
158
159

peekxc / splex / 6566946791

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous