05-ClassificationsSupervisees.qmd

---
jupyter: python3
from: markdown+emoji
execute:
  echo: true
  eval: true
  message: false
  warning: false
---

# Classifications d'images supervisées {#sec-chap05}


## :rocket: Préambule

Assurez-vous de lire ce préambule avant d'exécutez le reste du notebook.

### :dart: Objectifs
Dans ce chapitre, nous ferons une introduction générale à l'apprentissage automatique et abordons quelques techniques fondamentales. La librairie centrale utilisée dans ce chapitre sera [`sickit-learn`](https://scikit-learn.org/). Ce chapitre est aussi disponible sous la forme d'un notebook Python sur Google Colab:

[![](images/colab-badge.svg)](https://colab.research.google.com/github/sfoucher/TraitementImagesPythonVol1/blob/main/notebooks/05-ClassificationsSupervisees.ipynb){target="_blank"} 

### Librairies
Les librairies utilisées dans ce chapitre sont les suivantes:

* [SciPy](https://scipy.org/)

* [NumPy](https://numpy.org/) 

* [opencv-python · PyPI](https://pypi.org/project/opencv-python/)

* [scikit-image](https://scikit-image.org/)

* [Rasterio](https://rasterio.readthedocs.io/en/stable/)

* [xarray](https://docs.xarray.dev/en/stable/)

* [rioxarray](https://corteva.github.io/rioxarray/stable/index.html)

* [geopandas](https://geopandas.org)

* [scikit-learn](https://scikit-learn.org/)

Dans l'environnement Google Colab, seul `rioxarray` et `xrscipy` doit être installés:
```{python}
#| eval: false
%%capture
!pip install -qU matplotlib rioxarray xrscipy
```
Vérifier les importations nécessaires en premier:
```{python}
import numpy as np
import rioxarray as rxr
from scipy import signal
import xarray as xr
import rasterio
import xrscipy
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import geopandas
from shapely.geometry import Point
import pandas as pd
from numba import jit
from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import train_test_split
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.inspection import DecisionBoundaryDisplay
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis, QuadraticDiscriminantAnalysis
from sklearn.datasets import make_blobs, make_classification, make_gaussian_quantiles
```
### Images utilisées

Nous allons utilisez les images suivantes dans ce chapitre:
```{python}
#| eval: false
%%capture
import gdown

gdown.download('https://drive.google.com/uc?export=download&confirm=pbef&id=1a6Ypg0g1Oy4AJt9XWKWfnR12NW1XhNg_', output= 'RGBNIR_of_S2A.tif')
gdown.download('https://drive.google.com/uc?export=download&confirm=pbef&id=1a4PQ68Ru8zBphbQ22j0sgJ4D2quw-Wo6', output= 'landsat7.tif')
gdown.download('https://drive.google.com/uc?export=download&confirm=pbef&id=1_zwCLN-x7XJcNHJCH6Z8upEdUXtVtvs1', output= 'berkeley.jpg')
gdown.download('https://drive.google.com/uc?export=download&confirm=pbef&id=1dM6IVqjba6GHwTLmI7CpX8GP2z5txUq6', output= 'SAR.tif')
gdown.download('https://drive.google.com/uc?export=download&confirm=pbef&id=1aAq7crc_LoaLC3kG3HkQ6Fv5JfG0mswg', output= 'carte.tif')
```
Vérifiez que vous êtes capable de les lire :
```{python}
#| output: false

with rxr.open_rasterio('berkeley.jpg', mask_and_scale= True) as img_rgb:
    print(img_rgb)
with rxr.open_rasterio('RGBNIR_of_S2A.tif', mask_and_scale= True) as img_rgbnir:
    print(img_rgbnir)
with rxr.open_rasterio('SAR.tif', mask_and_scale= True) as img_SAR:
    print(img_SAR)
with rxr.open_rasterio('carte.tif', mask_and_scale= True) as img_carte:
    print(img_carte)
```

## Principes généraux

Une classification supervisée ou dirigée consiste à attribuer une étiquette (une classe) de manière automatique à chaque point d'un jeu de données. Cette classification peut se faire à l'aide d'une cascade de règles pré-établies (arbre de décision) ou à l'aide de techniques d'apprentissage automatique (*machine learning*). L'utilisation de règles pré-établies atteint vite une limite car ces règles doivent être fournies manuellement par un expert. Ainsi, l'avantage de l'apprentissage automatique est que les règles de décision sont dérivées automatiquement du jeu de données via une phase dite d’entraînement. On parle souvent de solutions générées par les données (*Data Driven Solutions*). Cet ensemble de règles est souvent appelé **modèle**. On visualise souvent ces règles sous la forme de *frontières de décisions* dans l'espace des données. Cependant, un des défis majeur de ce type de technique est d'être capable de produire des règles qui soient généralisables au-delà du jeu d’entraînement. 

Les classifications supervisées ou dirigées présupposent donc que nous avons à disposition **un jeu d’entraînement** déjà étiqueté. Celui-ci va nous permettre de construire un modèle. Afin que ce modèle soit représentatif et robuste, il nous faut assez de données d’entraînement. Les algorithmes d'apprentissage automatique sont très nombreux et plus ou moins complexes pouvant produire des frontières de décision très complexes et non linéaires. 

**curse of dimensionnality, capacité d'un modèle, sur-aprrentissage, sous-apprentissage**

### Comportement d'un modèle

Cet exemple tiré de [`sickit-learn`](https://scikit-learn.org/stable/auto_examples/model_selection/plot_underfitting_overfitting.html#sphx-glr-auto-examples-model-selection-plot-underfitting-overfitting-py) illustre les problèmes d'ajustement insuffisant ou **sous-apprentissage** (*underfitting*) et d'ajustement excessif ou **sur-apprentissage** (*overfitting*) et montre comment nous pouvons utiliser la régression linéaire avec un modèle polynomiale pour approximer des fonctions non linéaires. La @fig-overfitting montre la fonction que nous voulons approximer, qui est une partie de la fonction cosinus (couleur orange). En outre, les échantillons de la fonction réelle et les approximations de différents modèles sont affichés en bleu. Les modèles ont des caractéristiques polynomiales de différents degrés. Nous pouvons constater qu'une fonction linéaire (polynôme de degré 1) n'est pas suffisante pour s'adapter aux échantillons d'apprentissage. C'est ce qu'on appelle un sous-ajustement (underfitting) qui produit un biais systématique quelque soit les points d’entraînement. Un polynôme de degré 4 se rapproche presque parfaitement de la fonction réelle. Cependant, pour des degrés plus élevés, le modèle s'adaptera trop aux données d'apprentissage, c'est-à-dire qu'il apprendra le bruit des données d'apprentissage. Nous évaluons quantitativement le sur-apprentissage et le sous-apprentissage à l'aide de la validation croisée. Nous calculons l'erreur quadratique moyenne (EQM) sur l'ensemble de validation. Plus elle est élevée, moins le modèle est susceptible de se généraliser correctement à partir des données d'apprentissage.

```{python}
#| echo: false
#| label: fig-overfitting
#| fig-cap: "Exemples de sur et sous-apprentissage."
import matplotlib.pyplot as plt
import numpy as np

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import cross_val_score, cross_validate
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures


def true_fun(X):
    return np.cos(1.5 * np.pi * X+np.pi/2)


np.random.seed(0)
noise_level= 0.1
n_samples = 30
degrees = [1,4, 15]
#degrees= range(1,16)
X = np.sort(np.random.rand(n_samples))
y = true_fun(X) + np.random.randn(n_samples) * noise_level

X_test = np.sort(np.random.rand(10))
y_test = true_fun(X_test) + np.random.randn(int(10)) * noise_level

plt.figure(figsize=(14, 5))
results= []
for i in range(len(degrees)):
    ax = plt.subplot(1, len(degrees), i + 1)
    plt.setp(ax, xticks=(), yticks=())

    polynomial_features = PolynomialFeatures(degree=degrees[i], include_bias=False)
    linear_regression = LinearRegression()
    pipeline = Pipeline(
        [
            ("polynomial_features", polynomial_features),
            ("linear_regression", linear_regression),
        ]
    )
    pipeline.fit(X[:, np.newaxis], y)

    # Evaluate the models using crossvalidation
    scores = cross_validate(
        pipeline, X[:, np.newaxis], y, scoring="neg_mean_squared_error", cv=10, return_train_score=True
    )

    X_true = np.linspace(0, 1, 100)
    plt.plot(X_true, pipeline.predict(X_true[:, np.newaxis]), label="Modèle")
    plt.plot(X_true, true_fun(X_true), label="Vraie fonction")
    plt.scatter(X, y, edgecolor="b", s=20, label="Entr.")
    plt.scatter(X_test, y_test, edgecolor="g", s=20, label="Test")
    plt.xlabel("x")
    plt.ylabel("y")
    plt.xlim((0, 1))
    plt.ylim((-2, 2))
    plt.legend(loc="best")
    plt.title(
        "Degré {}\nErreur = {:.1e}(+/- {:.1e})".format(
            degrees[i], -scores['test_score'].mean(), scores['test_score'].std()
        ), fontsize='small'
    )
    results.append([degrees[i], -scores['train_score'].mean(), -scores['test_score'].mean(),scores['train_score'].std(),scores['test_score'].std()])
plt.show()
```

On constate aussi que sans les échantillons de validation, nous serions incapable de déterminer la situation de sur-apprentissage, l'erreur sur les points d’entraînement seul étant excellente pour un degré 15.

### Pipeline

La construction d'un modèle implique généralement toujours les mêmes étapes illustrées sur la figure [@fig-pipeline]:

1. La préparation des données implique parfois un pré-traitement afin de normaliser les données.

2. Partage des données en trois groupes: entraînement, validation et test

3. L'apprentissage du modèle sur l'ensemble d'entraînement. Cet apprentissage nécessite de déterminer les valeurs des hyper-paramètres du modèle par l'usager.

4. La validation du modèle sur l'ensemble de validation. Cette étape vise à vérifier que les hyper-paramètres du modèle sont adéquate.

5. Enfin le test du modèle sur un ensemble de donnée indépendant

```{mermaid}
%%| echo: false
%%| label: fig-pipeline
%%| fig-cap: "Étapes standards dans un entraînement."
flowchart TD
    A[fa:fa-database Données] --> B(fa:fa-gear Prétraitement)
    B --> C(fa:fa-folder-tree Partage des données) -.-> D(fa:fa-gears Entraînement)
    H[[Hyper-paramètres]] --> D
    D --> |Modèle| E>Validation]
    E --> |Modèle| G>Test]
    C -.-> E
    C -.-> G
```

### Construction d'un ensemble d’entraînement {#sec-05.02.02}

Les données d’entraînement vont permettre de construire un modèle. Ces données peuvent prendre des formes très variées mais on peut voir cela sous la forme d'un tableau $N \times D$:

1. La taille $N$ du jeu de donnée

2. Chaque entrée définit un échantillon ou un point dans un espace à plusieurs dimension.

3. Chaque échantillon est décrit par $D$ dimensions ou caractéristiques (*features*).

Une façon simple de construire un ensemble d’entraînement est d'échantillonner un produit existant. Nous allons utiliser la carte d'occupation des sols suivante qui contient 12 classes différentes.

```{python}
couleurs_classes= {'NoData': 'black', 'Commercial': 'yellow', 'Nuages': 'lightgrey', 
                    'Foret': 'darkgreen', 'Faible_végétation': 'green', 'Sol_nu': 'saddlebrown',
                  'Roche': 'dimgray', 'Route': 'red', 'Urbain': 'orange', 'Eau': 'blue', 'Tourbe': 'salmon', 'Végétation éparse': 'darkgoldenrod', 'Roche avec végétation': 'darkseagreen'}
nom_classes= [*couleurs_classes.keys()]
couleurs_classes= [*couleurs_classes.values()]
```
On peut visualiser la carte de la façon suivante:
```{python}
import matplotlib.pyplot as plt
import rioxarray as rxr
cmap_classes = ListedColormap(couleurs_classes)

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))
img_carte.squeeze().plot.imshow(cmap=cmap_classes, vmin=0, vmax=12)
ax.set_title("Carte d'occupation des sols", fontsize="small")
```

On peut facilement calculer la fréquence d’occurrence des 12 classes dans l'image à l'aide de `numpy`:

```{python}
img_carte= img_carte.squeeze() # nécessaire pour ignorer la dimension du canal
compte_classe = np.unique(img_carte.data, return_counts=True)
print(compte_classe)
```
La fréquence d'apparition de chaque classe varie grandement, on parle alors d'un **ensemble déséquilibré**. Ceci est très commun dans la plupart des ensembles d’entraînement, les classes n'apparaissent pas avec la même fréquence. 
```{python}
valeurs, comptes = compte_classe

# Create the histogram
plt.figure(figsize=(5, 3))
plt.bar(valeurs, comptes/comptes.sum()*100)
plt.xlabel("Classes")
plt.ylabel("%")
plt.title("Fréquences des classes", fontsize="small")
plt.xticks(range(len(nom_classes)), nom_classes, rotation=45, ha='right')
plt.show()
```

On peut échantillonner 100 points aléatoires pour chaque classe:

```{python}

img_carte= img_carte.squeeze()
class_counts = np.unique(img_carte.data, return_counts=True)

# Liste vide des points échantillonnées
sampled_points = []
class_labels= [] # contient les étiquettes des classes
for class_label in range(1,13): # pour chacune des 12 classes
  # On cherche tous les pixels pour cette étiquette
  class_pixels = np.argwhere(img_carte.data == class_label)

  # On se limite à 100 pixels par classe
  n_samples = min(100, len(class_pixels))

  # On les choisit les positions aléatoirement
  np.random.seed(0) # ceci permet de répliquer le tirage aléatoire
  sampled_indices = np.random.choice(len(class_pixels), n_samples, replace=False)

  # On prends les positions en lignes, colonnes
  sampled_pixels = class_pixels[sampled_indices]

  # On ajoute les points à la liste
  sampled_points.extend(sampled_pixels)
  class_labels.extend(np.array([class_label]*n_samples)[:,np.newaxis])

# Conversion en NumPy array
sampled_points = np.array(sampled_points)
class_labels = np.array(class_labels)
# On peut naviguer les points à l'aide de la géoréférence
transformer = rasterio.transform.AffineTransformer(img_carte.rio.transform())
transform_sampled_points= transformer.xy(sampled_points[:,0], sampled_points[:,1])

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))
img_carte.squeeze().plot.imshow(cmap=cmap_classes, vmin=0, vmax=12)
ax.scatter(transform_sampled_points[0], transform_sampled_points[1], c='w', s=1)  # Plot sampled points
ax.set_title("Carte d'occupation des sols avec les points échantillonnés", fontsize="small")
plt.show()
```

Une fois les points sélectionnés, il faut ajouter les valeurs des bandes provenant d'une image satellite. Pour cela, on peut utiliser la méthodes `sample()` de `rasterio`. Éventuellement, la librairie [`geopandas`](https://geopandas.org) permet de gérer les données d’entraînement sous la forme d'un tableau transportant aussi l'information de géoréférence. Afin de pouvoir classifier ces points, nous allons ajouter les valeurs radiométriques provenant de l'image Sentinel-2 à 4 bandes `RGBNIR_of_S2A.tif`. Ces valeurs seront stockées dans la colonne `value` sous la forme d'un vecteur en format `string`:

```{python}
points = [Point(xy) for xy in zip(transform_sampled_points[0], transform_sampled_points[1])]
gdf = geopandas.GeoDataFrame(range(1,len(points)+1), geometry=points, crs=img_carte.rio.crs)
coord_list = [(x, y) for x, y in zip(gdf["geometry"].x, gdf["geometry"].y)]
with rasterio.open('RGBNIR_of_S2A.tif') as src:
  gdf["value"] = [x for x in src.sample(coord_list)]
gdf['class']= class_labels
gdf.to_csv('sampling_points.csv') # sauvegarde sous forme d'un format csv
gdf.head()
```

## Analyse préliminaire des données

Une bonne pratique avant d'appliquer une technique d'apprentissage automatique est de regarder les caractéristiques de vos données:

1. Le nombre de dimensions (*features*)

2. Certaines dimensions sont informatives (discriminantes) et d'autres ne le sont pas

3. Le nombre classes

4. Le nombre de modes (*clusters*) par classes

5. Le nombre d'échantillons par classe

6. La forme des groupes 

7. La séparabilité des classes ou des groupes

Une manière d'évaluer la séparabilité de vos classes est d'appliquer des modèles Gaussiens sur chacune des classes. Le modèle Gaussien multivarié suppose que les données sont distribuées comme un nuage de points symétrique et unimodale. La distribution d'un point ${\bold x}$ appartenant à la classe $i$ est la suivante:

$$
P(\bold{x} | Classe=i) = \frac{1}{(2\pi)^{D/2} |\Sigma_i|^{1/2}}\exp\left(-\frac{1}{2} (\bold{x}-\bold{m}_i)^t \Sigma_k^{-1} (\bold{x}-\bold{m}_i)\right)
$$

La méthode [`QuadraticDiscriminantAnalysis`](https://scikit-learn.org/stable/modules/generated/sklearn.discriminant_analysis.QuadraticDiscriminantAnalysis.html) permet de calculer les paramètres des Gaussiennes multivariées pour chacune des classes.

On peut calculer une distance entre deux nuages Gaussiens avec la distance dites de Jeffries-Matusita (JM) basée sur la distance de Bhattacharyya $B$:

$$
JM_{ij}= 2(1-e^{-B})\\
B=\frac{1}{8}({\bold m}_i-\bold{m}_j)^t { \frac{\Sigma_i+\Sigma_j}{2} }(\bold{m}_i-\bold{m}_j)+\frac{1}{2}ln { \frac{|(\Sigma_i+\Sigma_j)/2|}{|\Sigma_i|^{1/2}|\Sigma_j|^{1/2}}}
$$

Cette distance présuppose que chaque classe $i$ est décrite par son centre $\bold{m}_i$ et de sa dispersion dans l'espace à $D$ dimensions mesurée par la matrice de covariance $\Sigma_i$. On peut en faire facilement une fonction Python à l'aide de `numpy`:
```{python}
def bhattacharyya_distance(m1, s1, m2, s2):
    # Calcul de la covariance moyenne
    s = (s1 + s2) / 2
    
    # Calcul du premier terme (différence des moyennes)
    m_diff = m1 - m2
    term1 = np.dot(np.dot(m_diff.T, np.linalg.inv(s)), m_diff) / 8
    
    # Calcul du second terme (différence de covariances)
    term2 = 0.5 * np.log(np.linalg.det(s) / np.sqrt(np.linalg.det(s1) * np.linalg.det(s2)))
    
    return term1 + term2

def jeffries_matusita_distance(m1, s1, m2, s2):
    B = bhattacharyya_distance(m1, s1, m2, s2)
    return 2 * (1 - np.exp(-B))
```

La figure ci-dessous illustre différentes situations avec des données artificielles:

```{python}
#| echo: false
#| warning: false
import matplotlib.pyplot as plt

from sklearn.datasets import make_blobs, make_classification, make_gaussian_quantiles


cmap_classes = ListedColormap( ['blue','yellow','red'])

np.random.seed(42)
plt.figure(figsize=(8, 8))
plt.subplots_adjust(bottom=0.05, top=0.9, left=0.05, right=0.95)

plt.subplot(321)

X1, Y1 = make_classification(
    n_features=2, n_redundant=0, n_informative=1, n_clusters_per_class=1
)
qda = QuadraticDiscriminantAnalysis(store_covariance=True)
qda.fit(X1, Y1)
jm=jeffries_matusita_distance(qda.means_[0], qda.covariance_[0], qda.means_[1], qda.covariance_[1])
plt.title(f"Une dimension informative, un mode par classe JM_12={jm:3.2f}", fontsize="small")
plt.scatter(X1[:, 0], X1[:, 1], marker="o", c=Y1, s=25, edgecolor="k", cmap=cmap_classes)

plt.subplot(322)

X1, Y1 = make_classification(
    n_features=2, n_redundant=0, n_informative=2, n_clusters_per_class=1
)
qda = QuadraticDiscriminantAnalysis(store_covariance=True)
qda.fit(X1, Y1)
jm=jeffries_matusita_distance(qda.means_[0], qda.covariance_[0], qda.means_[1], qda.covariance_[1])
plt.title(f"Deux dimensions informatives, un mode par classe JM_12={jm:3.2f}", fontsize="small")
plt.scatter(X1[:, 0], X1[:, 1], marker="o", c=Y1, s=25, edgecolor="k", cmap=cmap_classes)

plt.subplot(323)

X2, Y2 = make_classification(n_features=2, n_redundant=0, n_informative=2)
qda = QuadraticDiscriminantAnalysis(store_covariance=True)
qda.fit(X1, Y1)
jm=jeffries_matusita_distance(qda.means_[0], qda.covariance_[0], qda.means_[1], qda.covariance_[1])
plt.title(f"Deux dimensions informatives, deux modes par classe JM_12={jm:3.2f}", fontsize="small")
plt.scatter(X2[:, 0], X2[:, 1], marker="o", c=Y2, s=25, edgecolor="k", cmap=cmap_classes)

plt.subplot(324)

X1, Y1 = make_classification(
    n_features=2, n_redundant=0, n_informative=2, n_clusters_per_class=1, n_classes=3
)
qda = QuadraticDiscriminantAnalysis(store_covariance=True)
qda.fit(X1, Y1)
jm=jeffries_matusita_distance(qda.means_[0], qda.covariance_[0], qda.means_[1], qda.covariance_[1])
plt.title(f"Trois classes, deux dimensions informatives, un mode JM_12={jm:3.2f}", fontsize="small")
plt.scatter(X1[:, 0], X1[:, 1], marker="o", c=Y1, s=25, edgecolor="k", cmap=cmap_classes)

plt.subplot(325)

X1, Y1 = make_blobs(n_features=2, centers=3)
qda = QuadraticDiscriminantAnalysis(store_covariance=True)
qda.fit(X1, Y1)
jm=jeffries_matusita_distance(qda.means_[0], qda.covariance_[0], qda.means_[1], qda.covariance_[1])
plt.title(f"Trois classes JM_12={jm:3.2f}", fontsize="small")
plt.scatter(X1[:, 0], X1[:, 1], marker="o", c=Y1, s=25, edgecolor="k", cmap=cmap_classes)

plt.subplot(326)

X1, Y1 = make_gaussian_quantiles(n_features=2, n_classes=3)
qda = QuadraticDiscriminantAnalysis(store_covariance=True)
qda.fit(X1, Y1)
jm=jeffries_matusita_distance(qda.means_[0], qda.covariance_[0], qda.means_[1], qda.covariance_[1])
plt.title(f"Trois classes, Gaussiennes superposées JM_12={jm:3.2f}", fontsize="small")
plt.scatter(X1[:, 0], X1[:, 1], marker="o", c=Y1, s=25, edgecolor="k", cmap=cmap_classes)

plt.show()
```

On forme notre ensemble d'entrainement à partir du fichier `csv` de la section @sec-05.02.02. 

```{python}
df= pd.read_csv('sampling_points.csv')
# Extraire la colonne 'value'.
# 'value' est une chaîne de caractères représentation d'une liste de nombres.
# Nous devons la convertir en données numériques réelles.
X = df['value'].apply(lambda x: np.fromstring(x[1:-1], dtype=float, sep=' ')).to_list()

# on obtient une liste de numpy array  qu'il faut convertir en un numpy array 2D
X= np.array([row.tolist() for row in X])
idx= X.sum(axis=-1)>0 # on exclut certains points sans valeurs
X= X[idx,...]
y = df['class'].to_numpy()
y= y[idx]
class_labels = np.unique(y).tolist() # on cherche à savoir combien de classes uniques
n_classes = len(class_labels)
if max(class_labels) > n_classes: # il se peut que certaines classes soit absentes
  y_new= []
  for i,l in enumerate(class_labels):
    y_new.extend([i]*sum(y==l))
  y_new = np.array(y_new)

couleurs_classes2= [couleurs_classes[c] for c in np.unique(y).tolist()] # couleurs des classes
cmap_classes2 = ListedColormap(couleurs_classes2)
```

On peut faire une analyse de séparabilité sur notre ensemble d'entrainement de 10 classes. On obtient un tableau symmétrique de 10x10 valeurs. On peut observer des valeurs inférieures à 1 ce qui indique des séparabilités faibles entre ces classes sous l'hypothèse du modèle Gaussien:

```{python}
qda= QuadraticDiscriminantAnalysis(store_covariance=True)
qda.fit(X, y_new) # calcul des paramètres des distributions Gaussiennes
JM= []
classes= np.unique(y_new).tolist() # étiquettes uniques des classes
for cl1 in classes:
  for cl2 in classes:
    JM.append(jeffries_matusita_distance(qda.means_[cl1], qda.covariance_[cl1], qda.means_[cl2], qda.covariance_[cl2]))

JM= np.array(JM).reshape(len(classes),len(classes))
JM= pd.DataFrame(JM, index=classes, columns=classes)
JM.head(10)
```

## Méthodes non paramétriques

Les méthodes non paramétriques ne font pas d'hypothèses particulières sur les données. Un des inconvénients de ces modèles est que le nombre de paramètres du modèles augmente avec la taille des données.


### Méthode des parallélépipèdes {#sec-0511}

La méthode du parallélépipède est probablement la plus simple et consiste à délimiter directement le domaine des points d'une classe par une boite (un parallélépipède) à $D$ dimensions. Les limites de ces parallélépipèdes forment alors des frontières de décision manuelles qui vont permettre décider de la classe d'appartenance d'un nouveau point. Un des avantages de cette technique est que si un point n'est dans aucun parallélépipède alors on peut le laisser comme non classifié. Par contre, la construction de ces parallélépipèdes se complexifient grandement avec le nombre de bandes. À une dimension, deux paramètres, équivalents à un seuillage d'histogramme, sont suffisants. À deux dimensions, vous devez définir 4 segments par classe. Avec 3 bandes, vous devez définir 6 plans par classes et à D dimensions, D hyperplans à D-1 dimensions par classe. Le modèle ici est donc une suite de valeurs `min` et `max` pour chacune des bandes et des classes: 

```{python}
def parrallepiped_train(X_train, y_train):
  classes= np.unique(y_train).tolist()
  clf= []
  for cl in classes:
      data_cl= X_train[y_train == cl,...] # on cherche les données pour la classe courante
      
      limits=[]
      for b in range(data_cl.shape[1]):
        limits.append([data_cl[:,b].min(), data_cl[:,b].max()]) # on calcul le min et max pour chaque bande
      clf.append(np.array(limits))
  return clf
clf= parrallepiped_train(X, y_new)
```
La prédiction consiste à trouver pour chaque point la première limite qui est satisfaite. Notez qu'il n'y a pas de moyen de décider quelle est la meilleure classe si le point appartient à plusieurs classes.

```{python}
@jit(nopython=True)
def parrallepiped_predict(clf, X_test):
  y_pred= []
  for data in X_test:
    y_pred.append(np.nan)
    for cl, limits in enumerate(clf):
      inside= True
      for b,limit in enumerate(limits):
        inside = inside and (data[b] >= limit[0]) & (data[b] <= limit[1])
        if ~inside:
          break
      if inside:
        y_pred[-1]=cl
  return np.array(y_pred)
```
On peut appliquer ensuite le modèle sur l'image au complet. Les résultats sont assez mauvais, seule la classe eau en bleu semble être bien classifiée.
```{python}
data_image= img_rgbnir.to_numpy().transpose(1,2,0).reshape(img_rgbnir.shape[1]*img_rgbnir.shape[2],4)
y_image= parrallepiped_predict(clf, data_image)
y_image= y_image.reshape(img_rgbnir.shape[1],img_rgbnir.shape[2])

fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))
plt.imshow(y_image, cmap=cmap_classes2)
ax.set_title("Méthode des parrallélépipèdes", fontsize="small")
plt.show()
```

#### La malédiction de la haute dimension

Augmenter le nombre de dimension ou de caractéristiques des données permet de résoudre des problèmes complexes comme la classification d'image. Cependant, cela amène beaucoup de contraintes sur le volume des données. Supposons que nous avons N points occupant un segment linéaire de taille d. La densité de points est $N/d$. Si nous augmentons le nombre de dimension D, la densité de points va diminuer exponentiellement en $1/d^D$. Par conséquent, pour garder une densité constante et donc une bonne estimation des parallélépipèdes, il nous faudrait augmenter le nombre de points en puissance de D. Ceci porte le nom de la malédiction de la dimensionnalité (*dimensionality curse*). En résumé, l'espace vide augmente plus rapidement que le nombre de données d'entraînement et l'espace des données devient de plus en plus parcimonieux (*sparse*). Pour contrecarrer ce problème, on peut sélectionner les meilleures caractéristiques ou appliquer une réduction de dimension. 

### Plus proches voisins

La méthode des plus proches voisins (*K-Nearest-Neighbors* en Anglais) est certainement la plus simple des méthodes pour classifier des données. Elle consiste à comparer une nouvelle données avec ces voisins les plus proches en fonction d'une simple distance Euclidienne. Si une majorité de ces $K$ voisins appartiennent à une classe majoritaire alors cette classe est sélectionnée. Afin de permettre un vote majoritaire, on choisira un nombre impair pour la valeur de $K$. Mallgré sa simplicité, cette technique peut devenir assez demandante en terme de calcul pour un nombre important de points avec un nombre élevé de dimensions.

Reprensons l'ensemble d’entraînement formé à partir de notre image RGBNIR précédente:

```{python}
df= pd.read_csv('sampling_points.csv')
# Extraire la colonne 'value'.
# 'value' est une chaîne de caractères comme représentation d'une liste de valeurs.
# Nous devons la convertir en données numériques réelles.
X = df['value'].apply(lambda x: np.fromstring(x[1:-1], dtype=float, sep=' ')).to_list()

# on obtient une liste de numpy array  qu'il faut convertir en un numpy array 2D
X= np.array([row.tolist() for row in X])
idx= X.sum(axis=-1)>0 # il se peut qu'il y ait des valeurs erronées
X= X[idx,...]
y = df['class'].to_numpy()
y= y[idx]
class_labels = np.unique(y).tolist() # on cherche à savoir combien de classes uniques
n_classes = len(class_labels)
if max(class_labels) > n_classes: # il se peut que certaines classes soit absentes
  y_new= []
  for i,l in enumerate(class_labels):
    y_new.extend([i]*sum(y==l))
  y_new = np.array(y_new)
```

Il est important de faire précéder la méthode K-NN par une normalisation des données de façon à ce que chaque caractéristique soit de moyenne 0 et d'écart-type égale à 1 (on dit parfois que l'on blanchit les données). Cette normalisation permet à ce que chaque dimension ait le même poids dans le calcul des distances entre points. Cette opération porte le nom de `StandardScaler` dans `scikit-learn`. On peut alors former un pipeline de traitement combinant les deux opérations:

```{python}
clf = Pipeline(
    steps=[("scaler", StandardScaler()), ("knn", KNeighborsClassifier(n_neighbors=1))]
)
```

Avant d'effectuer un entraînement, on met généralement une portion des données pour valider les performances:

```{python}
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=0)
```
On peut visualiser les frontières de décision du K-NN pour différentes valeurs de $K$ lorsque seulement deux bandes sont utilisées (Rouge et proche infra-rouge ici):

```{python}
#| echo: false
#| warning: false

_, axs = plt.subplots(nrows=2,ncols=2, figsize=(9,9))
couleurs_classes2= [couleurs_classes[c] for c in np.unique(y).tolist()]
cmap_classes2 = ListedColormap(couleurs_classes2)
for ax, K in zip(axs.flatten(), [1,3,7,15]):
    clf.set_params(knn__weights='distance', knn__n_neighbors = K).fit(X_train[:,2:4], y_train)
    y_pred = clf.predict(X_test[:,2:4])
    print("Number of mislabeled points out of a total %d points : %d"
      % (X_test.shape[0], (y_test != y_pred).sum()))
    disp = DecisionBoundaryDisplay.from_estimator(
        clf,
        X[:,2:4],
        response_method="predict",
        plot_method="contourf",
        shading="auto",
        alpha=0.5,
        ax=ax,
        cmap= cmap_classes2,
    )
    scatter = disp.ax_.scatter(X_test[:, 2], X_test[:, 3], c=y_test, edgecolors="k", cmap=cmap_classes2)
    disp.ax_.set_xlabel('Rouge', fontsize="small")
    disp.ax_.set_ylabel('Proche infra-rouge', fontsize="small")
    _ = disp.ax_.set_title(
        f"K={clf[-1].n_neighbors} erreur={(y_test != y_pred).sum()/X_test.shape[0]*100:3.1f}%", fontsize="small"
    )

plt.show()
```
On peut voir comment les différentes frontières de décision se forment dans l'espace des bandes Rouge-NIR. L'augmentation de K rend ces frontières plus complexes et le calcul plus long.

```{python}
clf.set_params(knn__weights='distance', knn__n_neighbors = 7).fit(X_train, y_train)
y_pred = clf.predict(X_test)
print("Nombre de points misclassifiés sur %d points : %d"
  % (X_test.shape[0], (y_test != y_pred).sum()))
```

L'application du modèle (la prédiction) peut se faire sur toute l'image en transposant l'image sous forme d'une matrice avec Largeur x Hauteur lignes et 4 colonnes:

```{python}
data_image= img_rgbnir.to_numpy().transpose(1,2,0).reshape(img_rgbnir.shape[1]*img_rgbnir.shape[2],4)
y_classe= clf.predict(data_image)
y_classe= y_classe.reshape(img_rgbnir.shape[1],img_rgbnir.shape[2])
```

```{python}
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))
plt.imshow(y_classe, cmap=cmap_classes2)
ax.set_title("Carte d'occupation des sols avec K-NN", fontsize="small")
plt.show()
```

### Méthodes par arbre de décision

## Méthodes paramétriques

Les méthodes paramétriques se basent sur des modélisations statistiques des données pour permettre une classification. Contraitement au méthodes non paramétriques, elles ont un nombre fixe de paramètres qui ne dépend pas de la taille du jeu de données. Par contre, des hypothèses sont faites sur le comportement statistique des données. La classification consiste alors à trouver la classe la plus vraisemblable dont le modèle statistique décrit le mieux les valeurs observées. L'ensemble d’entraînement permettra alors de calculer les paramètres de chaque Gaussienne pour chacune des classes d'intérêt. 

### Méthode Bayésienne naïve

La méthode Bayésienne naïve Gaussienne consiste faire des hypothèses simplificatrices sur les données, en particulier l'indépendance des données et des dimensions. Ceci permet un calcul plus simple.
```{python}
from sklearn.naive_bayes import GaussianNB
gnb = GaussianNB()
y_pred = gnb.fit(X_train, y_train).predict(X_test)
print("Nombre de points erronés sur %d points : %d"
      % (X_test.shape[0], (y_test != y_pred).sum()))
```

```{python}
#| echo: false
#| warning: false
_, axs = plt.subplots(nrows=1,ncols=1, figsize=(4, 4))
couleurs_classes2= [couleurs_classes[c] for c in np.unique(y).tolist()]
cmap_classes2 = ListedColormap(couleurs_classes2)

#clf.set_params(knn__weights='distance', knn__n_neighbors = K).fit(X_train[:,2:4], y_train)
gnb= gnb.fit(X_train[:,2:4], y_train)
disp = DecisionBoundaryDisplay.from_estimator(
    gnb,
    X[:,2:4],
    response_method="predict",
    plot_method="contourf",
    shading="auto",
    alpha=0.5,
    ax=axs,
    cmap= cmap_classes2,
)
scatter = disp.ax_.scatter(X_test[:, 2], X_test[:, 3], c=y_test, edgecolors="k", cmap=cmap_classes2)
disp.ax_.set_xlabel('Rouge', fontsize="small")
disp.ax_.set_ylabel('NIR', fontsize="small")
_ = disp.ax_.set_title(
    f"Bayésien naif", fontsize="small"
)

plt.show()
```
On peut observer que les frontières de décision sont beaucoup plus régulières que pour K-NN.


```{python}
gnb.fit(X_train, y_train)
y_pred = gnb.predict(X_test)
print("Nombre de points misclassifiés sur %d points : %d"
  % (X_test.shape[0], (y_test != y_pred).sum()))
```

De la même manière, la prédiction peut s'appliquer sur toute l'image:

```{python}
#| echo: false
#| warning: false
data_image= img_rgbnir.to_numpy().transpose(1,2,0).reshape(img_rgbnir.shape[1]*img_rgbnir.shape[2],4)
y_classe= gnb.predict(data_image)
y_classe= y_classe.reshape(img_rgbnir.shape[1],img_rgbnir.shape[2])
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))
plt.imshow(y_classe, cmap=cmap_classes2)
ax.set_title("Carte d'occupation des sols avec la méthode Bayésienne naive", fontsize="small")
plt.show()
```
### Analyse Discriminante Quadratique (ADQ)

L'analyse discriminante quadratique peut-être vue comme une généralisation de l'approche Bayésienne naive qui suppose des modèles Gaussiens indépendants pour chaque dimension et chaque point. Ici, on va considérer un modèle Gaussien multivarié.

```{python}
qda = QuadraticDiscriminantAnalysis(store_covariance=True)
qda.fit(X_train, y_train)
y_pred = qda.predict(X_test)
print("Nombre de points misclassifiés sur %d points : %d"
  % (X_test.shape[0], (y_test != y_pred).sum()))
```


Les Gaussiennes multivariées peuvent être visualiser sous forme d'éllipses décrivant le domaine des valeurs de chaque classe:
```{python}
#| echo: false
#| warning: false

import matplotlib as mpl
colors= [couleurs_classes[c] for c in np.unique(y).tolist()]
def make_ellipses(gmm, bands, ax):
    for n, color in enumerate(colors):
        
        covariances = gmm.covariance_[n][np.ix_(bands, bands)]
        
        v, w = np.linalg.eigh(covariances)
        u = w[0] / np.linalg.norm(w[0])
        angle = np.arctan2(u[1], u[0])
        angle = 180 * angle / np.pi  # convert to degrees
        v = 2.0 * np.sqrt(2.0) * np.sqrt(v)
        #print(v)
        ell = mpl.patches.Ellipse(
            gmm.means_[n, :2], v[0], v[1], angle=180 + angle, color=color
        )
        ell.set_clip_box(ax.bbox)
        ell.set_alpha(0.5)
        ax.add_artist(ell)
        ax.set_aspect("equal", "datalim")

plt.figure(figsize=(9, 3))
#plt.subplots_adjust(
#    bottom=0.01, top=0.95, hspace=0.15, wspace=0.05, left=0.01, right=0.99
#)
noms_bandes= ['Bleu','Vert','Rouge','NIR']
for index, b in enumerate([0,1,3]):
  h = plt.subplot(1, 3, index+1)
  bands= [2,b]
  make_ellipses(qda, bands, h)
  for n, color in enumerate(colors):
          data = X[y_new == n,...]
          plt.scatter(
              data[:, bands[0]], data[:, bands[1]], s=0.8, color=color
          )
  h.set_xlabel(noms_bandes[2], fontsize="small")
  h.set_ylabel(noms_bandes[b], fontsize="small")
```

De la même manière, la prédiction peut s'appliquer sur toute l'image:

```{python}
#| echo: false
#| warning: false
data_image= img_rgbnir.to_numpy().transpose(1,2,0).reshape(img_rgbnir.shape[1]*img_rgbnir.shape[2],4)
y_classe= qda.predict(data_image)
y_classe= y_classe.reshape(img_rgbnir.shape[1],img_rgbnir.shape[2])
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))
plt.imshow(y_classe, cmap=cmap_classes2)
ax.set_title("Carte d'occupation des sols avec la méthode ADQ", fontsize="small")
plt.show()
```

### Réseaux de neurones

Les réseaux de neurones artificiels (RNA) ont connu un essor très important depuis les années 2010 avec des approches dites profondes. Ces aspects seront surtout abordés dans le tome 2 consacré à l'intelligence artificielle. On abordera ici seulement le perceptron simple et le perceptron multi-couches (MLP).

Le perceptron est l'unité de base d'un RNA et consiste en N connections, une unité de calcul (le neurone) avec une fonction d'activation et une sortie. Le perceptron ne permet de construire que des frontières de décision linéaires.

Le perceptron multi-couches est un réseau dense (*fully connected*) avec des couches cachées entre la couche d'entrée et la couche de sortie. qui permet de construire des frontières de décision beaucoup plus complexes via une hiérarchie de frontières de décision.

Ces réseaux sont entraînés via des techniques itératives d'optimisation de type descente en gradient avec une correction des paramètres (les poids) à l'aide de la rétro-propagation de l'erreur. L'erreur est mesurée via une fonction de coût que l'on cherche à réduire.