examples/errvarexample_freyberg.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib.patches import Rectangle as rect\n",
    "import matplotlib.pyplot as plt\n",
    "import warnings\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import matplotlib as mpl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model background\n",
    "Here is an example based on the model of Freyberg, 1988.  The synthetic model is a 2-dimensional MODFLOW model with 1 layer,  40 rows, and 20 columns.  The model has 2 stress periods: an initial steady-state stress period used for calibration, and a 5-year transient stress period.  The calibration period uses the recharge and well flux of Freyberg, 1988; the last stress period use 25% less recharge and 25% more pumping.\n",
    "\n",
    "The inverse problem has 761 parameters: hydraulic conductivity of each active model cell, calibration and forecast period recharge multipliers, storage and specific yield, calibration and forecast well flux for each of the six wells, and river bed conductance for each 40 cells with river-type boundary conditions.  The inverse problem has 12 head obseravtions, measured at the end of the steady-state calibration period.  The forecasts of interest include the sw-gw exchange flux during both stress periods (observations named ```sw_gw_0``` and ``sw_gw_1``), and the water level in well cell 6 located in at row 28 column 5 at the end of the stress periods (observations named ```or28c05_0``` and ```or28c05_1```).  The forecasts are included in the Jacobian matrix as zero-weight observations. The model files, pest control file and previously-calculated jacobian matrix are in the `freyberg/` folder\n",
    "\n",
    "\n",
    "Freyberg, David L. \"AN EXERCISE IN GROUND‐WATER MODEL CALIBRATION AND PREDICTION.\" Groundwater 26.3 (1988): 350-360."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import flopy\n",
    "\n",
    "# load the model\n",
    "model_ws = os.path.join(\"Freyberg\",\"extra_crispy\")\n",
    "ml = flopy.modflow.Modflow.load(\"freyberg.nam\",model_ws=model_ws)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Because this model is old -- it predates flopy's modelgrid implementation. \n",
    "# And because modelgrid has been implemented without backward compatability \n",
    "# the modelgrid object is not constructed properly. \n",
    "# - We will use some sneaky pyemu to get things to how they should be \n",
    "import pyemu\n",
    "sr = pyemu.helpers.SpatialReference.from_namfile(\n",
    "    os.path.join(model_ws, ml.namefile), \n",
    "    delc=ml.dis.delc, \n",
    "    delr=ml.dis.delr\n",
    ")\n",
    "ml.modelgrid.set_coord_info(\n",
    "    xoff=sr.xll,\n",
    "    yoff=sr.yll,\n",
    "    angrot=sr.rotation,\n",
    "    proj4=sr.proj4_str,\n",
    "    merge_coord_info=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot some model attributes\n",
    "fig = plt.figure(figsize=(10,10))\n",
    "ax = plt.subplot(111,aspect=\"equal\")\n",
    "ml.upw.hk.plot(axes=ax,colorbar=\"K m/d\",alpha=0.3)\n",
    "ml.wel.plot(axes=ax)  # flopy possibly now only plots BCs in black\n",
    "ml.riv.plot(axes=ax)\n",
    "\n",
    "# plot obs locations\n",
    "obs = pd.read_csv(os.path.join(\"Freyberg\",\"misc\",\"obs_rowcol.dat\"),\n",
    "                  delim_whitespace=True)\n",
    "obs_x = [ml.modelgrid.xcellcenters[r-1,c-1] for r,c in obs.loc[:,[\"row\",\"col\"]].values]\n",
    "obs_y = [ml.modelgrid.ycellcenters[r-1,c-1] for r,c in obs.loc[:,[\"row\",\"col\"]].values]\n",
    "ax.scatter(obs_x,obs_y,marker='.',label=\"obs\" )\n",
    "\n",
    "#plot names on the pumping well locations\n",
    "wel_data = ml.wel.stress_period_data[0]\n",
    "wel_x = ml.modelgrid.xcellcenters[wel_data[\"i\"],wel_data[\"j\"]]\n",
    "wel_y = ml.modelgrid.ycellcenters[wel_data[\"i\"],wel_data[\"j\"]]\n",
    "for i,(x,y) in enumerate(zip(wel_x,wel_y)):\n",
    "    ax.text(x,y,\"{0}  \".format(i+1),ha=\"right\",va=\"center\", font=dict(size=15), color='r')\n",
    "\n",
    "ax.set_ylabel(\"y\")\n",
    "ax.set_xlabel(\"x\")\n",
    "\n",
    "ax.add_patch(rect((0,0),0,0,label=\"well\",ec=\"none\",fc=\"r\"))\n",
    "ax.add_patch(rect((0,0),0,0,label=\"river\",ec=\"none\",fc=\"g\"))\n",
    "\n",
    "ax.legend(bbox_to_anchor=(1.5,1.0),frameon=False)\n",
    "plt.savefig(\"domain.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using `pyemu`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyemu"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First create a linear_analysis object.  We will use `ErrVar`  derived type, which replicates the behavior of the `PREDVAR` suite of PEST as well as `ident_par` utility.  We pass it the name of the jacobian matrix file.  Since we don't pass an explicit argument for `parcov` or `obscov`, `pyemu` attempts to build them from the parameter bounds and observation weights in a pest control file (.pst) with the same base case name as the jacobian.  Since we are interested in forecast uncertainty as well as parameter uncertainty, we also pass the names of the forecast sensitivity vectors we are interested in, which are stored in the jacobian as well.  Note that the `forecasts` argument can be a mixed list of observation names, other jacobian files or PEST-compatible ASCII matrix files.  Remember you can pass a filename to the `verbose` argument to write log file.\n",
    "\n",
    "Since most groundwater model history-matching analyses focus on adjusting hetergeneous hydraulic properties and not boundary condition elements, let's identify the well flux and recharge parameters as `omitted` in the error variance analysis.  We can conceptually think of this action as excluding these parameters from the history-matching process. Later we will explicitly calculate the penalty for not adjusting these parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the list of forecast names from the pest++ argument\n",
    "# in the pest control file\n",
    "jco = os.path.join(\"Freyberg\",\"freyberg.jcb\")\n",
    "pst = pyemu.Pst(jco.replace(\"jcb\",\"pst\"))\n",
    "omitted = [pname for pname in pst.par_names if \\\n",
    "           pname.startswith(\"wf\") or pname.startswith(\"rch\")]\n",
    "forecasts = pst.pestpp_options[\"forecasts\"].split(',')\n",
    "la = pyemu.ErrVar(jco=jco,verbose=\"errvar_freyberg.log\",\n",
    "                  omitted_parameters=omitted)\n",
    "print(\"observations, parameters found in jacobian:\",la.jco.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameter identifiability\n",
    "The `errvar` dervied type exposes a method to get a `pandas` dataframe of parameter identifiability information.  Recall that parameter identifiability is expressed as $d_i = \\Sigma(\\mathbf{V}_{1i})^2$, where $d_i$ is the parameter identifiability, which ranges from 0 (not identified by the data) to 1 (full identified by the data), and $\\mathbf{V}_1$ are the right singular vectors corresonding to non-(numerically) zero singular values.  First let's look at the singular spectrum of $\\mathbf{Q}^{\\frac{1}{2}}\\mathbf{J}$, where $\\mathbf{Q}$ is the cofactor matrix and $\\mathbf{J}$ is the jacobian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = la.qhalfx.s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pylab as plt\n",
    "figure = plt.figure(figsize=(10, 5))\n",
    "ax = plt.subplot(111)\n",
    "ax.plot(s.x)\n",
    "ax.set_title(\"singular spectrum\")\n",
    "ax.set_ylabel(\"power\")\n",
    "ax.set_xlabel(\"singular value\")\n",
    "ax.set_xlim(0,20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the singluar spectrum decays rapidly (not uncommon) and that we can really only support about 12 right singular vectors even though we have 700+ parameters in the inverse problem.  \n",
    "\n",
    "Let's get the identifiability dataframe at 12 singular vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the method is passed the number of singular vectors to include in V_1\n",
    "ident_df = la.get_identifiability_dataframe(12) \n",
    "ident_df.sort_values(by=\"ident\",ascending=False).iloc[0:10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the indentifiability:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = ident_df.sort_values(by=\"ident\",ascending=False).iloc[0:20].\\\n",
    "     loc[:,\"ident\"].plot(kind=\"bar\",figsize=(10,10))\n",
    "ax.set_ylabel(\"identifiability\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Forecast error variance \n",
    "\n",
    "Now let's explore the error variance of the forecasts we are interested in.  We will use an extended version of the forecast error variance equation:   \n",
    "\n",
    "$\\sigma_{s - \\hat{s}}^2 = \\underbrace{\\textbf{y}_i^T({\\bf{I}} - {\\textbf{R}})\\boldsymbol{\\Sigma}_{{\\boldsymbol{\\theta}}_i}({\\textbf{I}} - {\\textbf{R}})^T\\textbf{y}_i}_{1} + \\underbrace{{\\textbf{y}}_i^T{\\bf{G}}\\boldsymbol{\\Sigma}_{\\mathbf{\\epsilon}}{\\textbf{G}}^T{\\textbf{y}}_i}_{2} + \\underbrace{{\\bf{p}}\\boldsymbol{\\Sigma}_{{\\boldsymbol{\\theta}}_o}{\\bf{p}}^T}_{3}$\n",
    "\n",
    "Where term 1 is the null-space contribution, term 2 is the solution space contribution and term 3 is the model error term (the penalty for not adjusting uncertain parameters).  Remember the well flux and recharge parameters that we marked as omitted?  The consequences of that action can now be explicitly evaluated.  See Moore and Doherty (2005) and White and other (2014) for more explanation of these terms.  Note that if you don't have any `omitted_parameters`, the only terms 1 and 2 contribute to error variance\n",
    "\n",
    "First we need to create a list (or numpy ndarray) of the singular values we want to test.  Since we have 12 data, we only need to test up to $13$ singular values because that is where the action is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sing_vals = np.arange(13)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `ErrVar` derived type exposes a method to get a multi-index pandas dataframe with each of the terms of the error variance equation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errvar_df = la.get_errvar_dataframe(sing_vals)\n",
    "errvar_df.iloc[0:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "errvar_df[[\"first\"]].to_latex(\"sw_gw_0.tex\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plot the error variance components for each forecast:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = {\"first\": 'g', \"second\": 'b', \"third\": 'c'}\n",
    "max_idx = 19\n",
    "idx = sing_vals[:max_idx]\n",
    "for ipred, pred in enumerate(forecasts):\n",
    "    pred = pred.lower()\n",
    "    fig = plt.figure(figsize=(10, 10))\n",
    "    ax = plt.subplot(111)\n",
    "    ax.set_title(pred)\n",
    "    first = errvar_df[(\"first\", pred)][:max_idx]\n",
    "    second = errvar_df[(\"second\", pred)][:max_idx]\n",
    "    third = errvar_df[(\"third\", pred)][:max_idx]\n",
    "    ax.bar(idx, first, width=1.0, edgecolor=\"none\", \n",
    "           facecolor=colors[\"first\"], label=\"first\",bottom=0.0)\n",
    "    ax.bar(idx, second, width=1.0, edgecolor=\"none\", \n",
    "           facecolor=colors[\"second\"], label=\"second\", bottom=first)\n",
    "    ax.bar(idx, third, width=1.0, edgecolor=\"none\", \n",
    "           facecolor=colors[\"third\"], label=\"third\", \n",
    "           bottom=second+first)\n",
    "    ax.set_xlim(-1,max_idx+1)\n",
    "    ax.set_xticks(idx+0.5)\n",
    "    ax.set_xticklabels(idx)\n",
    "    #if ipred == 2:\n",
    "    ax.set_xlabel(\"singular value\")\n",
    "    ax.set_ylabel(\"error variance\")\n",
    "    ax.legend(loc=\"upper right\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we see the trade off between getting a good fit to push down the null-space (1st) term and the penalty for overfitting (the rise of the solution space (2nd) term)).  The sum of the first two terms in the \"appearent\" error variance (e.g. the uncertainty that standard analyses would yield) without considering the contribution from the omitted parameters.  You can verify this by checking prior uncertainty from the Schur's complement notebook against the zero singular value result using only terms 1 and 2. Note that the top of the green bar is the limit of traditional uncertainty/error variance analysis: accounting for parameter and observation\n",
    "\n",
    "We also see the added penalty for not adjusting the well flux and recharge parameters.  For the water level at the end of the calibration period forecast (``or28c05_0``), the fact the we have left parameters out doesn't matter - the parameter compensation associated with fixing uncertain model inputs can be \"calibrated out\" beyond 2 singular values.  For the water level forecast during forecast period (``or28c05_1``), the penalty for fixed parameters persists -it s nearly constant over the range of singular values.  \n",
    "\n",
    "For ``sw_gw_0``, the situation is much worse: not only are we greatly underestimating uncertainty by omitting parameters, worse, calibration increases the uncertainty for this forecast because the adjustable parametres are compensating for the omitted, uncertaint parameters in ways that are damanaging to the forecast. \n",
    "\n",
    "For the forecast period sw-gw exchange (``sw_gw_1``), calibration doesn't help or hurt - this forecast depend entirely on null space parameter components.  But treating the recharge and well pumpage as \"fixed\" (omitted) results in greatly underestimated uncertainty.     \n",
    "\n",
    "Let's check the ```errvar``` results against the results from ```schur```. This is simple with ```pyemu```, we simply  cast the ```errvar``` type to a ```schur``` type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "schur = la.get(astype=pyemu.Schur)\n",
    "schur_prior = schur.prior_forecast\n",
    "schur_post = schur.posterior_forecast\n",
    "print(\"{0:10s} {1:>12s} {2:>12s} {3:>12s} {4:>12s}\"\n",
    "      .format(\"forecast\",\"errvar prior\",\"errvar min\",\n",
    "              \"schur prior\", \"schur post\"))\n",
    "for ipred, pred in enumerate(forecasts):\n",
    "    first = errvar_df[(\"first\", pred)][:max_idx]\n",
    "    second = errvar_df[(\"second\", pred)][:max_idx]  \n",
    "    min_ev = np.min(first + second)\n",
    "    prior_ev = first[0] + second[0]\n",
    "    prior_sh = schur_prior[pred]\n",
    "    post_sh = schur_post[pred]\n",
    "    print(\"{0:12s} {1:12.6f} {2:12.6f} {3:12.6} {4:12.6f}\"\n",
    "          .format(pred,prior_ev,min_ev,prior_sh,post_sh))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the prior from ```schur``` class matches the two-term ```errvar``` result at zero singular values.  We also see, as expected, the posterior from ```schur``` is slightly lower than the two-term ```errvar``` result.  This shows us that the \"appearent\" uncertainty in these predictions, as found through application of Bayes equation, is being under estimated because if the ill effects of the omitted parameters."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}