# -*- coding: utf8 -*-
import numpy as np
import time

#metoda eliminacji Gaussa
#www.algorytm.org

#definiowanie przykładowych danych

dataset = 2
if dataset==1:
    A = np.matrix("10, -7, 0;0, -0.1, 6; 0, 2.5, 5")
    b = np.matrix("6;5.8;0")
    x = np.matrix([-22/31.,-58/31.,29/31.]).T
elif dataset==2:
    A = np.matrix("10, -7, 0, 5, 12, 1;11, -0.1, 6, 0, 1, 0; 0, 21.5, 5, 12, 0, 6; 0, 21.5, 3, 12, 2, 6; 0, 14.5, 5, 12, 8, 6; 9, 1.5, 1, 2, 4, 1")
    b = np.matrix("6;5.8;0;6;1;5")
    x = np.matrix([-22/31.,-58/31.,29/31.]).T

print A
print b

def gauss(_A, _b):
    """ Najprostsza metoda gaussa. """
    n = _A.shape[0]
    #eliminacja zmiennych
    for i in range(0,n-1): #dla każdej kolumny
        for k in range(i+1,n): #wszyztkie wiersze
            l=(_A[k,i]/_A[i,i])
            _A[k] = _A[k]-_A[i]*l
            _b[k] = _b[k]-_b[i]*l
    #postępowanie odwrotne TODO to nie jest ładne. To nie jest wydajne.
    _x = [0]*n
    _x[n-1] = _b.item(n-1)/_A.item(n-1,n-1)
    for k in range(n-2,-1,-1):
        _x[k] = (_b[k]-sum(_A[k]*np.matrix(_x).T)).item(0)/_A.item(k,k) #TODO na pewno da się to zrobić inaczej. Nie jest to wydajne.
    return np.matrix(_x).T

def swap_rows(_A, _b, rows):
    """metoda zamienia ze sobą dwa wiersze w układzie. Działa tylko dla 2 wymiarów."""
    for mat in (_A, _b):
        tmp = mat[rows[0]].copy()
        mat[rows[0]] = mat[rows[1]]
        mat[rows[1]] = tmp


def swap_column(_A, _b):
    """metoda zamienia ze sobą dwa wiersze w układzie. Niezaimplementowano."""
    pass
    
def gauss_stable(_A, _b):
    """ Gauss z zaimplementowanym częściowym wyborem elementu głównego. """
    n = _A.shape[0]
    #eliminacja zmiennych
    for i in range(0,n-1): #jedziemy po przekątnej.
        maxrow = i + abs(_A.T[i,i:]).argmax() #prawdopodobnie wąskie gardło metody.
        currow = i
        if maxrow != currow: swap_rows(_A,_b,(maxrow, currow))
        for k in range(i+1,n): #wszyztkie wiersze
            l=(_A[k,i]/_A[i,i])
            _A[k] = _A[k]-_A[i]*l
            _b[k] = _b[k]-_b[i]*l
    #postępowanie odwrotne TODO to nie jest ładne. To nie jest wydajne.
    _x = [0]*n
    _x[n-1] = _b.item(n-1)/_A.item(n-1,n-1)
    for k in range(n-2,-1,-1):
        _x[k] = (_b[k]-sum(_A[k]*np.matrix(_x).T)).item(0)/_A.item(k,k) #TODO na pewno da się to zrobić inaczej. Nie jest to wydajne.
    return np.matrix(_x).T

def get_time(f, count):
    t0 = time.clock()
    for i in range(0, count):
        f()
    return time.clock() - t0

print "START"
print gauss(A.copy(),b.copy())

if __name__ == "__main__":
    print "-----rozwiązanie ze skryptu:"
    print x
    print "-----rozwiązanie wyliczone w sposób bezpośredni:"
    print np.linalg.inv(A)*b
    #~ print "Czy rozwiązania się pokrywają: ", True if sum(x-np.linalg.inv(A)*b).A1 < 1e-10 else False
    print "-----rozwiązanie z gaussa:"
    print gauss(A.copy(),b.copy())
    print "Czy rozwiązania się pokrywają: ", True if sum(gauss(A.copy(),b.copy())-np.linalg.inv(A)*b).A1 < 1e-10 else False
    print "-----rozwiązanie z gaussa z częściowym wyborem elementu głównego:"
    print gauss_stable(A.copy(),b.copy())
    print "Czy rozwiązania się pokrywają: ", True if sum(gauss_stable(A.copy(),b.copy())-np.linalg.inv(A)*b).A1 < 1e-10 else False

    print "----testowanie czasowe:"
    print "numpy odwrócenie macierzy A:", get_time(lambda :np.linalg.inv(A)*b, 10000)
    print "numpy solve A:", get_time(lambda :np.linalg.solve(A,b), 10000)
    print "liczenie moim gaussem:", get_time(lambda :gauss(A.copy(),b.copy()), 10000)
    print "liczenie moim gaussem stabilnym:", get_time(lambda :gauss_stable(A.copy(),b.copy()), 10000)

    print "co niestety muszę napisać... mój gauss jest teraz niestabilny numerycznie, niewydajny i brzydki w kodzie. Zrobiłem to źle i przyznaje się do tego." + """
        dla A[6,6]:
        solve 0.31s
        gauss 14.8s
        gauss_stable 18.9
        ciekaw jestem czy solve nie czituje przez jakąś pamięć podręczną.
    """