__author__ = 'true'
import numpy as np
import scipy as sp


def lu_decomposition(a):
    m, n = a.shape
    l = np.zeros((m, m))
    u = np.zeros((m, m))
    for i in range(n):
        for j in range(n):
            u[i][j] = a[i][j]

    for i in range(n):
        for j in range(i, n):
            try:
                l[j][i] = u[j][i] / u[i][i]
            except ZeroDivisionError:
                print("Деление на ноль!")

    for k in range(1, n):
        for i in range(k - 1, n):
            for j in range(i, n):
                try:
                    l[j][i] = u[j][i] / u[i][i]
                except ZeroDivisionError:
                    print("Деление на ноль")
        for i in range(k, n):
            for j in range(k - 1, n):
                u[i][j] -= l[i][k - 1] * u[k - 1][j]

    return l, u


def lu_solve(a, b):
    if a.shape[0] != len(b):
        raise ValueError("Размеры входных массивов не совпадают!")
    l, u = lu_decomposition(a)
    n = a.shape[0]
    y = np.zeros(n)
    for i in range(n):
        y[i] = b[i]
        for j in range(i):
            y[i] -= y[j] * l[i][j]
    x = np.zeros(n)
    for i in reversed(range(n)):
        x[i] = y[i]
        for j in range(i + 1, n):
            x[i] -= u[i][j] * x[j]
        x[i] /= u[i][i]
    return x


# A = np.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]])
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

Q, R = np.linalg.qr(A)

print("numpy QR")
print("Q:")
print(Q)
print("R:")
print(R)


def house(x):
    n = len(x)
    mu = np.linalg.norm(x, 2)
    result = np.zeros(n)
    np.copyto(result, x)
    if mu != 0:
        beta = x[0] + np.sign(x[0]) * mu
        result[1:n] = result[1:n] / beta
    result[0] = 1
    print("HOUSE:", result)
    return result


def row_house(a, v):
    beta = -2 / np.dot(v.T, v)
    print("BETA", beta)
    w = beta * np.dot(a.T, v)
    print("W:", w)
    print("ADDITION:", np.dot(v, w.T))
    return a + np.dot(v, w.T)


m, n = A.shape

V = np.zeros(m)
for j in range(n):
    V[j:m] = house(A[j:m, j])
    A[j:m, j:n] = row_house(A[j:m, j:n], V[j:m])
    if j < m:
        A[j + 1:m, j] = V[j + 1:m]
    print("j:",j)
    print("A:")
    print(A)

print("\nГолуб:")
print(A)