# Computing resolution properties#

This example demonstrates how block and smooth TEM modelling is done and how to access resolution matrices and kernels.

```
```

We import the necessary libraries

```import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.physics import em
from pygimli.viewer.mpl import drawModel1D
```

Next, we define a time vector from 100µs to 10ms and assume a loop of 100x100m. We use the TDEMBlockModelling operator to generate synthetic data for a three-layer case and plot it. Note that the output is apparent resistivity. We could now setup a Marquard-type block inversion as done in the VES example.

```t = np.logspace(-4, -2, 20)
fopkw = dict(times=t, txArea=100**2)
synth = [50, 100, 300, 30, 300]
fopSynth = em.TDEMBlockModelling(nLayers=(len(synth)+1)//2, **fopkw)
data = fopSynth(synth)
plt.loglog(data, t, "x-")
plt.grid()
```

Instead, we want to do a smooth (Occam) style inversion where the layer thicknesses are known. For this, we use the TDEMSmoothModelling operator and logarithmic transforms for both model and data.

```thk = np.logspace(0.7, 1.7, 15)
print(sum(thk))
fop = em.TDEMSmoothModelling(thk=thk, **fopkw)
inv = pg.Inversion(fop=fop)
inv.dataTrans = pg.trans.TransLog()
inv.modelTrans = pg.trans.TransLog()
model = inv.run(data, relativeError=0.03, verbose=True)
```
```302.43843861107877
fop: <pygimli.physics.em.tdem.TDEMSmoothModelling object at 0x7ff746ed6cf0>
Data transformation: <pygimli.core._pygimli_.RTransLog object at 0x7ff747d11fd0>
Model transformation: <pygimli.core._pygimli_.RTransLog object at 0x7ff747d133d0>
min/max (data): 52.5/102
min/max (error): 3%/3%
min/max (start model): 68.23/68.23
--------------------------------------------------------------------------------
inv.iter 0 ... chi² =   55.52
--------------------------------------------------------------------------------
inv.iter 1 ... chi² =   23.79 (dPhi = 56.76%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 2 ... chi² =    8.46 (dPhi = 60.27%) lam: 20.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² =    0.41 (dPhi = 77.27%) lam: 20.0

################################################################################
#                  Abort criterion reached: chi² <= 1 (0.41)                   #
################################################################################
```

The inverted model is plotted along with the synthetic model, it resembles a smoothed image of it.

```fig, ax = plt.subplots()
drawModel1D(ax, model=synth, label="synthetic")
drawModel1D(ax, plot="semilogx", thickness=thk, values=model, label="inverted")
ax.grid(which='minor')
_ = ax.legend()
```

We are now interested in the resolution properties of the inverse problem. We call the function resolution matrix and obtain both the model resolution and the data resolution matrices. For only one, there are the functions modelResolutionMatrix and dataResolutionMatrix.

```from pygimli.frameworks.resolution import resolutionMatrix
RM, RD = resolutionMatrix(inv, returnRD=True)
z = np.hstack([np.cumsum(thk)-thk/2, sum(thk)])
```

We display the model resolution matrix using imshow.

```fig, ax = plt.subplots()
im = ax.imshow(RM, vmin=-0.3, vmax=0.3, cmap="RdBu_r")
ticks = np.arange(0, len(z), 2)
labels = [f"{zi:.0f}" for zi in z[ticks]]
ax.set_xticks(ticks)
ax.set_xticklabels(labels)
ax.set_yticks(ticks)
ax.set_yticklabels(labels)
cb = plt.colorbar(im)
```

Except the last layer, the highest values occur on the main diagonal with maximum values at depths of about 100-200m. We can plot this formal resolution

```fig, ax = plt.subplots()
drawModel1D(ax, plot="plot", thickness=thk, values=np.diag(RM), label="resolution")
```

A single resolution kernel can be obtained by extracting a column of the matrix. For large-scale inverse problems, the model size might be too big to compute the whole matrix. One can then compute single resolution kernels by solving an inverse problem using modelResolutionKernel.

```from pygimli.frameworks.resolution import modelResolutionKernel
fig, ax = plt.subplots()
nr = 10
drawModel1D(ax, plot="plot", thickness=thk, values=RM[:, nr], label="column")
knl = modelResolutionKernel(inv, nr)
drawModel1D(ax, plot="plot", thickness=thk, values=knl, label="kernel", ls="--")
ax.legend()
ax.hlines(z[nr], *ax.get_xlim())
```
```<matplotlib.collections.LineCollection object at 0x7ff74c2546d0>
```

We also have a look at the data resolution matrix telling us about the inter-dependency and information content of the individual data.

```im = plt.imshow(RD, vmin=-0.3, vmax=0.3, cmap="RdBu_r")
plt.colorbar(im)
```
```<matplotlib.colorbar.Colorbar object at 0x7ff74d369c60>
```

As expected, the neighboring time gates are highly correlated. The most important gates are the lowest and highest one as they define the first and last layers and cannot be replaced by others.

Total running time of the script: (1 minutes 22.766 seconds)

Gallery generated by Sphinx-Gallery