Shit and cum

piss_and_cry.py

import numpy as np


def gauss(A, B):
    n = len(A)
    
    def argmax(A):
        i, m = 0, A[0]
        for j, a in enumerate(A):
            if a > m:
                i, m = j, a
        return i

    def forward_elimination(M):
        for i in range(n):
            column_max = argmax(A[i:,i]) + i # pick max element in fixed column within all rows
            M[[i, column_max]] = M[[column_max, i]] # swap first and column_max rows
            M[i] /= M[i, i] # normalize first row
            for j in range(i + 1, n):
                M[j] -= M[i] * M[j, i] # eliminate i-th variable
        return M

    def back_substitution(M):
        for i in range(n-2, -1, -1):
            for j in range(i + 1, n):
                M[i] -= M[j] * M[i, j]
        return M
        
    B = np.expand_dims(B, axis=0) # necessary for concatenate
    M = np.concatenate((A, B.T), axis=1)
    return back_substitution(forward_elimination(M))[:, -1] # return answer only

def simple_iterative(A, b, epsilon):
    def norm_1_matrix(A):
        n = len(A)
        return max([sum([abs(A[i, j]) for i in range(n)]) for j in range(n)])

    def norm_1_vec(v):
        return sum([abs(x) for x in v])

    def is_positive_definite(A):
        try:
            np.linalg.cholesky(A)
        except np.linalg.LinAlgError:
            return False
        return True

    if not is_positive_definite(A):
        A, b = np.dot(A.T, A), np.dot(A.T, b)
    mu = 1 / norm_1_matrix(A)
    B, c = np.identity(len(A)) - mu*A, mu*b
    x0 = c
    while 1:
        x1 = np.dot(B, x0) + c
        if norm_1_vec(x0 - x1) < epsilon:
            return x1
        x0 = x1

# bad conditioned matrix, simple iterative should have a very poor accuracy
def bad(sz, eps, n):
    assert sz > 0
    val = eps * n
    A = np.full((sz, sz), val)
    for i in range(sz):
        A[i, i] += 1
        for j in range(i + 1, sz):
            A[i, j] = -1 - val
    b = np.full(sz, -1)
    b[-1] = 1
    return A, b

# well conditioned matrix, simple iterative should give a great accuracy
def good(n):
    A = np.array([
            [n+2, 1, 1],
            [1, n+4, 1],
            [1, 1, n+6]
        ], 'float')
    b = np.array([n+4, n+6, n+8])
    return A, b


def pretty_print(A, b, eps):
    print("Matrix:", A, sep='\n')
    print("Coefficient column:", b)
    print("Gaussian method:", gauss(A, b))
    print("Simple iteration method: ", simple_iterative(A, b, eps))
    print()

        
def test():
    n, eps, = 16, 10**(-6)

    A, b = good(n)
    pretty_print(A, b, eps)
    
    for i in range(2, 6):
        A_, b_ = bad(i, 10**(-1 - i), n)
        pretty_print(A_, b_, eps)


if __name__ == "__main__":
    test()