initial import of pyGIMLI notebook to do a VES inversion

VES.ipynb

@@ -0,0 +1,301 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Checkout www.pygimli.org for more examples"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "# VES inversion for a blocky model\n",
+    "\n",
+    "This tutorial shows how an built-in forward operator is used for inversion.\n",
+    "A DC 1D (VES) modelling is used to generate data, noisify and invert them.\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We import numpy, matplotlib and the 1D plotting function\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import pygimli as pg\n",
+    "from pygimli.physics import VESManager"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "some definitions before (model, data and error)\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "#ab2 = np.logspace(-0.5, 2.5, 40)  # AB/2 distance (current electrodes)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "define a synthetic model and do a forward simulatin including noise\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "#synres = [100., 500., 30.]  # synthetic resistivity\n",
+    "#synthk = [0.5, 3.5]  # synthetic thickness (nlay-th layer is infinite)"
+   ]
+  },
+  {
+   "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,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "#synthModel = synthk + synres  # concatenate thickness and resistivity\n",
+    "#ves = VESManager()\n",
+    "#rhoa, err = ves.simulate(synthModel, ab2=ab2, mn2=0.5,\n",
+    "#                         noiseLevel=0.03, seed=1337)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# use data from Bausenberg\n",
+    "ab2 = [1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 4.0, 5.0, 6.0, 7.5, 10.0, 12.5, 15.0, 20.0, 25.0, 30.0, 40.0, 50.0, 60.0, 75.0, 100.]\n",
+    "rhoa = [64.0, 70., 77., 89., 103., 115., 152., 179., 206., 244., 305., 360., 410., 489., 540., 588., 672., 605., 535., 433., 283.]\n",
+    "err = np.ones(len(rhoa))*0.02\n",
+    "err[-4:] = 0.05\n",
+    "mn2 = np.ones(len(rhoa))*0.5"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(ab2)\n",
+    "print(rhoa)\n",
+    "print(err)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "nlay = 3\n",
+    "ves = VESManager()\n",
+    "ves.invert(rhoa, err, ab2=ab2, mn2=mn2,\n",
+    "           nLayers=nlay, lam=1000, lambdaFactor=0.8, cTpye=2, verbose=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "show estimated & synthetic models and data with model response in 2 subplots\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(ncols=2, figsize=(8, 6))  # two-column figure\n",
+    "#ves.showModel(synthModel, ax=ax[0], label=\"synth\", plot=\"semilogy\", zmax=20)\n",
+    "ves.showModel(ves.model, ax=ax[0], label=\"model\", zmax=50, color=\"C1\")\n",
+    "ves.showData(rhoa, ax=ax[1], label=\"data\", color=\"C0\", marker=\"o\")\n",
+    "out = ves.showData(ves.inv.response, ax=ax[1], label=\"model\", color=\"C1\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We are interested in the model uncertaincies and through model covariance\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from pygimli.frameworks.resolution import modelCovariance\n",
+    "var, MCM = modelCovariance(ves.inv)\n",
+    "pg.info(var)\n",
+    "fig, ax = plt.subplots()\n",
+    "im = ax.imshow(MCM, vmin=-1, vmax=1, cmap=\"bwr\")\n",
+    "plt.colorbar(im)\n",
+    "labels = [rf'$d_{i+1}$' for i in range(nlay-1)] + \\\n",
+    "    [rf'$\\rho_{i+1}$' for i in range(nlay)]\n",
+    "plt.xticks(np.arange(nlay*2-1), labels)\n",
+    "_ = plt.yticks(np.arange(nlay*2-1), labels)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The model covariance matrix delivers variances and a scaled (dimensionless)\n",
+    "correlation matrix. The latter show the interdependency of the parameters\n",
+    "among each other. The first and last resistivity is best resolved, also the\n",
+    "first layer thickness. The remaining resistivities and thicknesses are highly\n",
+    "correlated. The variances can be used as error bars in the model plot.\n",
+    "\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "thk = ves.model[:nlay-1]\n",
+    "res = ves.model[nlay-1:]\n",
+    "z = np.cumsum(thk)\n",
+    "mid = np.hstack([z - thk/2, z[-1]*1.1])\n",
+    "resmean = np.sqrt(res[:-1]*res[1:])\n",
+    "fig, ax = plt.subplots()\n",
+    "#ves.showModel(synthModel, ax=ax, label=\"synth\", plot=\"semilogy\", zmax=20)\n",
+    "ves.showModel(ves.model, ax=ax, label=\"model\", zmax=50)\n",
+    "ax.errorbar(res, mid, marker=\"*\", ls=\"None\", xerr=res*var[nlay-1:])\n",
+    "ax.errorbar(resmean, z, marker=\"*\", ls=\"None\", yerr=thk*var[:nlay-1])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "print(res)\n",
+    "print(res*var[nlay-1:])\n",
+    "print(thk)\n",
+    "print(thk*var[:nlay-1])"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "pyGIMLI",
+   "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.11.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

environment.yml

@@ -0,0 +1,8 @@
+name: pyGIMLI
+channels:
+  - conda-forge
+  - defaults
+dependencies:
+  - pygimli>=1.5.0
+  - jupyterlab
+  - suitesparse=5.10.1