import numpy as np
from ex1 import resolve_up, resolve_down


def est_symetrique_definie_positive(A) -> bool:
    """
    Vérifie que A = A^t et que toutes les valeurs propres soit positive
    :param A:
    :return:
    """

    # Symétrie (avec tolérance pour les flottants)
    sym = np.allclose(A, A.T)
    # Définie positive : toutes les valeurs propres > 0
    valeurs_propres = np.linalg.eigvals(A)
    definie_positive = np.all(valeurs_propres > 0)
    return sym and definie_positive


def gauss(A, b):
    """
    Pas de condition, mais pas optimal non plus.
    :param A: Une matrice
    :param b: Un vecteur
    :return: A triangularisée et b modifié en conséquence.
    """
    n = len(A)
    for k in range(n): # k = pivot
        for i in range(k+1, n): # i = ligne en dessous du pivot
            g = A[i, k] / A[k, k] # multiplication pour éliminer A[i, k]
            A[i, k:] -= g * A[k, k:] # row operation on A
            b[i] -= g * b[k] # same row operation on b

    return A, b # matrice triangularisée et le vecteur b modifié



def cholesky(A):
    if not est_symetrique_definie_positive(A):
        raise ValueError('Matrice non symétrique définie positive')

    m = A.shape[0]  # Taille de la matrice (m x m)
    L = np.zeros_like(A, dtype=np.float32)  # Matrice L initialisée à zéro

    for j in range(m):  # Parcourt chaque colonne j de 0 à m-1
        for i in range(j, m):  # Parcourt chaque ligne i de j à m-1
            s = sum(L[i, k] * L[j, k] for k in range(j))  # Somme des produits L[i,k] * L[j,k] pour k < j

            if i == j:  # Élément diagonal
                L[i, j] = np.sqrt(A[i, i] - s)  # Calcul de la racine carrée pour la diagonale
            else:  # Élément hors diagonale
                L[i, j] = (A[i, j] - s) / L[j, j]  # Calcul des éléments hors diagonale
    return L


if __name__ == '__main__':
    A = np.array([
        [1, 1, 1, 1],
        [1, 5, 5, 5],
        [1, 5, 14, 14],
        [1, 5, 14, 15]
    ], dtype=np.float64)

    b = np.array([1, 0, 0, 0], dtype=np.float64)

    A_triangle, b_gauss = gauss(A, b)
    print(f"A triangulsarisée:\n {A_triangle}")

    x = resolve_up(A_triangle, b_gauss)
    print(f"x={x}")

    A = np.array([
        [4, 1, 1, 1],
        [1, 5, 2, 1],
        [1, 2, 6, 2],
        [1, 1, 2, 7]
    ])

    L = cholesky(A)
    print(f"L: {L}")
    x = resolve_down(L, b)
    print(f'x: {x}')