# Les matrices peuvent être représentées par des objets de type numpy.array. 
# Créér un array peut se faire à partir d'une liste de liste:
A = np.array([[10,6,6],
              [8,2,4],
              [8,4,2]])

# Accéder à $a_{ij}$ se fait par indexation:
A[0,1] # => 6

# Il existe des constructeurs pour des matrices particulières:
I = np.identity(3) # Matrice identité de taille 3
N = np.zeros((3,3)) #Matrice de zeros de taille 3
X,Y = np.meshgrid([1,2,3],[1,2,3]) # Matrices donnant le produit cartésien de deux vecteurs. 

x = np.array([1,2,3])
# Le produit matriciel pour les np.array peut se faire avec la méthode .dot:
A.dot(x) # Seule possibilité pour Python < 3.5
# Ou l'opérateur @:
A @ x #Python >= 3.5
assert np.all(A.dot(x) ==  A @ x)

### Vecteurs lignes, vecteurs colonnes... ###
# Le produit matriciel n'est défini que pour des matrice de dimensions compatibles. 
# Exemple, le produit: A * B necessite que B ait autant de colonnes que A ait de lignes. 
# Le résultat est une matrice carrée. 
# (M,N) * (K,M) -> (M,M)

# Par défaut les array numpy de dimension 1 sont orientés à chaque opération pour que le produit matriciel fonctionne:
x = np.array([1,2,3])
print("value: {}, dim: {}, shape: {}, A@x: {}, x@A {}".format(x, x.ndim, x.shape, A@x, x@A))

# Ainsi, Pour deux arrays de dimension 1, @ est le produit scalaire.
x@x # équivalent à <x,x> ou x^T * x si x est un vecteur colonne. 

# La syntaxe np.newaxis permet d'augmenter sa dimension. Maintenant, l'orientation ligne ou colonne est importante.
b = x[:,np.newaxis] # `b` est un vecteur colonne (3 lignes, 1 colonne)
print("value: {}, dim: {}, shape: {}, A@b: {}".format(b, b.ndim, b.shape, A@b))

c = x[np.newaxis,:] # `c` est un vecteur ligne (1 ligne, 3 colonnes)
print("value: {}, dim: {}, shape: {}, b@A:{}".format(c, c.ndim, c.shape, c@A))

c@b #=> un scalaire (produit scalaire des deux vecteurs <b,c>).
b@c #=> une matrice 3x3. 

A = np.array([[1.2,0.46],
              [.45,0.8]])
x = [0.7,0.3]

# Voyons, graphiquement, l'effet d'un produit matriciel sur un point:
fig, ax = plt.subplots(figsize=(8,8))
ax.scatter(*x, label=r'$\vec{x}$')
ax.scatter(*A@x, label=r'$A \circ \vec{x}$')
ax.scatter(0,0, label=r'0',color='k')
ax.arrow(0,0,*x,color='C0')
ax.arrow(0,0,*A@x,color='C1')

ax.set(xlim=(0,1.3*np.max(A@x)), ylim=(0,1.3*np.max(A@x)))
ax.legend()
ax.grid(1)

lspace = np.linspace(1,10,5)
fig, ax = plt.subplots(figsize=(10,10))

# Meshgrid donne le produit cartésien des deux vecteurs.
# Et donc de la grille sur laquelle on veut appliquer A.
X,Y = np.meshgrid(lspace,lspace)

# np.vectorize permet de donner à une fonction quelquonque (ici le lambda)
# des capacités de broadcasting (application à tous les éléments).
# Par défaut vectorize attend une fonction qui prend des scalaires et retourne des scalaires.
# Si ce n'est pas le cas il faut préciser la signature. Ici la fonction prend deux scalaires "(),()" et
# retourne un array de taille 2 "->(2)"
f = np.vectorize(lambda x,y: A@[x,y], signature='(),()->(2)')
U = f(X,Y) # On applique la fonction à tous les éléments de X,Y ! 

# L'attribut `flat` d'un array renvoie une vue de l'array "applatie", c'est à dire de dimension 1.
ax.scatter(X.flat, Y.flat, color='C0', label=r'$\vec{x}$')
ax.scatter(U[:,:,0].flat,  U[:,:,1].flat, color='C1', label=r'$A \circ \vec{x}$')
ax.quiver(X, Y, U[:,:,0]-X, U[:,:,1]-Y, angles='xy', units='xy', scale=1, color='C2', alpha=0.4)
ax.set(xlim=(0,np.max(U)),ylim=(0,np.max(U)))
ax.legend()
ax.grid()

def display_transform(A, ax=None, title=None, name='A'):
    """Affiche l'effet de la matrice A sur le carré unité"""
    if ax is None:
        fig, ax = plt.subplots(figsize=(6,6))
    if title is not None:
        ax.set(title=title)
        
    # unit square
    square = [(0,0),(0,1),(1,1),(1,0)]

    # Plot the canoncial basis and its image 
    arrow ={'linewidth':2, 'head_width':0.05, 'length_includes_head':True}
    ax.arrow(0, 0, 0, 1, color='k', label='e', **arrow)
    ax.arrow(0, 0, 1, 0, color='k', **arrow)
    ax.arrow(0, 0, *A@np.array([0,1]), color='C0' , **arrow)
    ax.arrow(0, 0, *A@np.array([1,0]), color='C0', label=name+'e', **arrow)
        
    # Plot the unit square and its image
    Asquare = [A@x for x in square]
    p1=ax.add_patch(Polygon(square, facecolor='k', alpha=0.3))
    p2=ax.add_patch(Polygon(Asquare, facecolor='C0', alpha=0.3))
    
    ## Plot a unit circle
    u = np.linspace(-1,1,1000)
    ax.plot(u, np.sqrt(1-u**2),color='k',alpha=0.1)
    ax.plot(u, -np.sqrt(1-u**2),color='k',alpha=0.1)
    
    # Plot a grid 
    ax.grid(1)
    
    # Informative legend
    tex =  ('${}'.format(name)
            +r' = \left [ \begin{array}{cc}'
            + r"{0[0][0]:2} & {0[0][1]:2} \\ {0[1][0]:2} & {0[1][1]:2}".format(A) 
            + r'  \end{array} \right ] $')
    ax.legend((p1,p2),('$e$','${} \circ e$'.format(name)), title=tex)
    
    # Set axis limits
    ax.set(xlim=(np.floor(np.min([np.min(Asquare),-1])),
                 np.ceil(np.max([np.max(Asquare),1]))), 
           ylim=(np.floor(np.min([np.min(Asquare),-1])),
                 np.ceil(np.max([np.max(Asquare),1]))))
    return ax

display_transform(A);

B = np.array([[1,1],
              [0,1]])  # Cisaillement 
C = np.array([[0,-1],
             [1,0]]) # Rotation de 180 degrés.

fig, ax = plt.subplots(1,4,figsize=(24,6))
display_transform(B, ax[0], title='Cisaillement', name='B')
display_transform(C, ax[1], title='Rotation', name='C') 
display_transform(B@C, ax[2], title='Cisaillement puis rotation', name=r'B\circ C')
display_transform(C@B, ax[3], title='Rotation puis cisaillement', name=r'C\circ B');

invA = np.linalg.inv(A) # Inverse numériquement la matrice. 

fig, ax = plt.subplots(1,3,figsize=(24,8))
display_transform(A, ax[0])
display_transform(invA, ax[1], name='A^{-1}')
display_transform(A@invA, ax[2], name='AA^{-1}'); # Remarquez l'erreur numérique ce n'est pas exactement I

# Regardons quelques matrices et leur déterminant, quand le volume du carré unité passe par 0. 
fig, axes = plt.subplots(1,5, figsize=(5*6,6))
theta = np.linspace(0, np.pi)
j = 0
dots = []
detlist = []
for i,t in enumerate(theta): 
    # Cette famille de matrice ne change pas le second vecteur de la base canonique,
    # et fait décrire le cercle unité au premier.
    M = np.array([[np.cos(t),0],
                  [np.sin(t),1]])
    # np.linalg.det permet de calculer le déterminant. 
    detlist.append(np.linalg.det(M))
    
    if i>1 and i%(len(theta)//4) == 0:
        t = np.round(t, 2)
        M = np.array([[np.cos(t),0],
              [np.sin(t),1]])
        j+=1
        display_transform(M, axes[j], title='t={}'.format(t), name='M')
        dots.append((t,np.linalg.det(M)))
        
axes[0].plot(theta, detlist)
axes[0].scatter(*zip(*dots))
axes[0].hlines(0, theta.min(),theta.max())
axes[0].vlines(0,-1,1)

def display_eigenvectors(A, ax=None):
    """Display eigenvectors of the matrix A"""
    if ax is None:
        ax = plt.gca()

    # Compute eigenvectors
    eigval, eigvect = np.linalg.eig(A)
    
    # Plot eigenvectors
    # If they have complex component, it is impossible to plot this way and a 
    # type error will be raised.
    arrow ={'linewidth':2, 'head_width':0.05, 'length_includes_head':True}
    try:
        ax.arrow(0, 0, eigval[0]*eigvect[0,0], eigval[0]*eigvect[1,0], color='C1', label='e', **arrow)
        u = ax.arrow(0, 0,  eigval[1]*eigvect[0,1],eigval[1]* eigvect[1,1], color='C1', **arrow)
    except TypeError:
        pass
    
    xlim = ax.get_xlim()
    ylim = ax.get_ylim()
    xmax = 10*np.max([ax.get_xlim(),ax.get_ylim()])
    try:
        ax.plot([0,eigvect[0,0]*xmax],[0,eigvect[1,0]*xmax], color='C1',alpha=.3)
        ax.plot([0,eigvect[0,1]*xmax],[0,eigvect[1,1]*xmax], color='C1',alpha=.3)
        ax.plot([0,-eigvect[0,0]*xmax],[0,-eigvect[1,0]*xmax], color='C1',alpha=.3)
        ax.plot([0,-eigvect[0,1]*xmax],[0,-eigvect[1,1]*xmax], color='C1',alpha=.3)
    except TypeError:
        pass
    
    ax.set(xlim=xlim, ylim=ylim)
    return ax

display_transform(A)
display_eigenvectors(A)