examples/pyemu_ensemble_basics.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"
   ]
  },
  {
   "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\"),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": [
    "The plot shows the Freyberg (1988) model domain.  The colorflood is the hydraulic conductivity ($\\frac{m}{d}$).  Red and green cells coorespond to well-type and river-type boundary conditions. Blue dots indicate the locations of water levels used for calibration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using `pyemu`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyemu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pst = pyemu.Pst(os.path.join(\"Freyberg\",\"freyberg.pst\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Drawing from the prior\n",
    "Now we need a prior-realized ``ParameterEnsemble``, which stores a ``pandas.DataFrame`` under the hood.\n",
    "\n",
    "### ```draw```\n",
    "\n",
    "The ``ParameterEnsemble`` class has several ``draw`` type methods to generate stochastic values from (multivariate) (log) gaussian, uniform and triangular distributions.  Much of what we do is predicated on the gaussian distribution, so let's use that here.\n",
    "\n",
    "The gaussian draw accepts a `cov` arg which can be a `pyemu.Cov` instance.  If this isn't passed, then the draw method constructs a diagonal covariance matrix from the parameter bounds (assuming a certain number of standard deviations represented by the distance between the bounds - the `sigma_range` argument)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,num_reals=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``draw`` also accepts a ``num_reals`` argument to specify the number of draws to make:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pe.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that these ``draw`` methods use initial parameter values in the control file (the `Pst.parameter_data.parval`1 attribute) the $\\boldsymbol{\\mu}$  (mean) prior parameter vector.  To change that, we need to update the parameter values in the control file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pst.parrep(pst.filename.replace(\".pst\",\".par\"))\n",
    "pst.parameter_data.parval1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,num_reals=200)\n",
    "pe.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## plotting\n",
    "\n",
    "Since ```ParameterEnsemble``` stores a ```pandas.DataFrame```, it has all the cool methods and attributes we all love.  Let's compare the results of drawing from a uniform vs a gaussian distribution.  The actual dataframe is stored under the private attribute `ParameterEnsemble._df`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pe = pyemu.ParameterEnsemble.from_uniform_draw(pst=pst,num_reals=1000)\n",
    "ax = plt.subplot(111)\n",
    "pe._df.loc[:,\"rch_1\"].plot(kind=\"hist\",bins=20,ax=ax,alpha=0.5)\n",
    "pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,num_reals=1000)\n",
    "pe._df.loc[:,\"rch_1\"].plot(kind=\"hist\",bins=20,ax=ax,alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gaussian histo go beyond the parameter bound - bad times.  Luckily, `ParameterEnsemble` includes an `enforce` method to apply parameter bounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pe = pyemu.ParameterEnsemble.from_uniform_draw(pst=pst,num_reals=1000)\n",
    "ax = plt.subplot(111)\n",
    "pe._df.loc[:,\"rch_1\"].plot(kind=\"hist\",bins=20,ax=ax,alpha=0.5)\n",
    "pe = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,num_reals=1000)\n",
    "pe.enforce(how=\"reset\")\n",
    "pe._df.loc[:,\"rch_1\"].plot(kind=\"hist\",bins=20,ax=ax,alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## bayes linear monte carlo\n",
    "\n",
    "we can use the bayes linear posterior parameter covariance matrix (aka Schur compliment) to \"precondition\" the realizations using linear algebra so that they hopefully yield a lower phi.  The trick is we just need to pass this\n",
    "posterior covariance matrix to the draw method.  Note this covariance matrix is the second moment of the posterior (under the FOSM assumptions) and the final parameter values is the first moment (which we `parrep`'ed into the `Pst` earlier)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the list of forecast names from the pest++ argument in the pest control file\n",
    "jco = os.path.join(\"Freyberg\",\"freyberg.jcb\")\n",
    "sc = pyemu.Schur(jco=jco)\n",
    "pe_post = pyemu.ParameterEnsemble.from_gaussian_draw(pst=pst,cov=sc.posterior_parameter, num_reals=1000)\n",
    "pe_post.enforce()\n",
    "ax = plt.subplot(111)\n",
    "pe._df.loc[:,\"rch_1\"].plot(kind=\"hist\",bins=20,ax=ax,alpha=0.5)\n",
    "pe_post._df.loc[:,\"rch_1\"].plot(kind=\"hist\",bins=20,ax=ax,alpha=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Now we just need to run this preconditioned ensemble to validate the FOSM assumptions (that the realizations do yield an acceptably low phi and that the relation between parameters and forecasts is linear)\n"
   ]
  }
 ],
 "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
}