{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "eba00088",
   "metadata": {},
   "source": [
    "**Work in progress**\n",
    "\n",
    "Contenu rédigé à partir des diapositives coursTransformers.pdf, faute de TP\n",
    "préexistant pour cette séance. Susceptible d'être ajusté.\n",
    "\n",
    "# XIII.2 - Travaux pratiques - Prédiction de séries temporelles avec un Transformer\n",
    "\n",
    "La séance précédente, nous avons entraîné un réseau récurrent (GRU) pour prédire la température future à partir des paramètres météorologiques passés. L'objectif de cette séance est de reprendre exactement le même problème, mais en remplaçant le réseau récurrent par un **Transformer** : nous implémenterons d'abord un mécanisme d'auto-attention « à la main », puis nous nous en servirons pour construire un petit encodeur Transformer avec PyTorch, que nous entraînerons et comparerons au GRU de la semaine dernière.\n",
    "\n",
    "Les Transformers ont été introduits pour le traitement du langage (traduction automatique), mais l'architecture ne fait aucune hypothèse sur la nature des données en entrée : tout ce dont elle a besoin est une **séquence de vecteurs**. Une séquence de mesures météorologiques au cours du temps convient donc tout aussi bien qu'une séquence de mots — c'est cette généralité que nous allons exploiter ici, uniquement pour des séries temporelles (nous ne ferons pas de traitement du langage dans ce TP).\n",
    "\n",
    "## Rappel : préparation des données\n",
    "\n",
    "Nous reprenons le jeu de données météorologique et la mise en forme *many-to-one* de la fin de la séance précédente (prédiction de la température à l'instant $T=6$ à partir des $M=720$ observations passées, sous-échantillonnées d'un facteur $6$, soit des séquences de longueur $120$). Pour le détail de ces choix, se reporter au TP précédent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6550f3b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import urllib.request\n",
    "\n",
    "if not os.path.exists('meteo_dataset.csv'):\n",
    "    urllib.request.urlretrieve('https://cedric.cnam.fr/vertigo/Cours/ml2/docs/meteo_dataset.csv', 'meteo_dataset.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8f4456c",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import torch\n",
    "from torch import nn\n",
    "from torch.utils.data import Dataset, DataLoader\n",
    "\n",
    "data = pd.read_csv('meteo_dataset.csv')\n",
    "\n",
    "split_fraction = 0.80\n",
    "split_idx = int(split_fraction * len(data))\n",
    "\n",
    "train_data = data.loc[0:split_idx-1].values\n",
    "val_data = data.loc[split_idx:].values\n",
    "\n",
    "past_length = 720\n",
    "future_step = 6\n",
    "step_size = 6\n",
    "sequence_length = int(past_length / step_size)\n",
    "batch_size = 256\n",
    "\n",
    "class FutureTimeseriesDataset(Dataset):\n",
    "    \"\"\" Fenêtres glissantes de longueur `sequence_length` (échantillonnées tous les\n",
    "        `sampling_rate` pas de temps), chacune associée à une unique étiquette (celle de\n",
    "        `targets` à `future_step` pas de temps après la fin de la fenêtre). Voir le TP précédent. \"\"\"\n",
    "    def __init__(self, data, targets, sequence_length, sampling_rate, future_step):\n",
    "        self.data = torch.as_tensor(data, dtype=torch.float32)\n",
    "        self.targets = torch.as_tensor(targets, dtype=torch.float32)\n",
    "        self.sequence_length = sequence_length\n",
    "        self.sampling_rate = sampling_rate\n",
    "        self.window_span = sequence_length * sampling_rate\n",
    "        self.future_step = future_step\n",
    "\n",
    "    def __len__(self):\n",
    "        return len(self.data) - self.window_span - self.future_step\n",
    "\n",
    "    def __getitem__(self, idx):\n",
    "        x = self.data[idx: idx + self.window_span: self.sampling_rate]\n",
    "        y = self.targets[idx + self.window_span + self.future_step]\n",
    "        return x, y\n",
    "\n",
    "dataset_train = FutureTimeseriesDataset(\n",
    "    data=train_data, targets=train_data[:, 1],\n",
    "    sequence_length=sequence_length, sampling_rate=step_size, future_step=future_step,\n",
    ")\n",
    "train_loader = DataLoader(dataset_train, batch_size=batch_size, shuffle=True)\n",
    "\n",
    "dataset_val = FutureTimeseriesDataset(\n",
    "    data=val_data, targets=val_data[:, 1],\n",
    "    sequence_length=sequence_length, sampling_rate=step_size, future_step=future_step,\n",
    ")\n",
    "val_loader = DataLoader(dataset_val, batch_size=batch_size, shuffle=False)\n",
    "\n",
    "input_dim = train_data.shape[1]\n",
    "print(f\"Longueur des séquences : {sequence_length}, dimension des observations : {input_dim}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c5aaa5b",
   "metadata": {},
   "source": [
    "## Pourquoi un Transformer pour une série temporelle ?\n",
    "\n",
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Pourquoi ne peut-on pas paralléliser le calcul des états cachés d'un réseau récurrent le long d'une séquence ? En quoi l'architecture Transformer permet-elle de lever cette limite ? On pourra aussi discuter de la taille du contexte pris en compte par chaque pas de temps dans les deux architectures.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8eb74f6f",
   "metadata": {},
   "source": [
    "## Attention : implémentation « à la main »\n",
    "\n",
    "Comme pour la régression logistique et le perceptron des séances précédentes, nous allons commencer par implémenter nous-mêmes le calcul au cœur du Transformer — l'attention — avant de nous appuyer sur PyTorch pour la suite.\n",
    "\n",
    "Pour chaque pas de temps $i$ d'une séquence, on définit trois vecteurs obtenus par projection linéaire de l'observation $x_i$ :\n",
    "\n",
    "- une **query** $q_i = x_i W^Q$ : ce que ce pas de temps cherche comme information,\n",
    "- une **key** $k_i = x_i W^K$ : ce que ce pas de temps a à offrir,\n",
    "- une **value** $v_i = x_i W^V$ : le contenu qui sera effectivement agrégé.\n",
    "\n",
    "En pratique, on empile les $q_i$, $k_i$, $v_i$ de toute la séquence dans des matrices $Q$, $K$, $V$ (une ligne par pas de temps), ce qui permet de calculer l'attention de tous les pas de temps simultanément :\n",
    "\n",
    "$$\n",
    "\\text{Attention}(Q, K, V) = \\text{softmax}\\left(\\frac{QK^T}{\\sqrt{d_k}}\\right) V\n",
    "$$\n",
    "\n",
    "où $d_k$ est la dimension des vecteurs *key* (la division par $\\sqrt{d_k}$ évite que les scores ne deviennent trop grands en valeur absolue, ce qui saturerait le softmax).\n",
    "\n",
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Que représente la matrice $QK^T$ avant l'application du softmax ? Quelles sont ses dimensions si $Q$ et $K$ sont de dimensions $(L, d_k)$, $L$ étant la longueur de la séquence ? Qu'obtient-on après avoir appliqué le softmax **le long de chaque ligne** ?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "817b10c6",
   "metadata": {},
   "source": [
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Compléter la fonction `scaled_dot_product_attention` ci-dessous, qui implémente l'équation (attention). On rappelle que le softmax d'une matrice, appliqué le long de l'axe -1 (les colonnes), s'écrit `torch.exp(x) / torch.sum(torch.exp(x), dim=-1, keepdim=True)`.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e0601032",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "def scaled_dot_product_attention(Q, K, V):\n",
    "    \"\"\" Entrées :\n",
    "        - Q : matrice des queries, (L, d_k)\n",
    "        - K : matrice des keys, (L, d_k)\n",
    "        - V : matrice des values, (L, d_v)\n",
    "\n",
    "        Renvoie :\n",
    "        - Z : sortie de l'attention, (L, d_v)\n",
    "        - alpha : poids d'attention, (L, L)\n",
    "    \"\"\"\n",
    "    # À compléter\n",
    "    ...\n",
    "    return Z, alpha"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b1bf6c6",
   "metadata": {},
   "source": [
    "Vérifions notre implémentation sur une vraie fenêtre de notre jeu de données météorologique. Pour rester simple, on projette directement les observations (dimension 7) vers des queries, keys et values de petite dimension à l'aide de matrices aléatoires fixées (nous n'avons pas encore appris $W^Q$, $W^K$, $W^V$ par descente de gradient, ce qui sera le rôle de PyTorch dans la suite) :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e515e175",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_train = train_data  # alias, comme dans le TP précédent\n",
    "\n",
    "torch.manual_seed(0)\n",
    "d_qkv = 8\n",
    "\n",
    "window = torch.as_tensor(x_train[:sequence_length], dtype=torch.float32)  # une fenêtre de 120 pas de temps, 7 variables\n",
    "# Standardisation (moyenne nulle, écart-type unitaire) : les grandeurs physiques\n",
    "# brutes (pression ~1000, densité ~1300...) ont des échelles très différentes et\n",
    "# feraient déborder l'exponentielle du softmax si on les projetait telles quelles.\n",
    "window = (window - window.mean(dim=0)) / window.std(dim=0)\n",
    "\n",
    "W_Q = torch.randn(input_dim, d_qkv)\n",
    "W_K = torch.randn(input_dim, d_qkv)\n",
    "W_V = torch.randn(input_dim, d_qkv)\n",
    "\n",
    "Q = window @ W_Q\n",
    "K = window @ W_K\n",
    "V = window @ W_V\n",
    "\n",
    "Z, alpha = scaled_dot_product_attention(Q, K, V)\n",
    "print(f\"Z : {Z.shape}, alpha : {alpha.shape}, somme de la première ligne de alpha : {alpha[0].sum():.3f}\")\n",
    "\n",
    "plt.figure(figsize=(6, 5))\n",
    "plt.imshow(alpha.numpy(), cmap='viridis')\n",
    "plt.colorbar(label=\"Poids d'attention\")\n",
    "plt.xlabel(\"Pas de temps (key)\")\n",
    "plt.ylabel(\"Pas de temps (query)\")\n",
    "plt.title(\"Matrice d'attention (poids aléatoires, non entraînés)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "497f4e8a",
   "metadata": {},
   "source": [
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Les poids $W^Q$, $W^K$, $W^V$ utilisés ci-dessus sont aléatoires, pas appris. À quoi vous attendez-vous à ce que ressemble la matrice d'attention obtenue ? Est-ce cohérent avec la figure obtenue ?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55233395",
   "metadata": {},
   "source": [
    "## Bloc Transformer avec PyTorch\n",
    "\n",
    "En pratique, on n'utilise pas une seule tête d'attention mais plusieurs en parallèle (*multi-head attention*) : chaque tête peut ainsi apprendre à se concentrer sur des aspects différents de la séquence. Réimplémenter le multi-head attention à la main serait long ; nous utiliserons la couche `nn.MultiheadAttention`, qui implémente exactement l'équation (attention) (avec plusieurs têtes, et les projections $W^Q$, $W^K$, $W^V$ apprises).\n",
    "\n",
    "Un bloc encodeur Transformer complet (voir le cours) combine cette attention avec une normalisation par couche (*layer normalization*), des connexions résiduelles, et un petit réseau *feed-forward* appliqué indépendamment à chaque pas de temps :\n",
    "\n",
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Compléter la méthode `forward` du bloc Transformer ci-dessous. On veillera à :\n",
    "- appliquer l'attention multi-tête `self.att` en *self-attention*, c'est-à-dire avec la même séquence `inputs` comme query, key et value (``nn.MultiheadAttention`` renvoie un couple ``(sortie, poids d'attention)``),\n",
    "- ajouter une connexion résiduelle entre `inputs` et la sortie de l'attention, puis normaliser (`self.layernorm1`),\n",
    "- appliquer le réseau *feed-forward* `self.ffn` (déjà défini), puis à nouveau une connexion résiduelle et une normalisation (`self.layernorm2`).\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "522b51a0",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "class TransformerBlock(nn.Module):\n",
    "    def __init__(self, d_model, num_heads, ff_dim):\n",
    "        super().__init__()\n",
    "        self.att = nn.MultiheadAttention(embed_dim=d_model, num_heads=num_heads, batch_first=True)\n",
    "        self.ffn = nn.Sequential(\n",
    "            nn.Linear(d_model, ff_dim),\n",
    "            nn.ReLU(),\n",
    "            nn.Linear(ff_dim, d_model),\n",
    "        )\n",
    "        self.layernorm1 = nn.LayerNorm(d_model, eps=1e-6)\n",
    "        self.layernorm2 = nn.LayerNorm(d_model, eps=1e-6)\n",
    "\n",
    "    def forward(self, inputs):\n",
    "        # À compléter\n",
    "        ...\n",
    "        return out2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6681a136",
   "metadata": {},
   "source": [
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Contrairement aux couches récurrentes ou convolutives, l'attention ne fait aucune différence entre un pas de temps et un autre : permuter l'ordre des pas de temps en entrée donnerait exactement la même sortie (à la permutation près). Pourtant, l'ordre temporel est évidemment important pour une série temporelle. Comment le Transformer original résout-il ce problème ? Quelle solution retient-on en pratique, d'après le cours ?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "237117a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "class PositionalEmbedding(nn.Module):\n",
    "    def __init__(self, sequence_length, d_model):\n",
    "        super().__init__()\n",
    "        self.position_embedding = nn.Embedding(num_embeddings=sequence_length, embedding_dim=d_model)\n",
    "\n",
    "    def forward(self, inputs):\n",
    "        positions = torch.arange(inputs.shape[1], device=inputs.device)\n",
    "        return inputs + self.position_embedding(positions)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3e36b48",
   "metadata": {},
   "source": [
    "## Construction et entraînement du modèle complet\n",
    "\n",
    "Nous pouvons maintenant assembler le modèle complet : une projection linéaire des observations vers la dimension $d_{\\text{model}}$ du Transformer, l'ajout de l'encodage positionnel, un bloc Transformer, puis une réduction de la séquence en un unique vecteur (par moyenne sur l'axe temporel, l'équivalent PyTorch de `GlobalAveragePooling1D`) suivie d'une couche de sortie pour la régression.\n",
    "\n",
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Compléter la définition du modèle ci-dessous (une classe héritant de `nn.Module`). On utilisera `d_model = 32`, `num_heads = 4` et `ff_dim = 32` pour le bloc Transformer, et on réduira la séquence de sortie du bloc en faisant la moyenne sur l'axe temporel (`.mean(dim=1)`) avant la couche de sortie. On créera aussi un optimiseur `Adam` (pas d'apprentissage $0{,}001$) et l'erreur quadratique moyenne (`nn.MSELoss`) comme fonction de coût, comme pour le GRU de la semaine dernière.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee0c96ae",
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "d_model = 32\n",
    "num_heads = 4\n",
    "ff_dim = 32\n",
    "learning_rate = 0.001\n",
    "\n",
    "class TransformerModel(nn.Module):\n",
    "    def __init__(self, input_dim, sequence_length, d_model, num_heads, ff_dim):\n",
    "        super().__init__()\n",
    "        # À compléter\n",
    "        ...\n",
    "\n",
    "    def forward(self, x):\n",
    "        # À compléter\n",
    "        ...\n",
    "\n",
    "model = TransformerModel(input_dim, sequence_length, d_model, num_heads, ff_dim)\n",
    "criterion = ...  # À compléter\n",
    "optimizer = ...  # À compléter\n",
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7ffe13d",
   "metadata": {},
   "source": [
    "Entraînons ce modèle, comme pour le GRU, pendant 5 époques :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8918a9ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs = 5\n",
    "\n",
    "for epoch in range(epochs):\n",
    "    model.train()\n",
    "    for x, y in train_loader:\n",
    "        y_pred = model(x)\n",
    "        loss = criterion(y_pred, y)\n",
    "        optimizer.zero_grad()\n",
    "        loss.backward()\n",
    "        optimizer.step()\n",
    "\n",
    "    model.eval()\n",
    "    with torch.no_grad():\n",
    "        val_losses = [criterion(model(x), y).item() for x, y in val_loader]\n",
    "    print(f\"Epoch {epoch}: loss (train) = {loss.item():.4f}, loss (val) = {sum(val_losses)/len(val_losses):.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cc3236f",
   "metadata": {},
   "source": [
    "## Évaluation et comparaison avec le GRU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "745b07f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_plot(sequence, prediction, ground_truth, delta=future_step/step_size, title=None):\n",
    "    plt.title(title)\n",
    "    plt.plot(sequence.flatten(), \".-\", label=\"Température passée\")\n",
    "    plt.plot(len(sequence) + delta, prediction, \"go\", label=\"Prédiction\")\n",
    "    plt.plot(len(sequence) + delta, ground_truth, \"rx\", label=\"Vraie température\")\n",
    "    plt.legend()\n",
    "    plt.xlim(-1, len(sequence) + delta + 3)\n",
    "    plt.xlabel(\"Pas de temps\")\n",
    "    plt.show()\n",
    "\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    x, y = next(iter(val_loader))\n",
    "    y_pred = model(x)\n",
    "    for i in range(5):\n",
    "        show_plot(x[i][:, 1].numpy(), y_pred[i].item(), y[i].item(), title=\"Prédiction de la température (Transformer)\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7949637",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean_absolute_error = nn.L1Loss()\n",
    "error = 0\n",
    "\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    for x, y in val_loader:\n",
    "        y_pred = model(x)\n",
    "        error += mean_absolute_error(y_pred, y).item()\n",
    "\n",
    "print(f\"Erreur absolue moyenne (Transformer) : {error/len(val_loader):.5f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a423a6c",
   "metadata": {},
   "source": [
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "Comparer l'erreur absolue moyenne obtenue ici à celle du GRU de la semaine dernière, sur la même tâche de prédiction. Le Transformer fait-il mieux ? Comparer aussi le nombre de paramètres des deux modèles (calculé avec `sum(p.numel() for p in model.parameters())`) ainsi que, qualitativement, le temps d'entraînement d'une époque. Ce résultat vous semble-t-il cohérent avec la taille du jeu de données utilisé ?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b3ecb9c",
   "metadata": {},
   "source": [
    "<div class=\"admonition-question admonition\">\n",
    "\n",
    "**Question**\n",
    "\n",
    "*À faire chez vous pour approfondir, si vous avez le temps.*\n",
    "\n",
    "Empiler plusieurs blocs Transformer (par exemple 2 ou 3) à la suite les uns des autres, comme le prévoit l'architecture originale. Les performances s'améliorent-elles ? Vous pouvez également essayer de remplacer le `GlobalAveragePooling1D` par la sortie du seul dernier pas de temps de la séquence (`transformer_out[:, -1, :]`) : quel effet cela a-t-il ?\n",
    "\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
