{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Smoothing parameters and time-lapse inversion\n",
    "In this notebook, we will see what it the impact of the different smoothing parameters and how we can use them to perform time-lapse inversion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-21T23:26:32.697117Z",
     "iopub.status.busy": "2021-01-21T23:26:32.695952Z",
     "iopub.status.idle": "2021-01-21T23:26:33.607431Z",
     "shell.execute_reply": "2021-01-21T23:26:33.608418Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sys\n",
    "sys.path.append('../src')\n",
    "from emagpy import Problem\n",
    "datadir = '../src/examples/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": "false"
   },
   "source": [
    "## Smoothing parameters\n",
    "There are three different smoothing parameters:\n",
    "- `alpha`: linked to the vertical smoothing\n",
    "- `beta`: linked to lateral smoothing\n",
    "- `gamma`: linked to smoothing with the first inverted survey (usefull for time-lapse inversion)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-21T23:26:33.612614Z",
     "iopub.status.busy": "2021-01-21T23:26:33.611388Z",
     "iopub.status.idle": "2021-01-21T23:26:33.864547Z",
     "shell.execute_reply": "2021-01-21T23:26:33.865474Z"
    }
   },
   "outputs": [],
   "source": [
    "k = Problem()\n",
    "k.createSurvey(datadir + 'cover-crop/coverCrop.csv')\n",
    "k.surveys[0].df = k.surveys[0].df[:30] # truncate dataset\n",
    "k.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": "false"
   },
   "source": [
    "### Vertical smoothing (`alpha`)\n",
    "The `alpha` parameter regulates the vertical smoothing of the profile. It is needed to regularize the inversion.\n",
    "It plays an important role in smooth inversion with many layers such as used with `Problem.invertGN()`.\n",
    "To determine the optimal value of `alpha` a common method is to plot the so-called 'L-curve' that show the size of the model and data misfit as a function of alpha. It can be plotted using the `Problem.lcurve()` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-21T23:26:33.869575Z",
     "iopub.status.busy": "2021-01-21T23:26:33.868965Z",
     "iopub.status.idle": "2021-01-21T23:26:34.717451Z",
     "shell.execute_reply": "2021-01-21T23:26:34.718372Z"
    }
   },
   "outputs": [],
   "source": [
    "k.lcurve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": "false"
   },
   "source": [
    "From the L-curve, we look for the alpha value that offers the best trade-off between the model and data misfit. This value is usually taken at the maximum curvature point of the curve. In this case it would corresponds to 0.07. This value can then be used in the call to `Problem.invert(alpha=alpha)`.\n",
    "\n",
    "The graphs below shows the impact of different values of `alpha`. Larger alpha values will give a smoother model, while smaller alpha values can show very high/low EC values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-21T23:26:34.722649Z",
     "iopub.status.busy": "2021-01-21T23:26:34.721394Z",
     "iopub.status.idle": "2021-01-21T23:26:37.223044Z",
     "shell.execute_reply": "2021-01-21T23:26:37.224028Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))\n",
    "for i, alpha in enumerate([0, 0.07, 0.5]):\n",
    "    ax = axs[i]\n",
    "    k.invert(alpha=alpha)\n",
    "    k.showResults(ax=ax)\n",
    "    ax.set_title('alpha = {:.2f}'.format(alpha))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": "false"
   },
   "source": [
    "### Lateral smoothing (`beta`)\n",
    "We can see that the data present some similarity between continguous measurements. To force lateral smoothing we can use the `beta` parameters. It will penalize the objective function with the difference between two continuous profiles. Note that this can be dangerous in some conditions, for instance if the first profile is badly inverted.\n",
    "\n",
    "The figure below show the effect of lateral smoothing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-21T23:26:37.228418Z",
     "iopub.status.busy": "2021-01-21T23:26:37.227200Z",
     "iopub.status.idle": "2021-01-21T23:26:39.423142Z",
     "shell.execute_reply": "2021-01-21T23:26:39.424080Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))\n",
    "for i, beta in enumerate([0, 0.1, 0.5]):\n",
    "    ax = axs[i]\n",
    "    k.invert(beta=beta)\n",
    "    k.showResults(ax=ax)\n",
    "    ax.set_title('beta = {:.2f}'.format(beta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": "false"
   },
   "source": [
    "### Reference survey smoothing (`gamma`) or time-lapse inversion\n",
    "`gamma` controls the smoothing between the first survey and subsequent survey. This way a quasi time-lapse inversion can be performed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-01-21T23:26:39.428447Z",
     "iopub.status.busy": "2021-01-21T23:26:39.427196Z",
     "iopub.status.idle": "2021-01-21T23:26:39.570827Z",
     "shell.execute_reply": "2021-01-21T23:26:39.568867Z"
    }
   },
   "outputs": [],
   "source": [
    "k = Problem()\n",
    "k.createTimeLapseSurvey([datadir + 'timelapse-wheat/170316.csv',\n",
    "                        datadir + 'timelapse-wheat/170427.csv',\n",
    "                        datadir + 'timelapse-wheat/170516.csv'])\n",
    "for s in k.surveys:\n",
    "    s.df = s.df[:20]\n",
    "\n",
    "# show apparent values collected\n",
    "fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))\n",
    "for i in range(len(k.surveys)):\n",
    "    ax = axs[i]\n",
    "    k.show(index=i, ax=ax)\n",
    "ax.set_ylim([12, 42])\n",
    "k.setInit(depths0=[0.7], fixedDepths=[False])\n",
    "bounds = [(0.1, 1),(2, 60),(2, 60)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inversion with gamma = 0 (no smoothing between surveys)\n",
    "k.invert(gamma=0, bnds=bounds)\n",
    "fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))\n",
    "for i in range(3):\n",
    "    ax = axs[i]\n",
    "    k.showResults(index=i, ax=ax, maxDepth=1.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# inversion with constrain to the first survey\n",
    "k.invert(gamma=0.1, bnds=bounds)\n",
    "fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))\n",
    "for i in range(3):\n",
    "    ax = axs[i]\n",
    "    k.showResults(index=i, ax=ax, maxDepth=1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": "false"
   },
   "source": [
    "### Observations\n",
    "In the second inversion with between survey constrain we can clearly see that the depths obtained during the first inversion constrain the depth of the other survey. We can also notice that the last two profiles are badly inverted in the first survey and so this propagate for the other subsequent surveys."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": "false"
   },
   "source": [
    "## Inverting change in ECa\n",
    "An other way of producing time-lapse inversion is to actually invert for the change in ECa. This can be done using the CS function and the `Problem.invertGN()` function. See appendix of Whalley et al. (2017) for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = Problem()\n",
    "k.createTimeLapseSurvey([datadir + 'timelapse-wheat/170316.csv',\n",
    "                        datadir + 'timelapse-wheat/170427.csv',\n",
    "                        datadir + 'timelapse-wheat/170516.csv'])\n",
    "for s in k.surveys:\n",
    "    s.df = s.df[:20]\n",
    "    \n",
    "k.computeApparentChange()\n",
    "# show change in apparent values collected\n",
    "fig, axs = plt.subplots(1, 3, sharex=True, figsize=(12,3))\n",
    "for i in range(len(k.surveys)):\n",
    "    ax = axs[i]\n",
    "    k.show(index=i, ax=ax)\n",
    "    if i > 0:\n",
    "        ax.set_ylabel('Change in EC [mS/m]')\n",
    "fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k.setInit(depths0=np.linspace(0.2, 2, 10)) # use a smooth model\n",
    "k.invertGN(alpha=0.1)\n",
    "fig, axs = plt.subplots(1, 3, sharex=True, sharey=True, figsize=(12,3))\n",
    "for i in range(3):\n",
    "    ax = axs[i]\n",
    "    if i == 0: # reference survey\n",
    "        k.showResults(index=i, ax=ax)\n",
    "    else:\n",
    "        k.showResults(index=i, ax=ax, cmap='bwr_r', vmin=-8, vmax=2)\n",
    "        fig.axes[-1].set_ylabel('EC change [mS/m]')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}