import numpy as np

# Исходные данные

As = [
    np.array([
        [-7.8, -2.4, -0.3, -3.1],
        [-8.0,  4.5, -9.2,  2.3],
        [-4.8,  5.5,  1.7, -9.5],
        [-2.5,  0.3, -8.8,  7.5]
    ]),
    np.array([
        [-8.5, -5.1, -1.3, -6.8],
        [ 0.7, -9.4,  3.1,  1.7],
        [-9.2, -2.2,  7.7, -4.9],
        [-3.2, -5.4,  3.2,  8.2]
    ]),
    np.array([
        [-6.6, -8.8,  6.6, -8.1],
        [-5.7, -9.3, -1.4,  5.2],
        [-1.0, -1.2,  4.1, -5.9],
        [ 4.3,  0.8, -7.2, -6.6]
    ])
]

x = np.array([1, 2, 3, 4])
bs = [A.dot(x) for A in As]


# вспомогательные функции

def _get_column_base(M, col):
    """Ищём главный элемент столбца"""
    base_i = 0
    base = 0
    for row in range(col, M.shape[0]):
        if abs(base) < abs(M[row][col]):
            base = M[row][col]
            base_i = row
    return base, base_i


def _swap_rows(M, i, j):
    if i == j:
        return
    temp = np.copy(M[i])
    M[i] = M[j]
    M[j] = temp
    # M[[i, j]] = M[[j, i]]


def _permutate_b(P, b):
    b = np.copy(b)
    n = len(P)
    for i in range(n):
        b[i], b[P[i]] = b[P[i]], b[i]
    return b


# основные функции

def LUP_factor(A):
    M = np.copy(A)
    n = M.shape[0]
    P = np.arange(n)
    for col in range(n - 1):
        base, base_i = _get_column_base(M, col)
        _swap_rows(M, base_i, col)
        P[col] = P[base_i]
        for row in range(col + 1, n):
            k = M[row, col] / base
            M[row, col] = k
            for i in range(col + 1, n):
                M[row, i] -= M[col, i] * k
    return M, P


"""
1) Решить заданную систему линейных алгебраических уравнений
методом LU разложения с выбором главного элемента по столбцу.
Ax=b
"""

print('Задание №1', end="\n---\n")


def solve_LU(M, b):
    LU, P = LUP_factor(M)
    b = _permutate_b(P, b)
    n = M.shape[0]
    y = np.zeros((n,))
    x = np.zeros((n,))
    # TODO Is it possible to combine x, y?

    # Ly = b => y
    for i in range(n):
        left_sum = 0
        for j in range(i):
            left_sum += LU[i][j] * y[j]
        y[i] = b[i] - left_sum

    # Ux = y => x
    for i in range(n - 1, -1, -1):
        left_sum = 0
        for j in range(i + 1, n):
            left_sum += LU[i][j] * x[j]
        x[i] = (y[i] - left_sum) / LU[i][i]
    return x


for i, (A, b) in enumerate(zip(As, bs), start=1):
    print(f'Решение системы №{i}: ', solve_LU(A, b))

"""
2) Используя полученное ранее LU разложение вычислить обратную матрицу к исходной.
"""

print('\nЗадание №2', end="\n---\n")


def inverse(M):
    n = M.shape[0]
    cols = []
    for k in range(n):
        b = np.eye(1, n, k=k)[0]  # правая часть, где x_i = 0, кроме x_k = 1
        cols.append(solve_LU(M, b))
    M_inv = np.array(cols).transpose()
    return M_inv


for i, A in enumerate(As, start=1):
    print(f'Обратная матрица для A_{i}:')
    print(inverse(A), end='\n\n')


"""
3) Вычислить определитель и найти число обусловленности исходной матри-цы
в кубической, октаэдрической и евклидовой нормах.
"""

print('\nЗадание №3', end="\n---\n")


def _sign_det(P):
    P = np.copy(P)
    pc = 0
    for i in range(len(P)):
        if P[i] != i:
            pc += 1
            P[P[i]] = P[i]
    return (1, -1)[pc % 2]


def det_LU(LU, P):
    n = LU.shape[0]
    det = LU[0, 0]
    for i in range(1, n):
        det *= LU[i, i]
    return det * _sign_det(P)


def C_norm(M):  # np.linalg.norm(A1, ord=1)
    max_el = 0
    n = M.shape[0]
    for j in range(n):
        el_sum = 0
        for i in range(n):
            el_sum += abs(M[i][j])
        if el_sum > max_el:
            max_el = el_sum
    return max_el


def O_norm(M):  # np.linalg.norm(A1, ord=np.inf)
    max_el = 0
    n = M.shape[0]
    for i in range(n):
        el_sum = 0
        for j in range(n):
            el_sum += abs(M[i][j])
        if el_sum > max_el:
            max_el = el_sum
    return max_el


def E_norm(M):  # np.linalg.norm(A1, ord='fro')
    s = 0
    n = M.shape[0]
    for i in range(n):
        for j in range(n):
            s += M[i][j]**2
    return np.sqrt(s)


def cond(A, norm):
    A_inv = inverse(A)
    return norm(A) * norm(A_inv)


for i, A in enumerate(As, start=1):
    LU, P = LUP_factor(A)
    print(f'det(A{i}) =', det_LU(LU, P))
    print(f'cond_C(A{i}) =', cond(A, norm=C_norm))
    print(f'cond_O(A{i}) =', cond(A, norm=O_norm))
    print(f'cond_E(A{i}) =', cond(A, norm=E_norm))
    print()