---
jupytext:
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.0
kernelspec:
  name: python3
  display_name: Python 3
---

# NumPy — algèbre linéaire et opérations avancées

```{code-cell} python
:tags: [hide-input]

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
import seaborn as sns

sns.set_theme(style="whitegrid", palette="muted", font_scale=1.1)
```

## Manipulation de formes

La capacité à remodeler des tableaux sans copier leurs données est l'une des forces de NumPy. La plupart des opérations de transformation de forme retournent des **vues** — des objets qui partagent le même bloc mémoire que le tableau d'origine. Modifier la vue modifie donc aussi le tableau source.

```{code-cell} python
a = np.arange(24)
print("a :", a)
print("shape :", a.shape)

# reshape : changer la forme sans copier les données
b = a.reshape(4, 6)
print("\nreshape(4, 6) :\n", b)

c = a.reshape(2, 3, 4)
print("\nreshape(2, 3, 4) :\n", c)

# L'argument -1 laisse NumPy inférer la dimension manquante
d = a.reshape(6, -1)    # 6 lignes, nombre de colonnes inféré = 4
print("\nreshape(6, -1) :\n", d)
```

```{code-cell} python
# ravel et flatten : mise à plat
a = np.array([[1, 2, 3], [4, 5, 6]])

vue = a.ravel()          # vue (si possible) — modification répercutée sur a
copie = a.flatten()      # toujours une copie

print("ravel()   :", vue)
print("flatten() :", copie)

# Modifier la vue modifie le tableau source
vue[0] = 99
print("a après modification de vue :", a)
```

```{code-cell} python
# transpose et T
M = np.arange(12).reshape(3, 4)
print("M :\n", M)
print("M.T :\n", M.T)

# Pour un tableau 3D, transpose accepte les axes dans l'ordre voulu
T3 = np.ones((2, 3, 4))
print("\nT3.shape :", T3.shape)
print("T3.transpose(2, 0, 1).shape :", T3.transpose(2, 0, 1).shape)
```

```{code-cell} python
# np.newaxis : ajouter une dimension de taille 1
v = np.array([1, 2, 3])          # forme (3,)
print("v.shape         :", v.shape)
print("v[:, np.newaxis] :", v[:, np.newaxis].shape)  # (3, 1) — vecteur colonne
print("v[np.newaxis, :] :", v[np.newaxis, :].shape)  # (1, 3) — vecteur ligne

# squeeze : supprimer les dimensions de taille 1
a = np.zeros((1, 3, 1, 4))
print("\nBefore squeeze :", a.shape)
print("After squeeze  :", np.squeeze(a).shape)    # (3, 4)
print("squeeze(axis=2):", np.squeeze(a, axis=2).shape)  # (1, 3, 4)
```

```{note}
`np.newaxis` est simplement `None` — les deux notations sont interchangeables. Il est cependant fortement conseillé d'utiliser `np.newaxis` pour la lisibilité, car l'intention est immédiatement claire. L'ajout de dimensions avec `np.newaxis` est un idiome fréquent pour rendre deux tableaux compatibles au broadcasting : transformer un vecteur de forme `(n,)` en `(n, 1)` ou `(1, n)` avant une opération.
```

## Empilement et découpage

Assembler plusieurs tableaux en un seul, ou au contraire diviser un tableau en morceaux, sont des opérations très fréquentes lors du prétraitement de données.

```{code-cell} python
a = np.array([[1, 2], [3, 4]])
b = np.array([[5, 6], [7, 8]])

# Concaténation selon un axe existant
print("concatenate axis=0 :\n", np.concatenate([a, b], axis=0))  # empiler verticalement
print("concatenate axis=1 :\n", np.concatenate([a, b], axis=1))  # côte à côte
```

```{code-cell} python
# vstack, hstack, dstack — raccourcis expressifs
v = np.array([1, 2, 3])
w = np.array([4, 5, 6])

print("vstack :\n", np.vstack([v, w]))   # (2, 3)
print("hstack :", np.hstack([v, w]))     # (6,)
print("column_stack :\n", np.column_stack([v, w]))  # (3, 2)
```

```{code-cell} python
# np.stack : créer une NOUVELLE dimension
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

print("stack axis=0 :\n", np.stack([a, b], axis=0))  # (2, 3)
print("stack axis=1 :\n", np.stack([a, b], axis=1))  # (3, 2)
```

```{code-cell} python
# np.split : découper un tableau
M = np.arange(24).reshape(4, 6)
parties = np.split(M, 2, axis=0)   # 2 parties égales selon les lignes
print("split axis=0 :")
for p in parties:
    print(p)

colonnes = np.hsplit(M, 3)         # 3 parties égales selon les colonnes
print("\nhsplit en 3 :")
for c in colonnes:
    print(c)
```

## Algèbre linéaire avec `np.linalg`

L'algèbre linéaire est au cœur de la quasi-totalité des algorithmes de machine learning et de deep learning. NumPy fournit un sous-module `np.linalg` complet, qui délègue les calculs aux bibliothèques BLAS et LAPACK — des implémentations hautement optimisées, souvent accélérées matériellement.

```{code-cell} python
A = np.array([[2.0, 1.0, -1.0],
              [-3.0, -1.0, 2.0],
              [-2.0, 1.0, 2.0]])
b = np.array([8.0, -11.0, -3.0])

# Résolution du système linéaire Ax = b
x = np.linalg.solve(A, b)
print("Solution Ax = b :", x)
print("Vérification A @ x :", A @ x)
```

```{code-cell} python
# Produit matriciel : l'opérateur @ (Python 3.5+)
A = np.random.default_rng(0).standard_normal((3, 4))
B = np.random.default_rng(1).standard_normal((4, 5))
C = A @ B
print("A @ B shape :", C.shape)

# Autres opérations linalg
M = np.array([[4.0, 2.0], [1.0, 3.0]])
print("\nDéterminant :", np.linalg.det(M))
print("Inverse :\n", np.linalg.inv(M))
print("Norme (Frobenius) :", np.linalg.norm(M))
print("Rang :", np.linalg.matrix_rank(M))
```

```{code-cell} python
# Valeurs propres et vecteurs propres
M = np.array([[3.0, 1.0], [1.0, 3.0]])
valeurs, vecteurs = np.linalg.eig(M)
print("Valeurs propres :", valeurs)
print("Vecteurs propres (en colonnes) :\n", vecteurs)

# Vérification : M @ v = lambda * v
for i in range(len(valeurs)):
    v = vecteurs[:, i]
    print(f"M @ v{i} = {M @ v},  λ{i} * v{i} = {valeurs[i] * v}")
```

```{admonition} Décomposition en valeurs singulières (SVD)
:class: tip
La **décomposition en valeurs singulières** (*Singular Value Decomposition*, SVD) factorise toute matrice M de taille (m, n) en trois matrices : M = U Σ Vᵀ, où U est orthogonale (m, m), Σ est diagonale (m, n) avec des valeurs non négatives décroissantes appelées **valeurs singulières**, et Vᵀ est orthogonale (n, n). La SVD est fondamentale en data science : elle sous-tend l'ACP (*Analyse en Composantes Principales*), la recommandation collaborative, la compression d'images et la pseudo-inverse.
```

```{code-cell} python
# SVD complète
M = np.random.default_rng(42).standard_normal((4, 6))
U, sigma, Vt = np.linalg.svd(M, full_matrices=True)
print("M shape :", M.shape)
print("U shape :", U.shape)
print("sigma   :", sigma)
print("Vt shape:", Vt.shape)

# Reconstruction exacte
Sigma = np.zeros_like(M)
np.fill_diagonal(Sigma, sigma)
M_rec = U @ Sigma @ Vt
print("\nErreur de reconstruction :", np.max(np.abs(M - M_rec)))
```

```{code-cell} python
:tags: [hide-input]

# Illustration : compression d'image par SVD (approximation de rang k)
rng = np.random.default_rng(7)

# Créer une "image" synthétique (gradient + bruit structuré)
x = np.linspace(0, 1, 64)
y = np.linspace(0, 1, 64)
X, Y = np.meshgrid(x, y)
image = np.sin(4 * np.pi * X) * np.cos(3 * np.pi * Y) + 0.3 * rng.standard_normal((64, 64))

U, sigma, Vt = np.linalg.svd(image, full_matrices=False)

rangs = [1, 3, 8, 20, 64]
fig, axes = plt.subplots(1, len(rangs), figsize=(15, 3.5))

for ax, k in zip(axes, rangs):
    approx = U[:, :k] @ np.diag(sigma[:k]) @ Vt[:k, :]
    ax.imshow(approx, cmap='RdBu_r', aspect='auto',
              vmin=image.min(), vmax=image.max())
    ax.set_title(f"Rang k = {k}", fontsize=11, fontweight='bold')
    ax.set_xticks([])
    ax.set_yticks([])

fig.suptitle("Décomposition SVD — reconstruction progressive d'une image synthétique",
             fontsize=13, fontweight='bold', y=1.04)
plt.tight_layout()
plt.show()
```

## Statistiques

NumPy fournit un ensemble complet de fonctions statistiques descriptives. La plupart acceptent un argument `axis` pour calculer la statistique le long d'une dimension donnée plutôt que sur l'ensemble du tableau.

```{code-cell} python
rng = np.random.default_rng(0)
data = rng.standard_normal((5, 4))   # 5 observations, 4 variables
print("Données :\n", np.round(data, 3))
print()

# Statistiques globales
print("Moyenne globale :", np.round(data.mean(), 4))
print("Écart-type global:", np.round(data.std(), 4))
print()

# Statistiques par axe
print("Moyenne par colonne (axis=0) :", np.round(data.mean(axis=0), 4))
print("Moyenne par ligne  (axis=1) :", np.round(data.mean(axis=1), 4))
print()

# Autres statistiques
print("Médiane :", np.round(np.median(data), 4))
print("Percentiles [25, 50, 75] :", np.round(np.percentile(data, [25, 50, 75]), 4))
print("Minimum :", data.min(), "| Maximum :", data.max())
print("Argmin  :", data.argmin(), "| Argmax  :", data.argmax())
```

```{code-cell} python
# Corrélation et histogramme
x = rng.standard_normal(100)
y = 0.7 * x + 0.3 * rng.standard_normal(100)

corrcoef = np.corrcoef(x, y)
print("Matrice de corrélation :\n", np.round(corrcoef, 4))

counts, bins = np.histogram(x, bins=10)
print("\nHistogramme de x :")
print("Comptages :", counts)
print("Bords des bacs :", np.round(bins, 2))
```

## Entrées et sorties

NumPy propose plusieurs fonctions pour persister des tableaux sur disque. Le choix du format dépend du contexte : interopérabilité avec d'autres outils, taille des données, nécessité de compression.

```{code-cell} python
import tempfile
import os

# Créer un tableau à sauvegarder
a = np.arange(12).reshape(3, 4)

with tempfile.TemporaryDirectory() as tmpdir:
    # Format binaire .npy (un seul tableau)
    path_npy = os.path.join(tmpdir, "tableau.npy")
    np.save(path_npy, a)
    a_chargé = np.load(path_npy)
    print("Chargé depuis .npy :\n", a_chargé)

    # Format .npz (plusieurs tableaux, avec compression optionnelle)
    b = np.linspace(0, 1, 6)
    path_npz = os.path.join(tmpdir, "multi.npz")
    np.savez(path_npz, matrice=a, vecteur=b)
    données = np.load(path_npz)
    print("\nClés dans .npz :", list(données.keys()))
    print("matrice :\n", données["matrice"])
    print("vecteur :", données["vecteur"])

    # Format texte .csv
    path_csv = os.path.join(tmpdir, "matrice.csv")
    np.savetxt(path_csv, a, delimiter=",", fmt="%d", header="c0,c1,c2,c3")
    a_txt = np.loadtxt(path_csv, delimiter=",", skiprows=1)
    print("\nChargé depuis .csv :\n", a_txt)
```

```{note}
Le format `.npy` est le plus rapide à lire et écrire car il ne fait aucune conversion : les octets du tableau sont écrits directement sur disque avec un en-tête minimal décrivant le `dtype` et la `shape`. Pour des données volumineuses, préférer `np.savez_compressed` (format `.npz` avec compression `zlib`) ou les formats de haut niveau comme **HDF5** (via `h5py`) ou **Zarr**, qui supportent l'accès partiel aux données et la lecture en parallèle.
```

## Performance : vues, contiguïté et `np.einsum`

Maîtriser les mécanismes de performance de NumPy permet d'écrire du code qui tire pleinement parti du matériel, notamment pour des calculs intensifs sur de grands tableaux.

**Vues vs copies.** Une vue partage la même zone mémoire que le tableau source. Les slices de base créent des vues ; le *fancy indexing* crée des copies. On peut vérifier si deux tableaux partagent la même mémoire avec `np.shares_memory(a, b)`.

```{code-cell} python
a = np.arange(12).reshape(3, 4)

vue = a[1:, ::2]        # slice → vue
copie = a[[0, 2], :]    # fancy indexing → copie

print("Vue partage la mémoire :", np.shares_memory(a, vue))
print("Copie partage la mémoire :", np.shares_memory(a, copie))
```

**Contiguïté mémoire.** Un tableau est dit **C-contigu** si les éléments sont ordonnés ligne par ligne (ordre C), et **Fortran-contigu** s'ils sont ordonnés colonne par colonne. NumPy crée par défaut des tableaux C-contigus. Les opérations qui traversent le tableau dans l'ordre de sa contiguïté sont plus rapides, car elles génèrent moins de défauts de cache.

```{code-cell} python
a = np.arange(12).reshape(3, 4)
print("C-contigu :", a.flags['C_CONTIGUOUS'])
print("F-contigu :", a.flags['F_CONTIGUOUS'])

a_f = np.asfortranarray(a)
print("\nAprès asfortranarray :")
print("C-contigu :", a_f.flags['C_CONTIGUOUS'])
print("F-contigu :", a_f.flags['F_CONTIGUOUS'])
```

**`np.einsum` — contraction tensorielle.** `einsum` est une interface générale pour exprimer des sommes de produits sur des tableaux multidimensionnels via la **notation d'Einstein**. Elle permet d'exprimer en une seule ligne des opérations aussi variées que la transposition, le produit matriciel, la trace, les produits externes ou les contractions tensorielles.

```{code-cell} python
A = np.random.default_rng(0).standard_normal((3, 4))
B = np.random.default_rng(1).standard_normal((4, 5))

# Produit matriciel : A_{ik} B_{kj} → C_{ij}
C_einsum = np.einsum('ik,kj->ij', A, B)
C_matmul = A @ B
print("Erreur max A@B vs einsum :", np.max(np.abs(C_einsum - C_matmul)))

# Trace d'une matrice
M = np.random.default_rng(2).standard_normal((4, 4))
print("Trace via einsum :", np.einsum('ii', M))
print("Trace via linalg :", np.trace(M))

# Produit externe de deux vecteurs
u = np.array([1.0, 2.0, 3.0])
v = np.array([4.0, 5.0])
print("Produit externe u⊗v :\n", np.einsum('i,j->ij', u, v))
```

```{admonition} Quand utiliser einsum ?
:class: note
`np.einsum` est particulièrement utile lorsque l'opération souhaitée n'a pas d'équivalent direct dans NumPy, ou lorsque l'enchaînement de plusieurs opérations (`matmul`, `sum`, `transpose`) produirait des tableaux intermédiaires inutiles. Par exemple, le calcul de la trace du produit AB, `np.trace(A @ B)`, alloue la matrice entière AB avant de sommer la diagonale. L'équivalent `np.einsum('ij,ji->', A, B)` effectue le calcul sans allouer AB. Sur de grands tableaux, l'argument `optimize=True` permet à NumPy de chercher le chemin de contraction optimal.
```

## Résumé

Ce chapitre a approfondi NumPy au-delà des bases, en couvrant les opérations avancées qui constituent le quotidien du travail en data science :

- La **manipulation des formes** — `reshape`, `ravel`, `flatten`, `transpose`, `np.newaxis`, `squeeze` — permet de réorganiser les données sans copie, en exploitant le mécanisme des vues et des `strides`.
- L'**empilement et le découpage** — `concatenate`, `stack`, `vstack`, `hstack`, `split` — permettent d'assembler et de désassembler des tableaux pour constituer des ensembles de données ou traiter des sous-blocs.
- Le module **`np.linalg`** fournit toute l'algèbre linéaire nécessaire : résolution de systèmes (`solve`), déterminant, inverse, valeurs propres (`eig`) et décomposition SVD. L'opérateur `@` est la notation recommandée pour le produit matriciel.
- Les **fonctions statistiques** — `mean`, `std`, `var`, `median`, `percentile`, `corrcoef`, `histogram` — acceptent toutes un argument `axis` pour opérer sur une dimension précise.
- Les formats **`.npy` et `.npz`** sont les formats natifs de NumPy, rapides et fiables pour persister des tableaux. `np.savetxt` et `np.loadtxt` servent pour l'interopérabilité avec les outils texte.
- Comprendre la distinction **vue vs copie** et la **contiguïté mémoire** est essentiel pour éviter des bugs subtils et maximiser les performances. `np.einsum` offre une syntaxe élégante pour les contractions tensorielles complexes.

Dans le chapitre suivant, nous passerons à la visualisation avec Matplotlib et Seaborn, pour donner une forme graphique aux données que NumPy nous permet de manipuler si efficacement.