notebooks/live-skytem-inversion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4d97d6c7-795f-4d02-b2fe-8bde20eb7a60",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# Geophysical inversions with `SimPEG`\n",
    "## An example with airborne electromagnetic data\n",
    "Presenters:\n",
    "Joseph Capriotti, Lindsey Heagy, Seogi Kang"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ccc4efb8-9556-40a3-bd95-ff79295ed020",
   "metadata": {
    "editable": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": "skip"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from utilities.gex_parser import parse_gex_file\n",
    "# If on google colab you will likely need to do instead:\n",
    "# from gex_parser import parse_gex_file\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Start with just importing the time domain module of simpeg, and a utility to plot a layered model.\n",
    "from SimPEG.electromagnetics import time_domain as tdem\n",
    "from SimPEG.utils import plot_1d_layer_model\n",
    "\n",
    "import discretize"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "99341123-e0b6-4dd8-9a0d-c03fe8260533",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Reading in the Data\n",
    "\n",
    "First a bit about the skytem system. There are in general two sets of data:\n",
    "* The low-moment data which is more sensitive to the near surface\n",
    "* The high-moment data, which is (relatively) more sensitive to deeper structures.\n",
    "\n",
    "The system has a hexagonal loop transmitter, with a $\\frac{\\partial \\vec{B}}{\\partial t}$ sensor."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46afbbed-2656-4575-83e5-fdbc2a4a9e87",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "There are two files that represent our data:\n",
    "* A configuration file. Usually a `.gex` extension.\n",
    "* The processed data file.\n",
    "  * Simple CSV file $\\rightarrow$ use `pandas`\n",
    " \n",
    "We've provided a simple parser to read the information in the configuration file into python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7653a6cc-1175-4f7a-989b-0d6968d12dc1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# gex_file = "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f71f8f2-2eb6-413e-ac85-1ddf4b217806",
   "metadata": {
    "editable": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "gex_file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cedba3b2-876a-4a82-8493-f96395dee827",
   "metadata": {},
   "source": [
    "Can simply use `pandas.read_csv` to read in the data file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e184a7d-f2f1-444d-a194-014e7ed73906",
   "metadata": {
    "editable": true,
    "scrolled": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_file = pd.read_csv('../data/MCWD3_dat.xyz')\n",
    "# you can alternatively use pandas to directly fetch this file from the internet (if on colab)\n",
    "# data_file = pd.read_csv(\"https://github.com/simpeg/segns-2024-tutorial/raw/main/data/MCWD3_dat.xyz\")\n",
    "data_file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "798f63dc-158e-4076-ac73-fcd8fa3454fc",
   "metadata": {},
   "source": [
    "inside this data file, missing data are marked with the value 9999, let's replace them with NaNs instead (which the plotting utilities will then ignore)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c250733-360b-4512-acc2-f73b74f3c306",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the data file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a074a99-eed6-4bf0-99c6-b69e34f6055d",
   "metadata": {},
   "source": [
    "Let's plot the location of all the stations in this survey"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05dc7e16-3146-4ddd-8c75-903ca69a4ac3",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.scatter(data_file.UTMX, data_file.UTMY, s=0.5)\n",
    "plt.xlabel('Easting (m)')\n",
    "plt.ylabel('Northing (m)');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a408e7d-8328-44ee-b3e4-3a008e764868",
   "metadata": {},
   "source": [
    "A few bits of processing have already been done to this data set, in particular, they have been normalized for source current strength, normalized by the transmitter area, and have also been spatially averaged along flight lines to form \"stations\".\n",
    "\n",
    "Now we need a bit of information from the `.gex` file to tell us the number of high moment and low moment time gates for this survey"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e955391-8fe4-4a6e-ae5d-f59bf20d608b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# n_lm_gates = \n",
    "# n_hm_gates ="
   ]
  },
  {
   "cell_type": "markdown",
   "id": "009b6aff-35c2-4659-9b84-659d4f725e00",
   "metadata": {},
   "source": [
    "The data file will have `n_lm_gates` of db/dt data, `n_lm_gates` of relative error estimates, then `n_hm_gates` of the high moment db/dt data, followed by `n_hm_gates` relative error estimates."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cfbcd72-2e93-4142-9dac-67a9b959709b",
   "metadata": {},
   "source": [
    "### Zooming in to a single station to work with"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1639f8b7-f867-41de-8e13-95bf65d9f85d",
   "metadata": {},
   "source": [
    "Now let's look at a single station, along some randomly chosen line.\n",
    "\n",
    "I can use panda's to easily group all of the data by a common line number\n",
    "and then get that group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63468d8e-fa81-49c0-8a98-ec73263eb890",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# line_no\n",
    "# line_grouping\n",
    "# line"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1fc4012-c3d9-4888-a2f9-e1fe3a9fb9c2",
   "metadata": {},
   "source": [
    "Let's look at all the data along that line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab16db74-7a91-404a-ae7c-1f545237e70f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "lm_data = line.iloc[:, 9:9+n_lm_gates]\n",
    "hm_data = line.iloc[:, 9+2*n_lm_gates:9+2*n_lm_gates + n_hm_gates]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c765091a-fb27-4df0-9c50-2375dcd7e3dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the line data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f99c2a1-14f9-4283-bf7e-59748831e9fa",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# select a single sounding along that line, and grab all of the data associated with that station\n",
    "# Record 5590\n",
    "# station = \n",
    "\n",
    "station_lm_data = station.iloc[0, 9:9+n_lm_gates].to_numpy()\n",
    "station_lm_std = station.iloc[0, 9+n_lm_gates:9+2*n_lm_gates].to_numpy()\n",
    "station_hm_data = station.iloc[0, 9+2*n_lm_gates:9+2*n_lm_gates + n_hm_gates].to_numpy()\n",
    "station_hm_std = station.iloc[0,  9+2*n_lm_gates + n_hm_gates:].to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd833038-50e0-4c76-881c-ebd375b3c854",
   "metadata": {},
   "source": [
    "Now we have numpy arrays of the recorded data for both the low and high moment data sets, along with their relative standard deviations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3c1aaa4-4f59-4636-b9ae-7f1dd8dcb0ea",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.scatter(data_file.UTMX, data_file.UTMY, s=0.5)\n",
    "# add the line and station highlights\n",
    "\n",
    "plt.xlabel('Easting (m)')\n",
    "plt.ylabel('Northing (m)');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a0578df-ea98-47ae-80c7-8d8ff9db94a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.semilogy(line.UTMY, lm_data)\n",
    "plt.semilogy(line.UTMY, hm_data)\n",
    "# show where the station is along the line using axvline\n",
    "plt.xlabel('Northing');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f340ea12-4eee-4cc1-8351-75f806f374e6",
   "metadata": {},
   "source": [
    "In order to model the data, we still need a few things pieces of information.\n",
    "* The location in time of the measurements\n",
    "* The transmitter waveforms\n",
    "* The transmitter shape\n",
    "* The receiver location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "792242bf-90bf-4751-9d15-d5cddf53e986",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# In general the location of the gates is\n",
    "# gate_time = gate_centers + shift + delay\n",
    "# The Gex file also tells us which gates the respective moments use\n",
    "gate_centers = gex_file['General']['GateTimes']['center']\n",
    "\n",
    "channel_info = gex_file['Channel1']\n",
    "shift = channel_info['GateTimeShift']\n",
    "delay = channel_info['MeaTimeDelay']\n",
    "first_gate = channel_info['RemoveInitialGates']\n",
    "last_gate = channel_info['NoGates']\n",
    "# lm_times =\n",
    "\n",
    "channel_info = gex_file['Channel2']\n",
    "shift = channel_info['GateTimeShift']\n",
    "delay = channel_info['MeaTimeDelay']\n",
    "first_gate = channel_info['RemoveInitialGates']\n",
    "last_gate = channel_info['NoGates']\n",
    "# hm_times ="
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21c0cff3-a2df-49c7-b7fe-60ae2f69090c",
   "metadata": {},
   "source": [
    "Some of the station's data are NaN's (which represent no data), Let's remove those points from our consideration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dec09ae6-f51b-4da4-b52a-c0ff06198608",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "lm_good_data = ~np.isnan(station_lm_data)\n",
    "hm_good_data = ~np.isnan(station_hm_data)\n",
    "\n",
    "lm_times = lm_times[lm_good_data]\n",
    "hm_times = hm_times[hm_good_data]\n",
    "\n",
    "station_lm_data = station_lm_data[lm_good_data]\n",
    "station_lm_std = station_lm_std[lm_good_data]\n",
    "station_hm_data = station_hm_data[hm_good_data]\n",
    "station_hm_std = station_hm_std[hm_good_data]\n",
    "\n",
    "# and of course reset our counts too\n",
    "n_lm_gates = len(lm_times)\n",
    "n_hm_gates = len(hm_times)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79ff21e1-08d3-452a-9600-466c546fac01",
   "metadata": {},
   "source": [
    "and finally plot some decay curves.\n",
    "\n",
    "This is the data we will try to match."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1083dca9-07dc-4d55-9c84-3b352020b982",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# loglog plots of time and data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a453795-c69a-49df-8492-0d2286bdf10c",
   "metadata": {},
   "source": [
    "### Setting up the survey"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed986b97-33ca-4bac-af55-e11148b17201",
   "metadata": {},
   "source": [
    "Now we want represent all of the configuration with `SimPEG` objects.\n",
    "\n",
    "SimPEG's general structure to hold this information is:\n",
    "* Create `Receiver`s (that observe data)\n",
    "* Create `Source`s (and attach receivers to those sources)\n",
    "* Gather all the sources into a `Survey`\n",
    "\n",
    "In the time-domain EM module:\n",
    "\n",
    "* `Receivers`\n",
    "    * Observe at a specific location in space.\n",
    "    * Observe a specific component (Here it will be the vertical component of $\\frac{d\\vec{B}}{dt}$)\n",
    "    * At specific times (the time-gate centers).\n",
    " \n",
    "* `Sources`\n",
    "    * Have specific types (Line Current, Magnetic Dipole, etc.)\n",
    "    * Have a location (or set of locations for Line current)\n",
    "    * Have waveforms\n",
    "    * Are \"listened to\" by receivers.\n",
    "\n",
    "So let's start by grabbing the information we will need from the gex file.\n",
    "\n",
    "We already have the time gate locations. So, let's grab all of the information related to the transmitter.\n",
    "\n",
    "Specifically here, the waveforms of the low and high moment sources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9e56a8b-46b8-4662-843a-4aa7cb420419",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "waves = gex_file['General']['Waveforms']\n",
    "# lm_wave_time\n",
    "# lm_wave_form\n",
    "# hm_wave_time\n",
    "# hm_wave_form\n",
    "# plot them in two subplots"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dee437f1-9b0b-43f0-bf11-30089f030320",
   "metadata": {},
   "source": [
    "The shape of the transmitter , and the offset of the receiver coil."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6436252e-b840-4543-a891-052e010b974a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# We pad the locations of the Tx points with 0, to expand it from a 2D (x,y) pair to a 3D (x, y, z) pair.\n",
    "# We also add another row because we will need to close the transmitter loop.\n",
    "# get the tx shape, pad it and close the loop\n",
    "# tx_shape = \n",
    "\n",
    "# get the receiver offset\n",
    "# rx_offset = \n",
    "\n",
    "# let's plot them"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "968e79a6-2eed-43df-b2cd-be8f97ec39b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# tx_shape cw or ccw?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59ed04b8-f452-4224-9a96-4c4fe291d5ce",
   "metadata": {},
   "source": [
    "Let's give the transmitter and receiver explicit locations by adding the UTM coordinates of the station to the transmitter shape and the receiver offset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ac109d2-2a5f-47bc-b6ea-a7f2ab833310",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# give them an absolute position (mostly just need to set the height here)\n",
    "# tx_loc = tx_shape + ...\n",
    "# rx_loc = rx_offset + ...\n",
    "# also let's grab the area\n",
    "# tx_area = "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3837d38d-4ca6-4282-9ba1-d8218b9fec44",
   "metadata": {},
   "source": [
    "We are often asked how to deal with both the low and high moment data in SimPEG. The answer is actually fairly straightforward as long as you understand the structure of a SimPEG survey, we model them as two seperate sources (who just so happen to be in the same location)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3dcdb315-63f0-4a95-b1d5-7c1eb4c36b5a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Low moment:\n",
    "# rx_lm PointMagneticFluxTimeDerivative\n",
    "\n",
    "# wave_lm PiecewiseLinearWaveform\n",
    "\n",
    "# src_lm LineCurrent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8905e18d-4d69-4057-83e2-433cf1ca07c7",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# high moment\n",
    "# rx_hm\n",
    "\n",
    "# wave_hm\n",
    "\n",
    "# src_hm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5516c9c4-a009-4d27-921d-315c269b4dde",
   "metadata": {},
   "source": [
    "Now we can set up a survey that has two sources!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcf2685a-1f57-4334-99e0-cec6c12171b3",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# srv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c335963d-4268-4d75-965c-f7e71b3a2a19",
   "metadata": {},
   "source": [
    "Now we're ready to set up the forward simulation, of which we will use (for this tutorial), a Layered simulation.\n",
    "\n",
    "It has two physical properties:\n",
    "* The conductivity of each layer (`sigma`)\n",
    "* The thickness of each layer (`thicknesses`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25f0d5c5-d9f1-4069-8c4f-5c6ef8f93756",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#first a simple test\n",
    "thicknesses = []\n",
    "conductivities = [1E-2]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "902f26ee-e701-4000-bb69-ee190aa74cf6",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sim = tdem.Simulation1DLayered(srv, sigma=conductivities, thicknesses=thicknesses)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a725c677-75f4-45a4-9f4c-7fed337fe50a",
   "metadata": {},
   "source": [
    "We can set these properties on the simulation and create data. The only requirement is that the conductivity array must be one longer than the thickness array.\n",
    "\n",
    "**note**: For this simulation the conductivities and thicknesses are defined from the surface downwards."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2c348cf4-5cc2-4e4b-88c7-f366e67761b6",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# sim.thicknesses =\n",
    "# sim.sigma = \n",
    "pre = sim.dpred(None)  # no \"model\" (will get to this more later).\n",
    "\n",
    "pre_lm = pre[:n_lm_gates]\n",
    "pre_hm = pre[n_lm_gates:]\n",
    "\n",
    "plt.loglog(lm_times, -pre_lm)  # negative here to account for simpeg convention of z+ up\n",
    "plt.loglog(hm_times, -pre_hm)\n",
    "plt.loglog(lm_times, station_lm_data*tx_area, color='C0', marker='.', linestyle='')  # un-normalize by area\n",
    "plt.loglog(hm_times, station_hm_data*tx_area, color='C1', marker='.', linestyle='')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82e06706-f60b-4382-b7bb-7f67d760b704",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "# Setup an inversion\n",
    "\n",
    "An inversion has many pieces that we must all setup for simpeg.\n",
    "\n",
    "There are:\n",
    "* `Map`: Objects that tell simpeg what to invert for.\n",
    "* `Data`: a container for the observed data, its standard deviation, and the survey.\n",
    "* `ObjectiveFunctions`, which describe the function we minimize to perform the inversion.\n",
    "    * `DataMisfit`: How we measure the fitness of a model to the observed data.\n",
    "    * `Regularization`: A measure of the simplicity of a model.\n",
    "* `Minimization` Routines: Which method we use to iteratively minize the objective function.\n",
    "* `InverseProblem`: Defines the optimization problem and the minimization routine.\n",
    "* `Directives`: Operations that adjust parameters in the objective function as the inversion proceeds.\n",
    "* `Inversion`: Groups together all of the above to actually run an inversion.\n",
    "\n",
    "We will start by describing parametric inversions for simple 1D models, where the number of layers we invert for is much less than the number of data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f23f060-2af8-43f6-99fd-80420aceebb0",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Maps\n",
    "This is how we tell simpeg what to invert for.\n",
    "\n",
    "Maps define how we go from our inversion model to the physical properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c44e8423-f20d-41a8-a709-1e2d0ef28176",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from SimPEG import maps"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46dcdc32-a1e8-4b8a-9262-40ec4957c7ed",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "for example, this is a simple map that will transform its input according to:\n",
    "$$out = e^{in}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b29e91c-e542-4be5-bcdb-1b89de1caf82",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# exp_map ="
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caca0e71-b1a9-4129-8d1b-7ea05f728504",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "This map is particularly useful for solving for conductivity for two reasons:\n",
    "1) Generally conductivity varies on a logarithmic scale\n",
    "2) We do not need to handle positivity constraints while minimizing\n",
    "\n",
    "We use maps in `SimPEG` to tell simulations what we will be inverting for."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8091876-c007-44e5-8f57-d48a0313ff72",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Data\n",
    "We need to create a container for the data, things to remember about here:\n",
    "\n",
    "* The processed data were normalized by the transmitter area.\n",
    "* There is a sign difference in the convention for SimPEG and the processed data\n",
    "* The standard deviations are actually relative errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2538e797-3ba5-4e27-bfeb-d384beaf4e42",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from SimPEG import data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "372d912c-0a37-4d19-a213-69c3d9dc0e5b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# dobs\n",
    "# rel_err\n",
    "\n",
    "# data_container"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "942faa76-ad60-4ce9-a27c-498eef207d06",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "Since the data container knows about the survey, we can directly index it\n",
    "with sources and receivers to retrieve the data associated with that pair.\n",
    "\n",
    "This is particularly useful to ensure that you have passed your observed data in the order that `SimPEG` expects it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c9537b4-d19d-4c47-ba89-86a5caf505dd",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_container[src_hm, rx_hm]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7835cca0-8417-4e72-897d-5b4fc7ced99d",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "## Parametric Halfspace inversion.\n",
    "\n",
    "With a few key component concepts out of the way, Let's see how we can setup the layered simulation to solve for the best fitting halfspace.\n",
    "\n",
    "Let's setup a new simulation that does this, we can re-use our previous survey object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2e846e1-d411-493c-b8c8-dddd365bdd50",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "sim_inv1 = tdem.Simulation1DLayered(srv, sigmaMap=exp_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aba24a2a-2580-4ac0-a159-29cf5bfdaa6a",
   "metadata": {},
   "source": [
    "When we assign a model (or pass a model to a function), the simulation uses the physical property maps to translate the model to the physical properties"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36c7ac2d-9a94-4ae3-879a-59c4db42297a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# try setting the model and looking at sigma"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f71d3fdb-2042-40c2-bc04-d24fbbc65d1a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "### Objective functions\n",
    "\n",
    "SimPEG treats inversions as minimization functions, as such we need an objective function to minimize. For parametric inversions a data misfit is usually sufficient.\n",
    "\n",
    "$$\n",
    "\\phi = \\phi_d + \\beta \\phi_m\n",
    "$$\n",
    "\n",
    "We will start by defining just a simple data misfit function\n",
    "* It needs to know how to compute data\n",
    "* The data set you want to compare against.\n",
    "\n",
    "It needs to be able to evaluate:\n",
    "$$\n",
    "\\phi_d = |W_d ( d_{obs} - F(m))|^2\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e164c4a7-0d2e-46b3-a6ee-0b290626fc92",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from SimPEG.data_misfit import L2DataMisfit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88e14b5a-3bbb-408a-a8e2-3939a7f95e96",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# phi_d_1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "907c97d3-c17f-4552-83e0-d678e5bf7e76",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "We can evaluate the data misfit of a particular model by calling this object. (It also knows how to do derivative and (approximate) Hessian operations needed for minimization)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fafa2c5d-eee4-480e-b470-50b8d63e678c",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# we can evaluate this function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75c52a31-0d5c-4f63-8cf7-392a288add8f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "If we wanted to know the best fitting halfspace for this data we could simply minimze this function. Let's grab a Gauss Newton minimizer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26527625-3792-45dc-b970-1b7250398440",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from SimPEG import optimization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f06b099-2ec9-45c4-b7eb-724707a96520",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "minimizer = opt = optimization.InexactGaussNewton(\n",
    "    maxIter=10, maxIterLS=20, maxIterCG=10, tolCG=1e-3\n",
    ")\n",
    "# Here, Inexact means we are going to use CG to solve for the step direction."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f676bb1-b89b-4b99-80a0-5523be73dc59",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "Define our inverse problem to be minimized:\n",
    "One odd quirk of a SimPEG minimization is that we must pass a regularization object to our inverse problem. However we can set `beta=0` to ignore it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fb39c04-f952-4ade-9393-479d8ef8dae5",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from SimPEG import regularization, inverse_problem\n",
    "\n",
    "empty_reg = regularization.Smallness(discretize.TensorMesh([1]))\n",
    "# note we needed to pass a mesh that had a single cell\n",
    "# in it because our model has 1 value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "558acc60-b1a0-4287-89ae-17021f7de229",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "#inv_prob_1\n",
    "\n",
    "# sets up, phi_d + 0 * phi_m, meaning the minimizer is only going to act on the data misfit term."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "915f8919-eefc-46d3-bfdc-3ab929c59c0a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from SimPEG import inversion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0ad33b0-891c-44c6-becc-af4c35b1d2d8",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# inv1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3711219b-14b4-484b-b368-20c5478e0880",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Run inversion\n",
    "recovered_model = inv1.run(m_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8da3ad9d-0d74-418b-b606-d4c5a5c68508",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# What was our recovered best fitting halfspace?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cfd94bc-3f9c-48c1-bafc-d9dd02f7e3d2",
   "metadata": {},
   "source": [
    "There's another useful way of computing data from a simulation, using the function\n",
    "`make_synthetic_data` which returns a `Data` object instead of just the `numpy` array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c3cefdf-a8f4-4320-8c88-797d0edc6503",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# This function returns a data object\n",
    "# (so we can easily index it with receivers for plotting)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae669815-ce21-4e61-9f7a-c406a193a609",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def plot_data(data_obj):\n",
    "    plt.loglog(lm_times, -data_obj[src_lm, rx_lm])\n",
    "    plt.loglog(hm_times, -data_obj[src_hm, rx_hm])\n",
    "    \n",
    "    plt.loglog(\n",
    "        lm_times, -data_container[src_lm, rx_lm], color='C0', marker='x', linestyle=\"\"\n",
    "    )\n",
    "    plt.loglog(\n",
    "        hm_times, -data_container[src_hm, rx_hm], color='C1', marker='x', linestyle=\"\"\n",
    "    )\n",
    "plot_data(data_pre)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37d511d2-5cba-42da-9183-7568b02a6ad4",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "clearly a 1D model fits this reasonably well (in general), can we do better?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d0f1c3d-f650-414f-8469-dd5ffe4e34c7",
   "metadata": {},
   "source": [
    "### Multiple layers?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0eea3f98-2055-46ed-a8ae-a703b7ac5eb0",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "How do we solve for a multi (but few) layered model? We need to tell simpeg to invert for both conductivity and thickness.\n",
    "\n",
    "Since we also need to ensure thicknesses are positive, lets use another `ExpMap` for them, but now our `model` consists of both values representing conductivity and thicknesses, so we need to tell SimPEG which parts correspond to each of them. There is a simple helper class to construct these called a `Wires`. It is very simple it says the first `X` values correspond to part 1, the next `Y` correspong to part 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "438eeb5c-003a-4f67-b467-cb1971da9f3f",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_layers = 2\n",
    "# wire_map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2f5881f3-c149-4b28-9a47-d45bd1247de2",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "It contains to projection maps which split the model into its conductivities and thickness portions, we still need to use the `exp_map` to say we want to work in the log domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "053863d4-7a9e-4501-9828-b221cae506ac",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# setup with exponent map\n",
    "# sigma_map\n",
    "# thick_map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6334a912-1b2f-4ce3-8723-c642401fb780",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# let's use the best fitting half-space as our initial model\n",
    "# setup initial model\n",
    "# m_sigma_0_2\n",
    "# m_h_0_2\n",
    "m_0_2 = np.r_[m_sigma_0_2, m_h_0_2]\n",
    "\n",
    "sim_inv2 = tdem.Simulation1DLayered(srv, sigmaMap=sigma_map, thicknessesMap=thick_map)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39bb1129-7bd6-461d-b64c-a28daa07eb3c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The simulation knows how to go from the model to the physical properties\n",
    "sim_inv2.model = m_0_2\n",
    "# sim_inv2. ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37ac4ac7-e316-427c-8111-768425fd8af8",
   "metadata": {},
   "source": [
    "Now we can set up the same pieces as the half-space inversion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a6ef5f3-cb03-4967-932b-99ed94297346",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "phi_d_2 = L2DataMisfit(data=data_container, simulation=sim_inv2)\n",
    "\n",
    "# create all of the same components of the inversion problem as before\n",
    "minimizer2 = opt = optimization.InexactGaussNewton(\n",
    "    maxIter=10, maxIterLS=20, maxIterCG=10, tolCG=1e-3\n",
    ")\n",
    "\n",
    "# again create a reg that has the same input size as the model vector, and then turn it off by setting `beta=0`\n",
    "empty_reg = regularization.Smallness(discretize.TensorMesh([len(m_0_2)]))\n",
    "inv_prob_2 = inverse_problem.BaseInvProblem(\n",
    "    phi_d_2, reg=empty_reg, opt=minimizer, beta=0.0\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "277f0ae4-c671-41eb-af48-4066fa347b57",
   "metadata": {},
   "source": [
    "And run it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a851bd4e-5d54-442f-97c4-355a6be4fd26",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "inv2 = inversion.BaseInversion(inv_prob_2, [])\n",
    "\n",
    "# Run inversion\n",
    "recovered_model_2 = inv2.run(m_0_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfda91e2-e622-46e6-8266-a46b27f52aaf",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_pre_2 = sim_inv2.make_synthetic_data(recovered_model_2)\n",
    "sim_inv2.sigma, sim_inv2.thicknesses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6fc988c-88da-4abe-86c4-f807a0cbb6f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = plot_1d_layer_model(sim_inv2.thicknesses, sim_inv2.sigma)\n",
    "ax.axvline(sim_inv1.sigma, color='C1')\n",
    "ax.set_ylim([50, 0])\n",
    "ax.set_xlim([1E-2, 1E0]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bb7f174-91d5-4a41-83ac-2a511dd950b2",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_data(data_pre_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "040bbfc8-e8c7-413d-9427-e7678400c940",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "We could keep going this route to add more and more layers, but at some point we start to enter the realm of underdetermined problems. E.g. when the number of model parameters exceed the number of observed data. For this we need to think about the regularization objective function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "658f8521-eade-44a9-b75e-8f3a26edc24a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "slide"
    },
    "tags": []
   },
   "source": [
    "## Regularized inversion\n",
    "\n",
    "Let's keep building off of the pieces we had before. Let's discretize our model into many different layers. But this time we will only invert for the conductivity of each layer. For this we need to set up a mesh. This should be a 1D mesh, and it's single dimension would represent depth, aka we define z+ down for this example (Most simpeg is z+ up though)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ec59d6a-beae-4e6a-bc35-c0a56603af24",
   "metadata": {},
   "source": [
    "We're going to define it with a single cell at 1m thickness, then increase the sizes geometrically by a factor of 1.05, for a total of 64 cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f37eb6d2-781a-4a9d-849b-faa57553f7cd",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# We can setup a mesh discretization using some shorthand:\n",
    "# h = [dh, (dh, n_1), (dh, n_2, expansion_factor)\n",
    "# h\n",
    "# mesh"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4230c6f9-30d2-4ae6-a4d5-f6087010d5fd",
   "metadata": {},
   "source": [
    "This time we only need to setup a simulation that inverts for conductivity, but we give it the thickness of each layer.\n",
    "Which would be the `h[0]` property on the mesh, but truncated to 63 cells. (The last layer is interpreted as infinity depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4b5f8d7-8a10-4ac6-b838-4c48d6dd6fdf",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# sim_reg"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bd5f744f-69f4-4e29-8b09-416313a08928",
   "metadata": {},
   "source": [
    "Define yet another data misfit measure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f485549e-9d66-403a-99af-63f295beef3a",
   "metadata": {},
   "outputs": [],
   "source": [
    "phi_d_reg = L2DataMisfit(data=data_container, simulation=sim_reg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7216c7dd-4dbf-48d8-9b0a-5750e54bf672",
   "metadata": {},
   "source": [
    "Now we need to also create a meaningful regularization function.\n",
    "\n",
    "We're going to use something like:\n",
    "$$\n",
    "\\phi_m = \\alpha_v \\int_V m^2 dV + \\alpha_z \\int_V (\\frac{\\partial m}{\\partial z})^2 dV\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7a69975-25be-4253-8bdc-481fed21164a",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# This class creates a function that measures the smallness\n",
    "# and the smoothness of the model.\n",
    "# reg WeightedLeastSquares"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c92e22f-b881-482f-95d4-2c47c67615d3",
   "metadata": {},
   "source": [
    "This class has some good default choices for $\\alpha_s$ and $\\alpha_x$\n",
    "\n",
    "$$\\alpha_s=1$$\n",
    "$$\\alpha_z = min(mesh.h)^2$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69ffcbdb-247d-4e07-9645-e31dae279fdb",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# m_0_reg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f7001f2-8c39-4115-920f-5c53ed4f18ed",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "\n",
    "# create all of the same components of the inversion problem as before\n",
    "minimizer_reg = opt = optimization.InexactGaussNewton(\n",
    "    maxIter=10, maxIterLS=20, maxIterCG=10, tolCG=1e-3\n",
    ")\n",
    "\n",
    "# inv_prob_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44c979e5-d738-4d15-9ac5-a441201520d1",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "inv3 = inversion.BaseInversion(inv_prob_3, [])\n",
    "\n",
    "# Run inversion\n",
    "recovered_model_3 = inv3.run(m_0_reg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c92eeb9e-0563-4a3f-83ed-65078905f260",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "data_pre_3 = sim_reg.make_synthetic_data(recovered_model_3)\n",
    "ax = plot_1d_layer_model(sim_inv2.thicknesses, sim_inv2.sigma)\n",
    "ax.axvline(sim_inv1.sigma, color='C1')\n",
    "plot_1d_layer_model(sim_reg.thicknesses, sim_reg.sigma, ax=ax, color='C2')\n",
    "ax.set_ylim([100, 0])\n",
    "ax.set_xlim([1E-2, 1E0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe0528e3-ab60-4206-bbd3-4f332fcbadeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_data(data_pre_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4d8a0c0-e6d1-4af1-865e-68b2ddfd2302",
   "metadata": {},
   "source": [
    "If we knew our noise level, We usually look for a model such that tit's $\\phi_d(m) <= nD/2$. This comes from the expected value of a chi^2 disrtibutionl."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fe1cc40-1b99-4434-9288-5b6f9d0a132b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "This process of manually adjusting the regularization parameter is a bit tiresome. In SimPEG we can use `Directive`s to automate this process for us.\n",
    "\n",
    "Say we want to decrease the regularization parameter by a factor of 5 every two iterations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "213d1744-3984-4a83-9f6d-b1e6f687e7cd",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from SimPEG import directives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "234c56be-b981-407f-a700-2f6fe6cb8acf",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# beta_cooler BetaSchedule"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a4564fb-ff68-4836-8184-04a7a165a643",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "we also would likely want to keep track of the model and the function values at each iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a35c5d5e-1212-45e5-be2b-a77d43b5b21e",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# save_dict SaveOutputDictEveryIteration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72ef54c8-c877-4f53-bedf-4a5e448278aa",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "source": [
    "We combine these into a list, and pass them to the inversion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93654801-2790-484f-b84c-42882e82bb65",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# create all of the same components of the inversion problem as before\n",
    "# increase the number of iterations from 10 to 20, maybe do a few more CG iterations\n",
    "minimizer_reg = opt = optimization.InexactGaussNewton(\n",
    "    maxIter=10, maxIterLS=20, maxIterCG=10, tolCG=1e-3\n",
    ")\n",
    "\n",
    "# choose a starting beta\n",
    "inv_prob_4 = inverse_problem.BaseInvProblem(\n",
    "    phi_d_reg, reg=reg, opt=minimizer_reg, beta=10\n",
    ")\n",
    "\n",
    "inv4 = inversion.BaseInversion(inv_prob_4, [save_dict, beta_cooler])\n",
    "\n",
    "# Run inversion\n",
    "recovered_model_4 = inv4.run(m_0_reg)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22aaa446-16c7-4037-9117-dda0d41718ac",
   "metadata": {},
   "source": [
    "The returned model here is just the last model from the inversion run, it might not necessarily be the best one though..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "78264d1b-ef41-4f72-81de-81e6a1486b3f",
   "metadata": {},
   "source": [
    "The `save_dict.outDict` dictionary is indexed per iteration, returning another dictionary containing the parameters and function evaluations. Let's collect all of the $\\phi_d$, $\\phi_m$ and $\\beta$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29513afe-85bd-4c39-8577-f5f376d0702b",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "n_iter = save_dict.opt.iter\n",
    "phi_ds = [save_dict.outDict[i]['phi_d'] for i in range(1, n_iter)]\n",
    "phi_ms = [save_dict.outDict[i]['phi_m'] for i in range(1, n_iter)]\n",
    "betas =  [save_dict.outDict[i]['beta'] for i in range(1, n_iter)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "383e7c41-7ded-4f79-8d58-33ba9481adb2",
   "metadata": {},
   "source": [
    "We can make a few plots to investigate how these function evaluations change per iteration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36b27ed3-cc43-4a0a-b08b-7aff531a1b14",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Make some plots of phi_d, phi_m and beta\n",
    "# choose a \"final\" model\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ae33479-be31-456f-8e71-458d10bac6e4",
   "metadata": {},
   "source": [
    "And what does the model and data look like at this point?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c71501c-1a07-4041-b0f2-e6b4eb807119",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# m_final\n",
    "d_final = sim_reg.make_synthetic_data(m_final)\n",
    "ax = plot_1d_layer_model(sim_reg.thicknesses, exp_map * m_final)\n",
    "plot_1d_layer_model(sim_inv2.thicknesses, sim_inv2.sigma, ax=ax)\n",
    "ax.set_ylim([100, 0])\n",
    "ax.set_xlim([1E-2, 1E0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba404fc7-e7fb-4539-930c-b9e3926af9ba",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "plot_data(d_final)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27c72d17-991c-453a-babe-1c42cd67421c",
   "metadata": {},
   "source": [
    "**note** The overall fit looks fairly close."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05ef7c09-7783-4f45-8e39-41031df55d06",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "subslide"
    },
    "tags": []
   },
   "source": [
    "You might ask..\n",
    "\"What is the noise level corresponding to this regularization parameter?\"\n",
    "There's a few ways you can define it, let's compare the relative sizes of the data residual vector to the observed data vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35704fe8-d9f5-40ea-bd31-3a1567f9c739",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": "fragment"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# rel_diff"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}