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, 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)

Running inversion with SimPEG v0.22.2
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.83e+01  2.00e+05  0.00e+00  2.00e+05    2.42e+06      0
   1  1.83e+01  3.40e+01  4.59e+01  8.76e+02    7.54e-06      0
------------------------- STOP! -------------------------
0 : |fc-fOld| = 1.9912e+05 <= tolF*(1+|f0|) = 2.0000e+04
0 : |xc-x_last| = 4.0098e+00 <= tolX*(1+|x0|) = 1.0000e-01
1 : |proj(x-g)-x|    = 7.5419e-06 <= tolG          = 1.0000e-01
1 : |proj(x-g)-x|    = 7.5419e-06 <= 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()
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.22.2
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.5281823872507709, 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.16e-02  2.00e+05  0.00e+00  2.00e+05    2.42e+06      0
geophys. misfits: 2.0 (target 20.0 [True]) | smallness misfit: 4663.6 (target: 100.0 [False])
mref changed in  26  places
Beta cooling evaluation: progress: [2.]; minimum progress targets: [160000.]
Warming alpha_pgi to favor clustering:  10.23238713518336
   1  3.16e-02  1.95e+00  5.13e+02  1.82e+01    1.92e+01      0
geophys. misfits: 3.2 (target 20.0 [True]) | smallness misfit: 910.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: [3.2]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  63.141625636935046
   2  3.16e-02  3.24e+00  6.55e+02  2.39e+01    4.36e+01      0
geophys. misfits: 4.7 (target 20.0 [True]) | smallness misfit: 478.4 (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: [4.7]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  270.2127421572115
   3  3.16e-02  4.67e+00  1.42e+03  4.96e+01    2.69e+02      0
geophys. misfits: 7.4 (target 20.0 [True]) | smallness misfit: 394.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: [7.4]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  725.5914794315696
   4  3.16e-02  7.45e+00  3.00e+03  1.02e+02    1.45e+03      0
geophys. misfits: 15.2 (target 20.0 [True]) | smallness misfit: 341.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: [15.2]; minimum progress targets: [22.]
Warming alpha_pgi to favor clustering:  953.3968727578674
   5  3.16e-02  1.52e+01  3.40e+03  1.23e+02    3.53e+03      0
geophys. misfits: 20.3 (target 20.0 [False]) | smallness misfit: 321.9 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.3]; minimum progress targets: [22.]
   6  3.16e-02  2.03e+01  3.22e+03  1.22e+02    4.49e+03      0   Skip BFGS
geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.]
   7  3.16e-02  2.02e+01  3.22e+03  1.22e+02    4.43e+03      0
geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.]
   8  3.16e-02  2.02e+01  3.22e+03  1.22e+02    4.44e+03      0
geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.]
   9  3.16e-02  2.02e+01  3.22e+03  1.22e+02    4.43e+03      0
geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.]
  10  3.16e-02  2.02e+01  3.22e+03  1.22e+02    4.43e+03      0
geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.]
  11  3.16e-02  2.02e+01  3.22e+03  1.22e+02    4.43e+03      0
geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.]
  12  3.16e-02  2.02e+01  3.22e+03  1.22e+02    4.43e+03      0
geophys. misfits: 20.2 (target 20.0 [False]) | smallness misfit: 322.0 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.2]; minimum progress targets: [22.]
  13  3.16e-02  2.02e+01  3.22e+03  1.22e+02    4.43e+03      0
geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.]
  14  3.16e-02  2.07e+01  3.18e+03  1.21e+02    7.96e+00      0   Skip BFGS
geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.]
  15  3.16e-02  2.07e+01  3.18e+03  1.21e+02    6.78e+00      0   Skip BFGS
geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.]
  16  3.16e-02  2.07e+01  3.18e+03  1.21e+02    1.18e+00      0
geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.]
  17  3.16e-02  2.07e+01  3.18e+03  1.21e+02    9.20e-01      0   Skip BFGS
geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.]
  18  3.16e-02  2.07e+01  3.18e+03  1.21e+02    6.82e-01      0
geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.]
  19  3.16e-02  2.07e+01  3.18e+03  1.21e+02    5.17e-01      0   Skip BFGS
geophys. misfits: 20.7 (target 20.0 [False]) | smallness misfit: 318.5 (target: 100.0 [False])
mref changed in  0  places
Beta cooling evaluation: progress: [20.7]; minimum progress targets: [22.]
  20  3.16e-02  2.07e+01  3.18e+03  1.21e+02    2.80e-01      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 3.6151e-08 <= tolF*(1+|f0|) = 2.0000e+04
1 : |xc-x_last| = 7.4356e-06 <= tolX*(1+|x0|) = 1.0000e-01
0 : |proj(x-g)-x|    = 2.8026e-01 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 2.8026e-01 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =      20    <= iter          =     20
------------------------- DONE! -------------------------

Estimated memory usage: 166 MB