{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effect of noise and height on the inversion\n",
    "This study has been demonstrated in the EMagPy paper [McLachlan et al. (2021)](https://doi.org/10.1016/j.cageo.2020.104561).\n",
    "\n",
    "All inversions are performed with the ROPE solver on a two-layer model with a varying depth. (a) Inversion with 0% noise with device on the ground. (b) Inversion with 5% noise on the ground. (c) Inversion with 0% noise at 1 m above the ground (d) Inversion with 5% noise at 1 m above the ground. The red line represents the true interface between the two layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import sys\n",
    "sys.path.append('../src') # add path where emagpy is\n",
    "from emagpy import Problem\n",
    "\n",
    "datadir = '../src/examples/'\n",
    "\n",
    "letters = ['a','b','c','d','e','f','g','h','i','j']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters for the synthetic model\n",
    "nlayer = 2 # number of layer\n",
    "npos = 20 # number of sampling positions\n",
    "conds = np.ones((npos, nlayer))*[20, 100]\n",
    "x = np.linspace(0.1, 2, npos)[:,None]\n",
    "depths = 0.65 + 0.2 * np.sin(x*np.pi*2) # wave\n",
    "coils0 = ['VCP1.48f10000h0', 'VCP2.82f10000h0', 'VCP4.49f10000h0',\n",
    "          'HCP1.48f10000h0', 'HCP2.82f10000h0', 'HCP4.49f10000h0']\n",
    "coils1 = ['VCP1.48f10000h1', 'VCP2.82f10000h1', 'VCP4.49f10000h1',\n",
    "          'HCP1.48f10000h1', 'HCP2.82f10000h1', 'HCP4.49f10000h1']\n",
    "coils = [coils0, coils0, coils1, coils1]\n",
    "noises = [0, 0.05, 0, 0.05]\n",
    "ks = []\n",
    "# generate ECa using forward model\n",
    "for i in range(4):\n",
    "    k = Problem()\n",
    "    k.setModels([depths], [conds])\n",
    "    _ = k.forward(forwardModel='FSlin', coils=coils[i], noise=noises[i])\n",
    "    ks.append(k)\n",
    "\n",
    "# invert\n",
    "for k in ks:\n",
    "    k.setInit(depths0=np.array([0.5]), fixedDepths=[False])\n",
    "    k.invert(forwardModel='FSlin', method='ROPE', regularization='l1',\n",
    "             bnds=[(0.05, 2.5), (5, 150), (5, 150)], rep=1000, njobs=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# graph of inverted values\n",
    "titles = ['(a) at 0 m (no noise)', '(b) at 0 m (with 5% noise)',\n",
    "          '(c) at 1 m (no noise)', '(d) at 1 m (with 5% noise)']\n",
    "fig, axs = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(14,3))\n",
    "for i in range(4):\n",
    "    ax = axs[i]\n",
    "    ks[i].showResults(ax=ax, vmin=0, vmax=120, maxDepth=2, errorbar=True)\n",
    "    rmseDepths = np.sqrt(np.sum((ks[i].depths[0][:,0] - depths[:,0])**2)/len(depths[:,0]))\n",
    "    ax.set_title('{:s} \\nRMSE$_{{depth}}$={:.2f}'.format(titles[i], rmseDepths))\n",
    "    ax.step(np.arange(depths.shape[0])+0.5, -depths, 'r', where='post') # true depth\n",
    "    if i < 3:\n",
    "        fig.axes[-1].remove()\n",
    "    if i > 0:\n",
    "        ax.set_ylabel('')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# graph of apparent values\n",
    "titles = ['at 0 m (no noise)', 'at 0 m (with 5% noise)',\n",
    "          'at 1 m (no noise)', 'at 1 m (with 5% noise)']\n",
    "fig, axs = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(14,3))\n",
    "for i in range(4):\n",
    "    ax = axs[i]\n",
    "    ks[i].show(ax=ax, vmin=10, vmax=100)\n",
    "    ax.set_title('({:s}) {:s}'.format(letters[i], titles[i]))\n",
    "    if i > 0:\n",
    "        ax.get_legend().remove()\n",
    "        ax.set_ylabel('')"
   ]
  }
 ],
 "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}