In the previous notebooks, we explored how to discretize the 1D MT problem, how to design an appropriate mesh, and explored the an example of non-uniquness. In this notebook, we will put those concepts to use and set up a 1D inversion.

The aim of this notebook is to highlight the elements we use to set up and solve the inverse problem, in the notebook that follows this, we will dive futher into the impact of various parameter choices (eg. the trade-off parameter $\beta$, the stopping criteria, noise model and regulariztion parameters).

We will use deterministic approach and pose the inverse problem as an optimization problem of the form

$$ \min_{\mathbf{m}} \phi(\mathbf{m}) = \phi_d(\mathbf{m}) + \beta\phi_m(\mathbf{m}) $$where:

- $\mathbf{m}$ is our model - the vector of numbers that describes our earth model. It is what we are solving for in the inversion
- $\phi_d(\mathbf{m})$ is the data misfit
- $\phi_m(\mathbf{m})$ is the regularization
- $\beta$ is a trade-off parameter

The data misfit, $\phi_d(\mathbf{m})$, is often taken to be a weighted $\ell_2$-norm, where the weights capture the noise model (eg. we want to assign higher weights and do a good job fitting data that we are confident are less noisy, and assign less weight / influence to data that are noisy). The $\ell_2$ norm is the correct norm to choose when noise is Gaussian (or approximately Gaussian, or if you have no additional information and assume it is Gaussian). An $\ell_2$ data misfit is captured mathematically by

$$ \phi_d(\mathbf{m}) = \frac{1}{2}\|\mathbf{W_d} (\mathcal{F}(\mathbf{m}) - \mathbf{d}^{\text{obs}})\|^2 $$where

- $\mathbf{W_d} = \text{diag}(1/\boldsymbol{\epsilon})$ is a diagonal matrix with diagonal entries $1/\epsilon_i$, where $\epsilon_i$ is an estimated standard deviation of the $i$th datum.
- $\mathcal{F}$ is the forward modelling operator that simulates the 1D MT response given a model
- $\mathbf{m}$ is the model: the set of parameters which describe the physical properties of the layered earth
- $\mathbf{d}^{\text{obs}}$ is the observed data

(The factor of $1/2$ and the choice of using the squared-norm is a matter of convienence. We will be using gradient based optimzation methods, so it is easier to take derivatives of a norm squared than a norm.)

The inverse problem is an ill posed problem. There are multiple (actually infinitely many!) models that can fit the data. There are a couple ways to observe this:

Mathematically: If we start by thinking about a linear problem $\mathbf{G}\mathbf{m} = \mathbf{d}$, the matrix $\mathbf{G}$ is wide, so it is not directly invertible (eg. see Matt Hall's Linear Inversion Tutorial). Here, we are dealing with a non-linear system of equations, but the principle is the same.

An Example: In the Forward Modelling and Nonuniqueness notebook, we used forward modelling to demonstrate non-uniquness with the conductivity-thickness product of a single, conductive layer, and this is a very simple model compared to most geologic settings!

Thus, to choose from the infinitely many solutions and arrive at a sensible one, we employ a regularization: $\phi_m$. Tikhonov regularization, which again employs $\ell_2$-norms, is a standard choice (It has a few nice features: it is convex and easy to differentiate). It takes the form: $$ \phi_m(\mathbf{m}) = \frac{1}{2}\big(\alpha_s\|\mathbf{W_s} (\mathbf{m} - \mathbf{m}_{\text{ref}})\|^2 + \alpha_z\|\mathbf{W_z} (\mathbf{m})\|^2 \big) $$

The first term is often referred to as the "smallness" as it measures the "size" of the model (in the $\ell_2$ sense). The matrix $\mathbf{W_s}$ is generally taken to be a diagonal matrix that may contain information about the length scales of the model or be used to weight the relative importance of various parameters in the model. The scalar $\alpha_s$ weights the relative importance of this term in the regularization. Notice that we include a reference model ($\mathbf{m}_{\text{ref}}$. Often this is defined as a constant value, but if more information is known about the background, that can be used to construct a reference model. Note that saying "I am not going to use a reference model" means that you are actually using $\mathbf{m}_{\text{ref}} = 0$, this is important to realize... in the inversion we demonstrate here, our `model`

will be $\mathbf{m} = \text{log}(\sigma)$. If we set $\mathbf{m}_{\text{ref}} = 0$, then we are favoring models close to 1 S/m - which is quite conductive!

The second term is often referred to as the "smoothness". The matrix $\mathbf{W_z}$ approximate the derivative of the model with respect to depth, and is hense a measure of how "smooth" the model is. The term $\alpha_z$ weights its relative importance in the regularization.

Although we pose the inverse problem as an optimization problem, we aren't necessarily going to solve to a true minimum. There are always uncertainties with the data, so there is no point trying to drive $\phi_d$ to its true minimium - we would end up introducing structures in the model to fit the data (we will explore this here). So another choice that has to be made when setting up an inversion is where to stop.

In this notebook, we will walk through how to set up an inversion and discuss one of the knobs:

- The $\beta$ knob: How well should we fit the data verses the regularization? Which $\beta$ value provides us the right balance? eg. When do we achieve the target misfit $\phi_d^*$?
- The $\alpha$ knob: Tuning the
**regularization**($\alpha_s$, $\alpha_z$): how*smooth*or*small*(in the $\ell_2$-sense) do we want our model to be?

We will also point out other parameters you can explore the impact of, including,

- starting model $\mathbf{m}_0$
- reference model $\mathbf{m}_{\text{ref}}$
- designing a $\beta$-cooling schedule

In the next notebook we will dive further into the $\alpha$ knob.

**A Note**

We will employ second-order optimization methods, meaning we have to be able to both simulate data and compute derivatives. We will take those steps for granted in this notebook, but if you are curious to look under the hood and see how sensitivities are computed for this problem, have a look at the Sensitivities Notebook.

In [1]:

```
from scipy.constants import mu_0
import numpy as np
import matplotlib.pyplot as plt
from SimPEG import (
Mesh, Maps, SolverLU, DataMisfit, Regularization,
Optimization, InvProblem, Inversion, Directives, Utils
)
try:
from pymatsolver import PardisoSolver as Solver
except importError:
from SimPEG import SolverLU as Solver
from MT1D import MT1DProblem, MT1DSurvey, MT1DSrc, ZxyRx, Survey, AppResPhaRx
%matplotlib inline
```

Now, we define our model parameters and survey setup. This includes defining

- location of receivers
- frequencies we are sampling
- type of dat we are measuring (impedance or apparent resistivity and phase)

Similar to the model shown in 1_MT1D_NumericalSetup.ipynb and 2_MT1D_ForwardModellingAndNonuniqueness. We will consider a model which consists of 5 units:

In [2]:

```
layer_tops = np.r_[0., -600., -1991., -5786., -9786.] # in m
rho_layers = np.r_[250., 25, 100., 10., 25.]
```

In [3]:

```
rxloc = np.r_[0.]
frequency = np.logspace(-3, 2, 25)
```

Next, we set up a `survey`

object. Here we are following the SimPEG approach and define

- a receiver object: given fields found from a simulation, it has the machinery to compute data. An impedance receiver
`ZxyRx`

calculates the ratio of $E_x$ and $H_y$ from calculated electric and magnetic fields (both real and imaginary components)

- a source object: this provides the right-hand-side for our MT problem. This is a plane wave source with at the specified frequency.

These are combined in a `survey`

.

In [4]:

```
# Create a receiver object
rx = ZxyRx(
rxloc, # location of the receiver
component="both", # measure both the real and imaginary components of the impedance (alternatively "real" / "imag")
frequency=frequency
)
# create a plane wave source
src = MT1DSrc([rx])
# define a survey
survey = MT1DSurvey([src])
```

In the Forward Modelling and Nonuniqueness notebook, we discussed how to design a mesh that extends sufficiently far and has fine enough cells near the surface to accurately simulate the MT response across the frequency range of interest. We have wrapped up that knowledge in the utility function `setMesh`

.

Since most of the geologic units we are considering are on the order of $100 \Omega m$, we will use a conductivity value of $10^{-2}$ S/m for creating the mesh.

In [5]:

```
max_depth_core = 15000.
mesh = survey.setMesh(
sigma=0.01, # approximate conductivity of the background
max_depth_core=max_depth_core, # extent of the core region of the mesh
ncell_per_skind=10, # number of cells per the smallest skin depth
n_skind=2, # number of skin depths that the mesh should extend to ensure the lowest-frequency fields have decayed
core_meshType = "log", # cell spacings in the core region of the mesh ("linear" or "log")
max_hz_core=1000. # If using a logarithmic core mesh, what is the maximum cell size?
)
```

In [6]:

```
# plot the mesh
fig, ax = plt.subplots(1,1, figsize=(8, 3))
mesh.plotGrid(centers=True, faces=True, ax=ax)
ax.legend(["centers", "faces"])
ax.invert_xaxis() # so that the surface is on our left hand side
ax.set_xlabel('z (m)')
ax.grid(which="both", linewidth=0.5)
```

In the first notebook, we set up the machinery to solve the 1D MT problem. Here, we will use that functionality. Please see Sensitivities Notebook for further discussion on how `MT1DProblems`

is set up.

In [7]:

```
prob = MT1DProblem(
mesh, # The mesh contains the geometry, grids, etc necessary for constructing the discrete PDE system
sigmaMap=Maps.ExpMap(mesh), # in the inversion, we want to invert for log-conductivity (enforces postivity, electrical conductivity tends to vary logarithmically)
verbose=False, # print information as we are setting up and solving
Solver=Solver # solver to employ for solving Ax = b
)
# tell the problem and survey about each other so we can construct our matrix system
# and right hand-side
prob.pair(survey)
```

Physical properties are defined at cell centers, so there should be a `sigma`

value for every cell center. Above, we defined our model as resisvitities (in units of $\Omega m$), but for the inversion, we want to work in conductivities (in units of S/m), so we take the reciprocal ($\sigma = 1/\rho$)

In [8]:

```
# start with nans so we can do a check to make sure all
# layer conductivities have been properly assigned
rho = np.ones(mesh.nC) * np.nan
# loop over each layer in the model and assign to mesh
for layer_top, rho_layer in zip(layer_tops, rho_layers):
inds = mesh.vectorCCx < layer_top
rho[inds] = rho_layer
sigma = 1./rho
```

In [9]:

```
# plot the model
fig, ax = plt.subplots(1, 1, figsize=(8, 3))
# trickery to plot from node to node rather than at cell centers
z = np.repeat(mesh.vectorCCx[1:], 2, axis=0)
z = np.r_[mesh.vectorCCx[0], z, mesh.vectorCCx[-1]]
sigma_plt = np.repeat(sigma, 2, axis=0)
ax.loglog(-z, sigma_plt, lw=2)
ax.set_ylabel("Conductivity (S/m)", fontsize=13)
ax.set_xlabel("Depth (m)", fontsize=13)
ax.grid(True, which='both', linewidth=0.4)
ax.set_ylim(2e-3, 2e-1)
ax.set_xlim(0, max_depth_core*2.)
```

Out[9]:

Based on the conductivity structure plotted above, we will create synthetic data. These data will later be used in the inversion.

In the inversion, the "model" that we will invert for is in $\log(\sigma)$. There are a couple of reasons for this: the electrical conductivity of earth materials can vary over several orders of magnitude, and electrical conductivity is always positive, so inverting for $\log(\sigma)$ enforces that.

In SimPEG, the method `dpred`

of the `survey`

class solves the PDE and computes the data required by the receivers.

In [10]:

```
mtrue = np.log(sigma) # since our "model" is log conductivity, we take the log
dtrue = survey.dpred(mtrue) # these are clean data (no noise yet.)
```

Add noise to generate "observed" data

In [11]:

```
np.random.seed(1) # set a seed to the results are reproducable
std = 0.1 # standard deviation of the noise (10%)
# add noise
uncert = std * np.abs(dtrue)
noise = uncert * np.random.randn(survey.nFreq*2)
survey.dobs = dtrue + noise
```

Although we will work with real and impaginary impedance values as our data in the inversion, it is a bit more intuitive to plot the data in terms of apparent resistivity (which hase units of resistivity, $\Omega m$) and phase

$$ \rho_a = \frac{1}{\mu_0\omega} \big|Z_{xy}\big|^2 $$$$ \phi = \tan^{-1}\left(\frac{\text{Im}(Z_{xy})}{\text{Re}(Z_{xy})}\right) $$In [12]:

```
def omega(frequency):
"""
angular frequency
"""
return 2*np.pi*frequency
def appres_phase_from_data(data, frequency):
"""
Compute apparent resistivity and phase given impedances (real and imaginary components)
and the frequency.
"""
# data are arranged (Zxy_real, Zxy_imag) for each frequency
Zxy_real = data.reshape((survey.nFreq, 2))[:,0]
Zxy_imag = data.reshape((survey.nFreq, 2))[:,1]
Zxy = Zxy_real+1j*Zxy_imag
# compute apparent resistivity and phase from complex impedance
app_res = abs(Zxy)**2 / (mu_0*omega(frequency))
phase = np.rad2deg(np.arctan(Zxy_imag / Zxy_real))
return app_res, phase
```

In [13]:

```
# fetch the apparent resistivity and phase for the clean (true)
# and noisy (obs) data
app_res_true, phase_true = appres_phase_from_data(dtrue, frequency)
app_res_obs, phase_obs = appres_phase_from_data(survey.dobs, frequency)
```

In [14]:

```
fig, ax = plt.subplots(2, 1, figsize=(8,6))
# apparent resistivity
ax[0].loglog(frequency, app_res_true, '-k', lw=1, label="clean")
ax[0].loglog(frequency, app_res_obs, '-o', lw=1, label="noisy")
ax[0].set_ylabel("$\\rho_a \ (\Omega m)$", fontsize = 14)
ax[0].set_ylim([10, 3e2])
# phase
ax[1].semilogx(frequency, phase_true, '-k', lw=1, label="clean")
ax[1].semilogx(frequency, phase_obs, '-o', lw=1, label="noisy")
ax[1].set_ylabel("Phase (degree)")
ax[1].set_ylim([30, 80])
for a in ax:
a.grid(True, which='both', linewidth=0.4)
a.legend()
a.set_xlim(frequency.max(), frequency.min())
a.set_xlabel("Frequency (Hz)")
plt.tight_layout()
```

Our setup of the inversion follows the SimPEG framework.

The "inversion implementation" consists of 8 modules, 3 of which we have already been working with. This is a very brief overview. For more details, see the SimPEG docs

`Mesh`

: mesh geometry and differential operators`Problem`

: physics engine. contains the machinery to construct and solve the PDE`Survey`

: sources and receivers`Data Misfit`

: measure of how "far" the predicted data are from the observed data`Regularization`

: Regularization on the model. Here we use Tikhonov regularization`InvProb`

: statement of the inverse problem, brings together the data misfit and regularization to define an objective function that we minimize in the inversion`Optimization`

: algorithm we will use to perform the optimization in the inversion, here, we choose a gradient based approach`Inversion`

: Bring everything together to run the inversion. This includes`directives`

which are instructions on updates that should be made during the course of the inversion (eg. cooling beta) and stopping criteria (eg. target misfit)

In [15]:

```
def run_MT1Dinv(
prob, # 1D MT problem
survey, # 1D MT survey with sources and receivers
m0, # starting model
mref=None, # reference model
alpha_s=1., # smallness weight
alpha_z=1., # smoothness weight
beta0=1e1, # trade off parameter
use_betaest=False, # estimate an initial beta based on the the data misfit and regularization terms
beta0_ratio=None , # if we are estimating beta, how much should we favor the regularization?
coolingFactor=1.5, # cooling factor
coolingRate=1, # cool beta after this many iterations
use_target=True, # stop the inversion at the target misfit?
):
# if the starting model is not defined, use the reference model
if mref is None:
mref = m0
# Data misfit
dmisfit = DataMisfit.l2_DataMisfit(survey)
dmisfit.W = 1./uncert
# Regularization
reg = Regularization.Simple(
prob.mesh, alpha_s=alpha_s, alpha_x=alpha_z, mref=mref
) # since we are in 1D, we work with the first dimension
# Optimization
opt = Optimization.InexactGaussNewton(maxIter=35, LSshorten=0.05)
# Statement of the inverse problem
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
# Inversion Directives
beta = Directives.BetaSchedule(
coolingFactor=coolingFactor, coolingRate=coolingRate
)
invProb.beta = beta0
invProgress = Directives.SaveOutputEveryIteration(save_txt=False)
target = Directives.TargetMisfit()
directives = [beta, invProgress]
if use_target:
directives.append(target)
if use_betaest:
if beta0_ratio is None:
beta0_ratio = 1.
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=beta0_ratio)
directives.append(betaest)
# assemble in an inversion
inv = Inversion.BaseInversion(invProb, directiveList=directives)
prob.counter = opt.counter = Utils.Counter()
opt.remember('xc')
# run the inversion
mopt = inv.run(m0)
xc = opt.recall("xc")
return xc, invProgress, reg
```

If the noise is Gaussian, then the sum of squares (our data misfit) is a Chi-squared distribution, which has an expected value of $N_\text{data}$ (in our case, we divide this by two to match our definition of $\phi_d$). Thus, the ideal choice of $\beta$ is one that gives us $$\phi_d^* \approx \frac{1}{2} N_\text{data}$$

You can choose to stop the inversion once it achieves the target misfit (`use_target=True`

). For this example, we want to also demonstrate what happens when you over-fit the data, so we will not stop it at the target misfit, but run the inversion until it hits the maximum number of iterations.

In this example, we will run an inversion where we cool (reduce the value of) `beta`

every iteration. In practice, `beta`

cooling strategy is often applied in non-linear inversions. The data misfit contribution objective function is non-convex, meaning that it has more than one local minimum. The Tikhonov regularization is convex. Having a large `beta`

ensures that at least locally, the function is convex, so we can make progress using a gradient-based optimization step. However, keeping `beta`

very high means that we do a good job fitting the regularization, and not necessarily fitting the data. Adopting a cooling schedule is one strategy for trying to address these competing challenges: by starting with a large beta and gradually reducing it, we make progress on a locally convex problem and then gradually allow more structure to the model by gradually reducing the value of `beta`

. Typically multiple model updates are made for each beta-value (eg. you can increase the `coolingRate`

to 3 to have 3 model update steps taken for each value of beta). The `coolingFactor`

is the factor by which we reduce beta at the next beta-iteration (eg. `beta[i+1] = beta[i]/coolingFactor`

).

The `use_betaest`

selects if we estimate a first beta or prescribe it. The directive which estimates the first `beta`

value approximates the relative magnitudes of both $\phi_d$ and $\phi_m$ (eg. if $\phi_d \sim 10$ and $\phi_m \sim 1$, then to balance their influence in the objective function, we would choose $\beta = 10$). The `beta_factor`

is a scaling factor so that we can specify if we would like the influence of $\phi_m$ to be `beta_factor`

times larger than $\phi_d$ (eg. if $\phi_d \sim 10$ and $\phi_m \sim 1$ and we use `beta_factor = 2`

, then the initial beta would be 20).

For the following example, we will gradually reduce $\beta$ at every iteration - this allows us to plot out a "typical" Tikhonov or trade-off curve of $\phi_d$ vs. $\phi_m$. We encourage you to explore different beta-cooling strategies.

Here, we perform a "smooth" inversion, that is, we define $\alpha_s \ll \alpha_z$, and penalize large gradients in the model. What happens if you reverse the values of `alpha_s`

and `alpha_z`

and instead perform a "small" inversion?

For both a starting and reference model, we use the background conductivity of $10^{-2}$ S/m. Since in the inversion,the model is $\log(\sigma)$, we take the logarithm. What happens when you change these?

In [16]:

```
sigma_ref = 1e-2 # reference conductivity
sigma_0 = 1e-2 # starting conductivity
# translate the starting and reference model to log-conductivity
mref = np.log(sigma_ref)*np.ones(mesh.nC)
m0 = np.log(sigma_0)*np.ones(mesh.nC)
xc, invProgress, reg = run_MT1Dinv(
prob,
survey,
m0=m0, # starting model
mref=mref, # reference model
alpha_s=1e-2, # smallness contribution
alpha_z=1., # smoothness contribution to the regularization
use_betaest=True, # estimate the initial beta
beta0_ratio=10., # starting contribution of regularization 10x larger than the data misfit
coolingFactor=1.5, # how much should we reduce beta at every beta-iteration
coolingRate=1, # reduce beta after each `coolingRate` iterations
use_target=False # stop the inversion at the target misfit?
)
```

In [17]:

```
np.savetxt('mt_inv.txt', np.vstack([mesh.vectorCCx, xc[12]]).T)
```

In the following plots, we show

- $\phi_d$ vs. $\beta$
- $\phi_m$ vs. $\beta$
- $\phi_d$ vs. $\phi_m$

Notice that as $\beta$ decreases, $\phi_d$ decreases and $\phi_m$ increases (as we decrease the importance of the model regularization, we introduce more structure in the model and do a better job fitting the data)

In all plots, the star indicates the iteration where we achieved target misfit.

In [18]:

```
invProgress.plot_tikhonov_curves()
```

In [19]:

```
invProgress.plot_misfit_curves()
```

In what follows, we will generate an interactive plots that show

(a) the true and recovered model

(b) the Tikhonov curve with the target misfit plotted as a star, and the current iteration as a dot

(c) predicted and observed apparent resistivity curve

(d) predicted and observed phase curve

What does our data fit look like at each iteration? When have we underfit the data? overfit the data? Let's scroll through the inversion results and see.

In [20]:

```
# grab beta, phi_d, phi_m from the inversion. They are lists, so we convert them to arrays
beta = np.array(invProgress.beta)
phi_d = np.array(invProgress.phi_d)
phi_m = np.array(invProgress.phi_m)
# find the iteration where we achieved target misfit
i_target = invProgress.i_target
print("The inversion reached target misfit after {} iterations".format(i_target))
```

In [21]:

```
from ipywidgets import interact, IntSlider, ToggleButtons, fixed
```

In [22]:

```
def view_1Dinversion_results(iteration, scale="linear", color="C0", alpha=1, ax=None, plotTrue=True):
if ax is None:
fig, ax = plt.subplots(2, 2, figsize=(12, 6))
ax = ax.flatten()
# get the apparent resistivity and phase data for this iteration
dpred = survey.dpred(xc[iteration])
app_res_pred, phase_pred = appres_phase_from_data(
dpred, frequency
)
# plot the true and recovered models
m_iter = xc[iteration]
sigtrue = np.repeat(sigma, 2, axis=0)
z = np.repeat(mesh.vectorCCx[1:], 2, axis=0)
z = np.r_[mesh.vectorCCx[0], z, mesh.vectorCCx[-1]]
if plotTrue is True:
ax[0].loglog(-z, sigtrue, 'k', lw=1, label='true')
ax[2].loglog(phi_m, phi_d, 'k-', lw=2)
ax[2].loglog(phi_m[i_target], phi_d[i_target], 'C3*', ms=16, label="target")
ax[1].loglog(frequency, app_res_obs, '.k-', linewidth=1, label='obs')
ax[3].semilogx(frequency, phase_obs, '.k-', linewidth=1, label="obs")
ax[0].loglog(-mesh.vectorCCx, np.exp(m_iter), color=color, alpha=alpha, lw=2, label='rec it {}'.format(iteration))
ax[0].set_ylabel("Conductivity (S/m)", fontsize = 14)
ax[0].set_xlabel("Depth (m)", fontsize = 14)
ax[0].set_ylim(2e-3, 2e-1)
ax[0].set_xlim((-z).min(), (-z).max())
ax[0].legend(loc=2)
# plot the Tikhonov curve
ax[2].loglog(phi_m[iteration], phi_d[iteration], color+'o', ms=8, label="it {}".format(iteration))
ax[2].set_xlim(phi_m.min(), phi_m.max())
ax[2].set_xlabel("$\phi_m$", fontsize = 14)
ax[2].set_ylabel("$\phi_d$", fontsize = 14)
ax[2].set_xscale(scale)
ax[2].legend()
# add iteration and beta values
ax[2].text(
phi_m[iteration], phi_d[iteration]*1.25,
"$\\beta$ = {:1.1e}".format(
beta[iteration]
), fontsize=11, #, xy=(0, 0), xytext=(0, 0)
)
# plot the apparend resistivity data
ax[1].loglog(frequency, app_res_pred, color, alpha=alpha, label='pred, it {}'.format(iteration))
ax[1].set_xlim(frequency.max(), frequency.min())
ax[1].set_ylim(10, 3e2)
ax[1].set_xlabel("Frequency (Hz)", fontsize = 14)
ax[1].set_ylabel("$\\rho_a \ (\Omega m)$", fontsize = 14)
ax[1].legend(loc=1)
# plot the phase
ax[3].semilogx(frequency, phase_pred, color, alpha=alpha, label="pred, it {}".format(iteration))
ax[3].set_xlim(frequency.max(), frequency.min())
ax[3].set_xlabel("Frequency (Hz)", fontsize = 14)
ax[3].set_ylabel("Phase (degree)", fontsize = 14)
ax[3].legend(loc=1)
for a, title in zip(ax, ['(a)', '(c)', '(b)', '(d)']):
a.set_title(title, fontsize=14)
if title != "c":
a.grid(True, which='both', linewidth=0.4)
```

In [23]:

```
# for ipywidets, all kwarg inputs must be specified, and plt.show() needs to be used
def interact_view_inversion(iteration, scale):
view_1Dinversion_results(iteration, scale)
plt.tight_layout()
plt.show()
interact(
interact_view_inversion,
iteration=IntSlider(min=0, max=beta.shape[0]-1, step=1, value=0),
scale=ToggleButtons(options=["linear", "log"], value="log"),
)
```

Out[23]:

Here is the figure from the article.

In [24]:

```
fig, ax = plt.subplots(2, 2, figsize=(12, 6), dpi=400)
ax = ax.flatten()
for it, col, alpha, plotTrue in zip(
[5, 9, 21], [ "C1", "C0", "C2"], [0.5, 1., 0.5], [True, False, False]
):
view_1Dinversion_results(
iteration=it, scale="log", color=col, alpha=alpha, ax=ax, plotTrue=plotTrue
)
plt.tight_layout()
```