Note
Go to the end to download the full example code.
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
survey = prob.make_synthetic_data(mtrue, relative_error=std, add_noise=True)
# 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, cg_maxiter=50, cg_rtol=1e-3)
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)
Running inversion with SimPEG v0.25.0
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
0 1.88e+01 2.00e+05 0.00e+00 2.00e+05 0 inf inf
1 1.88e+01 3.06e+02 4.18e+01 1.09e+03 2.52e+06 0 13 6.85e-04 1.72e+03
2 1.88e+01 3.21e+01 4.48e+01 8.74e+02 1.72e+03 0 22 1.58e-04 2.73e-01 Skip BFGS
3 1.88e+00 3.97e+00 4.78e+01 9.39e+01 5.11e+02 0 32 2.57e-04 1.31e-01 Skip BFGS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 2.2448e+01 <= tolF*(1+|f0|) = 2.0000e+04
0 : |xc-x_last| = 1.8099e-01 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x| = 5.1111e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 5.1111e+02 <= 1e3*eps = 1.0000e-02
0 : maxIter = 10 <= iter = 3
------------------------- 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, cg_maxiter=50, cg_rtol=1e-3)
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()
addmref = directives.PGI_AddMrefInSmooth(verbose=True)
# 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[0]):
axes[0].plot(prob.G[i, :])
axes[0].set_title("Columns of matrix G")
axes[1].hist(mtrue, bins=20, linewidth=3.0, density=True, color="k")
axes[1].set_xlabel("Model value")
axes[1].set_xlabel("Occurence")
axes[1].hist(mnormal, bins=20, density=True, color="b")
axes[1].hist(mcluster, bins=20, density=True, color="r")
axes[1].legend(["Mtrue Hist.", "L2 Model Hist.", "PGI Model Hist."])
axes[2].plot(mesh.cell_centers_x, mtrue, color="black", linewidth=3)
axes[2].plot(mesh.cell_centers_x, mnormal, color="blue")
axes[2].plot(mesh.cell_centers_x, mcluster, "r-")
axes[2].plot(mesh.cell_centers_x, invProb.reg.objfcts[0].reference_model, "r--")
axes[2].legend(("True Model", "L2 Model", "PGI Model", "Learned Mref"))
axes[2].set_ylim([-2, 2])
plt.show()
Running inversion with SimPEG v0.25.0
Alpha scales: [np.float64(0.5272875999563319), np.float64(0.0)]
<class 'simpeg.regularization.pgi.PGIsmallness'>
================================================= Projected GNCG =================================================
# beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment
-----------------------------------------------------------------------------------------------------------------
0 3.28e-02 2.00e+05 0.00e+00 2.00e+05 0 inf inf
1 3.28e-02 3.09e+01 3.47e+02 4.23e+01 2.52e+06 0 14 4.47e-04 1.13e+03
geophys. misfits: 30.9 (target 20.0 [False]) | smallness misfit: 2957.6 (target: 100.0 [False])
mref changed in 25 places
Beta cooling evaluation: progress: [30.9]; minimum progress targets: [160000.]
2 3.28e-02 5.11e-01 7.88e+01 3.10e+00 1.12e+03 0 29 9.63e-04 1.08e+00
geophys. misfits: 0.5 (target 20.0 [True]) | smallness misfit: 3944.4 (target: 100.0 [False])
mref changed in 3 places
Beta cooling evaluation: progress: [0.5]; minimum progress targets: [24.7]
Warming alpha_pgi to favor clustering: 39.12906831118989
3 3.28e-02 2.29e+00 3.15e+02 1.26e+01 7.13e+01 0 50 2.31e+00 1.65e+02
geophys. misfits: 2.3 (target 20.0 [True]) | smallness misfit: 540.1 (target: 100.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: [2.3]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering: 341.5125412107303
4 3.28e-02 5.34e+00 1.31e+03 4.84e+01 2.66e+02 0 50 5.97e+00 1.59e+03
geophys. misfits: 5.3 (target 20.0 [True]) | smallness misfit: 343.0 (target: 100.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: [5.3]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering: 1280.1350865738573
5 3.28e-02 2.40e+01 3.54e+03 1.40e+02 1.66e+03 0 50 2.94e+00 4.89e+03
geophys. misfits: 24.0 (target 20.0 [False]) | smallness misfit: 265.2 (target: 100.0 [False])
mref changed in 0 places
Beta cooling evaluation: progress: [24.]; minimum progress targets: [22.]
Decreasing beta to counter data misfit increase.
6 1.64e-02 1.03e+01 4.11e+03 7.77e+01 4.90e+03 0 50 5.70e-01 2.79e+03
geophys. misfits: 10.3 (target 20.0 [True]) | smallness misfit: 309.7 (target: 100.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: [10.3]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering: 2491.805533412284
7 1.64e-02 2.30e+01 6.82e+03 1.35e+02 2.81e+03 0 50 1.70e+00 4.76e+03
geophys. misfits: 23.0 (target 20.0 [False]) | smallness misfit: 267.7 (target: 100.0 [False])
mref changed in 0 places
Beta cooling evaluation: progress: [23.]; minimum progress targets: [22.]
Decreasing beta to counter data misfit increase.
8 8.21e-03 1.06e+01 7.85e+03 7.51e+01 4.78e+03 0 50 6.37e-01 3.04e+03
geophys. misfits: 10.6 (target 20.0 [True]) | smallness misfit: 309.2 (target: 100.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: [10.6]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering: 4686.890701947216
9 8.21e-03 1.97e+01 1.31e+04 1.27e+02 3.05e+03 0 50 1.54e+00 4.68e+03
geophys. misfits: 19.7 (target 20.0 [True]) | smallness misfit: 276.7 (target: 100.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: [19.7]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering: 4759.385038141183
10 8.21e-03 2.05e+01 1.32e+04 1.29e+02 4.68e+03 0 50 9.84e-01 4.61e+03
geophys. misfits: 20.5 (target 20.0 [False]) | smallness misfit: 273.9 (target: 100.0 [False])
mref changed in 0 places
Beta cooling evaluation: progress: [20.5]; minimum progress targets: [22.]
11 8.21e-03 2.16e+01 1.30e+04 1.29e+02 4.61e+03 0 50 9.78e-01 4.50e+03
geophys. misfits: 21.6 (target 20.0 [False]) | smallness misfit: 271.0 (target: 100.0 [False])
mref changed in 0 places
Beta cooling evaluation: progress: [21.6]; minimum progress targets: [22.]
12 8.21e-03 2.18e+01 1.30e+04 1.29e+02 4.50e+03 0 50 1.01e+00 4.53e+03
geophys. misfits: 21.8 (target 20.0 [False]) | smallness misfit: 270.4 (target: 100.0 [False])
mref changed in 0 places
Beta cooling evaluation: progress: [21.8]; minimum progress targets: [22.]
13 8.21e-03 2.21e+01 1.29e+04 1.28e+02 4.53e+03 0 50 2.81e-03 1.27e+01 Skip BFGS
geophys. misfits: 22.1 (target 20.0 [False]) | smallness misfit: 268.2 (target: 100.0 [False])
mref changed in 0 places
Beta cooling evaluation: progress: [22.1]; minimum progress targets: [22.]
Decreasing beta to counter data misfit increase.
14 4.10e-03 9.81e+00 1.50e+04 7.13e+01 2.84e+02 0 50 2.15e+00 6.12e+02
geophys. misfits: 9.8 (target 20.0 [True]) | smallness misfit: 311.6 (target: 100.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: [9.8]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering: 9705.776963575992
15 4.10e-03 2.27e+01 2.60e+04 1.30e+02 6.81e+02 0 50 1.40e-01 9.54e+01
geophys. misfits: 22.7 (target 20.0 [False]) | smallness misfit: 266.6 (target: 100.0 [False])
mref changed in 0 places
Beta cooling evaluation: progress: [22.7]; minimum progress targets: [22.]
Decreasing beta to counter data misfit increase.
16 2.05e-03 9.70e+00 3.04e+04 7.20e+01 3.06e+02 0 50 7.23e-03 2.21e+00
geophys. misfits: 9.7 (target 20.0 [True]) | smallness misfit: 311.3 (target: 100.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: [9.7]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering: 20003.699403770454
17 2.05e-03 2.36e+01 5.30e+04 1.32e+02 3.34e+02 0 50 7.15e-02 2.38e+01
geophys. misfits: 23.6 (target 20.0 [False]) | smallness misfit: 264.4 (target: 100.0 [False])
mref changed in 0 places
Beta cooling evaluation: progress: [23.6]; minimum progress targets: [22.]
Decreasing beta to counter data misfit increase.
18 1.03e-03 1.00e+01 6.21e+04 7.37e+01 2.98e+02 0 50 9.99e-03 2.98e+00
geophys. misfits: 10.0 (target 20.0 [True]) | smallness misfit: 309.6 (target: 100.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: [10.]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering: 39819.73501690069
19 1.03e-03 2.35e+01 1.06e+05 1.32e+02 3.20e+02 0 50 1.63e-02 5.22e+00
geophys. misfits: 23.5 (target 20.0 [False]) | smallness misfit: 264.8 (target: 100.0 [False])
mref changed in 0 places
Beta cooling evaluation: progress: [23.5]; minimum progress targets: [22.]
Decreasing beta to counter data misfit increase.
20 5.13e-04 9.99e+00 1.24e+05 7.34e+01 2.97e+02 0 48 8.30e-04 2.47e-01
geophys. misfits: 10.0 (target 20.0 [True]) | smallness misfit: 309.9 (target: 100.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: [10.]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering: 79701.88625681559
------------------------- STOP! -------------------------
1 : |fc-fOld| = 5.9135e+01 <= tolF*(1+|f0|) = 2.0000e+04
1 : |xc-x_last| = 4.1282e-02 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x| = 3.2211e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 3.2211e+02 <= 1e3*eps = 1.0000e-02
1 : maxIter = 20 <= iter = 20
------------------------- DONE! -------------------------
Total running time of the script: (0 minutes 15.167 seconds)
Estimated memory usage: 322 MB