{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "48e48d7f-053e-4a1b-8f95-f2c083c949cd",
      "metadata": {},
      "source": [
        "# Tutorial 02: TENAX vs SMEV (extreme value method)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "24266dc5-298b-47a7-8884-6b391dd78400",
      "metadata": {},
      "source": [
        "Import all libraries needed for running TENAX and SMEV"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "9fb356f4-980b-473d-bdc7-1d4b0a3d8d3a",
      "metadata": {},
      "outputs": [],
      "source": [
        "from importlib.resources import files\n",
        "import numpy as np\n",
        "import pandas as pd\n",
        "import matplotlib.pyplot as plt\n",
        "# Import pyTENAX\n",
        "from pyTENAX import tenax, smev, plotting"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d075c82d-500c-46a9-ab7c-bdd65283d44b",
      "metadata": {},
      "source": [
        "Let's initiate TENAX class and SMEV class with given setup."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "2f471a9b-98b3-450e-b326-60ad6117cd08",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Initiate TENAX class with customized setup\n",
        "S = tenax.TENAX(\n",
        "    return_period=[\n",
        "        2,\n",
        "        5,\n",
        "        10,\n",
        "        20,\n",
        "        50,\n",
        "        100,\n",
        "        200,\n",
        "    ],\n",
        "    durations=[10, 60, 180, 360, 720, 1440], #durations are in minutes and they refer to depth of rainfall within given duration\n",
        "    time_resolution=10,  # time resolution in minutes\n",
        "    left_censoring=[0, 0.90], # left censoring threshold \n",
        "    alpha=0.05, #dependence of shape on T depends on statistical significance at the alpha-level.\n",
        "    min_rain = 0.1, #minimum rainfall depth threshold\n",
        ")\n",
        "\n",
        "# Initiate SMEV class with customized setup following TENAX \n",
        "S_SMEV = smev.SMEV(\n",
        "    return_period=S.return_period,\n",
        "    durations=S.durations,\n",
        "    time_resolution=5,  # time resolution in minutes\n",
        "    min_rain=0.1,\n",
        "    storm_separation_time=24,\n",
        "    min_event_duration=30,\n",
        "    left_censoring=[S.left_censoring[1], 1],\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "1641d7d1-9551-4417-b039-c9b1d10a5660",
      "metadata": {},
      "source": [
        "Load same test data as in Tutorial 01."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b43a02dd-66fc-4007-ba4d-f16f10247962",
      "metadata": {},
      "outputs": [],
      "source": [
        "# Load precipitation data\n",
        "# Create input path file for the test file\n",
        "file_path_input = files('pyTENAX.res').joinpath('prec_data_Aadorf.parquet')\n",
        "# Load data from csv file\n",
        "data = pd.read_parquet(file_path_input)\n",
        "# Convert 'prec_time' column to datetime, if it's not already\n",
        "data[\"prec_time\"] = pd.to_datetime(data[\"prec_time\"])\n",
        "# Set 'prec_time' as the index\n",
        "data.set_index(\"prec_time\", inplace=True)\n",
        "name_col = \"prec_values\"  # name of column containing data to extract\n",
        "\n",
        "# load temperature data\n",
        "file_path_temperature = files('pyTENAX.res').joinpath('temp_data_Aadorf.parquet')\n",
        "t_data = pd.read_parquet(file_path_temperature)\n",
        "# Convert 'temp_time' column to datetime if it's not already in datetime format\n",
        "t_data[\"temp_time\"] = pd.to_datetime(t_data[\"temp_time\"])\n",
        "# Set 'temp_time' as the index\n",
        "t_data.set_index(\"temp_time\", inplace=True)\n",
        "temp_name_col = \"temp_values\""
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d716d3c5-4bff-4465-bd59-a97e9d9c242a",
      "metadata": {},
      "source": [
        "## Repeat the preprocessing from Tutorial 01.\n",
        "We once again focus only on 10-minute rainfall depth."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f4240b96-48dd-4767-ac21-644fd9e82fe3",
      "metadata": {},
      "outputs": [],
      "source": [
        "data = S.remove_incomplete_years(data, name_col)\n",
        "\n",
        "# get data from pandas to numpy array\n",
        "df_arr = np.array(data[name_col])\n",
        "df_dates = np.array(data.index)\n",
        "df_arr_t_data = np.array(t_data[temp_name_col])\n",
        "df_dates_t_data = np.array(t_data.index)\n",
        "\n",
        "# extract indexes of ordinary events\n",
        "# these are time-wise indexes =>returns list of np arrays with np.timeindex\n",
        "idx_ordinary = S.get_ordinary_events(data=df_arr, \n",
        "                                     dates=df_dates, \n",
        "                                     name_col=name_col,\n",
        "                                     check_gaps=False)\n",
        "\n",
        "# get ordinary events by removing too short events\n",
        "# returns boolean array, dates of OE in TO, FROM format, and count of OE in each years\n",
        "arr_vals, arr_dates, n_ordinary_per_year = S.remove_short(idx_ordinary)\n",
        "\n",
        "# assign ordinary events values by given durations, values are in depth per duration, NOT in intensity mm/h\n",
        "dict_ordinary, dict_AMS = S.get_ordinary_events_values(data=df_arr, \n",
        "                                                       dates=df_dates, \n",
        "                                                       arr_dates_oe=arr_dates)\n",
        "\n",
        "dict_ordinary, _, n_ordinary_per_year = S.associate_vars(dict_ordinary, \n",
        "                                                         df_arr_t_data, \n",
        "                                                         df_dates_t_data)\n",
        "\n",
        "# Your data (P, T arrays) and threshold thr=3.8\n",
        "P = dict_ordinary[\"10\"][\"ordinary\"].to_numpy()  # Replace with your actual data\n",
        "T = dict_ordinary[\"10\"][\"T\"].to_numpy()  # Replace with your actual data\n",
        "blocks_id = dict_ordinary[\"10\"][\"year\"].to_numpy()  # Replace with your actual data\n",
        "# Number of threshold\n",
        "thr = dict_ordinary[\"10\"][\"ordinary\"].quantile(S.left_censoring[1])"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "bdf70a0b-a461-4201-a6df-8378918fd4ac",
      "metadata": {},
      "source": [
        "## Repeat TENAX model (same as Tutorial 01)  \n",
        "Note: This doesn't run bootstrapping for TENAX uncertainty."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "84b1fa4f-abd1-4655-8026-f4625d9e8b11",
      "metadata": {},
      "outputs": [],
      "source": [
        "eT = np.arange(\n",
        "    np.min(T)-4, np.max(T) + 4, 1\n",
        ")  # define T values to calculate distributions. +4 to go beyond graph end\n",
        "\n",
        "# magnitude model\n",
        "F_phat, loglik, _, _ = S.magnitude_model(P, T, thr)\n",
        "\n",
        "# temperature model\n",
        "g_phat = S.temperature_model(T)\n",
        "\n",
        "# Sampling intervals for the Montecarlo\n",
        "Ts = np.arange(np.min(T) - S.temp_delta, \n",
        "               np.max(T) + S.temp_delta, \n",
        "               S.temp_res_monte_carlo)\n",
        "\n",
        "#  mean n of ordinary events\n",
        "n = n_ordinary_per_year.sum() / len(n_ordinary_per_year)\n",
        "\n",
        "# estimates return levels using MC samples\n",
        "RL, _, P_check = S.model_inversion(F_phat, g_phat, n, Ts)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5904f039-9287-4d1b-a814-880b14fb35fb",
      "metadata": {},
      "source": [
        "## SMEV model"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "908256db-c93b-494f-9a55-ab4683b9b41b",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "# estimate shape and  scale parameters of weibull distribution\n",
        "smev_shape, smev_scale = S_SMEV.estimate_smev_parameters(P, S_SMEV.left_censoring)\n",
        "# estimate return period (quantiles) with SMEV\n",
        "smev_RL = S_SMEV.smev_return_values(\n",
        "    S_SMEV.return_period, smev_shape, smev_scale, n.item()\n",
        ")\n",
        "\n",
        "smev_RL_unc = S_SMEV.smev_bootstrap_uncertainty(P, blocks_id, S.niter_smev, n.item())"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ec6dae24-d55d-4f03-b82e-ba175a44ed33",
      "metadata": {},
      "source": [
        "## Plot results TENAX vs SMEV \n",
        "Note, we excluded TENAX uncertainty here as it would take too long to run."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "54fb390c-44e2-45dd-8f78-2ab96299bec2",
      "metadata": {
        "trusted": true
      },
      "outputs": [],
      "source": [
        "AMS = dict_AMS[\"10\"]  # yet the annual maxima\n",
        "plotting.TNX_FIG_valid(AMS,\n",
        "                       S.return_period, \n",
        "                       RL=RL, \n",
        "                       smev_RL=smev_RL, \n",
        "                       RL_unc=[], \n",
        "                       smev_RL_unc=smev_RL_unc)\n",
        "plt.title(\"fig 4\")\n",
        "plt.ylabel(\"10-minute precipitation (mm)\")\n",
        "plt.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.2))\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4f0e577f-98a4-4f92-b92f-27523d2d9c65",
      "metadata": {},
      "source": [
        "### *Plot above represents how well TENAX fits extreme value method which does not use temperature as covariate.*"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python (Pyodide)",
      "language": "python",
      "name": "python"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "python",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.8"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}