Joint PGI of Gravity + Magnetic on an Octree mesh without petrophysical information
%matplotlib inlineThis tutorial shows through a joint inversion of Gravity and Magnetic data on an Octree mesh how to use the PGI framework introduced in Astic & Oldenburg (2019) and Astic et al. (2021) to make geologic assumptions and learn a suitable petrophysical distribution when no quantitative petrophysical information is available.
Thibaut Astic, Douglas W. Oldenburg,
A framework for petrophysically and geologically guided geophysical inversion
using a dynamic Gaussian mixture model prior, Geophysical Journal International,
Volume 219, Issue 3, December 2019, Pages 1989–2012, DOI:
10.1093/gji/ggz389 <https://doi.org/10.1093/gji/ggz389>_.
Thibaut Astic, Lindsey J. Heagy, Douglas W Oldenburg,
Petrophysically and geologically guided multi-physics inversion using a dynamic
Gaussian mixture model, Geophysical Journal International,
Volume 224, Issue 1, January 2021, Pages 40-68, DOI: 10.1093/gji/ggaa378 <https://doi.org/10.1093/gji/ggaa378>_.
Import modules¶
from discretize import TreeMesh
from discretize.utils import active_from_xyz
import matplotlib.pyplot as plt
import numpy as np
import simpeg.potential_fields as pf
from simpeg import (
data_misfit,
directives,
inverse_problem,
inversion,
maps,
optimization,
regularization,
utils,
)
from simpeg.utils import io_utilsSetup¶
# Load Mesh
mesh_file = io_utils.download(
"https://storage.googleapis.com/simpeg/pgi_tutorial_assets/mesh_tutorial.ubc"
)
mesh = TreeMesh.read_UBC(mesh_file)
# Load True geological model for comparison with inversion result
true_geology_file = io_utils.download(
"https://storage.googleapis.com/simpeg/pgi_tutorial_assets/geology_true.mod"
)
true_geology = mesh.read_model_UBC(true_geology_file)
# Plot true geology model
fig, ax = plt.subplots(1, 4, figsize=(20, 4))
ticksize, labelsize = 14, 16
for _, axx in enumerate(ax):
axx.set_aspect(1)
axx.tick_params(labelsize=ticksize)
mesh.plot_slice(
true_geology,
normal="X",
ax=ax[0],
ind=-17,
clim=[0, 2],
pcolor_opts={"cmap": "inferno_r"},
grid=True,
)
mesh.plot_slice(
true_geology,
normal="Y",
ax=ax[1],
clim=[0, 2],
pcolor_opts={"cmap": "inferno_r"},
grid=True,
)
geoplot = mesh.plot_slice(
true_geology,
normal="Z",
ax=ax[2],
clim=[0, 2],
ind=-10,
pcolor_opts={"cmap": "inferno_r"},
grid=True,
)
geocb = plt.colorbar(geoplot[0], cax=ax[3], ticks=[0, 1, 2])
geocb.set_label(
"True geology model\n(classification/density/mag. susc.)", fontsize=labelsize
)
geocb.set_ticklabels(
["BCKGRD (0 g/cc; 0 SI)", "PK (-0.8 g/cc; 5e-3 SI)", "VK (-0.2 g/cc; 2e-2 SI)"]
)
geocb.ax.tick_params(labelsize=ticksize)
ax[3].set_aspect(10)
plt.show()
# Load geophysical data
data_grav_file = io_utils.download(
"https://storage.googleapis.com/simpeg/pgi_tutorial_assets/gravity_data.obs"
)
data_grav = io_utils.read_grav3d_ubc(data_grav_file)
data_mag_file = io_utils.download(
"https://storage.googleapis.com/simpeg/pgi_tutorial_assets/magnetic_data.obs"
)
data_mag = io_utils.read_mag3d_ubc(data_mag_file)
# plot data and mesh
fig, ax = plt.subplots(2, 2, figsize=(15, 10))
ax = ax.reshape(-1)
plt.gca().set_aspect("equal")
plt.gca().set_xlim(
[
data_mag.survey.receiver_locations[:, 0].min(),
data_mag.survey.receiver_locations[:, 0].max(),
],
)
plt.gca().set_ylim(
[
data_mag.survey.receiver_locations[:, 1].min(),
data_mag.survey.receiver_locations[:, 1].max(),
]
)
mesh.plot_slice(
np.ones(mesh.nC),
normal="Z",
ind=int(-10),
grid=True,
pcolor_opts={"cmap": "Greys"},
ax=ax[0],
)
mm = utils.plot2Ddata(
data_grav.survey.receiver_locations,
-data_grav.dobs,
ax=ax[0],
level=True,
nx=20,
ny=20,
dataloc=True,
ncontour=12,
shade=True,
contourOpts={"cmap": "Blues_r", "alpha": 0.8},
levelOpts={"colors": "k", "linewidths": 0.5, "linestyles": "dashed"},
)
ax[0].set_aspect(1)
ax[0].set_title(
"Gravity data values and locations,\nwith mesh and geology overlays", fontsize=16
)
plt.colorbar(mm[0], cax=ax[2], orientation="horizontal")
ax[2].set_aspect(0.05)
ax[2].set_title("mGal", fontsize=16)
mesh.plot_slice(
np.ones(mesh.nC),
normal="Z",
ind=int(-10),
grid=True,
pcolor_opts={"cmap": "Greys"},
ax=ax[1],
)
mm = utils.plot2Ddata(
data_mag.survey.receiver_locations,
data_mag.dobs,
ax=ax[1],
level=True,
nx=20,
ny=20,
dataloc=True,
ncontour=11,
shade=True,
contourOpts={"cmap": "Reds", "alpha": 0.8},
levelOpts={"colors": "k", "linewidths": 0.5, "linestyles": "dashed"},
)
ax[1].set_aspect(1)
ax[1].set_title(
"Magnetic data values and locations,\nwith mesh and geology overlays", fontsize=16
)
plt.colorbar(mm[0], cax=ax[3], orientation="horizontal")
ax[3].set_aspect(0.05)
ax[3].set_title("nT", fontsize=16)
# overlay true geology model for comparison
indz = -9
indslicezplot = mesh.gridCC[:, 2] == mesh.cell_centers_z[indz]
for i in range(2):
utils.plot2Ddata(
mesh.gridCC[indslicezplot][:, [0, 1]],
true_geology[indslicezplot],
nx=200,
ny=200,
contourOpts={"alpha": 0},
clim=[0, 2],
ax=ax[i],
level=True,
ncontour=2,
levelOpts={"colors": "k", "linewidths": 2, "linestyles": "--"},
method="nearest",
)
plt.subplots_adjust(hspace=-0.25, wspace=0.1)
plt.show()
# Load Topo
topo_file = io_utils.download(
"https://storage.googleapis.com/simpeg/pgi_tutorial_assets/CDED_Lake_warp.xyz"
)
topo = np.genfromtxt(topo_file, skip_header=1)
# find the active cells
actv = active_from_xyz(mesh, topo, "CC")
# Create active map to go from reduce set to full
ndv = np.nan
actvMap = maps.InjectActiveCells(mesh, actv, ndv)
nactv = int(actv.sum())
# Create simulations and data misfits
# Wires mapping
wires = maps.Wires(("den", actvMap.nP), ("sus", actvMap.nP))
gravmap = actvMap * wires.den
magmap = actvMap * wires.sus
idenMap = maps.IdentityMap(nP=nactv)
# Grav problem
simulation_grav = pf.gravity.simulation.Simulation3DIntegral(
survey=data_grav.survey,
mesh=mesh,
rhoMap=wires.den,
active_cells=actv,
engine="choclo",
)
dmis_grav = data_misfit.L2DataMisfit(data=data_grav, simulation=simulation_grav)
# Mag problem
simulation_mag = pf.magnetics.simulation.Simulation3DIntegral(
survey=data_mag.survey,
mesh=mesh,
chiMap=wires.sus,
active_cells=actv,
engine="choclo",
)
dmis_mag = data_misfit.L2DataMisfit(data=data_mag, simulation=simulation_mag)file already exists, new file is called /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/mesh_tutorial.ubc
Downloading https://storage.googleapis.com/simpeg/pgi_tutorial_assets/mesh_tutorial.ubc
saved to: /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/mesh_tutorial.ubc
Download completed!
file already exists, new file is called /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/geology_true.mod
Downloading https://storage.googleapis.com/simpeg/pgi_tutorial_assets/geology_true.mod
saved to: /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/geology_true.mod
Download completed!
file already exists, new file is called /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/gravity_data.obs
Downloading https://storage.googleapis.com/simpeg/pgi_tutorial_assets/gravity_data.obs
saved to: /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/gravity_data.obs
Download completed!
file already exists, new file is called /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/magnetic_data.obs
Downloading https://storage.googleapis.com/simpeg/pgi_tutorial_assets/magnetic_data.obs
saved to: /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/magnetic_data.obs
Download completed!
file already exists, new file is called /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/CDED_Lake_warp.xyz
Downloading https://storage.googleapis.com/simpeg/pgi_tutorial_assets/CDED_Lake_warp.xyz
saved to: /home/ssoler/git/user-tutorials/notebooks/12-pgi-inversion/CDED_Lake_warp.xyz
Download completed!
Create a joint Data Misfit
# Joint data misfit
dmis = 0.5 * dmis_grav + 0.5 * dmis_mag
# initial model
m0 = np.r_[-1e-4 * np.ones(actvMap.nP), 1e-4 * np.ones(actvMap.nP)]Inversion with no petrophysical information about the means¶
In this scenario, we do not know the true petrophysical signature of each rock unit. We thus make geologic assumptions to design a coupling term and perform a multi-physics inversion. in addition to a neutral background, we assume that one rock unit is only less dense, and the third one is only magnetic. As we do not know their mean petrophysical values. We start with an initial guess (-1 g/cc) for the updatable mean density-contrast value of the less dense unit (with a fixed susceptibility of 0 SI). The magnetic-contrasting unit’s updatable susceptibility is initialized at a value of 0.1 SI (with a fixed 0 g/cc density contrast). We then let the algorithm learn a suitable set of means under the set constrained (fixed or updatable value), through the kappa argument, denoting our confidences in each initial mean value (high confidence: fixed value; low confidence: updatable value).
Create a petrophysical GMM initial guess¶
The GMM is our representation of the petrophysical and geological information. Here, we focus on the petrophysical aspect, with the means and covariances of the physical properties of each rock unit. To generate the data above, the PK unit was populated with a density contrast of -0.8 g/cc and a magnetic susceptibility of 0.005 SI. The properties of the HK unit were set at -0.2 g/cc and 0.02 SI. But here, we assume we do not have this information. Thus, we start with initial guess for the means and confidences kappa such that one unit is only less dense and one unit is only magnetic, both embedded in a neutral background. The covariances matrices are set so that we assume petrophysical noise levels of around 0.05 g/cc and 0.001 SI for both unit. The background unit is set at a fixed null contrasts (0 g/cc 0 SI) with a petrophysical noise level of half of the above.
gmmref = utils.WeightedGaussianMixture(
n_components=3, # number of rock units: bckgrd, PK, HK
mesh=mesh, # inversion mesh
actv=actv, # actv cells
covariance_type="diag", # diagonal covariances
)
# required: initialization with fit
# fake random samples, size of the mesh
# number of physical properties: 2 (density and mag.susc)
rng = np.random.default_rng(seed=518936)
gmmref.fit(rng.normal(size=(nactv, 2)))
# set parameters manually
# set phys. prop means for each unit
gmmref.means_ = np.c_[
[0.0, 0.0], # BCKGRD density contrast and mag. susc
[-1, 0.0], # PK
[0, 0.1], # HK
].T
# set phys. prop covariances for each unit
gmmref.covariances_ = np.array(
[[6e-04, 3.175e-07], [2.4e-03, 1.5e-06], [2.4e-03, 1.5e-06]]
)
# important after setting cov. manually: compute precision matrices and cholesky
gmmref.compute_clusters_precisions()
# set global proportions; low-impact as long as not 0 or 1 (total=1)
gmmref.weights_ = np.r_[0.9, 0.075, 0.025]
# Plot the 2D GMM
ax = gmmref.plot_pdf(flag2d=True, plotting_precision=250)
ax[0].set_xlabel("Density contrast [g/cc]")
ax[0].set_ylim([0, 5])
ax[2].set_ylabel("magnetic Susceptibility [SI]")
ax[2].set_xlim([0, 100])
plt.show()
Inverse problem with no mean information¶
# Create PGI regularization
# Sensitivity weighting
wr_grav = np.sum(simulation_grav.G**2.0, axis=0) ** 0.5 / (mesh.cell_volumes[actv])
wr_grav = wr_grav / np.max(wr_grav)
wr_mag = np.sum(simulation_mag.G**2.0, axis=0) ** 0.5 / (mesh.cell_volumes[actv])
wr_mag = wr_mag / np.max(wr_mag)
# create joint PGI regularization with smoothness
reg = regularization.PGI(
gmmref=gmmref,
mesh=mesh,
wiresmap=wires,
maplist=[idenMap, idenMap],
active_cells=actv,
alpha_pgi=1.0,
alpha_x=1.0,
alpha_y=1.0,
alpha_z=1.0,
alpha_xx=0.0,
alpha_yy=0.0,
alpha_zz=0.0,
# use the classification of the initial model (here, all background unit)
# as initial reference model
reference_model=utils.mkvc(
gmmref.means_[gmmref.predict(m0.reshape(actvMap.nP, -1))]
),
weights_list=[wr_grav, wr_mag], # weights each phys. prop. by correct sensW
)
# Directives
# Add directives to the inversion
# ratio to use for each phys prop. smoothness in each direction:
# roughly the ratio of range of each phys. prop.
alpha0_ratio = np.r_[
1e-2 * np.ones(len(reg.objfcts[1].objfcts[1:])),
1e-2 * 100.0 * np.ones(len(reg.objfcts[2].objfcts[1:])),
]
Alphas = directives.AlphasSmoothEstimate_ByEig(alpha0_ratio=alpha0_ratio, verbose=True)
# initialize beta and beta/alpha_s schedule
beta = directives.BetaEstimate_ByEig(beta0_ratio=1e-4)
betaIt = directives.PGI_BetaAlphaSchedule(
verbose=True,
coolingFactor=2.0,
tolerance=0.2,
progress=0.2,
)
# geophy. and petro. target misfits
targets = directives.MultiTargetMisfits(
verbose=True,
chiSmall=0.5, # ask for twice as much clustering (target value is /2)
)
# add learned mref in smooth once stable
MrefInSmooth = directives.PGI_AddMrefInSmooth(
wait_till_stable=True,
verbose=True,
)
# update the parameters in smallness (L2-approx of PGI)
update_smallness = directives.PGI_UpdateParameters(
update_gmm=True, # update the GMM each iteration
kappa=np.c_[ # confidences in each mean phys. prop. of each cluster
1e10
* np.ones(
2
), # fixed background at 0 density, 0 mag. susc. (high confidences of 1e10)
[
0,
1e10,
], # density-contrasting cluster: updatable density mean, fixed mag. susc.
[
1e10,
0,
], # magnetic-contrasting cluster: fixed density mean, updatable mag. susc.
].T,
)
# pre-conditioner
update_Jacobi = directives.UpdatePreconditioner()
# iteratively balance the scaling of the data misfits
scaling_init = directives.ScalingMultipleDataMisfits_ByEig(chi0_ratio=[1.0, 100.0])
scale_schedule = directives.JointScalingSchedule(verbose=True)
# Create inverse problem
# Optimization
# set lower and upper bounds
lowerbound = np.r_[-2.0 * np.ones(actvMap.nP), 0.0 * np.ones(actvMap.nP)]
upperbound = np.r_[0.0 * np.ones(actvMap.nP), 1e-1 * np.ones(actvMap.nP)]
opt = optimization.ProjectedGNCG(
maxIter=30,
lower=lowerbound,
upper=upperbound,
maxIterLS=20,
cg_maxiter=100,
cg_rtol=1e-4,
)
# create inverse problem
invProb = inverse_problem.BaseInvProblem(dmis, reg, opt)
inv = inversion.BaseInversion(
invProb,
# directives: evaluate alphas (and data misfits scales) before beta
directiveList=[
Alphas,
scaling_init,
beta,
update_smallness,
targets,
scale_schedule,
betaIt,
MrefInSmooth,
update_Jacobi,
],
)
# Invert
pgi_model_no_info = inv.run(m0)
# Plot the result with full petrophysical information
density_model_no_info = gravmap * pgi_model_no_info
magsus_model_no_info = magmap * pgi_model_no_info
learned_gmm = reg.objfcts[0].gmm
quasi_geology_model_no_info = actvMap * reg.objfcts[0].compute_quasi_geology_model()
fig, ax = plt.subplots(3, 4, figsize=(15, 10))
for _, axx in enumerate(ax):
for _, axxx in enumerate(axx):
axxx.set_aspect(1)
axxx.tick_params(labelsize=ticksize)
indx = 15
indy = 17
indz = -9
# geology model
mesh.plot_slice(
quasi_geology_model_no_info,
normal="X",
ax=ax[0, 0],
clim=[0, 2],
ind=indx,
pcolor_opts={"cmap": "inferno_r"},
)
mesh.plot_slice(
quasi_geology_model_no_info,
normal="Y",
ax=ax[0, 1],
clim=[0, 2],
ind=indy,
pcolor_opts={"cmap": "inferno_r"},
)
geoplot = mesh.plot_slice(
quasi_geology_model_no_info,
normal="Z",
ax=ax[0, 2],
clim=[0, 2],
ind=indz,
pcolor_opts={"cmap": "inferno_r"},
)
geocb = plt.colorbar(geoplot[0], cax=ax[0, 3], ticks=[0, 1, 2])
geocb.set_ticklabels(["BCK", "PK", "VK"])
geocb.set_label("Quasi-Geology model\n(Rock units classification)", fontsize=16)
ax[0, 3].set_aspect(10)
# gravity model
mesh.plot_slice(
density_model_no_info,
normal="X",
ax=ax[1, 0],
clim=[-1, 0],
ind=indx,
pcolor_opts={"cmap": "Blues_r"},
)
mesh.plot_slice(
density_model_no_info,
normal="Y",
ax=ax[1, 1],
clim=[-1, 0],
ind=indy,
pcolor_opts={"cmap": "Blues_r"},
)
denplot = mesh.plot_slice(
density_model_no_info,
normal="Z",
ax=ax[1, 2],
clim=[-1, 0],
ind=indz,
pcolor_opts={"cmap": "Blues_r"},
)
dencb = plt.colorbar(denplot[0], cax=ax[1, 3])
dencb.set_label("Density contrast\nmodel (g/cc)", fontsize=16)
ax[1, 3].set_aspect(10)
# magnetic model
mesh.plot_slice(
magsus_model_no_info,
normal="X",
ax=ax[2, 0],
clim=[0, 0.025],
ind=indx,
pcolor_opts={"cmap": "Reds"},
)
mesh.plot_slice(
magsus_model_no_info,
normal="Y",
ax=ax[2, 1],
clim=[0, 0.025],
ind=indy,
pcolor_opts={"cmap": "Reds"},
)
susplot = mesh.plot_slice(
magsus_model_no_info,
normal="Z",
ax=ax[2, 2],
clim=[0, 0.025],
ind=indz,
pcolor_opts={"cmap": "Reds"},
)
suscb = plt.colorbar(susplot[0], cax=ax[2, 3])
suscb.set_label("Magnetic susceptibility\nmodel (SI)", fontsize=16)
ax[2, 3].set_aspect(10)
# overlay true geology model for comparison
indslicexplot = mesh.gridCC[:, 0] == mesh.cell_centers_x[indx]
indsliceyplot = mesh.gridCC[:, 1] == mesh.cell_centers_y[indy]
indslicezplot = mesh.gridCC[:, 2] == mesh.cell_centers_z[indz]
for i in range(3):
for j, (plane, indd) in enumerate(
zip([[1, 2], [0, 2], [0, 1]], [indslicexplot, indsliceyplot, indslicezplot])
):
utils.plot2Ddata(
mesh.gridCC[indd][:, plane],
true_geology[indd],
nx=100,
ny=100,
contourOpts={"alpha": 0},
clim=[0, 2],
ax=ax[i, j],
level=True,
ncontour=2,
levelOpts={"colors": "grey", "linewidths": 2, "linestyles": "--"},
method="nearest",
)
# plot the locations of the cross-sections
for i in range(3):
ax[i, 0].plot(
mesh.cell_centers_y[indy] * np.ones(2), [-300, 500], c="k", linestyle="dotted"
)
ax[i, 0].plot(
[
data_mag.survey.receiver_locations[:, 1].min(),
data_mag.survey.receiver_locations[:, 1].max(),
],
mesh.cell_centers_z[indz] * np.ones(2),
c="k",
linestyle="dotted",
)
ax[i, 0].set_xlim(
[
data_mag.survey.receiver_locations[:, 1].min(),
data_mag.survey.receiver_locations[:, 1].max(),
],
)
ax[i, 1].plot(
mesh.cell_centers_x[indx] * np.ones(2), [-300, 500], c="k", linestyle="dotted"
)
ax[i, 1].plot(
[
data_mag.survey.receiver_locations[:, 0].min(),
data_mag.survey.receiver_locations[:, 0].max(),
],
mesh.cell_centers_z[indz] * np.ones(2),
c="k",
linestyle="dotted",
)
ax[i, 1].set_xlim(
[
data_mag.survey.receiver_locations[:, 0].min(),
data_mag.survey.receiver_locations[:, 0].max(),
],
)
ax[i, 2].plot(
mesh.cell_centers_x[indx] * np.ones(2),
[
data_mag.survey.receiver_locations[:, 1].min(),
data_mag.survey.receiver_locations[:, 1].max(),
],
c="k",
linestyle="dotted",
)
ax[i, 2].plot(
[
data_mag.survey.receiver_locations[:, 0].min(),
data_mag.survey.receiver_locations[:, 0].max(),
],
mesh.cell_centers_y[indy] * np.ones(2),
c="k",
linestyle="dotted",
)
ax[i, 2].set_xlim(
[
data_mag.survey.receiver_locations[:, 0].min(),
data_mag.survey.receiver_locations[:, 0].max(),
],
)
ax[i, 2].set_ylim(
[
data_mag.survey.receiver_locations[:, 1].min(),
data_mag.survey.receiver_locations[:, 1].max(),
],
)
plt.tight_layout()
plt.show()
# Plot the learned 2D GMM
fig = plt.figure(figsize=(10, 10))
ax0 = plt.subplot2grid((4, 4), (3, 1), colspan=3)
ax1 = plt.subplot2grid((4, 4), (0, 1), colspan=3, rowspan=3)
ax2 = plt.subplot2grid((4, 4), (0, 0), rowspan=3)
ax = [ax0, ax1, ax2]
learned_gmm.plot_pdf(flag2d=True, ax=ax, padding=1, plotting_precision=100)
ax[0].set_xlabel("Density contrast [g/cc]")
ax[0].set_ylim([0, 5])
ax[2].set_xlim([0, 50])
ax[2].set_ylabel("magnetic Susceptibility [SI]")
ax[1].scatter(
density_model_no_info[actv],
magsus_model_no_info[actv],
c=quasi_geology_model_no_info[actv],
cmap="inferno_r",
edgecolors="k",
label="recovered PGI model",
alpha=0.5,
)
ax[0].hist(density_model_no_info[actv], density=True, bins=50)
ax[2].hist(magsus_model_no_info[actv], density=True, bins=50, orientation="horizontal")
ax[1].scatter(
[0, -0.8, -0.02],
[0, 0.005, 0.02],
label="True petrophysical means",
cmap="inferno_r",
c=[0, 1, 2],
marker="v",
edgecolors="k",
s=200,
)
ax[1].legend()
plt.show()
Running inversion with SimPEG v0.25.0
Alpha scales: [np.float64(9293029.056720901), np.float64(0.0), np.float64(7011631.545233198), np.float64(0.0), np.float64(13519145.111732943), np.float64(0.0), np.float64(979157391.2547016), np.float64(0.0), np.float64(720315070.5979666), np.float64(0.0), np.float64(1710373999.0876644), np.float64(0.0)]
<class 'simpeg.regularization.pgi.PGIsmallness'>
Initial data misfit scales: [0.98146488 0.01853512]
================================================= 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.56e-08 4.27e+06 9.99e+04 4.27e+06 0 inf inf
1 1.56e-08 4.08e+05 2.94e+09 4.08e+05 2.05e+02 0 100 4.13e-03 6.64e+03
geophys. misfits: 413393.8 (target 576.0 [False]); 126047.7 (target 576.0 [False]) | smallness misfit: 130560.2 (target: 11719.0 [False])
Beta cooling evaluation: progress: [413393.8 126047.7]; minimum progress targets: [3470231.3 572697.8]
mref changed in 4737 places
2 1.56e-08 1.83e+05 5.33e+09 1.83e+05 2.05e+01 0 100 6.35e-03 2.96e+03
geophys. misfits: 184995.6 (target 576.0 [False]); 58766.6 (target 576.0 [False]) | smallness misfit: 117051.8 (target: 11719.0 [False])
Beta cooling evaluation: progress: [184995.6 58766.6]; minimum progress targets: [330715. 100838.2]
mref changed in 2582 places
3 1.56e-08 7.00e+04 5.44e+09 7.01e+04 1.96e+01 0 100 1.48e-02 2.31e+03
geophys. misfits: 71131.9 (target 576.0 [False]); 9795.2 (target 576.0 [False]) | smallness misfit: 90884.6 (target: 11719.0 [False])
Beta cooling evaluation: progress: [71131.9 9795.2]; minimum progress targets: [147996.4 47013.3]
mref changed in 1039 places
4 1.56e-08 4.16e+03 5.72e+09 4.25e+03 1.83e+01 0 100 1.12e-02 8.53e+02
geophys. misfits: 4230.2 (target 576.0 [False]); 446.8 (target 576.0 [True]) | smallness misfit: 69477.8 (target: 11719.0 [False])
Updating scaling for data misfits by 1.2892515021215902
New scales: [0.98556331 0.01443669]
Beta cooling evaluation: progress: [4230.2 446.8]; minimum progress targets: [56905.5 7836.2]
mref changed in 567 places
5 1.56e-08 1.86e+02 5.77e+09 2.76e+02 1.79e+01 0 100 3.21e-02 5.44e+02
geophys. misfits: 183.9 (target 576.0 [True]); 349.1 (target 576.0 [True]) | smallness misfit: 63846.5 (target: 11719.0 [False])
Beta cooling evaluation: progress: [183.9 349.1]; minimum progress targets: [3384.1 691.2]
Warming alpha_pgi to favor clustering: 2.390829629260463
mref changed in 132 places
6 1.56e-08 1.81e+02 5.87e+09 2.73e+02 2.45e+01 0 100 1.66e-01 7.07e+02
geophys. misfits: 181.8 (target 576.0 [True]); 140.6 (target 576.0 [True]) | smallness misfit: 62391.0 (target: 11719.0 [False])
Beta cooling evaluation: progress: [181.8 140.6]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 8.685359106520854
mref changed in 35 places
7 1.56e-08 1.73e+02 6.34e+09 2.72e+02 3.77e+01 1 100 1.24e-01 3.05e+02
geophys. misfits: 154.3 (target 576.0 [True]); 1482.3 (target 576.0 [False]) | smallness misfit: 61319.8 (target: 11719.0 [False])
Updating scaling for data misfits by 3.732992756268878
New scales: [0.94815358 0.05184642]
Beta cooling evaluation: progress: [ 154.3 1482.3]; minimum progress targets: [691.2 691.2]
Decreasing beta to counter data misfit increase.
mref changed in 69 places
8 7.80e-09 2.19e+02 6.30e+09 2.68e+02 1.97e+01 4 100 7.61e-02 2.78e+03
geophys. misfits: 159.8 (target 576.0 [True]); 1302.1 (target 576.0 [False]) | smallness misfit: 61021.0 (target: 11719.0 [False])
Updating scaling for data misfits by 3.603862454545514
New scales: [0.8353769 0.1646231]
Beta cooling evaluation: progress: [ 159.8 1302.1]; minimum progress targets: [ 691.2 1185.8]
Decreasing beta to counter data misfit increase.
mref changed in 5 places
9 3.90e-09 2.83e+02 6.71e+09 3.09e+02 1.80e+01 2 100 5.55e-02 5.97e+03
geophys. misfits: 160.6 (target 576.0 [True]); 904.1 (target 576.0 [False]) | smallness misfit: 60199.1 (target: 11719.0 [False])
Updating scaling for data misfits by 3.5871586491331064
New scales: [0.58585689 0.41414311]
Beta cooling evaluation: progress: [160.6 904.1]; minimum progress targets: [ 691.2 1041.7]
Decreasing beta to counter data misfit increase.
mref changed in 81 places
10 1.95e-09 3.90e+02 7.16e+09 4.04e+02 1.78e+01 2 100 2.77e-02 6.11e+03
geophys. misfits: 174.5 (target 576.0 [True]); 694.1 (target 576.0 [False]) | smallness misfit: 60451.9 (target: 11719.0 [False])
Updating scaling for data misfits by 3.3008704977368515
New scales: [0.29999488 0.70000512]
Beta cooling evaluation: progress: [174.5 694.1]; minimum progress targets: [691.2 723.3]
Decreasing beta to counter data misfit increase.
mref changed in 67 places
11 9.76e-10 4.94e+02 7.27e+09 5.01e+02 1.77e+01 3 100 2.11e-02 6.58e+03
geophys. misfits: 159.2 (target 576.0 [True]); 638.0 (target 576.0 [False]) | smallness misfit: 60711.5 (target: 11719.0 [False])
Updating scaling for data misfits by 3.6173313169513883
New scales: [0.10592496 0.89407504]
Beta cooling evaluation: progress: [159.2 638. ]; minimum progress targets: [691.2 691.2]
mref changed in 29 places
12 9.76e-10 1.08e+02 7.19e+09 1.15e+02 1.76e+01 0 100 7.29e-03 5.06e+03
geophys. misfits: 103.2 (target 576.0 [True]); 108.9 (target 576.0 [True]) | smallness misfit: 62523.1 (target: 11719.0 [False])
Beta cooling evaluation: progress: [103.2 108.9]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 47.22597670195155
mref changed in 135 places
13 9.76e-10 2.78e+01 1.00e+10 3.76e+01 8.99e+00 0 100 1.29e-02 3.37e+03
geophys. misfits: 101.1 (target 576.0 [True]); 19.2 (target 576.0 [True]) | smallness misfit: 147752.9 (target: 11719.0 [False])
Beta cooling evaluation: progress: [101.1 19.2]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 844.2478637435512
mref changed in 678 places
14 9.76e-10 2.77e+01 1.96e+11 2.19e+02 1.70e+01 7 100 9.46e-02 6.44e+03
/t40array/ssoler/miniforge3/envs/simpeg-user-tutorials/lib/python3.13/site-packages/simpeg/utils/pgi_utils.py:1110: ConvergenceWarning: Initialization 10 did not converge. Try different init parameters, or increase max_iter, tol or check for degenerate data.
self.fit_predict(X, y, debug)
geophys. misfits: 101.3 (target 576.0 [True]); 19.0 (target 576.0 [True]) | smallness misfit: 147379.3 (target: 11719.0 [False])
Beta cooling evaluation: progress: [101.3 19. ]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 15192.79572692244
mref changed in 3 places
15 9.76e-10 2.35e+02 2.47e+12 2.65e+03 1.69e+01 0 100 8.65e-02 6.64e+03
geophys. misfits: 1145.9 (target 576.0 [False]); 126.9 (target 576.0 [True]) | smallness misfit: 55132.3 (target: 11719.0 [False])
Updating scaling for data misfits by 4.539612709965989
New scales: [0.34973208 0.65026792]
Beta cooling evaluation: progress: [1145.9 126.9]; minimum progress targets: [691.2 691.2]
Decreasing beta to counter data misfit increase.
mref changed in 916 places
16 4.88e-10 4.69e+02 5.66e+11 7.45e+02 1.09e+02 4 100 2.26e-02 2.72e+03
geophys. misfits: 1056.4 (target 576.0 [False]); 152.6 (target 576.0 [True]) | smallness misfit: 53389.0 (target: 11719.0 [False])
Updating scaling for data misfits by 3.775314592811509
New scales: [0.67001803 0.32998197]
Beta cooling evaluation: progress: [1056.4 152.6]; minimum progress targets: [916.7 691.2]
Decreasing beta to counter data misfit increase.
mref changed in 5 places
17 2.44e-10 5.15e+02 4.44e+11 6.24e+02 1.20e+02 1 100 2.22e-02 3.44e+03
geophys. misfits: 621.9 (target 576.0 [False]); 299.2 (target 576.0 [True]) | smallness misfit: 43673.2 (target: 11719.0 [False])
Updating scaling for data misfits by 1.9253104331545225
New scales: [0.79630424 0.20369576]
Beta cooling evaluation: progress: [621.9 299.2]; minimum progress targets: [845.1 691.2]
mref changed in 10 places
18 2.44e-10 4.24e+02 4.43e+11 5.32e+02 2.22e+01 1 100 2.83e-02 2.20e+03
geophys. misfits: 408.3 (target 576.0 [True]); 485.5 (target 576.0 [True]) | smallness misfit: 40043.9 (target: 11719.0 [False])
Beta cooling evaluation: progress: [408.3 485.5]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 19728.276766778865
mref changed in 11 places
19 2.44e-10 3.66e+02 5.42e+11 4.98e+02 1.74e+01 1 100 3.58e-02 2.92e+03
geophys. misfits: 302.8 (target 576.0 [True]); 613.7 (target 576.0 [False]) | smallness misfit: 37432.2 (target: 11719.0 [False])
Updating scaling for data misfits by 1.9022252177664045
New scales: [0.67267956 0.32732044]
Beta cooling evaluation: progress: [302.8 613.7]; minimum progress targets: [691.2 691.2]
mref changed in 7 places
20 2.44e-10 3.83e+02 5.24e+11 5.11e+02 1.69e+01 2 100 2.54e-02 3.62e+03
geophys. misfits: 278.3 (target 576.0 [True]); 597.9 (target 576.0 [False]) | smallness misfit: 36598.1 (target: 11719.0 [False])
Updating scaling for data misfits by 2.069544577770102
New scales: [0.49825023 0.50174977]
Beta cooling evaluation: progress: [278.3 597.9]; minimum progress targets: [691.2 691.2]
mref changed in 2 places
21 2.44e-10 3.03e+02 4.79e+11 4.20e+02 1.68e+01 1 100 2.21e-02 4.85e+03
geophys. misfits: 262.5 (target 576.0 [True]); 342.5 (target 576.0 [True]) | smallness misfit: 35050.1 (target: 11719.0 [False])
Beta cooling evaluation: progress: [262.5 342.5]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 38236.26604250299
mref changed in 3 places
22 2.44e-10 2.69e+02 8.51e+11 4.76e+02 1.78e+01 2 100 4.05e-02 6.59e+03
geophys. misfits: 273.6 (target 576.0 [True]); 263.5 (target 576.0 [True]) | smallness misfit: 33752.0 (target: 11719.0 [False])
Beta cooling evaluation: progress: [273.6 263.5]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 82047.96442282462
mref changed in 0 places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
23 2.44e-10 3.43e+02 1.14e+12 6.22e+02 2.15e+01 0 100 6.10e-02 8.87e+03
geophys. misfits: 416.4 (target 576.0 [True]); 270.4 (target 576.0 [True]) | smallness misfit: 28003.5 (target: 11719.0 [False])
Beta cooling evaluation: progress: [416.4 270.4]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 144140.09622404317
mref changed in 0 places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
24 2.44e-10 2.67e+02 1.97e+12 7.48e+02 5.68e+01 3 100 4.13e-02 5.30e+03
geophys. misfits: 425.7 (target 576.0 [True]); 109.8 (target 576.0 [True]) | smallness misfit: 27764.7 (target: 11719.0 [False])
Beta cooling evaluation: progress: [425.7 109.8]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 475606.8903409355
mref changed in 0 places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
25 2.44e-10 4.71e+02 4.67e+12 1.61e+03 6.60e+01 1 100 9.04e-02 1.04e+04
geophys. misfits: 608.9 (target 576.0 [False]); 333.2 (target 576.0 [True]) | smallness misfit: 24261.2 (target: 11719.0 [False])
Updating scaling for data misfits by 1.7288454822208763
New scales: [0.63191819 0.36808181]
Beta cooling evaluation: progress: [608.9 333.2]; minimum progress targets: [691.2 691.2]
mref changed in 0 places
26 2.44e-10 5.79e+02 3.80e+12 1.51e+03 8.02e+01 0 100 5.16e-02 5.82e+03
geophys. misfits: 790.2 (target 576.0 [False]); 215.7 (target 576.0 [True]) | smallness misfit: 22804.4 (target: 11719.0 [False])
Updating scaling for data misfits by 2.6702026691996608
New scales: [0.82092236 0.17907764]
Beta cooling evaluation: progress: [790.2 215.7]; minimum progress targets: [691.2 691.2]
Decreasing beta to counter data misfit increase.
mref changed in 0 places
27 1.22e-10 5.94e+02 4.29e+12 1.12e+03 1.23e+02 1 100 5.48e-02 3.00e+03
geophys. misfits: 626.5 (target 576.0 [False]); 444.0 (target 576.0 [True]) | smallness misfit: 22809.7 (target: 11719.0 [False])
Updating scaling for data misfits by 1.2972686495349859
New scales: [0.85605091 0.14394909]
Beta cooling evaluation: progress: [626.5 444. ]; minimum progress targets: [691.2 691.2]
mref changed in 0 places
28 1.22e-10 5.29e+02 4.56e+12 1.08e+03 4.98e+01 1 100 4.91e-02 8.67e+03
geophys. misfits: 560.1 (target 576.0 [True]); 343.1 (target 576.0 [True]) | smallness misfit: 23523.7 (target: 11719.0 [False])
Beta cooling evaluation: progress: [560.1 343.1]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 643825.1964467242
mref changed in 0 places
Add mref to Smoothness. Changes in mref happened in 0.0 % of the cells
29 1.22e-10 5.06e+02 6.08e+12 1.25e+03 1.67e+01 0 100 5.41e-02 4.88e+03
geophys. misfits: 562.2 (target 576.0 [True]); 171.2 (target 576.0 [True]) | smallness misfit: 23635.4 (target: 11719.0 [False])
Beta cooling evaluation: progress: [562.2 171.2]; minimum progress targets: [691.2 691.2]
Warming alpha_pgi to favor clustering: 1412829.8824339863
mref changed in 1 places
30 1.22e-10 8.33e+02 1.04e+13 2.10e+03 3.65e+01 0 100 1.72e-01 6.92e+03
geophys. misfits: 829.4 (target 576.0 [False]); 854.2 (target 576.0 [False]) | smallness misfit: 22642.3 (target: 11719.0 [False])
Beta cooling evaluation: progress: [829.4 854.2]; minimum progress targets: [691.2 691.2]
Decreasing beta to counter data misfit increase.
mref changed in 0 places
------------------------- STOP! -------------------------
1 : |fc-fOld| = 6.6429e+02 <= tolF*(1+|f0|) = 4.2707e+05
0 : |xc-x_last| = 3.2552e-01 <= tolX*(1+|x0|) = 1.0153e-01
0 : |proj(x-g)-x| = 1.0735e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 1.0735e+02 <= 1e3*eps = 1.0000e-02
1 : maxIter = 30 <= iter = 30
------------------------- DONE! -------------------------
