models.mlpmultivariate.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f77bc74d-58a1-4c14-b477-f52d28f2a869",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| default_exp models.mlpmultivariate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15392f6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "12fa25a4",
   "metadata": {},
   "source": [
    "# MLPMultivariate\n",
    "> One of the simplest neural architectures are Multi Layer Perceptrons (`MLP`) composed of stacked Fully Connected Neural Networks trained with backpropagation. Each node in the architecture is capable of modeling non-linear relationships granted by their activation functions. Novel activations like Rectified Linear Units (`ReLU`) have greatly improved the ability to fit deeper networks overcoming gradient vanishing problems that were associated with `Sigmoid` and `TanH` activations. For the forecasting task the last layer is changed to follow a auto-regression problem. This version is multivariate, indicating that it will predict all time series of the forecasting problem jointly. <br><br>**References**<br>-[Rosenblatt, F. (1958). \"The perceptron: A probabilistic model for information storage and organization in the brain.\"](https://psycnet.apa.org/record/1959-09865-001)<br>-[Fukushima, K. (1975). \"Cognitron: A self-organizing multilayered neural network.\"](https://pascal-francis.inist.fr/vibad/index.php?action=getRecordDetail&idt=PASCAL7750396723)<br>-[Vinod Nair, Geoffrey E. Hinton (2010). \"Rectified Linear Units Improve Restricted Boltzmann Machines\"](https://www.cs.toronto.edu/~fritz/absps/reluICML.pdf)<br>"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "e6036ce9",
   "metadata": {},
   "source": [
    "![Figure 1. Three layer MLP with autorregresive inputs.](imgs_models/mlp.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2508f7a9-1433-4ad8-8f2f-0078c6ed6c3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "from fastcore.test import test_eq\n",
    "from nbdev.showdoc import show_doc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44065066-e72a-431f-938f-1528adef9fe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "\n",
    "from neuralforecast.losses.pytorch import MAE\n",
    "from neuralforecast.common._base_multivariate import BaseMultivariate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce70cd14-ecb1-4205-8511-fecbd26c8408",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "class MLPMultivariate(BaseMultivariate):\n",
    "    \"\"\" MLPMultivariate\n",
    "\n",
    "    Simple Multi Layer Perceptron architecture (MLP) for multivariate forecasting. \n",
    "    This deep neural network has constant units through its layers, each with\n",
    "    ReLU non-linearities, it is trained using ADAM stochastic gradient descent.\n",
    "    The network accepts static, historic and future exogenous data, flattens \n",
    "    the inputs and learns fully connected relationships against the target variables.\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `h`: int, forecast horizon.<br>\n",
    "    `input_size`: int, considered autorregresive inputs (lags), y=[1,2,3,4] input_size=2 -> lags=[1,2].<br>\n",
    "    `n_series`: int, number of time-series.<br>\n",
    "    `stat_exog_list`: str list, static exogenous columns.<br>\n",
    "    `hist_exog_list`: str list, historic exogenous columns.<br>\n",
    "    `futr_exog_list`: str list, future exogenous columns.<br>\n",
    "    `n_layers`: int, number of layers for the MLP.<br>\n",
    "    `hidden_size`: int, number of units for each layer of the MLP.<br>\n",
    "    `loss`: PyTorch module, instantiated train loss class from [losses collection](https://nixtla.github.io/neuralforecast/losses.pytorch.html).<br>\n",
    "    `valid_loss`: PyTorch module=`loss`, instantiated valid loss class from [losses collection](https://nixtla.github.io/neuralforecast/losses.pytorch.html).<br>\n",
    "    `max_steps`: int=1000, maximum number of training steps.<br>\n",
    "    `learning_rate`: float=1e-3, Learning rate between (0, 1).<br>\n",
    "    `num_lr_decays`: int=-1, Number of learning rate decays, evenly distributed across max_steps.<br>\n",
    "    `early_stop_patience_steps`: int=-1, Number of validation iterations before early stopping.<br>\n",
    "    `val_check_steps`: int=100, Number of training steps between every validation loss check.<br>\n",
    "    `batch_size`: int=32, number of different series in each batch.<br>\n",
    "    `step_size`: int=1, step size between each window of temporal data.<br>\n",
    "    `scaler_type`: str='identity', type of scaler for temporal inputs normalization see [temporal scalers](https://nixtla.github.io/neuralforecast/common.scalers.html).<br>\n",
    "    `random_seed`: int=1, random_seed for pytorch initializer and numpy generators.<br>\n",
    "    `num_workers_loader`: int=os.cpu_count(), workers to be used by `TimeSeriesDataLoader`.<br>\n",
    "    `drop_last_loader`: bool=False, if True `TimeSeriesDataLoader` drops last non-full batch.<br>\n",
    "    `alias`: str, optional,  Custom name of the model.<br>\n",
    "    `optimizer`: Subclass of 'torch.optim.Optimizer', optional, user specified optimizer instead of the default choice (Adam).<br>\n",
    "    `optimizer_kwargs`: dict, optional, list of parameters used by the user specified `optimizer`.<br>\n",
    "    `**trainer_kwargs`: int,  keyword trainer arguments inherited from [PyTorch Lighning's trainer](https://pytorch-lightning.readthedocs.io/en/stable/api/pytorch_lightning.trainer.trainer.Trainer.html?highlight=trainer).<br>    \n",
    "    \"\"\"\n",
    "    # Class attributes\n",
    "    SAMPLING_TYPE = 'multivariate'\n",
    "    \n",
    "    def __init__(self,\n",
    "                 h,\n",
    "                 input_size,\n",
    "                 n_series,\n",
    "                 futr_exog_list = None,\n",
    "                 hist_exog_list = None,\n",
    "                 stat_exog_list = None,\n",
    "                 num_layers = 2,\n",
    "                 hidden_size = 1024,\n",
    "                 loss = MAE(),\n",
    "                 valid_loss = None,\n",
    "                 max_steps: int = 1000,\n",
    "                 learning_rate: float = 1e-3,\n",
    "                 num_lr_decays: int = -1,\n",
    "                 early_stop_patience_steps: int =-1,\n",
    "                 val_check_steps: int = 100,\n",
    "                 batch_size: int = 32,\n",
    "                 step_size: int = 1,\n",
    "                 scaler_type: str = 'identity',\n",
    "                 random_seed: int = 1,\n",
    "                 num_workers_loader: int = 0,\n",
    "                 drop_last_loader: bool = False,\n",
    "                 optimizer = None,\n",
    "                 optimizer_kwargs = None,\n",
    "                 **trainer_kwargs):\n",
    "\n",
    "        # Inherit BaseMultivariate class\n",
    "        super(MLPMultivariate, self).__init__(h=h,\n",
    "                                  input_size=input_size,\n",
    "                                  n_series=n_series,\n",
    "                                  futr_exog_list=futr_exog_list,\n",
    "                                  hist_exog_list=hist_exog_list,\n",
    "                                  stat_exog_list=stat_exog_list,\n",
    "                                  loss=loss,\n",
    "                                  valid_loss=valid_loss,\n",
    "                                  max_steps=max_steps,\n",
    "                                  learning_rate=learning_rate,\n",
    "                                  num_lr_decays=num_lr_decays,\n",
    "                                  early_stop_patience_steps=early_stop_patience_steps,\n",
    "                                  val_check_steps=val_check_steps,\n",
    "                                  batch_size=batch_size,\n",
    "                                  step_size=step_size,\n",
    "                                  scaler_type=scaler_type,\n",
    "                                  num_workers_loader=num_workers_loader,\n",
    "                                  drop_last_loader=drop_last_loader,\n",
    "                                  random_seed=random_seed,\n",
    "                                  optimizer=optimizer,\n",
    "                                  optimizer_kwargs=optimizer_kwargs,\n",
    "                                  **trainer_kwargs)\n",
    "\n",
    "        # Architecture\n",
    "        self.num_layers = num_layers\n",
    "        self.hidden_size = hidden_size\n",
    "\n",
    "        self.futr_input_size = len(self.futr_exog_list)\n",
    "        self.hist_input_size = len(self.hist_exog_list)\n",
    "        self.stat_input_size = len(self.stat_exog_list)\n",
    "\n",
    "        input_size_first_layer = n_series * (input_size + self.hist_input_size * input_size + \\\n",
    "                                 self.futr_input_size*(input_size + h) + self.stat_input_size)\n",
    "\n",
    "        # MultiLayer Perceptron\n",
    "        layers = [nn.Linear(in_features=input_size_first_layer, out_features=hidden_size)]\n",
    "        for i in range(num_layers - 1):\n",
    "            layers += [nn.Linear(in_features=hidden_size, out_features=hidden_size)]\n",
    "        self.mlp = nn.ModuleList(layers)\n",
    "\n",
    "        # Adapter with Loss dependent dimensions\n",
    "        self.out = nn.Linear(in_features=hidden_size, \n",
    "                             out_features=h * self.loss.outputsize_multiplier * n_series)\n",
    "\n",
    "    def forward(self, windows_batch):\n",
    "\n",
    "        # Parse windows_batch\n",
    "        x             = windows_batch['insample_y']             #   [batch_size (B), input_size (L), n_series (N)]\n",
    "        hist_exog     = windows_batch['hist_exog']              #   [B, hist_exog_size (X), L, N]\n",
    "        futr_exog     = windows_batch['futr_exog']              #   [B, futr_exog_size (F), L + h, N]\n",
    "        stat_exog     = windows_batch['stat_exog']              #   [N, stat_exog_size (S)]\n",
    "\n",
    "        # Flatten MLP inputs [B, C, L+H, N] -> [B, C * (L+H) * N]\n",
    "        # Contatenate [ Y^1_t, ..., Y^N_t | X^1_{t-L},..., X^1_{t}, ..., X^N_{t} | F^1_{t-L},..., F^1_{t+H}, ...., F^N_{t+H} | S^1, ..., S^N ]\n",
    "        batch_size = x.shape[0]\n",
    "        x = x.reshape(batch_size, -1)\n",
    "        if self.hist_input_size > 0:\n",
    "            x = torch.cat(( x, hist_exog.reshape(batch_size, -1) ), dim=1)\n",
    "\n",
    "        if self.futr_input_size > 0:\n",
    "            x = torch.cat(( x, futr_exog.reshape(batch_size, -1) ), dim=1)\n",
    "\n",
    "        if self.stat_input_size > 0:\n",
    "            x = torch.cat(( x, stat_exog.reshape(batch_size, -1) ), dim=1)\n",
    "\n",
    "        for layer in self.mlp:\n",
    "             x = torch.relu(layer(x))\n",
    "        x = self.out(x)\n",
    "        \n",
    "        x = x.reshape(batch_size, self.h, -1)\n",
    "        forecast = self.loss.domain_map(x)\n",
    "\n",
    "        # domain_map might have squeezed the last dimension in case n_series == 1\n",
    "        # Note that this fails in case of a tuple loss, but Multivariate does not support tuple losses yet.\n",
    "        if forecast.ndim == 2:\n",
    "            return forecast.unsqueeze(-1)\n",
    "        else:\n",
    "            return forecast"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfc06a06",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(MLPMultivariate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a23696b",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(MLPMultivariate.fit, name='MLPMultivariate.fit')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8475d33",
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(MLPMultivariate.predict, name='MLPMultivariate.predict')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bf909e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "import logging\n",
    "import warnings\n",
    "\n",
    "from neuralforecast import NeuralForecast\n",
    "from neuralforecast.utils import AirPassengersPanel, AirPassengersStatic\n",
    "from neuralforecast.losses.pytorch import MAE, MSE, RMSE, MAPE, SMAPE, MASE, relMSE, QuantileLoss, MQLoss, DistributionLoss,PMM, GMM, NBMM, HuberLoss, TukeyLoss, HuberQLoss, HuberMQLoss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7ee8d15",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "# Test losses\n",
    "logging.getLogger(\"pytorch_lightning\").setLevel(logging.ERROR)\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "Y_train_df = AirPassengersPanel[AirPassengersPanel.ds<AirPassengersPanel['ds'].values[-12]].reset_index(drop=True) # 132 train\n",
    "Y_test_df = AirPassengersPanel[AirPassengersPanel.ds>=AirPassengersPanel['ds'].values[-12]].reset_index(drop=True) # 12 test\n",
    "\n",
    "AirPassengersStatic_single = AirPassengersStatic[AirPassengersStatic[\"unique_id\"] == 'Airline1']\n",
    "Y_train_df_single = Y_train_df[Y_train_df[\"unique_id\"] == 'Airline1']\n",
    "Y_test_df_single = Y_test_df[Y_test_df[\"unique_id\"] == 'Airline1']\n",
    "\n",
    "losses = [MAE(), MSE(), RMSE(), MAPE(), SMAPE(), MASE(seasonality=12), relMSE(y_train=Y_train_df), QuantileLoss(q=0.5), MQLoss(), DistributionLoss(distribution='Bernoulli'), DistributionLoss(distribution='Normal'), DistributionLoss(distribution='Poisson'), DistributionLoss(distribution='StudentT'), DistributionLoss(distribution='NegativeBinomial'), DistributionLoss(distribution='Tweedie'), PMM(), GMM(), NBMM(), HuberLoss(), TukeyLoss(), HuberQLoss(q=0.5), HuberMQLoss()]\n",
    "valid_losses = [MAE(), MSE(), RMSE(), MAPE(), SMAPE(), MASE(seasonality=12), relMSE(y_train=Y_train_df), QuantileLoss(q=0.5), MQLoss(), DistributionLoss(distribution='Bernoulli'), DistributionLoss(distribution='Normal'), DistributionLoss(distribution='Poisson'), DistributionLoss(distribution='StudentT'), DistributionLoss(distribution='NegativeBinomial'), DistributionLoss(distribution='Tweedie'), PMM(), GMM(), NBMM(), HuberLoss(), TukeyLoss(), HuberQLoss(q=0.5), HuberMQLoss()]\n",
    "\n",
    "for loss, valid_loss in zip(losses, valid_losses):\n",
    "    try:\n",
    "        model = MLPMultivariate(h=12, \n",
    "                        input_size=24,\n",
    "                        n_series=2,\n",
    "                        loss = loss,\n",
    "                        valid_loss = valid_loss,\n",
    "                        scaler_type='robust',\n",
    "                        learning_rate=1e-3,\n",
    "                        max_steps=2,\n",
    "                        val_check_steps=10,\n",
    "                        early_stop_patience_steps=2,\n",
    "                        )\n",
    "\n",
    "        fcst = NeuralForecast(models=[model], freq='M')\n",
    "        fcst.fit(df=Y_train_df, static_df=AirPassengersStatic, val_size=12)\n",
    "        forecasts = fcst.predict(futr_df=Y_test_df)\n",
    "    except Exception as e:\n",
    "        assert str(e) == f\"{loss} is not supported in a Multivariate model.\"\n",
    "\n",
    "\n",
    "# Test n_series = 1\n",
    "model = MLPMultivariate(h=12, \n",
    "                    input_size=24,\n",
    "                    n_series=1,\n",
    "                    loss = MAE(),\n",
    "                    scaler_type='robust',\n",
    "                    learning_rate=1e-3,\n",
    "                    max_steps=2,\n",
    "                    val_check_steps=10,\n",
    "                    early_stop_patience_steps=2,\n",
    "                )\n",
    "fcst = NeuralForecast(models=[model], freq='M')\n",
    "fcst.fit(df=Y_train_df_single, static_df=AirPassengersStatic_single, val_size=12)\n",
    "forecasts = fcst.predict(futr_df=Y_test_df_single)        "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c3e4e0f",
   "metadata": {},
   "source": [
    "## Usage Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72b60ba0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#| eval: false\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pytorch_lightning as pl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from neuralforecast import NeuralForecast\n",
    "# from neuralforecast.models import MLP\n",
    "from neuralforecast.losses.pytorch import MAE\n",
    "from neuralforecast.tsdataset import TimeSeriesDataset\n",
    "from neuralforecast.utils import AirPassengers, AirPassengersPanel, AirPassengersStatic\n",
    "\n",
    "Y_train_df = AirPassengersPanel[AirPassengersPanel.ds<AirPassengersPanel['ds'].values[-12]] # 132 train\n",
    "Y_test_df = AirPassengersPanel[AirPassengersPanel.ds>=AirPassengersPanel['ds'].values[-12]].reset_index(drop=True) # 12 test\n",
    "\n",
    "model = MLPMultivariate(h=12, \n",
    "            input_size=24,\n",
    "            n_series=2,\n",
    "            loss = MAE(),\n",
    "            scaler_type='robust',\n",
    "            learning_rate=1e-3,\n",
    "            max_steps=200,\n",
    "            val_check_steps=10,\n",
    "            early_stop_patience_steps=2)\n",
    "\n",
    "fcst = NeuralForecast(\n",
    "    models=[model],\n",
    "    freq='M'\n",
    ")\n",
    "fcst.fit(df=Y_train_df, static_df=AirPassengersStatic, val_size=12)\n",
    "forecasts = fcst.predict(futr_df=Y_test_df)\n",
    "\n",
    "Y_hat_df = forecasts.reset_index(drop=False).drop(columns=['unique_id','ds'])\n",
    "plot_df = pd.concat([Y_test_df, Y_hat_df], axis=1)\n",
    "plot_df = pd.concat([Y_train_df, plot_df])\n",
    "\n",
    "plot_df = plot_df[plot_df.unique_id=='Airline1'].drop('unique_id', axis=1)\n",
    "plt.plot(plot_df['ds'], plot_df['y'], c='black', label='True')\n",
    "plt.plot(plot_df['ds'], plot_df['MLPMultivariate'], c='blue', label='median')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.plot()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}