###########################################
###########################################
## TP 04 : Calcul matriciel
###########################################
###########################################

## Importations
import numpy as np
import numpy.linalg as la
import random as rd

###########################################
## Exercice 1 : la fonction np.array
###########################################

A = np.array([[1,2,3,4],[11,12,13,14],[21,22,23,24]])
assert A.shape == (3,4)

A2 = A.reshape(2,6)
A3 = A.reshape(6,2)
# print(A)
# print(A2)
# print(A3)

# # 'reshape' est possible si on conserve le nombre de coefficients :
# # A.reshape(n,p) fonctionne ssi n*p = A.size

B = np.array([[1,11],[2,12],[3,13],[4,14],[5,15]])
assert B.shape == (5,2)

B2 = B.reshape(2,5)
B3 = B.reshape(10,1)

L = np.array([[1,2,3,4,5]])
C = np.array([[6],[7],[8],[9],[10]])
assert L.shape == (1,5)
assert C.shape == (5,1)

###########################################
## Exercice 2 : coefficients Aij définis par des formules
###########################################

# avec des boucles 'for' imbriquées
A = np.zeros((3,4))
for i in range(3) :
    for j in range(4) :
        A[i,j] = i+2*j
print('\n A = ')
print(A)

B = np.zeros((3,3))
for i in range(3) :
    for j in range(3) :
        B[i,j] = i+j
print('\n B = ')
print(B)

C = np.zeros((4,4))
for i in range(4) :
    for j in range(4) :
        C[i,j] = i*j
print('\n C = ')
print(C)

D = np.zeros((3,4))
for i in range(3) :
    for j in range(4) :
        D[i,j] = i-j
print('\n D = ')
print(D)

E = np.zeros((5,5))
for i in range(5) :
    for j in range(5) :
        E[i,j] = i**2+j**2
print('\n E = ')
print(E)

F = np.zeros((4,3))
for i in range(4) :
    for j in range(3) :
        F[i,j] = (-1)**(i+j)*(i+j)
print('\n F = ')
print(F)

# avec des listes définies 'en compréhension'
A2 = np.array([[i+2*j for j in range(4)] for i in range(3)])
B2 = np.array([[i+j for j in range(3)] for i in range(3)])
C2 = np.array([[i*j for j in range(4)] for i in range(4)])
D2 = np.array([[i-j for j in range(4)] for i in range(3)])
E2 = np.array([[i**2+j**2 for j in range(5)] for i in range(5)])
F2 = np.array([[(-1)**(i+j)*(i+j) for j in range(3)] for i in range(4)])

###########################################
## Exercice 3 : extraction de coefficients, lignes ou colonnes
###########################################

A = np.array([[i+2*j for j in range(5)] for i in range(4)])

a12 = A[0,1] # extraction d'un coefficient
L1 = A[:1] # extraction d'une ligne
Cn = np.transpose(np.array([A[:,-1]])) # extraction d'une colonne
# remarque : A[:,j] renvoie la liste des coefficients de la colonne d'indice j de A

# print(A)
# print(a12)
# print(L1)
# print(Cn)

sousA = A[0:3] # ne conserve dans A que les lignes d'indices de 0 à 2
tsousA = np.transpose(sousA) # on transpose pour pouvoir sélectionner des colonnes
exA = tsousA[1:4] # ne conserve que les colonnes d'indices 1 à 3
extraitA = np.transpose(exA) # on retranspose

# print('A =')
# print(A)
# print('\n Sous-matrice = ')
# print(extraitA)

B = np.array([[i**2+2*j+1 for j in range(5)] for i in range(5)])

Diag = []
for i in range(5) :
    Diag.append(B[i,i])

overDiag = []
for i in range(4) :
    overDiag.append(B[i,i+1])

# print('B =')
# print(B)
# print('\n éléments diagonaux :')
# print(Diag)
# print('\n éléments au-dessus de la diagonale :')
# print(overDiag)

C = np.array([[2*i + 3*j + 1 for j in range(4)] for i in range(6)])
L2 = C[1]
L4 = C[3]
C1 = np.transpose(np.array([C[:,0]]))
C3 = np.transpose(np.array([C[:,2]]))

# print('C =')
# print(C)
# print('\n ligne 2 de C :')
# print(L2)
# print('\n ligne 4 de C :')
# print(L4)
# print('\n colonne 1 de C')
# print(C1)
# print('\n colonne 3 de C')
# print(C3)

D = np.array([[2*i - 3*j + 1 for j in range(7)] for i in range(7)])

Dl = np.array([D[2*i] for i in range(4) ])
tDl = np.transpose(Dl)
tDlc = np.array([tDl[2*i] for i in range(4)])
tamisD = np.transpose(tDlc)

# print('D =')
# print(D)
# print('\n Matrice-extraite des lignes et colonnes impaires :')
# print(tamisD)

###########################################
## Exercice 4 : fonctions zeros, ones, eye...
###########################################

Nulle3 = np.zeros((3,3))
U3 = np.ones((3,3))
I4 = np.eye(4)

E = np.zeros((5,5))
for j in range(5) :
    E[0,j] = 1
    E[4,j] = 1
for i in range(5) :
    E[i,0] = 1
    E[i,4] = 1

D = 2*I4

L = [-1 for _ in range(5)]
F = np.eye(6) + np.diag(L,1) + np.diag(L,-1)

###########################################
## Exercice 5 : Produit matriciel
###########################################

# 1a Rappel mathématique
# formule du produit matriciel :
# Si A de taille (n,p) et B de taille (p,q),
# alors C = A x B existe et est de taille (n,q).
# Pour tout (i,j) avec i dans |[ 1,n ]| et j dans |[ 1,q ]|, on a
# C_ij = SOMME_{k dans |[ 1,p ]|} A_ik x B_kj


## 1b Fonction multipliant des matrices de tailles compatibles
def mult(A,B) :
    (n,p) = A.shape
    (r,q) = B.shape
    if p != r :
        return 'tailles incompatibles'
    C = np.zeros((n,q))
    for i in range(n) :
        for j in range(q) :
            S = 0
            for k in range(p) :
                S += A[i,k]*B[k,j]
            C[i,j] = S
    return C

## 2 comparaison mult et np.dot
A = np.array([[1,5,2],[4,-1,1]])
B = np.array([[2,-1],[5,4],[1,-2]])

print('Produit AB avec la fonction np.dot :')
print(np.dot(A,B))
print('\n Produit AB avec la fonction mult :')
print(mult(A,B))
print('\n Produit BA avec la fonction np.dot :')
print(np.dot(B,A))
print('\n Produit AB avec la fonction mult :')
print(mult(B,A))

## 3 puissance d'une matrice carrée
C = np.array([[1,5,2],[4,-1,1],[0,1,-1]])
print('C² avec np.dot :')
print(np.dot(C,C))
print('\n C² avec mult :')
print(mult(C,C))

## 4 matrices aléatoires
randA = np.array([[rd.randint(-5,5) for _ in range(3)] for _ in range(3)])
randB = np.array([[rd.randint(-5,5) for _ in range(3)] for _ in range(3)])

print('AB avec np.dot :')
print(np.dot(randA,randB))
print('\n AB avec mult :')
print(mult(randA,randB))

## 5 Fonction calculant A,B,C

def mult3(A,B,C) :
    (n,p) = A.shape
    (r,q) = B.shape
    (s,t) = C.shape
    if p != r or q != s :
        return 'tailles incompatibles'
    AB = mult(A,B)
    return mult(AB,C)

randA = np.array([[rd.randint(-5,5) for _ in range(2)] for _ in range(2)])
randB = np.array([[rd.randint(-5,5) for _ in range(3)] for _ in range(2)])
randC = np.array([[rd.randint(-5,5) for _ in range(2)] for _ in range(3)])

print('ABC calculé avec np.dot :')
print(np.dot(np.dot(randA,randB),randC))
print('\n ABC calculé avec mult3 :')
print(mult3(randA,randB,randC))

###########################################
## Exercice 6 : Puissances de matrices carrées
###########################################

# 1a : fonction puissance, version itérative (boucle for)

def puissance(A,n) :
    (p,q) = A.shape
    if p != q :
        return 'matrice non carée !!'
    I = np.eye(p)
    for _ in range(n) :
        I = mult(I,A)
    return I

# 1b : fonction puissance, version récursive

def puissanceRecursif(A,n) :
    (p,q) = A.shape
    if p != q :
        return 'matrice non carée !!'
    if n == 0 :
        return np.eye(p)
    return mult(A,puissanceRecursif(A,n-1))

## 2 : comparaison avec la.matrix_power

randA = np.array([[rd.randint(-5,5) for _ in range(2)] for _ in range(2)])
print('A**3 avec la fonction puissance :')
print(puissance(randA,3))
print('\n A**3 avec la.matrix_power :')
print(la.matrix_power(randA,3))

## 3 : vérification que pour tout n, I_3**n = I_3

for n in range(6) :
    print('\n I3**'+str(n)+' = ')
    print(puissance(np.eye(3),n))

## 4 : vérification pour les puissances d'une matrice diagonale

D = np.diag([2,3,4])

print("Puissance d'une matrice diagonale avec la fonction 'puissance'")
print(puissance(D,5))
print("\n Puissance d'une matrice diagonale avec la formule du cours")
print(np.diag([2**5,3**5,4**5]))

## 5 : vérification pour une matrice de rotation

B = np.array([[0,1],[-1,0]])

print("B**4 avec la fonction 'puissance'")
print(puissance(B,4))
print("\n B**4 avec 'matrix_power'")
print(la.matrix_power(B,4))


###########################################
## Exercice 7 : Rang, inverse et éléments propres d'une matrice carrée
###########################################

# 1 : calcul du rang d'une matrice

A = np.array([[3*i+j+1 for j in range(3)] for i in range(3)])
B = np.eye(3)

print('rg(A) = '+str(la.matrix_rank(A)))
print('rg(B) = '+str(la.matrix_rank(B)))

# 2 : calcul de l'inverse d'une matrice

C = np.array([[1,2,3],[0,-1,2],[4,0,-1]])
Cinv = la.inv(C)

print('\n inverse de C = ')
print(Cinv)

print('\n C*C^{-1} = ')
print(mult(C,Cinv))

# 3 rang, valeurs propres

E = np.array([[1,2,3],[0,1,4],[5,6,0]])

print('\n rg(E) = '+str(la.matrix_rank(E)))

spectre = ''
SP = la.eigvals(E)
for elt in SP :
    spectre = spectre + str(elt) + ', '
print('\n Sp(E) = {'+spectre+'}')

print('\n card(Sp(E)) = 3 donc E est diagonalisable.')

# 4 Cas des matrices nulles et de 1

print('\n La rang de la matrice nulle de taille 4 est : '+str(la.matrix_rank(np.zeros((4,4)))))

print('\n La rang de la matrice "un" de taille 4 est : '+str(la.matrix_rank(np.ones((4,4)))))