From 96f67f6bc8eceb6c6d82d63f892683524cbb7779 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 25 May 2024 17:42:11 +0200 Subject: [PATCH 01/34] pkg: made `.venv` in `.gitignore` more general --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index bdd5bb4..7eb4c23 100644 --- a/.gitignore +++ b/.gitignore @@ -83,7 +83,7 @@ celerybeat-schedule .env # virtualenv -.venv +.venv* venv/ ENV/ From 3deb05ca115a11e1d7ad419a9e3b7223a6ea680d Mon Sep 17 00:00:00 2001 From: MothNik Date: Fri, 7 Jun 2024 21:40:00 +0200 Subject: [PATCH 02/34] wip: [11] refactored algorithm 1 to handle multiple right-hand sides --- src/pentapy/core.py | 7 ++ src/pentapy/solver.pxd | 2 +- src/pentapy/solver.pyx | 272 ++++++++++++++++++++++++++++++++++------- 3 files changed, 237 insertions(+), 44 deletions(-) diff --git a/src/pentapy/core.py b/src/pentapy/core.py index 067189d..9a55780 100644 --- a/src/pentapy/core.py +++ b/src/pentapy/core.py @@ -77,7 +77,14 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): else: mat_flat = create_banded(mat, col_wise=False, dtype=np.double) rhs = np.asarray(rhs, dtype=np.double) + single_rhs = rhs.ndim == 1 + if single_rhs: + rhs = rhs[:, np.newaxis] + try: + if single_rhs: + return penta_solver1(mat_flat, rhs).ravel() + return penta_solver1(mat_flat, rhs) except ZeroDivisionError: warnings.warn("pentapy: PTRANS-I not suitable for input-matrix.") diff --git a/src/pentapy/solver.pxd b/src/pentapy/solver.pxd index e7c471e..05d249f 100644 --- a/src/pentapy/solver.pxd +++ b/src/pentapy/solver.pxd @@ -1,4 +1,4 @@ # cython: language_level=3 -cdef double[:] c_penta_solver1(double[:, :] mat_flat, double[:] rhs) +cdef double[:, :] c_penta_solver1(double[:, :] mat_flat, double[:, :] rhs) cdef double[:] c_penta_solver2(double[:, :] mat_flat, double[:] rhs) diff --git a/src/pentapy/solver.pyx b/src/pentapy/solver.pyx index 469b074..c59226b 100644 --- a/src/pentapy/solver.pyx +++ b/src/pentapy/solver.pyx @@ -1,14 +1,23 @@ -# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=True +# cython: language_level=3, boundscheck=True, wraparound=False, cdivision=True + """ This is a solver linear equation systems with a penta-diagonal matrix, -implemented in cython. +implemented in Cython. + """ + +### Imports ### + import numpy as np cimport numpy as np +from libc.stdint cimport int64_t, uint64_t -def penta_solver1(double[:, :] mat_flat, double[:] rhs): +### Main Python Interface ### + + +def penta_solver1(double[:, :] mat_flat, double[:, :] rhs): return np.asarray(c_penta_solver1(mat_flat, rhs)) @@ -16,56 +25,233 @@ def penta_solver2(double[:, :] mat_flat, double[:] rhs): return np.asarray(c_penta_solver2(mat_flat, rhs)) -cdef double[:] c_penta_solver1(double[:, :] mat_flat, double[:] rhs): - - cdef int mat_j = mat_flat.shape[1] - - cdef double[:] result = np.zeros(mat_j) +### Solver Algorithm 1 ### - cdef double[:] al = np.zeros(mat_j) - cdef double[:] be = np.zeros(mat_j) - cdef double[:] ze = np.zeros(mat_j) - cdef double[:] ga = np.zeros(mat_j) - cdef double[:] mu = np.zeros(mat_j) - cdef int i +cdef double[:, :] c_penta_solver1(double[:, :] mat_flat, double[:, :] rhs): + """ + Solves the pentadiagonal system of equations ``Ax = b`` with the matrix ``A`` and + the right-hand side ``b`` by - mu[0] = mat_flat[2, 0] - al[0] = mat_flat[1, 0] / mu[0] - be[0] = mat_flat[0, 0] / mu[0] - ze[0] = rhs[0] / mu[0] + - factorizing the matrix ``A`` into auxiliary coefficients and a unit upper + triangular matrix ``U`` + - transforming the right-hand side into a vector ``zeta`` + - solving the system of equations ``Ux = zeta`` by backward substitution - ga[1] = mat_flat[3, 1] - mu[1] = mat_flat[2, 1] - al[0] * ga[1] - al[1] = (mat_flat[1, 1] - be[0] * ga[1]) / mu[1] - be[1] = mat_flat[0, 1] / mu[1] - ze[1] = (rhs[1] - ze[0] * ga[1]) / mu[1] + """ - for i in range(2, mat_j-2): - ga[i] = mat_flat[3, i] - al[i-2] * mat_flat[4, i] - mu[i] = mat_flat[2, i] - be[i-2] * mat_flat[4, i] - al[i-1] * ga[i] - al[i] = (mat_flat[1, i] - be[i-1] * ga[i]) / mu[i] - be[i] = mat_flat[0, i] / mu[i] - ze[i] = (rhs[i] - ze[i-2] * mat_flat[4, i] - ze[i-1] * ga[i]) / mu[i] + cdef uint64_t mat_n_rows = mat_flat.shape[1] + cdef uint64_t rhs_n_cols = rhs.shape[1] + cdef uint64_t iter_col + cdef double[::, ::1] result = np.empty(shape=(mat_n_rows, rhs_n_cols)) + cdef double[::, ::1] mat_factorized = np.empty(shape=(mat_n_rows, 5)) - ga[mat_j-2] = mat_flat[3, mat_j-2] - al[mat_j-4] * mat_flat[4, mat_j-2] - mu[mat_j-2] = mat_flat[2, mat_j-2] - be[mat_j-4] * mat_flat[4, mat_j-2] - al[mat_j-3] * ga[mat_j-2] - al[mat_j-2] = (mat_flat[1, mat_j-2] - be[mat_j-3] * ga[mat_j-2]) / mu[mat_j-2] + # first, the matrix is factorized + c_penta_factorize_algo1( + mat_flat, + mat_n_rows, + mat_factorized, + ) - ga[mat_j-1] = mat_flat[3, mat_j-1] - al[mat_j-3] * mat_flat[4, mat_j-1] - mu[mat_j-1] = mat_flat[2, mat_j-1] - be[mat_j-3] * mat_flat[4, mat_j-1] - al[mat_j-2] * ga[mat_j-1] + # then, all the right-hand sides are solved + for iter_col in range(rhs_n_cols): + c_solve_penta_from_factorize_algo_1( + mat_n_rows, + mat_factorized, + rhs[::, iter_col], + result[::, iter_col], + ) - ze[mat_j-2] = (rhs[mat_j-2] - ze[mat_j-4] * mat_flat[4, mat_j-2] - ze[mat_j-3] * ga[mat_j-2]) / mu[mat_j-2] - ze[mat_j-1] = (rhs[mat_j-1] - ze[mat_j-3] * mat_flat[4, mat_j-1] - ze[mat_j-2] * ga[mat_j-1]) / mu[mat_j-1] - - # Backward substitution - result[mat_j-1] = ze[mat_j-1] - result[mat_j-2] = ze[mat_j-2] - al[mat_j-2] * result[mat_j-1] + return result - for i in range(mat_j-3, -1, -1): - result[i] = ze[i] - al[i] * result[i+1] - be[i] * result[i+2] - return result +cdef void c_penta_factorize_algo1( + double[:, :] mat_flat, + uint64_t mat_n_rows, + double[::, ::1] mat_factorized, +): + """ + Factorizes the pentadiagonal matrix ``A`` into + + - auxiliary coefficients ``e``, ``mu`` and ``gamma`` for the transformation of the + right-hand side + - a unit upper triangular matrix with the main diagonals ``alpha`` and ``beta`` + for the following backward substitution. Its unit main diagonal is implicit. + + They are overwriting the memoryview ``mat_factorized`` as follows: + + ```bash + [[ * mu_0 * al_0 be_0 ] + [ * mu_1 ga_1 al_1 be_1 ] + [ e_2 mu_2 ga_2 al_2 be_2 ] + ... + [ e_i mu_i ga_i al_i be_i ] + ... + [ e_{n-2} mu_{n-2} ga_{n-2} al_{n-2} * ] + [ e_{n-1} mu_{n-1} ga_{n-1} * * ]] + ``` + + where the entries marked with ``*`` are not used by design, but overwritten with + zeros. + + """ + + ### Variable declarations ### + + cdef uint64_t iter_row + cdef double mu_i, ga_i, e_i + cdef double al_i, al_i_minus_1, al_i_plus_1 + + ### Factorization ### + + # First row + mu_i = mat_flat[2, 0] + al_i_minus_1 = mat_flat[1, 0] / mu_i + be_i_minus_1 = mat_flat[0, 0] / mu_i + + mat_factorized[0, 0] = 0.0 + mat_factorized[0, 1] = mu_i + mat_factorized[0, 2] = 0.0 + mat_factorized[0, 3] = al_i_minus_1 + mat_factorized[0, 4] = be_i_minus_1 + + # Second row + ga_i = mat_flat[3, 1] + mu_i = mat_flat[2, 1] - al_i_minus_1 * ga_i + al_i = (mat_flat[1, 1] - be_i_minus_1 * ga_i) / mu_i + be_i = mat_flat[0, 1] / mu_i + + mat_factorized[1, 0] = 0.0 + mat_factorized[1, 1] = mu_i + mat_factorized[1, 2] = ga_i + mat_factorized[1, 3] = al_i + mat_factorized[1, 4] = be_i + + # Central rows + for iter_row in range(2, mat_n_rows-2): + e_i = mat_flat[4, iter_row] + ga_i = mat_flat[3, iter_row] - al_i_minus_1 * e_i + mu_i = mat_flat[2, iter_row] - be_i_minus_1 * e_i - al_i * ga_i + + al_i_plus_1 = (mat_flat[1, iter_row] - be_i * ga_i) / mu_i + al_i_minus_1 = al_i + al_i = al_i_plus_1 + + be_i_plus_1 = mat_flat[0, iter_row] / mu_i + be_i_minus_1 = be_i + be_i = be_i_plus_1 + + mat_factorized[iter_row, 0] = e_i + mat_factorized[iter_row, 1] = mu_i + mat_factorized[iter_row, 2] = ga_i + mat_factorized[iter_row, 3] = al_i + mat_factorized[iter_row, 4] = be_i + + # Second to last row + e_i = mat_flat[4, mat_n_rows-2] + ga_i = mat_flat[3, mat_n_rows-2] - al_i_minus_1 * e_i + mu_i = mat_flat[2, mat_n_rows-2] - be_i_minus_1 * e_i - al_i * ga_i + al_i_plus_1 = (mat_flat[1, mat_n_rows-2] - be_i * ga_i) / mu_i + + mat_factorized[mat_n_rows-2, 0] = e_i + mat_factorized[mat_n_rows-2, 1] = mu_i + mat_factorized[mat_n_rows-2, 2] = ga_i + mat_factorized[mat_n_rows-2, 3] = al_i_plus_1 + mat_factorized[mat_n_rows-2, 4] = 0.0 + + # Last Row + e_i = mat_flat[4, mat_n_rows-1] + ga_i = mat_flat[3, mat_n_rows-1] - al_i * e_i + mu_i = mat_flat[2, mat_n_rows-1] - be_i * e_i - al_i_plus_1 * ga_i + + mat_factorized[mat_n_rows-1, 0] = e_i + mat_factorized[mat_n_rows-1, 1] = mu_i + mat_factorized[mat_n_rows-1, 2] = ga_i + mat_factorized[mat_n_rows-1, 3] = 0.0 + mat_factorized[mat_n_rows-1, 4] = 0.0 + + return + + +cdef void c_solve_penta_from_factorize_algo_1( + uint64_t mat_n_rows, + double[::, ::1] mat_factorized, + double[::] rhs_single, + double[::] result_view, +): + """ + Solves the pentadiagonal system of equations ``Ax = b`` with the factorized + unit upper triangular matrix ``U`` and the right-hand side ``b``. + It overwrites the right-hand side ``b`` first with the transformed vector ``zeta`` + and then with the solution vector ``x`` for ``Ux = zeta``. + + """ + + ### Variable declarations ### + + cdef int64_t iter_row + cdef double ze_i, ze_i_minus_1, ze_i_plus_1 + + ### Transformation ### + + # first, the right-hand side is transformed into the vector ``zeta`` + # First row + + ze_i_minus_1 = rhs_single[0] / mat_factorized[0, 1] + result_view[0] = ze_i_minus_1 + + # Second row + ze_i = (rhs_single[1] - ze_i_minus_1 * mat_factorized[1, 2]) / mat_factorized[1, 1] + result_view[1] = ze_i + + # Central rows + for iter_row in range(2, mat_n_rows-2): + ze_i_plus_1 = ( + rhs_single[iter_row] + - ze_i_minus_1 * mat_factorized[iter_row, 0] + - ze_i * mat_factorized[iter_row, 2] + ) / mat_factorized[iter_row, 1] + ze_i_minus_1 = ze_i + ze_i = ze_i_plus_1 + result_view[iter_row] = ze_i_plus_1 + + # Second to last row + ze_i_plus_1 = ( + rhs_single[mat_n_rows-2] + - ze_i_minus_1 * mat_factorized[mat_n_rows-2, 0] + - ze_i * mat_factorized[mat_n_rows-2, 2] + ) / mat_factorized[mat_n_rows-2, 1] + ze_i_minus_1 = ze_i + ze_i = ze_i_plus_1 + result_view[mat_n_rows-2] = ze_i_plus_1 + + # Last row + ze_i_plus_1 = ( + rhs_single[mat_n_rows-1] + - ze_i_minus_1 * mat_factorized[mat_n_rows-1, 0] + - ze_i * mat_factorized[mat_n_rows-1, 2] + ) / mat_factorized[mat_n_rows-1, 1] + result_view[mat_n_rows-1] = ze_i_plus_1 + + ### Backward substitution ### + + # The solution vector is calculated by backward substitution that overwrites the + # right-hand side vector with the solution vector + ze_i -= mat_factorized[mat_n_rows-2, 3] * ze_i_plus_1 + result_view[mat_n_rows-2] = ze_i + + for iter_row in range(mat_n_rows-3, -1, -1): + result_view[iter_row] -= ( + mat_factorized[iter_row, 3] * ze_i + + mat_factorized[iter_row, 4] * ze_i_plus_1 + ) + ze_i_plus_1 = ze_i + ze_i = result_view[iter_row] + + return + + +### Solver Algorithm 2 ### cdef double[:] c_penta_solver2(double[:, :] mat_flat, double[:] rhs): From 328d9c90492ec4dd10a3a62f50917547c63082eb Mon Sep 17 00:00:00 2001 From: MothNik Date: Fri, 7 Jun 2024 21:40:29 +0200 Subject: [PATCH 03/34] feat: [11] added doctest runs to `pytest` --- pytest.ini | 2 ++ 1 file changed, 2 insertions(+) create mode 100644 pytest.ini diff --git a/pytest.ini b/pytest.ini new file mode 100644 index 0000000..2bed0f3 --- /dev/null +++ b/pytest.ini @@ -0,0 +1,2 @@ +[pytest] +addopts = --doctest-modules \ No newline at end of file From 90bdeb143f8677427d4255e332a00d83b51a40cf Mon Sep 17 00:00:00 2001 From: MothNik Date: Fri, 7 Jun 2024 22:33:59 +0200 Subject: [PATCH 04/34] tests: [11] replaced single tools test by scalable parametrized tools test that test more different cases --- tests/test_tools.py | 205 ++++++++++++++++++++++++++++++++++++++++++++ tests/util_funcs.py | 175 +++++++++++++++++++++++++++++++++++++ 2 files changed, 380 insertions(+) create mode 100644 tests/test_tools.py create mode 100644 tests/util_funcs.py diff --git a/tests/test_tools.py b/tests/test_tools.py new file mode 100644 index 0000000..879a713 --- /dev/null +++ b/tests/test_tools.py @@ -0,0 +1,205 @@ +""" +This test suite implements the test for the ``tools`` module of the ``pentapy`` package. + +""" + +### Imports ### + +import warnings +from typing import Optional, Tuple, Type + +import numpy as np +import pentapy as pp +import pytest +import util_funcs as uf +from pentapy.tools import _check_penta + +warnings.simplefilter("always") + +### Constants ### + +SEED = 19_031_977 +N_ROWS = [ + 3, + 4, + 5, + 10, + 11, + 25, + 26, + 50, + 51, + 100, + 101, + 250, + 251, + 500, + 501, + 1_000, + 1_001, + 10_000, + 10_001, +] + +### Tests ### + + +@pytest.mark.parametrize("offset", [0, 1, 2, -1, -2]) +@pytest.mark.parametrize("n_rows", N_ROWS) +def test_diag_indices(n_rows: int, offset: int) -> None: + """ + Tests the generation of the diagonal indices via the function + ``pentapy.diag_indices``. + + """ + + # the diagonal indices are obtained with NumPy and pentapy + row_idxs_ref, col_idxs_ref = uf.get_diag_indices(n=n_rows, offset=offset) + row_idxs, col_idxs = pp.diag_indices(n=n_rows, offset=offset) + + # the diagonal indices are compared + assert np.array_equal(row_idxs_ref, row_idxs) + assert np.array_equal(col_idxs_ref, col_idxs) + + +@pytest.mark.parametrize("copy", [True, False]) +@pytest.mark.parametrize("with_shift", [True, False]) +@pytest.mark.parametrize("col_wise", [True, False]) +@pytest.mark.parametrize("n_rows", N_ROWS) +def test_penta_generators( + n_rows: int, + col_wise: bool, + with_shift: bool, + copy: bool, +) -> None: + """ + Tests the generation of pentadiagonal matrices where the matrix. + + """ + + # a reference matrix is initialised + mat_full_ref = uf.gen_rand_penta_matrix_dense_int( + n_rows=n_rows, + seed=SEED, + with_pentapy_indices=False, + ) + + # then, it is turned into a banded matrix ... + mat_banded = pp.create_banded(mat_full_ref, col_wise=col_wise) + + # ... which is maybe shifted + # Case 1: copied shift + if with_shift and copy: + mat_banded = pp.shift_banded(mat_banded, col_to_row=col_wise, copy=True) + col_wise = not col_wise + + # Case 2: in-place shift + if with_shift and not copy: + mat_banded = pp.shift_banded(mat_banded, col_to_row=col_wise, copy=False) + col_wise = not col_wise + + # ... from which a full matrix is created again + mat_full = pp.create_full(mat_banded, col_wise=col_wise) + + # the matrices are compared + assert np.array_equal(mat_full_ref, mat_full) + + +@pytest.mark.parametrize( + "shape, exception", + [ + ((5, 5), None), # Valid 2D Array with 5 rows and 5 rows + ((5, 2), ValueError), # 2D Array with 5 rows but only 2 columns + ((2, 5), ValueError), # 2D Array with 2 rows but 5 columns + ((5,), ValueError), # 1D Array + ], +) +def test_create_banded_raises( + shape: Tuple[int, ...], + exception: Optional[Type[Exception]], +) -> None: + """ + Test if the function ``pentapy.create_banded`` raises the expected exceptions. + + """ + + # the test matrix is initialised + np.random.seed(SEED) + mat = np.random.rand(*shape) + + # Case 1: no exception should be raised + if exception is None: + pp.create_banded(mat) + return + + # Case 2: an exception should be raised + with pytest.raises(exception): + pp.create_banded(mat) + + +@pytest.mark.parametrize( + "shape, exception", + [ + ((5, 5), None), # Valid 2D Array with 5 bands and 5 columns + ((5, 10), None), # Valid 2D Array with 5 bands and 10 columns + ((5, 3), None), # 2D Array with 5 bands and the minimum number of columns + ((6, 20), ValueError), # 2D Array does not have 5 bands + ((4, 30), ValueError), # 2D Array does not have 5 bands + ((5, 1), ValueError), # 2D Array with 5 bands but too little columns + ((5, 2), ValueError), # 2D Array with 5 bands but too little columns + ((5,), ValueError), # 1D Array + ], +) +def test_create_full_raises( + shape: Tuple[int, ...], + exception: Optional[Type[Exception]], +) -> None: + """ + Test if the function ``pentapy.create_full`` raises the expected exceptions. + + """ + + # the test matrix is initialised + np.random.seed(SEED) + mat = np.random.rand(*shape) + + # Case 1: no exception should be raised + if exception is None: + pp.create_full(mat) + return + + # Case 2: an exception should be raised + with pytest.raises(exception): + pp.create_full(mat) + + +@pytest.mark.parametrize( + "shape, exception", + [ + ((5, 3), None), # Valid 2D Array with 5 bands and 3 rows + ((5, 2), ValueError), # 2D Array with 5 bands but less than 3 rows + ((4, 3), ValueError), # 2D Array with less than 5 bands + ((5,), ValueError), # 1D Array + ], +) +def test_check_penta( + shape: Tuple[int, ...], + exception: Optional[Type[Exception]], +) -> None: + """ + Test if the function ``pentapy.tools._check_penta`` raises the expected exceptions. + + """ + + # the test matrix is initialised + np.random.seed(SEED) + mat = np.random.rand(*shape) + + # Case 1: no exception should be raised + if exception is None: + _check_penta(mat) + return + + # Case 2: an exception should be raised + with pytest.raises(exception): + _check_penta(mat) diff --git a/tests/util_funcs.py b/tests/util_funcs.py new file mode 100644 index 0000000..2402d25 --- /dev/null +++ b/tests/util_funcs.py @@ -0,0 +1,175 @@ +""" +This test suite implements the utility functions for testing the ``pentapy`` package. + +""" + +### Imports ### + +from functools import partial +from typing import Tuple + +import numpy as np +import pentapy as pp + +### Utility Functions ### + + +def get_diag_indices( + n: int, + offset: int, +) -> Tuple[np.ndarray, np.ndarray]: + """ + Computes the row and column indices of the diagonal of a matrix ``mat``. + + This answer is based on the Stack Overflow answer that can be found at: + https://stackoverflow.com/a/18081653/14814813 + + Doctests + -------- + >>> # Setting up a test matrix + >>> n_rows = 5 + >>> mat = np.arange(start=0, stop=n_rows * n_rows).reshape(n_rows, n_rows) + + >>> # Getting the main diagonal indices + >>> row_idxs, col_idxs = get_diag_indices(n=n_rows, offset=0) + >>> row_idxs + array([0, 1, 2, 3, 4]) + >>> col_idxs + array([0, 1, 2, 3, 4]) + >>> mat[row_idxs, col_idxs] + array([ 0, 6, 12, 18, 24]) + + >>> # Getting the first upper diagonal indices + >>> row_idxs, col_idxs = get_diag_indices(n=n_rows, offset=1) + >>> row_idxs + array([0, 1, 2, 3]) + >>> col_idxs + array([1, 2, 3, 4]) + >>> mat[row_idxs, col_idxs] + array([ 1, 7, 13, 19]) + + >>> # Getting the second upper diagonal indices + >>> row_idxs, col_idxs = get_diag_indices(n=n_rows, offset=2) + >>> row_idxs + array([0, 1, 2]) + >>> col_idxs + array([2, 3, 4]) + >>> mat[row_idxs, col_idxs] + array([ 2, 8, 14]) + + >>> # Getting the first lower diagonal indices + >>> row_idxs, col_idxs = get_diag_indices(n=n_rows, offset=-1) + >>> row_idxs + array([1, 2, 3, 4]) + >>> col_idxs + array([0, 1, 2, 3]) + >>> mat[row_idxs, col_idxs] + array([ 5, 11, 17, 23]) + + >>> # Getting the second lower diagonal indices + >>> row_idxs, col_idxs = get_diag_indices(n=n_rows, offset=-2) + >>> row_idxs + array([2, 3, 4]) + >>> col_idxs + array([0, 1, 2]) + >>> mat[row_idxs, col_idxs] + array([10, 16, 22]) + + """ + + row_idxs, col_idxs = np.diag_indices(n=n, ndim=2) + if offset < 0: + row_idx_from = -offset + row_idx_to = None + col_idx_from = 0 + col_idx_to = offset + elif offset > 0: + row_idx_from = 0 + row_idx_to = -offset + col_idx_from = offset + col_idx_to = None + else: + row_idx_from = None + row_idx_to = None + col_idx_from = None + col_idx_to = None + + return ( + row_idxs[row_idx_from:row_idx_to], + col_idxs[col_idx_from:col_idx_to], + ) + + +def gen_rand_penta_matrix_dense_int( + n_rows: int, + seed: int, + with_pentapy_indices: bool, +) -> np.ndarray: + """ + Generates a random dense pentadiagonal matrix with shape ``(n_rows, n_rows)`` and + data type ``int64``. + + Doctests + -------- + >>> # Generating a random pentadiagonal matrix with NumPy indices + >>> n_rows = 5 + >>> seed = 19_031_977 + >>> with_pentapy_indices = False + + >>> mat_no_pentapy = gen_rand_penta_matrix_dense_int( + ... n_rows=n_rows, + ... seed=seed, + ... with_pentapy_indices=with_pentapy_indices + ... ) + >>> mat_no_pentapy + array([[117, 499, 43, 0, 0], + [378, 149, 857, 353, 0], + [285, 769, 767, 229, 484], + [ 0, 717, 214, 243, 877], + [ 0, 0, 410, 611, 79]], dtype=int64) + + >>> # Generating a random pentadiagonal matrix with pentapy indices + >>> mat_with_pentapy = gen_rand_penta_matrix_dense_int( + ... n_rows=n_rows, + ... seed=seed, + ... with_pentapy_indices=True + ... ) + >>> mat_with_pentapy + array([[117, 499, 43, 0, 0], + [378, 149, 857, 353, 0], + [285, 769, 767, 229, 484], + [ 0, 717, 214, 243, 877], + [ 0, 0, 410, 611, 79]], dtype=int64) + + >>> # Checking if the two matrices are equal + >>> np.array_equal(mat_no_pentapy, mat_with_pentapy) + True + + """ + + # first, a matrix of zeros is initialised ... + mat = np.zeros((n_rows, n_rows), dtype=np.int64) + # ... together with a partially specified random vector generator + # NOTE: this ensures consistent random numbers for both cases + gen_rand_int = partial(np.random.randint, low=1, high=1_000) + + # then, the diagonal index function is obtained + diag_idx_func = get_diag_indices + if with_pentapy_indices: + diag_idx_func = pp.diag_indices + + # then, the diagonals are filled with random integers + np.random.seed(seed=seed) + for offset in range(-2, 3): + row_idxs, col_idxs = diag_idx_func(n=n_rows, offset=offset) + mat[row_idxs, col_idxs] = gen_rand_int(size=n_rows - abs(offset)) + + return mat + + +### Doctests ### + +if __name__ == "__main__": # pragma: no cover + import doctest + + doctest.testmod() From 50806411c1087036776576f7c4e68a05d0f3b6a6 Mon Sep 17 00:00:00 2001 From: MothNik Date: Fri, 7 Jun 2024 22:35:43 +0200 Subject: [PATCH 05/34] tests: [11] removed version from coverage --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index 7b0aec6..4c400f8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -103,6 +103,7 @@ max-line-length = 120 "*examples*", "*tests*", "*paper*", + "pentapy/src/pentapy/_version.py", ] [tool.coverage.report] From 2f576092dfc0bb52e4ab268248e0673dc0eda667 Mon Sep 17 00:00:00 2001 From: MothNik Date: Fri, 7 Jun 2024 23:02:53 +0200 Subject: [PATCH 06/34] lint: [11] fixed block comment lint error --- src/pentapy/solver.pyx | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/pentapy/solver.pyx b/src/pentapy/solver.pyx index c59226b..bf075da 100644 --- a/src/pentapy/solver.pyx +++ b/src/pentapy/solver.pyx @@ -6,7 +6,7 @@ implemented in Cython. """ -### Imports ### +# Imports import numpy as np @@ -14,7 +14,7 @@ cimport numpy as np from libc.stdint cimport int64_t, uint64_t -### Main Python Interface ### +# Main Python Interface def penta_solver1(double[:, :] mat_flat, double[:, :] rhs): @@ -25,7 +25,7 @@ def penta_solver2(double[:, :] mat_flat, double[:] rhs): return np.asarray(c_penta_solver2(mat_flat, rhs)) -### Solver Algorithm 1 ### +# Solver Algorithm 1 cdef double[:, :] c_penta_solver1(double[:, :] mat_flat, double[:, :] rhs): @@ -96,13 +96,13 @@ cdef void c_penta_factorize_algo1( """ - ### Variable declarations ### + # Variable declarations cdef uint64_t iter_row cdef double mu_i, ga_i, e_i cdef double al_i, al_i_minus_1, al_i_plus_1 - ### Factorization ### + # Factorization # First row mu_i = mat_flat[2, 0] @@ -187,12 +187,12 @@ cdef void c_solve_penta_from_factorize_algo_1( """ - ### Variable declarations ### + # Variable declarations cdef int64_t iter_row cdef double ze_i, ze_i_minus_1, ze_i_plus_1 - ### Transformation ### + # Transformation # first, the right-hand side is transformed into the vector ``zeta`` # First row @@ -233,7 +233,7 @@ cdef void c_solve_penta_from_factorize_algo_1( ) / mat_factorized[mat_n_rows-1, 1] result_view[mat_n_rows-1] = ze_i_plus_1 - ### Backward substitution ### + # Backward substitution # The solution vector is calculated by backward substitution that overwrites the # right-hand side vector with the solution vector @@ -251,7 +251,7 @@ cdef void c_solve_penta_from_factorize_algo_1( return -### Solver Algorithm 2 ### +# Solver Algorithm 2 cdef double[:] c_penta_solver2(double[:, :] mat_flat, double[:] rhs): From 8db6ba4841b697725a60ba7d205b7db59f4b3d79 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 10:34:23 +0200 Subject: [PATCH 07/34] style: [11] improved headlines --- src/pentapy/solver.pyx | 16 ++++++++-------- tests/test_tools.py | 6 +++--- tests/util_funcs.py | 10 +++++----- 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/src/pentapy/solver.pyx b/src/pentapy/solver.pyx index bf075da..c8a22a2 100644 --- a/src/pentapy/solver.pyx +++ b/src/pentapy/solver.pyx @@ -14,7 +14,7 @@ cimport numpy as np from libc.stdint cimport int64_t, uint64_t -# Main Python Interface +# === Main Python Interface === def penta_solver1(double[:, :] mat_flat, double[:, :] rhs): @@ -25,7 +25,7 @@ def penta_solver2(double[:, :] mat_flat, double[:] rhs): return np.asarray(c_penta_solver2(mat_flat, rhs)) -# Solver Algorithm 1 +# === Solver Algorithm 1 === cdef double[:, :] c_penta_solver1(double[:, :] mat_flat, double[:, :] rhs): @@ -96,13 +96,13 @@ cdef void c_penta_factorize_algo1( """ - # Variable declarations + # === Variable declarations === cdef uint64_t iter_row cdef double mu_i, ga_i, e_i cdef double al_i, al_i_minus_1, al_i_plus_1 - # Factorization + # === Factorization === # First row mu_i = mat_flat[2, 0] @@ -187,12 +187,12 @@ cdef void c_solve_penta_from_factorize_algo_1( """ - # Variable declarations + # === Variable declarations === cdef int64_t iter_row cdef double ze_i, ze_i_minus_1, ze_i_plus_1 - # Transformation + # === Transformation === # first, the right-hand side is transformed into the vector ``zeta`` # First row @@ -233,7 +233,7 @@ cdef void c_solve_penta_from_factorize_algo_1( ) / mat_factorized[mat_n_rows-1, 1] result_view[mat_n_rows-1] = ze_i_plus_1 - # Backward substitution + # === Backward substitution === # The solution vector is calculated by backward substitution that overwrites the # right-hand side vector with the solution vector @@ -251,7 +251,7 @@ cdef void c_solve_penta_from_factorize_algo_1( return -# Solver Algorithm 2 +# === Solver Algorithm 2 === cdef double[:] c_penta_solver2(double[:, :] mat_flat, double[:] rhs): diff --git a/tests/test_tools.py b/tests/test_tools.py index 879a713..2f54c48 100644 --- a/tests/test_tools.py +++ b/tests/test_tools.py @@ -3,7 +3,7 @@ """ -### Imports ### +# === Imports === import warnings from typing import Optional, Tuple, Type @@ -16,7 +16,7 @@ warnings.simplefilter("always") -### Constants ### +# === Constants === SEED = 19_031_977 N_ROWS = [ @@ -41,7 +41,7 @@ 10_001, ] -### Tests ### +# === Tests === @pytest.mark.parametrize("offset", [0, 1, 2, -1, -2]) diff --git a/tests/util_funcs.py b/tests/util_funcs.py index 2402d25..594e75d 100644 --- a/tests/util_funcs.py +++ b/tests/util_funcs.py @@ -3,7 +3,7 @@ """ -### Imports ### +# === Imports === from functools import partial from typing import Tuple @@ -11,7 +11,7 @@ import numpy as np import pentapy as pp -### Utility Functions ### +# === Utility Functions === def get_diag_indices( @@ -25,7 +25,7 @@ def get_diag_indices( https://stackoverflow.com/a/18081653/14814813 Doctests - -------- + ======-- >>> # Setting up a test matrix >>> n_rows = 5 >>> mat = np.arange(start=0, stop=n_rows * n_rows).reshape(n_rows, n_rows) @@ -110,7 +110,7 @@ def gen_rand_penta_matrix_dense_int( data type ``int64``. Doctests - -------- + ======-- >>> # Generating a random pentadiagonal matrix with NumPy indices >>> n_rows = 5 >>> seed = 19_031_977 @@ -167,7 +167,7 @@ def gen_rand_penta_matrix_dense_int( return mat -### Doctests ### +# === Doctests === if __name__ == "__main__": # pragma: no cover import doctest From 1cdd525f8a93611e7c369253b724572388841819 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 10:35:27 +0200 Subject: [PATCH 08/34] wip: [11] formatted --- tests/util_funcs.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/util_funcs.py b/tests/util_funcs.py index 594e75d..040bd04 100644 --- a/tests/util_funcs.py +++ b/tests/util_funcs.py @@ -9,6 +9,7 @@ from typing import Tuple import numpy as np + import pentapy as pp # === Utility Functions === From ae656cf94f38b4651545161e848443a97930b53f Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 11:13:44 +0200 Subject: [PATCH 09/34] test: [11] created conditioned banded matrix creator for testing purposes --- tests/util_funcs.py | 199 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 197 insertions(+), 2 deletions(-) diff --git a/tests/util_funcs.py b/tests/util_funcs.py index 040bd04..ca8f2be 100644 --- a/tests/util_funcs.py +++ b/tests/util_funcs.py @@ -9,9 +9,15 @@ from typing import Tuple import numpy as np +from scipy import sparse as sprs import pentapy as pp +# === Constants === + +_MIN_DIAG_VAL = 1e-3 + + # === Utility Functions === @@ -26,7 +32,7 @@ def get_diag_indices( https://stackoverflow.com/a/18081653/14814813 Doctests - ======-- + -------- >>> # Setting up a test matrix >>> n_rows = 5 >>> mat = np.arange(start=0, stop=n_rows * n_rows).reshape(n_rows, n_rows) @@ -111,7 +117,7 @@ def gen_rand_penta_matrix_dense_int( data type ``int64``. Doctests - ======-- + -------- >>> # Generating a random pentadiagonal matrix with NumPy indices >>> n_rows = 5 >>> seed = 19_031_977 @@ -168,6 +174,195 @@ def gen_rand_penta_matrix_dense_int( return mat +def gen_conditioned_rand_penta_matrix_dense( + n_rows: int, + seed: int, + ill_conditioned: bool, +) -> np.ndarray: + """ + Generates a well- or ill-conditioned random banded pentadiagonal matrix with shape + ``(n_rows, n_rows)``. + + This is achieved as follows: + - a fake LDU decomposition is generated where ``L`` and ``U`` are unit lower and + upper triangular matrices, respectively, and ``D`` is a diagonal matrix + - the matrix is then reconstructed by multiplying the three matrices and converting + the result to a banded matrix + + If ``D`` does not have any zeros or values of small magnitude compared to the + largest value, the matrix should be well-conditioned. + Otherwise, it is ill-conditioned. + + Doctests + -------- + >>> # Imports + >>> from scipy.linalg import bandwidth + + >>> # 1) Generating a super small well-conditioned random pentadiagonal matrix + >>> n_rows = 3 + >>> seed = 19_031_977 + + >>> mat = gen_conditioned_rand_penta_matrix_dense( + ... n_rows=n_rows, + ... seed=seed, + ... ill_conditioned=False, + ... ) + >>> mat + array([[ 0.92453713, 0.28308514, -0.09972199], + [-0.09784268, 0.2270634 , -0.1509019 ], + [-0.23431267, 0.00468463, 0.22991003]]) + >>> # its bandwidth is computed and should be equal to 2 + >>> bandwidth(mat) + (2, 2) + >>> # its condition number is computed and values below 1e10 can be considered good + >>> np.linalg.cond(mat) + 4.976880305142543 + + >>> # 2) Generating a super small ill-conditioned random pentadiagonal matrix + >>> mat = gen_conditioned_rand_penta_matrix_dense( + ... n_rows=n_rows, + ... seed=seed, + ... ill_conditioned=True, + ... ) + >>> mat + array([[ 0.92453713, 0.28308514, -0.09972199], + [-0.09784268, 0.2270634 , -0.1509019 ], + [-0.23431267, 0.00468463, -0.02273771]]) + >>> # its bandwidth is computed and should be equal to 2 + >>> bandwidth(mat) + (2, 2) + >>> # its condition number is computed and its value should be close to the + >>> # reciprocal floating point precision, i.e., ~1e16 + >>> np.linalg.cond(mat) + 1.493156437173682e+17 + + >>> # 3) Generating a small well-conditioned random pentadiagonal matrix + >>> n_rows = 7 + + >>> mat = gen_conditioned_rand_penta_matrix_dense( + ... n_rows=n_rows, + ... seed=seed, + ... ill_conditioned=False, + ... ) + >>> np.round(mat, 2) + array([[ 0.92, -0.72, 0.73, 0. , 0. , 0. , 0. ], + [ 0.83, -0.02, 1.08, 0.41, 0. , 0. , 0. ], + [-0.58, 0.13, -0.13, -0.37, 0.18, 0. , 0. ], + [ 0. , -0.07, -0.58, 0.46, -0.31, 0.28, 0. ], + [ 0. , 0. , 0.43, 0.13, 0.39, -0.1 , -0.15], + [ 0. , 0. , 0. , 0.06, -0.14, 0.4 , 0.28], + [ 0. , 0. , 0. , 0. , -0.14, 0.36, 0.53]]) + >>> # its bandwidth is computed and should be equal to 2 + >>> bandwidth(mat) + (2, 2) + >>> # its condition number is computed and values below 1e10 can be considered good + >>> np.linalg.cond(mat) + 42.4847446467131 + + >>> # 4) Generating a small ill-conditioned random pentadiagonal matrix + >>> mat = gen_conditioned_rand_penta_matrix_dense( + ... n_rows=n_rows, + ... seed=seed, + ... ill_conditioned=True, + ... ) + >>> np.round(mat, 2) + array([[ 0.92, -0.72, 0.73, 0. , 0. , 0. , 0. ], + [ 0.83, -0.02, 1.08, 0.41, 0. , 0. , 0. ], + [-0.58, 0.13, -0.13, -0.37, 0.18, 0. , 0. ], + [ 0. , -0.07, -0.58, 0.46, -0.31, 0.28, 0. ], + [ 0. , 0. , 0.43, 0.13, 0.39, -0.1 , -0.15], + [ 0. , 0. , 0. , 0.06, -0.14, 0.4 , 0.28], + [ 0. , 0. , 0. , 0. , -0.14, 0.36, 0.28]]) + >>> # its bandwidth is computed and should be equal to 2 + >>> bandwidth(mat) + (2, 2) + >>> # its condition number is computed and its value should be close to the + >>> # reciprocal floating point precision, i.e., ~1e16 + >>> np.linalg.cond(mat) + 1.1079218802103074e+17 + + >>> # 5) Generating a large well-conditioned random pentadiagonal matrix + >>> n_rows = 1_000 + + >>> mat = gen_conditioned_rand_penta_matrix_dense( + ... n_rows=n_rows, + ... seed=seed, + ... ill_conditioned=False, + ... ) + >>> # its bandwidth is computed and should be equal to 2 + >>> bandwidth(mat) + (2, 2) + >>> # its condition number is computed and values below 1e10 can be considered good + >>> np.linalg.cond(mat) + 9570.995402466417 + + >>> # 6) Generating a large ill-conditioned random pentadiagonal matrix + >>> mat = gen_conditioned_rand_penta_matrix_dense( + ... n_rows=n_rows, + ... seed=seed, + ... ill_conditioned=True, + ... ) + >>> # its bandwidth is computed and should be equal to 2 + >>> bandwidth(mat) + (2, 2) + >>> # its condition number is computed and its value should be close to the + >>> # reciprocal floating point precision, i.e., ~1e16 + >>> np.linalg.cond(mat) + 5.058722571393928e+17 + + """ + + # first, the fake diagonal matrix is generated whose entries are strictly + # positive and sorted in descending order + np.random.seed(seed=seed) + d_diag = np.flip(np.sort(np.random.rand(n_rows))) + + # the conditioning is achieved by manipulating the smallest diagonal entry + # Case 1: well-conditioned matrix + if not ill_conditioned: + # here, the smallest diagonal entry is set to a value that is enforced to have + # a minimum magnitude + d_diag = np.maximum(d_diag, _MIN_DIAG_VAL) + + # Case 2: ill-conditioned matrix + else: + # here, the smallest diagonal entry is set to a value that is numerically zero + # compared to the largest entry + d_diag[n_rows - 1] = 0.1 * np.finfo(np.float64).eps * d_diag[0] + + # ... followed by a unit lower triangular matrix with 2 sub-diagonals, but here + # the entries may be negative ... + diagonals = [ + 1.0 - 2.0 * np.random.rand(n_rows - 2), + 1.0 - 2.0 * np.random.rand(n_rows - 1), + np.ones(n_rows), + ] + l_mat = sprs.diags( + diagonals=diagonals, + offsets=[-2, -1, 0], # type: ignore + shape=(n_rows, n_rows), + format="csr", + dtype=np.float64, + ) + + # ... and an upper triangular matrix with 2 super-diagonals + diagonals = [ + np.ones(n_rows), + 1.0 - 2.0 * np.random.rand(n_rows - 1), + 1.0 - 2.0 * np.random.rand(n_rows - 2), + ] + u_mat = sprs.diags( + diagonals=diagonals, + offsets=[0, 1, 2], # type: ignore + shape=(n_rows, n_rows), + format="csr", + dtype=np.float64, + ) + + # finally, the matrix is reconstructed by multiplying the three matrices + return (l_mat.multiply(d_diag[np.newaxis, ::]).dot(u_mat)).toarray() + + # === Doctests === if __name__ == "__main__": # pragma: no cover From 26752f292f569e96b0a7b39d5a4fb95e695b6f72 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 11:44:29 +0200 Subject: [PATCH 10/34] tests: [11] add doctested reference solver; made ill-conditioning more severe --- tests/util_funcs.py | 85 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 72 insertions(+), 13 deletions(-) diff --git a/tests/util_funcs.py b/tests/util_funcs.py index ca8f2be..8ef8efa 100644 --- a/tests/util_funcs.py +++ b/tests/util_funcs.py @@ -9,6 +9,7 @@ from typing import Tuple import numpy as np +from scipy import linalg as spla from scipy import sparse as sprs import pentapy as pp @@ -195,9 +196,6 @@ def gen_conditioned_rand_penta_matrix_dense( Doctests -------- - >>> # Imports - >>> from scipy.linalg import bandwidth - >>> # 1) Generating a super small well-conditioned random pentadiagonal matrix >>> n_rows = 3 >>> seed = 19_031_977 @@ -212,7 +210,7 @@ def gen_conditioned_rand_penta_matrix_dense( [-0.09784268, 0.2270634 , -0.1509019 ], [-0.23431267, 0.00468463, 0.22991003]]) >>> # its bandwidth is computed and should be equal to 2 - >>> bandwidth(mat) + >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and values below 1e10 can be considered good >>> np.linalg.cond(mat) @@ -229,7 +227,7 @@ def gen_conditioned_rand_penta_matrix_dense( [-0.09784268, 0.2270634 , -0.1509019 ], [-0.23431267, 0.00468463, -0.02273771]]) >>> # its bandwidth is computed and should be equal to 2 - >>> bandwidth(mat) + >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and its value should be close to the >>> # reciprocal floating point precision, i.e., ~1e16 @@ -253,7 +251,7 @@ def gen_conditioned_rand_penta_matrix_dense( [ 0. , 0. , 0. , 0.06, -0.14, 0.4 , 0.28], [ 0. , 0. , 0. , 0. , -0.14, 0.36, 0.53]]) >>> # its bandwidth is computed and should be equal to 2 - >>> bandwidth(mat) + >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and values below 1e10 can be considered good >>> np.linalg.cond(mat) @@ -274,7 +272,7 @@ def gen_conditioned_rand_penta_matrix_dense( [ 0. , 0. , 0. , 0.06, -0.14, 0.4 , 0.28], [ 0. , 0. , 0. , 0. , -0.14, 0.36, 0.28]]) >>> # its bandwidth is computed and should be equal to 2 - >>> bandwidth(mat) + >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and its value should be close to the >>> # reciprocal floating point precision, i.e., ~1e16 @@ -290,7 +288,7 @@ def gen_conditioned_rand_penta_matrix_dense( ... ill_conditioned=False, ... ) >>> # its bandwidth is computed and should be equal to 2 - >>> bandwidth(mat) + >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and values below 1e10 can be considered good >>> np.linalg.cond(mat) @@ -303,12 +301,12 @@ def gen_conditioned_rand_penta_matrix_dense( ... ill_conditioned=True, ... ) >>> # its bandwidth is computed and should be equal to 2 - >>> bandwidth(mat) + >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and its value should be close to the >>> # reciprocal floating point precision, i.e., ~1e16 >>> np.linalg.cond(mat) - 5.058722571393928e+17 + 1.7137059583101745e+19 """ @@ -326,9 +324,8 @@ def gen_conditioned_rand_penta_matrix_dense( # Case 2: ill-conditioned matrix else: - # here, the smallest diagonal entry is set to a value that is numerically zero - # compared to the largest entry - d_diag[n_rows - 1] = 0.1 * np.finfo(np.float64).eps * d_diag[0] + # here, the smallest diagonal entry is set to zero + d_diag[n_rows - 1] = 0.0 # ... followed by a unit lower triangular matrix with 2 sub-diagonals, but here # the entries may be negative ... @@ -363,6 +360,68 @@ def gen_conditioned_rand_penta_matrix_dense( return (l_mat.multiply(d_diag[np.newaxis, ::]).dot(u_mat)).toarray() +def solve_penta_matrix_dense_scipy( + mat: np.ndarray, + rhs: np.ndarray, +) -> np.ndarray: + """ + Solves a pentadiagonal matrix system using SciPy's banded solver. + + Doctests + -------- + >>> # Setting up a small test matrix and right-hand side + >>> n_rows = 5 + >>> seed = 19_031_977 + + >>> mat = gen_conditioned_rand_penta_matrix_dense( + ... n_rows=n_rows, + ... seed=seed, + ... ill_conditioned=False, + ... ) + >>> rhs = np.random.rand(n_rows, 5) + + >>> # Solving the system using SciPy's banded solver + >>> sol = solve_penta_matrix_dense_scipy(mat=mat, rhs=rhs) + >>> np.round(sol, 2) + array([[-2.16, -0.36, 0.72, 0.23, -0.2 ], + [ 4.07, 1.3 , 0.81, 1.31, 0.48], + [ 4.05, 0.33, 2.19, 1.22, 0.58], + [-1.9 , -0.79, 1.02, -0.39, 1.02], + [ 6.31, 1.81, 1.29, 1.41, 0.37]]) + + >>> # the solution is checked by verifying that the residual is close to zero + >>> np.max(np.abs(mat @ sol - rhs)) <= np.finfo(np.float64).eps * n_rows + True + + >>> # Setting up a large test matrix and right-hand side + >>> n_rows = 1_000 + + >>> mat = gen_conditioned_rand_penta_matrix_dense( + ... n_rows=n_rows, + ... seed=seed, + ... ill_conditioned=False, + ... ) + >>> rhs = np.random.rand(n_rows, 5) + + >>> # Solving the system using SciPy's banded solver + >>> sol = solve_penta_matrix_dense_scipy(mat=mat, rhs=rhs) + >>> # the solution is checked by verifying that the residual is close to zero + >>> np.max(np.abs(mat @ sol - rhs)) <= np.finfo(np.float64).eps * n_rows + True + + """ + + # first, the matrix is converted to LAPACK banded storage format + mat_banded = pp.create_banded(mat=mat, col_wise=True) + + # then, the system is solved using SciPy's banded solver + return spla.solve_banded( + l_and_u=(2, 2), + ab=mat_banded, + b=rhs, + ) + + # === Doctests === if __name__ == "__main__": # pragma: no cover From a0c8849e05d96236e1a234af00095def29848efb Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 12:05:50 +0200 Subject: [PATCH 11/34] feat: [11] made error messages informative --- src/pentapy/tools.py | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/src/pentapy/tools.py b/src/pentapy/tools.py index 3db7126..ca80d27 100644 --- a/src/pentapy/tools.py +++ b/src/pentapy/tools.py @@ -172,10 +172,14 @@ def create_banded(mat, up=2, low=2, col_wise=True, dtype=None): """ mat = np.asanyarray(mat) if mat.ndim != 2: - msg = "create_banded: matrix has to be 2D" + msg = f"create_banded: matrix has to be 2D, got {mat.ndim}D" raise ValueError(msg) + if mat.shape[0] != mat.shape[1]: - msg = "create_banded: matrix has to be n x n" + msg = ( + f"create_banded: matrix has to be n x n, " + f"got {mat.shape[0]} x {mat.shape[1]}" + ) raise ValueError(msg) size = mat.shape[0] @@ -246,14 +250,23 @@ def create_full(mat, up=2, low=2, col_wise=True): """ mat = np.asanyarray(mat) if mat.ndim != 2: - msg = "create_full: matrix has to be 2D" + msg = f"create_full: matrix has to be 2D, got {mat.ndim}D" raise ValueError(msg) + if mat.shape[0] != up + low + 1: - msg = "create_full: matrix has wrong count of bands" + msg = ( + f"create_full: matrix has wrong count of bands, required " + f"{up} + {low} + 1 = {up + low + 1}, got {mat.shape[0]} bands" + ) raise ValueError(msg) + if mat.shape[1] < max(up, low) + 1: - msg = "create_full: matrix has to few information" + msg = ( + f"create_full: matrix has to few information, required " + f"{max(up, low) + 1} columns, got {mat.shape[1]} columns" + ) raise ValueError(msg) + size = mat.shape[1] mat_full = np.diag(mat[up]) if col_wise: @@ -266,16 +279,17 @@ def create_full(mat, up=2, low=2, col_wise=True): mat_full[diag_indices(size, up - i)] = mat[i, : -(up - i)] for i in range(low): mat_full[diag_indices(size, -(low - i))] = mat[-i - 1, (low - i) :] + return mat_full def _check_penta(mat): if mat.ndim != 2: - msg = "pentapy: matrix has to be 2D" + msg = f"pentapy: matrix has to be 2D, got {mat.ndim}D" raise ValueError(msg) if mat.shape[0] != 5: - msg = "pentapy: matrix needs 5 bands" + msg = f"pentapy: matrix needs 5 bands, got {mat.shape[0]} bands" raise ValueError(msg) if mat.shape[1] < 3: - msg = "pentapy: matrix needs at least 3 rows" + msg = f"pentapy: matrix needs at least 3 rows, got {mat.shape[1]} rows" raise ValueError(msg) From ee6943a2fae9d3fb676bc56e20c90f41b2b4d29b Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 13:23:46 +0200 Subject: [PATCH 12/34] fix: [23] disabled cdivision to fix error handling on Python side --- src/pentapy/solver.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pentapy/solver.pyx b/src/pentapy/solver.pyx index c8a22a2..091288c 100644 --- a/src/pentapy/solver.pyx +++ b/src/pentapy/solver.pyx @@ -1,4 +1,4 @@ -# cython: language_level=3, boundscheck=True, wraparound=False, cdivision=True +# cython: language_level=3, boundscheck=True, wraparound=False, cdivision=False """ This is a solver linear equation systems with a penta-diagonal matrix, From 028484372fcce5d445695cae86f9db884bbe4290 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 13:28:14 +0200 Subject: [PATCH 13/34] feat/refactor: [11] enabled multipe right-hand sides for solver I; improved import chain; improved code readability --- src/pentapy/core.py | 58 ++++++++++++++++++++++++++++++--------------- 1 file changed, 39 insertions(+), 19 deletions(-) diff --git a/src/pentapy/core.py b/src/pentapy/core.py index 9a55780..2393122 100644 --- a/src/pentapy/core.py +++ b/src/pentapy/core.py @@ -5,8 +5,8 @@ import numpy as np +from pentapy import tools as ptools from pentapy.solver import penta_solver1, penta_solver2 -from pentapy.tools import _check_penta, create_banded, shift_banded def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): @@ -66,40 +66,60 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): result : :class:`numpy.ndarray` Solution of the equation system """ + if solver in [1, "1", "PTRANS-I"]: if is_flat and index_row_wise: mat_flat = np.asarray(mat, dtype=np.double) - _check_penta(mat_flat) + ptools._check_penta(mat_flat) elif is_flat: mat_flat = np.array(mat, dtype=np.double) - _check_penta(mat_flat) - shift_banded(mat_flat, copy=False) + ptools._check_penta(mat_flat) + ptools.shift_banded(mat_flat, copy=False) else: - mat_flat = create_banded(mat, col_wise=False, dtype=np.double) + mat_flat = ptools.create_banded(mat, col_wise=False, dtype=np.double) + rhs = np.asarray(rhs, dtype=np.double) + + # Special case: Early exit when the matrix has only 3 rows/columns + # NOTE: this avoids memory leakage in the Cython-solver that will iterate over + # at least 4 rows/columns no matter what + if mat_flat.shape[1] == 3: + return np.linalg.solve( + a=ptools.create_full(mat_flat, col_wise=False), + b=rhs, + ) + + # if there is only a single right-hand side, it has to be reshaped to a 2D array + # NOTE: this has to be reverted at the end single_rhs = rhs.ndim == 1 + rhs_og_shape = rhs.shape if single_rhs: rhs = rhs[:, np.newaxis] try: + # if there was only a 1D right-hand side, the result has to be flattened if single_rhs: return penta_solver1(mat_flat, rhs).ravel() return penta_solver1(mat_flat, rhs) + except ZeroDivisionError: warnings.warn("pentapy: PTRANS-I not suitable for input-matrix.") - return np.full_like(rhs, np.nan) + return np.full(shape=rhs_og_shape, fill_value=np.nan) + elif solver in [2, "2", "PTRANS-II"]: if is_flat and index_row_wise: mat_flat = np.asarray(mat, dtype=np.double) - _check_penta(mat_flat) + ptools._check_penta(mat_flat) elif is_flat: mat_flat = np.array(mat, dtype=np.double) - _check_penta(mat_flat) - shift_banded(mat_flat, copy=False) + ptools._check_penta(mat_flat) + ptools.shift_banded(mat_flat, copy=False) else: - mat_flat = create_banded(mat, col_wise=False, dtype=np.double) + mat_flat = ptools.create_banded(mat, col_wise=False, dtype=np.double) + rhs = np.asarray(rhs, dtype=np.double) + try: return penta_solver2(mat_flat, rhs) except ZeroDivisionError: @@ -113,12 +133,12 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): raise ValueError(msg) from imp_err if is_flat and index_row_wise: mat_flat = np.array(mat) - _check_penta(mat_flat) - shift_banded(mat_flat, col_to_row=False, copy=False) + ptools._check_penta(mat_flat) + ptools.shift_banded(mat_flat, col_to_row=False, copy=False) elif is_flat: mat_flat = np.asarray(mat) else: - mat_flat = create_banded(mat) + mat_flat = ptools.create_banded(mat) return solve_banded((2, 2), mat_flat, rhs) elif solver in [4, "4", "spsolve"]: # pragma: no cover try: @@ -129,12 +149,12 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): raise ValueError(msg) from imp_err if is_flat and index_row_wise: mat_flat = np.array(mat) - _check_penta(mat_flat) - shift_banded(mat_flat, col_to_row=False, copy=False) + ptools._check_penta(mat_flat) + ptools.shift_banded(mat_flat, col_to_row=False, copy=False) elif is_flat: mat_flat = np.asarray(mat) else: - mat_flat = create_banded(mat) + mat_flat = ptools.create_banded(mat) size = mat_flat.shape[1] M = sps.spdiags(mat_flat, [2, 1, 0, -1, -2], size, size, format="csc") return spsolve(M, rhs, use_umfpack=False) @@ -153,12 +173,12 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): raise ValueError(msg) from imp_err if is_flat and index_row_wise: mat_flat = np.array(mat) - _check_penta(mat_flat) - shift_banded(mat_flat, col_to_row=False, copy=False) + ptools._check_penta(mat_flat) + ptools.shift_banded(mat_flat, col_to_row=False, copy=False) elif is_flat: mat_flat = np.asarray(mat) else: - mat_flat = create_banded(mat) + mat_flat = ptools.create_banded(mat) size = mat_flat.shape[1] M = sps.spdiags(mat_flat, [2, 1, 0, -1, -2], size, size, format="csc") return spsolve(M, rhs, use_umfpack=True) From c32bec2744b08d455627197a56eca6927bfa1a0b Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 13:28:40 +0200 Subject: [PATCH 14/34] tests: [11] added shape check to doctest --- tests/util_funcs.py | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/tests/util_funcs.py b/tests/util_funcs.py index 8ef8efa..be9f3c1 100644 --- a/tests/util_funcs.py +++ b/tests/util_funcs.py @@ -209,7 +209,9 @@ def gen_conditioned_rand_penta_matrix_dense( array([[ 0.92453713, 0.28308514, -0.09972199], [-0.09784268, 0.2270634 , -0.1509019 ], [-0.23431267, 0.00468463, 0.22991003]]) - >>> # its bandwidth is computed and should be equal to 2 + >>> # it has to be square and its bandwidth is computed and should be equal to 2 + >>> mat.shape[0] == mat.shape[1] + True >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and values below 1e10 can be considered good @@ -226,7 +228,9 @@ def gen_conditioned_rand_penta_matrix_dense( array([[ 0.92453713, 0.28308514, -0.09972199], [-0.09784268, 0.2270634 , -0.1509019 ], [-0.23431267, 0.00468463, -0.02273771]]) - >>> # its bandwidth is computed and should be equal to 2 + >>> # it has to be square and its bandwidth is computed and should be equal to 2 + >>> mat.shape[0] == mat.shape[1] + True >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and its value should be close to the @@ -250,7 +254,9 @@ def gen_conditioned_rand_penta_matrix_dense( [ 0. , 0. , 0.43, 0.13, 0.39, -0.1 , -0.15], [ 0. , 0. , 0. , 0.06, -0.14, 0.4 , 0.28], [ 0. , 0. , 0. , 0. , -0.14, 0.36, 0.53]]) - >>> # its bandwidth is computed and should be equal to 2 + >>> # it has to be square and its bandwidth is computed and should be equal to 2 + >>> mat.shape[0] == mat.shape[1] + True >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and values below 1e10 can be considered good @@ -271,7 +277,9 @@ def gen_conditioned_rand_penta_matrix_dense( [ 0. , 0. , 0.43, 0.13, 0.39, -0.1 , -0.15], [ 0. , 0. , 0. , 0.06, -0.14, 0.4 , 0.28], [ 0. , 0. , 0. , 0. , -0.14, 0.36, 0.28]]) - >>> # its bandwidth is computed and should be equal to 2 + >>> # it has to be square and its bandwidth is computed and should be equal to 2 + >>> mat.shape[0] == mat.shape[1] + True >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and its value should be close to the @@ -287,7 +295,9 @@ def gen_conditioned_rand_penta_matrix_dense( ... seed=seed, ... ill_conditioned=False, ... ) - >>> # its bandwidth is computed and should be equal to 2 + >>> # it has to be square and its bandwidth is computed and should be equal to 2 + >>> mat.shape[0] == mat.shape[1] + True >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and values below 1e10 can be considered good @@ -300,7 +310,9 @@ def gen_conditioned_rand_penta_matrix_dense( ... seed=seed, ... ill_conditioned=True, ... ) - >>> # its bandwidth is computed and should be equal to 2 + >>> # it has to be square and its bandwidth is computed and should be equal to 2 + >>> mat.shape[0] == mat.shape[1] + True >>> spla.bandwidth(mat) (2, 2) >>> # its condition number is computed and its value should be close to the From e7d92c86e67ad2cc77a9e293e9d35c46f81ac624 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 13:29:57 +0200 Subject: [PATCH 15/34] tests: [11] added extensive parametrized tests for solver I that also cover the edge cases --- tests/test_solver_1.py | 146 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 146 insertions(+) create mode 100644 tests/test_solver_1.py diff --git a/tests/test_solver_1.py b/tests/test_solver_1.py new file mode 100644 index 0000000..df00f9d --- /dev/null +++ b/tests/test_solver_1.py @@ -0,0 +1,146 @@ +""" +Test suite for testing the pentadiagonal solver based on Algorithm PTRANS-I. + +""" + +# === Imports === + +from typing import Literal + +import numpy as np +import pentapy as pp +import pytest +import util_funcs as uf + +# === Constants === + +SEED = 19_031_977 +N_ROWS = [ + 3, + 4, + 5, + 10, + 11, + 25, + 26, + 50, + 51, + 100, + 101, + 250, + 251, + 500, + 501, + 1_000, + 1_001, + 10_000, + 10_001, +] +REF_WARNING = "pentapy: PTRANS-I not suitable for input-matrix." + +# === Tests === + + +@pytest.mark.parametrize("induce_error", [False, True]) +@pytest.mark.parametrize("solver_alias", [1]) # "1", "PTRANS-I"]) +@pytest.mark.parametrize("input_layout", ["full", "banded_row_wise", "banded_col_wise"]) +@pytest.mark.parametrize("n_rhs", [None, 1, 10]) +@pytest.mark.parametrize("n_rows", N_ROWS) +def test_penta_solver1( + n_rows: int, + n_rhs: int, + input_layout: Literal["full", "banded_row_wise", "banded_col_wise"], + solver_alias: Literal[1, "1", "PTRANS-I"], + induce_error: bool, +) -> None: + """ + Tests the pentadiagonal solver based on Algorithm PTRANS-I when starting from + different input layouts, number of right-hand sides, number of rows, and also + when inducing an error by making the first diagonal element zero. + It has to be ensured that the edge case of ``n_rows = 3`` is also covered. + + """ + + # first, a random pentadiagonal matrix is generated + mat_full = uf.gen_conditioned_rand_penta_matrix_dense( + n_rows=n_rows, + seed=SEED, + ill_conditioned=False, + ) + + # an error is induced by setting the first diagonal element to zero + if induce_error: + # the induction of the error is only possible if the matrix does not have + # only 3 rows + if n_rows == 3: + pytest.skip( + "Only 3 rows, cannot induce error because this will not go into " + "PTRANS-I, but NumPy" + ) + + mat_full[0, 0] = 0.0 + + # the right-hand side is generated + np.random.seed(SEED) + if n_rhs is not None: + rhs = np.random.rand(n_rows, n_rhs) + result_shape = (n_rows, n_rhs) + else: + rhs = np.random.rand(n_rows) + result_shape = (n_rows,) + + # the matrix is converted to the desired layout + if input_layout == "full": + mat = mat_full + kwargs = dict(is_flat=False) + + elif input_layout == "banded_row_wise": + mat = pp.create_banded(mat_full, col_wise=False) + kwargs = dict( + is_flat=True, + index_row_wise=True, + ) + + elif input_layout == "banded_col_wise": + mat = pp.create_banded(mat_full, col_wise=True) + kwargs = dict( + is_flat=True, + index_row_wise=False, + ) + + else: + raise ValueError(f"Invalid input layout: {input_layout}") + + # the solution is computed + # Case 1: in case of an error, a warning has to be issued and the result has to + # be NaN + if induce_error: + with pytest.warns(UserWarning, match=REF_WARNING): + sol = pp.solve( + mat=mat, + rhs=rhs, + solver=solver_alias, # type: ignore + **kwargs, + ) + assert sol.shape == result_shape + assert np.isnan(sol).all() + + return + + # Case 2: in case of no error, the solution can be computed without any issues + sol = pp.solve( + mat=mat, + rhs=rhs, + solver=solver_alias, # type: ignore + **kwargs, + ) + assert sol.shape == result_shape + + # if no error was induced, the reference solution is computed with SciPy + sol_ref = uf.solve_penta_matrix_dense_scipy( + mat=mat_full, + rhs=rhs, + ) + + # the solutions are compared + assert np.allclose(sol, sol_ref) From 647012f05677306e911ef0edcfa1162785efbafe Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 14:11:28 +0200 Subject: [PATCH 16/34] tests: [11] added more intermediate sizes for tests --- tests/test_solver_1.py | 4 ++++ tests/test_tools.py | 7 ++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/tests/test_solver_1.py b/tests/test_solver_1.py index df00f9d..78956f6 100644 --- a/tests/test_solver_1.py +++ b/tests/test_solver_1.py @@ -33,6 +33,10 @@ 501, 1_000, 1_001, + 2500, + 2501, + 5_000, + 5_001, 10_000, 10_001, ] diff --git a/tests/test_tools.py b/tests/test_tools.py index 2f54c48..cabac61 100644 --- a/tests/test_tools.py +++ b/tests/test_tools.py @@ -9,9 +9,10 @@ from typing import Optional, Tuple, Type import numpy as np -import pentapy as pp import pytest import util_funcs as uf + +import pentapy as pp from pentapy.tools import _check_penta warnings.simplefilter("always") @@ -37,6 +38,10 @@ 501, 1_000, 1_001, + 2500, + 2501, + 5_000, + 5_001, 10_000, 10_001, ] From 8ab514f9404c803b4f6f368a5c21f607729d4201 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 15:51:36 +0200 Subject: [PATCH 17/34] feat/fix: [11] added cython annotations to build process; fixed wrong f-string; fixed typo --- setup.py | 9 ++++++--- src/pentapy/.gitignore | 2 ++ 2 files changed, 8 insertions(+), 3 deletions(-) create mode 100644 src/pentapy/.gitignore diff --git a/setup.py b/setup.py index 8081e05..fc8648c 100644 --- a/setup.py +++ b/setup.py @@ -1,15 +1,18 @@ -"""pentapy: A toolbox for pentadiagonal matrizes.""" +"""pentapy: A toolbox for pentadiagonal matrices.""" import os +import Cython.Compiler.Options import numpy as np from Cython.Build import cythonize from setuptools import Extension, setup +Cython.Compiler.Options.annotate = True + # cython extensions CY_MODULES = [ Extension( - name=f"pentapy.solver", + name="pentapy.solver", sources=[os.path.join("src", "pentapy", "solver.pyx")], include_dirs=[np.get_include()], define_macros=[("NPY_NO_DEPRECATED_API", "NPY_1_7_API_VERSION")], @@ -17,7 +20,7 @@ ] setup( - ext_modules=cythonize(CY_MODULES), + ext_modules=cythonize(CY_MODULES, nthreads=1, annotate=True), package_data={"pentapy": ["*.pxd"]}, # include pxd files include_package_data=False, # ignore other files zip_safe=False, diff --git a/src/pentapy/.gitignore b/src/pentapy/.gitignore new file mode 100644 index 0000000..53cc0d6 --- /dev/null +++ b/src/pentapy/.gitignore @@ -0,0 +1,2 @@ +# Cython html files +*.html \ No newline at end of file From 5086959a21a974f697aa3212ecfb22d8671248d4 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 15:52:50 +0200 Subject: [PATCH 18/34] refactor/doc: fixed missing variable declarations; added clarifying comments --- src/pentapy/solver.pyx | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/pentapy/solver.pyx b/src/pentapy/solver.pyx index 091288c..9d45c90 100644 --- a/src/pentapy/solver.pyx +++ b/src/pentapy/solver.pyx @@ -73,10 +73,11 @@ cdef void c_penta_factorize_algo1( """ Factorizes the pentadiagonal matrix ``A`` into - - auxiliary coefficients ``e``, ``mu`` and ``gamma`` for the transformation of the - right-hand side - - a unit upper triangular matrix with the main diagonals ``alpha`` and ``beta`` - for the following backward substitution. Its unit main diagonal is implicit. + - auxiliary coefficients ``e``, ``mu`` and ``gamma`` (``ga``) for the transformation + of the right-hand side + - a unit upper triangular matrix with the main diagonals ``alpha``(``al``) and + ``beta`` (``be``) for the following backward substitution. Its unit main + diagonal is implicit. They are overwriting the memoryview ``mat_factorized`` as follows: @@ -99,8 +100,9 @@ cdef void c_penta_factorize_algo1( # === Variable declarations === cdef uint64_t iter_row - cdef double mu_i, ga_i, e_i - cdef double al_i, al_i_minus_1, al_i_plus_1 + cdef double mu_i, ga_i, e_i # mu, gamma, e + cdef double al_i, al_i_minus_1, al_i_plus_1 # alpha + cdef double be_i, be_i_minus_1, be_i_plus_1 # beta # === Factorization === From f062b886c3b8218460f651fc742f75a6df305744 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 17:26:08 +0200 Subject: [PATCH 19/34] style: [11] made - signs better readable --- src/pentapy/solver.pyx | 64 +++++++++++++++++++++--------------------- 1 file changed, 32 insertions(+), 32 deletions(-) diff --git a/src/pentapy/solver.pyx b/src/pentapy/solver.pyx index 9d45c90..ea15048 100644 --- a/src/pentapy/solver.pyx +++ b/src/pentapy/solver.pyx @@ -130,7 +130,7 @@ cdef void c_penta_factorize_algo1( mat_factorized[1, 4] = be_i # Central rows - for iter_row in range(2, mat_n_rows-2): + for iter_row in range(2, mat_n_rows - 2): e_i = mat_flat[4, iter_row] ga_i = mat_flat[3, iter_row] - al_i_minus_1 * e_i mu_i = mat_flat[2, iter_row] - be_i_minus_1 * e_i - al_i * ga_i @@ -150,27 +150,27 @@ cdef void c_penta_factorize_algo1( mat_factorized[iter_row, 4] = be_i # Second to last row - e_i = mat_flat[4, mat_n_rows-2] - ga_i = mat_flat[3, mat_n_rows-2] - al_i_minus_1 * e_i - mu_i = mat_flat[2, mat_n_rows-2] - be_i_minus_1 * e_i - al_i * ga_i - al_i_plus_1 = (mat_flat[1, mat_n_rows-2] - be_i * ga_i) / mu_i + e_i = mat_flat[4, mat_n_rows - 2] + ga_i = mat_flat[3, mat_n_rows - 2] - al_i_minus_1 * e_i + mu_i = mat_flat[2, mat_n_rows - 2] - be_i_minus_1 * e_i - al_i * ga_i + al_i_plus_1 = (mat_flat[1, mat_n_rows - 2] - be_i * ga_i) / mu_i - mat_factorized[mat_n_rows-2, 0] = e_i - mat_factorized[mat_n_rows-2, 1] = mu_i - mat_factorized[mat_n_rows-2, 2] = ga_i - mat_factorized[mat_n_rows-2, 3] = al_i_plus_1 - mat_factorized[mat_n_rows-2, 4] = 0.0 + mat_factorized[mat_n_rows - 2, 0] = e_i + mat_factorized[mat_n_rows - 2, 1] = mu_i + mat_factorized[mat_n_rows - 2, 2] = ga_i + mat_factorized[mat_n_rows - 2, 3] = al_i_plus_1 + mat_factorized[mat_n_rows - 2, 4] = 0.0 # Last Row - e_i = mat_flat[4, mat_n_rows-1] - ga_i = mat_flat[3, mat_n_rows-1] - al_i * e_i - mu_i = mat_flat[2, mat_n_rows-1] - be_i * e_i - al_i_plus_1 * ga_i + e_i = mat_flat[4, mat_n_rows - 1] + ga_i = mat_flat[3, mat_n_rows - 1] - al_i * e_i + mu_i = mat_flat[2, mat_n_rows - 1] - be_i * e_i - al_i_plus_1 * ga_i - mat_factorized[mat_n_rows-1, 0] = e_i - mat_factorized[mat_n_rows-1, 1] = mu_i - mat_factorized[mat_n_rows-1, 2] = ga_i - mat_factorized[mat_n_rows-1, 3] = 0.0 - mat_factorized[mat_n_rows-1, 4] = 0.0 + mat_factorized[mat_n_rows - 1, 0] = e_i + mat_factorized[mat_n_rows - 1, 1] = mu_i + mat_factorized[mat_n_rows - 1, 2] = ga_i + mat_factorized[mat_n_rows - 1, 3] = 0.0 + mat_factorized[mat_n_rows - 1, 4] = 0.0 return @@ -207,7 +207,7 @@ cdef void c_solve_penta_from_factorize_algo_1( result_view[1] = ze_i # Central rows - for iter_row in range(2, mat_n_rows-2): + for iter_row in range(2, mat_n_rows - 2): ze_i_plus_1 = ( rhs_single[iter_row] - ze_i_minus_1 * mat_factorized[iter_row, 0] @@ -219,30 +219,30 @@ cdef void c_solve_penta_from_factorize_algo_1( # Second to last row ze_i_plus_1 = ( - rhs_single[mat_n_rows-2] - - ze_i_minus_1 * mat_factorized[mat_n_rows-2, 0] - - ze_i * mat_factorized[mat_n_rows-2, 2] - ) / mat_factorized[mat_n_rows-2, 1] + rhs_single[mat_n_rows - 2] + - ze_i_minus_1 * mat_factorized[mat_n_rows - 2, 0] + - ze_i * mat_factorized[mat_n_rows - 2, 2] + ) / mat_factorized[mat_n_rows - 2, 1] ze_i_minus_1 = ze_i ze_i = ze_i_plus_1 - result_view[mat_n_rows-2] = ze_i_plus_1 + result_view[mat_n_rows - 2] = ze_i_plus_1 # Last row ze_i_plus_1 = ( - rhs_single[mat_n_rows-1] - - ze_i_minus_1 * mat_factorized[mat_n_rows-1, 0] - - ze_i * mat_factorized[mat_n_rows-1, 2] - ) / mat_factorized[mat_n_rows-1, 1] - result_view[mat_n_rows-1] = ze_i_plus_1 + rhs_single[mat_n_rows - 1] + - ze_i_minus_1 * mat_factorized[mat_n_rows - 1, 0] + - ze_i * mat_factorized[mat_n_rows - 1, 2] + ) / mat_factorized[mat_n_rows - 1, 1] + result_view[mat_n_rows - 1] = ze_i_plus_1 # === Backward substitution === # The solution vector is calculated by backward substitution that overwrites the # right-hand side vector with the solution vector - ze_i -= mat_factorized[mat_n_rows-2, 3] * ze_i_plus_1 - result_view[mat_n_rows-2] = ze_i + ze_i -= mat_factorized[mat_n_rows - 2, 3] * ze_i_plus_1 + result_view[mat_n_rows - 2] = ze_i - for iter_row in range(mat_n_rows-3, -1, -1): + for iter_row in range(mat_n_rows - 3, -1, -1): result_view[iter_row] -= ( mat_factorized[iter_row, 3] * ze_i + mat_factorized[iter_row, 4] * ze_i_plus_1 From 15c787e19c2396f18eb16d4a7717d5739b516632 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 18:56:44 +0200 Subject: [PATCH 20/34] doc/refactor: [11] improved docs and comments of algorithm I; removed `uint` --- src/pentapy/solver.pyx | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/src/pentapy/solver.pyx b/src/pentapy/solver.pyx index ea15048..451eee9 100644 --- a/src/pentapy/solver.pyx +++ b/src/pentapy/solver.pyx @@ -11,7 +11,7 @@ implemented in Cython. import numpy as np cimport numpy as np -from libc.stdint cimport int64_t, uint64_t +from libc.stdint cimport int64_t # === Main Python Interface === @@ -40,12 +40,16 @@ cdef double[:, :] c_penta_solver1(double[:, :] mat_flat, double[:, :] rhs): """ - cdef uint64_t mat_n_rows = mat_flat.shape[1] - cdef uint64_t rhs_n_cols = rhs.shape[1] - cdef uint64_t iter_col + # === Variable declarations === + + cdef int64_t mat_n_rows = mat_flat.shape[1] + cdef int64_t rhs_n_cols = rhs.shape[1] + cdef int64_t iter_col cdef double[::, ::1] result = np.empty(shape=(mat_n_rows, rhs_n_cols)) cdef double[::, ::1] mat_factorized = np.empty(shape=(mat_n_rows, 5)) + # === Solving the system of equations === + # first, the matrix is factorized c_penta_factorize_algo1( mat_flat, @@ -67,7 +71,7 @@ cdef double[:, :] c_penta_solver1(double[:, :] mat_flat, double[:, :] rhs): cdef void c_penta_factorize_algo1( double[:, :] mat_flat, - uint64_t mat_n_rows, + int64_t mat_n_rows, double[::, ::1] mat_factorized, ): """ @@ -82,14 +86,14 @@ cdef void c_penta_factorize_algo1( They are overwriting the memoryview ``mat_factorized`` as follows: ```bash - [[ * mu_0 * al_0 be_0 ] - [ * mu_1 ga_1 al_1 be_1 ] - [ e_2 mu_2 ga_2 al_2 be_2 ] + [[ * mu_0 * al_0 be_0 ] + [ * mu_1 ga_1 al_1 be_1 ] + [ e_2 mu_2 ga_2 al_2 be_2 ] ... [ e_i mu_i ga_i al_i be_i ] - ... - [ e_{n-2} mu_{n-2} ga_{n-2} al_{n-2} * ] - [ e_{n-1} mu_{n-1} ga_{n-1} * * ]] + [ e_{n-3} mu_{n-3} ga_{n-3} al_{n-3} be_{n-3} ] ... + [ e_{n-2} mu_{n-2} ga_{n-2} al_{n-2} * ] + [ e_{n-1} mu_{n-1} ga_{n-1} * * ]] ``` where the entries marked with ``*`` are not used by design, but overwritten with @@ -99,7 +103,7 @@ cdef void c_penta_factorize_algo1( # === Variable declarations === - cdef uint64_t iter_row + cdef int64_t iter_row cdef double mu_i, ga_i, e_i # mu, gamma, e cdef double al_i, al_i_minus_1, al_i_plus_1 # alpha cdef double be_i, be_i_minus_1, be_i_plus_1 # beta @@ -176,7 +180,7 @@ cdef void c_penta_factorize_algo1( cdef void c_solve_penta_from_factorize_algo_1( - uint64_t mat_n_rows, + int64_t mat_n_rows, double[::, ::1] mat_factorized, double[::] rhs_single, double[::] result_view, From ef7ec51a6eceda76a54c39b83b912baeafe2f4bf Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 19:45:01 +0200 Subject: [PATCH 21/34] wip: [11] restructured chaotic alias comparisons --- src/pentapy/_models.py | 49 +++++++++++++++++ src/pentapy/core.py | 122 ++++++++++++++++++++++++----------------- 2 files changed, 121 insertions(+), 50 deletions(-) create mode 100644 src/pentapy/_models.py diff --git a/src/pentapy/_models.py b/src/pentapy/_models.py new file mode 100644 index 0000000..c75eb8c --- /dev/null +++ b/src/pentapy/_models.py @@ -0,0 +1,49 @@ +""" +Auxiliary models for the pentapy package. + +""" + +# === Imports === + +from enum import IntEnum +from typing import Dict + +# === Models === + + +class PentaSolverAliases(IntEnum): + """ + Defines all available solver aliases for pentadiagonal systems, namely + + - ``PTRANS_I``: The PTRANS-I algorithm + - ``PTRANS_II``: The PTRANS-II algorithm + - ``LAPACK``: Scipy's LAPACK solver :func:`scipy.linalg.solve_banded` + - ``SUPER_LU``: Scipy's SuperLU solver :func:`scipy.sparse.linalg.spsolve(..., use_umfpack=False)` + - ``UMFPACK``: Scipy's UMFpack solver :func:`scipy.sparse.linalg.spsolve(..., use_umfpack=True)` + + """ # noqa: E501 + + PTRANS_I = 1 + PTRANS_II = 2 + LAPACK = 3 + SUPER_LU = 4 + UMFPACK = 5 + + +# === Constants === + +_SOLVER_ALIAS_CONVERSIONS: Dict[str, PentaSolverAliases] = { + "1": PentaSolverAliases.PTRANS_I, + "ptrans-i": PentaSolverAliases.PTRANS_I, + "2": PentaSolverAliases.PTRANS_II, + "ptrans-ii": PentaSolverAliases.PTRANS_II, + "3": PentaSolverAliases.LAPACK, + "lapack": PentaSolverAliases.LAPACK, + "solve_banded": PentaSolverAliases.LAPACK, + "4": PentaSolverAliases.SUPER_LU, + "spsolve": PentaSolverAliases.SUPER_LU, + "5": PentaSolverAliases.UMFPACK, + "spsolve_umf": PentaSolverAliases.UMFPACK, + "umf": PentaSolverAliases.UMFPACK, + "umf_pack": PentaSolverAliases.UMFPACK, +} diff --git a/src/pentapy/core.py b/src/pentapy/core.py index 2393122..ccf0c8f 100644 --- a/src/pentapy/core.py +++ b/src/pentapy/core.py @@ -1,15 +1,48 @@ """The core module of pentapy.""" # pylint: disable=C0103, C0415, R0911, E0611 + +# === Imports === + import warnings +from typing import Literal import numpy as np +from pentapy import _models as pmodels +from pentapy import solver as psolver # type: ignore from pentapy import tools as ptools -from pentapy.solver import penta_solver1, penta_solver2 - -def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): +# === Solver === + + +def solve( + mat: np.ndarray, + rhs: np.ndarray, + is_flat: bool = False, + index_row_wise: bool = True, + solver: Literal[ + 1, + "1", + "PTRANS-I", + "ptrans-i", + 2, + "2", + "PTRANS-II", + "ptrans-ii", + 3, + "3", + "lapack", + 4, + "4", + "spsolve", + 5, + "5", + "spsolve_umf", + "umf", + "umf_pack", + ] = 1, +) -> np.ndarray: """ Solver for a pentadiagonal system. @@ -39,35 +72,39 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): Parameters ---------- - mat : :class:`numpy.ndarray` - The Matrix or the flattened Version of the pentadiagonal matrix. - rhs : :class:`numpy.ndarray` - The right hand side of the equation system. - is_flat : :class:`bool`, optional + mat : :class:`numpy.ndarray` of shape (m, m) or (5, m) + The full or flattened version of the pentadiagonal matrix. + rhs : :class:`numpy.ndarray` of shape (m,) or (m, n) + The right hand side(s) of the equation system. Its shape is preserved. + is_flat : :class:`bool`, default=False State if the matrix is already flattend. Default: ``False`` - index_row_wise : :class:`bool`, optional + index_row_wise : :class:`bool`, default=True State if the flattend matrix is row-wise flattend. Default: ``True`` - solver : :class:`int` or :class:`str`, optional + solver : :class:`int` or :class:`str`, default=1 Which solver should be used. The following are provided: - * ``[1, "1", "PTRANS-I"]`` : The PTRANS-I algorithm + * ``[1, "1", "PTRANS-I"]`` : The PTRANS-I algorithm (default) * ``[2, "2", "PTRANS-II"]`` : The PTRANS-II algorithm - * ``[3, "3", "lapack", "solve_banded"]`` : - scipy.linalg.solve_banded - * ``[4, "4", "spsolve"]`` : - The scipy sparse solver without umf_pack - * ``[5, "5", "spsolve_umf", "umf", "umf_pack"]`` : - The scipy sparse solver with umf_pack + * ``[3, "3", "lapack", "solve_banded"]`` : :func:`scipy.linalg.solve_banded` + * ``[4, "4", "spsolve"]`` : :func:`scipy.sparse.linalg.spsolve(..., use_umfpack=False)` + * ``[5, "5", "spsolve_umf", "umf", "umf_pack"]`` : :func:`scipy.sparse.linalg.spsolve(..., use_umfpack=False)` - Default: ``1`` + Strings are not case-sensitive. Returns ------- - result : :class:`numpy.ndarray` - Solution of the equation system + result : :class:`numpy.ndarray` of shape (m,) or (m, n) + Solution of the equation system with the same shape as ``rhs``. + """ - if solver in [1, "1", "PTRANS-I"]: + # first, the solver is converted to the internal name to avoid confusion + solver_inter = pmodels._SOLVER_ALIAS_CONVERSIONS[str(solver).lower()] + + if solver_inter in { + pmodels.PentaSolverAliases.PTRANS_I, + pmodels.PentaSolverAliases.PTRANS_II, + }: if is_flat and index_row_wise: mat_flat = np.asarray(mat, dtype=np.double) ptools._check_penta(mat_flat) @@ -97,35 +134,23 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): rhs = rhs[:, np.newaxis] try: + solver_func = ( + psolver.penta_solver1 + if solver_inter == pmodels.PentaSolverAliases.PTRANS_I + else psolver.penta_solver2 + ) + # if there was only a 1D right-hand side, the result has to be flattened if single_rhs: - return penta_solver1(mat_flat, rhs).ravel() + return solver_func(mat_flat, rhs).ravel() - return penta_solver1(mat_flat, rhs) + return solver_func(mat_flat, rhs) except ZeroDivisionError: warnings.warn("pentapy: PTRANS-I not suitable for input-matrix.") return np.full(shape=rhs_og_shape, fill_value=np.nan) - elif solver in [2, "2", "PTRANS-II"]: - if is_flat and index_row_wise: - mat_flat = np.asarray(mat, dtype=np.double) - ptools._check_penta(mat_flat) - elif is_flat: - mat_flat = np.array(mat, dtype=np.double) - ptools._check_penta(mat_flat) - ptools.shift_banded(mat_flat, copy=False) - else: - mat_flat = ptools.create_banded(mat, col_wise=False, dtype=np.double) - - rhs = np.asarray(rhs, dtype=np.double) - - try: - return penta_solver2(mat_flat, rhs) - except ZeroDivisionError: - warnings.warn("pentapy: PTRANS-II not suitable for input-matrix.") - return np.full_like(rhs, np.nan) - elif solver in [3, "3", "lapack", "solve_banded"]: # pragma: no cover + elif solver_inter == pmodels.PentaSolverAliases.LAPACK: # pragma: no cover try: from scipy.linalg import solve_banded except ImportError as imp_err: # pragma: no cover @@ -140,7 +165,8 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): else: mat_flat = ptools.create_banded(mat) return solve_banded((2, 2), mat_flat, rhs) - elif solver in [4, "4", "spsolve"]: # pragma: no cover + + elif solver_inter == pmodels.PentaSolverAliases.SUPER_LU: # pragma: no cover try: from scipy import sparse as sps from scipy.sparse.linalg import spsolve @@ -158,13 +184,8 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): size = mat_flat.shape[1] M = sps.spdiags(mat_flat, [2, 1, 0, -1, -2], size, size, format="csc") return spsolve(M, rhs, use_umfpack=False) - elif solver in [ - 5, - "5", - "spsolve_umf", - "umf", - "umf_pack", - ]: # pragma: no cover + + elif solver_inter == pmodels.PentaSolverAliases.UMFPACK: # pragma: no cover try: from scipy import sparse as sps from scipy.sparse.linalg import spsolve @@ -182,6 +203,7 @@ def solve(mat, rhs, is_flat=False, index_row_wise=True, solver=1): size = mat_flat.shape[1] M = sps.spdiags(mat_flat, [2, 1, 0, -1, -2], size, size, format="csc") return spsolve(M, rhs, use_umfpack=True) + else: # pragma: no cover msg = f"pentapy.solve: unknown solver ({solver})" raise ValueError(msg) From 94d110ce03f676fe9ba6643d1f5ef69788cf715b Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 20:15:26 +0200 Subject: [PATCH 22/34] feat: [11] finalized multiple right-hand side support (serial) on Cython level; fixed typos --- src/pentapy/solver.pxd | 4 +- src/pentapy/solver.pyx | 262 ++++++++++++++++++++++++++++++++++------- 2 files changed, 222 insertions(+), 44 deletions(-) diff --git a/src/pentapy/solver.pxd b/src/pentapy/solver.pxd index 05d249f..b16f8a0 100644 --- a/src/pentapy/solver.pxd +++ b/src/pentapy/solver.pxd @@ -1,4 +1,4 @@ # cython: language_level=3 -cdef double[:, :] c_penta_solver1(double[:, :] mat_flat, double[:, :] rhs) +cdef double[::, ::] c_penta_solver1(double[::, ::] mat_flat, double[::, ::] rhs) -cdef double[:] c_penta_solver2(double[:, :] mat_flat, double[:] rhs) +cdef double[::, ::] c_penta_solver2(double[::, ::] mat_flat, double[::, ::] rhs) diff --git a/src/pentapy/solver.pyx b/src/pentapy/solver.pyx index 451eee9..1494a25 100644 --- a/src/pentapy/solver.pyx +++ b/src/pentapy/solver.pyx @@ -1,4 +1,4 @@ -# cython: language_level=3, boundscheck=True, wraparound=False, cdivision=False +# cython: language_level=3, boundscheck=False, wraparound=False, cdivision=False """ This is a solver linear equation systems with a penta-diagonal matrix, @@ -17,18 +17,18 @@ from libc.stdint cimport int64_t # === Main Python Interface === -def penta_solver1(double[:, :] mat_flat, double[:, :] rhs): +def penta_solver1(double[::, ::] mat_flat, double[::, ::] rhs): return np.asarray(c_penta_solver1(mat_flat, rhs)) -def penta_solver2(double[:, :] mat_flat, double[:] rhs): +def penta_solver2(double[::, ::] mat_flat, double[::, ::] rhs): return np.asarray(c_penta_solver2(mat_flat, rhs)) # === Solver Algorithm 1 === -cdef double[:, :] c_penta_solver1(double[:, :] mat_flat, double[:, :] rhs): +cdef double[::, ::] c_penta_solver1(double[::, ::] mat_flat, double[::, ::] rhs): """ Solves the pentadiagonal system of equations ``Ax = b`` with the matrix ``A`` and the right-hand side ``b`` by @@ -70,7 +70,7 @@ cdef double[:, :] c_penta_solver1(double[:, :] mat_flat, double[:, :] rhs): cdef void c_penta_factorize_algo1( - double[:, :] mat_flat, + double[::, ::] mat_flat, int64_t mat_n_rows, double[::, ::1] mat_factorized, ): @@ -260,53 +260,231 @@ cdef void c_solve_penta_from_factorize_algo_1( # === Solver Algorithm 2 === -cdef double[:] c_penta_solver2(double[:, :] mat_flat, double[:] rhs): +cdef double[::, ::] c_penta_solver2(double[::, ::] mat_flat, double[::, ::] rhs): + """ + Solves the pentadiagonal system of equations ``Ax = b`` with the matrix ``A`` and + the right-hand side ``b`` by + + - factorizing the matrix ``A`` into auxiliary coefficients and a unit lower + triangular matrix ``L`` + - transforming the right-hand side into a vector ``omega`` + - solving the system of equations ``Lx = omega`` by backward substitution + + """ + + # Variable declarations + + cdef int64_t mat_n_rows = mat_flat.shape[1] + cdef int64_t rhs_n_cols = rhs.shape[1] + cdef int64_t iter_col + cdef double[::, ::1] result = np.empty(shape=(mat_n_rows, rhs_n_cols)) + cdef double[::, ::1] mat_factorized = np.empty(shape=(mat_n_rows, 5)) + + # first, the matrix is factorized + c_penta_factorize_algo2( + mat_flat, + mat_n_rows, + mat_factorized, + ) - cdef int mat_j = mat_flat.shape[1] + # then, all the right-hand sides are solved + for iter_col in range(rhs_n_cols): + c_solve_penta_from_factorize_algo_2( + mat_n_rows, + mat_factorized, + rhs[::, iter_col], + result[::, iter_col], + ) - cdef double[:] result = np.zeros(mat_j) + return result - cdef double[:] ps = np.zeros(mat_j) # psi - cdef double[:] si = np.zeros(mat_j) # sigma - cdef double[:] ph = np.zeros(mat_j) # phi - cdef double[:] ro = np.zeros(mat_j) # rho - cdef double[:] we = np.zeros(mat_j) # w +cdef void c_penta_factorize_algo2( + double[::, ::] mat_flat, + int64_t mat_n_rows, + double[::, ::1] mat_factorized, +): + """ + Factorizes the pentadiagonal matrix ``A`` into - cdef int i + - auxiliary coefficients ``psi`` (``ps``), ``rho`` and ``b`` for the transformation + of the right-hand side + - a unit lower triangular matrix with the main diagonals ``phi`` and ``sigma`` + (``si``) for the following forward substitution. Its unit main diagonal is + implicit. - ps[mat_j-1] = mat_flat[2, mat_j-1] - si[mat_j-1] = mat_flat[3, mat_j-1] / ps[mat_j-1] - ph[mat_j-1] = mat_flat[4, mat_j-1] / ps[mat_j-1] - we[mat_j-1] = rhs[mat_j-1] / ps[mat_j-1] + They are overwriting the memoryview ``mat_factorized`` as follows: - ro[mat_j-2] = mat_flat[1, mat_j-2] - ps[mat_j-2] = mat_flat[2, mat_j-2] - si[mat_j-1] * ro[mat_j-2] - si[mat_j-2] = (mat_flat[3, mat_j-2] - ph[mat_j-1] * ro[mat_j-2]) / ps[mat_j-2] - ph[mat_j-2] = mat_flat[4, mat_j-2] / ps[mat_j-2] - we[mat_j-2] = (rhs[mat_j-2] - we[mat_j-1] * ro[mat_j-2]) / ps[mat_j-2] + ```bash + [[ * * ps_0 rho_0 b_i ] + [ * si_1 ps_1 rho_1 b_1 ] + [ phi_2 si_2 ps_2 rho_2 b_2 ] + ... + [ phi_i si_i ps_i rho_i b_i ] + ... + [ phi_{n-3} si_{n-3} ps_{n-3} rho_{n-3} b_{n-3} ] + [ phi_{n-2} si_{n-2} ps_{n-2} rho_{n-2} * ] + [ phi_{n-1} si_{n-1} ps_{n-1} * * ]] + ``` - for i in range(mat_j-3, 1, -1): - ro[i] = mat_flat[1, i] - si[i+2] * mat_flat[0, i] - ps[i] = mat_flat[2, i] - ph[i+2] * mat_flat[0, i] - si[i+1] * ro[i] - si[i] = (mat_flat[3, i] - ph[i+1] * ro[i]) / ps[i] - ph[i] = mat_flat[4, i] / ps[i] - we[i] = (rhs[i] - we[i+2] * mat_flat[0, i] - we[i+1] * ro[i]) / ps[i] + where the entries marked with ``*`` are not used by design, but overwritten with + zeros. - ro[1] = mat_flat[1, 1] - si[3] * mat_flat[0, 1] - ps[1] = mat_flat[2, 1] - ph[3] * mat_flat[0, 1] - si[2] * ro[1] - si[1] = (mat_flat[3, 1] - ph[2] * ro[1]) / ps[1] + """ - ro[0] = mat_flat[1, 0] - si[2] * mat_flat[0, 0] - ps[0] = mat_flat[2, 0] - ph[2] * mat_flat[0, 0] - si[1] * ro[0] + # === Variable declarations === - we[1] = (rhs[1] - we[3] * mat_flat[0, 1] - we[2] * ro[1]) / ps[1] - we[0] = (rhs[0] - we[2] * mat_flat[0, 0] - we[1] * ro[0]) / ps[0] + cdef int64_t iter_row + cdef double ps_i, rho_i # psi, rho + cdef double si_i, si_i_minus_1, si_i_plus_1 # sigma + cdef double phi_i, phi_i_minus_1, phi_i_plus_1 # phi - # Foreward substitution - result[0] = we[0] - result[1] = we[1] - si[1] * result[0] + # === Factorization === - for i in range(2, mat_j): - result[i] = we[i] - si[i] * result[i-1] - ph[i] * result[i-2] + # First row + ps_i = mat_flat[2, mat_n_rows - 1] + si_i_plus_1 = mat_flat[3, mat_n_rows - 1] / ps_i + phi_i_plus_1 = mat_flat[4, mat_n_rows - 1] / ps_i - return result + mat_factorized[mat_n_rows - 1, 0] = phi_i_plus_1 + mat_factorized[mat_n_rows - 1, 1] = si_i_plus_1 + mat_factorized[mat_n_rows - 1, 2] = ps_i + mat_factorized[mat_n_rows - 1, 3] = 0.0 + mat_factorized[mat_n_rows - 1, 4] = 0.0 + + # Second row + rho_i = mat_flat[1, mat_n_rows-2] + ps_i = mat_flat[2, mat_n_rows-2] - si_i_plus_1 * rho_i + si_i = (mat_flat[3, mat_n_rows-2] - phi_i_plus_1 * rho_i) / ps_i + phi_i = mat_flat[4, mat_n_rows-2] / ps_i + + mat_factorized[mat_n_rows - 2, 0] = phi_i + mat_factorized[mat_n_rows - 2, 1] = si_i + mat_factorized[mat_n_rows - 2, 2] = ps_i + mat_factorized[mat_n_rows - 2, 3] = rho_i + mat_factorized[mat_n_rows - 2, 4] = 0.0 + + # Central rows + for iter_row in range(mat_n_rows-3, 1, -1): + b_i = mat_flat[0, iter_row] + rho_i = mat_flat[1, iter_row] - si_i_plus_1 * b_i + ps_i = mat_flat[2, iter_row] - phi_i_plus_1 * b_i - si_i * rho_i + si_i_minus_1 = (mat_flat[3, iter_row] - phi_i * rho_i) / ps_i + si_i_plus_1 = si_i + si_i = si_i_minus_1 + phi_i_minus_1 = mat_flat[4, iter_row] / ps_i + phi_i_plus_1 = phi_i + phi_i = phi_i_minus_1 + + mat_factorized[iter_row, 0] = phi_i + mat_factorized[iter_row, 1] = si_i + mat_factorized[iter_row, 2] = ps_i + mat_factorized[iter_row, 3] = rho_i + mat_factorized[iter_row, 4] = b_i + + # Second to last row + b_i = mat_flat[0, 1] + rho_i = mat_flat[1, 1] - si_i_plus_1 * b_i + ps_i = mat_flat[2, 1] - phi_i_plus_1 * b_i - si_i * rho_i + si_i_minus_1 = (mat_flat[3, 1] - phi_i * rho_i) / ps_i + si_i_plus_1 = si_i + si_i = si_i_minus_1 + + mat_factorized[1, 0] = 0.0 + mat_factorized[1, 1] = si_i + mat_factorized[1, 2] = ps_i + mat_factorized[1, 3] = rho_i + mat_factorized[1, 4] = b_i + + # Last row + b_i = mat_flat[0, 0] + rho_i = mat_flat[1, 0] - si_i_plus_1 * b_i + ps_i = mat_flat[2, 0] - phi_i * b_i - si_i * rho_i + + mat_factorized[0, 0] = 0.0 + mat_factorized[0, 1] = 0.0 + mat_factorized[0, 2] = ps_i + mat_factorized[0, 3] = rho_i + mat_factorized[0, 4] = b_i + + return + + +cdef void c_solve_penta_from_factorize_algo_2( + int64_t mat_n_rows, + double[::, ::1] mat_factorized, + double[::] rhs_single, + double[::] result_view, +): + """ + Solves the pentadiagonal system of equations ``Ax = b`` with the factorized + unit lower triangular matrix ``L`` and the right-hand side ``b``. + It overwrites the right-hand side ``b`` first with the transformed vector ``omega`` + and then with the solution vector ``x`` for ``Lx = omega``. + + """ + + # === Variable declarations === + + cdef int64_t iter_row + cdef double om_i, om_i_minus_1, om_i_minus_2 # omega + + # === Transformation === + + # first, the right-hand side is transformed into the vector ``omega`` + # First row + om_i_plus_1 = rhs_single[mat_n_rows-1] / mat_factorized[mat_n_rows - 1, 2] + result_view[mat_n_rows-1] = om_i_plus_1 + + # Second row + om_i = ( + rhs_single[mat_n_rows-2] + - om_i_plus_1 * mat_factorized[mat_n_rows - 2, 3] + ) / mat_factorized[mat_n_rows - 2, 2] + result_view[mat_n_rows-2] = om_i + + # Central rows + for iter_row in range(mat_n_rows-3, 1, -1): + om_i_minus_1 = ( + rhs_single[iter_row] + - om_i_plus_1 * mat_factorized[iter_row, 4] + - om_i * mat_factorized[iter_row, 3] + ) / mat_factorized[iter_row, 2] + om_i_plus_1 = om_i + om_i = om_i_minus_1 + result_view[iter_row] = om_i + + # Second to last row + om_i_minus_1 = ( + rhs_single[1] + - om_i_plus_1 * mat_factorized[1, 4] + - om_i * mat_factorized[1, 3] + ) / mat_factorized[1, 2] + om_i_plus_1 = om_i + om_i = om_i_minus_1 + result_view[1] = om_i + + # Last row + om_i_minus_1 = ( + rhs_single[0] + - om_i_plus_1 * mat_factorized[0, 4] + - om_i * mat_factorized[0, 3] + ) / mat_factorized[0, 2] + result_view[0] = om_i_minus_1 + + # === Forward substitution === + + # The solution vector is calculated by forward substitution that overwrites the + # right-hand side vector with the solution vector + om_i -= mat_factorized[1, 1] * om_i_minus_1 + result_view[1] = om_i + + for iter_row in range(2, mat_n_rows): + result_view[iter_row] = ( + result_view[iter_row] + - mat_factorized[iter_row, 0] * om_i_minus_1 + - mat_factorized[iter_row, 1] * om_i + ) + om_i_minus_1 = om_i + om_i = result_view[iter_row] + + return \ No newline at end of file From ba5945f3e1b8ba10862d6c4e7302aa3e0f4dec00 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 20:16:15 +0200 Subject: [PATCH 23/34] tests: [11] unified test for both pentapy algorithms, made tests a lot more extensive --- ...t_solver_1.py => test_solvers_internal.py} | 22 +++++++++++++------ 1 file changed, 15 insertions(+), 7 deletions(-) rename tests/{test_solver_1.py => test_solvers_internal.py} (82%) diff --git a/tests/test_solver_1.py b/tests/test_solvers_internal.py similarity index 82% rename from tests/test_solver_1.py rename to tests/test_solvers_internal.py index 78956f6..2b543cf 100644 --- a/tests/test_solver_1.py +++ b/tests/test_solvers_internal.py @@ -1,5 +1,6 @@ """ -Test suite for testing the pentadiagonal solver based on Algorithm PTRANS-I. +Test suite for testing the pentadiagonal solver based on either Algorithm PTRANS-I or +PTRANS-II. """ @@ -40,13 +41,17 @@ 10_000, 10_001, ] -REF_WARNING = "pentapy: PTRANS-I not suitable for input-matrix." +REF_WARNING_CONTENT = "not suitable for input-matrix." +SOLVER_ALIASES_PTRANS_I = [1, "1", "PTRANS-I", "ptrans-i"] +SOLVER_ALIASES_PTRANS_II = [2, "2", "PTRANS-II", "ptrans-ii"] # === Tests === @pytest.mark.parametrize("induce_error", [False, True]) -@pytest.mark.parametrize("solver_alias", [1]) # "1", "PTRANS-I"]) +@pytest.mark.parametrize( + "solver_alias", SOLVER_ALIASES_PTRANS_I + SOLVER_ALIASES_PTRANS_II +) @pytest.mark.parametrize("input_layout", ["full", "banded_row_wise", "banded_col_wise"]) @pytest.mark.parametrize("n_rhs", [None, 1, 10]) @pytest.mark.parametrize("n_rows", N_ROWS) @@ -72,17 +77,20 @@ def test_penta_solver1( ill_conditioned=False, ) - # an error is induced by setting the first diagonal element to zero + # an error is induced by setting the first or last diagonal element to zero if induce_error: # the induction of the error is only possible if the matrix does not have # only 3 rows if n_rows == 3: pytest.skip( "Only 3 rows, cannot induce error because this will not go into " - "PTRANS-I, but NumPy" + "PTRANS-I, but NumPy." ) - mat_full[0, 0] = 0.0 + if solver_alias in SOLVER_ALIASES_PTRANS_I: + mat_full[0, 0] = 0.0 + else: + mat_full[n_rows - 1, n_rows - 1] = 0.0 # the right-hand side is generated np.random.seed(SEED) @@ -119,7 +127,7 @@ def test_penta_solver1( # Case 1: in case of an error, a warning has to be issued and the result has to # be NaN if induce_error: - with pytest.warns(UserWarning, match=REF_WARNING): + with pytest.warns(UserWarning, match=REF_WARNING_CONTENT): sol = pp.solve( mat=mat, rhs=rhs, From 026e8eab05dfe8efffe3d0eff0348db5e2728c13 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 20:24:26 +0200 Subject: [PATCH 24/34] tests: [11] removed superficial unittest --- tests/test_pentapy.py | 115 ------------------------------------------ 1 file changed, 115 deletions(-) delete mode 100755 tests/test_pentapy.py diff --git a/tests/test_pentapy.py b/tests/test_pentapy.py deleted file mode 100755 index 558de5e..0000000 --- a/tests/test_pentapy.py +++ /dev/null @@ -1,115 +0,0 @@ -""" -This is the unittest for pentapy. -""" - -import unittest - -# import platform -import warnings - -import numpy as np - -import pentapy as pp - -warnings.simplefilter("always") - - -class TestPentapy(unittest.TestCase): - def setUp(self): - self.seed = 19031977 - self.size = 1000 - self.rand = np.random.RandomState(self.seed) - self.mat = (self.rand.rand(5, self.size) - 0.5) * 1e-5 - self.rhs = self.rand.rand(self.size) * 1e5 - - def test_tools(self): - self.mat_int = np.zeros((100, 100), dtype=int) - # fill bands of pentadiagonal matrix - self.mat_int[pp.diag_indices(100, 0)] = self.rand.randint(1, 1000, size=100) - self.mat_int[pp.diag_indices(100, 1)] = self.rand.randint(1, 1000, size=99) - self.mat_int[pp.diag_indices(100, 2)] = self.rand.randint(1, 1000, size=98) - self.mat_int[pp.diag_indices(100, -1)] = self.rand.randint(1, 1000, size=99) - self.mat_int[pp.diag_indices(100, -2)] = self.rand.randint(1, 1000, size=98) - # create banded - self.mat_int_col = pp.create_banded(self.mat_int) - self.mat_int_row = pp.create_banded(self.mat_int, col_wise=False) - # create full - self.mat_int_col_ful = pp.create_full(self.mat_int_col, col_wise=True) - self.mat_int_row_ful = pp.create_full(self.mat_int_row, col_wise=False) - # shifting - self.mat_shift_cr = pp.shift_banded(self.mat_int_col) - self.mat_shift_rc = pp.shift_banded(self.mat_int_row, col_to_row=False) - # in place shifting - self.mat_int_col_ip = pp.create_banded(self.mat_int) - self.mat_int_row_ip = pp.create_banded(self.mat_int, col_wise=False) - pp.shift_banded(self.mat_int_col_ip, copy=False) - pp.shift_banded(self.mat_int_row_ip, copy=False, col_to_row=False) - # checking - self.assertEqual(np.sum(self.mat_int > 0), 494) - self.assertTrue(np.array_equal(self.mat_int_col, self.mat_shift_rc)) - self.assertTrue(np.array_equal(self.mat_int_row, self.mat_shift_cr)) - self.assertTrue(np.array_equal(self.mat_int_col, self.mat_int_row_ip)) - self.assertTrue(np.array_equal(self.mat_int_row, self.mat_int_col_ip)) - self.assertTrue(np.array_equal(self.mat_int, self.mat_int_col_ful)) - self.assertTrue(np.array_equal(self.mat_int, self.mat_int_row_ful)) - - def test_solve1(self): - self.mat_col = pp.shift_banded(self.mat, col_to_row=False) - self.mat_ful = pp.create_full(self.mat, col_wise=False) - - sol_row = pp.solve(self.mat, self.rhs, is_flat=True, solver=1) - sol_col = pp.solve( - self.mat_col, - self.rhs, - is_flat=True, - index_row_wise=False, - solver=1, - ) - sol_ful = pp.solve(self.mat_ful, self.rhs, solver=1) - - diff_row = np.max(np.abs(np.dot(self.mat_ful, sol_row) - self.rhs)) - diff_col = np.max(np.abs(np.dot(self.mat_ful, sol_col) - self.rhs)) - diff_ful = np.max(np.abs(np.dot(self.mat_ful, sol_ful) - self.rhs)) - - diff_row_col = np.max(np.abs(self.mat_ful - pp.create_full(self.mat_col))) - self.assertAlmostEqual(diff_row * 1e-5, 0.0) - self.assertAlmostEqual(diff_col * 1e-5, 0.0) - self.assertAlmostEqual(diff_ful * 1e-5, 0.0) - self.assertAlmostEqual(diff_row_col * 1e5, 0.0) - - def test_solve2(self): - self.mat_col = pp.shift_banded(self.mat, col_to_row=False) - self.mat_ful = pp.create_full(self.mat, col_wise=False) - - sol_row = pp.solve(self.mat, self.rhs, is_flat=True, solver=2) - sol_col = pp.solve( - self.mat_col, - self.rhs, - is_flat=True, - index_row_wise=False, - solver=2, - ) - sol_ful = pp.solve(self.mat_ful, self.rhs, solver=2) - - diff_row = np.max(np.abs(np.dot(self.mat_ful, sol_row) - self.rhs)) - diff_col = np.max(np.abs(np.dot(self.mat_ful, sol_col) - self.rhs)) - diff_ful = np.max(np.abs(np.dot(self.mat_ful, sol_ful) - self.rhs)) - - diff_row_col = np.max(np.abs(self.mat_ful - pp.create_full(self.mat_col))) - self.assertAlmostEqual(diff_row * 1e-5, 0.0) - self.assertAlmostEqual(diff_col * 1e-5, 0.0) - self.assertAlmostEqual(diff_ful * 1e-5, 0.0) - self.assertAlmostEqual(diff_row_col * 1e5, 0.0) - - def test_error(self): - self.err_mat = np.array( - [[3, 2, 1, 0], [-3, -2, 7, 1], [3, 2, -1, 5], [0, 1, 2, 3]] - ) - self.err_rhs = np.array([6, 3, 9, 6]) - sol_2 = pp.solve(self.err_mat, self.err_rhs, is_flat=False, solver=2) - diff_2 = np.max(np.abs(np.dot(self.err_mat, sol_2) - self.err_rhs)) - self.assertAlmostEqual(diff_2, 0.0) - - -if __name__ == "__main__": - unittest.main() From f6fa4cbb28cb18b2d921bc4f6ad3243d539cdad0 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 20:25:26 +0200 Subject: [PATCH 25/34] misc: [11] fixed typos --- examples/README.rst | 4 ++-- paper/paper.md | 6 +++--- src/pentapy/core.py | 12 ++++++------ src/pentapy/tools.py | 22 +++++++++++----------- 4 files changed, 22 insertions(+), 22 deletions(-) diff --git a/examples/README.rst b/examples/README.rst index ea7bac4..a1b7ae0 100644 --- a/examples/README.rst +++ b/examples/README.rst @@ -88,7 +88,7 @@ If M is a full matrix, you call the following: X = pp.solve(M, Y) -If M is flattend in row-wise order you have to set the keyword argument ``is_flat=True``: +If M is flattened in row-wise order you have to set the keyword argument ``is_flat=True``: .. code-block:: python @@ -99,7 +99,7 @@ If M is flattend in row-wise order you have to set the keyword argument ``is_fla X = pp.solve(M, Y, is_flat=True) -If you got a col-wise flattend matrix you have to set ``index_row_wise=False``: +If you got a col-wise flattened matrix you have to set ``index_row_wise=False``: .. code-block:: python diff --git a/paper/paper.md b/paper/paper.md index 53c6f04..fa0fdc3 100644 --- a/paper/paper.md +++ b/paper/paper.md @@ -51,8 +51,8 @@ $$ Here, $d_i$ are the diagonal entries and $d_i^{(j)}$ represent the $j$-th minor diagonal. -Recently, @askar presented two algorithms to -solve the linear systems of equations for $X$, ``PTRANS-I`` and ``PTRANS-II``, +Recently, @askar presented two algorithms to +solve the linear systems of equations for $X$, ``PTRANS-I`` and ``PTRANS-II``, applying first transformation to a triangular matrix and then, respectively, backward and forward substitution. ``pentapy`` provides Cython [@cython] implementations of these algorithms and a set of tools to convert matrices to row-wise or @@ -73,7 +73,7 @@ The linear algebra solver of NumPy [@numpy] served as a standard reference, whic ``pentapy`` is designed to provide a fast solver for the special case of a pentadiagonal linear system. To the best of the author's knowledge, this package outperforms the current algorithms for solving pentadiagonal systems in Python. -The solver can handle different input formats of the coefficient matrix, i.e., a flattend matrix or a +The solver can handle different input formats of the coefficient matrix, i.e., a flattened matrix or a quadratic matrix. diff --git a/src/pentapy/core.py b/src/pentapy/core.py index ccf0c8f..ad704a6 100644 --- a/src/pentapy/core.py +++ b/src/pentapy/core.py @@ -46,8 +46,8 @@ def solve( """ Solver for a pentadiagonal system. - The matrix can be given as a full n x n matrix or as a flattend one. - The flattend matrix can be given in a row-wise flattend form:: + The matrix can be given as a full n x n matrix or as a flattened one. + The flattened matrix can be given in a row-wise flattened form:: [[Dup2[0] Dup2[1] Dup2[2] ... Dup2[N-2] 0 0 ] [Dup1[0] Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] 0 ] @@ -55,7 +55,7 @@ def solve( [0 Dlow1[1] Dlow1[2] ... Dlow1[N-2] Dlow1[N-1] Dlow1[N]] [0 0 Dlow2[2] ... Dlow2[N-2] Dlow2[N-2] Dlow2[N]]] - Or a column-wise flattend form:: + Or a column-wise flattened form:: [[0 0 Dup2[2] ... Dup2[N-2] Dup2[N-1] Dup2[N] ] [0 Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] Dup1[N] ] @@ -65,7 +65,7 @@ def solve( Dup1 and Dup2 are the first and second upper minor-diagonals and Dlow1 resp. Dlow2 are the lower ones. - If you provide a column-wise flattend matrix, you have to set:: + If you provide a column-wise flattened matrix, you have to set:: index_row_wise=False @@ -77,9 +77,9 @@ def solve( rhs : :class:`numpy.ndarray` of shape (m,) or (m, n) The right hand side(s) of the equation system. Its shape is preserved. is_flat : :class:`bool`, default=False - State if the matrix is already flattend. Default: ``False`` + State if the matrix is already flattened. Default: ``False`` index_row_wise : :class:`bool`, default=True - State if the flattend matrix is row-wise flattend. Default: ``True`` + State if the flattened matrix is row-wise flattened. Default: ``True`` solver : :class:`int` or :class:`str`, default=1 Which solver should be used. The following are provided: diff --git a/src/pentapy/tools.py b/src/pentapy/tools.py index ca80d27..4bb0b52 100644 --- a/src/pentapy/tools.py +++ b/src/pentapy/tools.py @@ -52,8 +52,8 @@ def shift_banded(mat, up=2, low=2, col_to_row=True, copy=True): Either from column-wise to row-wise storage or vice versa. - The Matrix has to be given as a flattend matrix. - Either in a column-wise flattend form:: + The Matrix has to be given as a flattened matrix. + Either in a column-wise flattened form:: [[0 0 Dup2[2] ... Dup2[N-2] Dup2[N-1] Dup2[N] ] [0 Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] Dup1[N] ] @@ -65,7 +65,7 @@ def shift_banded(mat, up=2, low=2, col_to_row=True, copy=True): col_to_row=True - Or in a row-wise flattend form:: + Or in a row-wise flattened form:: [[Dup2[0] Dup2[1] Dup2[2] ... Dup2[N-2] 0 0 ] [Dup1[0] Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] 0 ] @@ -98,7 +98,7 @@ def shift_banded(mat, up=2, low=2, col_to_row=True, copy=True): Returns ------- :class:`numpy.ndarray` - Shifted bandend matrix + Shifted banded matrix """ if copy: mat_flat = np.copy(mat) @@ -124,8 +124,8 @@ def shift_banded(mat, up=2, low=2, col_to_row=True, copy=True): def create_banded(mat, up=2, low=2, col_wise=True, dtype=None): """Create a banded matrix from a given quadratic Matrix. - The Matrix will to be returned as a flattend matrix. - Either in a column-wise flattend form:: + The Matrix will to be returned as a flattened matrix. + Either in a column-wise flattened form:: [[0 0 Dup2[2] ... Dup2[N-2] Dup2[N-1] Dup2[N] ] [0 Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] Dup1[N] ] @@ -137,7 +137,7 @@ def create_banded(mat, up=2, low=2, col_wise=True, dtype=None): col_wise=True - Or in a row-wise flattend form:: + Or in a row-wise flattened form:: [[Dup2[0] Dup2[1] Dup2[2] ... Dup2[N-2] 0 0 ] [Dup1[0] Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] 0 ] @@ -168,7 +168,7 @@ def create_banded(mat, up=2, low=2, col_wise=True, dtype=None): Returns ------- :class:`numpy.ndarray` - Bandend matrix + Banded matrix """ mat = np.asanyarray(mat) if mat.ndim != 2: @@ -202,8 +202,8 @@ def create_banded(mat, up=2, low=2, col_wise=True, dtype=None): def create_full(mat, up=2, low=2, col_wise=True): """Create a (n x n) Matrix from a given banded matrix. - The given Matrix has to be a flattend matrix. - Either in a column-wise flattend form:: + The given Matrix has to be a flattened matrix. + Either in a column-wise flattened form:: [[0 0 Dup2[2] ... Dup2[N-2] Dup2[N-1] Dup2[N] ] [0 Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] Dup1[N] ] @@ -215,7 +215,7 @@ def create_full(mat, up=2, low=2, col_wise=True): col_wise=True - Or in a row-wise flattend form:: + Or in a row-wise flattened form:: [[Dup2[0] Dup2[1] Dup2[2] ... Dup2[N-2] 0 0 ] [Dup1[0] Dup1[1] Dup1[2] ... Dup1[N-2] Dup1[N-1] 0 ] From 9b3e1224b205f50caec5fa136771680d952bab55 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 20:39:43 +0200 Subject: [PATCH 26/34] style/refactor: [11] made core code readable to humans and not only a machine --- src/pentapy/core.py | 55 +++++++++++++++++++++++++++------------------ 1 file changed, 33 insertions(+), 22 deletions(-) diff --git a/src/pentapy/core.py b/src/pentapy/core.py index ad704a6..4158d94 100644 --- a/src/pentapy/core.py +++ b/src/pentapy/core.py @@ -101,6 +101,7 @@ def solve( # first, the solver is converted to the internal name to avoid confusion solver_inter = pmodels._SOLVER_ALIAS_CONVERSIONS[str(solver).lower()] + # Case 1: the pentapy solvers if solver_inter in { pmodels.PentaSolverAliases.PTRANS_I, pmodels.PentaSolverAliases.PTRANS_II, @@ -150,29 +151,14 @@ def solve( warnings.warn("pentapy: PTRANS-I not suitable for input-matrix.") return np.full(shape=rhs_og_shape, fill_value=np.nan) + # Case 2: LAPACK's banded solver elif solver_inter == pmodels.PentaSolverAliases.LAPACK: # pragma: no cover try: from scipy.linalg import solve_banded except ImportError as imp_err: # pragma: no cover msg = "pentapy.solve: scipy.linalg.solve_banded could not be imported" raise ValueError(msg) from imp_err - if is_flat and index_row_wise: - mat_flat = np.array(mat) - ptools._check_penta(mat_flat) - ptools.shift_banded(mat_flat, col_to_row=False, copy=False) - elif is_flat: - mat_flat = np.asarray(mat) - else: - mat_flat = ptools.create_banded(mat) - return solve_banded((2, 2), mat_flat, rhs) - elif solver_inter == pmodels.PentaSolverAliases.SUPER_LU: # pragma: no cover - try: - from scipy import sparse as sps - from scipy.sparse.linalg import spsolve - except ImportError as imp_err: - msg = "pentapy.solve: scipy.sparse could not be imported" - raise ValueError(msg) from imp_err if is_flat and index_row_wise: mat_flat = np.array(mat) ptools._check_penta(mat_flat) @@ -181,17 +167,27 @@ def solve( mat_flat = np.asarray(mat) else: mat_flat = ptools.create_banded(mat) - size = mat_flat.shape[1] - M = sps.spdiags(mat_flat, [2, 1, 0, -1, -2], size, size, format="csc") - return spsolve(M, rhs, use_umfpack=False) - elif solver_inter == pmodels.PentaSolverAliases.UMFPACK: # pragma: no cover + # NOTE: since this is a general banded solver, the number of sub- and super- + # diagonals has to be provided + return solve_banded( + l_and_u=(2, 2), + ab=mat_flat, + b=rhs, + ) + + # Case 3: SciPy's sparse solver with or without UMFPACK + elif solver_inter in { + pmodels.PentaSolverAliases.SUPER_LU, + pmodels.PentaSolverAliases.UMFPACK, + }: try: from scipy import sparse as sps from scipy.sparse.linalg import spsolve except ImportError as imp_err: msg = "pentapy.solve: scipy.sparse could not be imported" raise ValueError(msg) from imp_err + if is_flat and index_row_wise: mat_flat = np.array(mat) ptools._check_penta(mat_flat) @@ -200,9 +196,24 @@ def solve( mat_flat = np.asarray(mat) else: mat_flat = ptools.create_banded(mat) + + # the solvers require a sparse left-hand side matrix, so this is created here + # NOTE: the UMFPACK solver will not be triggered for multiple right-hand sides + use_umfpack = solver_inter == pmodels.PentaSolverAliases.UMFPACK size = mat_flat.shape[1] - M = sps.spdiags(mat_flat, [2, 1, 0, -1, -2], size, size, format="csc") - return spsolve(M, rhs, use_umfpack=True) + M = sps.spdiags( + data=mat_flat, + diags=[2, 1, 0, -1, -2], + m=size, + n=size, + format="csc", + ) + + return spsolve( + A=M, + b=rhs, + use_umfpack=use_umfpack, + ) else: # pragma: no cover msg = f"pentapy.solve: unknown solver ({solver})" From 92abbbdaca93083334917c7db42fe63e63ed6c23 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 21:27:50 +0200 Subject: [PATCH 27/34] feat/refactor: [11] ensured all solvers behave the same in terms of output and warning behaviour; tested it altogether --- src/pentapy/core.py | 32 ++++++--- tests/test_solvers_external.py | 121 +++++++++++++++++++++++++++++++++ tests/test_solvers_internal.py | 2 +- 3 files changed, 145 insertions(+), 10 deletions(-) create mode 100644 tests/test_solvers_external.py diff --git a/src/pentapy/core.py b/src/pentapy/core.py index 4158d94..c701573 100644 --- a/src/pentapy/core.py +++ b/src/pentapy/core.py @@ -142,10 +142,11 @@ def solve( ) # if there was only a 1D right-hand side, the result has to be flattened + sol = solver_func(mat_flat, rhs) if single_rhs: - return solver_func(mat_flat, rhs).ravel() + sol = sol.ravel() - return solver_func(mat_flat, rhs) + return sol except ZeroDivisionError: warnings.warn("pentapy: PTRANS-I not suitable for input-matrix.") @@ -170,11 +171,17 @@ def solve( # NOTE: since this is a general banded solver, the number of sub- and super- # diagonals has to be provided - return solve_banded( - l_and_u=(2, 2), - ab=mat_flat, - b=rhs, - ) + # NOTE: LAPACK handles all the reshaping and flattening internally + try: + return solve_banded( + l_and_u=(2, 2), + ab=mat_flat, + b=rhs, + ) + + except np.linalg.LinAlgError: + warnings.warn("pentapy: LAPACK solver encountered singular matrix.") + return np.full(shape=rhs.shape, fill_value=np.nan) # Case 3: SciPy's sparse solver with or without UMFPACK elif solver_inter in { @@ -184,7 +191,7 @@ def solve( try: from scipy import sparse as sps from scipy.sparse.linalg import spsolve - except ImportError as imp_err: + except ImportError as imp_err: # pragma: no cover msg = "pentapy.solve: scipy.sparse could not be imported" raise ValueError(msg) from imp_err @@ -209,12 +216,19 @@ def solve( format="csc", ) - return spsolve( + sol = spsolve( A=M, b=rhs, use_umfpack=use_umfpack, ) + # NOTE: spsolve flattens column-vectors, thus their shape has to be restored + # NOTE: it already fills the result vector with NaNs if the matrix is singular + if rhs.ndim == 2 and 1 in rhs.shape: + sol = sol[::, np.newaxis] + + return sol + else: # pragma: no cover msg = f"pentapy.solve: unknown solver ({solver})" raise ValueError(msg) diff --git a/tests/test_solvers_external.py b/tests/test_solvers_external.py new file mode 100644 index 0000000..075dfb4 --- /dev/null +++ b/tests/test_solvers_external.py @@ -0,0 +1,121 @@ +""" +Test suite for testing the external solvers that can be called via pentapy. The tests +are not exhaustive and only check whether the solvers can be called and return a +solution. + +""" + +# === Imports === + +from typing import Literal + +import numpy as np +import pentapy as pp +import pytest +import util_funcs as uf + +# === Constants === + +SEED = 19_031_977 +N_ROWS = [ + 3, + 4, + 5, + 10, + 11, + 25, + 26, + 50, + 51, +] +REF_WARNING_CONTENT = "singular" +SOLVER_ALIASES_LAPACK = [3, "3", "lapack", "LaPaCk"] +SOLVER_ALIASES_SPSOLVE = [4, "4", "spsolve", "SpSoLvE"] + +# === Tests === + + +@pytest.mark.parametrize("induce_error", [False, True]) +@pytest.mark.parametrize("solver_alias", SOLVER_ALIASES_LAPACK + SOLVER_ALIASES_SPSOLVE) +@pytest.mark.parametrize("input_layout", ["full", "banded_row_wise", "banded_col_wise"]) +@pytest.mark.parametrize("n_rhs", [None, 1, 10]) +@pytest.mark.parametrize("n_rows", N_ROWS) +def test_external_solvers( + n_rows: int, + n_rhs: int, + input_layout: Literal["full", "banded_row_wise", "banded_col_wise"], + solver_alias: Literal[1, "1", "PTRANS-I"], + induce_error: bool, +) -> None: + """ + Tests the external bindings for solving pentadiagonal systems starting from + different input layouts, number of right-hand sides, number of rows, and when an + error is induced by a zero matrix. + It has to be ensured that the edge case of ``n_rows = 3`` is also covered. + + """ + + # first, a random pentadiagonal matrix is generated + mat_full = np.zeros(shape=(n_rows, n_rows)) + if not induce_error: + mat_full[::, ::] = uf.gen_conditioned_rand_penta_matrix_dense( + n_rows=n_rows, + seed=SEED, + ill_conditioned=False, + ) + + # the right-hand side is generated + np.random.seed(SEED) + if n_rhs is not None: + rhs = np.random.rand(n_rows, n_rhs) + result_shape = (n_rows, n_rhs) + else: + rhs = np.random.rand(n_rows) + result_shape = (n_rows,) + + # the matrix is converted to the desired layout + if input_layout == "full": + mat = mat_full + kwargs = dict(is_flat=False) + + elif input_layout == "banded_row_wise": + mat = pp.create_banded(mat_full, col_wise=False) + kwargs = dict( + is_flat=True, + index_row_wise=True, + ) + + elif input_layout == "banded_col_wise": + mat = pp.create_banded(mat_full, col_wise=True) + kwargs = dict( + is_flat=True, + index_row_wise=False, + ) + + else: + raise ValueError(f"Invalid input layout: {input_layout}") + + # the solution is computed + # Case 1: in case of an error, a warning has to be issued and the result has to + # be NaN + if induce_error: + with pytest.warns(UserWarning, match=REF_WARNING_CONTENT): + sol = pp.solve( + mat=mat, + rhs=rhs, + solver=solver_alias, # type: ignore + **kwargs, + ) + assert sol.shape == result_shape + assert np.isnan(sol).all() + + return + + # Case 2: in case of no error, the solution can be computed without any issues + sol = pp.solve( + mat=mat, + rhs=rhs, + solver=solver_alias, # type: ignore + **kwargs, + ) + assert sol.shape == result_shape diff --git a/tests/test_solvers_internal.py b/tests/test_solvers_internal.py index 2b543cf..cdc55fa 100644 --- a/tests/test_solvers_internal.py +++ b/tests/test_solvers_internal.py @@ -55,7 +55,7 @@ @pytest.mark.parametrize("input_layout", ["full", "banded_row_wise", "banded_col_wise"]) @pytest.mark.parametrize("n_rhs", [None, 1, 10]) @pytest.mark.parametrize("n_rows", N_ROWS) -def test_penta_solver1( +def test_pentapy_solvers( n_rows: int, n_rhs: int, input_layout: Literal["full", "banded_row_wise", "banded_col_wise"], From fb0cf2a342240e3b637774f2f3a8661424d1d515 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 21:28:23 +0200 Subject: [PATCH 28/34] fix: [11] fixed coverage typo --- src/pentapy/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pentapy/__init__.py b/src/pentapy/__init__.py index 655c428..3705064 100644 --- a/src/pentapy/__init__.py +++ b/src/pentapy/__init__.py @@ -44,7 +44,7 @@ try: from pentapy._version import __version__ -except ImportError: # pragma: nocover +except ImportError: # pragma: no cover # package is not installed __version__ = "0.0.0.dev0" From c2e8032007a172528adb596f667eb44aedc89c3b Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 21:28:50 +0200 Subject: [PATCH 29/34] lint: [11] linted cython files --- src/pentapy/solver.pyx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/pentapy/solver.pyx b/src/pentapy/solver.pyx index 1494a25..2da01db 100644 --- a/src/pentapy/solver.pyx +++ b/src/pentapy/solver.pyx @@ -104,9 +104,9 @@ cdef void c_penta_factorize_algo1( # === Variable declarations === cdef int64_t iter_row - cdef double mu_i, ga_i, e_i # mu, gamma, e - cdef double al_i, al_i_minus_1, al_i_plus_1 # alpha - cdef double be_i, be_i_minus_1, be_i_plus_1 # beta + cdef double mu_i, ga_i, e_i # mu, gamma, e + cdef double al_i, al_i_minus_1, al_i_plus_1 # alpha + cdef double be_i, be_i_minus_1, be_i_plus_1 # beta # === Factorization === @@ -334,9 +334,9 @@ cdef void c_penta_factorize_algo2( # === Variable declarations === cdef int64_t iter_row - cdef double ps_i, rho_i # psi, rho - cdef double si_i, si_i_minus_1, si_i_plus_1 # sigma - cdef double phi_i, phi_i_minus_1, phi_i_plus_1 # phi + cdef double ps_i, rho_i # psi, rho + cdef double si_i, si_i_minus_1, si_i_plus_1 # sigma + cdef double phi_i, phi_i_minus_1, phi_i_plus_1 # phi # === Factorization === @@ -426,7 +426,7 @@ cdef void c_solve_penta_from_factorize_algo_2( # === Variable declarations === cdef int64_t iter_row - cdef double om_i, om_i_minus_1, om_i_minus_2 # omega + cdef double om_i, om_i_minus_1, om_i_plus_1 # omega # === Transformation === @@ -487,4 +487,4 @@ cdef void c_solve_penta_from_factorize_algo_2( om_i_minus_1 = om_i om_i = result_view[iter_row] - return \ No newline at end of file + return From b4190f8144a24a122cf79c6724bb5ef938278dae Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 21:31:46 +0200 Subject: [PATCH 30/34] package: [11] made requirements file-based and dynamic to make facilitate development without compromising on build --- pyproject.toml | 33 ++++++--------------------------- requirements/all.txt | 2 ++ requirements/base.txt | 1 + requirements/check.txt | 4 ++++ requirements/doc.txt | 8 ++++++++ requirements/scipy.txt | 1 + requirements/test.txt | 4 ++++ requirements/umfpack | 1 + 8 files changed, 27 insertions(+), 27 deletions(-) create mode 100644 requirements/all.txt create mode 100644 requirements/base.txt create mode 100644 requirements/check.txt create mode 100644 requirements/doc.txt create mode 100644 requirements/scipy.txt create mode 100644 requirements/test.txt create mode 100644 requirements/umfpack diff --git a/pyproject.toml b/pyproject.toml index 4c400f8..b26ae54 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,7 +14,7 @@ name = "pentapy" authors = [{name = "Sebastian Müller", email = "info@geostat-framework.org"}] readme = "README.md" license = {text = "MIT"} -dynamic = ["version"] +dynamic = ["version", "dependencies", "optional-dependencies"] description = "pentapy: A toolbox for pentadiagonal matrizes." classifiers = [ "Development Status :: 5 - Production/Stable", @@ -35,32 +35,10 @@ classifiers = [ "Topic :: Scientific/Engineering", "Topic :: Utilities", ] -dependencies = ["numpy>=1.20.0"] -[project.optional-dependencies] -scipy = ["scipy"] -umfpack = ["scikit-umfpack"] -all = [ - "scipy", - "scikit-umfpack", -] -doc = [ - "m2r2>=0.2.8", - "scipy>=1.1.0", - "matplotlib>=3", - "perfplot<0.9", - "numpydoc>=1.1", - "sphinx>=7", - "sphinx-gallery>=0.8", - "sphinx-rtd-theme>=2", -] -test = ["pytest-cov>=3"] -check = [ - "black>=24,<25", - "isort[colors]", - "pylint", - "cython-lint", -] +[tool.setuptools.dynamic] +dependencies = {file = ["requirements/base.txt"]} +optional-dependencies = {scipy = {file = ["requirements/scipy.txt"]}, umfpack = {file = ["requirements/umfpack.txt"]}, all = {file = ["requirements/all.txt"]}, doc = {file = ["requirements/doc.txt"]}, test = {file = ["requirements/test.txt"]}, check = {file = ["requirements/check.txt"]}} [project.urls] Homepage = "https://github.com/GeoStat-Framework/pentapy" @@ -103,7 +81,8 @@ max-line-length = 120 "*examples*", "*tests*", "*paper*", - "pentapy/src/pentapy/_version.py", + "src/pentapy/_version.py", + "src/pentapy/__init__.py", ] [tool.coverage.report] diff --git a/requirements/all.txt b/requirements/all.txt new file mode 100644 index 0000000..be8d325 --- /dev/null +++ b/requirements/all.txt @@ -0,0 +1,2 @@ +scikit-umfpack +scipy \ No newline at end of file diff --git a/requirements/base.txt b/requirements/base.txt new file mode 100644 index 0000000..19b3787 --- /dev/null +++ b/requirements/base.txt @@ -0,0 +1 @@ +numpy>=1.20.0 \ No newline at end of file diff --git a/requirements/check.txt b/requirements/check.txt new file mode 100644 index 0000000..4af46fc --- /dev/null +++ b/requirements/check.txt @@ -0,0 +1,4 @@ +black>=24,<25 +isort[colors] +pylint +cython-lint \ No newline at end of file diff --git a/requirements/doc.txt b/requirements/doc.txt new file mode 100644 index 0000000..c49be85 --- /dev/null +++ b/requirements/doc.txt @@ -0,0 +1,8 @@ +m2r2>=0.2.8 +scipy>=1.1.0 +matplotlib>=3 +perfplot<0.9 +numpydoc>=1.1 +sphinx>=7 +sphinx-gallery>=0.8 +sphinx-rtd-theme>=2 \ No newline at end of file diff --git a/requirements/scipy.txt b/requirements/scipy.txt new file mode 100644 index 0000000..9c61c73 --- /dev/null +++ b/requirements/scipy.txt @@ -0,0 +1 @@ +scipy \ No newline at end of file diff --git a/requirements/test.txt b/requirements/test.txt new file mode 100644 index 0000000..2f8c0c7 --- /dev/null +++ b/requirements/test.txt @@ -0,0 +1,4 @@ +pytest>=8 +pytest-cov>=3 +pytest-xdist>=3 +scipy>=1.1.0 \ No newline at end of file diff --git a/requirements/umfpack b/requirements/umfpack new file mode 100644 index 0000000..a8630c1 --- /dev/null +++ b/requirements/umfpack @@ -0,0 +1 @@ +scikit-umfpack \ No newline at end of file From 4c1e3a0827ec42a78a31be65d4aba32077c3d7b8 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 21:32:11 +0200 Subject: [PATCH 31/34] feat: [11] updated chagelog --- CHANGELOG.md | 34 ++++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 1417c4c..b197caf 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,39 @@ All notable changes to **pentapy** will be documented in this file. +## [1.4.0] - 2024-06 + +See [#22](https://github.com/GeoStat-Framework/pentapy/pull/22) + +### Enhancements + +- added support for multiple right-hand sides (currently serial) +- improved error handling and added debug information to error messages + +### Changes + +- shotgun refactored and documented the Cython implementation of PTRANS-I and PTRANS-II for single and multiple right-hand sides support +- fully typed the function ``pentapy.solve`` +- made internal solver alias handling of ``pentapy.solve`` smarter, more robust, and removed all duplicate code +- gave all solvers a consistent interface +- made code in ``pentapy.core`` more human-readable and maintainable and added comments +- fixed typos in documentation + +### Bugfixes + +- fixed error handling in case of zero-division to trigger dead error handling branch (see [Issue 23](https://github.com/GeoStat-Framework/pentapy/issues/23)) +- fixed edge case error for row/column of 3 (see [Issue 24](https://github.com/GeoStat-Framework/pentapy/issues/24)) + +### Tests + +- transitioned from ``unittest``-based testing to fully ``pytest``-based testing with parametrized and parallelized exhaustive testing (see [Issue 25](https://github.com/GeoStat-Framework/pentapy/issues/25)) +- made actual tests more meaningful by comparing them to LAPACK as reference standard (see [Issue 25](https://github.com/GeoStat-Framework/pentapy/issues/25)) +- included external solver bindings accessible via ``pentapy.solve`` as part of the test suite +- increased true coverage (not line-hit coverage) close to 100% + +### Packaging + +- made dependency specification file-based and dynamic ## [1.3.0] - 2024-04 @@ -100,6 +133,7 @@ This is the first release of pentapy, a python toolbox for solving pentadiagonal The solver is implemented in cython, which makes it really fast. +[1.4.0]: https://github.com/GeoStat-Framework/pentapy/compare/v1.3.0...v1.4.0 [1.3.0]: https://github.com/GeoStat-Framework/pentapy/compare/v1.2.0...v1.3.0 [1.2.0]: https://github.com/GeoStat-Framework/pentapy/compare/v1.1.2...v1.2.0 [1.1.2]: https://github.com/GeoStat-Framework/pentapy/compare/v1.1.1...v1.1.2 From b4f79e5fd3b73e666c738d3756e223ef952930a0 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 21:53:28 +0200 Subject: [PATCH 32/34] fix: [11] fixed wrong coverage exclude --- src/pentapy/core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pentapy/core.py b/src/pentapy/core.py index c701573..456155e 100644 --- a/src/pentapy/core.py +++ b/src/pentapy/core.py @@ -153,7 +153,7 @@ def solve( return np.full(shape=rhs_og_shape, fill_value=np.nan) # Case 2: LAPACK's banded solver - elif solver_inter == pmodels.PentaSolverAliases.LAPACK: # pragma: no cover + elif solver_inter == pmodels.PentaSolverAliases.LAPACK: try: from scipy.linalg import solve_banded except ImportError as imp_err: # pragma: no cover From a2edcb4d7d458e6ffd2cc095d42312b4095032fc Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 22:34:59 +0200 Subject: [PATCH 33/34] doc: [11] improved wording for preserving shape --- src/pentapy/core.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/pentapy/core.py b/src/pentapy/core.py index 456155e..56e2a54 100644 --- a/src/pentapy/core.py +++ b/src/pentapy/core.py @@ -75,7 +75,8 @@ def solve( mat : :class:`numpy.ndarray` of shape (m, m) or (5, m) The full or flattened version of the pentadiagonal matrix. rhs : :class:`numpy.ndarray` of shape (m,) or (m, n) - The right hand side(s) of the equation system. Its shape is preserved. + The right hand side(s) of the equation system. Its shape determines the shape + of the output as they will be identical. is_flat : :class:`bool`, default=False State if the matrix is already flattened. Default: ``False`` index_row_wise : :class:`bool`, default=True From 007bc5a0f794bcabe6edab86b2f10038ce84a4d9 Mon Sep 17 00:00:00 2001 From: MothNik Date: Sat, 8 Jun 2024 22:36:04 +0200 Subject: [PATCH 34/34] fix: [11] fixed changelog --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b197caf..d13b435 100755 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,7 @@ All notable changes to **pentapy** will be documented in this file. ## [1.4.0] - 2024-06 -See [#22](https://github.com/GeoStat-Framework/pentapy/pull/22) +See [#26](https://github.com/GeoStat-Framework/pentapy/pull/26) ### Enhancements