I want to attract attention to my code

[quant12345]

Feature description

Hi, everyone!
A few weeks ago I made a pull request. Replacing loops and conditional statements in the jacobi iteration method function with vector operations numpy.

Which on average gives an acceleration of the calculation of the entire algorithm by about 50 times!
Attached scripts to test performance and graphics. But, for some reason, there is no response to my code.

I am attaching the code here: the function names have been changed to 'jacobi_old' and 'jacobi_new', the 'strictly_diagonally_dominant' function has remained unchanged. The code starts a loop, where n is the size of the generated arrays: coefficient, constant, init_val. For the convergence of the method, the diagonal elements of the coefficient array are multiplied by 10 so that there is a diagonal dominance of the matrix. At each iteration, the size of the array, the time spent by the old and the new algorithm(in seconds), the ratio of time are printed, and the resulting arrays are compared (are they the same).
If you have any questions, I am open to dialogue.

Code that runs the new and old code for comparison

"""
Jacobi Iteration Method - https://en.wikipedia.org/wiki/Jacobi_method
"""
from __future__ import annotations

import datetime

import numpy as np
from numpy import float64
from numpy.typing import NDArray


# Method to find solution of system of linear equations

def jacobi_old(
        coefficient_matrix: NDArray[float64],
        constant_matrix: NDArray[float64],
        init_val: list[int],
        iterations: int,
) -> list[float]:
    """
    Jacobi Iteration Method:
    An iterative algorithm to determine the solutions of strictly diagonally dominant
    system of linear equations

    4x1 +  x2 +  x3 =  2
     x1 + 5x2 + 2x3 = -6
     x1 + 2x2 + 4x3 = -4

    x_init = [0.5, -0.5 , -0.5]

    Examples:

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_old(coefficient, constant, init_val, iterations)
    [0.909375, -1.14375, -0.7484375]


    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_old(coefficient, constant, init_val, iterations)
    Traceback (most recent call last):
        ...
    ValueError: Coefficient matrix dimensions must be nxn but received 2x3

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_old(
    ...     coefficient, constant, init_val, iterations
    ... )  # doctest: +NORMALIZE_WHITESPACE
    Traceback (most recent call last):
        ...
    ValueError: Coefficient and constant matrices dimensions must be nxn and nx1 but
                received 3x3 and 2x1

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_old(
    ...     coefficient, constant, init_val, iterations
    ... )  # doctest: +NORMALIZE_WHITESPACE
    Traceback (most recent call last):
        ...
    ValueError: Number of initial values must be equal to number of rows in coefficient
                matrix but received 2 and 3

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 0
    >>> jacobi_old(coefficient, constant, init_val, iterations)
    Traceback (most recent call last):
        ...
    ValueError: Iterations must be at least 1
    """

    rows1, cols1 = coefficient_matrix.shape
    rows2, cols2 = constant_matrix.shape

    if rows1 != cols1:
        msg = f"Coefficient matrix dimensions must be nxn but received {rows1}x{cols1}"
        raise ValueError(msg)

    if cols2 != 1:
        msg = f"Constant matrix must be nx1 but received {rows2}x{cols2}"
        raise ValueError(msg)

    if rows1 != rows2:
        msg = (
            "Coefficient and constant matrices dimensions must be nxn and nx1 but "
            f"received {rows1}x{cols1} and {rows2}x{cols2}"
        )
        raise ValueError(msg)

    if len(init_val) != rows1:
        msg = (
            "Number of initial values must be equal to number of rows in coefficient "
            f"matrix but received {len(init_val)} and {rows1}"
        )
        raise ValueError(msg)

    if iterations <= 0:
        raise ValueError("Iterations must be at least 1")

    table: NDArray[float64] = np.concatenate(
        (coefficient_matrix, constant_matrix), axis=1
    )

    rows, cols = table.shape

    strictly_diagonally_dominant(table)

    # Iterates the whole matrix for given number of times
    for _ in range(iterations):
        new_val = []
        for row in range(rows):
            temp = 0
            for col in range(cols):
                if col == row:
                    denom = table[row][col]
                elif col == cols - 1:
                    val = table[row][col]
                else:
                    temp += (-1) * table[row][col] * init_val[col]
            temp = (temp + val) / denom
            new_val.append(temp)
        init_val = new_val

    return [float(i) for i in new_val]


def jacobi_new(
        coefficient_matrix: NDArray[float64],
        constant_matrix: NDArray[float64],
        init_val: list[float],
        iterations: int,
) -> list[float]:
    """
    Jacobi Iteration Method:
    An iterative algorithm to determine the solutions of strictly diagonally dominant
    system of linear equations

    4x1 +  x2 +  x3 =  2
     x1 + 5x2 + 2x3 = -6
     x1 + 2x2 + 4x3 = -4

    x_init = [0.5, -0.5 , -0.5]

    Examples:

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_new(coefficient, constant, init_val, iterations)
    [0.909375, -1.14375, -0.7484375]


    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_new(coefficient, constant, init_val, iterations)
    Traceback (most recent call last):
        ...
    ValueError: Coefficient matrix dimensions must be nxn but received 2x3

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_new(
    ...     coefficient, constant, init_val, iterations
    ... )  # doctest: +NORMALIZE_WHITESPACE
    Traceback (most recent call last):
        ...
    ValueError: Coefficient and constant matrices dimensions must be nxn and nx1 but
                received 3x3 and 2x1

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5]
    >>> iterations = 3
    >>> jacobi_new(
    ...     coefficient, constant, init_val, iterations
    ... )  # doctest: +NORMALIZE_WHITESPACE
    Traceback (most recent call last):
        ...
    ValueError: Number of initial values must be equal to number of rows in coefficient
                matrix but received 2 and 3

    >>> coefficient = np.array([[4, 1, 1], [1, 5, 2], [1, 2, 4]])
    >>> constant = np.array([[2], [-6], [-4]])
    >>> init_val = [0.5, -0.5, -0.5]
    >>> iterations = 0
    >>> jacobi_new(coefficient, constant, init_val, iterations)
    Traceback (most recent call last):
        ...
    ValueError: Iterations must be at least 1
    """

    rows1, cols1 = coefficient_matrix.shape
    rows2, cols2 = constant_matrix.shape

    if rows1 != cols1:
        msg = f"Coefficient matrix dimensions must be nxn but received {rows1}x{cols1}"
        raise ValueError(msg)

    if cols2 != 1:
        msg = f"Constant matrix must be nx1 but received {rows2}x{cols2}"
        raise ValueError(msg)

    if rows1 != rows2:
        msg = (
            "Coefficient and constant matrices dimensions must be nxn and nx1 but "
            f"received {rows1}x{cols1} and {rows2}x{cols2}"
        )
        raise ValueError(msg)

    if len(init_val) != rows1:
        msg = (
            "Number of initial values must be equal to number of rows in coefficient "
            f"matrix but received {len(init_val)} and {rows1}"
        )
        raise ValueError(msg)

    if iterations <= 0:
        raise ValueError("Iterations must be at least 1")

    table: NDArray[float64] = np.concatenate(
        (coefficient_matrix, constant_matrix), axis=1
    )

    rows, cols = table.shape

    strictly_diagonally_dominant(table)

    """
    # Iterates the whole matrix for given number of times
    for _ in range(iterations):
        new_val = []
        for row in range(rows):
            temp = 0
            for col in range(cols):
                if col == row:
                    denom = table[row][col]
                elif col == cols - 1:
                    val = table[row][col]
                else:
                    temp += (-1) * table[row][col] * init_val[col]
            temp = (temp + val) / denom
            new_val.append(temp)
        init_val = new_val
    """

    """
    denom - a list of values along the diagonal
    """
    denom = np.diag(coefficient_matrix)
    """
    val_last - values of the last column of the table array
    """
    val_last = table[:, -1]
    """
    masks - boolean mask of all strings without diagonal
    elements array coefficient_matrix
    """
    masks = ~np.eye(coefficient_matrix.shape[0], dtype=bool)
    """
    no_diag - coefficient_matrix array values without diagonal elements
    """
    no_diag = coefficient_matrix[masks].reshape(-1, rows - 1)
    """
    Here we get 'i_col' - these are the column numbers, for each row
    without diagonal elements, except for the last column.
    """
    i_row, i_col = np.where(masks)
    ind = i_col.reshape(-1, rows - 1)
    """
    'i_col' is converted to a two-dimensional list 'ind',
    which will be used to make selections from 'init_val'
    ('arr' array see below).
    """

    # Iterates the whole matrix for given number of times
    for _ in range(iterations):
        arr = np.take(init_val, ind)
        temp = np.sum((-1) * no_diag * arr, axis=1)
        new_val = (temp + val_last) / denom
        init_val = new_val

    return new_val.tolist()


# Checks if the given matrix is strictly diagonally dominant
def strictly_diagonally_dominant(table: NDArray[float64]) -> bool:
    """
    >>> table = np.array([[4, 1, 1, 2], [1, 5, 2, -6], [1, 2, 4, -4]])
    >>> strictly_diagonally_dominant(table)
    True

    >>> table = np.array([[4, 1, 1, 2], [1, 5, 2, -6], [1, 2, 3, -4]])
    >>> strictly_diagonally_dominant(table)
    Traceback (most recent call last):
        ...
    ValueError: Coefficient matrix is not strictly diagonally dominant
    """

    rows, cols = table.shape

    is_diagonally_dominant = True

    for i in range(0, rows):
        total = 0
        for j in range(0, cols - 1):
            if i == j:
                continue
            else:
                total += table[i][j]

        if table[i][i] <= total:
            raise ValueError("Coefficient matrix is not strictly diagonally dominant")

    return is_diagonally_dominant


for n in [100, 100, 200, 300, 400, 500]:# array size 'coefficient'
    coefficient = np.random.randint(low=-5, high=5, size=(n, n))
    coefficient[np.diag_indices_from(coefficient)] = 10 * n

    constant = np.random.randint(low=-5, high=5, size=(n, 1))

    init_val = np.random.choice([0.5, -0.5], n)

    now = datetime.datetime.now()
    result_old = np.round(jacobi_old(coefficient, constant, init_val, n), 5)
    time_old = datetime.datetime.now() - now

    now = datetime.datetime.now()
    result_new = np.round(jacobi_new(coefficient, constant, init_val, n), 5)
    time_new = datetime.datetime.now() - now

    word = ('array size n x n {0} time_old {1} time_new {2}'
            ' difference in time {3} comparison result array {4}'
            .format(n, time_old.total_seconds(), time_new.total_seconds(),
                    round(time_old / time_new, 2), np.all(result_new == result_old)))

    print(word)



# Test Cases
if __name__ == "__main__":
    import doctest

   doctest.testmod()

Output:

array size n x n 100 time_old 0.89429 time_new 0.018528 difference in time 48.27 comparison result array True
array size n x n 100 time_old 0.873833 time_new 0.020888 difference in time 41.83 comparison result array True
array size n x n 200 time_old 6.827408 time_new 0.138856 difference in time 49.17 comparison result array True
array size n x n 300 time_old 23.260927 time_new 0.474395 difference in time 49.03 comparison result array True
array size n x n 400 time_old 59.593184 time_new 1.084421 difference in time 54.95 comparison result array True
array size n x n 500 time_old 120.983446 time_new 2.177925 difference in time 55.55 comparison result array True

I also tried an array of size 1000: where the calculation time of the old algorithm is
12 minutes, and the new one is 17 seconds.

[rohan472000]

@quant12345 , close this as #8938 is merged