# Petrophysically guided inversion (PGI): Linear example#

We do a comparison between the classic least-squares inversion and our formulation of a petrophysically constrained inversion. We explore it through the UBC linear example.

## Tikhonov Inversion##

```import discretize as Mesh
import matplotlib.pyplot as plt
import numpy as np
from SimPEG import (
data_misfit,
directives,
inverse_problem,
inversion,
maps,
optimization,
regularization,
simulation,
utils,
)

# Random seed for reproductibility
np.random.seed(1)
# Mesh
N = 100
mesh = Mesh.TensorMesh([N])

# Survey design parameters
nk = 20
jk = np.linspace(1.0, 60.0, nk)
p = -0.25
q = 0.25

# Physics
def g(k):
return np.exp(p * jk[k] * mesh.cell_centers_x) * np.cos(
np.pi * q * jk[k] * mesh.cell_centers_x
)

G = np.empty((nk, mesh.nC))

for i in range(nk):
G[i, :] = g(i)

# True model
mtrue = np.zeros(mesh.nC)
mtrue[mesh.cell_centers_x > 0.2] = 1.0
mtrue[mesh.cell_centers_x > 0.35] = 0.0
t = (mesh.cell_centers_x - 0.65) / 0.25
indx = np.abs(t) < 1
mtrue[indx] = -(((1 - t ** 2.0) ** 2.0)[indx])

mtrue = np.zeros(mesh.nC)
mtrue[mesh.cell_centers_x > 0.3] = 1.0
mtrue[mesh.cell_centers_x > 0.45] = -0.5
mtrue[mesh.cell_centers_x > 0.6] = 0

# SimPEG problem and survey
prob = simulation.LinearSimulation(mesh, G=G, model_map=maps.IdentityMap())
std = 0.01

# Setup the inverse problem
reg = regularization.WeightedLeastSquares(mesh, alpha_s=1.0, alpha_x=1.0)
dmis = data_misfit.L2DataMisfit(data=survey, simulation=prob)
opt = optimization.ProjectedGNCG(maxIter=10, maxIterCG=50, tolCG=1e-4)
invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)
directiveslist = [
directives.BetaEstimate_ByEig(beta0_ratio=1e-5),
directives.BetaSchedule(coolingFactor=10.0, coolingRate=2),
directives.TargetMisfit(),
]

inv = inversion.BaseInversion(invProb, directiveList=directiveslist)
m0 = np.zeros_like(mtrue)

mnormal = inv.run(m0)
```
```SimPEG.InvProblem will set Regularization.reference_model to m0.
SimPEG.InvProblem will set Regularization.reference_model to m0.
SimPEG.InvProblem will set Regularization.reference_model to m0.

SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using the default solver Pardiso and no solver_opts.***

model has any nan: 0
=============================== Projected GNCG ===============================
#     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
0  1.86e+01  1.00e+05  0.00e+00  1.00e+05    1.26e+06      0
1  1.86e+01  2.17e+01  2.28e+01  4.46e+02    2.18e-05      0
------------------------- STOP! -------------------------
0 : |fc-fOld| = 9.9554e+04 <= tolF*(1+|f0|) = 1.0000e+04
0 : |xc-x_last| = 4.0057e+00 <= tolX*(1+|x0|) = 1.0000e-01
1 : |proj(x-g)-x|    = 2.1819e-05 <= tolG          = 1.0000e-01
1 : |proj(x-g)-x|    = 2.1819e-05 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      10    <= iter          =      1
------------------------- DONE! -------------------------
```

## Petrophysically constrained inversion ##

```# fit a Gaussian Mixture Model with n components
# on the true model to simulate the laboratory
# petrophysical measurements
n = 3
clf = utils.WeightedGaussianMixture(
mesh=mesh,
n_components=n,
covariance_type="full",
max_iter=100,
n_init=3,
reg_covar=5e-4,
)
clf.fit(mtrue.reshape(-1, 1))

# Petrophyically constrained regularization
reg = regularization.PGI(
gmmref=clf,
mesh=mesh,
alpha_pgi=1.0,
alpha_x=1.0,
)

# Optimization
opt = optimization.ProjectedGNCG(maxIter=20, maxIterCG=50, tolCG=1e-4)
opt.remember("xc")

# Setup new inverse problem
invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)

# directives
Alphas = directives.AlphasSmoothEstimate_ByEig(alpha0_ratio=10.0, verbose=True)
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-8)
betaIt = directives.PGI_BetaAlphaSchedule(
verbose=True,
coolingFactor=2.0,
warmingFactor=1.0,
tolerance=0.1,
update_rate=1,
progress=0.2,
)
targets = directives.MultiTargetMisfits(verbose=True)
petrodir = directives.PGI_UpdateParameters()

# Setup Inversion
inv = inversion.BaseInversion(
invProb, directiveList=[Alphas, beta, petrodir, targets, addmref, betaIt]
)

# Initial model same as for WeightedLeastSquares
mcluster = inv.run(m0)

# Final Plot
fig, axes = plt.subplots(1, 3, figsize=(12 * 1.2, 4 * 1.2))
for i in range(prob.G.shape):
axes.plot(prob.G[i, :])
axes.set_title("Columns of matrix G")

axes.hist(mtrue, bins=20, linewidth=3.0, density=True, color="k")
axes.set_xlabel("Model value")
axes.set_xlabel("Occurence")
axes.hist(mnormal, bins=20, density=True, color="b")
axes.hist(mcluster, bins=20, density=True, color="r")
axes.legend(["Mtrue Hist.", "L2 Model Hist.", "PGI Model Hist."])

axes.plot(mesh.cell_centers_x, mtrue, color="black", linewidth=3)
axes.plot(mesh.cell_centers_x, mnormal, color="blue")
axes.plot(mesh.cell_centers_x, mcluster, "r-")
axes.plot(mesh.cell_centers_x, invProb.reg.objfcts.reference_model, "r--")

axes.legend(("True Model", "L2 Model", "PGI Model", "Learned Mref"))
axes.set_ylim([-2, 2])

plt.show()
``` ```SimPEG.InvProblem will set Regularization.reference_model to m0.
SimPEG.InvProblem will set Regularization.reference_model to m0.

SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using the default solver Pardiso and no solver_opts.***

Alpha scales: [0.5249361360067218, 0.0]
<class 'SimPEG.regularization.pgi.PGIsmallness'>
model has any nan: 0
=============================== Projected GNCG ===============================
#     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
0  3.29e-02  1.00e+05  0.00e+00  1.00e+05    1.26e+06      0
geophys. misfits: 5.9 (target 10.0 [True]) | smallness misfit: 1922.4 (target: 50.0 [False])
mref changed in  26  places
Beta cooling evaluation: progress: [5.9] ; minimum progress targets: [80000.]
Warming alpha_pgi to favor clustering:  1.6860956791065138
1  3.29e-02  5.93e+00  4.68e+01  7.47e+00    2.89e+00      0
geophys. misfits: 6.0 (target 10.0 [True]) | smallness misfit: 913.5 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [6.] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  2.8177119557212387
2  3.29e-02  5.98e+00  4.47e+01  7.45e+00    8.10e-01      0   Skip BFGS
geophys. misfits: 6.1 (target 10.0 [True]) | smallness misfit: 663.3 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [6.1] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  4.646567486373403
3  3.29e-02  6.06e+00  5.28e+01  7.80e+00    1.12e+00      0   Skip BFGS
geophys. misfits: 6.2 (target 10.0 [True]) | smallness misfit: 462.7 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [6.2] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  7.542826304193983
4  3.29e-02  6.16e+00  6.12e+01  8.17e+00    1.31e+00      0
geophys. misfits: 6.3 (target 10.0 [True]) | smallness misfit: 313.5 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [6.3] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  12.032866446632116
5  3.29e-02  6.27e+00  6.95e+01  8.56e+00    8.89e+00      0
geophys. misfits: 6.4 (target 10.0 [True]) | smallness misfit: 212.0 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [6.4] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  18.85630634342623
6  3.29e-02  6.38e+00  7.80e+01  8.95e+00    3.90e+00      0
geophys. misfits: 6.5 (target 10.0 [True]) | smallness misfit: 146.9 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [6.5] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  29.016267549390637
7  3.29e-02  6.50e+00  8.68e+01  9.35e+00    2.56e+00      0
geophys. misfits: 6.6 (target 10.0 [True]) | smallness misfit: 106.8 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [6.6] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  43.7896331598517
8  3.29e-02  6.63e+00  9.63e+01  9.79e+00    3.14e+00      0
geophys. misfits: 6.8 (target 10.0 [True]) | smallness misfit: 82.2 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [6.8] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  64.70122388136394
9  3.29e-02  6.77e+00  1.07e+02  1.03e+01    3.96e+00      0
geophys. misfits: 6.9 (target 10.0 [True]) | smallness misfit: 66.5 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [6.9] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  93.35424462712366
10  3.29e-02  6.93e+00  1.19e+02  1.09e+01    4.86e+00      0
geophys. misfits: 7.1 (target 10.0 [True]) | smallness misfit: 56.0 (target: 50.0 [False])
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [7.1] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  131.15514751222793
11  3.29e-02  7.12e+00  1.33e+02  1.15e+01    5.89e+00      0
geophys. misfits: 7.3 (target 10.0 [True]) | smallness misfit: 48.6 (target: 50.0 [True])
All targets have been reached
mref changed in  0  places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
Beta cooling evaluation: progress: [7.3] ; minimum progress targets: [11.]
Warming alpha_pgi to favor clustering:  178.91791743902695
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.0000e+04
1 : |xc-x_last| = 2.6698e-02 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x|    = 5.8930e+00 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 5.8930e+00 <= 1e3*eps       = 1.0000e-02
0 : maxIter   =      20    <= iter          =     12
------------------------- DONE! -------------------------
```

Total running time of the script: ( 0 minutes 6.748 seconds)

Estimated memory usage: 19 MB

Gallery generated by Sphinx-Gallery