losses.numpy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| default_exp losses.numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NumPy Evaluation\n",
    "\n",
    "> NeuralForecast contains a collection NumPy loss functions aimed to be used during the models' evaluation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most important train signal is the forecast error, which is the difference between the observed value $y_{\\tau}$ and the prediction $\\hat{y}_{\\tau}$, at time $y_{\\tau}$:\n",
    "\n",
    "$$e_{\\tau} = y_{\\tau}-\\hat{y}_{\\tau} \\qquad \\qquad \\tau \\in \\{t+1,\\dots,t+H \\}$$\n",
    "\n",
    "The train loss summarizes the forecast errors in different evaluation metrics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "from typing import Optional, Union\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "from IPython.display import Image\n",
    "from nbdev.showdoc import show_doc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "WIDTH = 600\n",
    "HEIGHT = 300"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def _divide_no_nan(a: np.ndarray, b: np.ndarray) -> np.ndarray:\n",
    "    \"\"\"\n",
    "    Auxiliary funtion to handle divide by 0\n",
    "    \"\"\"\n",
    "    div = a / b\n",
    "    div[div != div] = 0.0\n",
    "    div[div == float('inf')] = 0.0\n",
    "    return div"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def _metric_protections(y: np.ndarray, y_hat: np.ndarray, \n",
    "                        weights: Optional[np.ndarray]) -> None:\n",
    "    assert (weights is None) or (np.sum(weights) > 0), 'Sum of weights cannot be 0'\n",
    "    assert (weights is None) or (weights.shape == y.shape),\\\n",
    "        f'Wrong weight dimension weights.shape {weights.shape}, y.shape {y.shape}'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Scale-dependent Errors\n",
    "\n",
    "These metrics are on the same scale as the data."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean Absolute Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def mae(y: np.ndarray, y_hat: np.ndarray,\n",
    "        weights: Optional[np.ndarray] = None,\n",
    "        axis: Optional[int] = None) -> Union[float, np.ndarray]:\n",
    "    \"\"\"Mean Absolute Error\n",
    "\n",
    "    Calculates Mean Absolute Error between\n",
    "    `y` and `y_hat`. MAE measures the relative prediction\n",
    "    accuracy of a forecasting method by calculating the\n",
    "    deviation of the prediction and the true\n",
    "    value at a given time and averages these devations\n",
    "    over the length of the series.\n",
    "\n",
    "    $$ \\mathrm{MAE}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}_{\\\\tau}) = \\\\frac{1}{H} \\\\sum^{t+H}_{\\\\tau=t+1} |y_{\\\\tau} - \\hat{y}_{\\\\tau}| $$\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `y`: numpy array, Actual values.<br>\n",
    "    `y_hat`: numpy array, Predicted values.<br>\n",
    "    `mask`: numpy array, Specifies date stamps per serie to consider in loss.<br>\n",
    "\n",
    "    **Returns:**<br>\n",
    "    `mae`: numpy array, (single value).    \n",
    "    \"\"\"\n",
    "    _metric_protections(y, y_hat, weights)\n",
    "    \n",
    "    delta_y = np.abs(y - y_hat)\n",
    "    if weights is not None:\n",
    "        mae = np.average(delta_y[~np.isnan(delta_y)], \n",
    "                         weights=weights[~np.isnan(delta_y)],\n",
    "                         axis=axis)\n",
    "    else:\n",
    "        mae = np.nanmean(delta_y, axis=axis)\n",
    "        \n",
    "    return mae"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(mae, title_level=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](imgs_losses/mae_loss.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean Squared Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def mse(y: np.ndarray, y_hat: np.ndarray, \n",
    "        weights: Optional[np.ndarray] = None,\n",
    "        axis: Optional[int] = None) -> Union[float, np.ndarray]:\n",
    "    \"\"\"  Mean Squared Error\n",
    "\n",
    "    Calculates Mean Squared Error between\n",
    "    `y` and `y_hat`. MSE measures the relative prediction\n",
    "    accuracy of a forecasting method by calculating the \n",
    "    squared deviation of the prediction and the true\n",
    "    value at a given time, and averages these devations\n",
    "    over the length of the series.\n",
    "\n",
    "    $$ \\mathrm{MSE}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}_{\\\\tau}) = \\\\frac{1}{H} \\\\sum^{t+H}_{\\\\tau=t+1} (y_{\\\\tau} - \\hat{y}_{\\\\tau})^{2} $$\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `y`: numpy array, Actual values.<br>\n",
    "    `y_hat`: numpy array, Predicted values.<br>\n",
    "    `mask`: numpy array, Specifies date stamps per serie to consider in loss.<br>\n",
    "\n",
    "    **Returns:**<br>\n",
    "    `mse`: numpy array, (single value).\n",
    "    \"\"\"\n",
    "    _metric_protections(y, y_hat, weights)\n",
    "\n",
    "    delta_y = np.square(y - y_hat)\n",
    "    if weights is not None:\n",
    "        mse = np.average(delta_y[~np.isnan(delta_y)],\n",
    "                         weights=weights[~np.isnan(delta_y)],\n",
    "                         axis=axis)\n",
    "    else:\n",
    "        mse = np.nanmean(delta_y, axis=axis)\n",
    "\n",
    "    return mse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(mse, title_level=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](imgs_losses/mse_loss.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Root Mean Squared Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def rmse(y: np.ndarray, y_hat: np.ndarray,\n",
    "         weights: Optional[np.ndarray] = None,\n",
    "         axis: Optional[int] = None) -> Union[float, np.ndarray]:\n",
    "    \"\"\" Root Mean Squared Error\n",
    "\n",
    "    Calculates Root Mean Squared Error between\n",
    "    `y` and `y_hat`. RMSE measures the relative prediction\n",
    "    accuracy of a forecasting method by calculating the squared deviation\n",
    "    of the prediction and the observed value at a given time and\n",
    "    averages these devations over the length of the series.\n",
    "    Finally the RMSE will be in the same scale\n",
    "    as the original time series so its comparison with other\n",
    "    series is possible only if they share a common scale.\n",
    "    RMSE has a direct connection to the L2 norm.\n",
    "\n",
    "    $$ \\mathrm{RMSE}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}_{\\\\tau}) = \\\\sqrt{\\\\frac{1}{H} \\\\sum^{t+H}_{\\\\tau=t+1} (y_{\\\\tau} - \\hat{y}_{\\\\tau})^{2}} $$\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `y`: numpy array, Actual values.<br>\n",
    "    `y_hat`: numpy array, Predicted values.<br>\n",
    "    `mask`: numpy array, Specifies date stamps per serie to consider in loss.<br>\n",
    "\n",
    "    **Returns:**<br>\n",
    "    `rmse`: numpy array, (single value).\n",
    "    \"\"\"\n",
    "    return np.sqrt(mse(y, y_hat, weights, axis))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(rmse, title_level=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](imgs_losses/rmse_loss.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Percentage errors\n",
    "\n",
    "These metrics are unit-free, suitable for comparisons across series."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean Absolute Percentage Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def mape(y: np.ndarray, y_hat: np.ndarray, \n",
    "         weights: Optional[np.ndarray] = None,\n",
    "         axis: Optional[int] = None) -> Union[float, np.ndarray]:\n",
    "    \"\"\" Mean Absolute Percentage Error\n",
    "\n",
    "    Calculates Mean Absolute Percentage Error  between\n",
    "    `y` and `y_hat`. MAPE measures the relative prediction\n",
    "    accuracy of a forecasting method by calculating the percentual deviation\n",
    "    of the prediction and the observed value at a given time and\n",
    "    averages these devations over the length of the series.\n",
    "    The closer to zero an observed value is, the higher penalty MAPE loss\n",
    "    assigns to the corresponding error.\n",
    "\n",
    "    $$ \\mathrm{MAPE}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}_{\\\\tau}) = \\\\frac{1}{H} \\\\sum^{t+H}_{\\\\tau=t+1} \\\\frac{|y_{\\\\tau}-\\hat{y}_{\\\\tau}|}{|y_{\\\\tau}|} $$\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `y`: numpy array, Actual values.<br>\n",
    "    `y_hat`: numpy array, Predicted values.<br>\n",
    "    `mask`: numpy array, Specifies date stamps per serie to consider in loss.<br>\n",
    "\n",
    "    **Returns:**<br>\n",
    "    `mape`: numpy array, (single value).\n",
    "    \"\"\"\n",
    "    _metric_protections(y, y_hat, weights)\n",
    "        \n",
    "    delta_y = np.abs(y - y_hat)\n",
    "    scale = np.abs(y)\n",
    "    mape = _divide_no_nan(delta_y, scale)\n",
    "    mape = np.average(mape, weights=weights, axis=axis)\n",
    "    \n",
    "    return mape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(mape, title_level=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](imgs_losses/mape_loss.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## SMAPE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def smape(y: np.ndarray, y_hat: np.ndarray,\n",
    "          weights: Optional[np.ndarray] = None,\n",
    "          axis: Optional[int] = None) -> Union[float, np.ndarray]:\n",
    "    \"\"\" Symmetric Mean Absolute Percentage Error\n",
    "\n",
    "    Calculates Symmetric Mean Absolute Percentage Error between\n",
    "    `y` and `y_hat`. SMAPE measures the relative prediction\n",
    "    accuracy of a forecasting method by calculating the relative deviation\n",
    "    of the prediction and the observed value scaled by the sum of the\n",
    "    absolute values for the prediction and observed value at a\n",
    "    given time, then averages these devations over the length\n",
    "    of the series. This allows the SMAPE to have bounds between\n",
    "    0% and 200% which is desirable compared to normal MAPE that\n",
    "    may be undetermined when the target is zero.\n",
    "\n",
    "    $$ \\mathrm{sMAPE}_{2}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}_{\\\\tau}) = \\\\frac{1}{H} \\\\sum^{t+H}_{\\\\tau=t+1} \\\\frac{|y_{\\\\tau}-\\hat{y}_{\\\\tau}|}{|y_{\\\\tau}|+|\\hat{y}_{\\\\tau}|} $$\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `y`: numpy array, Actual values.<br>\n",
    "    `y_hat`: numpy array, Predicted values.<br>\n",
    "    `mask`: numpy array, Specifies date stamps per serie to consider in loss.<br>\n",
    "\n",
    "    **Returns:**<br>\n",
    "    `smape`: numpy array, (single value).\n",
    "    \n",
    "    **References:**<br>\n",
    "    [Makridakis S., \"Accuracy measures: theoretical and practical concerns\".](https://www.sciencedirect.com/science/article/pii/0169207093900793)\n",
    "    \"\"\"\n",
    "    _metric_protections(y, y_hat, weights)\n",
    "        \n",
    "    delta_y = np.abs(y - y_hat)\n",
    "    scale = np.abs(y) + np.abs(y_hat)\n",
    "    smape = _divide_no_nan(delta_y, scale)\n",
    "    smape = 2 * np.average(smape, weights=weights, axis=axis)\n",
    "    \n",
    "    if isinstance(smape, float):\n",
    "        assert smape <= 2, 'SMAPE should be lower than 200'\n",
    "    else:\n",
    "        assert all(smape <= 2), 'SMAPE should be lower than 200'\n",
    "    \n",
    "    return smape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(smape, title_level=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Scale-independent Errors\n",
    "\n",
    "These metrics measure the relative improvements versus baselines."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mean Absolute Scaled Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def mase(y: np.ndarray, y_hat: np.ndarray, \n",
    "         y_train: np.ndarray,\n",
    "         seasonality: int,\n",
    "         weights: Optional[np.ndarray] = None,\n",
    "         axis: Optional[int] = None) -> Union[float, np.ndarray]:\n",
    "    \"\"\" Mean Absolute Scaled Error \n",
    "    Calculates the Mean Absolute Scaled Error between\n",
    "    `y` and `y_hat`. MASE measures the relative prediction\n",
    "    accuracy of a forecasting method by comparinng the mean absolute errors\n",
    "    of the prediction and the observed value against the mean\n",
    "    absolute errors of the seasonal naive model.\n",
    "    The MASE partially composed the Overall Weighted Average (OWA), \n",
    "    used in the M4 Competition.\n",
    "\n",
    "    $$ \\mathrm{MASE}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}_{\\\\tau}, \\\\mathbf{\\hat{y}}^{season}_{\\\\tau}) = \\\\frac{1}{H} \\sum^{t+H}_{\\\\tau=t+1} \\\\frac{|y_{\\\\tau}-\\hat{y}_{\\\\tau}|}{\\mathrm{MAE}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}^{season}_{\\\\tau})} $$\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `y`: numpy array, (batch_size, output_size), Actual values.<br>\n",
    "    `y_hat`: numpy array, (batch_size, output_size)), Predicted values.<br>\n",
    "    `y_insample`: numpy array, (batch_size, input_size), Actual insample Seasonal Naive predictions.<br>\n",
    "    `seasonality`: int. Main frequency of the time series; Hourly 24,  Daily 7, Weekly 52, Monthly 12, Quarterly 4, Yearly 1.        \n",
    "    `mask`: numpy array, Specifies date stamps per serie to consider in loss.<br>\n",
    "\n",
    "    **Returns:**<br>\n",
    "    `mase`: numpy array, (single value).\n",
    "    \n",
    "    **References:**<br>\n",
    "    [Rob J. Hyndman, & Koehler, A. B. \"Another look at measures of forecast accuracy\".](https://www.sciencedirect.com/science/article/pii/S0169207006000239)<br>\n",
    "    [Spyros Makridakis, Evangelos Spiliotis, Vassilios Assimakopoulos, \"The M4 Competition: 100,000 time series and 61 forecasting methods\".](https://www.sciencedirect.com/science/article/pii/S0169207019301128)\n",
    "    \"\"\"\n",
    "    delta_y = np.abs(y - y_hat)\n",
    "    delta_y = np.average(delta_y, weights=weights, axis=axis)\n",
    "\n",
    "    scale = np.abs(y_train[:-seasonality] - y_train[seasonality:])\n",
    "    scale = np.average(scale, axis=axis)\n",
    "\n",
    "    mase = delta_y / scale\n",
    "\n",
    "    return mase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(mase, title_level=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](imgs_losses/mase_loss.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relative Mean Absolute Error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def rmae(y: np.ndarray, \n",
    "         y_hat1: np.ndarray, y_hat2: np.ndarray, \n",
    "         weights: Optional[np.ndarray] = None,\n",
    "         axis: Optional[int] = None) -> Union[float, np.ndarray]:\n",
    "    \"\"\" RMAE\n",
    "            \n",
    "    Calculates Relative Mean Absolute Error (RMAE) between\n",
    "    two sets of forecasts (from two different forecasting methods).\n",
    "    A number smaller than one implies that the forecast in the \n",
    "    numerator is better than the forecast in the denominator.\n",
    "    \n",
    "    $$ \\mathrm{rMAE}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}_{\\\\tau}, \\\\mathbf{\\hat{y}}^{base}_{\\\\tau}) = \\\\frac{1}{H} \\sum^{t+H}_{\\\\tau=t+1} \\\\frac{|y_{\\\\tau}-\\hat{y}_{\\\\tau}|}{\\mathrm{MAE}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}^{base}_{\\\\tau})} $$\n",
    "    \n",
    "    **Parameters:**<br>\n",
    "    `y`: numpy array, observed values.<br>\n",
    "    `y_hat1`: numpy array. Predicted values of first model.<br>\n",
    "    `y_hat2`: numpy array. Predicted values of baseline model.<br>\n",
    "    `weights`: numpy array, optional. Weights for weighted average.<br>\n",
    "    `axis`: None or int, optional.Axis or axes along which to average a.<br> \n",
    "        The default, axis=None, will average over all of the elements of\n",
    "        the input array.\n",
    "    \n",
    "    **Returns:**<br>\n",
    "    `rmae`: numpy array or double.\n",
    "\n",
    "    **References:**<br>\n",
    "    [Rob J. Hyndman, & Koehler, A. B. \"Another look at measures of forecast accuracy\".](https://www.sciencedirect.com/science/article/pii/S0169207006000239)\n",
    "    \"\"\"\n",
    "    numerator = mae(y=y, y_hat=y_hat1, weights=weights, axis=axis)\n",
    "    denominator = mae(y=y, y_hat=y_hat2, weights=weights, axis=axis)\n",
    "    rmae = numerator / denominator\n",
    "\n",
    "    return rmae"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(rmae, title_level=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](imgs_losses/rmae_loss.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Probabilistic Errors\n",
    "\n",
    "These measure absolute deviation non-symmetrically, that produce under/over estimation."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quantile Loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def quantile_loss(y: np.ndarray, y_hat: np.ndarray, q: float = 0.5, \n",
    "                  weights: Optional[np.ndarray] = None,\n",
    "                  axis: Optional[int] = None) -> Union[float, np.ndarray]:\n",
    "    \"\"\" Quantile Loss\n",
    "\n",
    "    Computes the quantile loss between `y` and `y_hat`.\n",
    "    QL measures the deviation of a quantile forecast.\n",
    "    By weighting the absolute deviation in a non symmetric way, the\n",
    "    loss pays more attention to under or over estimation.\n",
    "    A common value for q is 0.5 for the deviation from the median (Pinball loss).\n",
    "\n",
    "    $$ \\mathrm{QL}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}^{(q)}_{\\\\tau}) = \\\\frac{1}{H} \\\\sum^{t+H}_{\\\\tau=t+1} \\Big( (1-q)\\,( \\hat{y}^{(q)}_{\\\\tau} - y_{\\\\tau} )_{+} + q\\,( y_{\\\\tau} - \\hat{y}^{(q)}_{\\\\tau} )_{+} \\Big) $$\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `y`: numpy array, Actual values.<br>\n",
    "    `y_hat`: numpy array, Predicted values.<br>\n",
    "    `q`: float, between 0 and 1. The slope of the quantile loss, in the context of quantile regression, the q determines the conditional quantile level.<br>\n",
    "    `mask`: numpy array, Specifies date stamps per serie to consider in loss.<br>\n",
    "\n",
    "    **Returns:**<br>\n",
    "    `quantile_loss`: numpy array, (single value).\n",
    "    \n",
    "    **References:**<br>\n",
    "    [Roger Koenker and Gilbert Bassett, Jr., \"Regression Quantiles\".](https://www.jstor.org/stable/1913643)\n",
    "    \"\"\"\n",
    "    _metric_protections(y, y_hat, weights)\n",
    "\n",
    "    delta_y = y - y_hat\n",
    "    loss = np.maximum(q * delta_y, (q - 1) * delta_y)\n",
    "\n",
    "    if weights is not None:\n",
    "        quantile_loss = np.average(loss[~np.isnan(loss)], \n",
    "                             weights=weights[~np.isnan(loss)],\n",
    "                             axis=axis)\n",
    "    else:\n",
    "        quantile_loss = np.nanmean(loss, axis=axis)\n",
    "        \n",
    "    return quantile_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(quantile_loss, title_level=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](imgs_losses/q_loss.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-Quantile Loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| export\n",
    "def mqloss(y: np.ndarray, y_hat: np.ndarray, \n",
    "           quantiles: np.ndarray, \n",
    "           weights: Optional[np.ndarray] = None,\n",
    "           axis: Optional[int] = None) -> Union[float, np.ndarray]:\n",
    "    \"\"\"  Multi-Quantile loss\n",
    "\n",
    "    Calculates the Multi-Quantile loss (MQL) between `y` and `y_hat`.\n",
    "    MQL calculates the average multi-quantile Loss for\n",
    "    a given set of quantiles, based on the absolute \n",
    "    difference between predicted quantiles and observed values.\n",
    "\n",
    "    $$ \\mathrm{MQL}(\\\\mathbf{y}_{\\\\tau},[\\\\mathbf{\\hat{y}}^{(q_{1})}_{\\\\tau}, ... ,\\hat{y}^{(q_{n})}_{\\\\tau}]) = \\\\frac{1}{n} \\\\sum_{q_{i}} \\mathrm{QL}(\\\\mathbf{y}_{\\\\tau}, \\\\mathbf{\\hat{y}}^{(q_{i})}_{\\\\tau}) $$\n",
    "\n",
    "    The limit behavior of MQL allows to measure the accuracy \n",
    "    of a full predictive distribution $\\mathbf{\\hat{F}}_{\\\\tau}$ with \n",
    "    the continuous ranked probability score (CRPS). This can be achieved \n",
    "    through a numerical integration technique, that discretizes the quantiles \n",
    "    and treats the CRPS integral with a left Riemann approximation, averaging over \n",
    "    uniformly distanced quantiles.    \n",
    "\n",
    "    $$ \\mathrm{CRPS}(y_{\\\\tau}, \\mathbf{\\hat{F}}_{\\\\tau}) = \\int^{1}_{0} \\mathrm{QL}(y_{\\\\tau}, \\hat{y}^{(q)}_{\\\\tau}) dq $$\n",
    "\n",
    "    **Parameters:**<br>\n",
    "    `y`: numpy array, Actual values.<br>\n",
    "    `y_hat`: numpy array, Predicted values.<br>\n",
    "    `quantiles`: numpy array,(n_quantiles). Quantiles to estimate from the distribution of y.<br>\n",
    "    `mask`: numpy array, Specifies date stamps per serie to consider in loss.<br>\n",
    "\n",
    "    **Returns:**<br>\n",
    "    `mqloss`: numpy array, (single value).\n",
    "    \n",
    "    **References:**<br>\n",
    "    [Roger Koenker and Gilbert Bassett, Jr., \"Regression Quantiles\".](https://www.jstor.org/stable/1913643)<br>\n",
    "    [James E. Matheson and Robert L. Winkler, \"Scoring Rules for Continuous Probability Distributions\".](https://www.jstor.org/stable/2629907)\n",
    "    \"\"\"\n",
    "    if weights is None: weights = np.ones(y.shape)\n",
    "        \n",
    "    _metric_protections(y, y_hat, weights)\n",
    "    n_q = len(quantiles)\n",
    "    \n",
    "    y_rep  = np.expand_dims(y, axis=-1)\n",
    "    error  = y_hat - y_rep\n",
    "    sq     = np.maximum(-error, np.zeros_like(error))\n",
    "    s1_q   = np.maximum(error, np.zeros_like(error))\n",
    "    mqloss = (quantiles * sq + (1 - quantiles) * s1_q)\n",
    "    \n",
    "    # Match y/weights dimensions and compute weighted average\n",
    "    weights = np.repeat(np.expand_dims(weights, axis=-1), repeats=n_q, axis=-1)\n",
    "    mqloss  = np.average(mqloss, weights=weights, axis=axis)\n",
    "\n",
    "    return mqloss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_doc(mqloss, title_level=3)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](imgs_losses/mq_loss.png)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Examples and Validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import unittest\n",
    "import torch as t \n",
    "import numpy as np\n",
    "\n",
    "from neuralforecast.losses.pytorch import (\n",
    "    MAE, MSE, RMSE,      # unscaled errors\n",
    "    MAPE, SMAPE,         # percentage errors\n",
    "    MASE,                # scaled error\n",
    "    QuantileLoss, MQLoss # probabilistic errors\n",
    ")\n",
    "\n",
    "from neuralforecast.losses.numpy import (\n",
    "    mae, mse, rmse,              # unscaled errors\n",
    "    mape, smape,                 # percentage errors\n",
    "    mase,                        # scaled error\n",
    "    quantile_loss, mqloss        # probabilistic errors\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#| hide\n",
    "# Test class for pytorch/numpy loss functions\n",
    "class TestLoss(unittest.TestCase):\n",
    "    def setUp(self):   \n",
    "        self.num_quantiles = np.random.randint(3, 10)\n",
    "        self.first_num = np.random.randint(1, 300)\n",
    "        self.second_num = np.random.randint(1, 300)\n",
    "        \n",
    "        self.y = t.rand(self.first_num, self.second_num)\n",
    "        self.y_hat = t.rand(self.first_num, self.second_num)\n",
    "        self.y_hat2 = t.rand(self.first_num, self.second_num)\n",
    "        self.y_hat_quantile = t.rand(self.first_num, self.second_num, self.num_quantiles)\n",
    "        \n",
    "        self.quantiles = t.rand(self.num_quantiles)\n",
    "        self.q_float = np.random.random_sample()\n",
    "\n",
    "    def test_mae(self):\n",
    "        mae_numpy   = mae(self.y, self.y_hat)\n",
    "        mae_pytorch = MAE()\n",
    "        mae_pytorch = mae_pytorch(self.y, self.y_hat).numpy()\n",
    "        self.assertAlmostEqual(mae_numpy, mae_pytorch, places=6)\n",
    "\n",
    "    def test_mse(self):\n",
    "        mse_numpy   = mse(self.y, self.y_hat)\n",
    "        mse_pytorch = MSE()\n",
    "        mse_pytorch = mse_pytorch(self.y, self.y_hat).numpy()\n",
    "        self.assertAlmostEqual(mse_numpy, mse_pytorch, places=6)\n",
    "\n",
    "    def test_rmse(self):\n",
    "        rmse_numpy   = rmse(self.y, self.y_hat)\n",
    "        rmse_pytorch = RMSE()\n",
    "        rmse_pytorch = rmse_pytorch(self.y, self.y_hat).numpy()\n",
    "        self.assertAlmostEqual(rmse_numpy, rmse_pytorch, places=6)\n",
    "\n",
    "    def test_mape(self):\n",
    "        mape_numpy   = mape(y=self.y, y_hat=self.y_hat)\n",
    "        mape_pytorch = MAPE()\n",
    "        mape_pytorch = mape_pytorch(y=self.y, y_hat=self.y_hat).numpy()\n",
    "        self.assertAlmostEqual(mape_numpy, mape_pytorch, places=6)\n",
    "\n",
    "    def test_smape(self):\n",
    "        smape_numpy   = smape(self.y, self.y_hat)\n",
    "        smape_pytorch = SMAPE()\n",
    "        smape_pytorch = smape_pytorch(self.y, self.y_hat).numpy()\n",
    "        self.assertAlmostEqual(smape_numpy, smape_pytorch, places=4)\n",
    "    \n",
    "    #def test_mase(self):\n",
    "    #    y_insample = t.rand(self.first_num, self.second_num)\n",
    "    #    seasonality = 24\n",
    "    #    # Hourly 24, Daily 7, Weekly 52\n",
    "    #    # Monthly 12, Quarterly 4, Yearly 1\n",
    "    #    mase_numpy   = mase(y=self.y, y_hat=self.y_hat,\n",
    "    #                        y_insample=y_insample, seasonality=seasonality)\n",
    "    #    mase_object  = MASE(seasonality=seasonality)\n",
    "    #    mase_pytorch = mase_object(y=self.y, y_hat=self.y_hat,\n",
    "    #                               y_insample=y_insample).numpy()\n",
    "    #    self.assertAlmostEqual(mase_numpy, mase_pytorch, places=2)\n",
    "\n",
    "    #def test_rmae(self):\n",
    "    #    rmae_numpy   = rmae(self.y, self.y_hat, self.y_hat2)\n",
    "    #    rmae_object  = RMAE()\n",
    "    #    rmae_pytorch = rmae_object(self.y, self.y_hat, self.y_hat2).numpy()\n",
    "    #    self.assertAlmostEqual(rmae_numpy, rmae_pytorch, places=4)\n",
    "\n",
    "    def test_quantile(self):\n",
    "        quantile_numpy = quantile_loss(self.y, self.y_hat, q = self.q_float)\n",
    "        quantile_pytorch = QuantileLoss(q = self.q_float)\n",
    "        quantile_pytorch = quantile_pytorch(self.y, self.y_hat).numpy()\n",
    "        self.assertAlmostEqual(quantile_numpy, quantile_pytorch, places=6)\n",
    "    \n",
    "    # def test_mqloss(self):\n",
    "    #     weights = np.ones_like(self.y)\n",
    "\n",
    "    #     mql_np_w = mqloss(self.y, self.y_hat_quantile, self.quantiles, weights=weights)\n",
    "    #     mql_np_default_w = mqloss(self.y, self.y_hat_quantile, self.quantiles)\n",
    "\n",
    "    #     mql_object = MQLoss(quantiles=self.quantiles)\n",
    "    #     mql_py_w = mql_object(y=self.y,\n",
    "    #                           y_hat=self.y_hat_quantile,\n",
    "    #                           mask=t.Tensor(weights)).numpy()\n",
    "        \n",
    "    #     print('self.y.shape', self.y.shape)\n",
    "    #     print('self.y_hat_quantile.shape', self.y_hat_quantile.shape)\n",
    "    #     mql_py_default_w = mql_object(y=self.y,\n",
    "    #                                   y_hat=self.y_hat_quantile).numpy()\n",
    "\n",
    "    #     weights[0,:] = 0\n",
    "    #     mql_np_new_w = mqloss(self.y, self.y_hat_quantile, self.quantiles, weights=weights)\n",
    "    #     mql_py_new_w = mql_object(y=self.y,\n",
    "    #                               y_hat=self.y_hat_quantile,\n",
    "    #                               mask=t.Tensor(weights)).numpy()\n",
    "\n",
    "    #     self.assertAlmostEqual(mql_np_w,  mql_np_default_w)\n",
    "    #     self.assertAlmostEqual(mql_py_w,  mql_py_default_w)\n",
    "    #     self.assertAlmostEqual(mql_np_new_w,  mql_py_new_w)\n",
    "    \n",
    "\n",
    "unittest.main(argv=[''], verbosity=2, exit=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "python3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}