import numpy as np
import matplotlib.pyplot as mp 
from math import *
import scipy.optimize as sp
import matplotlib.pyplot as plt
from pylab import *
from scipy import *
import scipy as sp 
import scipy.integrate as integ
from scipy.integrate import odeint
import time 
import timeit 
import random 
import datetime
##
def ProduitMatrice(A,B):
    
    if ( (A.shape[1]==B.shape[0]) and (B.shape[1]!=1) ) :
        P=zeros((len(A),len(A)))
        S=0
        for i in range(0,len(A)):
            for j in range(0,len(A)):
                for k in range(0,len(A)):
                    S=S+A[i,k]*B[k,j]
                P[i,j]=S
                S=0
    elif ( (A.shape[1]==B.shape[0]) and (B.shape[1]==1) ):
        P=zeros((len(A),1))
        S=0
        for i in range(0,len(A)):
            for j in range(0,len(A)):
                for k in range(0,len(A)):
                    S=S+A[i,k]*B[k,0]
                P[i,0]=S
                S=0
    else:
        print("Probleme : taille")
        return
    
    return P
##

# La fonction triangulaire_inférieure vérifie si une matrice carrée de taille NXN est triangulaire inférieure ou non 
# Cette fonction renvoie True si elle est triangulaire inférieure ou False sinon
# A est triangulaire inférieure si et seulement si A[i,j]=0 pour j>i
def triangulaire_inferieure(A):
    n=len(A)
    j=0
    test=True #on suppose que A est triangulaire inférieure
    while test==True and j<=n-1 :
        for i in range(0,j):
            if A[i,j]!=0:
                test=False
        j=j+1
    return test
  
##

# La fonction triangulaire_supérieure vérifie si une matrice carrée de taille NXN est triangulaire supérieure ou non 
# Cette fonction renvoie True si elle est triangulaire supérieure ou False sinon
# A est triangulaire supérieure si et seulement si A[i,j]=0 pour i>j

def triangulaire_superieure(A):
    n=len(A)
    i=0
    test=True #on suppose que A est triangulaire supérieure
    while test==True and i<=n-1 :
        for j in range(0,i):
            if A[i,j]!=0:
                test=False
        i=i+1
    return test
##
# La fonction diagonale vérifie si une matrice carrée de taille NXN est diagonale ou non 
# Cette fonction renvoie True si elle est diagonale ou False sinon
def diagonale(A):
    triang_inf=triangulaire_inferieure(A)
    triang_sup=triangulaire_superieure(A)
    if triang_inf==True and triang_sup==True :
        return True
    else :
        return False

##

# La fonction diago_inver verifie si la matrice A supposée diagonale de taille N X N est inversible ou non 
# A est inversible si et seulement si pour tout i dans [[1,N]] a[i,i] != 0

def diago_inver(A):
    verif=diagonale(A)
    i=0
    n=A.shape[0]
    test=True
    if verif==True :
        while  test==True and i<n : 
            if A[i,i]==0 :
                test=False
            else :
                i=i+1
    else :
        test=False            
        
    return test
##


#                         "METHODES DIRECTES" 

##
 
# La fonction RslMatDiag permet de résoudre le système linéaire à matrice diagonale inversible 

def RslMatDiag(A,b):
    verif=diago_inver(A)
    if verif==True : 
        n=A.shape[0]
        x=zeros([n,1]) 
        for i in range(0,n):
            x[i,0]=b[i,0]/A[i,i]
        return x
    else :
        print(" A n'est pas inversible !!! ")
            
##        
#La fonction RslTriInf permet de résoudre le système linéaire triangulaire inférieure inversible
#Resolution de A x = b   
#A est une matrice traingulaire inférieure (carree de taille N)
#b est un vecteur colonne (de taille N)

def RslTriInf(A,b):
    
    #Verification des dimensions de la matrice U et du vecteur b
    if A.shape[1]!=b.shape[0] or b.shape[1]!=1:
        print("Probleme : taille")
        return 
        
    if A.shape[0]!=A.shape[1]:
        print("Probleme : la matrice doit être carrée")
        return
        
    #On vérifie si la matrice A est triangulaire inférieure
    verif=triangulaire_inferieure(A)
    
    #Verification de l'inversibilite  de A 
    #A est inversible si A[i,i] != 0
    precision = 10**(-13) 
    i=0
    test=True
    while i<len(A) and test==True :
        if abs(A[i,i])<=precision:
            print("Attention : la matrice A n'est peut etre pas inversible")
            test=False
            return
        else :
            i=i+1
    #Calcul de x       
    if verif==True :
        n=A.shape[0]
        x=zeros([n,1])
        x[0,0]=b[0,0]/A[0,0]
        for i in range(1,n):
            s=0
            for j in range(0,i):
                s=s+A[i,j]*x[j,0]
            x[i,0]=(b[i,0]-s)/A[i,i]
        return x
    else :
        print(" A n'est pas triangulaire inférieure !!! ")
    

##
#La fonction RslTriSup permet de résoudre le système linéaire triangulaire supérieure inversible
#Resolution de A x = b   
#A est une matrice traingulaire superieure (carree de taille N)
#b est un vecteur colonne (de taille N)

def RslTriSup(A,b):
    #Verification des dimensions de la matrice U et du vecteur b
    if A.shape[1]!=b.shape[0] or b.shape[1]!=1:
        print("Probleme : taille")
        return 
        
    if A.shape[0]!=A.shape[1]:
        print("Probleme : la matrice doit être carrée")
        return
        
    #On vérifie si la matrice A est triangulaire inférieure
    verif=triangulaire_superieure(A)
    
    #Verification de l'inversibilite  de A 
    #A est inversible si A[i,i] != 0
    precision = 10**(-13) 
    i=0
    test=True
    while i<len(A) and test==True :
        if abs(A[i,i])<=precision:
            print("Attention : la matrice A n'est peut etre pas inversible")
            test=False
            return
        else :
            i=i+1
    #Calcul de x       
    if verif==True :
        n=len(A)
        x=zeros((n,1))
        for i in range(n-1,-1,-1):
            s=0
            for j in range(i,n):
                s=s+A[i,j]*x[j,0]
            x[i,0]=(b[i,0]-s)/A[i,i]
        return x
    else :
        print(" A n'est pas triangulaire supérieure !!! ")
    
##

#La fonction FactLU permet de calculer les matrices L et U dites matrice de factorisation LU associée
#à la matrice A telle que A = L U
"A est une matrice de Mn(K) dont les sous-matrices principales sont inversibles."
"L est une matrice de Mn(K) triangulaire inférieure avec L[i,i]=1"
"U est une matrice de Mn(K) triangulaire supérieure."

def FactLU(A):
    
    #Vérification si A est une matrice carrée ou non
    if A.shape[0]!=A.shape[1]:
        print("Probleme : la matrice doit être carrée")
        return
        
    n=len(A)
    U=zeros((n,n))
    L=eye(n,n)
    for i in range(0,n):
        for j in range(i,n):
            S1=0
            for k in range(0,i):
                S1=S1+L[i,k]*U[k,j]
            U[i,j]=A[i,j]-S1
           
        for j in range(i,n):
            S2=0
            for k in range(0,i):
                S2=S2+L[j,k]*U[k,i]
            L[j,i]=(A[j,i]-S2)/U[i,i]
        
    return L,U
        
##

#La fonction RSLFactLU permet de résoudre, par une factorisation LU le système linéaire A x = b
"A est une matrice de Mn(K) dont les sous-matrices principales sont inversibles définie positive."
"b : un vecteur de Rn"
"x : un vecteur de Rn"

def RSLFactLU(A,b):
    #Verification des dimensions de la matrice A et du vecteur b
    if A.shape[1]!=b.shape[0] or b.shape[1]!=1:
        print("Probleme : taille")
        return 
        
    if A.shape[0]!=A.shape[1]:
        print("Probleme : la matrice doit être carrée")
        return
        
    L,U=FactLU(A)
    y=RslTriInf(L,b)
    x=RslTriSup(U,y)
    return x
    
##  

#La fonction Cholesky permet de calculer la matrice B dites matrice de factorisation 
#positive de Cholesky associée à la matrice A telle que A = BB*
"A est une matrice de Mn hermitienne définie positive"

def Cholesky(A):
    
    #Vérification si A est une matrice carrée ou non
    if A.shape[0]!=A.shape[1]:
        print("Probleme : la matrice doit être carrée")
        return
    
    #Vérification si A est une matrice hermitienne ou non
    if Hermitienne(A)==False :
        print("Problème : la matrice n'est pas hermitienne")
        return
    #Vérification si A est une matrice définie positive ou non
    if DefiniePositive(A)==False :
        print("Problème : la matrice n'est pas définie positive")
        return
    
    #Si A est une matrice de Mn hermitienne définie positive
            
    n=len(A)
    B=zeros((n,n))
    
    if ( (DefiniePositive(A)==True) and (Hermitienne(A)==True) ) :
        for i in range(0,n):
            S1=0
            for j in range(0,i+1):
                S1=S1+abs(B[i,j]**2)
            B[i,i]=sqrt(A[i,i]-S1)
            
            for j in range(i,n):
                S2=0
                for k in range(0,i):
                    S2=S2+B[j,k]*conj(B[i,k])
                B[j,i]=(A[j,i]-S2)/B[i,i]
             
    return B
##  

#La fonction RSLCholesky permet de résoudre, par une factorisation de Cholesky positive, le système linéaire  A x = b
#ou A est une matrice de Mn hermitienne définie positive

def  RSLCholesky(A,b):
    
    #Verification des dimensions de la matrice A et du vecteur b
    if A.shape[1]!=b.shape[0] or b.shape[1]!=1:
        print("Probleme : taille")
        return 
        
    if A.shape[0]!=A.shape[1]:
        print("Probleme : la matrice doit être carrée")
        return
        
    #Vérification si A est une matrice hermitienne ou non
    if Hermitienne(A)==False :
        print("Problème : la matrice n'est pas hermitienne")
        return
        
    #Vérification si A est une matrice définie positive ou non
    if DefiniePositive(A)==False :
        print("Problème : la matrice n'est pas définie positive")
        return
        
    B=Cholesky(A)
    y=RslTriInf(B,b)
    U=MatriceAdjointe(B)
    x=RslTriSup(U,y)
    
    return x
    
##


#                           METHODES ITERATIVES

##   

# Resolution du systeme lineaire A x = b par la methode de Jacobi

" Parametres d'entree"
# A = matrice carree inversible de taille nXn 
# b = vecteur de taille nX1
# x0 = vecteur de taille nX1 (pour l initialisation)
# epsilon = precision souhaitee sur le residu
# kmax = nombre maximal d'iterations

"Parametres de sortie"
# x : vecteur de taille n
# n : nombre d iteration necessaire pour avoir un precision de epsilon sur le residu
# vecteur_residu : vecteur de taille le nombre d iterations qui contient l'evolution de la norme du residu


# Resolution du systeme lineaire A x = b par la methode de Jacobi

"Remarque"
"si A est à diagonale strictement dominante alors la méthode de Jacobi converge"

def RSLJacobi(A,b,x0,epsilon,kmax):
    
    #Verification des dimensions de la matrice A et du vecteur b
    if A.shape[1]!=b.shape[0] or b.shape[1]!=1:
        print("Probleme : taille")
        return 
        
    if A.shape[0]!=A.shape[1]:
        print("Probleme : la matrice doit être carrée")
        return
    #Verification de A[i,i] != 0
    i=0
    precision = 10**(-16) 
    while i<len(A):
        if abs(A[i,i])<=precision:
            print("Attention : la matrice A n'est peut etre pas inversible")
            return
        else :
            i=i+1
    #
    k=1
    x=x0
    Xp1=zeros((len(x0),1))# x(k+1)
    residu=linalg.norm( dot(A,x) - b )
    vecteur_residu=[]
    vecteur_residu.append(residu)
    
    while k<=kmax and residu>=epsilon:
    
        for i in range(0,len(A)):
            S=0
            for j in range(0,len(A)):
                if i != j:
                    S=S + A[i,j]*x[j,0]
            Xp1[i,0]= (b[i,0]-S)/A[i,i]
        
        x=Xp1 
        residu=linalg.norm( dot(A,x) - b )
        vecteur_residu.append(residu)
        k=k+1
    #
    
    if residu<=epsilon :
        err= abs( linalg.norm( np.linalg.solve(A, b) - x) )
        print("Erreur : ",err)
        return x ,vecteur_residu,k
    
    elif k==kmax or residu>epsilon :
        print("Probleme de convergence : ")
        print("Le nombre maximal d iterations est atteint ou le residu est vraiment tres grand")
        print("Il se peut que l'algorithme diverge")
        print("Norme du residu : ",residu, "nombre d'itération :",k)

##

# Resolution du systeme lineaire A x = b par la methode de Gauss-Seidel

" Parametres d'entree"
# A = matrice carree inversible de taille nXn 
# b = vecteur de taille nX1
# x0 = vecteur de taille nX1 (pour l initialisation)
# epsilon = precision souhaitee sur le residu
# kmax = nombre maximal d'iterations

"Parametres de sortie"
# x : vecteur de taille n
# n : nombre d iteration necessaire pour avoir un precision de epsilon sur le residu
# vecteur_residu : vecteur de taille le nombre d iterations qui contient l'evolution de la norme du residu


# Resolution du systeme lineaire A x = b par la methode de Gauss-Seidel

"Remarque"
"si A est à diagonale strictement dominante alors la méthode de Gauss-Seidel converge"

def RSLGaussSeidel(A,b,x0,epsilon,kmax):
    
    #Verification des dimensions de la matrice A et du vecteur b
    if A.shape[1]!=b.shape[0] or b.shape[1]!=1:
        print("Probleme : taille")
        return 
        
    if A.shape[0]!=A.shape[1]:
        print("Probleme : la matrice doit être carrée")
        return
    #Verification de A[i,i] != 0
    i=0
    precision = 10**(-16) 
    while i<len(A):
        if abs(A[i,i])<=precision:
            print("Attention : la matrice A n'est peut etre pas inversible")
            return
        else :
            i=i+1
    #
    k=1
    x=x0
    Xp1=zeros((len(x0),1))# x(k+1)
    residu=linalg.norm( dot(A,x) - b )
    vecteur_residu=[]
    vecteur_residu.append(residu)
    
    while k<=kmax and residu>=epsilon:
        p=x
        for i in range(0,len(A)):
            S=0
            for j in range(0,i):
                S= S + A[i,j]*x[j,0]
            for j in range(i+1,len(A)):
                S= S + A[i,j]*p[j,0]
            
            Xp1[i,0]=(b[i,0]-S)/A[i,i]
            x=Xp1
            
        residu=linalg.norm( dot(A,x) - b )
        vecteur_residu.append(residu)
        k=k+1
    
    #
    
    if residu<=epsilon :
        err= abs( linalg.norm( np.linalg.solve(A, b) - x) )
        print("Erreur : ",err)
        return x ,vecteur_residu,k
    
    elif k==kmax or residu>epsilon :
        print("Probleme de convergence : ")
        print("Le nombre maximal d iterations est atteint ou le residu est vraiment tres grand")
        print("Il se peut que l'algorithme diverge")
        print("Norme du residu : ",residu, "nombre d'itération :",k)       

##  

# Resolution du systeme lineaire A x = b par la methode de relaxation

" Parametres d'entree"
# A = matrice carree inversible de taille nXn 
# b = vecteur de taille nX1
# x0 = vecteur de taille nX1 (pour l initialisation)
# epsilon = precision souhaitee sur le residu
# kmax = nombre maximal d'iterations
# w = réel non nul

"Parametres de sortie"
# x : vecteur de taille n
# n : nombre d iteration necessaire pour avoir un precision de epsilon sur le residu
# vecteur_residu : vecteur de taille le nombre d iterations qui contient l'evolution de la norme du residu


# Resolution du systeme lineaire A x = b par la methode de relaxation

"Remarque"
" si  0 < w < 2 alors la méthode de relaxation converge"

def RSLRelaxation(A,b,x0,epsilon,kmax,w):
    
    #Verification des dimensions de la matrice A et du vecteur b
    if A.shape[1]!=b.shape[0] or b.shape[1]!=1:
        print("Probleme : taille")
        return 
        
    if A.shape[0]!=A.shape[1]:
        print("Probleme : la matrice doit être carrée")
        return
    #Verification de A[i,i] != 0
    i=0
    precision = 10**(-16) 
    while i<len(A):
        if abs(A[i,i])<=precision:
            print("Attention : la matrice A n'est peut etre pas inversible")
            return
        else :
            i=i+1
    #
    k=1
    x=x0
    Xp1=zeros((len(x0),1))# x(k+1)
    residu=linalg.norm( dot(A,x) - b )
    vecteur_residu=[]
    vecteur_residu.append(residu)
    
    while k<=kmax and residu>=epsilon:
        p=x
        for i in range(0,len(A)):
            S=0
            for j in range(0,i):
                S= S + A[i,j]*x[j,0]
            for j in range(i+1,len(A)):
                S= S + A[i,j]*p[j,0]
            
            Xp1[i,0]=( w*(b[i,0]-S) )/A[i,i] +(1-w)*p[i,0]
            x=Xp1
            
        residu=linalg.norm( dot(A,x) - b )
        vecteur_residu.append(residu)
        k=k+1
    
    #
    
    if residu<=epsilon :
        err= abs( linalg.norm( np.linalg.solve(A, b) - x) )
        print("Erreur : ",err)
        return x ,vecteur_residu,k
    
    elif k==kmax or residu>epsilon :
        print("Probleme de convergence : ")
        print("Le nombre maximal d iterations est atteint ou le residu est vraiment tres grand")
        print("Il se peut que l'algorithme diverge")
        print("Norme du residu : ",residu, "nombre d'itération :",k)   


##  
def tic():
    debut=time.clock()
    return debut
##
def toc(debut):
    fin=debut-time.clock() 
    return fin            

##        
def Norme(V,p):
    S=0
    for i in range(0,len(V)):
        S=S+( abs(V[i,0])**p )
    return S**(1/p)
##
def MatriceAdjointe(A):
    M=conj(A.T)
    return M
##
def DefiniePositive(x):
    return np.all(np.linalg.eigvals(x) > 0)
##
def Hermitienne(A):
    B=conj(A.T)
    T=find((B==A)==False)
    if len(T)==0:
        return True
    else:
        return False
##