colab_3_PredictiveMaintenance_workshop.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "586a8684-8061-4a7e-aed6-c61518987bc9",
   "metadata": {},
   "source": [
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f717c58-97d5-4468-a89b-20292cf35ad4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# setup darts\n",
    "!pip install \"darts @ git+https://github.com/unit8co/darts.git@master\" seaborn openpyxl catboost lightgbm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e43c00c-b78b-4c57-ac14-0e8a34387ae1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download data from data source to local\n",
    "import io\n",
    "import os\n",
    "import requests\n",
    "import shutil\n",
    "import zipfile\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62d1e582-c306-4aaf-b8cc-e484c00ea662",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = \"data\"\n",
    "pdm_dir = os.path.join(data_dir, \"pdm\")\n",
    "for dir_path in [data_dir, pdm_dir]:\n",
    "    if not os.path.exists(dir_path):\n",
    "        os.mkdir(dir_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fc8b339-a787-4a33-9b5b-c084e636aa51",
   "metadata": {},
   "source": [
    "## Data Download and Preprocessing - Predictive Maintenance - Wind Turbine Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fca327ee-a94a-48c3-ac01-44a62dee2d7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "FREQ = pd.tseries.frequencies.to_offset(\"10min\")\n",
    "\n",
    "def preprocess_data(df):\n",
    "    \"\"\"Preprocess the turbine signals\n",
    "    - drops duplicate time stamps (wrong times from daylight saving, no chance to identify correct times)\n",
    "    - resamples the data to a contiguous 10 Min DatetimeIndex\n",
    "    - flags missing dates with feature \"missin_date\"\n",
    "    - computes a generator/turbine on/off timer for each turbine.\n",
    "    \"\"\"\n",
    "    # drop duplicate time stamps (wrong times from daylight saving, no chance to\n",
    "    # identify correct times)\n",
    "    df = df.drop_duplicates([\"Turbine_ID\", \"Timestamp\"], keep=\"first\").reset_index(drop=True)\n",
    "    # remove time zone information\n",
    "    df[\"Timestamp\"] = pd.DatetimeIndex(df[\"Timestamp\"]).tz_localize(None)\n",
    "    df = df.sort_values(by=[\"Turbine_ID\", \"Timestamp\"])\n",
    "    df = df.set_index(\"Timestamp\")\n",
    "    df = df.groupby(\"Turbine_ID\").apply(compute_on_off_timer).reset_index(\"Turbine_ID\", drop=True)\n",
    "    return df.reset_index()\n",
    "\n",
    "\n",
    "def compute_on_off_timer(df):\n",
    "    \"\"\"Computes the time (number of 10 minute steps) since the generator last crossed the \n",
    "    1200 rotations per minute (RPM) mark.\n",
    "\n",
    "    This is a useful feature to let the model know how long that the turbine has been active/inactive,\n",
    "    and should ultimately improve generator temperature modelling.\n",
    "\n",
    "    - positive counter when the generator went from below 1200 RPM to above\n",
    "    - negative counter when the generator went from above 1200 RPM to below\n",
    "    - upper- and lower-bound by +/- 4 hours (4 * 6 (10 minutes steps) = 24)\n",
    "    \"\"\"\n",
    "    \n",
    "    df = df.ffill()\n",
    "    df = df.asfreq(FREQ)\n",
    "    df[\"missing_date\"] = df.isna().any(axis=1)\n",
    "    df = df.ffill()\n",
    "\n",
    "    # let's use a proxy for generator on if RPM >= 1200\n",
    "    df[\"state_on\"] = df[\"Gen_RPM_Avg\"] >= 1200\n",
    "    df[\"state_off\"] = df[\"Gen_RPM_Avg\"] < 1200\n",
    "\n",
    "    # count how many steps since last time the rotation speed went above threshold (and limit the counter)\n",
    "    b = (~df[\"state_on\"]).cumsum()\n",
    "    df[\"timer_on\"] = df.groupby(b)[\"state_on\"].cumsum()\n",
    "    df.loc[df[\"timer_on\"] > 4 * 6, \"timer_on\"] = 4 * 6\n",
    "\n",
    "    # count how many steps since last time the rotation speed went below threshold (and limit the counter)\n",
    "    b = (~df[\"state_off\"]).cumsum()\n",
    "    df[\"timer_off\"] = -(df.groupby(b)[\"state_off\"]).cumsum()\n",
    "    df.loc[df[\"timer_off\"] < -4 * 6, \"timer_off\"] = -4 * 6\n",
    "\n",
    "    # combine on and off timers\n",
    "    df[\"timer_on_off\"] = np.where(df[\"timer_on\"] > 0., df[\"timer_on\"], df[\"timer_off\"])\n",
    "\n",
    "    # remove the dedicated counters\n",
    "    return df.drop(columns=[\"timer_on\", \"timer_off\", \"state_off\"])\n",
    "\n",
    "\n",
    "def preprocess_anomalies(df):\n",
    "    df[\"Timestamp\"] = pd.DatetimeIndex(df[\"Timestamp\"]).tz_localize(None).round(FREQ.freqstr)\n",
    "    df[\"start\"] = df[\"Timestamp\"]\n",
    "    df[\"end\"] = df[\"start\"] + (6 * 24 - 1) * FREQ\n",
    "    df = df.sort_values(by=[\"Turbine_ID\", \"Timestamp\"]).reset_index(drop=True)\n",
    "    return df\n",
    "\n",
    "\n",
    "def download_wind_turbine():\n",
    "    # URL of the files\n",
    "    url_data = \"https://www.edp.com/sites/default/files/2023-04/Wind-Turbine-SCADA-signals-2016.xlsx\"\n",
    "    url_failures = \"https://www.edp.com/sites/default/files/2023-04/Historical-Failure-Logbook-2016.xlsx\"\n",
    "    \n",
    "    for url, preproc_fn in zip([url_data, url_failures], [preprocess_data, preprocess_anomalies]):\n",
    "        file_path = os.path.join(pdm_dir, url.split(\"/\")[-1])\n",
    "        csv_path = file_path.replace(\"xlsx\", \"csv\")\n",
    "        if not (os.path.exists(file_path) or os.path.exists(csv_path)):\n",
    "            # Send a GET request to download the zip file\n",
    "            response = requests.get(url)\n",
    "    \n",
    "            # Check if the request was successful\n",
    "            if response.status_code == 200:\n",
    "                # Save the zip file to the local drive\n",
    "                with open(file_path, \"wb\") as file:\n",
    "                    file.write(response.content)\n",
    "                print(\"File downloaded successfully.\")\n",
    "            else:\n",
    "                print(\"Failed to download.\")\n",
    "        else:\n",
    "            print(\"File already downloaded.\")\n",
    "    \n",
    "        # Extract the zip file\n",
    "        if os.path.exists(file_path):\n",
    "            print(\"Converting to csv..\")\n",
    "            df = pd.read_excel(file_path)\n",
    "            df = preproc_fn(df)\n",
    "            df.to_csv(csv_path, index=False)\n",
    "            print(\"Successfully converted to csv.\")\n",
    "            os.remove(file_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f18b2f96-5003-4b2b-bddf-1ca7aed76d52",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Downloading Wind Turbine Data..\")\n",
    "download_wind_turbine()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4edc3f1-bcb3-4d4d-9fca-8662ec827540",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52c5747e-fbbb-49a5-82f0-687dad2f25e6",
   "metadata": {},
   "source": [
    "## Helper Functions (no need to look at them for now)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "608cfc9c-97a3-472d-a706-d0f510235c4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from typing import List, Optional\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from darts import concatenate, TimeSeries\n",
    "from darts.dataprocessing.transformers import Scaler\n",
    "from darts.utils.data.tabularization import create_lagged_data\n",
    "from darts.utils.utils import generate_index\n",
    "\n",
    "\n",
    "FREQ = pd.tseries.frequencies.to_offset(\"10min\")\n",
    "DEV_LO, DEV_MD, DEV_HI, DEV_AGG = [\"dev_lo\", \"dev_md\", \"dev_hi\", \"dev_agg\"]\n",
    "ANOM_LO, ANOM_HI, ANOM_AGG = [\"anom_lo\", \"anom_hi\", \"anom_agg\"]\n",
    "\n",
    "\n",
    "def df_anom_to_series(df: pd.DataFrame):\n",
    "    \"\"\"Converts the wind turbine failure DataFrame (rows with anomaly start and end dates)\n",
    "    into a TimeSeries of binary anomalies.\"\"\"\n",
    "    # create time series with binary anomalies\n",
    "    time_index = generate_index(\n",
    "        start=pd.Timestamp(df[\"start\"].min()),\n",
    "        end=pd.Timestamp(df[\"end\"].max()),\n",
    "        freq=FREQ\n",
    "    ).tz_localize(None)\n",
    "\n",
    "    turbine_ids = sorted(df[\"Turbine_ID\"].unique())\n",
    "    df_anom = pd.DataFrame(\n",
    "        index=time_index,\n",
    "        data={t_id: 0. for t_id in turbine_ids},\n",
    "    )\n",
    "\n",
    "    for t_id in turbine_ids:\n",
    "        anom_times = df.loc[df[\"Turbine_ID\"] == t_id, [\"start\", \"end\"]]\n",
    "        for start, end in anom_times.values:\n",
    "            df_anom.loc[start:end, t_id] = 1.\n",
    "\n",
    "    return TimeSeries.from_dataframe(df_anom)\n",
    "\n",
    "\n",
    "def compute_anomalies(\n",
    "    model,\n",
    "    series: TimeSeries,\n",
    "    pred_series: TimeSeries,\n",
    "    quantiles: List[float],\n",
    "    min_value: float = 0.,\n",
    "    anom_window: int = 3,\n",
    "    min_anom_prob: float = 1.0,\n",
    "):\n",
    "    \"\"\"\n",
    "    Computes the anomalies based on the predicted normal operating range and the actual values of the feature.\n",
    "\n",
    "    - Ignores points which were only slightly outside the interval (threshold `min_value`)\n",
    "    - Scans the residuals in fixed size windows and counts how many points were out-of bounds in each window\n",
    "      (`anom_window`)\n",
    "    - For each window we can set a minimum out-of-bounds probability below which we do not consider the window as\n",
    "      anomalous (`min_anom_prob`)\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    model\n",
    "        A trained probabilistic forecasting model that predicts three quantiles.\n",
    "    series\n",
    "        The actual target series values\n",
    "    pred_series\n",
    "        The predicted target series quantile values\n",
    "    quantiles\n",
    "         The predicted quantiles.\n",
    "    min_value\n",
    "        Ignore points where the actual values are less than `min_value` outside the predicted interval.\n",
    "    anom_window\n",
    "        The window size to detect anomalies on.\n",
    "    min_anom_prob\n",
    "        For each window it is the minimum value for the fraction between the number of points that are\n",
    "        out-of-bounds and `anom_window`.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    anom_pred\n",
    "        the residuals (\"dev_lo\", ...) and final anomaly flags (\"anom_lo\", ...) as a `TimeSeries`.\n",
    "    df_anom_pred\n",
    "        a `DataFrame` where each row represents an anomaly, including the start, end date and some statistics\n",
    "    ql\n",
    "        the quantile loss for the historical forecast (as a metric how good the quantile predictions are, 0.\n",
    "        being the best score).\n",
    "    \"\"\"\n",
    "    n_quantiles = pred_series.n_components // series.n_components\n",
    "    # repeat target column so when can compute the residuals per quantile\n",
    "    series_ext = concatenate([series] * n_quantiles, axis=1)\n",
    "    residuals = model.residuals(\n",
    "        series=series_ext,\n",
    "        historical_forecasts=pred_series,\n",
    "        last_points_only=True,\n",
    "    )\n",
    "    # quantile loss (metric) per quantile\n",
    "    ql = quantile_loss(residuals, quantiles)\n",
    "\n",
    "    # ignore residuals where y_true was within high and low quantile\n",
    "    df = residuals.pd_dataframe(copy=True)\n",
    "    df.columns = [DEV_LO, DEV_MD, DEV_HI]\n",
    "\n",
    "    df = df.fillna(0.)\n",
    "    df.loc[df[DEV_LO] > -min_value, DEV_LO] = 0.\n",
    "    df.loc[df[DEV_HI] < min_value, DEV_HI] = 0.\n",
    "    df.loc[:, DEV_AGG] = df[DEV_LO] + df[DEV_HI]\n",
    "\n",
    "    # ignore residuals where the anomaly didn't last for a couple time steps\n",
    "    df = _find_anomaly_periods(df, anom_window, min_anom_prob)\n",
    "    # generate a dataframe where each row\n",
    "    anom = _compute_anomaly_table(df)\n",
    "    return TimeSeries.from_dataframe(df), anom, ql\n",
    "\n",
    "\n",
    "def _find_anomaly_periods(\n",
    "    df: pd.DataFrame,\n",
    "    anom_window: int = 3,\n",
    "    min_anom_prob: float = 1.0,\n",
    "):\n",
    "    # ignore residuals where the anomaly didn't last for a couple time steps\n",
    "    df_out = df.copy()\n",
    "    for col_dev, col_anom in zip([DEV_LO, DEV_HI, DEV_AGG], [ANOM_LO, ANOM_HI, ANOM_AGG]):\n",
    "        is_anomaly = df_out[col_dev] != 0\n",
    "\n",
    "        # windowed anomalies\n",
    "        windows_anom = create_lagged_data(\n",
    "            target_series=TimeSeries.from_series(is_anomaly),\n",
    "            lags=[i for i in range(-anom_window, 0)],\n",
    "            uses_static_covariates=False,\n",
    "            is_training=False,\n",
    "        )[0]\n",
    "        # windowing results in n - anom_window windows -> repeat the first window and prepend\n",
    "        windows_anom = np.concatenate([\n",
    "            np.zeros((anom_window - 1, anom_window, 1)),\n",
    "            windows_anom,\n",
    "        ])\n",
    "        # get the average number of anomalous steps per window\n",
    "        windows_anom = windows_anom.mean(axis=1)[:, 0]\n",
    "        # compute the actual anomalous time frames, with respect to probability\n",
    "        idx_anom = np.argwhere(windows_anom >= min_anom_prob)[:, 0]\n",
    "\n",
    "        # reset anomaly flags\n",
    "        windows_anom = np.zeros(shape=windows_anom.shape)\n",
    "        for i in range(anom_window):\n",
    "            windows_anom[idx_anom - i] = 1.\n",
    "\n",
    "        windows_anom = pd.Series(windows_anom, index=is_anomaly.index).astype(bool)\n",
    "        # add binary anomalies\n",
    "        df_out.loc[:, col_anom] = windows_anom\n",
    "        # remove deviations from too short anomalies\n",
    "        df_out.loc[~windows_anom, col_dev] = 0.0\n",
    "    return df_out\n",
    "\n",
    "\n",
    "def _compute_anomaly_table(df: pd.DataFrame) -> pd.DataFrame:\n",
    "    \"\"\"Computes start and end dates for each anomaly in `df`. The returned DataFrame has columns:\n",
    "        - \"start\": anomaly time stamp\n",
    "        - \"end\": anomaly end time stamp\n",
    "        - \"n_steps\": how many steps the anomaly lasted\n",
    "        - \"sum\": sum of all deviations (y_true above the high and/or below the low predicted quantile)\n",
    "        - \"name\" \"anom_lo\", \"anom_hi\" or \"anom_agg\"\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    df\n",
    "        the pandas DataFrame containing anomaly columns `columns`\n",
    "    \"\"\"\n",
    "    interval_dfs = []\n",
    "    for col in [ANOM_LO, ANOM_HI, ANOM_AGG]:\n",
    "        # make [0, 0, 1, 1, 0] groupable -> [0, 0, 1, 1, 2]\n",
    "        blocks = df[col].diff().ne(0).cumsum()\n",
    "        # remove groups with initial zeros, so that we can group only the ones-groups\n",
    "        blocks = blocks[df[col] > 0]\n",
    "        # group all ones and get the start and end dates (index)\n",
    "        anom_df = df.loc[df[col] > 0]\n",
    "        interval_df = (\n",
    "            anom_df[col].index.to_frame()\n",
    "            .groupby(blocks, sort=True)\n",
    "            .agg([\"min\", \"max\", \"size\"])\n",
    "            .rename(columns={\"min\": \"start\", \"max\": \"end\", \"size\": \"n_steps\"})\n",
    "            .reset_index(drop=True)\n",
    "        )\n",
    "        # drop the time tag from columns\n",
    "        interval_df.columns = interval_df.columns.droplevel(0)\n",
    "        # deviance column\n",
    "        dev_col = col.replace(\"anom\", \"dev\")\n",
    "        dev_df = (\n",
    "            anom_df[[dev_col]]\n",
    "            .groupby(blocks, sort=True)\n",
    "            .agg(\"sum\")\n",
    "            .reset_index(drop=True)\n",
    "            .rename(columns={dev_col: \"dev_tot\"})\n",
    "        )\n",
    "        interval_df = pd.concat([interval_df, dev_df], axis=1)\n",
    "        interval_df[\"name\"] = col\n",
    "        interval_dfs.append(interval_df)\n",
    "    return pd.concat(interval_dfs, ignore_index=True)\n",
    "\n",
    "\n",
    "def quantile_loss(residuals: TimeSeries, quantiles: List[float]):\n",
    "    \"\"\"Computes the quantile loss per predicted quantile.\"\"\"\n",
    "    errors = residuals.values(copy=False)\n",
    "    qs = np.array(quantiles)\n",
    "    losses = 2.0 * np.maximum((qs - 1) * errors, qs * errors)\n",
    "    return pd.DataFrame(np.nanmean(losses, axis=0), index=quantiles).T\n",
    "\n",
    "\n",
    "def plot_predicted_anomalies(\n",
    "    df: pd.DataFrame,\n",
    "    series: TimeSeries,\n",
    "    covs: TimeSeries,\n",
    "    hist_fc: TimeSeries,\n",
    "    anom_true: TimeSeries,\n",
    "    anom_pred: TimeSeries,\n",
    "    turbine_id: str,\n",
    "    max_plots: int = 5,\n",
    "):\n",
    "    \"\"\"For the first `max_plots` predicted anomaly in `df`, it plots the preceding and following 3 days of:\n",
    "\n",
    "    - actual target values\n",
    "    - historical forecast boundaries\n",
    "    - predicted anomaly as a green capsule\n",
    "    - actual anomaly as a red capsule (if there is any in the time frame)\n",
    "    - and the covariates scaled to a value range (0, 1).\n",
    "    \"\"\"\n",
    "\n",
    "    for idx, (start, end) in enumerate(df[[\"start\", \"end\"]].values):\n",
    "        start = pd.Timestamp(start)\n",
    "        end = pd.Timestamp(end)\n",
    "        p_start = start - FREQ * 6 * 24 * 3\n",
    "        p_end = end + FREQ * 6 * 24 * 3\n",
    "        fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(12, 9.6), sharex=True)\n",
    "\n",
    "        series[p_start:p_end].plot(ax=ax1)\n",
    "        plot_intervals(hist_fc[p_start:p_end], ax=ax1, plot_med=True, alpha=0.5)\n",
    "        plot_actual_anomalies(anom_true[turbine_id][p_start:p_end] * 100, ax=ax1)\n",
    "        plot_actual_anomalies(anom_pred[ANOM_AGG][p_start:p_end] * 100, ax=ax1)\n",
    "        anom_pred[DEV_AGG][p_start:p_end].plot(ax=ax1)\n",
    "\n",
    "        Scaler().fit_transform(covs)[p_start:p_end].plot(ax=ax2)\n",
    "        ax1.set_title(f\"Anom start: {start.round(FREQ)}, end: {end.round(FREQ)}\")\n",
    "        plt.show()\n",
    "        if idx == max_plots - 1:\n",
    "            break\n",
    "\n",
    "\n",
    "def plot_intervals(\n",
    "    series: TimeSeries,\n",
    "    ax=None,\n",
    "    plot_med: bool = True,\n",
    "    alpha: float=0.25,\n",
    "    c: Optional[str] = None\n",
    "):\n",
    "    \"\"\"Plot historical quantile forecasts as intervals.\"\"\"\n",
    "\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "\n",
    "    vals = series.values(copy=False)\n",
    "    if plot_med:\n",
    "        median_p = ax.plot(\n",
    "            series.time_index,\n",
    "            vals[:, 1],\n",
    "            label=series.columns[1]\n",
    "        )\n",
    "    else:\n",
    "        median_p = ax.plot([], [])\n",
    "    color_used = c or median_p[0].get_color()\n",
    "\n",
    "    ax.fill_between(\n",
    "        series.time_index,\n",
    "        vals[:, 0],\n",
    "        vals[:, -1],\n",
    "        color=color_used,\n",
    "        alpha=alpha,\n",
    "    )\n",
    "    ax.legend()\n",
    "\n",
    "\n",
    "def plot_actual_anomalies(\n",
    "    series: TimeSeries,\n",
    "    ax=None,\n",
    "    alpha: float = 0.25,\n",
    "    c: Optional[str] = None\n",
    "):\n",
    "    \"\"\"Plots a binary anomaly `series` as capsule.\"\"\"\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots()\n",
    "\n",
    "    vals = series.values(copy=False)\n",
    "    vals_zero = np.zeros(vals.shape)\n",
    "    if c is None:\n",
    "        empty_p = ax.plot([], [])\n",
    "        color_used = empty_p[0].get_color()\n",
    "    else:\n",
    "        color_used = c\n",
    "\n",
    "    ax.fill_between(\n",
    "        series.time_index,\n",
    "        vals_zero[:, 0],\n",
    "        vals[:, 0],\n",
    "        color=color_used,\n",
    "        alpha=alpha,\n",
    "        label=series.columns[0] + \"_anom\",\n",
    "    )\n",
    "    ax.legend()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b853461-7c77-4bab-b324-8743d8ae2cb5",
   "metadata": {},
   "source": [
    "# Predictive Maintenance With Darts\n",
    "\n",
    "In this notebook, we'll show how to perform a predictive maintenance task using Darts on the example of Wind Turbine Failures. We will look at:\n",
    "- which features to use (and not use) for modelling\n",
    "- how to setup and train a model to predict the expected normal operating range of a signal\n",
    "- how to detect anomalies from deviations between the actual signal values and the predicted operating range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dfe57351-d648-4911-8f62-da65bcc5e0b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from darts import TimeSeries"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb1caded-1990-44b7-a9f7-53df15908fcd",
   "metadata": {},
   "source": [
    "## Data Loading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "270c7738-8285-4611-98f7-164bdec6b129",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "\n",
    "# the data has a 10 minute frequency\n",
    "FREQ = pd.tseries.frequencies.to_offset(\"10min\")\n",
    "\n",
    "data_dir = os.path.join(\"data\", \"pdm\")\n",
    "fpath_data = os.path.join(data_dir, f\"Wind-Turbine-SCADA-signals-2016.csv\")\n",
    "fpath_failures = os.path.join(data_dir, f\"Historical-Failure-Logbook-2016.csv\")\n",
    "\n",
    "# get the turbine signal data\n",
    "df_data = pd.read_csv(fpath_data)\n",
    "df_data.index = pd.DatetimeIndex(df_data[\"Timestamp\"])\n",
    "df_data = df_data.drop(columns=\"Timestamp\")\n",
    "\n",
    "# get the failure logbook\n",
    "df_anom_true = pd.read_csv(fpath_failures)\n",
    "df_anom_true[\"start\"] = pd.to_datetime(df_anom_true[\"start\"])\n",
    "df_anom_true[\"end\"] = pd.to_datetime(df_anom_true[\"end\"])\n",
    "\n",
    "# convert the logbook into a TimeSeries of binary anomalies (using helper funtion `df_anom_to_series`)\n",
    "anom_true = df_anom_to_series(df_anom_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "199526a6-ef49-42de-93f5-656b51c871fc",
   "metadata": {},
   "source": [
    "## Data Investigation\n",
    "\n",
    "### Generator Anomalies\n",
    "\n",
    "Let's look at the anomalies of Turbine 6 and 7."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca86e903-4162-45f9-892f-5a137e264f25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "\n",
    "df_anom_true[df_anom_true[\"Turbine_ID\"] == \"T06\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bd55c33-f832-4a6e-9f9e-c6daf0804af5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "\n",
    "df_anom_true[df_anom_true[\"Turbine_ID\"] == \"T07\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ac5c2da-e1a8-4077-8e58-e8413b9eca5d",
   "metadata": {},
   "source": [
    "It seems there were many anomalies related to the generator. There were 5 anomalies related to the generator in turbine 6 and 1 in turbine 7. Also, most of the anomalies indicate that there was an issue with high temperatures, or malfunctioning temperature sensors. \n",
    "\n",
    "For this exercise, we'll only focus on anomalies related to the generator. Specifically, let's try to model one of the generator temperature signals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad380093-7ad6-4964-9d5f-f0564679abc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "\n",
    "# select only generator speicifc anomalies\n",
    "df_anom_true_gen = df_anom_true[df_anom_true[\"Component\"].str.contains(\"GENERATOR\")].reset_index(drop=True)\n",
    "\n",
    "# convert the logbook into a TimeSeries of binary anomalies\n",
    "anom_true_gen = df_anom_to_series(df_anom_true_gen)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b97bc8d7-3f4e-49cc-bd93-3f2f43eb99dd",
   "metadata": {},
   "source": [
    "### Target Variable Analysis (Generator Bearing Temperature)\n",
    "\n",
    "We select the target and covariates features as discussed in the presentation.\n",
    "\n",
    "The generator bearing temperature is the target feature/signal that we want to model. The location of the generator bearing can be found in [hello](#External-Feature-Analysis).\n",
    "\n",
    "To make the code a bit easier, we already extract the covariates here as well, and refer to turbine 6 as `train` and turbine 7 as `test`.\n",
    "\n",
    "Now let's extract the signal data per turbine and create a TimeSeries from it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43401d28-b6e1-4112-b739-992733749553",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "\n",
    "# train and test turbine IDs\n",
    "turbine_id_train, turbine_id_test = \"T06\", \"T07\"\n",
    "\n",
    "# target feature (we explain more later)\n",
    "tg_col = 'Gen_Bear_Temp_Avg'\n",
    "\n",
    "# external features (we explain more later)\n",
    "cov_cols = [\n",
    "    \"Gen_RPM_Max\",\n",
    "    \"Nac_Temp_Avg\",\n",
    "    \"timer_on_off\",\n",
    "    \"missing_date\",\n",
    "    \"Prod_LatestAvg_ActPwrGen0\",\n",
    "    \"Prod_LatestAvg_ActPwrGen1\",\n",
    "]\n",
    "\n",
    "# extract turbine specific signal data\n",
    "df_train = df_data.loc[df_data[\"Turbine_ID\"] == turbine_id_train, cov_cols + [tg_col]].copy()\n",
    "df_test = df_data.loc[df_data[\"Turbine_ID\"] == turbine_id_test, cov_cols + [tg_col]].copy()\n",
    "\n",
    "# extract the TimeSeries\n",
    "train_raw = TimeSeries.from_dataframe(df_train, fill_missing_dates=True, freq=FREQ)\n",
    "test_raw = TimeSeries.from_dataframe(df_test, fill_missing_dates=True, freq=FREQ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3ea669c-f96c-48d5-bfe2-442d316f94c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "\n",
    "def plot_series(train, test, anom_true):\n",
    "    fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(12, 6), sharex=True)\n",
    "    train[tg_col].plot(ax=ax1)\n",
    "    (anom_true[turbine_id_train] * 150).plot(label=\"is_anomaly\", alpha=0.6, ax=ax1)\n",
    "    ax1.set_title(\"Anomalies Train Set (Turbine 6)\")\n",
    "    \n",
    "    test[tg_col].plot(ax=ax2)\n",
    "    (anom_true[turbine_id_test] * 150).plot(label=\"is_anomaly\", alpha=0.6, ax=ax2)\n",
    "    ax2.set_title(\"Anomalies Test Set (Turbine 7)\")\n",
    "    fig.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "plot_series(train_raw, test_raw, anom_true_gen)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ca28b84-d917-46dc-8f2b-8b42efa2336b",
   "metadata": {},
   "source": [
    "> We reference the anomalies here as 1,2, ... from left to right\n",
    "\n",
    "**Turbine 6:**\n",
    "- Anomalies 1 and 5 in Turbine 6 were generator replacements. They are followed by a short time when the turbine was inactive (missing dates).\n",
    "- Anomaly 2 was a temperature sensor failure, but it seems it could have been detected earlier\n",
    "- Anomalies 3 (high temperature generator error) and 4 (refrigeration system and temperature sensor replaced) seem to be less obvious from a visual inspection\n",
    "\n",
    "**Turbine 7:**\n",
    "- There was only one anomaly related to the generator (high temperature in generator bearing). It also clearly visible, the predictive maintenance model should have an easy time finding this.\n",
    "\n",
    "And now a quick glance at the distribution of the temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de6559e0-1546-4ad4-a4d3-60e693ed2e59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(12, 6))\n",
    "df_targets = pd.concat([df_train[tg_col], df_test[tg_col]], axis=1)\n",
    "df_targets.columns += [\"_\" + turbine_id_train, \"_\" + turbine_id_test]\n",
    "df_targets.boxplot(ax=ax1)\n",
    "df_targets.plot.hist(bins=100, alpha=0.5, density=True, ax=ax2)\n",
    "df_targets.median()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e3868c2-2605-40fd-bb06-f35b3c7ab597",
   "metadata": {},
   "source": [
    "From the above, we can already spot the outliers at 200°C, which are related to high temperature errors (and usually followed by a replacement of the sensor)\n",
    "\n",
    "Also, we could estimate here that normal temperature range lies between roughly 20 - 100°C. We'll later see how to make the model learn this.\n",
    "\n",
    "If we went with only the distribution itself, we could only mark outliers as anomalies. However, this is not enough to detect anomalies 1, 3, and 4 in turbine 6. For these, we would have to analyze the temperature in relation to the other signals/measurements.\n",
    "\n",
    "And that is exactly what we want to do in this exercise. We want to build a predictive maintenance model, which estimates the normal operating range of the temperature at any time `t` based on some input signals."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd9233b0-6a89-4339-b60e-f5ce2b665d23",
   "metadata": {},
   "source": [
    "### External Feature Analysis\n",
    "\n",
    "Predictive maintenance is usually performed by analyzing device behavior in the most recent past (post hoc). In case there is abnormal device behavior, the model/tool should raise an alarm before serious damage occurs.\n",
    "\n",
    "Because it's done post hoc, we can model for example the generator temperature at all times `t` with actual measurement data from other sensors at times `t` (in regular forecasting, we wouldn't know these values as they lie in the future). \n",
    "\n",
    "For this example, we want to estimate the normal operating range for the generator bearing temperature only based on **uni-directional causal signals**.\n",
    "\n",
    "> Uni-directional causal: We'll only use external features (covariates) that can cause the temperature to change but which are not themselves affected by a change in temperature.\n",
    ">\n",
    "> Example of a **good feature**: the rotational speed of the generator shaft. The kinetic energy causes heat buildup in the generator. The rotational speed is also related with the amount of generated power (heat source).\n",
    "> \n",
    "> Example of a **bad feature**: A temperature sensor that is located close to the one we're modelling, would most likely also be affected by a temperature anomaly in the generator. Using such signals as model input would make the anomalous behavior predictable, which is not what we want. We only want the normal operating range.\n",
    "\n",
    "\n",
    "The pre-selected target and external features / measurements (covariates) are listed and shown in images below:\n",
    "\n",
    "- **Gen_Bearing_Temp_Avg [°C]**: (target) Average generator bearing temperature.\n",
    "- **Gen_RPM_Max [RPM]**: Maximum rotations per minute in the 10 minutes interval.\n",
    "- **Nac_Temp_Avg [°C]**: Average temperature of the nacelle (the housing).\n",
    "- **timer_on_off [-]**: Generated feature; counts the number of steps since the last time Gen_RPM_Max crossed the 1200 RPM mark (positive if it went from below 1200 RPM to above, and negative otherwise). It has an upper- and lower threshold of +/- 24 (=4 hours in 10 minute steps). It can be seen as a proxy for heat build-up and stagnation over time.\n",
    "- **missing_date [-]**: Generated feature; whether the date was missing in the dataset.\n",
    "- **Prod_LatestAvg_ActPwerGen0/1 [W]**: Average Power Production by generator 0/1.\n",
    "\n",
    "![covariates-illustration](./images/wind-turbine-components.jpg)\n",
    "![covariates-illustration](./images/wind_turb_gen_bearing.jpg)\n",
    "\n",
    "Let's investigate the covariates for turbine 6. We'll show:\n",
    "- Covariates over time against the temperature (we scale the covariates to a value range of (0, 1) to make them better visible)\n",
    "- Correlation Heatmap with the temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d60a0289-edbb-4ca6-ba37-be54b908d2da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "\n",
    "from darts.dataprocessing.transformers import Scaler\n",
    "import seaborn as sb\n",
    "\n",
    "times = slice(1000, 2000)\n",
    "\n",
    "def plot_feature_analysis(series, tg_col, cov_cols):\n",
    "    fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(12, 6), sharex=True)\n",
    "    \n",
    "    # target over time\n",
    "    series[tg_col][times].plot(ax=ax1)\n",
    "    ax1.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    \n",
    "    # scale covariates to (0, 1) for better visibility\n",
    "    Scaler().fit_transform(series[cov_cols])[times].plot(ax=ax2)\n",
    "    ax1.set_title(\"Temperature\")\n",
    "    ax2.set_title(\"External Features\")\n",
    "    ax2.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "    fig.tight_layout()\n",
    "    plt.show()\n",
    "    \n",
    "    # correlation heatmap\n",
    "    fig, ax = plt.subplots()\n",
    "    df_feats = pd.concat([series[tg_col].pd_dataframe(), series[cov_cols].pd_dataframe()], axis=1) \n",
    "    dataplot = sb.heatmap(\n",
    "        df_feats.corr(), cmap=\"RdBu\", annot=True, ax=ax, vmin=-1, vmax=1)\n",
    "    ax.set_title(\"Correlation heatmap\")\n",
    "    plt.show()\n",
    "\n",
    "plot_feature_analysis(train_raw, tg_col, cov_cols)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5bbfcd88-196d-4e5d-8ac6-34997f195273",
   "metadata": {},
   "source": [
    "All covariates show positive correlation with the generator bearing temperature.\n",
    "The nacelle temperature (e.g. changes in temperature during the day and throughout the year) and power generation have the highest positive correlation.\n",
    "\n",
    "Features Gen_RPM_Max and Prod_LatestAvg_ActPwrGen0/1 are relatively noisy due to frequent fluctuations. We can apply a moving average filter to reduce noise and potentially increase the correlation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41318f4a-060e-4e71-b6cf-038772c88e5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ===> CODE REQUIRED <===\n",
    "\n",
    "def filter_signals(series, cols, filter_window):\n",
    "    \"\"\"\n",
    "    Write a function that applies a moving average filter to the selected signals for noise reduction\n",
    "\n",
    "    - create a MovingAverageFilter with the window size `filter_window`\n",
    "    - and then filter only the columns `cols` of `series` (hint: you can select specific columns with series[cols])\n",
    "\n",
    "    Documentation: https://unit8co.github.io/darts/generated_api/darts.models.filtering.moving_average_filter.html#darts.models.filtering.moving_average_filter.MovingAverageFilter\n",
    "    \"\"\"\n",
    "\n",
    "    # < --- START OF YOUR CODE --- >\n",
    "    \n",
    "    # create the filter with a predifined window size\n",
    "    ma_filter = # your code here\n",
    "    cols_filtered = # your code here\n",
    "    \n",
    "    # < --- END OF YOUR CODE --- >\n",
    "    \n",
    "    # the filter renames the columns so we have to revert that\n",
    "    cols_other = series.columns.drop(cols).tolist()\n",
    "    renamed_cols = [f\"rolling_mean_{filter_window}_\" + col for col in filter_cols]\n",
    "    cols_filtered = cols_filtered.with_columns_renamed(renamed_cols, cols)\n",
    "    \n",
    "    # return a new series with the filtered signals\n",
    "    return series[cols_other].stack(cols_filtered)\n",
    "\n",
    "# filter the noisy features...\n",
    "filter_cols = [\n",
    "    \"Gen_RPM_Max\", \n",
    "    \"Prod_LatestAvg_ActPwrGen0\",\n",
    "    \"Prod_LatestAvg_ActPwrGen1\",\n",
    "]\n",
    "# ... with a 2 hour window = 2 * 6 (10 minute steps) = 12\n",
    "filter_window = 12\n",
    "\n",
    "train_raw_filt = filter_signals(train_raw, filter_cols)\n",
    "test_raw_filt = filter_signals(test_raw, filter_cols)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32ada528-1122-49f0-9e32-7e74acda7dc3",
   "metadata": {},
   "source": [
    "And we plot it again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f1f6fcdf-6591-45ef-b2ca-6fa8b7e416e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ===> CODE REQUIRED <===\n",
    "\n",
    "\"\"\"\n",
    "Plot the feature analysis again, this time with the filtered series, and filtered columns\n",
    "\"\"\"\n",
    "\n",
    "plot_feature_analysis(\n",
    "    # your code\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da2e1ad2-9010-440c-8df0-9d4da150e50b",
   "metadata": {},
   "source": [
    "Indeed, it did increase the correlation.\n",
    "\n",
    "Great, now we are ready to start with preparing the data for the modelling!\n",
    "\n",
    "## Dataset Preparation for Training\n",
    "\n",
    "### Anomaly Removal\n",
    "We want to estimate the normal operating range of the generator temperature. For this, we must train the model on an anomaly-free data. This can be done in different ways:\n",
    "\n",
    "- Split TimeSeries into several series that do not contain anomaly time frames\n",
    "- Generate a binary anomaly flag, mark anomalous time frames in training TimeSeries as `True`, and at prediction time as `False`. Use this feature as model input.\n",
    "\n",
    "Let's go with the latter option. Let's create a function that \"removes\" the known anomalies, including the preceding and following 7 days (2 week buffer window). \n",
    "\n",
    "> We also remove the non-generator-specific anomalies, since we have no guarantee that the generator was not affected by them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "579893d1-f751-4aaa-80a6-9700e17eef5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "\n",
    "def remove_anomalies(\n",
    "    series: TimeSeries,\n",
    "    tg_col: str,\n",
    "    anomalies_df: pd.DataFrame,\n",
    "    turbine_id: str,\n",
    "):\n",
    "    \"\"\"\n",
    "    Returns a new TimeSeries with removed anomalies (including the last and following 7 days). \n",
    "    Removed anomalies meaning:\n",
    "\n",
    "    - all covariate values are replaced with `np.nan`\n",
    "    - flag \"missing_date\" is set to `True`\n",
    "    - target values are set to `0.`\n",
    "    \"\"\"\n",
    "    df = series.pd_dataframe(copy=True)\n",
    "    one_week = FREQ * 6 * 24 * 7\n",
    "\n",
    "    # mask all anomalous time frames\n",
    "    anom_mask = np.zeros(len(series)).astype(bool)\n",
    "    anom_turb = anomalies_df[anomalies_df[\"Turbine_ID\"] == turbine_id].reset_index(drop=True)\n",
    "    for start in anom_turb[\"start\"]:\n",
    "        anom_mask = anom_mask | ((series.time_index >= start - one_week) & (series.time_index <= start + one_week))\n",
    "\n",
    "    df.loc[anom_mask, :] = np.nan\n",
    "    df.loc[anom_mask, tg_col] = 0.\n",
    "    df.loc[anom_mask, \"missing_date\"] = 1.\n",
    "    return TimeSeries.from_dataframe(df)\n",
    "\n",
    "# data cleaning (remove anomalies)\n",
    "train = remove_anomalies(\n",
    "    series=train_raw_filt,\n",
    "    tg_col=tg_col,\n",
    "    anomalies_df=df_anom_true,\n",
    "    turbine_id=turbine_id_train,\n",
    ")\n",
    "test = remove_anomalies(\n",
    "    series=test_raw_filt,\n",
    "    tg_col=tg_col,\n",
    "    anomalies_df=df_anom_true,\n",
    "    turbine_id=turbine_id_test,\n",
    ")\n",
    "\n",
    "plot_series(train, test, anom_true)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a7cdd97-e94b-4461-a19d-278824d6e782",
   "metadata": {},
   "source": [
    "### Target and Covariates Extraction for Training and Validation Set"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9c9aa62-61d7-48a8-9d23-ea762c96b191",
   "metadata": {},
   "source": [
    "Now we extract the target series and covariates to train the model.\n",
    "We'll use `CatBoostModel` which support a validation set for early stopping to reduce overfitting.\n",
    "\n",
    "We use the entire year 2016 of the anomaly-free Turbine 6 data for training. For the validation set, we use the first 10k points of Turbine 7 (roughly two months)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05000456-8c4a-48cd-8f11-86c9a6d0acd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ===> CODE REQUIRED <===\n",
    "\n",
    "\"\"\"Extract the target and covariates for both the train and thest turbine\n",
    "\n",
    "Hint: you can extract columns from a series with `series[list_of_columns]`\n",
    "\"\"\"\n",
    "\n",
    "# extract target variable from `train` and `test`\n",
    "tg_train = #\n",
    "tg_test = #\n",
    "\n",
    "# extract covariates from `train` and `test`\n",
    "covs_train = #\n",
    "covs_test = #\n",
    "\n",
    "\n",
    "\"\"\"Extract the first 10'000 values from the test target and covariates, that we will use ase a validation set.\n",
    "\n",
    "Hint: you can slice a series with `series[start_idx:end_idx]`\n",
    "\"\"\"\n",
    "\n",
    "# for early stopping extract values from `tg_test` ad `covs_test`\n",
    "tg_val = #\n",
    "covs_val = # "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26af3745-4d1d-4e33-bc36-612aefe12ca0",
   "metadata": {},
   "source": [
    "## Model Setup And Training\n",
    "\n",
    "> **Reminder:** Because our predictive maintenance algorithm is done post hoc, we can model the generator temperature at all times `t` with actual measurement data from other sensors at times `t` (in regular forecasting, we wouldn't know these values as they lie in the future).\n",
    "\n",
    "\n",
    "To perform this task, we set up the model as follows:\n",
    "\n",
    "- We train the model to predict only one step `t` at a time -> `output_chunk_length=1`\n",
    "- We train a probabilistic model to forecast temperature quantiles (0.05, 0.50, 0.95) which we will use a proxy for the normal range (`likelihood=\"quantile\", quantiles=[0.05, 0.50, 0.95]`)\n",
    "- We **don't use** the past of the temperature as a model input (because of the causal flow) -> `lags=None`\n",
    "- We only use external features (future covariates) because we can use information at time `t` to predict the temperature at time `t` -> `lags_future_covariates`\n",
    "- Additionally we add some calendar information (the month and hour of the day) as predictors for some seasonal patterns (e.g. the month for the different seasons, hour for day and night) -> `add_encoders`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0131d40c-bb20-4971-8720-e647ae561336",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# ===> CODE REQUIRED <===\n",
    "\n",
    "\"\"\"\n",
    "We set up the CatBoostModel for the predictive maintenance task:\n",
    "All the parameters you have to fill in are described in the model documentation: \n",
    "- https://unit8co.github.io/darts/generated_api/darts.models.forecasting.catboost_model.html#darts.models.forecasting.catboost_model.CatBoostModel\n",
    "\n",
    "1) add calendar information as input features to the model. These will be generated automatically by the model.\n",
    "\n",
    "Hint: add a cyclic enconding of \"month\", and \"hour\" as future covariates.\n",
    "\"\"\"\n",
    "add_encoders = {\n",
    "    \"cyclic\": {\n",
    "        \"future\": [\"month\", \"hour\"]\n",
    "    }\n",
    "}\n",
    "\n",
    "\"\"\"\n",
    "2) make the model not use any information from the target variable\n",
    "\"\"\"\n",
    "\n",
    "lags = # your code here\n",
    "\n",
    "\"\"\"\n",
    "3) for the future covariates, define component specific lags:\n",
    "- Specify a default lag with key \"default_lags\". It should use 24 steps from the past, and 1 step from the future\n",
    "- Additionally, add specific lags for the calendar month features:\n",
    "  - \"darts_enc_fc_cyc_month_sin\": use only 1 step from the future\n",
    "  - \"darts_enc_fc_cyc_month_cos\": use only 1 step from the future\n",
    "\"\"\"\n",
    "\n",
    "lags_future_covariates = {\n",
    "    # your code here\n",
    "}\n",
    "\n",
    "\"\"\"\n",
    "4) make the model probabilistic with Quantile Regression (likelihood)\n",
    "And give the quantiles (0.05, 0.50, 0.95) that the model should predict\n",
    "\"\"\"\n",
    "\n",
    "# define the quantile to be fitted\n",
    "likelihood = # use quantile likelihood\n",
    "quantiles = # give a list of three quantiles\n",
    "\n",
    "# create the model\n",
    "from darts.models import CatBoostModel\n",
    "model = CatBoostModel(\n",
    "    lags=None,\n",
    "    lags_future_covariates=lags_future_covariates,\n",
    "    output_chunk_length=1,\n",
    "    add_encoders=add_encoders,\n",
    "    likelihood=\"quantile\",\n",
    "    quantiles=quantiles,\n",
    "    random_state=42,\n",
    ")\n",
    "\n",
    "\"\"\"\n",
    "Train the model with a validation set to prevent overfitting.\n",
    "You have to use `tg_train`, `covs_train`, `tg_val`, `covs_val`\n",
    "\n",
    "Documentation: https://unit8co.github.io/darts/generated_api/darts.models.forecasting.catboost_model.html#darts.models.forecasting.catboost_model.CatBoostModel.fit\n",
    "\"\"\"\n",
    "model.fit(\n",
    "    # your code here\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "140249a3-37d1-433f-b9e5-0014ae543760",
   "metadata": {},
   "source": [
    "## Model Evaluation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "baa74b77-fcbd-4acd-aecf-d360b5d4e99f",
   "metadata": {},
   "source": [
    "### Anomaly-Free Turbine 6 Forecasts\n",
    "\n",
    "Now we should check how the model performed. First, we take a look at the anomaly-free (cleaned) Turbine 6 data from 2016.\n",
    "Ideally, if we let our model forecast the entire year, we shouldn't detect any major anomalous behavior.\n",
    "\n",
    "We generate the historical forecasts over the year 2016 and show the predictions for the entire year, and the first couple of days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9636325d-de72-41a6-8baf-3887b516a5f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ===> CODE REQUIRED <===\n",
    "\n",
    "# < --- START OF YOUR CODE --- >\n",
    "\n",
    "\"\"\"\n",
    "Define the historical forecast parameters. It should:\n",
    "- Perform a one step forecast (forecast_horizon)\n",
    "- Move by one step for the next prediction (stride)\n",
    "- return only the last points per prediction (last_points_only)\n",
    "- predict the quantile \"likelihood\" parameters directly (predict_likelihood_parameters)\n",
    "\"\"\"\n",
    "\n",
    "hist_fc_params = {\n",
    "    \"forecast_horizon\": # your code here,\n",
    "    \"stride\":,\n",
    "    \"last_points_only\":,\n",
    "    \"retrain\":,\n",
    "    \"predict_likelihood_parameters\":,\n",
    "}\n",
    "\n",
    "\"\"\"\n",
    "Generate the historical forecasts using series `tg_train`, future covariates `covs_train`, \n",
    "and all parameters from `hist_fc_params`.\n",
    "\n",
    "Hint: you can feed the parameters from a dict to a function like this:\n",
    "`some_function(**hist_fc_params)`\n",
    "\n",
    "The historical forecast documentation is here: https://unit8co.github.io/darts/generated_api/darts.models.forecasting.catboost_model.html#darts.models.forecasting.catboost_model.CatBoostModel.historical_forecasts\n",
    "\"\"\"\n",
    "\n",
    "hist_fc_train = # your code here\n",
    "\n",
    "# < --- END OF YOUR CODE --- >\n",
    "\n",
    "\n",
    "# we added a helper function `plot_intervals` to plot the quantile predictions as intervals\n",
    "\n",
    "def plot_hist_fcs(series, hist_fc):\n",
    "    fig, (ax1, ax2) = plt.subplots(nrows=2, figsize=(12, 6))\n",
    "    series.plot(ax=ax1)\n",
    "    plot_intervals(hist_fc, ax=ax1)\n",
    "    ax1.set_title(\"Entire Year 2016 Forecasts\")\n",
    "    \n",
    "    series[times].plot(ax=ax2)\n",
    "    plot_intervals(hist_fc[times], ax=ax2)\n",
    "    ax2.set_title(\"First Days 2016 Forecasts\")\n",
    "    \n",
    "    fig.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "plot_hist_fcs(tg_train, hist_fc_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c6ce6c3-5b1c-45fa-867c-d944424764ff",
   "metadata": {},
   "source": [
    "The graph shows the actual temperature in black, predicted median temperature in blue (quantile 0.50), and the estimated normal operating range as a light blue interval (upper bound = quantile 0.95, lower bound = quantile 0.05).\n",
    "\n",
    "These predictions look rather nice, but it's also the data that the model has been trained on. Nevertheless, let's continue the investigation.\n",
    "\n",
    "Let's also remove predictions at times when the flag \"missing_date\" was `True`, since these could also not be forecasted properly in the real-world scenario."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "612f47bb-3b88-4562-9093-7f7ee86016eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "def postprocess_preds(series: TimeSeries, covs: TimeSeries):\n",
    "    \"\"\"Return a new series, where the values have bin replaced with \n",
    "    `np.nan` where flag 'is_missing' is `True`.\n",
    "    \"\"\"\n",
    "    vals = series.values(copy=True)\n",
    "    is_missing = covs[\"missing_date\"].slice_intersect_values(series, copy=False)[:, 0, 0]\n",
    "    vals[is_missing == 1.] = np.nan\n",
    "    return series.with_values(vals)\n",
    "\n",
    "hist_fc_train = postprocess_preds(hist_fc_train, covs_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba47db98-a26d-45ca-b77c-615977842612",
   "metadata": {},
   "source": [
    "### Residuals Computation\n",
    "What we are interested in is: how often and by how much was the actual temperature outside the predicted normal operating range?\n",
    "\n",
    "For this, we can compute the residuals per predicted quantile."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "661d6689-b339-4814-b255-90ad1ed387c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---> No Code Required <---\n",
    "# repeat target column so when can compute the residuals per quantile\n",
    "from darts import concatenate\n",
    "train_extended = concatenate([tg_train] * 3, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b7c508f9-1c29-4a25-9cdb-5bbe2b52c205",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ===> CODE REQUIRED <===\n",
    "\n",
    "\"\"\"\n",
    "Compute the residuals between `train_extended` and `hist_fc_train`. \n",
    "\n",
    "Documentation: https://unit8co.github.io/darts/generated_api/darts.models.forecasting.catboost_model.html#darts.models.forecasting.catboost_model.CatBoostModel.residuals\n",
    "\"\"\"\n",
    "\n",
    "residuals = # you code here\n",
    "\n",
    "ax = residuals[times].plot(lw=1)\n",
    "ax.set_title(\"Residuals per quantile for the first days\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb002734-6549-418a-a34b-7c7fb123811a",
   "metadata": {},
   "source": [
    "The residuals are computed as `y_true - y_pred`, so we can find out when the actual temperature was above or below the interval when:\n",
    "\n",
    "- actual temperature is below if `q0.05 > 0`\n",
    "- actual temperature is above if `q0.95 < 0`\n",
    "\n",
    "From now on, we only focus on the times when one of these two conditions was met.\n",
    "\n",
    "Also this was the last code cell which requires your input. From here on we'll go through the rest together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4ac4d2c-db40-4faf-b653-bc82843c373f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = residuals.pd_dataframe(copy=True)\n",
    "\n",
    "# rename the residuals to something shorter\n",
    "dev_low, def_med, def_high, dev_agg = [\"dev_low\", \"def_med\", \"def_high\", \"dev_agg\"]\n",
    "df.columns = [dev_low, def_med, def_high]\n",
    "\n",
    "# let's drop the median prediction\n",
    "df = df[[dev_low, def_high]]\n",
    "\n",
    "# fill missing dates predictions with zeros.\n",
    "df = df.fillna(0.)\n",
    "\n",
    "# ignore residuals where y_true was within high and low quantile\n",
    "# points below the lower bound\n",
    "df.loc[df[dev_low] > 0., dev_low] = 0.\n",
    "# points above the upper bound\n",
    "df.loc[df[def_high] < 0., def_high] = 0.\n",
    "# points above or below the interval\n",
    "df.loc[:, dev_agg] = df[dev_low] + df[def_high]\n",
    "\n",
    "# compute the share of anomalous time steps for each condition\n",
    "(df != 0).mean() * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6943ace9-c0c0-44d1-bce7-f71565cc10bb",
   "metadata": {},
   "source": [
    "So, 3.5% of actual temperatures were above the interval, 3.8% were below, giving a total 7.2% out of bounds.\n",
    "\n",
    "This looks fine, as our quantile interval should include roughly 90% of all actual values.\n",
    "\n",
    "Let's quickly look at the residual analysis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d8360be-97d5-40d9-a478-f24cc6aa223c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# replace 0. values with np.nan to ignore them for the residuals analysis\n",
    "df_ = df[[dev_agg]].copy()\n",
    "df_[df_ == 0.] = np.nan\n",
    "\n",
    "from darts.utils.statistics import plot_residuals_analysis\n",
    "plot_residuals_analysis(TimeSeries.from_dataframe(df_), num_bins=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d1fca84-92da-4c2f-9a36-ffdd49a1912b",
   "metadata": {},
   "source": [
    "From January until November 2016, the residuals seem rather stable. Interestingly, from November onwards, the residuals tend to be much larger. We either miss some important feature/measurement that can explain these deviations, or the temperature measurements in this time frame are actually showing an anomalous behavior.\n",
    "\n",
    "> **Tip:** In such cases it is strongly recommended to communicate with the engineers/operators/... to find an explanation for it.\n",
    "\n",
    "### Anomaly Detection\n",
    "\n",
    "Okay, now we can start defining what we consider as an anomaly. A single point being slightly outside the interval can always happen by chance. Instead, let's try to find some significant anomalies. For this, we can apply the following:\n",
    "- we could ignore some points which were only slightly outside the interval (threshold `min_value`)\n",
    "- we scan our residuals in fixed size windows and count how many points were out-of bounds in each window (`anom_window`)\n",
    "- for each window we can set a minimum out-of-bounds probability below which we do not consider the window as anomalous (`min_anom_prob`)\n",
    "\n",
    "Since this is a non-trivial operation, we provide a helper function `compute_anomalies()` for it.\n",
    "The function performs all steps from above, including the residuals computation. It returns\n",
    "\n",
    "- `anom_pred`: the residuals (\"dev_lo\", ...) and final anomaly flags (\"anom_lo\", ...) as a `TimeSeries`.\n",
    "- `df_anom_pred`: a `DataFrame` where each row represents an anomaly, including the start, end date and some statistics\n",
    "- `ql`: the quantile loss for the historical forecast (as a metric of how good the quantile predictions are; with zero being the best score)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf229bb0-b87b-4ee0-8f6b-f4ad0760756a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# anomaly detection settings\n",
    "min_value = 2.  # minimum 2 degrees outside\n",
    "anom_window = 6*24  # a one day window\n",
    "min_anom_prob = 0.25  # minimum 25% of points should be out of bounds in a window\n",
    "\n",
    "anom_train_pred, df_anom_train_pred, ql_train = compute_anomalies(\n",
    "    model=model,\n",
    "    series=tg_train,\n",
    "    pred_series=hist_fc_train,\n",
    "    quantiles=quantiles,\n",
    "    min_value=min_value,\n",
    "    anom_window=anom_window,\n",
    "    min_anom_prob=min_anom_prob,\n",
    ")\n",
    "anom_train_pred[[\"anom_lo\", \"anom_hi\", \"anom_agg\"]].pd_dataframe().mean(axis=0) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a28f2ac-3a1c-4464-91b7-0cbfa3043c95",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_anom_train_pred = df_anom_train_pred[df_anom_train_pred[\"name\"] == \"anom_agg\"]\n",
    "df_anom_train_pred = df_anom_train_pred.sort_values(by=\"dev_tot\", ascending=False).reset_index(drop=True)\n",
    "df_anom_train_pred"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e062edb1-a310-479e-96f8-79af739981ad",
   "metadata": {},
   "source": [
    "With this setting, we only marked a single time period as anomalous. It also lies at the end of 2016 where we identified the larger residual values before.\n",
    "\n",
    "> The `min_anom_prob` value can be calibrated by comparing the predicted anomalies on the anomaly-free series (should give roughly zero anomalies) with the predicted anomalies on the anomalous series (should find roughly as many anomalies as there are true anomalies)\n",
    "\n",
    "And now we visualize the predicted anomaly with helper function `plot_predicted_anomalies()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9bd86ae5-2320-4c78-b10e-693f646c75f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_predicted_anomalies(\n",
    "    df_anom_train_pred,\n",
    "    tg_train,\n",
    "    covs_train,\n",
    "    hist_fc_train,\n",
    "    anom_true,\n",
    "    anom_train_pred,\n",
    "    turbine_id_train,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e0f5d0e-deca-4e56-8509-70e99e4001bf",
   "metadata": {},
   "source": [
    "This graph shows us:\n",
    "- Top: actual temperature (black), predicted normal operating range (blue), predicted anomaly window (green capsule), residual value if actual temperature was out of bounds (green), actual anomaly (red, only if there was an actual anomaly)\n",
    "- Bottom: the covariates scaled to a value range (0, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5725b3df-4cb3-4ffa-8bf4-61d9c9ce8886",
   "metadata": {},
   "source": [
    "### Turbine 6 Forecasts On Anomalous Data\n",
    "\n",
    "All in all, our model seems properly configured now on the training data with the removed anomalies. \n",
    "\n",
    "Now, let's see how the model performs when we add back the anomalous time frames into the Turbine 6 data.\n",
    "\n",
    "For simplicity, we make use of our helper functions directly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b79f2acf-730c-47c6-aebb-732babc4273c",
   "metadata": {},
   "outputs": [],
   "source": [
    "tg_train_anom = train_raw_filt[tg_col]\n",
    "covs_train_anom = train_raw_filt[cov_cols]\n",
    "\n",
    "hist_fc_train_anom = model.historical_forecasts(\n",
    "    series=tg_train_anom, future_covariates=covs_train_anom, **hist_fc_params\n",
    ")\n",
    "\n",
    "hist_fc_train_anom = postprocess_preds(hist_fc_train_anom, covs_train_anom)\n",
    "anom_train_anom_pred, df_anom_train_anom_pred, ql_train_anom = compute_anomalies(\n",
    "    model=model,\n",
    "    series=tg_train_anom,\n",
    "    pred_series=hist_fc_train_anom,\n",
    "    quantiles=quantiles,\n",
    "    min_value=min_value,\n",
    "    anom_window=anom_window,\n",
    "    min_anom_prob=min_anom_prob,\n",
    ")\n",
    "\n",
    "df_anom_train_anom_pred = df_anom_train_anom_pred[df_anom_train_anom_pred[\"name\"] == \"anom_agg\"]\n",
    "df_anom_train_anom_pred = df_anom_train_anom_pred.sort_values(by=\"dev_tot\", ascending=False).reset_index(drop=True)\n",
    "df_anom_train_anom_pred.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c5deab5-9cc8-40b4-b244-568bcd892203",
   "metadata": {},
   "source": [
    "It looks like now we predict more anomalies with much larger total deviation. Let's check them visually:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3272b1b-14d5-4d30-951e-fd0bde75ea8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_predicted_anomalies(\n",
    "    df_anom_train_anom_pred,\n",
    "    tg_train_anom,\n",
    "    covs_train_anom,\n",
    "    hist_fc_train_anom,\n",
    "    anom_true_gen,\n",
    "    anom_train_anom_pred,\n",
    "    turbine_id_train,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d91b59d-a365-495b-8439-fe0dbbfb5bf0",
   "metadata": {},
   "source": [
    "This looks great, as we can see actual anomalies (red capsules) in almost all the examples! And now plot actual vs predict anomalies throughout 2016."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dcf9c9e8-ed9a-4dca-922c-41ab7e1b2013",
   "metadata": {},
   "outputs": [],
   "source": [
    "anom_true_gen[turbine_id_train].plot(label=\"actual anomalies\")\n",
    "anom_train_anom_pred[\"anom_agg\"].plot(alpha=0.5, label=\"predicted anomalies\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "242f1489-5395-4715-b796-a79d404ab4b1",
   "metadata": {},
   "source": [
    "The model correctly detected all anomalies except number 3 from the left! The last predicted anomaly is the one that we saw already on the anomaly-free input data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80442529-381b-4595-adf6-215995a6bb6b",
   "metadata": {},
   "source": [
    "### Turbine 7 Forecasts On Anomalous Data\n",
    "\n",
    "As a last step, let's see how the model performs in detecting anomalies on turbine 7.\n",
    "Remember, this **model has not been trained on this turbine at all** (only the first 2 months were used for early stopping, not learning)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05813114-0160-46e5-adec-8dcfab6bbce6",
   "metadata": {},
   "outputs": [],
   "source": [
    "tg_test_anom = test_raw_filt[tg_col]\n",
    "covs_test_anom = test_raw_filt[cov_cols]\n",
    "\n",
    "hist_fc_test_anom = model.historical_forecasts(\n",
    "    series=tg_test_anom, future_covariates=covs_test_anom, **hist_fc_params\n",
    ")\n",
    "\n",
    "hist_fc_test_anom = postprocess_preds(hist_fc_test_anom, covs_test_anom)\n",
    "anom_test_anom_pred, df_anom_test_anom_pred, ql_test_anom = compute_anomalies(\n",
    "    model=model,\n",
    "    series=tg_test_anom,\n",
    "    pred_series=hist_fc_test_anom,\n",
    "    quantiles=quantiles,\n",
    "    min_value=min_value,\n",
    "    anom_window=anom_window,\n",
    "    min_anom_prob=min_anom_prob,\n",
    ")\n",
    "\n",
    "df_anom_test_anom_pred = df_anom_test_anom_pred[df_anom_test_anom_pred[\"name\"] == \"anom_agg\"]\n",
    "df_anom_test_anom_pred = df_anom_test_anom_pred.sort_values(by=\"dev_tot\", ascending=False).reset_index(drop=True)\n",
    "df_anom_test_anom_pred.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ed22eea-3a81-460c-a4aa-3cca9bb8adf0",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_predicted_anomalies(\n",
    "    df_anom_test_anom_pred,\n",
    "    tg_test_anom,\n",
    "    covs_test_anom,\n",
    "    hist_fc_test_anom,\n",
    "    anom_true_gen,\n",
    "    anom_test_anom_pred,\n",
    "    turbine_id_test,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8709efa-8803-4eff-b6e3-77ddce6a51a5",
   "metadata": {},
   "source": [
    "The most significant predicted anomaly (the first one in order) is also an actual anomaly (and the only actual generator anomaly).\n",
    "All other predicted anomalies have much lower total deviation (residuals). This is an excellent result, considering that the model has not been trained on this turbine."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55bbd177-de17-4e7e-be16-0e5ca71e2c9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "anom_true_gen[turbine_id_test].plot(label=\"actual anomalies\")\n",
    "anom_test_anom_pred[\"anom_agg\"].plot(alpha=0.5, label=\"predicted anomalies\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4d39b96-f5eb-4f45-87e4-c0cb5d2b1528",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "\n",
    "We showed how to perform a predictive maintenance task using Darts on the example of Wind Turbine Failures. We looked at:\n",
    "- which features to use (and not use) for modelling\n",
    "- how to setup and train a model to predict the expected normal operating range of a signal\n",
    "- how to detect anomalies from deviations between the actual signal values and the predicted operating range\n",
    "\n",
    "Since this was an open-source dataset, we did not have the chance to communicate with the engineers/operators of the turbines.\n",
    "To further the results, it would be crucial to further investigate with them:\n",
    "- if there are other relevant features for modelling the temperature signal\n",
    "- if there is an explanation why the residuals are higher at the end of the year\n",
    "- ...\n",
    "\n",
    "And in general, to improve the results: We only modeled one temperature signal in this example. For a real predictive maintenance pipeline, we would:\n",
    "- have dedicated models for different signals\n",
    "- perform the anomaly detection for each signal\n",
    "- from these results, we can draw conclusions on origin of anomaly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ac5b26f-b6f9-448a-95d8-5ddf330258fa",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.13"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}