Tutoriel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction au Reservoir Computing et à ReservoirPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import reservoirpy as rpy\n",
    "from reservoirpy.nodes import Reservoir, Ridge\n",
    "\n",
    "rpy.verbosity(0)\n",
    "# Rendre les résultats déterministes\n",
    "rpy.set_seed(2341)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Un exemple simple: Mackey-Glass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import de la série Mackey-Glass\n",
    "from reservoirpy.datasets import mackey_glass\n",
    "\n",
    "timeseries = mackey_glass(n_timesteps=2_000, tau=17, seed=260_418)\n",
    "\n",
    "print(type(timeseries), timeseries.shape)\n",
    "\n",
    "plt.figure(figsize=(6, 2))\n",
    "plt.title(\"Série Mackey-Glass\")\n",
    "plt.plot(timeseries)\n",
    "plt.xlim(0, 600)\n",
    "plt.xlabel(\"Pas de temps\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Séparation X/Y et train/test\n",
    "from reservoirpy.datasets import to_forecasting\n",
    "\n",
    "X_train, X_test, Y_train, Y_test = to_forecasting(timeseries=timeseries, test_size=400, forecast=20)\n",
    "\n",
    "print(\n",
    "    X_train.shape,\n",
    "    X_test.shape,\n",
    "    Y_train.shape,\n",
    "    Y_test.shape,\n",
    ")\n",
    "cutoff = X_train.shape[0]\n",
    "\n",
    "\n",
    "plt.figure(figsize=(15, 5))\n",
    "plt.plot(np.arange(0, 1580), X_train, label=\"Données d'entraînement\")\n",
    "plt.plot(np.arange(0, 1580), Y_train, label=\"Objectif d'entraînement\")\n",
    "plt.plot(np.arange(1580, 1980), X_test, label=\"Données de test\")\n",
    "plt.plot(np.arange(1580, 1980), Y_test, label=\"Objectif de test\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Créer le modèle\n",
    "reservoir = Reservoir(\n",
    "    units=500, # Nombre de neurones\n",
    "    sr=1.25,   # Rayon spectral de la matrice de poids\n",
    "    lr=0.3,    # Taux de décharge\n",
    "    input_scaling=1. # Facteur multiplicatif de l'entrée\n",
    ")\n",
    "readout = Ridge(ridge=1e-4) # Régularisation: pénalise le sur-apprentissage\n",
    "model = reservoir >> readout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Entraîner le modèle sur les données d'entraînement\n",
    "model.fit(X_train, Y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Évaluer le modèle sur les données de test\n",
    "Y_pred = model.run(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Affichage des résultats\n",
    "plt.figure(figsize=(6, 3))\n",
    "plt.plot(Y_test, color=\"blue\", linewidth=2, label=\"Réel\")\n",
    "plt.plot(Y_pred, color=\"red\", label=\"Prédit\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Auscultons le réservoir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "activity = reservoir.run(X_test)\n",
    "\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.plot(activity[:200, :20]) # 200 pas de temps, 20 neurones\n",
    "plt.title(\"Activité du réservoir\")\n",
    "plt.xlabel(\"Pas de temps\")\n",
    "plt.ylabel(\"Activité des neurones\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prédiction univariée (Consommation électrique sur le réseau français)\n",
    "\n",
    "Lancé en 2011, [www.data.gouv.fr](www.data.gouv.fr) est un site gouvernemental français qui met à disposition plus de 45 000 jeux de données sous licence libre.\n",
    "\n",
    "RTE (Réseau de Transport de l'Électricité) propose au grand public des informations historiques et [en temps réel](https://www.rte-france.com/eco2mix/la-consommation-delectricite-en-france) sur la consommation électrique sur le réseau en France.\n",
    "\n",
    "Le jeu de données suivant, ainsi mis à disposition par RTE sur la plateforme [www.data.gouv.fr](www.data.gouv.fr), est une série temporelle qui présente la consommation électrique (en MW) du réseau par demi-heure, de 2012 à 2022.\n",
    "\n",
    "Avec près de 200 000 pas de temps, il peut être intéressant d'évaluer des tâches de prédiction à plusieurs échelles de temps (en rééchantillonnant, par exemple via `np.reshape` puis `np.mean`).\n",
    "\n",
    "__Lien de téléchargement :__ https://www.data.gouv.fr/fr/datasets/consommation-quotidienne-brute/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "from matplotlib.dates import DateFormatter\n",
    "\n",
    "# chargement du dataset\n",
    "df_conso_rte = pd.read_csv('./datasets/conso_electricite_rte.csv')\n",
    "# conversion de la colonne date_heure au bon type\n",
    "df_conso_rte[\"date_heure\"] = pd.to_datetime(df_conso_rte[\"date_heure\"], utc=True)\n",
    "# Affichage\n",
    "plt.figure(figsize=(12, 4))\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.title(\"Évolution totale\")\n",
    "plt.plot(df_conso_rte[\"date_heure\"], df_conso_rte[\"consommation_brute_electricite_rte\"] / 1e3)\n",
    "plt.xlabel(\"Date\")\n",
    "plt.ylabel(\"Conso (GW)\")\n",
    "\n",
    "ax = plt.subplot(1,2,2)\n",
    "plt.title(\"Évolution hebdomadaire\")\n",
    "plt.plot(df_conso_rte[\"date_heure\"][:2*24*7], df_conso_rte[\"consommation_brute_electricite_rte\"][:2*24*7] / 1e3)\n",
    "plt.xlabel(\"Date\")\n",
    "plt.ylabel(\"Conso (GW)\")\n",
    "ax.xaxis.set_major_formatter(DateFormatter(\"%m-%d\"))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pré-traitement\n",
    "# Resampling journalier\n",
    "consommation_demi_heure = df_conso_rte[\"consommation_brute_electricite_rte\"].to_numpy()\n",
    "consommation_jour = consommation_demi_heure[:-3].reshape((-1, 2*24, 1)).mean(axis=1)\n",
    "\n",
    "# On vérifie que les dimensions du tableau sont bonnes\n",
    "print(f\"Dimensions des données: {consommation_jour.shape}\")\n",
    "\n",
    "# Normalisation\n",
    "conso_min, conso_max = np.min(consommation_jour), np.max(consommation_jour)\n",
    "consommation_jour = (consommation_jour - conso_min) / (conso_max - conso_min) - 0.5\n",
    "\n",
    "# Forecast et séparation train/test\n",
    "X_train, X_test, Y_train, Y_test = to_forecasting(timeseries=consommation_jour, forecast=5, test_size=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Création du modèle\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Entraînement du modèle\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test du modèle\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Évaluation des performances du modèle\n",
    "from reservoirpy.observables import nrmse\n",
    "print(nrmse(y_true=Y_test, y_pred=Y_pred))\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(Y_test, color=\"black\", linewidth=3)\n",
    "plt.plot(Y_pred, color=\"red\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Visualiser l'activité du réservoir\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exploration des hyper-paramètres\n",
    "\n",
    "Vous pouvez consulter le [tutoriel ReservoirPy sur la recherche d'hyper-paramètres](https://github.com/reservoirpy/reservoirpy/blob/master/tutorials/4-Understand_and_optimize_hyperparameters.ipynb) pour plus de détails :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Définition de la fonction objectif\n",
    "from reservoirpy.observables import nrmse, rsquare\n",
    "\n",
    "def objective(dataset, config, *, input_scaling, N, sr, lr, ridge, seed):\n",
    "    x_train, y_train, x_test, y_test = dataset\n",
    "    seed = int(seed)\n",
    "    \n",
    "    # Construction du modèle à partir des hyper-paramètres explorés\n",
    "    reservoir = Reservoir(\n",
    "        units=N, \n",
    "        sr=sr, \n",
    "        lr=lr, \n",
    "        input_scaling=input_scaling, \n",
    "        seed=seed,\n",
    "    )\n",
    "    readout = Ridge(ridge=ridge)\n",
    "\n",
    "    model = reservoir >> readout\n",
    "\n",
    "    y_pred = model.fit(x_train, y_train) \\\n",
    "                        .run(x_test)\n",
    "\n",
    "    return {'loss': nrmse(y_test, y_pred, norm_value=np.ptp(x_train)),\n",
    "            'r2': rsquare(y_test, y_pred)}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "hyperopt_config = {\n",
    "    \"exp\": \"recherche_1\",    # nom de l'expérience\n",
    "    \"hp_max_evals\": 200,     # nombre d'itérations\n",
    "    \"hp_method\": \"random\",   # méthode d'exploration\n",
    "    \"seed\": 42,              # seed des paramètres sélectionnés\n",
    "    \"hp_space\": {            # définition de l'espace exploré\n",
    "        \"N\": [\"choice\", 500],             # réservoirs à 500 neurons\n",
    "        \"sr\": [\"loguniform\", 1e-2, 10],   # rayon spectral exploré log-uniformément entre 0.01 et 10\n",
    "        \"lr\": [\"loguniform\", 1e-3, 1],    # idem avec le taux de décharge, de 1e-3 à 1\n",
    "        \"input_scaling\": [\"choice\", 1.0], # on fixe l'input scaling à 1\n",
    "        \"ridge\": [\"loguniform\", 1e-8, 1e1], # régularisation, log-exploré entre 10^-8 à 1\n",
    "        \"seed\": [\"uniform\", 1, 10000]       # seed d'initialisation du réservoir\n",
    "    }\n",
    "}\n",
    "\n",
    "# we precautionously save the configuration in a JSON file\n",
    "# each file will begin with a number corresponding to the current experimentation run number.\n",
    "with open(f\"{hyperopt_config['exp']}.config.json\", \"w+\") as f:\n",
    "    json.dump(hyperopt_config, f)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from reservoirpy.hyper import research\n",
    "\n",
    "dataset = (X_train, Y_train, X_test, Y_test)\n",
    "best = research(objective, dataset, f\"{hyperopt_config['exp']}.config.json\", \".\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from reservoirpy.hyper import plot_hyperopt_report\n",
    "fig = plot_hyperopt_report(hyperopt_config[\"exp\"], (\"lr\", \"sr\", \"ridge\"), metric=\"r2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification de séries temporelles: FordA\n",
    "\n",
    "On prend ici un jeu de données de classification :\n",
    "\n",
    "Le problème de classification consiste à diagnostiquer si un certain symptôme existe ou non dans un sous-système automobile. Chaque cas se compose de 500 mesures du bruit du moteur et d'une classification.\n",
    "http://timeseriesclassification.com/description.php?Dataset=FordA\n",
    "\n",
    "[Guide de la classification avec ReservoirPy](https://github.com/reservoirpy/reservoirpy/blob/master/tutorials/5-Classification-with-RC.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "with open(\"./datasets/FordA.p\", 'rb') as file:\n",
    "    fordA = pickle.load(file)\n",
    "\n",
    "print(\n",
    "    type(fordA),\n",
    "    fordA[0].shape, # X_train\n",
    "    fordA[1].shape, # Y_train\n",
    "    fordA[2].shape, # X_test\n",
    "    fordA[3].shape, # Y_test\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pré-traitement\n",
    "\n",
    "Pour ReservoirPy, toutes les données sont des séries temporelles. Ce qui n'est pas le cas ici pour les données à prédire en classification. Réglons ça."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train = fordA[0]\n",
    "Y_train = fordA[1]\n",
    "X_test = fordA[2]\n",
    "Y_test = fordA[3]\n",
    "\n",
    "Y_train = np.repeat(Y_train.reshape(-1, 1, 1), 500, axis=1)\n",
    "Y_test = np.repeat(Y_test.reshape(-1, 1, 1), 500, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Création d'un modèle et entraînement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Création du modèle\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Entraîner le modèle à classifier\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lancer le modèle sur l'ensemble de test\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Décoder les prédictions\n",
    "# Il faut avoir 1 sortie par série !\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Affichage des résultats\n",
    "from sklearn.metrics import classification_report, ConfusionMatrixDisplay\n",
    "\n",
    "print(classification_report(y_true=Y_test, y_pred=Y_pred))\n",
    "\n",
    "_ = ConfusionMatrixDisplay.from_predictions(y_true=Y_test, y_pred=Y_pred)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exploration des hyper-paramètres"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Définition de la fonction objectif\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configuration de la recherche\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Lancement de la recherche\n",
    "#TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Affichage de l'exploration\n",
    "#TODO"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "py311",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.1.-1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}