|
| 1 | +import scipy |
| 2 | +from toponetx import CellComplex |
| 3 | +from toponetx.utils import incidence_to_adjacency |
| 4 | + |
| 5 | + |
| 6 | +def cc_incidence_matrix( |
| 7 | + cell_complex: CellComplex, |
| 8 | + rank: int, |
| 9 | + signed: bool = True, |
| 10 | + index: bool = False, |
| 11 | +) -> scipy.sparse.csr_matrix | tuple[dict, dict, scipy.sparse.csr_matrix]: |
| 12 | + """ |
| 13 | + Same function as toponet but returns the boundary matrix of the cell complex in the correct order. |
| 14 | + """ |
| 15 | + nodelist = ( |
| 16 | + cell_complex.nodes |
| 17 | + ) # order simplices as they appear in the cell complex |
| 18 | + if rank == 0: |
| 19 | + A = scipy.sparse.lil_matrix((0, len(nodelist))) |
| 20 | + if index: |
| 21 | + node_index = {node: i for i, node in enumerate(nodelist)} |
| 22 | + if signed: |
| 23 | + return {}, node_index, A.asformat("csr") |
| 24 | + return {}, node_index, abs(A.asformat("csr")) |
| 25 | + |
| 26 | + if signed: |
| 27 | + return A.asformat("csr") |
| 28 | + return abs(A.asformat("csr")) |
| 29 | + # For rank 1 and 2, we need to have the indices of the vertices. |
| 30 | + node_index = {node: i for i, node in enumerate(nodelist)} |
| 31 | + # edgelist contains edges composed by the indices of the vertices they contain sorted, |
| 32 | + # this is, with the induced orientation by the order of the vertices in the cell complex. |
| 33 | + edgelist = [ |
| 34 | + sorted((node_index[e[0]], node_index[e[1]])) |
| 35 | + for e in cell_complex.edges |
| 36 | + ] |
| 37 | + if rank == 1: |
| 38 | + A = scipy.sparse.lil_matrix((len(nodelist), len(edgelist))) |
| 39 | + for ei, e in enumerate(edgelist): |
| 40 | + (ui, vi) = e[ |
| 41 | + :2 |
| 42 | + ] # Note that the indices are sorted, so we are orienting the edges |
| 43 | + # by the order of the nodes in the cell complex given by their indices |
| 44 | + A[ui, ei] = -1 |
| 45 | + A[vi, ei] = 1 |
| 46 | + if index: |
| 47 | + edge_index = {edge: i for i, edge in enumerate(edgelist)} |
| 48 | + if signed: |
| 49 | + return node_index, edge_index, A.asformat("csr") |
| 50 | + return node_index, edge_index, abs(A.asformat("csr")) |
| 51 | + if signed: |
| 52 | + return A.asformat("csr") |
| 53 | + return abs(A.asformat("csr")) |
| 54 | + if rank == 2: |
| 55 | + A = scipy.sparse.lil_matrix((len(edgelist), len(cell_complex.cells))) |
| 56 | + edge_index = { |
| 57 | + tuple(edge): i for i, edge in enumerate(edgelist) |
| 58 | + } # oriented edges indices |
| 59 | + for celli, cell in enumerate(cell_complex.cells): |
| 60 | + edge_visiting_dic = {} # this dictionary is cell dependent |
| 61 | + # mainly used to handle the cell complex non-regular case |
| 62 | + for edge in cell.boundary: |
| 63 | + edge_w_indices = (node_index[edge[0]], node_index[edge[1]]) |
| 64 | + ei = edge_index[tuple(sorted(edge_w_indices))] |
| 65 | + if ei not in edge_visiting_dic: |
| 66 | + if edge_w_indices in edge_index: |
| 67 | + edge_visiting_dic[ei] = 1 |
| 68 | + else: |
| 69 | + edge_visiting_dic[ei] = -1 |
| 70 | + else: |
| 71 | + if edge in edge_index: |
| 72 | + edge_visiting_dic[ei] = edge_visiting_dic[ei] + 1 |
| 73 | + else: |
| 74 | + edge_visiting_dic[ei] = edge_visiting_dic[ei] - 1 |
| 75 | + |
| 76 | + A[ei, celli] = edge_visiting_dic[ |
| 77 | + ei |
| 78 | + ] # this will update everytime we visit this edge for non-regular cell complexes |
| 79 | + # the regular case can be handled more efficiently : |
| 80 | + # if edge in edge_index: |
| 81 | + # A[ei, celli] = 1s |
| 82 | + # else: |
| 83 | + # A[ei, celli] = -1 |
| 84 | + if index: |
| 85 | + cell_index = { |
| 86 | + c.elements: i for i, c in enumerate(cell_complex.cells) |
| 87 | + } |
| 88 | + if signed: |
| 89 | + return edge_index, cell_index, A.asformat("csr") |
| 90 | + return edge_index, cell_index, abs(A.asformat("csr")) |
| 91 | + |
| 92 | + if signed: |
| 93 | + return A.asformat("csr") |
| 94 | + return abs(A.asformat("csr")) |
| 95 | + raise ValueError(f"Only dimensions 0, 1 and 2 are supported, got {rank}.") |
| 96 | + |
| 97 | + |
| 98 | +def hodge_laplacian_matrix( |
| 99 | + cell_complex: CellComplex, |
| 100 | + rank: int, |
| 101 | + signed: bool = True, |
| 102 | +) -> scipy.sparse.csr_matrix: |
| 103 | + assert ( |
| 104 | + cell_complex.dim >= rank >= 0 |
| 105 | + ) # No negative dimensional Hodge Laplacian |
| 106 | + if cell_complex.dim > rank >= 0: |
| 107 | + up_laplacian = up_laplacian_matrix(cell_complex, rank, True) |
| 108 | + else: |
| 109 | + up_laplacian = None |
| 110 | + if cell_complex.dim >= rank > 0: |
| 111 | + down_laplacian = down_laplacian_matrix(cell_complex, rank, True) |
| 112 | + else: |
| 113 | + down_laplacian = None |
| 114 | + if up_laplacian is not None and down_laplacian is not None: |
| 115 | + hodge_laplacian = up_laplacian + down_laplacian |
| 116 | + elif up_laplacian is not None: |
| 117 | + hodge_laplacian = up_laplacian |
| 118 | + elif down_laplacian is not None: |
| 119 | + hodge_laplacian = down_laplacian |
| 120 | + else: |
| 121 | + # Dimension is 0 because, if the dimension of the cell complex is one or higher, |
| 122 | + # we have at least lower laplacians. Also, if the dimension is greater than zero, |
| 123 | + # we have at least upper laplacians for all dimensions except for cell_complex.dim, for |
| 124 | + # which we have lower laplacian. |
| 125 | + hodge_laplacian = scipy.sparse.coo_matrix( |
| 126 | + (len(cell_complex.nodes), len(cell_complex.nodes)) |
| 127 | + ) |
| 128 | + if not signed: |
| 129 | + hodge_laplacian = abs(hodge_laplacian) |
| 130 | + return hodge_laplacian |
| 131 | + |
| 132 | + |
| 133 | +def up_laplacian_matrix( |
| 134 | + cell_complex: CellComplex, |
| 135 | + rank: int, |
| 136 | + signed: bool = True, |
| 137 | +) -> scipy.sparse.csr_matrix: |
| 138 | + """ |
| 139 | + Same function as toponet but returns the upper laplacian of the cell complex in the correct order. |
| 140 | + """ |
| 141 | + |
| 142 | + if cell_complex.dim > rank >= 0: |
| 143 | + B_next = cc_incidence_matrix(cell_complex, rank + 1) |
| 144 | + L_up = B_next @ B_next.transpose() |
| 145 | + else: |
| 146 | + raise ValueError( |
| 147 | + f"Rank should be larger or equal than 0 and <= {cell_complex.dim - 1} (maximal dimension cells-1), got {rank}" |
| 148 | + ) |
| 149 | + if not signed: |
| 150 | + L_up = abs(L_up) |
| 151 | + return L_up |
| 152 | + |
| 153 | + |
| 154 | +def down_laplacian_matrix( |
| 155 | + cell_complex: CellComplex, rank: int, signed: bool = True, weight=None |
| 156 | +) -> scipy.sparse.csr_matrix: |
| 157 | + """ |
| 158 | + Same function as toponet but returns the lower laplacian of the cell complex in the correct order. |
| 159 | + """ |
| 160 | + if weight is not None: |
| 161 | + raise ValueError("`weight` is not supported in this version") |
| 162 | + |
| 163 | + if cell_complex.dim >= rank > 0: |
| 164 | + B = cc_incidence_matrix(cell_complex, rank) |
| 165 | + L_down = B.transpose() @ B |
| 166 | + else: |
| 167 | + raise ValueError( |
| 168 | + f"Rank should be larger or equal than 1 and <= {cell_complex.dim} (maximal dimension cells), got {rank}." |
| 169 | + ) |
| 170 | + if not signed: |
| 171 | + L_down = abs(L_down) |
| 172 | + return L_down |
| 173 | + |
| 174 | + |
| 175 | +def lower_adjacency( |
| 176 | + cell_complex: CellComplex, dim: int, s: int = 1 |
| 177 | +) -> scipy.sparse.spmatrix: |
| 178 | + # A cell is neighbor of itself and all the other cells appearing in the lower hodge laplacian. |
| 179 | + if dim == 0: |
| 180 | + return scipy.sparse.coo_matrix( |
| 181 | + (len(cell_complex.nodes), len(cell_complex.nodes)) |
| 182 | + ) |
| 183 | + else: |
| 184 | + B = cc_incidence_matrix(cell_complex, dim, signed=False) |
| 185 | + A = incidence_to_adjacency(B, s=s) |
| 186 | + return A.tocoo() |
| 187 | + |
| 188 | + |
| 189 | +def upper_adjacency( |
| 190 | + cell_complex: CellComplex, dim: int, s: int = 1 |
| 191 | +) -> scipy.sparse.spmatrix: |
| 192 | + if cell_complex.dim == dim: |
| 193 | + match dim: |
| 194 | + case 0: |
| 195 | + return scipy.sparse.coo_matrix( |
| 196 | + (len(cell_complex.nodes), len(cell_complex.nodes)) |
| 197 | + ) |
| 198 | + case 1: |
| 199 | + return scipy.sparse.coo_matrix( |
| 200 | + (len(cell_complex.edges), len(cell_complex.edges)) |
| 201 | + ) |
| 202 | + case 2: |
| 203 | + return scipy.sparse.coo_matrix( |
| 204 | + (len(cell_complex.cells), len(cell_complex.cells)) |
| 205 | + ) |
| 206 | + else: |
| 207 | + # A cell is neighbor of itself and all the other cells appearing in the upper hodge laplacian. |
| 208 | + B_T = cc_incidence_matrix( |
| 209 | + cell_complex, dim + 1, signed=False |
| 210 | + ).transpose() |
| 211 | + A = incidence_to_adjacency(B_T, s=s) |
| 212 | + return A.tocoo() |
0 commit comments