move LUDecomposition to legacy.py

Author: mozman

profiling/banded_matrix_solver.py

@@ -8,9 +8,10 @@
     Matrix,
     BandedMatrixLU,
     banded_matrix,
-    LUDecomposition,
+    NumpySolver,
 )
 
+
 CWD = pathlib.Path("~/Desktop/Outbox").expanduser()
 if not CWD.exists():
     CWD = pathlib.Path(".")
@@ -28,10 +29,10 @@ def random_matrix(shape, m1: int, m2: int):
     return m
 
 
-def profile_LU_matrix_solver(count: int, A: Matrix, B: Matrix):
+def profile_numpy_matrix_solver(count: int, A: Matrix, B: Matrix):
     for _ in range(count):
-        lu = LUDecomposition(A)
-        lu.solve_matrix(B)
+        lu = NumpySolver(A.matrix)
+        lu.solve_matrix(B.matrix)
 
 
 def profile_banded_matrix_solver(count, A: Matrix, B: Matrix):
@@ -76,11 +77,11 @@ def profile(func, *args):
                     )
                 )
             )
-            t0 = profile(profile_LU_matrix_solver, REPEAT, A, B)
+            t0 = profile(profile_numpy_matrix_solver, REPEAT, A, B)
             t1 = profile(profile_banded_matrix_solver, REPEAT, A, B)
             factor = t0 / t1
             print(
-                f"Matrix {size}x{size}, m1={m1}, m2={m2}, {REPEAT}x: Standard LU {t0:0.3f}s Banded LU {t1:0.3}s factor: x{factor:.1f}"
+                f"Matrix {size}x{size}, m1={m1}, m2={m2}, {REPEAT}x: Numpy {t0:0.3f}s Banded LU {t1:0.3}s factor: x{factor:.1f}"
             )
             writer.writerow(
                 [

profiling/linalg.py

@@ -5,7 +5,6 @@
 
 from ezdxf.math.linalg import (
     Matrix,
-    LUDecomposition,
     numpy_matrix_solver,
     numpy_vector_solver,
 )
@@ -14,6 +13,7 @@
     gauss_matrix_solver,
     gauss_jordan_solver,
     gauss_jordan_inverse,
+    LUDecomposition,
 )
 
 

src/ezdxf/math/bspline.py

@@ -22,7 +22,6 @@
     Iterator,
     Sequence,
     TYPE_CHECKING,
-    Union,
     Optional,
 )
 import math
@@ -41,7 +40,7 @@
 )
 from ezdxf.math import linalg
 from ezdxf.lldxf.const import DXFValueError
-from ezdxf import PYPY
+
 
 if TYPE_CHECKING:
     from ezdxf.math import (
@@ -51,11 +50,6 @@
         Bezier4P,
     )
 
-# Acceleration of banded diagonal matrix solver kicks in at:
-# N=15 for CPython on Windows and Linux
-# N=60 for pypy3 on Windows and Linux
-USE_BANDED_MATRIX_SOLVER_CPYTHON_LIMIT = 15
-USE_BANDED_MATRIX_SOLVER_PYPY_LIMIT = 60
 
 __all__ = [
     # High level functions:
@@ -551,24 +545,26 @@ def double_knots(n: int, p: int, t: Sequence[float]) -> list[float]:
     return u
 
 
-def _get_best_solver(matrix: Union[list, linalg.Matrix], degree: int):
-    """Returns best suited linear equation solver depending on matrix
-    configuration and python interpreter.
+def _get_best_solver(matrix: list| linalg.Matrix, degree: int) -> linalg.Solver:
+    """Returns best suited linear equation solver depending on matrix configuration and 
+    python interpreter.
     """
-    A = matrix if isinstance(matrix, linalg.Matrix) else linalg.Matrix(matrix=matrix)
-    if PYPY:
-        limit = USE_BANDED_MATRIX_SOLVER_PYPY_LIMIT
-    else:
-        limit = USE_BANDED_MATRIX_SOLVER_CPYTHON_LIMIT
-    if A.nrows < limit:  # use default equation solver
-        return linalg.LUDecomposition(A.matrix)
+    # v1.2: added NumpySolver
+    #   Acceleration of banded diagonal matrix solver is still a thing but only for 
+    #   really big matrices N > 30
+    # PyPy has no advantages when using the NumpySolver
+    if not isinstance(matrix, linalg.Matrix):
+        matrix =linalg.Matrix(matrix)
+
+    if matrix.nrows < 30:  # use default equation solver
+        return linalg.NumpySolver(matrix.matrix)
     else:
         # Theory: band parameters m1, m2 are at maximum degree-1, for
         # B-spline interpolation and approximation:
         # m1 = m2 = degree-1
         # But the speed gain is not that big and just to be sure:
-        m1, m2 = linalg.detect_banded_matrix(A, check_all=False)
-        A = linalg.compact_banded_matrix(A, m1, m2)
+        m1, m2 = linalg.detect_banded_matrix(matrix, check_all=False)
+        A = linalg.compact_banded_matrix(matrix, m1, m2)
         return linalg.BandedMatrixLU(A, m1, m2)
 
 
@@ -606,7 +602,8 @@ def unconstrained_global_bspline_interpolation(
     )
     N = Basis(knots=knots, order=degree + 1, count=len(fit_points))
     solver = _get_best_solver([N.basis_vector(t) for t in t_vector], degree)
-    control_points = solver.solve_matrix(fit_points)
+    mat_B = np.array(fit_points, dtype=np.float64)
+    control_points = solver.solve_matrix(mat_B)
     return Vec3.list(control_points.rows()), knots
 
 

src/ezdxf/math/legacy.py

@@ -2,22 +2,18 @@
 # License: MIT License
 # legacy.py code is replaced by numpy - exists just for benchmarking, testing and nostalgia
 from __future__ import annotations
-from typing import Iterable
+from typing import Iterable, cast
 from itertools import repeat
-from .linalg import Matrix, MatrixData, NDArray
+import reprlib
+
+from .linalg import Matrix, MatrixData, NDArray, Solver
 
 __all__ = [
     "gauss_vector_solver",
     "gauss_matrix_solver",
     "gauss_jordan_solver",
     "gauss_jordan_inverse",
     "LUDecomposition",
-    "tridiagonal_vector_solver",
-    "tridiagonal_matrix_solver",
-    "detect_banded_matrix",
-    "compact_banded_matrix",
-    "BandedMatrixLU",
-    "banded_matrix",
 ]
 
 
@@ -264,4 +260,158 @@ def gauss_jordan_inverse(A: MatrixData) -> Matrix:
     else:
         matrix_a = list(A)
     nrows = len(matrix_a)
-    return gauss_jordan_solver(matrix_a, list(repeat([0.0], nrows)))[0]
+    return gauss_jordan_solver(matrix_a, list(repeat([0.0], nrows)))[0]
+
+
+class LUDecomposition(Solver):
+    """Represents a `LU decomposition`_ matrix of A, raise :class:`ZeroDivisionError`
+    for a singular matrix.
+
+    This algorithm is a little bit faster than the `Gauss-Elimination`_
+    algorithm using CPython and much faster when using pypy.
+
+    The :attr:`LUDecomposition.matrix` attribute gives access to the matrix data
+    as list of rows like in the :class:`Matrix` class, and the :attr:`LUDecomposition.index`
+    attribute gives access to the swapped row indices.
+
+    Args:
+        A: matrix [[a11, a12, ..., a1n], [a21, a22, ..., a2n], [a21, a22, ..., a2n],
+            ... [an1, an2, ..., ann]]
+
+    raises:
+        ZeroDivisionError: singular matrix
+
+    """
+
+    __slots__ = ("matrix", "index", "_det")
+
+    def __init__(self, A: MatrixData | NDArray):
+        lu: MatrixData = copy_float_matrix(A)
+        n: int = len(lu)
+        det: float = 1.0
+        index: list[int] = []
+
+        # find max value for each row, raises ZeroDivisionError for singular matrix!
+        scaling: list[float] = [1.0 / max(abs(v) for v in row) for row in lu]
+
+        for k in range(n):
+            big: float = 0.0
+            imax: int = k
+            for i in range(k, n):
+                temp: float = scaling[i] * abs(lu[i][k])
+                if temp > big:
+                    big = temp
+                    imax = i
+
+            if k != imax:
+                for col in range(n):
+                    temp = lu[imax][col]
+                    lu[imax][col] = lu[k][col]
+                    lu[k][col] = temp
+
+                det = -det
+                scaling[imax] = scaling[k]
+
+            index.append(imax)
+            for i in range(k + 1, n):
+                temp = lu[i][k] / lu[k][k]
+                lu[i][k] = temp
+                for j in range(k + 1, n):
+                    lu[i][j] -= temp * lu[k][j]
+
+        self.index: list[int] = index
+        self.matrix: MatrixData = lu
+        self._det: float = det
+
+    def __str__(self) -> str:
+        return str(self.matrix)
+
+    def __repr__(self) -> str:
+        return f"{self.__class__} {reprlib.repr(self.matrix)}"
+
+    @property
+    def nrows(self) -> int:
+        """Count of matrix rows (and cols)."""
+        return len(self.matrix)
+
+    def solve_vector(self, B: Iterable[float]) -> list[float]:
+        """Solves the linear equation system given by the nxn Matrix A . x = B,
+        right-hand side quantities as vector B with n elements.
+
+        Args:
+            B: vector [b1, b2, ..., bn]
+
+        Returns:
+            vector as list of floats
+
+        """
+        X: list[float] = [float(v) for v in B]
+        lu: MatrixData = self.matrix
+        index: list[int] = self.index
+        n: int = self.nrows
+        ii: int = 0
+
+        if len(X) != n:
+            raise ValueError(
+                "Item count of vector B has to be equal to matrix row count."
+            )
+
+        for i in range(n):
+            ip: int = index[i]
+            sum_: float = X[ip]
+            X[ip] = X[i]
+            if ii != 0:
+                for j in range(ii - 1, i):
+                    sum_ -= lu[i][j] * X[j]
+            elif sum_ != 0.0:
+                ii = i + 1
+            X[i] = sum_
+
+        for row in range(n - 1, -1, -1):
+            sum_ = X[row]
+            for col in range(row + 1, n):
+                sum_ -= lu[row][col] * X[col]
+            X[row] = sum_ / lu[row][row]
+        return X
+
+    def solve_matrix(self, B: MatrixData | NDArray) -> Matrix:
+        """Solves the linear equation system given by the nxn Matrix A . x = B,
+        right-hand side quantities as nxm Matrix B.
+
+        Args:
+            B: matrix [[b11, b12, ..., b1m], [b21, b22, ..., b2m],
+                ... [bn1, bn2, ..., bnm]]
+
+        Returns:
+            matrix as :class:`Matrix` object
+
+        """
+        if not isinstance(B, Matrix):
+            matrix_b = Matrix(matrix=[list(row) for row in B])
+        else:
+            matrix_b = cast(Matrix, B)
+
+        if matrix_b.nrows != self.nrows:
+            raise ValueError("Row count of self and matrix B has to match.")
+
+        return Matrix(
+            matrix=[self.solve_vector(col) for col in matrix_b.cols()]
+        ).transpose()
+
+    def inverse(self) -> Matrix:
+        """Returns the inverse of matrix as :class:`Matrix` object,
+        raise :class:`ZeroDivisionError` for a singular matrix.
+
+        """
+        return self.solve_matrix(Matrix.identity(shape=(self.nrows, self.nrows)).matrix)
+
+    def determinant(self) -> float:
+        """Returns the determinant of matrix, raises :class:`ZeroDivisionError`
+        if matrix is singular.
+
+        """
+        det: float = self._det
+        lu: MatrixData = self.matrix
+        for i in range(self.nrows):
+            det *= lu[i][i]
+        return det