From a62525c9bf8d9ff95393d52d244e9e6d7417a840 Mon Sep 17 00:00:00 2001 From: Jonathan Demaeyer Date: Wed, 10 Jun 2026 22:18:48 +0200 Subject: [PATCH 1/3] Added a block inversion algorithm for numerical case --- layercake/utils/matrix.py | 149 +++++++++++++++++++++++++++----------- 1 file changed, 108 insertions(+), 41 deletions(-) diff --git a/layercake/utils/matrix.py b/layercake/utils/matrix.py index db78e95..f5ab018 100644 --- a/layercake/utils/matrix.py +++ b/layercake/utils/matrix.py @@ -9,67 +9,105 @@ """ import warnings +import numpy as np +from numpy.linalg import LinAlgError + from sympy import MutableSparseMatrix from sympy.matrices.exceptions import NonInvertibleMatrixError def block_matrix_inverse(P, blocks_extent, simplify=True): - """Function to invert a symbolic matrix :math:`P` divided by blocks. + """Function to invert a symbolic or numeric matrix :math:`P` divided by blocks. Parameters ---------- - P: ~sympy.matrices.immutable.ImmutableSparseMatrix or ~sympy.matrices.mutable.MutableSparseMatrix + P: ~sympy.matrices.immutable.ImmutableSparseMatrix or ~sympy.matrices.mutable.MutableSparseMatrix or ~numpy.ndarray The block matrix to invert. blocks_extent: list(tuple) The extent of each block, as a list of 2-tuple. simplify: bool, optional Try to simplify the inverse. + Only applies to symbolic matrices. Default to `True`. Warnings -------- - If the `simplify` argument is set to `False`, this function will not check if the + In the symbolic case, if the `simplify` argument is set to `False`, this function will not check if the matrix :math:`P` is invertible ! """ - if simplify: - if P.det().simplify() == 0: - raise NonInvertibleMatrixError('block_matrix_inverse: The provided matrix is not invertible.') - else: - warnings.warn(f'Inverting a symbolic matrix without checking that it is invertible. ' - 'Be cautious about the result.') - be = blocks_extent.copy() - PP = P - B_list = list() - C_list = list() - Am1_list = list() - while len(be) >= 2: - rest_block_extent = (be[1][0], be[-1][1]) - Ainv, B, C, D, PPP = _block_matrix_inverse_2x2(PP, (be[0], rest_block_extent)) - - B_list.append(B) - C_list.append(C) - Am1_list.append(Ainv) - - PP = PPP - ret = be[0][1] - be = [(p[0]-ret, p[1]-ret) for p in be[1:]] - - PP = PP.adjugate() / PP.det() - for Am1, B, C in zip(Am1_list[::-1], B_list[::-1], C_list[::-1]): - Ashape = Am1.shape[0] - Pshape = PP.shape[0] - mat_shape = Ashape + Pshape - Pp1 = MutableSparseMatrix(mat_shape, mat_shape, {}) - Pp1[slice(0, Ashape), slice(0, Ashape)] = Am1 + Am1 @ B @ PP @ C @ Am1 - Pp1[slice(0, Ashape), slice(Ashape, mat_shape)] = - Am1 @ B @ PP - Pp1[slice(Ashape, mat_shape), slice(0, Ashape)] = - PP @ C @ Am1 - Pp1[slice(Ashape, mat_shape), slice(Ashape, mat_shape)] = PP - PP = Pp1 - - if simplify: - return PP.simplify() - else: + if isinstance(P, np.ndarray): + if np.linalg.det(P) == 0: + raise LinAlgError('block_matrix_inverse: The provided matrix is not invertible.') + be = blocks_extent.copy() + PP = P + B_list = list() + C_list = list() + Am1_list = list() + while len(be) >= 2: + rest_block_extent = (be[1][0], be[-1][1]) + Ainv, B, C, D, PPP = _block_matrix_inverse_2x2_num(PP, (be[0], rest_block_extent)) + + B_list.append(B) + C_list.append(C) + Am1_list.append(Ainv) + + PP = PPP + ret = be[0][1] + be = [(p[0] - ret, p[1] - ret) for p in be[1:]] + + PP = np.linalg.inv(PP) + for Am1, B, C in zip(Am1_list[::-1], B_list[::-1], C_list[::-1]): + Ashape = Am1.shape[0] + Pshape = PP.shape[0] + mat_shape = Ashape + Pshape + Pp1 = np.zeros((mat_shape, mat_shape)) + Pp1[slice(0, Ashape), slice(0, Ashape)] = Am1 + Am1 @ B @ PP @ C @ Am1 + Pp1[slice(0, Ashape), slice(Ashape, mat_shape)] = - Am1 @ B @ PP + Pp1[slice(Ashape, mat_shape), slice(0, Ashape)] = - PP @ C @ Am1 + Pp1[slice(Ashape, mat_shape), slice(Ashape, mat_shape)] = PP + PP = Pp1 + return PP + else: + if simplify: + if P.det().simplify() == 0: + raise NonInvertibleMatrixError('block_matrix_inverse: The provided matrix is not invertible.') + else: + warnings.warn(f'Inverting a symbolic matrix without checking that it is invertible. ' + 'Be cautious about the result.') + be = blocks_extent.copy() + PP = P + B_list = list() + C_list = list() + Am1_list = list() + while len(be) >= 2: + rest_block_extent = (be[1][0], be[-1][1]) + Ainv, B, C, D, PPP = _block_matrix_inverse_2x2(PP, (be[0], rest_block_extent)) + + B_list.append(B) + C_list.append(C) + Am1_list.append(Ainv) + + PP = PPP + ret = be[0][1] + be = [(p[0]-ret, p[1]-ret) for p in be[1:]] + + PP = PP.adjugate() / PP.det() + for Am1, B, C in zip(Am1_list[::-1], B_list[::-1], C_list[::-1]): + Ashape = Am1.shape[0] + Pshape = PP.shape[0] + mat_shape = Ashape + Pshape + Pp1 = MutableSparseMatrix(mat_shape, mat_shape, {}) + Pp1[slice(0, Ashape), slice(0, Ashape)] = Am1 + Am1 @ B @ PP @ C @ Am1 + Pp1[slice(0, Ashape), slice(Ashape, mat_shape)] = - Am1 @ B @ PP + Pp1[slice(Ashape, mat_shape), slice(0, Ashape)] = - PP @ C @ Am1 + Pp1[slice(Ashape, mat_shape), slice(Ashape, mat_shape)] = PP + PP = Pp1 + + if simplify: + return PP.simplify() + else: + return PP def _block_matrix_inverse_2x2(P, blocks_extent): be = blocks_extent @@ -81,6 +119,16 @@ def _block_matrix_inverse_2x2(P, blocks_extent): return Ainv, B, C, D, D - C @ Ainv @ B +def _block_matrix_inverse_2x2_num(P, blocks_extent): + be = blocks_extent + A = P[slice(*be[0]), slice(*be[0])] + B = P[slice(*be[0]), slice(*be[1])] + C = P[slice(*be[1]), slice(*be[0])] + D = P[slice(*be[1]), slice(*be[1])] + Ainv = np.linalg.inv(A) + + return Ainv, B, C, D, D - C @ Ainv @ B + if __name__ == '__main__': A = MutableSparseMatrix(3, 3, {}) @@ -102,3 +150,22 @@ def _block_matrix_inverse_2x2(P, blocks_extent): bl = [(0, 3), (3, 6), (6, 9)] Ginv = block_matrix_inverse(G, bl) + A = np.zeros((3, 3)) + A[0,0] = 1 + A[1,1] = 3 + A[2,2] = 2 + C = A /2 + D = 3 * A + G = np.zeros((9, 9)) + G[:3, :3] = A + G[:3, 3:6] = -C + G[:3, 6:] = D + G[3:6, :3] = C + G[3:6, 3:6] = C + G[3:6, 6:] = - D + G[6:, :3] = - C + G[6:, 3:6] = A + G[6:, 6:] = D + bl = [(0, 3), (3, 6), (6, 9)] + Ginv_num = block_matrix_inverse(G, bl) + From 2756307808ec5e1b56ce847e73ae699792d89c2c Mon Sep 17 00:00:00 2001 From: Jonathan Demaeyer Date: Thu, 11 Jun 2026 15:16:23 +0200 Subject: [PATCH 2/3] Matrix inversion algorithm is now configurable --- layercake/bakery/cake.py | 4 +++- layercake/bakery/layers.py | 5 +++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/layercake/bakery/cake.py b/layercake/bakery/cake.py index cc5f0f7..3e19a8b 100644 --- a/layercake/bakery/cake.py +++ b/layercake/bakery/cake.py @@ -55,6 +55,7 @@ def __init__(self): self.layers = list() self._lhs_inversion_in_layer = True self._simplify_after_LHS_inversion = True + self._matrix_inv = np.linalg.inv def add_layer(self, layer): """Add a layer object to the cake. @@ -65,6 +66,7 @@ def add_layer(self, layer): Layer object to add to the cake. """ layer._cake_order = len(self.layers) + layer._matrix_inv = self._matrix_inv self.layers.append(layer) layer._cake = self for equation in layer.equations: @@ -262,7 +264,7 @@ def tensor(self): if not self._lhs_inversion_in_layer: try: lhs_mat_inverted = np.zeros((self.ndim + 1, self.ndim + 1)) - lhs_mat_inverted[1:, 1:] = np.linalg.inv(lhs_mat.todense()[1:, 1:]) + lhs_mat_inverted[1:, 1:] = self._matrix_inv(lhs_mat.todense()[1:, 1:]) tensor = sp.COO(np.tensordot(lhs_mat_inverted, tensor.to_coo(), 1)) except LinAlgError: raise LinAlgError(f'The left-hand side of the cake is not invertible with the provided basis.') diff --git a/layercake/bakery/layers.py b/layercake/bakery/layers.py index f76a58f..0493d1b 100644 --- a/layercake/bakery/layers.py +++ b/layercake/bakery/layers.py @@ -63,6 +63,7 @@ def __init__(self, name=''): self._lhs_inverted = False self._lhs_mat = None self._simplify_after_LHS_inversion = True + self._matrix_inv = np.linalg.inv @property def _cake_first_index(self): @@ -298,7 +299,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i self._lhs_mat[lhs_order:lhs_order + ndim, ofield_order:ofield_order + ondim] + lhs_term.inner_products.todense() else: try: - lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = np.linalg.inv(eq.lhs_inner_products_addition.todense()) + lhs_mat_inverted[lhs_order:lhs_order + ndim, lhs_order:lhs_order + ndim] = self._matrix_inv(eq.lhs_inner_products_addition.todense()) self._lhs_inverted = True except LinAlgError: raise LinAlgError(f'The left-hand side of the equation {eq} is not invertible with the provided basis.') @@ -345,7 +346,7 @@ def compute_tensor(self, numerical=True, compute_inner_products=False, compute_i lhs_order += ndim if self._lhs_inversion and not self._lhs_inverted: try: - lhs_mat_inverted[1:, 1:] = np.linalg.inv(self._lhs_mat.todense()[1:, 1:]) + lhs_mat_inverted[1:, 1:] = self._matrix_inv(self._lhs_mat.todense()[1:, 1:]) self._lhs_inverted = True except LinAlgError: raise LinAlgError(f'The left-hand side of the layer {self} is not invertible with the provided basis.') From 894dc86a719d899556f7e52d6442987c4ae7b44e Mon Sep 17 00:00:00 2001 From: Jonathan Demaeyer Date: Fri, 12 Jun 2026 10:36:41 +0200 Subject: [PATCH 3/3] Preparing for v1.0.6 beta release --- documentation/source/conf.py | 2 +- layercake/__init__.py | 2 +- pyproject.toml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/documentation/source/conf.py b/documentation/source/conf.py index c8fc71d..2d358b7 100644 --- a/documentation/source/conf.py +++ b/documentation/source/conf.py @@ -19,7 +19,7 @@ project = 'LayerCake' copyright = '2025-2026, Jonathan Demaeyer and Oisín Hamilton' author = 'Jonathan Demaeyer and Oisín Hamilton' -release = 'v1.0.5-beta' +release = 'v1.0.6-beta' version = release # -- General configuration --------------------------------------------------- diff --git a/layercake/__init__.py b/layercake/__init__.py index 338471d..5899849 100644 --- a/layercake/__init__.py +++ b/layercake/__init__.py @@ -13,4 +13,4 @@ 'OperatorTerm', 'ProductOfTerms', 'AdditionOfTerms', 'LinearTerm', 'ConstantTerm', 'Equation', 'Laplacian', 'D', 'Layer', 'Cake'] -__version__ = '1.0.5b0' +__version__ = '1.0.6b0' diff --git a/pyproject.toml b/pyproject.toml index c96c38e..2b5c218 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -7,7 +7,7 @@ requires = ["setuptools", "wheel"] [project] requires-python = ">=3.10" name = "layercake_model" -version = "1.0.5-beta" +version = "1.0.6-beta" description = "A framework to design systems of partial differential equations (PDEs), and convert them to ordinary differential equations (ODEs) via Galerkin-type expansions. " readme = "README.md" authors = [