LRP_UNETs_rank_predictors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d545b75d",
   "metadata": {},
   "source": [
    "### Analyse de the performance of the network when applying the LRPcomposite"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1427bd4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Python ≥3.5 is required\n",
    "import sys\n",
    "assert sys.version_info >= (3, 5)\n",
    "\n",
    "# Scikit-Learn ≥0.20 is required\n",
    "import sklearn\n",
    "assert sklearn.__version__ >= '0.20'\n",
    "\n",
    "from sklearn.impute import SimpleImputer\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.compose import ColumnTransformer\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.metrics import confusion_matrix, classification_report, ConfusionMatrixDisplay, precision_score, recall_score, roc_auc_score, roc_curve\n",
    "from sklearn.utils import class_weight\n",
    "\n",
    "# TensorFlow ≥2.0 is required\n",
    "import tensorflow_addons as tfa\n",
    "import tensorflow as tf\n",
    "assert tf.__version__ >= '2.0'\n",
    "\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras import layers, regularizers\n",
    "\n",
    "print(\"Num GPUs Available: \", len(tf.config.list_physical_devices('GPU')))\n",
    "\n",
    "# Common imports\n",
    "import os\n",
    "import glob\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import xarray as xr\n",
    "import dask\n",
    "import datetime\n",
    "import math\n",
    "import pickle\n",
    "import pathlib\n",
    "import hashlib\n",
    "import seaborn\n",
    "dask.config.set({'array.slicing.split_large_chunks': False})\n",
    "\n",
    "# To make this notebook's output stable across runs\n",
    "np.random.seed(42)\n",
    "\n",
    "# Config matplotlib\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.ticker import StrMethodFormatter\n",
    "\n",
    "mpl.rc('axes', labelsize=14)\n",
    "mpl.rc('xtick', labelsize=12)\n",
    "mpl.rc('ytick', labelsize=12)\n",
    "\n",
    "# Dotenv\n",
    "from dotenv import dotenv_values\n",
    "\n",
    "# Custom utils\n",
    "from utils.utils_data import *\n",
    "from utils.utils_ml import *\n",
    "from utils.utils_resnet import *\n",
    "from utils.utils_plot import *\n",
    "from utils.DNN_models import *\n",
    "from utils.rank_relevances import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a538e73-d554-464b-8b56-9f2690110750",
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaml\n",
    "conf = yaml.safe_load(open(\"config.yaml\"))\n",
    "PRECIP_XTRM = 0.95 # Percentile (threshold) for the extremes\n",
    "PRECIP_DATA = 'ERA5-low' "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42e51a71-ad4d-45ac-a82c-cb7179ef0bbe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# load data\n",
    "# test samples\n",
    "dg_train_X = np.array(xr.open_dataarray('tmp/data/dg_train_X.nc'))\n",
    "dg_train_Y_xtrm = xr.open_dataarray('tmp/data/dg_train_Y_xtrm0.95th.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d20c633d-7d43-476c-bb0f-9cb5c12127e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "dg_test_X = np.array(xr.open_dataarray('tmp/data/dg_test_X.nc'))\n",
    "dg_test_Y = np.array(xr.open_dataarray('tmp/data/dg_test_Y.nc'))\n",
    "dg_test_Y_xtrm = xr.open_dataarray('tmp/data/dg_test_Y_xtrm0.95th.nc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "014720a6-68ce-4b01-a281-de8041a801dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# coordinates\n",
    "lons_x = np.load('tmp/data/lons_y.npy')\n",
    "lats_y = np.load('tmp/data/lats_y.npy')\n",
    "# test-times\n",
    "#times = np.arange(np.datetime64('2016-01-01'), np.datetime64('2022-01-01'))\n",
    "#times = pd.to_datetime(times)\n",
    "# train-times\n",
    "times = np.arange(np.datetime64('1979-01-01'), np.datetime64('2006-01-01'))\n",
    "times = pd.to_datetime(times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b35b19f3-566b-423d-a12f-7a711242f558",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_bool= dg_train_Y_xtrm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "df23810f-9b01-423e-95e0-db87fc9f89b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "models = ['UNET1','UNET2','UNET3','UNET4']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d09e07f2-ee3d-4494-b8e7-df5ac80afec8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# load weights\n",
    "lrp_all = []\n",
    "list_df = []\n",
    "    \n",
    "for m_id in models:\n",
    "    print('LRP',m_id)\n",
    "    lrp = np.load(f'tmp/LRP/lrpcomp_train_DNN_{m_id}.npy')\n",
    "    lrp_all.append(lrp)\n",
    "    localrel_avg, localrel_max  = getmap_localrel(lrp, conf['i_shape'], conf['varnames'], y_bool, lats_y, lons_x, times, icrop=3)\n",
    "    rel_varmeans = localrel_avg.mean(axis=(0,1))\n",
    "    df_sortedvars = [conf['varnames'][i] for i in np.argsort(rel_varmeans)[::-1]]\n",
    "    df_sortval = [rel_varmeans[i] for i in np.argsort(rel_varmeans)[::-1]]\n",
    "    data={'Model': m_id,'Variable' : df_sortedvars, 'Values' : df_sortval}\n",
    "    df = pd.DataFrame(data)\n",
    "    list_df.append(df)\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ada545b-329a-4862-9ea6-33542e7fe872",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_all = pd.concat([list_df[1], list_df[3], list_df[2], list_df[0]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b42ffe3-0962-490b-91e6-4e43742e3f4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "dfpivot = df_all.pivot(columns='Model')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b09181f-bd1d-4335-8449-aa2f2677d69f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf857054-e736-4be7-aa2c-c7210871bbfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_all.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2defec4-627f-4403-9481-b3547bd3600b",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(14,6))\n",
    "\n",
    "b= sns.barplot(ax=ax, x = \"Variable\", y = \"Values\", hue = \"Model\", data = df_all)\n",
    "plt.tick_params(labelsize=7)\n",
    "#plt.show()\n",
    "ax.set_ylabel(None)\n",
    "ax.set_xlabel(None)\n",
    "ax.tick_params(axis='x', rotation=90)\n",
    "b.tick_params(labelsize=15)\n",
    "\n",
    "plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.2f}')) \n",
    "# save figure \n",
    "fname = 'Ranking_predictors_test_UNETs'\n",
    "plt.savefig('figures/' + fname + '.pdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f2f0f8d-fda0-49e9-b324-5fb6f6cf7a2d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# select UNET2 as the best network (higher scores)\n",
    "rel = lrp_all[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "819ca964-3279-4d77-8a54-7b2cdd00fbc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "pixelrel_avg, pixelrel_max = getmap_rel(rel, conf['i_shape'], y_bool, lons_x, lats_y, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77619c26-999e-4acc-85af-8041e8c50cae",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_xr_rel(rel, lats_y,lons_x, vnames, fname, cmap='Reds', vmin=None, vmax=None, vcenter=None, plot=True):\n",
    "    \n",
    "    \n",
    "    mx= xr.DataArray(rel, dims=[\"lat\", \"lon\", \"variable\"],\n",
    "                  coords=dict(lat = lats_y, \n",
    "            lon = lons_x, variable= vnames ))\n",
    "    \n",
    "      \n",
    "    if vcenter is None:\n",
    "    \n",
    "        g = mx.plot.pcolormesh(\"lon\", \"lat\", col=\"variable\", col_wrap=6, robust=True, cmap=cmap,\n",
    "        yincrease = False, extend='both',vmin=vmin, vmax=vmax,\n",
    "        figsize=(12, 8),  cbar_kwargs={\"orientation\": \"vertical\", \"shrink\": 0.9, \"aspect\": 50})\n",
    "    \n",
    "    else:\n",
    "        \n",
    "        if mx.min() == 0:\n",
    "            norm = mcolors.TwoSlopeNorm(vcenter=0, vmax=mx.max())\n",
    "        \n",
    "        else:\n",
    "        \n",
    "            norm = mcolors.TwoSlopeNorm(vmin=mx.min(), vcenter=0, vmax=mx.max())\n",
    "        \n",
    "  \n",
    "        g = mx.plot.pcolormesh(\"lon\", \"lat\", col=\"variable\", col_wrap=6, robust=True, cmap=cmap,\n",
    "        yincrease = False, extend='both', norm = norm,\n",
    "        figsize=(14, 14),  cbar_kwargs={\"orientation\": \"vertical\", \"shrink\": 0.9, \"aspect\": 50})\n",
    "        \n",
    "    #figsize=(14, 12)\n",
    "    for ax, title in zip(g.axes.flat, vnames):\n",
    "\n",
    "        world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "        world.boundary.plot(ax=ax, lw=1, color='k')\n",
    "        ax.set_xlim(min(lons_x), max(lons_x))\n",
    "        ax.set_title(title)\n",
    "        ax.set_ylim(min(lats_y), max(lats_y))\n",
    "        \n",
    "    # To control the space\n",
    "    plt.subplots_adjust(right=0.8, wspace=0.1, hspace=0.3)\n",
    "    if plot:\n",
    "        #plt.tight_layout()\n",
    "        plt.savefig('figures/' + fname + '.pdf')\n",
    "    else:\n",
    "        \n",
    "        plt.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84a5f81a-740d-4cbf-92c8-82cd135260df",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_xr_rel(pixelrel_avg, lats_y, lons_x, conf['varnames'], 'relevances_train_UNET2_pixelwise')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a9fbc62-0413-4c3a-aef4-30b614ba4ffa",
   "metadata": {},
   "outputs": [],
   "source": [
    "#pixelrel_avg_u4, pixelrel_max_4 = getmap_rel(lrp_all[3], conf['i_shape'], y_bool, lons_x, lats_y, False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0213825-8832-4a2c-a6bd-ee371ac1628d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#plot_xr_rel(pixelrel_avg_u4, lats_y, lons_x, conf['varnames'], 'relevances_train_UNET4_pixelwise')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47b7f72a-0c7c-46eb-bebe-7f27147fa803",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "63fd5069d213b44bf678585dea6b12cceca9941eaf7f819626cde1f2670de90d"
  },
  "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}