{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a3c5787b",
   "metadata": {},
   "source": [
    "Visualization of the results of the model\n",
    "=================================\n",
    "\n",
    "\n",
    "A GAMCR model has two main methods to produce results:\n",
    "\n",
    "* `.predict_streamflow(matJ)`: to get the hydrograph\n",
    "\n",
    "* `.predict_transfer_function(X)`: to get the transfer functions (RRDs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8239dc8a",
   "metadata": {},
   "source": [
    "# 1. Simple visualization\n",
    "\n",
    "Below we provide is simple code to show the observed and predicted hydrographs, and the predicted average NRF transfer function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "614cac2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import GAMCR\n",
    "\n",
    "site= 'flashy'\n",
    "# You should change this line specifying the path of your model\n",
    "model.load_model(os.path.join(root,'{0}/{0}_best_model.pkl'.format(site)))\n",
    "X, matJ, y, timeyear, dates = model.load_data(os.path.join(root, '{0}/data/'.format(site)), max_files=20)\n",
    "\n",
    "yhat = model.predict_streamflow(matJ)\n",
    "plt.figure()\n",
    "plt.plot(yhat, label='estimation')\n",
    "plt.plot(y, label='ground truth', linestyle='--')\n",
    "plt.plot([],[], color='white', label='NSE: {0}'.format(GAMCR.nse(y, yhat)))\n",
    "plt.title(site)\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "a = np.mean(model.gam._modelmat(X), axis=0)\n",
    "smooth_P = np.kron( GAMCR.build_custom_matrix(model.L), np.dot(a.reshape(-1,1), a.reshape(1,-1)) )\n",
    "\n",
    "coeffs = model.gam.get_coeffs().reshape(-1)\n",
    "print( coeffs.T @ smooth_P @ coeffs)\n",
    "\n",
    "plt.figure()\n",
    "H = model.predict_transfer_function(X)\n",
    "plt.figure()\n",
    "plt.plot(np.array([i/24 for i in range(24*3)]), H.mean(axis=0)[:24*3], label='estimation')\n",
    "plt.title(site)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f182dec",
   "metadata": {},
   "source": [
    "# 2. Advanced analysis of the results\n",
    "\n",
    "By launching the script `compute_statistics.py`, you can compute different information on the learned transfer functions, such as:\n",
    "- the global average NRF/RDD/weighted RRD\n",
    "- the average NRF/RRD/weighted RRD over some ensembles (you can stratify either by precipitation intensity, antecendent wetness or by both)\n",
    "- the area, mean, peak and peak lag of the transfer function over different ensembles for the precipitation intensity\n",
    "\n",
    "## Note that if you want to compare GAMCR with ERRA, you should first launch the script ERRA_experiments.R located in the folder ./experiments/ERRA/\n",
    "\n",
    "If ERRA data is available, the ensembles will be automatically set to the ones used by ERRA to allow a fair comparison. Otherwise, we can define yourself the ensembles.\n",
    "\n",
    "For ease of use, we save the global (and stratified) NRF/RDD/weighted RRD in a CSV file `NRF_RRD.csv`, located in the `results` folder of the site under study. In this file, the columns named `NRF`, `RRD`, and `weighted_avg_RRD` correspond to the average over all time points. The columns named `NRF_group_{k}`, `RRD_group_{k}`, and `weighted_avg_RRD_group_{k}` correspond to the average over the time points for the group with index k. To determine the precipitation and/or antecedent wetness ranges corresponding to a group, you should open the file `groups.csv`.\n",
    "\n",
    "### ➡️ Run the Notebook: `figure_real_data_paper.ipynb`\n",
    "\n",
    "A set of tools has already been implemented in GAMCR. By opening the notebook `figures_paper.ipynb` in the `data_and_visualization` folder, you can view the visualizations available using GAMCR."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0770635f",
   "metadata": {},
   "source": [
    "![Screen shot of NRF averaged over different ranges of precipitation intensity.](../../_static/nb_visualization.png)"
   ]
  }
 ],
 "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}