# spy03_13_18_ajustement_lineaire.py
import numpy as np
import matplotlib.pyplot as plt

print('------ 3.13 ---------------')
print("Droite passant par deux points donnés.")
#  b + -2*a = 1 
#  b +  3*a = 5 
# (1  -2) _ (1)
# (1   3) ¯ (5) 
matA = np.array([[1, -2], [1, 3]], np.float64)
vecb = np.array([1, 5], np.float64)
v_sol = np.linalg.lstsq(matA, vecb, rcond=None)[0] # Solution approchée = solution exacte ici
print("matA=\n", matA) 
print("v_sol=", v_sol)
print("matA*v_sol=", np.dot(matA, v_sol)) # Vérification.
b = v_sol[0]
a = v_sol[1]

x1 = matA[:,1] # Abscisses des points donnés
y1 = vecb # Ordonnées des points donnés
x2 = np.array([-3, 4])
y2 = a*x2 + b
plt.plot(x1, y1, 'o')
plt.plot(x2, y2, '-')
plt.title('droite passant par les deux points.')
plt.xlabel('x')
plt.ylabel('y = a*x + b')
plt.text(-3, 5.5, 'a = ' + '{:10.6f}'.format(a), color='black', fontsize=12)
plt.text(-3, 5.0, 'b = ' + '{:10.6f}'.format(b), color='black', fontsize=12)
plt.show()

print('------ 3.14 ---------------')
print("Droite passant au mieux par 4 points donnés.")
points = np.array([[-2, 1], [3, 5], [8, 9], [13, 12]], np.float64) # Données
matA = np.array([np.ones(points.shape[0]), points[:,0]]).T
vecb = points[:,1]
v_sol = np.linalg.lstsq(matA, vecb, rcond=None)[0] # Solution approchée
print("matA=\n", matA) 
print("Solution = ", v_sol)
print("matA*v_sol=", np.dot(matA, v_sol)) # Vérification.
b = v_sol[0]
a = v_sol[1]

xx = np.array([-3, 14]) # Abscisses des extrémités de la droite
yy = a*xx + b # Ordonnées des extrémités de la droite
plt.plot(points[:,0], points[:,1], 'o') # Affiche les points
plt.plot(xx, yy, '-') # Trace la droite
plt.title('droite passant "au mieux" par les points.')
plt.xlabel('x')
plt.ylabel('y = a*x + b')
plt.text(-2.5, 12, 'a = ' + '{:10.6f}'.format(a), color='black', fontsize=12)
plt.text(-2.5, 11, 'b = ' + '{:10.6f}'.format(b), color='black', fontsize=12)
plt.axis([-3, 14, 0, 14]) # Limites des axes
plt.show()

print('------ 3.15 ---------------')
print("Droite passant au mieux par 5 points donnés.")
points = np.array([[-2, 13], [0, 10], [2, 8], [4, 7], [6, 5]], np.float64) # Données
matA = np.array([np.ones(points.shape[0]), points[:,0]]).T
print("A =\n", matA)

vecb = points[:,1]
v_sol = np.linalg.lstsq(matA, vecb, rcond=None)[0] # Solution approchée
print("Solution = ", v_sol)
b = v_sol[0]
a = v_sol[1]

xx = np.array([-3, 7]) # Abscisses des extrémités de la droite
yy = a*xx + b # Ordonnées des extrémités de la droite
plt.plot(points[:,0], points[:,1], 'o') # Affiche les points
plt.plot(xx, yy, '-') # Trace la droite
plt.title('droite passant "au mieux" par les points.')
plt.xlabel('x')
plt.ylabel('y = a*x + b')
plt.text(-2.5, 2, 'a = ' + '{:10.6f}'.format(a), color='black', fontsize=12)
plt.text(-2.5, 1, 'b = ' + '{:10.6f}'.format(b), color='black', fontsize=12)
plt.axis([-3, 7, 0, 14]) # Limites des axes
plt.show()

print('------ 3.16 ---------------')
print("Détermination d'un polynôme d'interpolation.")
ax = np.array([-2, -1, 0, 1, 2], np.float64) # Donnée
ay = np.sin(ax) # Donnée

# Création de la matrice
matA = np.zeros((ax.shape[0], ax.shape[0]), np.float64)
matA[:, 0] = np.ones(ax.shape[0]) # Première colonne de la matrice

# Remplissage des autres colonnes de la matrice
for jj in range(1, ax.shape[0]):
    matA[:, jj] = ax ** jj

print("A =\n", matA)

# Détermination des coefficients du polynôme
coefs = np.linalg.lstsq(matA, ay, rcond=None)[0] # Coefficients du polynôme
print("Solution = ", coefs)

# Graphique, les points, le sinus et le polynôme.
xx = np.linspace(-2.5, 2.5, 21,endpoint=True) # Abscisses des points pour tracer le polynôme
yy = coefs[0] + coefs[1]*xx + coefs[2]*xx**2 + coefs[3]*xx**3 + coefs[4]*xx**4
plt.plot(ax, ay, 'o') # Affiche les points
plt.plot(xx, yy, '-') # Trace le polynôme
plt.plot(xx, np.sin(xx), '-') # Trace le sinus
plt.title('Polynôme passant par les points.')
plt.xlabel('x')
plt.ylabel('y = sin(x)')
plt.legend(['points', 'y=poly(x)', 'y=sin(x)'], loc='upper left')
plt.axis([-2.5, 2.5, -1.5, 1.5]) # Limites des axes
plt.show()

print('------ 3.16 bis ---------------')
ax = np.array([-2, -1, 0, 1, 2], np.float64) # Donnée
ay = np.sin(ax) # Donnée

coefs = np.polyfit(ax, ay, 4) # Coefficients du polynôme de degré 4
print("Solution = ", coefs)
coefs = np.flip(coefs) # Change l'ordre des nombres dans le vecteur
print("Solution = ", coefs)

# Graphique, les points, le polynôme et le sinus.
xx = np.linspace(-2.5, 2.5, 81, endpoint=True) # Abscisses des points pour tracer le polynôme
yy = coefs[0] + coefs[1]*xx + coefs[2]*xx**2 + coefs[3]*xx**3 + coefs[4]*xx**4
plt.plot(ax, ay, 'o') # Affiche les points
plt.plot(xx, yy, '-') # Trace le polynôme
plt.plot(xx, np.sin(xx), '-') # Trace le sinus
plt.title('Polynôme passant par les points.')
plt.xlabel('x')
plt.ylabel('y = sin(x)')
plt.legend(['points', 'y=poly(x)', 'y=sin(x)'], loc='upper left')
plt.axis([-2.5, 2.5, -1.5, 1.5]) # Limites des axes
plt.show()

print('------ 3.17 ---------------')
print("Polynôme passant au mieux par des points donnés.")
ax = np.linspace(-2.5, 2.5, 21, np.float64) # Donnée
ay = np.sin(ax) # Donnée
nDeg = 17
coefs = np.polyfit(ax, ay, nDeg) # Coefficients du polynôme de degré 3
coefs = np.flip(coefs) # Change l'ordre des nombres dans le vecteur
print("Solution = ", coefs)

# Graphique, les points, le polynôme et le sinus.
xx = np.linspace(-2.5, 2.5, 81, np.float64) # Abscisses des points pour tracer le polynôme
yy = coefs[0]
for jj in range(1,nDeg+1): yy += coefs[jj]*xx**jj
plt.plot(ax, ay, 'o') # Affiche les points
plt.plot(xx, yy, '-') # Trace le polynôme
plt.plot(xx, np.sin(xx), '-') # Trace le sinus
plt.title('Polynôme passant au mieux par les points.')
plt.xlabel('x')
plt.ylabel('y = sin(x)')
plt.legend(['points', 'y=poly(x)', 'y=sin(x)'], loc='upper left')
plt.axis([-2.5, 2.5, -1.5, 1.5]) # Limites des axes
plt.show()

ydiff = yy - np.sin(xx)
plt.plot(xx, ydiff, '-') # Trace la différence entrel le polynôme et le sinus
plt.title('Différence : "Polynôme passant au mieux par les points" - Sinus.')
plt.xlabel('x')
plt.ylabel('y = différence')
plt.xlim(-2.5, 2.5) # Limites de l'axe X c.f. https://python.developpez.com/tutoriels/graphique-2d/matplotlib/#LII-D
fig = plt.gcf(); fig.set_size_inches(8, 5) # autre manière de faire
plt.show()


print('------ 3.18 ---------------')
print("Approximation d'une fonction paire par une somme de cosinus.")
ax = np.linspace(-np.pi, np.pi, 35, np.float64) # Donnée
ay = np.exp(-5 * ax**2) # Donnée
nbCoef = 15

# Création de la matrice
matA = np.zeros((ax.shape[0], nbCoef), np.float64)
matA[:, 0] = np.ones(ax.shape[0]) # Première colonne de la matrice

# Remplissage des autres colonnes de la matrice
for jj in range(1, nbCoef):
    matA[:, jj] = np.cos(jj * ax)
print("Taille de la matrice =", matA.shape)
#print("A=\n", matA)  # Pour info.

# Détermination des coefficients
coefs = np.linalg.lstsq(matA, ay, rcond=None)[0] # Coefficients
print("Solution = ", coefs)

# Graphique, les points, la fonction et son approximation.
xx = np.linspace(-np.pi, np.pi, 81,endpoint=True) # Abscisses des points pour tracer les fonctions
yy = coefs[0]
for jj in range(1, nbCoef): yy += coefs[jj] * np.cos(jj*xx)
plt.plot(ax, ay, 'o') # Affiche les points
plt.plot(xx, yy, '-') # Trace l'approximation
plt.plot(xx, np.exp(-5 * xx**2), '-') # Trace la courbe à approximer
plt.title('Approximation de : y = exp(-5 * x**2).')
plt.xlabel('x [rad]')
plt.ylabel('y = approx. et courbe')
plt.legend(['points', 'y=approx(x)', 'y=exp(-5 * x**2)'], loc='upper left')
plt.axis([-3.5, 3.5, -0.1, 1.1]) # Limites des axes
plt.show()

ydiff = yy - np.exp(-5 * xx**2)
plt.plot(xx, ydiff, '-') # Trace la différence entrel la fonction et son approximation
plt.title('Approximation - fonction')
plt.xlabel('x [rad]')
plt.ylabel('y = différence')
plt.xlim(-3.5, 3.5) # Limites de l'axe X c.f. https://python.developpez.com/tutoriels/graphique-2d/matplotlib/#LII-D
fig = plt.gcf(); fig.set_size_inches(8, 5) # autre manière de faire
plt.show()