7105318798

Committed 05 Dec 2023 07:00PM UTC coverage: 83.945% (-0.1%) from 84.072%

Build # 7105318798

Build Type

push

github

Committed by

peekxc

Commit Message

Updates

Run Details

3 of 4 new or added lines in 2 files covered. (75.0%)

1 existing line in 1 file now uncovered.

800 of 953 relevant lines covered (83.95%)

2.52 hits per line

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

96.92

/src/splex/sparse.py

import numpy as np 
from numpy.typing import ArrayLike
from scipy.sparse import coo_array, spmatrix
from collections.abc import Sized
from array import array
from .generics import *
from .predicates import *
from .Simplex import Simplex

from more_itertools import flatten, collapse, chunked, unique_everseen 

# def iterable2sequence(S: Iterable[Hashable], dtype): # dtype = (np.uint16, 3)
#   """Converts an arbitrary iterable of homogenous types to Sequence"""
#   assert dtype is not None
#   from hirola import HashTable
#   S_arr = np.fromiter(S, dtype=dtype)
#   h = HashTable(len(S_arr)*1.15, dtype=dtype)
#   h.add(S_arr)

#   I,J in combinations(S_arr.T, S_arr.shape[1]-1)
#   S_arr[]
# qr = np.array(rank_combs(S), dtype=np.uint64)

def _fast_boundary(S: Iterable[SimplexConvertible], F: Iterable[SimplexConvertible], dtype) -> spmatrix:
  assert len(dtype) == 2
  from hirola import HashTable
  S_arr = np.fromiter(S, dtype=dtype) if dtype[1] > 1 else np.fromiter(collapse(S), dtype=dtype[0])
  if S_arr.ndim > 1:
    S_arr.sort(axis=1) ## Ensure lexicographical order
  if dtype[1] == 1:
    return coo_array((0, len(S_arr)), dtype=dtype[0])
  else:
    m, d = S_arr.shape
    F_arr = np.fromiter(F, dtype=(S_arr.dtype, d-1)) if d > 2 else np.fromiter(collapse(F), dtype=S_arr.dtype)
    if m == 0:
      return coo_array((F_arr.shape[0], 0), dtype=dtype[0])
    ## Use euler characteristic bound
    tbl_sz = max(max((3*m)*1.25, 3), F_arr.shape[0]*1.25)
    h = HashTable(tbl_sz, dtype=(S_arr.dtype, d-1)) if d > 2 else HashTable(tbl_sz, dtype=S_arr.dtype)
    h.add(F_arr)
    ind = np.arange(d).astype(int)
    I = np.empty(m*d, dtype=S_arr.dtype)
    J = np.empty(m*d, dtype=S_arr.dtype)
    X = np.empty(m*d, dtype=int)
    ## Column-wise assignnment
    for i, idx in enumerate(combinations(ind, len(ind)-1)):
      f_ind = h[S_arr[:,idx]] if d > 2 else np.ravel(h[S_arr[:,idx]]) ## get rows in (0,1)
      I[i*m:(i+1)*m] = f_ind                                          ## row indices
      J[i*m:(i+1)*m] = np.fromiter(range(m), dtype=S_arr.dtype)       ## col indices
      X[i*m:(i+1)*m] = np.repeat((-1)**i, m)                          ## always choose lex order (?)
    # X = np.fromiter(flatten([(-1)**np.argsort(I[J == j]) for j in range(m)]), dtype=int)
    # X = np.empty(m*d, dtype=int)
    # for j in range(m):
    #   X[J==j] = (-1)**np.argsort(I[J == j]) ## this is wrong 
    D = coo_array((X, (I,J)), shape=(len(h.keys), m))
    return D


def _boundary(S: Iterable[SimplexConvertible], F: Optional[Sequence[SimplexConvertible]] = None):
  ## Load faces. If not given, by definition, the given p-simplices contain their boundary faces.
  if F is None: 
    assert is_repeatable(S), "Simplex iterable must be repeatable (a generator is not sufficient!)"
    F = list(faces(S))
  
  ## Ensure faces 'F' is indexable
  assert isinstance(F, Sequence), "Faces must be a valid Sequence (supporting .index(*) with SimplexLike objects!)"

  ## Build the boundary matrix from the sequence
  m = 0
  I,J,X = [],[],[] # row, col, data 
  face_dict = { s : i for i,s in enumerate(F) }
  for (j,s) in enumerate(S):
    if dim(s) > 0:
      # I.extend([F.index(f) for f in faces(s, dim(s)-1)]) ## TODO: this is painfully slow
      I.extend([face_dict[f] for f in faces(s, dim(s)-1)])
      J.extend(repeat(j, dim(s)+1))
      X.extend(islice(cycle([1,-1]), dim(s)+1))
    m += 1
  D = coo_array((X, (I,J)), shape=(len(F), m)).tolil(copy=False)
  return D 

## TODO: investigate whether to make a 'ChainLike' for extension with 'BoundaryChain'?
## Or maybe its fine to just rely on something with __index__
## class Chain(SupportsIndex[int], [Coefficient], Protocol):
##    def __init__(orientation = ...)
##    def __iter__() -> tuple[int, Coefficient]:
##    def __add__(Chain) -> Chain
##    def __iadd__(Chain) -> None
##    def pivot() -> (int, coefficient) [for homology]
##    def index(SimplexConvertible) -> int 
## something like chain(s: SimplexLike, c: ComplexLike, oriented: bool) -> ndarray, or generator (index, value)
## complexLike could be overloaded to handle .index(), checked if Sequence[SimplexLike] or if just Container[int] (+implying)
def boundary_matrix(K: Union[ComplexLike, FiltrationLike], p: Optional[Union[int, tuple]] = None):
  """
  Constructs a sparse boundary matrix of a given simplicial object _K_

  Parameters: 
    K: simplicial complex (optionally filtered) or ComplexLike. 
    p: dimension of the p-chains to form the columns. 
  
  Returns: 
    D: sparse matrix representing either the full or p-th boundary matrix (as List-of-Lists format)
  """
  if isinstance(p, tuple):
    return (boundary_matrix(K, pi) for pi in p)
  else: 
    assert p is None or isinstance(p, Integral) and p >= 0, "p must be non-negative integer, or None"
    assert isinstance(K, ComplexLike) or isinstance(K, FiltrationLike), f"Unknown input type '{type(K)}'"
    if p is None:
      simplices = list(faces(K)) # to ensure repeatable
      D = _boundary(simplices)
    else:
      p_simplices = map(Simplex, faces(K, p=p))
      p_faces = list(map(Simplex, faces(K, p=p-1)))
      # D = _boundary(p_simplices, p_faces)
      D = _fast_boundary(p_simplices, p_faces, dtype=(np.uint16, p+1))
      # face_gen = list(unique_everseen(chunked(collapse([s.boundary() for s in faces(K, p)]), p)))
      # face_ranks = rank_combs(face_gen, n=card(K,0), order="lex")
      # if p == 1:
      #   face_gen = collapse(face_gen)
      #   p_faces = np.fromiter(face_gen, dtype=np.uint16)
      #   p_faces = p_faces[np.argsort(face_ranks)] # lex order
      # elif p > 1: 
      #   p_faces = np.fromiter(face_gen, dtype=(np.uint16, p))
      #   p_faces = p_faces[np.argsort(face_ranks)] # lex order
      # else: 
      #   p_faces = face_ranks
      # D = _fast_boundary(faces(K, p=p), p_faces, dtype=(np.uint16, p+1))
      # K_shape = card(K,p-1), card(K,p)
      # if D.shape != K_shape:
      #   D = D.reshape(K_shape) # handles degenerate cases, otherwise should error
    return D

## It's enough to 
def coboundary_matrix():
  pass 
  # D1.shape[1]-np.flip(D1.col)-1, D1.shape[0]-np.flip(D1.row)-1

1	import numpy as np	3✔
2	from numpy.typing import ArrayLike	3✔
3	from scipy.sparse import coo_array, spmatrix	3✔
4	from collections.abc import Sized	3✔
5	from array import array	3✔
6	from .generics import *	3✔
7	from .predicates import *	3✔
8	from .Simplex import Simplex	3✔
9
10	from more_itertools import flatten, collapse, chunked, unique_everseen	3✔
11
12	# def iterable2sequence(S: Iterable[Hashable], dtype): # dtype = (np.uint16, 3)
13	# """Converts an arbitrary iterable of homogenous types to Sequence"""
14	# assert dtype is not None
15	# from hirola import HashTable
16	# S_arr = np.fromiter(S, dtype=dtype)
17	# h = HashTable(len(S_arr)*1.15, dtype=dtype)
18	# h.add(S_arr)
19
20	# I,J in combinations(S_arr.T, S_arr.shape[1]-1)
21	# S_arr[]
22	# qr = np.array(rank_combs(S), dtype=np.uint64)
23
24	def _fast_boundary(S: Iterable[SimplexConvertible], F: Iterable[SimplexConvertible], dtype) -> spmatrix:	3✔
25	assert len(dtype) == 2	3✔
26	from hirola import HashTable	3✔
27	S_arr = np.fromiter(S, dtype=dtype) if dtype[1] > 1 else np.fromiter(collapse(S), dtype=dtype[0])	3✔
28	if S_arr.ndim > 1:	3✔
29	S_arr.sort(axis=1) ## Ensure lexicographical order	3✔
30	if dtype[1] == 1:	3✔
31	return coo_array((0, len(S_arr)), dtype=dtype[0])	3✔
32	else:
33	m, d = S_arr.shape	3✔
34	F_arr = np.fromiter(F, dtype=(S_arr.dtype, d-1)) if d > 2 else np.fromiter(collapse(F), dtype=S_arr.dtype)	3✔
35	if m == 0:	3✔
36	return coo_array((F_arr.shape[0], 0), dtype=dtype[0])	3✔
37	## Use euler characteristic bound
38	tbl_sz = max(max((3m)1.25, 3), F_arr.shape[0]*1.25)	3✔
39	h = HashTable(tbl_sz, dtype=(S_arr.dtype, d-1)) if d > 2 else HashTable(tbl_sz, dtype=S_arr.dtype)	3✔
40	h.add(F_arr)	3✔
41	ind = np.arange(d).astype(int)	3✔
42	I = np.empty(m*d, dtype=S_arr.dtype)	3✔
43	J = np.empty(m*d, dtype=S_arr.dtype)	3✔
44	X = np.empty(m*d, dtype=int)	3✔
45	## Column-wise assignnment
46	for i, idx in enumerate(combinations(ind, len(ind)-1)):	3✔
47	f_ind = h[S_arr[:,idx]] if d > 2 else np.ravel(h[S_arr[:,idx]]) ## get rows in (0,1)	3✔
48	I[im:(i+1)m] = f_ind ## row indices	3✔
49	J[im:(i+1)m] = np.fromiter(range(m), dtype=S_arr.dtype) ## col indices	3✔
50	X[im:(i+1)m] = np.repeat((-1)**i, m) ## always choose lex order (?)	3✔
51	# X = np.fromiter(flatten([(-1)**np.argsort(I[J == j]) for j in range(m)]), dtype=int)
52	# X = np.empty(m*d, dtype=int)
53	# for j in range(m):
54	# X[J==j] = (-1)**np.argsort(I[J == j]) ## this is wrong
55	D = coo_array((X, (I,J)), shape=(len(h.keys), m))	3✔
56	return D	3✔
57
58
59	def _boundary(S: Iterable[SimplexConvertible], F: Optional[Sequence[SimplexConvertible]] = None):	3✔
60	## Load faces. If not given, by definition, the given p-simplices contain their boundary faces.
61	if F is None:	3✔
62	assert is_repeatable(S), "Simplex iterable must be repeatable (a generator is not sufficient!)"	3✔
63	F = list(faces(S))	3✔
64
65	## Ensure faces 'F' is indexable
66	assert isinstance(F, Sequence), "Faces must be a valid Sequence (supporting .index(*) with SimplexLike objects!)"	3✔
67
68	## Build the boundary matrix from the sequence
69	m = 0	3✔
70	I,J,X = [],[],[] # row, col, data	3✔
71	face_dict = { s : i for i,s in enumerate(F) }	3✔
72	for (j,s) in enumerate(S):	3✔
73	if dim(s) > 0:	3✔
74	# I.extend([F.index(f) for f in faces(s, dim(s)-1)]) ## TODO: this is painfully slow
75	I.extend([face_dict[f] for f in faces(s, dim(s)-1)])	3✔
76	J.extend(repeat(j, dim(s)+1))	3✔
77	X.extend(islice(cycle([1,-1]), dim(s)+1))	3✔
78	m += 1	3✔
79	D = coo_array((X, (I,J)), shape=(len(F), m)).tolil(copy=False)	3✔
80	return D	3✔
81
82	## TODO: investigate whether to make a 'ChainLike' for extension with 'BoundaryChain'?
83	## Or maybe its fine to just rely on something with __index__
84	## class Chain(SupportsIndex[int], [Coefficient], Protocol):
85	## def __init__(orientation = ...)
86	## def __iter__() -> tuple[int, Coefficient]:
87	## def __add__(Chain) -> Chain
88	## def __iadd__(Chain) -> None
89	## def pivot() -> (int, coefficient) [for homology]
90	## def index(SimplexConvertible) -> int
91	## something like chain(s: SimplexLike, c: ComplexLike, oriented: bool) -> ndarray, or generator (index, value)
92	## complexLike could be overloaded to handle .index(), checked if Sequence[SimplexLike] or if just Container[int] (+implying)
93	def boundary_matrix(K: Union[ComplexLike, FiltrationLike], p: Optional[Union[int, tuple]] = None):	3✔
94	"""
95	Constructs a sparse boundary matrix of a given simplicial object _K_
96
97	Parameters:
98	K: simplicial complex (optionally filtered) or ComplexLike.
99	p: dimension of the p-chains to form the columns.
100
101	Returns:
102	D: sparse matrix representing either the full or p-th boundary matrix (as List-of-Lists format)
103	"""
104	if isinstance(p, tuple):	3✔
105	return (boundary_matrix(K, pi) for pi in p)	×
106	else:
107	assert p is None or isinstance(p, Integral) and p >= 0, "p must be non-negative integer, or None"	3✔
108	assert isinstance(K, ComplexLike) or isinstance(K, FiltrationLike), f"Unknown input type '{type(K)}'"	3✔
109	if p is None:	3✔
110	simplices = list(faces(K)) # to ensure repeatable	3✔
111	D = _boundary(simplices)	3✔
112	else:
113	p_simplices = map(Simplex, faces(K, p=p))	3✔
114	p_faces = list(map(Simplex, faces(K, p=p-1)))	3✔
115	# D = _boundary(p_simplices, p_faces)
116	D = _fast_boundary(p_simplices, p_faces, dtype=(np.uint16, p+1))	3✔
117	# face_gen = list(unique_everseen(chunked(collapse([s.boundary() for s in faces(K, p)]), p)))
118	# face_ranks = rank_combs(face_gen, n=card(K,0), order="lex")
119	# if p == 1:
120	# face_gen = collapse(face_gen)
121	# p_faces = np.fromiter(face_gen, dtype=np.uint16)
122	# p_faces = p_faces[np.argsort(face_ranks)] # lex order
123	# elif p > 1:
124	# p_faces = np.fromiter(face_gen, dtype=(np.uint16, p))
125	# p_faces = p_faces[np.argsort(face_ranks)] # lex order
126	# else:
127	# p_faces = face_ranks
128	# D = _fast_boundary(faces(K, p=p), p_faces, dtype=(np.uint16, p+1))
129	# K_shape = card(K,p-1), card(K,p)
130	# if D.shape != K_shape:
131	# D = D.reshape(K_shape) # handles degenerate cases, otherwise should error
132	return D	3✔
133
134	## It's enough to
135	def coboundary_matrix():	3✔
UNCOV 136	pass	×
137	# D1.shape[1]-np.flip(D1.col)-1, D1.shape[0]-np.flip(D1.row)-1

peekxc / splex / 7105318798

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