MADS_calibration_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0629b047",
   "metadata": {},
   "source": [
    "# Model-data calibration for the DVM-DOS-TEM model using [MADS](https://github.com/madsjulia/Mads.jl) julia"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4330206",
   "metadata": {},
   "source": [
    "This example illustrates the usage of the MADS (Model Analysis & Decision Support) package. Here, we use MADS package for the data-model calibration, i.e., minimize the difference between observations and the corresponding modeled outputs. In this example, we wrap the DVM-DOS-TEM model with MADS to calibrate model parameters by comparing model outputs with observations from the Murphy Dome (black spruce) site near Fairbanks, Alaska. This example is for STEP1 only and assumes that a user has a background on six calibration steps and has all the required packages installed. Note that this example runs only from the docker container."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1495d0d2",
   "metadata": {},
   "source": [
    "## STEP1-MD1 \n",
    "In STEP1, we calibrate parameter `cmax` for four plant functional types (PFT) by tweaking `cmax` values in the `parameters/cmt_calparbgc.txt` file. Our target values are Gross Primary Productivity (GPP) values. The observed GPP values can be found in `calibration/calibration_targets.py`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d88cada",
   "metadata": {},
   "outputs": [],
   "source": [
    "# STEP1 (MD1)\n",
    "# parameters: cmax\n",
    "# targets: GPP[1:4], Nitrogen unlimited mode\n",
    "\n",
    "import Mads\n",
    "import PyCall\n",
    "@show pwd()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8c5e318",
   "metadata": {},
   "source": [
    "For **STEP1**, we setup the nitrogen unlimited mode. This mode can be enabled by setup the `calib_mode` eqaul to `'GPPAllIgnoringNitrogen'`. Below we embedded the python code inside the julia notebook. `PyCall` package provides the communication between `python` and `julia` languages. A similar simplified example that uses python from julia can be found on the MADS Github [page](https://github.com/madsjulia/Mads.jl/blob/master/notebooks/model_diagnostics_python/model_diagnostics_python.jl). Most of the utility functions required to run the calibration example live in `TEM.py`, which is located in `work/scripts` folder. The community number type (CMT) corresponding to black spruce forest is `cmtnum=1`. The code below uses two functions: `run_TEM` runs the model and grabs the outputs, and `get_targets` grabs the correspoding target values (in this case GPPs) from `calibration/calibration_targets.py` file. Typically vegetation carbon (above ground carbon) takes less time to equilibrate, which is why we use 100 years of pre-run and 200 of equilibrium run, which should suffice for the model to have a rich equilibrium state for vegetation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "652efb93",
   "metadata": {},
   "outputs": [],
   "source": [
    "PyCall.py\"\"\"\n",
    "\n",
    "import sys,os\n",
    "sys.path.append(os.path.join('/work','scripts'))\n",
    "import TEM\n",
    "\n",
    "def run_TEM(x):\n",
    "    \n",
    "    # update param files\n",
    "    for j in range(len(dvmdostem.params)):\n",
    "        dvmdostem.params[j]['val']=x[j]   \n",
    "    \n",
    "    dvmdostem.clean()             # clean results from previous run\n",
    "    dvmdostem.setup(calib=True)   # setup for a new run\n",
    "    dvmdostem.update_params()     # save updated paramters into the cmt_calparbgc.txt\n",
    "    dvmdostem.run()               # run the model\n",
    "\n",
    "    return dvmdostem.get_calibration_outputs()[:4]\n",
    "\n",
    "def get_targets():\n",
    "    return dvmdostem.get_calibration_outputs(calib=True)[:4]\n",
    "\n",
    "# istantiate the model\n",
    "dvmdostem=TEM.TEM_model()\n",
    "# set the path the input data  \n",
    "dvmdostem.site='/data/input-catalog/cru-ts40_ar5_rcp85_mri-cgcm3_MurphyDome_10x10'\n",
    "# set the path to the ouput folder\n",
    "dvmdostem.work_dir='/data/workflows/MD1'\n",
    "# set the mode\n",
    "dvmdostem.calib_mode='GPPAllIgnoringNitrogen'\n",
    "# set the run setup\n",
    "dvmdostem.opt_run_setup='--pr-yrs 100 --eq-yrs 200 --sp-yrs 0 --tr-yrs 0 --sc-yrs 0'\n",
    "# define parameters which will participate in calibration \n",
    "dvmdostem.set_params(cmtnum=1, params=['cmax','cmax','cmax','cmax'], \\\n",
    "                               pftnums=[0,1,2,3])\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0a18b1b",
   "metadata": {},
   "source": [
    "Define the initial guess (IG) vector. Typically comes from the `cmt_calparbgc.txt`. These IG values need to be updated during the calibration process. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fba10d8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#initial_guess=[385.0, 115.0, 201.0, 95.0]\n",
    "initial_guess=[381.2, 113.9, 207.1, 93.3]\n",
    "# the y_init line runs the run_TEM functions in the /data/workflows/MD1 folder\n",
    "# y_init corresponds to the modeled GPP outputs\n",
    "y_init=PyCall.py\"run_TEM\"(initial_guess)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2d7f280",
   "metadata": {},
   "source": [
    "Let's grab the target GPP values and get their length. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19553fcc",
   "metadata": {},
   "outputs": [],
   "source": [
    "targets=PyCall.py\"get_targets\"()\n",
    "obs_time=1:length(targets)\n",
    "targets"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dde86150",
   "metadata": {},
   "source": [
    "Ths function below serves as a broker between `MADS` julia and python's `run_TEM` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a22e583",
   "metadata": {},
   "outputs": [],
   "source": [
    "function TEM_pycall(parameters::AbstractVector)\n",
    "        predictions = PyCall.py\"run_TEM\"(parameters)\n",
    "        return predictions\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db8e9ceb",
   "metadata": {},
   "source": [
    "Below we set the calibration problem (`createproblem`) to minimize the difference between GPP observed and GPP modeled. `paramdist` sets up the range for the calibrated parameters (see `paramkey`). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f80d800d",
   "metadata": {},
   "outputs": [],
   "source": [
    "md = Mads.createproblem(\n",
    "    initial_guess,               #IG vector to start the calibration\n",
    "    targets,                     #obsevations\n",
    "    TEM_pycall;                  #GPP modeled\n",
    "    paramkey=[\"cmax0\",\"cmax1\",\"cmax2\",\"cmax3\"],\n",
    "    #381.2, 113.9, 207.1, 93.3\n",
    "    paramdist=[\"Uniform(370, 390)\",\"Uniform(110, 120)\",\"Uniform(200, 210)\",\"Uniform(90, 100)\"],\n",
    "    obstime=obs_time,            #length of the observation vector\n",
    "    obsweight=[100,100,100,100], #weight for target values\n",
    "    problemname=\"STEP1-MD1\")     "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "532b1140",
   "metadata": {},
   "source": [
    "Dislpay parameter values and observations (targets). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82e9cdc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "Mads.showparameters(md)\n",
    "Mads.showobservations(md)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5285a717",
   "metadata": {},
   "source": [
    "Calibrate functions will run the calibration algorithm. In this case, by default, it is using [Levenberg-Marquardt Algorithm](https://en.wikipedia.org/wiki/Levenberg%E2%80%93Marquardt_algorithm) to minimize the difference between modeled and observed GPPs. `tolOF` sets the tolerance number for the minimization function. `OF` stands for the root mean square difference between observations and modeled outputs.  `tolOFcount` sets the number of times after which it stops if there are no improvements in the `OF` number. This step could take a while to finish.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6635e259",
   "metadata": {},
   "outputs": [],
   "source": [
    "calib_param, calib_information = Mads.calibrate(md, tolOF=0.01, tolOFcount=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a50a6a38",
   "metadata": {},
   "source": [
    "Save the plot of the match between observation and modeled outputs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73a66964",
   "metadata": {},
   "outputs": [],
   "source": [
    "Mads.plotmatches(md, calib_param, \n",
    "    xtitle=\"# of observations\", ytitle=\"GPP\",filename=\"STEP1-MD1-matchplot.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d58fe47",
   "metadata": {},
   "source": [
    "Display the match plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b221292",
   "metadata": {},
   "outputs": [],
   "source": [
    "Mads.display(\"STEP1-MD1-matchplot.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4bd6ca36",
   "metadata": {},
   "source": [
    "Quick sensitivity analysis will help better understand the dependencies between parameters and targets. If you run the calibration, you should have `calib_param` ready. If not, use the second line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5bddd9b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "localsa = Mads.localsa(md; filename=\"model_diagnostics.png\", par=collect(values(calib_param)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfe2568b",
   "metadata": {},
   "outputs": [],
   "source": [
    "localsa = Mads.localsa(md; filename=\"model_diagnostics.png\", par=initial_guess)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1361b566",
   "metadata": {},
   "outputs": [],
   "source": [
    "[Mads.getparamlabels(md) localsa[\"stddev\"]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "153f3bc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "Mads.display(\"model_diagnostics-jacobian.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a199ccc",
   "metadata": {},
   "outputs": [],
   "source": [
    "Mads.display(\"model_diagnostics-eigenmatrix.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f67b7419",
   "metadata": {},
   "outputs": [],
   "source": [
    "Mads.display(\"model_diagnostics-eigenvalues.png\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60354b85",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.7.2",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}