optim: moved mobius inversion in point_measure

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):

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(

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
 

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

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
 

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