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.86e+01  2.00e+05  0.00e+00  2.00e+05                         0           inf          inf
   1  1.86e+01  3.19e+02  4.23e+01  1.11e+03    2.48e+06      0      13       7.32e-04     1.81e+03
   2  1.86e+01  3.61e+01  4.54e+01  8.82e+02    1.81e+03      0      22       4.55e-04     8.25e-01
   3  1.86e+00  1.19e+01  4.80e+01  1.01e+02    5.09e+02      0      50       1.83e+01     9.34e+03   Skip BFGS
------------------------- STOP! -------------------------
1 : |fc-fOld| = 1.9315e+01 <= tolF*(1+|f0|) = 2.0000e+04
0 : |xc-x_last| = 1.4864e-01 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x|    = 5.0886e+02 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 5.0886e+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.5382521687173576), 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.22e-02  2.00e+05  0.00e+00  2.00e+05                         0           inf          inf
   1  3.22e-02  3.45e+01  3.52e+02  4.59e+01    2.48e+06      0      14       2.94e-04     7.29e+02
geophys. misfits: 34.5 (target 20.0 [False]) | smallness misfit: 2975.6 (target: 100.0 [False])
mref changed in  26  places
Beta cooling evaluation: progress: [34.5]; minimum progress targets: [160000.]
   2  3.22e-02  7.56e+00  5.60e+01  9.37e+00    7.29e+02      0      43       8.14e-04     5.93e-01
geophys. misfits: 7.6 (target 20.0 [True]) | smallness misfit: 2097.4 (target: 100.0 [False])
mref changed in  1  places
Beta cooling evaluation: progress: [7.6]; minimum progress targets: [27.6]
Warming alpha_pgi to favor clustering:  2.643834989913313
   3  3.22e-02  7.74e+00  7.41e+01  1.01e+01    2.62e+00      0      50       3.30e+01     8.65e+01
geophys. misfits: 7.7 (target 20.0 [True]) | smallness misfit: 1222.5 (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: [7.7]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  6.828472989300841
   4  3.22e-02  7.90e+00  1.05e+02  1.13e+01    8.66e+01      0      50       7.05e-01     6.11e+01
geophys. misfits: 7.9 (target 20.0 [True]) | smallness misfit: 676.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: [7.9]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  17.288939607168714
   5  3.22e-02  8.20e+00  1.51e+02  1.31e+01    6.16e+01      0      50       1.94e+00     1.19e+02
geophys. misfits: 8.2 (target 20.0 [True]) | smallness misfit: 449.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: [8.2]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  42.175135379185726
   6  3.22e-02  8.17e+00  2.15e+02  1.51e+01    1.20e+02      0      50       2.67e+00     3.22e+02
geophys. misfits: 8.2 (target 20.0 [True]) | smallness misfit: 267.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: [8.2]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  103.29858992275079
   7  3.22e-02  8.87e+00  3.15e+02  1.90e+01    3.23e+02      0      50       3.68e+00     1.19e+03
geophys. misfits: 8.9 (target 20.0 [True]) | smallness misfit: 189.8 (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: [8.9]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  232.97822307645592
   8  3.22e-02  9.11e+00  5.15e+02  2.57e+01    1.19e+03      0      50       7.98e-01     9.49e+02
geophys. misfits: 9.1 (target 20.0 [True]) | smallness misfit: 165.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.1]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  511.35516163272905
   9  3.22e-02  1.06e+01  9.02e+02  3.96e+01    9.54e+02      0      50       1.90e+00     1.82e+03
geophys. misfits: 10.6 (target 20.0 [True]) | smallness misfit: 150.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: [10.6]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  969.0743313683653
  10  3.22e-02  1.31e+01  1.48e+03  6.08e+01    1.82e+03      0      50       2.66e+00     4.85e+03
geophys. misfits: 13.1 (target 20.0 [True]) | smallness misfit: 138.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: [13.1]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  1477.194388636353
  11  3.22e-02  1.68e+01  2.04e+03  8.26e+01    4.85e+03      0      50       1.07e+00     5.19e+03
geophys. misfits: 16.8 (target 20.0 [True]) | smallness misfit: 129.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: [16.8]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  1760.4240943412674
  12  3.22e-02  1.86e+01  2.34e+03  9.41e+01    5.19e+03      0      50       1.16e+00     6.01e+03   Skip BFGS
geophys. misfits: 18.6 (target 20.0 [True]) | smallness misfit: 125.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: [18.6]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  1889.3671426041958
  13  3.22e-02  2.00e+01  2.46e+03  9.92e+01    6.01e+03      0      50       1.17e+00     7.04e+03   Skip BFGS
geophys. misfits: 20.0 (target 20.0 [True]) | smallness misfit: 122.8 (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: [20.]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  1894.0370339354456
  14  3.22e-02  2.01e+01  2.46e+03  9.93e+01    7.04e+03      0      50       8.89e-01     6.26e+03
geophys. misfits: 20.1 (target 20.0 [False]) | smallness misfit: 122.6 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.1]; minimum progress targets: [22.]
  15  3.22e-02  2.05e+01  2.40e+03  9.81e+01    6.26e+03      0      50       3.22e-02     2.02e+02   Skip BFGS
geophys. misfits: 20.5 (target 20.0 [False]) | smallness misfit: 119.8 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.5]; minimum progress targets: [22.]
  16  3.22e-02  2.10e+01  2.39e+03  9.80e+01    2.02e+02      0      50       8.99e-01     1.81e+02
geophys. misfits: 21.0 (target 20.0 [False]) | smallness misfit: 119.0 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [21.]; minimum progress targets: [22.]
  17  3.22e-02  2.10e+01  2.39e+03  9.80e+01    1.81e+02      0      50       1.02e+00     1.84e+02   Skip BFGS
geophys. misfits: 21.0 (target 20.0 [False]) | smallness misfit: 118.9 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [21.]; minimum progress targets: [22.]
  18  3.22e-02  2.10e+01  2.39e+03  9.80e+01    1.84e+02      0      50       2.72e-01     5.01e+01
geophys. misfits: 21.0 (target 20.0 [False]) | smallness misfit: 118.9 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [21.]; minimum progress targets: [22.]
  19  3.22e-02  2.11e+01  2.39e+03  9.80e+01    5.01e+01      0      50       1.93e-01     9.65e+00   Skip BFGS
geophys. misfits: 21.1 (target 20.0 [False]) | smallness misfit: 118.8 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [21.1]; minimum progress targets: [22.]
  20  3.22e-02  2.11e+01  2.38e+03  9.80e+01    9.65e+00      0      50       1.04e+00     1.01e+01
geophys. misfits: 21.1 (target 20.0 [False]) | smallness misfit: 118.8 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [21.1]; minimum progress targets: [22.]
------------------------- STOP! -------------------------
1 : |fc-fOld| = 5.7172e-05 <= tolF*(1+|f0|) = 2.0000e+04
1 : |xc-x_last| = 1.7692e-04 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x|    = 1.0053e+01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 1.0053e+01 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      20    <= iter          =     20
------------------------- DONE! -------------------------
Total running time of the script: (0 minutes 13.171 seconds)
Estimated memory usage: 319 MB
