#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Jun  7 12:09:51 2024

@author: mbpadmin
"""

import numpy as np
import matplotlib.pyplot as plt

## Exercice 1

# Q2. 
def matrice_fibo(n):
    A = np.zeros((n, n))
    a, b = 0, 1
    for k in range(1, n): # on remplit les indices i+j <= n
        for i in range(k+1):
            A[i, k-i] = A[k-i, i] = b
        a, b = b, a+b
    for k in range(n, 2*n-1): # on remplit les indices i+j > n
        for i in range(k-n+1, n):
            A[i, k-i] = A[k-i, i] = b
        a, b = b, a+b   
    return A

def matrice_fibo2(n):
    F = [0, 1]
    for k in range(2, 2*n - 1):
        F.append(F[-2] + F[-1])
    A = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            A[i, j] = F[i + j]
    return A

"""
Il est aussi possible d'écrire une première fonction calculant la suite de 
Fibonacci et une seconde appelant celle-ci pour remplir la matrice. Si l'on 
procède ainsi, il faut renvoyer la liste (ou le dictionnaire) des valeurs de 
F[k] pour toutes les valeurs de k nécessaires et non simplement F[n], ce qui 
ferait recalculer plein de fois tous les éléments de la suite depuis le départ.

Bien sûr, il faut en priorité éviter une complexité exponentielle en utilisant
naïvement la définition récursive de F[n].

Le plus simple est de procéder linéairement (fonction fibo). On peut aussi 
utiliser la mémoïsation (fonction fibo_memo), mais c'est plus long à écrire et
il faut être sûr de ne pas se tromper.
"""

def fibo(n): 
    dic = {0:0, 1:1}
    for i in range(2, n+1):
        dic[i] = dic[i-1] + dic[i-2]
    return dic

def fibo_memo(n):
    dic = {0:0, 1:1}
    def remplissage(k):
        if k in dic:
            return dic[k]
        dic[k] = remplissage(k-1) + remplissage(k-2)
        return dic[k]
    remplissage(n)
    return dic

def matrice_fibo3(n):
    dic = fibo(2*n - 2)
    M = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            M[i, j] = dic[i + j]
    return M

A, B, C = matrice_fibo(10), matrice_fibo2(10), matrice_fibo3(10)
assert np.array_equal(A, B)
assert np.array_equal(B, C)

"""
Python n'affiche rien à la compilation du assert: il est content (donc 
l'affirmation est juste).
"""
      
# Q3.
def spectre(n):
    return np.linalg.eigvals(matrice_fibo(n)) 


# Q4.
R = [np.linalg.matrix_rank(matrice_fibo(n)) for n in range(2, 10)]

def vp_nonnulles(n):
    sp = spectre(n)    
    return np.sum((abs(sp) > 1e-8))

Alpha, Beta = [], []
for n in range(2, 30):
    sp = spectre(n).real 
    sp.sort() # tri de la liste des valeurs propres
    Alpha.append(sp[0])
    Beta.append(sp[-1])



"""
Exercice 2
"""

# Q1.
def fadeev(A):
    p = A.shape[0] # récupère la taille de la matrice A
    k, Ak, L = 1, A, [1]
    while k <= p:
        alphak = np.trace(Ak)/k
        L.append(- alphak)
        Ak = A.dot(Ak - alphak * np.eye(p))
        k += 1
    return L

# Q2   
A = np.array([1, -4, 1, 1, -1, 1, 0, 1, 0, 0, -2, 1, 0, 0, -1, 3]).reshape(4, 4)
p1, p2 = fadeev(A), np.poly(A)

"""
Expérience sur les références partagées

B = A
B = B.dot(A)
print(f"A = {A}"), print(f'B = {B}')

B = A
C = A.copy()
B[0, 0] = -1
print(f'B = {B}'), print(f'A = {A}')
C = A.copy()
C[0,0] = -2
print(f'A = {A}'), print(f'C = {C}')
"""



"""
Exercice 3
"""
# Question 1
def diagonale_dominante_stricte(M):
    n = np.shape(M)[0]
    for i in range(n):
        if 2 * abs(M[i, i]) <= np.sum(abs(M[i, :])):
            return False
    return True

def diagonale_dominante_stricte2(M):
    n = np.shape(M)[0]
    test, i = True, 0
    while test and i < n:
        m, S, j = 2 * abs(M[i, i]), 0, 0
        while S < m and j < n:
            S += abs(M[i, j])
            j += 1
        test = (S < m)
        i += 1
    return test

# Question 2
A = np.array([5, -1, 2, 1, -4, -2, 3, -1, 5]).reshape(3, 3)
X = np.linalg.solve(A, np.array([1, 1, 1]))

"""
Dans la console : 
diagonale_dominante_stricte(A)
X
A.dot(X)
"""

# Question 6
def methode_iterative(A, B, epsilon):
    n = np.shape(A)[0]
    Dinv = np.diag([1/A[i, i] for i in range(n)])
    AA = abs(A)
    g = max([np.sum(AA[i, :])/AA[i, i] - 1 for i in range(n)])
    assert g < 1, "A n'est pas à diagonale strictement dominante" 
    def une_iteration(X):
        return (np.eye(n) - np.dot(Dinv, A)).dot(X) + Dinv.dot(B)
    X = np.array([1, 0, 0])
    Y = une_iteration(X)
    p = np.log(epsilon*(1-g) / np.max(abs(Y - X))) / np.log(g)
    for k in range(1, int(p)):
        Y = une_iteration(Y)
    return Y



"""
Exercice 4
"""
# Question 3
def GramSchmidt(base):
    BON, dim = [], len(base)
    for k in range(len(base)):
        E, Xk = np.zeros(dim), base[k]
        for i in range(k):
            E += np.vdot(Xk, BON[i]) * BON[i]
        V = Xk - E
        BON.append(V / np.linalg.norm(V)) # ou V / np.sqrt(np.vdot(V, V))
    return BON

B = [np.random.random(100) for _ in range(100)]
BON = GramSchmidt(B)
Gram = np.array([[np.vdot(BON[i], BON[j]) for i in range(100)] for j in range(100)])

"""
Dans la console :
np.amax(abs(Gram - np.eye(100)))
"""

# Question 4
def QR(A):
    """
    On n'écrit pas de documentation à l'oral !
    
    Entrée : une matrice inversible 
    Sortie : un couple (Q, R) avec Q orthogonale et R triangulaire supérieure
             à diagonale strictement positive telle que A = QR.
    """
    n = np.shape(A)[0]
    BON = GramSchmidt([A[:, k] for k in range(n)])
    R = np.zeros((n, n))
    for i in range(n):
        for j in range(i, n):
            R[i, j] = np.vdot(BON[i], A[:, j])
    return np.array(BON).T, R


def verif(Q, R):
    n, p = R.shape
    q, r = Q.shape
    if not n == p == q == r:
        raise ValueError("Pour Centrale, c'est pas dans la poche...")
    Zero = Q.dot(Q.T) - np.eye(Q.shape[0])
    if np.amax(abs(Zero)) > 1e-9 or R[0, 0] <= 0:
        return False
    for i in range(1, n):
        if R[i, i] <= 0:
            return False
        for j in range(i):
            if abs(R[i, j]) > 1e-9:
                return False
    return True   

A = np.random.random((20, 20)) # test sur une matrice aléatoire
Q, R = QR(A)
print(verif(Q, R))


"""
Exercice 5
"""
# Question 3.a.
from numpy.polynomial import Polynomial as pol

def matrice(F, G):
    d = G.degree()
    M = np.zeros((d, d))
    Can = [pol([0] * k + [1]) for k in range(d + 1)]
    for k in range(d):
        img = ((F * Can[k]) % G).coef
        M[: len(img), k] = img
    return M

F, G = pol([0, 1]), pol([0, -2, -1, 2, 1])
A = matrice(F, G)
P = np.poly(A)

# Question 3.b.
vp = np.linalg.eig(A)[1]
L = []
plt.figure()
X = np.linspace(-3, 2, 50)
plt.grid()
for j in range(4):
    f = pol(vp[:, j])
    L.append(f.roots())
    plt.axis([-3.1, 2.1, -6, 6])
    plt.title(f'Les vecteurs propres de $g$')
    plt.plot(X, f(X))
plt.show()


# Question 3.c.
FF = [pol([1]), pol([-1, 1]), pol([2,3,2]), pol([1,3,-1,1]), pol([2,1,0,1,4])]
vec_pro = []
for F in FF:
    vp = np.linalg.eig(matrice(F, G))[1]
    L = []
    for j in range(4):
        f = pol(vp[:, j])
        L.append(f.roots())
    vec_pro.append(L)


"""
Exercice 6
"""