{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Checkout www.pygimli.org for more examples\n%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nVES inversion for a smooth model\n================================\n\nThis tutorial shows how an built-in forward operator is used for an Occam type\n(smoothness-constrained) inversion with fixed regularization (most natural).\nA direct current (DC) one-dimensional (1D) VES (vertical electric sounding)\nmodelling operator is used to generate data, add noise and inversion.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We import numpy numerics, mpl plotting, pygimli and the 1D plotting function\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nnp.random.seed(1337)\n\nimport matplotlib.pyplot as plt\n\nimport pygimli as pg\nfrom pygimli.viewer.mpl import drawModel1D"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "##############################################################################\n initialize the forward modelling operator and compute synthetic noisy data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "synres = [100., 500., 20., 800.]  # synthetic resistivity\nsynthk = [4, 6, 10]  # synthetic thickness (lay layer is infinite)\nab2 = np.logspace(-1, 2, 25)  # 0.1 to 100 in 25 steps (8 points per decade)\nfBlock = pg.core.DC1dModelling(len(synres), ab2, ab2/3)\nrhoa = fBlock(synthk+synres)\n# The data are noisified using a\nerrPerc = 3.  # relative error of 3 percent\nrhoa = rhoa * (np.random.randn(len(rhoa)) * errPerc / 100. + 1.)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The forward operator can be called by f.response(model) or simply f(model)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "thk = np.logspace(-0.5, 0.5, 30)\nf = pg.core.DC1dRhoModelling(thk, ab2, ab2/3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create some transformations used for inversion\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "transRho = pg.trans.TransLogLU(1, 1000)  # lower and upper bound\ntransRhoa = pg.trans.TransLog()  # log transformation also for data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set up inversion\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "inv = pg.core.Inversion(rhoa, f, transRhoa, transRho, False)  # data vector, f, ...\n# The transformations can also be omitted and set individually by\n# inv.setTransData(transRhoa)\n# inv.setTransModel(transRho)\ninv.setRelativeError(errPerc / 100.0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create a homogeneous starting model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "model = pg.Vector(len(thk)+1, np.median(rhoa))  # uniform values\ninv.setModel(model)  #"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Set pretty large regularization strength and run inversion\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "print(\"inversion with lam=200\")\ninv.setLambda(100)\nres100 = inv.run()  # result is a pg.Vector, but compatible to numpy array\nprint('rrms={:.2f}%, chi^2={:.3f}'.format(inv.relrms(), inv.chi2()))\n# Decrease the regularization (smoothness) and start (from old result)\nprint(\"inversion with lam=20\")\ninv.setLambda(10)\nres10 = inv.run()  # result is a pg.Vector, but compatible to numpy array\nprint('rrms={:.2f}%, chi^2={:.3f}'.format(inv.relrms(), inv.chi2()))\n# We now optimize lambda such that data are fitted within noise (chi^2=1)\nprint(\"chi^2=1 optimized inversion\")\nresChi = inv.runChi1()  # ends up in a lambda of about 3\nprint(\"optimized lambda value:\", inv.getLambda())\nprint('rrms={:.2f}%, chi^2={:.3f}'.format(inv.relrms(), inv.chi2()))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Show model (inverted and synthetic) as well as model response and data\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(ncols=2, figsize=(8, 6))  # two-column figure\ndrawModel1D(ax[0], synthk, synres, color='b', label='synthetic',\n            plot='semilogx')\ndrawModel1D(ax[0], thk, res100, color='g', label=r'$\\lambda$=100')\ndrawModel1D(ax[0], thk, res10, color='c', label=r'$\\lambda$=10')\ndrawModel1D(ax[0], thk, resChi, color='r',\n            label=r'$\\chi$={:.2f}'.format(inv.getLambda()))\nax[0].grid(True, which='both')\nax[0].legend(loc='best')\nax[1].loglog(rhoa, ab2, 'rx-', label='measured')\nax[1].loglog(inv.response(), ab2, 'b-', label='fitted')\nax[1].set_ylim((max(ab2), min(ab2)))\nax[1].grid(True, which='both')\nax[1].set_xlabel(r'$\\rho_a$ [$\\Omega$m]')\nax[1].set_ylabel('AB/2 [m]')\nax[1].yaxis.set_label_position('right')\nax[1].yaxis.set_ticks_position('right')\nax[1].legend(loc='best')\nplt.show()"
      ]
    }
  ],
  "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": 0
}