{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Package R \"vortex\"\n",
    "\n",
    "using Vectorial exOgenous variables to foRecast Time-sErieX.\n",
    "\n",
    "Ce package permet de prévoir des courbes de PM10 (par exemple), en se basant sur l'historique des valeurs mais aussi des variables exogènes (par exemple la météo).\n",
    "\n",
    "Fonctions principales :\n",
    "\n",
    " * <code>getData</code> : charge un jeu de données en mémoire\n",
    " * <code>getForecast</code> : prédit les lendemains aux indices demandés\n",
    "\n",
    "Diverses méthodes permettent ensuite d'analyser les performances : <code>getError</code>, <code>plotXYZ</code> : voir la section \"see also\" dans <code>?plotError</code>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Chargement de la librairie (après compilation, \"R CMD INSTALL ppmfun/\")\n",
    "library(vortex)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note : sur la base de nos dernières expériences, on considère que \n",
    "\n",
    " * on ne touche pas à la fenêtre obtenue par la fonction <code>optimize</code> ;\n",
    " * on oublie la méthode consistant à prédire forme et niveau de manière complètement déconnectée : il faut relier les deux.\n",
    "\n",
    "### Acquisition des données\n",
    "\n",
    "Compte-tenu de la nature hétérogène des données utilisées $-$ fonctionnelles pour les PM10, vectorielles pour les variables exogènes $-$, celles-ci sont encapsulées (comme des listes) dans un objet de type *Data*. En interne, la $i^{eme}$ cellule correspondant aux données disponibles au $i^{eme}$ jour à l'heure $H$ de prédiction choisie (1h00, 8h00 ou 14h00) : c'est-à-dire les valeurs des PM10 de $H-24h$ à $H-1h$, ainsi que les variables météo prédites pour la période de 1h à 0h du jour courant (sauf si on prédit à 0h : on prend alors les valeurs mesurées de la veille).\n",
    "\n",
    "Méthodes d'un objet de classe \"Data\" : elles prennent comme argument \"index\", qui est un index entier ; mais une fonction de conversion existe : <code>dateIndexToInteger</code>.\n",
    "\n",
    " * <code>getTime</code> : suite des date+heure\n",
    " * <code>getCenteredSerie</code> : série centrée\n",
    " * <code>getLevel</code> : niveau\n",
    " * <code>getSerie</code> : série *non* centrée\n",
    " * <code>getExoHat</code> : variables exogènes prévues\n",
    " * <code>getExoDm1</code> : variables exogènes mesurées la veille\n",
    "\n",
    "Exemple :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Voir ?getData pour les arguments\n",
    "data = getData(ts_data=\"data/pm10_mesures_H_loc.csv\", exo_data=\"data/meteo_extra_noNAs.csv\",\n",
    "  input_tz = \"Europe/Paris\", working_tz=\"Europe/Paris\", predict_at=\"07\")\n",
    "data$getLevel(10) #niveau du jour 10\n",
    "data$getExoHat(17) #météo prévue pour le jour 18"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prédiction\n",
    "\n",
    "Deux types de prévisions du prochain bloc de $24h$ sont à distinguer :\n",
    "\n",
    " * prévision de la forme (centrée) ;\n",
    " * prévision du saut d'une fin de série au début de la suivante.\n",
    "\n",
    "Il faut ainsi préciser à la fois une méthode de prévision de forme (\"Average\", \"Persistence\" et \"Neighbors\" sont implémentées), et une méthode de prédiction de saut (\"Zero\", \"Persistence\" ou \"Neighbors\"). On détaille surtout la méthode à voisins ci-après, les autres étant des approches naïves que l'on peut considérer comme des références à améliorer.\n",
    "\n",
    " 1. **Préparation des données** : fenêtrage si demandé (paramètre \"memory\"), recherche des paires de jours consécutifs sans valeurs manquantes.\n",
    " 2. **Optimisation des paramètres d'échelle** : via la fonction <code>optimize</code> minimisant la somme des 45 dernières erreurs jounalières RMSE, sur des jours similaires.\n",
    " 3. **Prédiction finale** : une fois le (ou les, si \"simtype\" vaut \"mix\") paramètre d'échelle $h$ déterminé, les similarités sont évaluées sur les variables exogènes et/ou endogènes, sous la forme $s(i,j) = \\mbox{exp}\\left(-\\frac{\\mbox{dist}^2(i,j)}{h^2}\\right)$. La formule indiquée plus haut dans le rapport est alors appliquée.\n",
    "\n",
    "Détail technique : la sortie de la méthode <code>getForecast</code> est un objet de type Forecast, encapsulant les séries prévues ainsi que tous les paramètres optimisés par la méthode \"Neighbors\". Fonctions disponibles (argument \"index\" comme pour les fonctions sur Data) :\n",
    "\n",
    " * <code>getSerie</code> : série prévue (sans les information de temps)\n",
    " * <code>getParams</code> : liste des paramètres (poids, fenêtre, ...)\n",
    " * <code>getIndexInData</code> : indice du jour où s'effectue la prévision relativement au jeu de données\n",
    "\n",
    "### Calcul des erreurs\n",
    "\n",
    "Pour chacun des instants à prévoir jusqu'à minuit du jour courant (ou avant : argument *horizon*), on calcule l'erreur moyenne sur tous les instants similaires du passé. Deux types d'erreurs sont considérées :\n",
    "\n",
    " * l'erreur \"abs\" égale à la valeur absolue moyenne entre la mesure et la prédiction ;\n",
    " * l'erreur \"MAPE\" égale à l'erreur absolue normalisée par la mesure.\n",
    "\n",
    "### Expériences numériques"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=9, repr.plot.height=3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_endo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"endo\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_exo = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"exo\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_mix = predictPM10(data, 2200, 2230, 0,0, \"Neighbors\", \"Neighbors\", simtype=\"mix\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p = list(p_endo, p_exo, p_mix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "yrange_MAPE = range(p_mix$errors$MAPE, p_endo$errors$MAPE, p_exo$errors$MAPE)\n",
    "yrange_abs = range(p_mix$errors$abs, p_endo$errors$abs, p_exo$errors$abs)\n",
    "yrange_RMSE = range(p_mix$errors$RMSE, p_endo$errors$RMSE, p_exo$errors$RMSE)\n",
    "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
    "\n",
    "par(mfrow=c(1,3))\n",
    "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
    "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
    "{\n",
    "  for (j in 1:3) #model (mix,endo,exo)\n",
    "  {\n",
    "    xlab = if (j>1) \"\" else \"Temps\"\n",
    "    ylab = if (j>1) \"\" else \"Erreur\"\n",
    "    main = if (j>1) \"\" else titles[i]\n",
    "    plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=main, xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
    "    if (j<3)\n",
    "        par(new=TRUE)\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ne tenir compte que des similarités sur les variables exogènes semble conduire à l'erreur la plus faible."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_nn = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_np = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Persistence\", sameSeaon=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_nz = predictPM10(data, 2200, 2230, 0, 0, \"Neighbors\", \"Zero\", sameSeaon=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_pp = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Persistence\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_pz = predictPM10(data, 2200, 2230, 0, 0, \"Persistence\", \"Zero\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p = list(p_nn, p_np, p_nz, p_pp, p_pz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "yrange_MAPE = range(p_nn$errors$MAPE, p_nz$errors$MAPE, p_np$errors$MAPE, p_pp$errors$MAPE, p_pz$errors$MAPE)\n",
    "yrange_abs = range(p_nn$errors$abs, p_nz$errors$abs, p_np$errors$abs, p_pp$errors$abs, p_pz$errors$abs)\n",
    "yrange_RMSE = range(p_nn$errors$RMSE, p_nz$errors$RMSE, p_np$errors$RMSE, p_pp$errors$RMSE, p_pz$errors$RMSE)\n",
    "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
    "\n",
    "par(mfrow=c(1,3))\n",
    "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
    "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
    "{\n",
    "  for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
    "  {\n",
    "    xlab = if (j>1) \"\" else \"Temps\"\n",
    "    ylab = if (j>1) \"\" else \"Erreur\"\n",
    "    main = if (j>1) \"\" else titles[i]\n",
    "    plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
    "    if (j<5)\n",
    "      par(new=TRUE)\n",
    "  }\n",
    "}\n",
    "        \n",
    "\n",
    "p = list(p_nn_epandage, p_nn_nonpollue, p_nn_chauffage)\n",
    "forecasts_2 = lapply(1:length(data), function(index) ( if (is.na(p[[2]]$forecasts[[index]][1])) rep(NA,24) else p[[2]]$forecasts[[index]]$pred ) )\n",
    "e1 = getErrors(data, forecasts_1, 17)\n",
    " \n",
    "e = list(e1,e2,e3)\n",
    "yrange_MAPE = range(e1$MAPE, e2$MAPE, e3$MAPE)\n",
    "yrange_abs = range(e1$abs, e2$abs, e3$abs)\n",
    "yrange_RMSE = range(e1$RMSE, e2$RMSE, e3$RMSE)\n",
    "ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
    "\n",
    "par(mfrow=c(1,3))\n",
    "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
    "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
    "{\n",
    "  for (j in 1:3) #model (nn,np,nz,pp,pz)\n",
    "  {\n",
    "    xlab = if (j>1) \"\" else \"Temps\"\n",
    "    ylab = if (j>1) \"\" else \"Erreur\"\n",
    "    main = if (j>1) \"\" else titles[i]\n",
    "    plot(e[[j]][[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
    "    if (j<3)\n",
    "      par(new=TRUE)\n",
    "  }\n",
    "}\n",
    "\n",
    "par(mfrow=c(1,2))\n",
    "#p[[i]]$forecasts[[index]]\n",
    "#futurs des blocs du passé pour le jour 2290 ::\n",
    "futurs = lapply(1:length(p[[1]]$forecasts[[2290]]$indices),\n",
    "       function(index) ( data[[ p[[1]]$forecasts[[2290]]$indices[index]+1 ]]$pm10 ) )\n",
    "#vrai futur (en rouge), vrai jour (en noir)\n",
    "r_min = min( sapply( 1:length(futurs), function(index) ( min(futurs[[index]] ) ) ) )\n",
    "r_max = max( sapply( 1:length(futurs), function(index) ( max(futurs[[index]] ) ) ) )\n",
    "for (i in 1:length(futurs))\n",
    "{\n",
    "    plot(futurs[[i]], col=1, ylim=c(r_min,r_max), type=\"l\")\n",
    "    if (i<length(futurs))\n",
    "        par(new=TRUE)\n",
    "}\n",
    "\n",
    "plot(data[[2290]]$pm10, ylim=c(r_min, r_max), col=1, type=\"l\")\n",
    "    par(new=TRUE)\n",
    "plot(data[[2291]]$pm10, ylim=c(r_min, r_max), col=2, type=\"l\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Meilleurs results: nn et nz (np moins bon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#%%TODO: analyse sur les trois périodes indiquées par Michel ; simtype==\"exo\" par defaut\n",
    "#16/03/2015 2288\n",
    "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width=9, repr.plot.height=6)\n",
    "plot(p_nn_epandage$errors$abs, type=\"l\", col=1, main=\"Erreur absolue\", xlab=\"Temps\",\n",
    "     ylab=\"Erreur\", ylim=range(p_nn_epandage$errors$abs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#19/01/2015 2232\n",
    "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#23/02/2015 2267\n",
    "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"mix\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot(p_nn_nonpollue$errors$abs, type=\"l\", col=2, main=\"Erreur absolue\", xlab=\"Temps\",\n",
    "     ylab=\"Erreur\", ylim=range(p_nn_nonpollue$errors$abs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "library(ppmfun)\n",
    "data = getData(\"local\", \"7h\")\n",
    "p_nn_epandage = predictPM10(data, 2287, 2293, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")\n",
    "p_nn_nonpollue = predictPM10(data, 2266, 2272, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")\n",
    "p_nn_chauffage = predictPM10(data, 2231, 2237, 0, 0, \"Neighbors\", \"Neighbors\", sameSeaon=TRUE, simtype=\"endo\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p = list(p_nn_epandage, p_nn_nonpollue, p_nn_chauffage)\n",
    "#yrange_MAPE = range(p[[1]]$errors$MAPE, p[[2]]$errors$MAPE, p[[3]]$errors$MAPE)\n",
    "#yrange_abs = range(p[[1]]$errors$abs, p[[2]]$errors$abs, p[[3]]$errors$abs)\n",
    "#yrange_RMSE = range(p[[1]]$errors$RMSE, p[[2]]$errors$RMSE, p[[3]]$errors$RMSE)\n",
    "#ranges = list(yrange_MAPE,yrange_abs,yrange_RMSE)\n",
    "print(p[[1]]$forecasts[[2290]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "par(mfrow=c(1,3))\n",
    "titles = paste(\"Erreur\",c(\"MAPE\",\"abs\",\"RMSE\"))\n",
    "for (i in 1:3) #error type (MAPE,abs,RMSE)\n",
    "{\n",
    "  for (j in 1:5) #model (nn,np,nz,pp,pz)\n",
    "  {\n",
    "    xlab = if (j>1) \"\" else \"Temps\"\n",
    "    ylab = if (j>1) \"\" else \"Erreur\"\n",
    "    main = if (j>1) \"\" else titles[i]\n",
    "    plot(p[[j]]$errors[[i]], type=\"l\", col=j, main=titles[i], xlab=xlab, ylab=ylab, ylim=ranges[[i]])\n",
    "    if (j<5)\n",
    "      par(new=TRUE)\n",
    "  }\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bilan\n",
    "\n",
    "TODO"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}