{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quickstart: the wrapper for VS-Lite in R"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Expected time to run through: 5 mins**\n",
    "\n",
    "\n",
    "In this tutorial, we demonstrate how to run the R package of VS-Lite with PyVSL."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import PyVSL\n",
    "import pandas as pd\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the below, we check the wrapper function first, and let's pay attention to the units of the input T and P."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\u001b[0;31mSignature:\u001b[0m\n",
       "\u001b[0mPyVSL\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mVSL_R\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0msyear\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0meyear\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mphi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mT1\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mT2\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m23\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mM1\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.01\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mM2\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.05\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mMmax\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.76\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mMmin\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.01\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0malph\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.093\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mm_th\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m4.886\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mmu_th\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m5.8\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mrootd\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1000\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mM0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0msubstep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mI_0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mI_f\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m12\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mhydroclim\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'P'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m    \u001b[0mRlib_path\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n",
       "\u001b[0;34m\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
       "\u001b[0;31mDocstring:\u001b[0m\n",
       "VS-Lite tree-ring PSM\n",
       "\n",
       "Porting the VSLiteR code (https://github.com/fzhu2e/VSLiteR),\n",
       "which is a fork of the original version by Suz Tolwinski-Ward (https://github.com/suztolwinskiward/VSLiteR)\n",
       "NOTE: Need to install the R library following the below steps:\n",
       "    1. install.packages(\"devtools\")\n",
       "    2. library(devtools)\n",
       "    3. install_github(\"fzhu2e/VSLiteR\")\n",
       "\n",
       "Args:\n",
       "    syear (int): Start year of simulation.\n",
       "    eyear (int): End year of simulation.\n",
       "    phi (float): Latitude of site (in degrees N).\n",
       "    T (array): temperature timeseries in deg C, with length of 12*Nyrs\n",
       "    P (array): precipitation timeseries in mm/month, with length of 12*Nyrs\n",
       "    # original param. in R package: T (12 x Nyrs) Matrix of ordered mean monthly temperatures (in degEes C).\n",
       "    # original param. in R package: P (12 x Nyrs) Matrix of ordered accumulated monthly precipitation (in mm).\n",
       "    T1: Lower temperature threshold for growth to begin (scalar, deg. C).\n",
       "    T2: Upper temperature threshold for growth sensitivity to temp (scalar, deg. C).\n",
       "    M1: Lower moisture threshold for growth to begin (scalar, v.v).\n",
       "    M2: Upper moisture threshold for growth sensitivity to moisture (scalar, v/v).\n",
       "    Mmax: Scalar maximum soil moisture held by the soil (in v/v).\n",
       "    Mmin: Scalar minimum soil moisture (for error-catching) (in v/v).\n",
       "    alph: Scalar runoff parameter 1 (in inverse months).\n",
       "    m_th: Scalar runoff parameter 3 (unitless).\n",
       "    mu_th: Scalar runoff parameter 2 (unitless).\n",
       "    rootd: Scalar root/\"bucket\" depth (in mm).\n",
       "    M0: Initial value for previous month's soil moisture at t = 1 (in v/v).\n",
       "    substep: Use leaky bucket code with sub-monthly time-stepping? (TRUE/FALSE)\n",
       "    I_0: lower bound of integration window (months before January in NH)\n",
       "    I_f: upper bound of integration window (months after January in NH)\n",
       "    hydroclim: Switch; value is either \"P\" (default) or \"M\" depending on whether the second input climate variable\n",
       "        is precipitation, in which case soil moisture isestimated using the Leaky Bucket model of the CPC,\n",
       "        or soil moisture, in which case the inputs are used directly to compute the growth response.\n",
       "\n",
       "Returns:\n",
       "    res (dict): a dictionary with several lists, keys including:\n",
       "        trw: tree ring width\n",
       "        gT, gM, gE, M, potEv, sample.mean.width, sample.std.width\n",
       "\n",
       "Reference:\n",
       "    + The original VSLiteR code by Suz Tolwinski-Ward (https://github.com/suztolwinskiward/VSLiteR)\n",
       "    + Tolwinski-Ward, S.E., M.N. Evans, M.K. Hughes, K.J. Anchukaitis, An efficient forward model of the climate controls\n",
       "        on interannual variation in tree-ring width, Climate Dynamics, doi:10.1007/s00382-010-0945-5, (2011).\n",
       "    + Tolwinski-Ward, S.E., K.J. Anchukaitis, M.N. Evans, Bayesian parameter estimation and interpretation for an\n",
       "        intermediate model of tree-ring width , Climate of the Past, doi:10.5194/cpd-9-615-2013, (2013).\n",
       "    + Tolwinski-Ward, S.E., M.P. Tingley, M.N. Evans, M.K. Hughes, D.W. Nychka, Probabilistic reconstructions of localtemperature and\n",
       "        soil moisture from tree-ring data with potentially time-varying climatic response, Climate Dynamics, doi:10.1007/s00382-014-2139-z, (2014).\n",
       "\u001b[0;31mFile:\u001b[0m      ~/Github/PyVSL/PyVSL/wrapper.py\n",
       "\u001b[0;31mType:\u001b[0m      function\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "PyVSL.VSL_R?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Load test data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First of all, we load the test data (can be downloaded with this [link](https://github.com/fzhu2e/PyVSL/raw/main/docsrc/tutorial/data/test_T_P.pkl)), which contains T and P at a grid point of the CCSM4 past1000 simulation.\n",
    "Note that the unit for T is degree C and that for P is accumulated monthly precipitation (mm/month)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dict = pd.read_pickle('./data/test_T_P.pkl')\n",
    "time = data_dict['time']\n",
    "T = data_dict['T'] # unit: degC\n",
    "P = data_dict['P'] # unit: mm/month"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 850.04109589  850.12328767  850.20273973 ... 1850.7890411  1850.8739726\n",
      " 1850.95616438]\n",
      "-31.395203 16.369354\n",
      "1.3947426 362.0332\n"
     ]
    }
   ],
   "source": [
    "print(time)\n",
    "print(np.min(T), np.max(T))\n",
    "print(np.min(P), np.max(P))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the wrapper"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to run the wrapper to generate simulated TRW."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 167 ms, sys: 24.5 ms, total: 192 ms\n",
      "Wall time: 190 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "res = PyVSL.VSL_R(\n",
    "        850, 1850, 45,                  # the starting year and ending year of the input T & P, along with the latitude\n",
    "        T, P,\n",
    "        T1=1, T2=15, M1=0.01, M2=0.05,  # parameters of the thresholds\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['trw', 'gT', 'gM', 'gE', 'M', 'potEv', 'sample.mean.width', 'sample.std.width', 'trw_org'])\n"
     ]
    }
   ],
   "source": [
    "print(res.keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.03242287  1.62074875  2.00288833 ... -0.87266271  0.07929194\n",
      " -1.20604461]\n"
     ]
    }
   ],
   "source": [
    "print(res['trw'])  # this is the simulated TRW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}