#!/usr/bin/python3
#
# Fonctions utiles pour la décomposition QR

import numpy as np

# retourne la norme du vecteur v
def norme(v):
    return np.sqrt(np.sum(v**2))

# retourne la projection du vecteur w sur le vecteur v
def projeter(v, w, epsilon=1.0E-10):
    vnorme = norme(v)
    if vnorme < epsilon:   # projection sur le vecteur nul = vecteur nul par déf.
        return np.zeros(v.shape[0])
    resultat = np.copy(v)  # le vecteur résultat
    resultat *= v @ w      # multiplier par le produit scalaire v.w
    resultat /= vnorme**2  # diviser par |v|²
    return resultat

# retourne une matrice dont les colonnes sont les colonnes orthonormalisées de a
def gram_schmidt(A, epsilon=1.0E-10):
    n = A.shape[0]
    resultat = np.copy(A)       # matrice pour mettre les résultats
    for i in range(n):          # orthonormalisation de Gram-Schmidt modifiée:
        v = resultat[:, i]           # i-ème colonne
        vnorme = norme(v)         
        if vnorme >= epsilon:    #  normaliser (sauf les vecteurs zéro)
            v /= vnorme
        for j in range(i+1, n): # des vecteurs suivants, soustraire la projection
            resultat[:, j] -= projeter(v, resultat[:, j], epsilon)
    return resultat

# retourne q et r tel que q est orthogonale, r est triangulaire supérieure, et qr=a
def decomp_QR(A, epsilon=1.0E-10):
    Q = gram_schmidt(A, epsilon)
    R = Q.T @ A
    return Q, R
    
# tester    
if __name__ == '__main__':
    A = np.array(range(16), dtype=float).reshape((4,4))
    n = A.shape[0]
    for i in range(n): # symétriser
        for j in range(i, n):
            A[i, j] = A[j, i]
    print(A)
    print(gram_schmidt(A))
    Q, R = decomp_QR(A)
    print(Q)
    print(R)
    print(Q @ R)