From 0dda850db94d774e80844cd99620f01b3400990a Mon Sep 17 00:00:00 2001 From: David Loiseaux Date: Wed, 8 Jan 2025 16:26:22 +0100 Subject: [PATCH] optim: moved mobius inversion in point_measure --- multipers/ml/signed_betti.py | 48 -------------------- multipers/ml/signed_measures.py | 1 + multipers/point_measure.pyx | 51 ++++++++++++++++++++++ multipers/simplex_tree_multi.pyx.tp | 3 +- multipers/slicer.pyx.tp | 3 +- multipers/tests/old_test_rank_invariant.py | 2 +- multipers/tests/test_signed_betti.py | 2 +- 7 files changed, 56 insertions(+), 54 deletions(-) delete mode 100644 multipers/ml/signed_betti.py diff --git a/multipers/ml/signed_betti.py b/multipers/ml/signed_betti.py deleted file mode 100644 index e3bfeac..0000000 --- a/multipers/ml/signed_betti.py +++ /dev/null @@ -1,48 +0,0 @@ -## This code was initially written by Luis Scoccola -import numpy as np -# from scipy.sparse import coo_array -# from scipy.ndimage import convolve1d - - -def signed_betti(hilbert_function:np.ndarray, threshold:bool=False): - n = hilbert_function.ndim - # hilbert_function = hilbert_function if inplace else np.copy(hilbert_function) - # zero out the "end" of the Hilbert function - if threshold: - for dimension in range(n): - slicer = tuple([slice(None) if i != dimension else -1 for i in range(n)]) - hilbert_function[slicer] = 0 - # weights = np.array([0, 1, -1], dtype=int) - for i in range(n): - # res = convolve1d(res, weights, axis=i, mode="constant", cval=0) - minus = tuple([slice(None) if j!=i else slice(0,-1) for j in range(n)]) - plus = tuple([slice(None) if j!=i else slice(1,None) for j in range(n)]) - hilbert_function[plus]-= hilbert_function[minus] - else: - return hilbert_function - -def rank_decomposition_by_rectangles(rank_invariant:np.ndarray, threshold:bool=False): - # takes as input the rank invariant of an n-parameter persistence module - # M : [0, ..., s_1 - 1] x ... x [0, ..., s_n - 1] ---> Vec - # on a grid with dimensions of sizes s_1, ..., s_n. The input is assumed to be - # given as a tensor of dimensions (s_1, ..., s_n, s_1, ..., s_n), so that, - # at index [i_1, ..., i_n, j_1, ..., j_n] we have the rank of the structure - # map M(i) -> M(j), where i = (i_1, ..., i_n) and j = (j_1, ..., j_n), and - # i <= j, meaning that i_1 <= j_1, ..., i_n <= j_n. - # NOTE : - # - About the input, we assume that, if not( i <= j ), then at index - # [i_1, ..., i_n, j_1, ..., j_n] we have a zero. - # - Similarly, the output at index [i_1, ..., i_n, j_1, ..., j_n] only - # makes sense when i <= j. For indices where not( i <= j ) the output - # may take arbitrary values and they should be ignored. - n = rank_invariant.ndim // 2 - if threshold: - # zero out the "end" - for dimension in range(n): - slicer = tuple( - [slice(None) for _ in range(n)] - + [slice(None) if i != dimension else -1 for i in range(n)] - ) - rank_invariant[slicer] = 0 - to_flip = tuple(range(n, 2 * n)) - return np.flip(signed_betti(np.flip(rank_invariant, to_flip)), to_flip) diff --git a/multipers/ml/signed_measures.py b/multipers/ml/signed_measures.py index 132074f..f56002d 100644 --- a/multipers/ml/signed_measures.py +++ b/multipers/ml/signed_measures.py @@ -11,6 +11,7 @@ import multipers as mp from multipers.grids import compute_grid as reduce_grid from multipers.ml.convolutions import available_kernels, convolution_signed_measures +from multipers.point_measure import signed_betti, rank_decomposition_by_rectangles class FilteredComplex2SignedMeasure(BaseEstimator, TransformerMixin): diff --git a/multipers/point_measure.pyx b/multipers/point_measure.pyx index f344dbf..f3b0d6b 100644 --- a/multipers/point_measure.pyx +++ b/multipers/point_measure.pyx @@ -29,6 +29,57 @@ ctypedef fused some_float: import cython cimport cython + + + +# from scipy.sparse import coo_array +# from scipy.ndimage import convolve1d + +def signed_betti(hilbert_function:np.ndarray, bool threshold=False): + cdef int n = hilbert_function.ndim + # zero out the "end" of the Hilbert function + if threshold: + for dimension in range(n): + slicer = tuple([slice(None) if i != dimension else -1 for i in range(n)]) + hilbert_function[slicer] = 0 + for i in range(n): + minus = tuple([slice(None) if j!=i else slice(0,-1) for j in range(n)]) + plus = tuple([slice(None) if j!=i else slice(1,None) for j in range(n)]) + hilbert_function[plus] -= hilbert_function[minus] + return hilbert_function + +def rank_decomposition_by_rectangles(rank_invariant:np.ndarray, bool threshold=False): + # takes as input the rank invariant of an n-parameter persistence module + # M : [0, ..., s_1 - 1] x ... x [0, ..., s_n - 1] ---> Vec + # on a grid with dimensions of sizes s_1, ..., s_n. The input is assumed to be + # given as a tensor of dimensions (s_1, ..., s_n, s_1, ..., s_n), so that, + # at index [i_1, ..., i_n, j_1, ..., j_n] we have the rank of the structure + # map M(i) -> M(j), where i = (i_1, ..., i_n) and j = (j_1, ..., j_n), and + # i <= j, meaning that i_1 <= j_1, ..., i_n <= j_n. + # NOTE : + # - About the input, we assume that, if not( i <= j ), then at index + # [i_1, ..., i_n, j_1, ..., j_n] we have a zero. + # - Similarly, the output at index [i_1, ..., i_n, j_1, ..., j_n] only + # makes sense when i <= j. For indices where not( i <= j ) the output + # may take arbitrary values and they should be ignored. + cdef int n = rank_invariant.ndim // 2 + if threshold: + # zero out the "end" + for dimension in range(n): + slicer = tuple( + [slice(None) for _ in range(n)] + + [slice(None) if i != dimension else -1 for i in range(n)] + ) + rank_invariant[slicer] = 0 + to_flip = tuple(range(n, 2 * n)) + return np.flip(signed_betti(np.flip(rank_invariant, to_flip)), to_flip) + + + + + + + @cython.boundscheck(False) @cython.wraparound(False) def integrate_measure( diff --git a/multipers/simplex_tree_multi.pyx.tp b/multipers/simplex_tree_multi.pyx.tp index 5b34b77..80a6ee4 100644 --- a/multipers/simplex_tree_multi.pyx.tp +++ b/multipers/simplex_tree_multi.pyx.tp @@ -68,8 +68,7 @@ from typing import Iterable,Literal,Optional from tqdm import tqdm import multipers.grids as mpg -from multipers.ml.signed_betti import rank_decomposition_by_rectangles -from multipers.point_measure import sparsify +from multipers.point_measure import signed_betti, rank_decomposition_by_rectangles, sparsify from warnings import warn diff --git a/multipers/slicer.pyx.tp b/multipers/slicer.pyx.tp index 260fa8f..dd899f9 100644 --- a/multipers/slicer.pyx.tp +++ b/multipers/slicer.pyx.tp @@ -26,8 +26,7 @@ from multipers.slicer cimport * from multipers.filtrations cimport * from multipers.filtration_conversions cimport * ## TODO: these two are not needed, remove that by updating rank code. -from multipers.point_measure import sparsify -from multipers.ml.signed_betti import rank_decomposition_by_rectangles +from multipers.point_measure import sparsify, rank_decomposition_by_rectangles import numpy as np cimport cython diff --git a/multipers/tests/old_test_rank_invariant.py b/multipers/tests/old_test_rank_invariant.py index d970cea..e1d0132 100644 --- a/multipers/tests/old_test_rank_invariant.py +++ b/multipers/tests/old_test_rank_invariant.py @@ -1,6 +1,6 @@ import multipers as mp import numpy as np -from multipers.ml.signed_betti import signed_betti +from multipers.point_measure import signed_betti from multipers.hilbert_function import hilbert_surface from multipers.euler_characteristic import euler_surface diff --git a/multipers/tests/test_signed_betti.py b/multipers/tests/test_signed_betti.py index 02bfcf6..df2fdde 100644 --- a/multipers/tests/test_signed_betti.py +++ b/multipers/tests/test_signed_betti.py @@ -1,5 +1,5 @@ import numpy as np -from multipers.ml.signed_betti import signed_betti, rank_decomposition_by_rectangles +from multipers.point_measure import signed_betti, rank_decomposition_by_rectangles # only tests rank functions with 1 and 2 parameters