2-linear-inversion-app.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoscilabs.inversion.LinearInversionDirect import LinearInversionDirectApp\n",
    "from ipywidgets import interact, FloatSlider, ToggleButtons, IntSlider, FloatText, IntText"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "XoBdBeGBCevcAMVqqZtM",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 1
     },
     "outputId": {
      "block": "nzHSDMZI4foQGe88I1D5",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 1
     }
    }
   },
   "outputs": [],
   "source": [
    "app = LinearInversionDirectApp()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "GAKaomYwcnkDpl2fpbAI",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 1
     }
    }
   },
   "source": [
    "# Linear Inversion App\n",
    "\n",
    "This app is based upon the inversion tutorial: \"INVERSION FOR APPLIED GEOPHYSICS\" by Oldenburg and Li (2005).\n",
    "\n",
    "Douglas W. Oldenburg and Yaoguo Li (2005) 5. Inversion for Applied Geophysics: A Tutorial. Near-Surface Geophysics: pp. 89-150.\n",
    "eISBN: 978-1-56080-171-9 \n",
    "print ISBN: 978-1-56080-130-6 \n",
    "https://doi.org/10.1190/1.9781560801719.ch5 \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "wzGxtmJjFucP0mYtZzVD",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 2
     }
    }
   },
   "source": [
    "## Purpose\n",
    "\n",
    "By using a simple decaying and oscillating kernel function, which emulates the physics of electromagnetic (EM) survey, we understand basic concepts of inverting data. Three items that we are going to explore are:\n",
    "\n",
    "- Step1: Create a model ($\\mathbf{m}$)\n",
    "- Step2: Generate a sensitivity kernel (or matrix), $\\mathbf{G}$\n",
    "- Step3: Simulate data ($\\mathbf{d} = \\mathbf{G}\\mathbf{m}$)\n",
    "- Step4: All three steps together\n",
    "- Step5: Invert the data, and explore inversion results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "TVYiH721oymtWnC2zQ1w",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 2
     }
    }
   },
   "source": [
    "## Forward problem\n",
    "\n",
    "\n",
    "Let $g_j(x)$ denote the kernel function for $j$th datum. With a given model $m(x)$, $j$th datum can be computed by solving following integral equation:\n",
    "\n",
    " $$ d_j = \\int_a^{b} g_j(x) m(x) dx $$\n",
    "\n",
    "where \n",
    "\n",
    "$$ g_j(x) = e^{p_jx} cos (2 \\pi q_jx) $$ \n",
    "\n",
    "By discretizing $g_j(x)$ we obtain"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "s9H8M9Mtzk7wkXBPF8iZ",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 2
     }
    }
   },
   "source": [
    "$$ \\mathbf{g}_j(\\mathbf{x}) = e^{p_j\\mathbf{x}} cos (2 \\pi q_j \\mathbf{x})$$\n",
    "\n",
    "where\n",
    "\n",
    "- $\\mathbf{g}_j$: $j$th row vector for the sensitivty matrix ($1 \\times M$)\n",
    "- $\\mathbf{x}$: model location ($1 \\times M$)\n",
    "- $p_j$: decaying constant (<0)\n",
    "- $q_j$: oscillating constant (>0)\n",
    "\n",
    "By stacking multiple rows of $\\mathbf{g}_j$, we obtain sensitivity matrix, $\\mathbf{G}$: \n",
    "\n",
    "\\begin{align}\n",
    "    \\mathbf{G} = \n",
    "    \\begin{bmatrix}\n",
    "        \\mathbf{g}_1\\\\\n",
    "        \\vdots\\\\\n",
    "        \\mathbf{g}_{N}\n",
    "    \\end{bmatrix}\n",
    "\\end{align}\n",
    "\n",
    "Here, the size of the matrix $\\mathbf{G}$ is $(N \\times M)$. \n",
    "Finally data, $\\mathbf{d}$, can be written as a linear equation:\n",
    "\n",
    "$$ \\mathbf{d} = \\mathbf{G}\\mathbf{m}$$\n",
    "\n",
    "where $\\mathbf{m}$ is an inversion model; this is a column vector ($M \\times 1$). \n",
    "\n",
    "In real measurments, there will be various noises source, and hence observation, $\\mathbf{d}^{obs}$, can be written as \n",
    "\n",
    "$$ \\mathbf{d}^{obs} = \\mathbf{G}\\mathbf{m} + \\mathbf{noise}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "UY71pbGgRaPr7PBxH3OT",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 2
     }
    }
   },
   "source": [
    "##  Step1: Create a model, $\\mathbf{m}$\n",
    "\n",
    "The model $m$ is a function defined on the interval (-2,2). Here we generate a model that is the sum of a: (a) background $m_{ref}$, (b) box car $m_1$ and (c) Gaussian $m_2$. The box car is defined by\n",
    "- `m$_{background}$` : amplitude of the background\n",
    "- `m1` : amplitude\n",
    "- `$m1_{center}$` : center\n",
    "- `$m1_{width}$` : width\n",
    "the Gaussian is defined by \n",
    "- `m2` : amplitude\n",
    "- `$m2_{center}$` : center\n",
    "- `$m2_{sigma}$` : width of Gaussian (as defined by a standard deviation)\n",
    "- `M`: # of model parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "BoaOYiDvY7yFMo9QvOpF",
      "project": "FesuZy9nERkRH9r67vBe",
      "version": 1
     },
     "outputId": {
      "block": "vuewWXXAWRcE2KNdwnEk",
      "project": "FesuZy9nERkRH9r67vBe",
      "version": 1
     }
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b8ec921437494d54ab190cabc88af603",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, continuous_update=False, description='m$_{background}$', max=2.0,…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "app.reset_to_defaults()\n",
    "Q_model = app.interact_plot_model()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "0zUv3gNqulptgg5RQABq",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 4
     }
    }
   },
   "source": [
    "##  Step2: Generate a sensitivity kernel (or matrix), $\\mathbf{G}$\n",
    "\n",
    "By using the following app, we explore each row vector of the sensitivity matrix, $\\mathbf{g}_j$. Parameters of the apps are:\n",
    "\n",
    "- `M`: # of model parameters\n",
    "- `N`: # of data\n",
    "- `p`: decaying constant (<0)\n",
    "- `q`: oscillating constant (>0)\n",
    "- `ymin`: maximum limit for y-axis\n",
    "- `ymax`: minimum limit for y-axis\n",
    "- `show_singular`: show singualr values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "x9PxRUeZbvYukM3fIzi0",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 27
     },
     "outputId": {
      "block": "xjP1IRGiznblm7UECCWU",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 27
     }
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7912e3d288204167ae717f5fa298c456",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(IntSlider(value=20, continuous_update=False, description='N', min=1), IntSlider(value=10…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "app.reset_to_defaults()\n",
    "Q_kernel = app.interact_plot_G()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "B3aFbgUWKTzmqttCPDu3",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 1
     }
    }
   },
   "source": [
    "## Step3: Simulate data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "fOWoINnfZAyspfxAKRnr",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 3
     }
    }
   },
   "source": [
    "The $j$-th datum is the inner product of the $j$-th kernel $g_j(x)$ and the model $m(x)$. In discrete form it can be written as the dot product of the vector $g_j$ and the model vector $m$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "xeKrzyj4YLEYLfraRG1y",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 1
     }
    }
   },
   "source": [
    "### $$ d_j = \\mathbf{g}_j \\mathbf{m} $$\n",
    "\n",
    "If there are $N$ data, these data can be written as a column vector, $\\mathbf{d}$:\n",
    "\n",
    "\\begin{align}\n",
    "    \\mathbf{d} = \\mathbf{G}\\mathbf{m} = \n",
    "    \\begin{bmatrix}\n",
    "        d_1\\\\\n",
    "        \\vdots\\\\\n",
    "        d_{N}\n",
    "    \\end{bmatrix}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "k5SdPFYbjSTR0Unuprld",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 1
     }
    }
   },
   "source": [
    "### Adding Noise\n",
    "\n",
    "Observational data are always contaminated with noise. Here we add Gaussian noise $N(0,\\epsilon)$ (zero mean and standard deviation $\\sigma$). Here we choose \n",
    "\n",
    "$$ \\epsilon = \\% |d| + \\text{floor} $$\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "lSiXnfA4hZvPpbkCW4st",
      "project": "FesuZy9nERkRH9r67vBe",
      "version": 1
     },
     "outputId": {
      "block": "NsHNyjWDNLM2ZFx0quKY",
      "project": "FesuZy9nERkRH9r67vBe",
      "version": 1
     }
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b666c7a20012473d9ee8a647df9ffc13",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(Checkbox(value=False, description='add_noise'), FloatText(value=0.0, description='percen…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "app.reset_to_defaults()\n",
    "Q_data = app.interact_plot_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "Kfas1vohMX9B1w8haaeF",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 1
     }
    }
   },
   "source": [
    "## Step4: All three steps together"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "gpeg40wOdccGxZMMBIxU",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 26
     },
     "outputId": {
      "block": "yFInpwVxV9MsRBuL5pNr",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 26
     }
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e7467249c0c14dd7af088033b97d9a7c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HBox(children=(VBox(children=(Button(description='Model', layout=Layout(height='30px', width='1…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "app.interact_plot_all_three_together()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "I5cRNN1g9FcrGvbH7UKL",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 1
     }
    }
   },
   "source": [
    "## Inverse Problem\n",
    "\n",
    "In the inverse problem we attempt to find the model $\\mathbf{m}$ that gave rise to the observational data $\\mathbf{d}^{obs}$. The inverse problem is formulated as an optimization problem: \n",
    "\n",
    "\n",
    "$$\\text{minimize} \\ \\ \\ \\phi(\\mathbf{m}) = \\phi_d(\\mathbf{m}) + \\beta \\phi_m(\\mathbf{m}) $$\n",
    "\n",
    "where \n",
    "\n",
    "- $\\phi_d$: data misfit\n",
    "- $\\phi_m$: model regularization\n",
    "- $\\beta$: trade-off (or Tikhonov) parameter  $0<\\beta<\\infty$\n",
    "\n",
    "Data misfit is defined as \n",
    "\n",
    "$$ \\phi_d = \\sum_{j=1}^{N}\\Big(\\frac{\\mathbf{g}_j\\mathbf{m}-d^{obs}_j}{\\epsilon_j}\\Big)^2$$\n",
    "\n",
    "where $\\epsilon_j$  is an estimate of the standard deviation of the $j$th datum.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "DPEqdTQucMMV6G5ONR7c",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 1
     }
    }
   },
   "source": [
    "The model regularization term, $\\phi_m$, can be written as \n",
    "\n",
    "$$ \\phi_m(\\mathbf{m}) = \\alpha_s \\int (\\mathbf{m}-\\mathbf{m}_{ref}) dx + \\alpha_x \\int (\\frac{d \\mathbf{m}}{dx}) dx$$\n",
    "\n",
    "The first term is referred to as the \"smallness\" term. Minimizing this generates a model that is close to a reference model $m_{ref}$. The second term penalizes roughness of the model. It is generically referred to as a \"flattest\" or \"smoothness\" term.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "j7hVPYIA73jbxv0bAinV",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 6
     }
    }
   },
   "source": [
    "## Step5: Invert the data, and explore inversion results\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "iooxa": {
     "id": {
      "block": "8k53qHfi1hcwT29TMxNC",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 3
     }
    }
   },
   "source": [
    "In the inverse problem we define parameters needed to evaluate the data misfit and the model regularization terms. We then deal with parameters associated with the inversion.\n",
    "\n",
    "### Parameters\n",
    "\n",
    "- `mode`: `Run` or `Explore`\n",
    "    - `Run`: Each click of the app, will run `n_beta` times of inversion\n",
    "    - `Explore`: Not running inversions, but explore result of the inversions\n",
    "\n",
    "\n",
    "- `noise option`: `error contaminated` or `clean data`\n",
    "\n",
    "#### Misfit\n",
    "- `percent`: percentage of the uncertainty (%)\n",
    "\n",
    "- `floor`: floor of the uncertainty (%)\n",
    "\n",
    "- `chifact`: chi factor for stopping criteria (when $\\phi_d^{\\ast}=N \\rightarrow$ `chifact=1`)\n",
    "\n",
    "#### Model norm\n",
    "- `mref`: reference model\n",
    "\n",
    "- `alpha_s`: $\\alpha_s$ for smallness\n",
    "\n",
    "- `alpha_x`: $\\alpha_x$ for smoothness\n",
    "\n",
    "#### Beta\n",
    "- `beta_min`: minimum $\\beta$\n",
    "\n",
    "- `beta_max`: maximum $\\beta$\n",
    "\n",
    "- `n_beta`: the number of $\\beta$\n",
    "\n",
    "#### Plotting options\n",
    "\n",
    "- `data`: `obs & pred` or `normalized misfit`\n",
    "    - `obs & pred`: show observed and predicted data\n",
    "    - `normalized misfit`: show normalized misfit\n",
    "\n",
    "\n",
    "- `tikhonov`: `phi_d & phi_m` or `phi_d vs phi_m`\n",
    "    - `phi_d & phi_m`: show $\\phi_d$ and $\\phi_m$ as a function of $\\beta$\n",
    "    - `phi_d vs phi_m`: show tikhonov curve\n",
    "    \n",
    "- `i_beta`: i-th $\\beta$ value\n",
    "\n",
    "- `scale`: `linear` or `log`\n",
    "    - `linear`: linear scale for plotting the third panel\n",
    "    - `log`: log scale for plotting the third panel     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### $$\\text{minimize} \\ \\ \\ \\phi(\\mathbf{m}) = \\phi_d(\\mathbf{m}) + \\beta \\phi_m(\\mathbf{m}) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "iooxa": {
     "id": {
      "block": "zYsKlEza7On1cZcJUO1x",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 21
     },
     "outputId": {
      "block": "XHsuSTVmGPtaRtV3WO9m",
      "project": "VNMrkxzChhdveZyf6lmb",
      "version": 21
     }
    },
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "87337f8a284040a2ad713455eae8302e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HBox(children=(RadioButtons(description='mode', options=('Run', 'Explore'), value='Run'), Radio…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "app.interact_plot_inversion()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "iooxa": {
   "id": {
    "block": "pwZt1IzpzYFhHsoFuTHA",
    "project": "FesuZy9nERkRH9r67vBe",
    "version": 1
   }
  },
  "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": 2
}