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
from numpy.polynomial.polynomial import polyval
from numpy import linalg as LA
from random import uniform

##

def coefficient_polynome_Lagrange(vecteurx, vecteurf):
    
    # Verification des dimensions du vecteur vecteurx et du vecteur vecteurf
    if vecteurx.shape[0]!=vecteurf.shape[0] or vecteurx.shape[1]!=vecteurf.shape[1]:
        print("Probleme : taille")
        return 
    
    n=len(vecteurx)
    A=zeros((n,n))
    
    for i in range(0,n):
        for j in range(0,n):
            A[i,j]=(vecteurx[i,0])**j
    # Résolution de systèmes linéaires A x = vecteurf
    x = np.linalg.solve(A, vecteurf)
    
    ConditionnementA = cond(A)
    
    return x , ConditionnementA
##

def representation_polynome_interpolation():
    
    # Abscisse des points a interpoler
    vecteurx = linspace(0, 3, 5)
    vecteurx=vecteurx.reshape( len(vecteurx) ,1) 
    # Ordonnees des points a interpoler (definies de maniere aleatoire)
    vecteurf = array([uniform(0,1) for i in range( len(vecteurx) )])
    vecteurf=vecteurf.reshape( len(vecteurf) ,1) 
    
    # Calcul des coefficients du polynome d interpolation de vecteurx, vecteurf
    vecteura , conditionnementA = coefficient_polynome_Lagrange(vecteurx, vecteurf)
    
    # Trace du polynome d interpolation
    # Recherche d un intervalle [a,b] contenant tous les points a interpoler
    L = max(vecteurx) - min(vecteurx) 
    a = min(vecteurx) - L/10 
    b = max(vecteurx) + L/10
   
    # Vecteur des points auquel on va calculer p(t) 
    vecteurt = linspace(0,4, 100)
    vecteurp = np.polyval(vecteura[range(len(vecteurx)-1,-1,-1)] , vecteurt)
    
    #TRACE
    fig = figure(figsize=(14,8))
    plt.plot(vecteurx, vecteurf,"o",label="Points d'interpolation")
    plt.plot(vecteurt, vecteurp,"--",label="Polynome d'interpolation")
    plt.grid(True)
    plt.legend()
    show()
     
##
def test_coefficient_polynome_Lagrange():
    
    vecteurConditionnement=zeros((20,1))
    for i in range(5,25):
        
        # Abscisse des points a interpoler
        vecteurx = linspace(0, 3, i)
        vecteurx=vecteurx.reshape( len(vecteurx) ,1) 
        
        # Ordonnees des points a interpoler (definies de maniere aleatoire)
        vecteurf = array([uniform(0,1) for i in range( len(vecteurx) )])
        vecteurf=vecteurf.reshape( len(vecteurf) ,1) 
        
        # Calcul des coefficients du polynome d interpolation de vecteurx, vecteurf
        vecteura , conditionnementA = coefficient_polynome_Lagrange(vecteurx, vecteurf)
        
        vecteurConditionnement[i-5,0]=conditionnementA
    
    #TRACE
    fig = figure(figsize=(14,8))
    n=np.arange(5,25)
    plt.semilogy(n,vecteurConditionnement)
    plt.grid(True)
    plt.xlabel("n")
    plt.ylabel("Conditionnement de A")
    plt.legend()
    show()

##    

#La fonction Lagrange permet de calculer le polynôme d'interpolation de Lagrange

def Lagrange(vecteurx,vecteury,t):
    
    # Verification des dimensions du vecteur vecteurx et du vecteur vecteurf
    #if vecteurx.shape[0]!=vecteurf.shape[0] or vecteurx.shape[1]!=vecteurf.shape[1]:
    #print("Probleme : taille")
    #    return 
    
    n=len(vecteurx)
    p=0
    for i in range(0,n):
        L=1
        for j in range(0,n):
            if i!=j :
                L= L * ( (t-vecteurx[j,0]) / (vecteurx[i,0]-vecteurx[j,0]) )
        p= p + vecteury[i,0] * L 
    
    return p
##

# f est une fonction continue sur [a,b]
# n nombre de points a interpoler
# a et b deux réels tq a<b  

def InterpolationLagrange_Pointsuniformes(f,n,a,b):
    
    x=np.linspace(a,b,100)
    
    #Points uniformes
    xi=np.linspace(a,b,n)
    xi=xi.reshape( len(xi) ,1) #redimensionnement 
    
    #Ordonnnees des points a interpoler 
    fi=f(xi)
    fi=fi.reshape( len(fi) ,1)#redimensionnement  
    
    #Calcul de pi(x)
    p=zeros((len(x),1))
    for i in range(0,len(x)):
        p[i,0]=Lagrange(xi,fi,x[i])
    
    #Trace
    fig = figure(figsize=(14,8))
    plt.plot(x,p,label="La fonction y=Pn(t) ")
    plt.plot(x,f(x),label="La fonction y=f(t) ")
    plt.plot(xi, fi,"o",label="Points uniformes")
    plt.grid(True)
    title("Polynomes d'interpolation de Lagrange")
    plt.legend()
    show()



##

# f est une fonction continue sur [a,b]
# n nombre de points a interpoler
# a et b deux réels tq a<b  

def InterpolationLagrange_PointsTchebyshev(f,n,a,b):
# f est une fonction continue sur [a,b]
# n nombre de points a interpoler
# a et b deux réels tq a<b      
    x=np.linspace(a,b,100)
    
    #Points de Tchebyshev
    xi=array([0.5*(a+b) +0.5*(b-a)*cos( i*pi/n ) for i in range(0,n+1)])    
    xi=xi.reshape( len(xi) ,1) #redimensionnement

    #Ordonnnees des points a interpoler 
    fi=f(xi)
    fi=fi.reshape( len(fi) ,1) #redimensionnement 
    
    #Calcul de pi(x)
    p=zeros((len(x),1))
    for i in range(0,len(x)):
        p[i,0]=Lagrange(xi,fi,x[i])
    
    #Trace
    fig = figure(figsize=(14,8))
    plt.plot(x,p,label="La fonction y=Pn(t) ")
    plt.plot(x,f(x),label="La fonction y=f(t) ")
    plt.plot(xi, fi,"o",label="Points Tchebyshev")
    plt.grid(True)
    title("Polynomes d'interpolation de Lagrange")
    plt.legend()
    show()
##

def PolyLi(i,x,t):
    
    n=len(x)
    Li=1
    for j in range(0,n):
        if i != j :
            Li= Li * ( (t-x[j,0])/(x[i,0]-x[j,0]) )
        
    return Li
    
##

def PolyLp(i,x):
    
    n=len(x)
    Lp=0
    for j in range(0,n):
        if i != j :
            Lp= Lp + ( 1/(x[i,0]-x[j,0]) )
            
    return Lp
    
##

def PolyA(i,x,t):
    
    A= ( 1- 2*PolyLp(i,x) * (t-x[i,0]) ) * ( PolyLi(i,x,t)**2 )
    
    return A
    
##

def PolyB(i,x,t):
    
    B= ( t-x[i,0] ) * ( PolyLi(i,x,t)**2 )
    
    return B
    
##

def Hermite(x,y,z,t):
    n=len(x)
    H=0
    for i in range(0,n):
        H= H + y[i,0]*PolyA(i,x,t) + z[i,0]*PolyB(i,x,t)
        
    return H

##

# f est une fonction de classe C1 sur [a,b]
# fp est la dérivée de f
# n nombre de points a interpoler
# a et b deux réels tq a<b  
 
def InterpolationHermite_PointsChebyshev(f,fp,n,a,b):
        
    x=np.linspace(a,b,100)
    
    #Points de Tchebyshev

    xi=array([0.5*(a+b) +0.5*(b-a)*cos( i*pi/n ) for i in range(0,n+1)])    
    xi=xi.reshape( len(xi) ,1) #redimensionnement

    #Ordonnnees des points a interpoler 
    yi=f(xi)
    yi=yi.reshape( len(yi) ,1) #redimensionnement
    
    zi=fp(xi)
    zi=zi.reshape( len(zi) ,1) #redimensionnement
    
    #Calcul de Hi(x)
    H=zeros((len(x),1))
    for i in range(0,len(x)):
        H[i,0]=Hermite(xi,yi,zi,x[i])
    
    #Trace
    fig = figure(figsize=(14,8))
    plt.plot(x,H,label="La fonction y=Hn(t) ")
    plt.plot(x,f(x),label="La fonction y=f(t) ")
    plt.plot(xi, yi,"o",label="Points Tchebyshev")
    ylim(-1,1)
    plt.grid(True)
    title("Polynomes d'interpolation de Lagrange-Hermite")
    plt.legend()
    show()
 
## 

# f est une fonction de classe C1 sur [a,b]
# fp est la dérivée de f
# n nombre de points a interpoler
# a et b deux réels tq a<b  
 
def InterpolationHermite_Pointsuniformes(f,fp,n,a,b):
        
    x=np.linspace(a,b,100)
    
    #Points uniformes
    xi=np.linspace(a,b,n)
    xi=xi.reshape( len(xi) ,1) #redimensionnement 

    #Ordonnnees des points a interpoler 
    yi=f(xi)
    yi=yi.reshape( len(yi) ,1) #redimensionnement
    
    zi=fp(xi)
    zi=zi.reshape( len(zi) ,1) #redimensionnement
    
    #Calcul de Hi(x)
    H=zeros((len(x),1))
    for i in range(0,len(x)):
        H[i,0]=Hermite(xi,yi,zi,x[i])
    
    #Trace
    fig = figure(figsize=(14,8))
    plt.plot(x,H,label="La fonction y=Hn(t) ")
    plt.plot(x,f(x),label="La fonction y=f(t) ")
    plt.plot(xi, yi,"o",label="Points uniformes")
    plt.grid(True)
    title("Polynomes d'interpolation de Lagrange-Hermite")
    ylim(-1,1)
    plt.legend()
    show() 
 
##