'''
Un developement de la methode du pivot de Gauss basee sur les
operations avec les lignes.  La fonction principale est
  pivotDeGauss
qui prend comme argument une matrice.

Il y a aussi une fonction qui calcule le determinant d'une
matrice (caree).
'''
import numpy as np

def U_ij(n, i, j, x):
    '''
    Retourne la matrice inversible U_ij = I_n+X, ou X a une seule 
    composante non nulle, X[i, j]=x.  
    Remrquez que le produit U_ij*A est la matrice A avec A[i,:],
    i.e. la ligne i, remplacee par A[i,:] + x*A[j,:].
    Des matrices de type U_ij seront utilisees pour modifier A
    en appliquant le pivot de Gauss avec la ligne j comme ligne
    du pivot (le pivot sera un element de cette ligne). 
    
    L'element x dependra de la matrice A qu'on voudra echelonner.
    '''
    P = np.eye(n)
    P[i, j] = x
    return P


def lElement_x(A, posPivot, ligne):
    '''
    A est la matrice qu'on veut echelonner ; posPivot est la
    position du pivot, i.e. une liste [i, j] ; le troisieme argument
    ligne (un entier) designe la ligne de A qui sera modifiee par
    un produit a gauche du type U*A.  
    '''
    lPivot, cPivot = posPivot
    pivot = A[lPivot, cPivot]
    return -A[ligne, cPivot]/pivot
    

def utilisationDUnPivot(A, posPivot):
    '''
    A est une matrice et posPivot la position d'un pivot.  La
    fonction retourne la matrice obtenue en utililsant la ligne du pivot.
    Le coefficient x avec lequel une ligne i sera modifiee en
    appliquant l'algorithme du pivot est renvoye par la fonction
    lElement_x.

    Remarque : Le pivot n'est pas change !
    '''
    lPivot, cPivot = posPivot
    n = A.shape[0]
    for i in range(n):
        if i != lPivot:
            U = U_ij(n, i, lPivot, lElement_x(A, posPivot, i))
            A = U.dot(A)
            # c = np.shape(M)
    return A


def maxEtPosition(A_2darray):
    '''
    Retourne un element maximal et sa position dans A_2darray.  L'argument
    A_2darray est un tableau ayant deux dimensions (i.e. une matrice).
    C'est cette position qui deviendra le prochain pivot dans le
    developement de l'algorithme de Gauss.
    '''
    indices = np.where(A_2darray == A_2darray.max())
    # print(indices, "\n")  # il faut comprendre la sortie de np.where !
    position = [indices[0][0], indices[1][0]]
    return A_2darray[position[0], position[1]], position


def zeroSurLignesEtColonnes(A, positions_list):
    '''
    Remplace par des zeros les elements correspondant aux lignes et 
    colonnes des positions de positions_list.  Enleve les signes des
    autres termes.
    '''
    B = abs(A)
    for pos in positions_list:
        B[pos[0], :] = 0
        B[:, pos[1]] = 0
    return B


def pivotDeGauss(A, eps=.1**14):
    '''
    Le pivot de Gauss en utilisant les lignes.  La fonction suit
    les operations qu'un humain ferrait.
    Le deuxieme argument controle l'egalite avec 0.
    La fonction renvoie la matrice "echelonnee" (sans changer la
    position des pivots lors de la construction) ainsi que la liste
    des positions des pivots utilises.
    '''
    posPivots_list = []
    tmp_array = abs(A)
    a, pPivot = maxEtPosition(tmp_array)
    while a>eps:  # a is always >= 0
        posPivots_list.append(pPivot)
        A = utilisationDUnPivot(A, pPivot)
        tmp_array = zeroSurLignesEtColonnes(A, posPivots_list)
        a, pPivot = maxEtPosition(tmp_array)
    return A, posPivots_list


def determinant(A):
    E, posPivots_list = pivotDeGauss(A)
    if np.shape(A)[0] != np.shape(A)[1]:
        return False
    if len(posPivots_list)==np.shape(A)[0]:
        d = 1.
        for p in posPivots_list:
            d = d*(-1)**(p[0]+p[1])*E[p[0], p[1]]
    else:
        d = 0.
    return d


def montrerMatrice(A, eps=.1**8):
    '''
    Transforme les coefficients dont la valeur absolue est <eps en 0
    et montre la matrice dans le terminal.
    '''
    m, n = np.shape(A)
    C = np.zeros((m, n))
    for i in range(m):
        for j in range(n):
            if abs(A[i, j])>eps:
               C[i, j] =  A[i, j]
    print(C)


####################################################################
####################################################################
### tests numeriques


A = np.array([[1, 2, 3, -4], [5, -7, 9, 11], [0, 0, 0, 2], [0, -1, 0, 3]])
B = np.array([[1, -8, 3], [-1, 4, 8], [0, 0, 0]])
C = np.array([[2, -1, 0, 1, 0, 0], [-1, 2, -1, 0, 1, 0], [0, -1, 2, 0, 0, 1]])

M = B

E, listePositionsPivots = pivotDeGauss(M)
print(f"\npour la matrice \n{M} \nune forme echelonee est")
montrerMatrice(E)

print(f"\nde plus sont determinant est {determinant(M)}")
print(np.linalg.det(M), " est le dereminant calcule par linalg")