Note
Go to the end to download the full example code.
Invert three surveys line for a thin plate using the surrogate model#
If you are running this notebook locally, make sure you’ve followed steps here to set up the environment. (This environment.yml file specifies a list of packages required to run the notebooks)
This notebook assumes that you have created a surrogate model by executing the following two notebooks: - Latin Hypercube Sampling - Surrogate model creation
# -------------------------------------------------------- #
# #
# Uncomment below to set up environment on "colab" #
# #
# -------------------------------------------------------- #
# !pip install -U cofi
# !pip install git+https://github.com/JuergHauser/PyP223.git
# !pip install smt
# !git clone https://github.com/inlab-geo/cofi-examples.git
# %cd cofi-examples/examples/vtem_max
# If this notebook is run locally PyP223 and smt need to be installed separately by uncommenting the following lines,
# that is by removing the # and the white space between it and the exclamation mark.
# !pip install git+https://github.com/JuergHauser/PyP223.git
# !pip install smt
import pickle
import numpy
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import cofi
import arviz
from vtem_max_forward_lib import (
problem_setup,
system_spec,
survey_setup,
ForwardWrapper,
plot_predicted_profile,
plot_transient,
plot_plate_faces,
plot_plate_faces_single
)
numpy.random.seed(42)
Background#
This example inverts three survey line of VTEM max data using the vertical component for a thin plate target. It thus becomes possible to invert for the easting,northing, depth of the plate reference point, the plate dip and plate azimuth. Solving the forward problem, that is calculating the objective function, usess the surrogate model that has been created by the Kriging approach applied to the latin hypercube samples of the objective function.
Problem definition#
tx_min = 115
tx_max = 281
tx_interval = 15
ty_min = 25
ty_max = 176
ty_interval = 75
tx_points = numpy.arange(tx_min, tx_max, tx_interval)
ty_points = numpy.arange(ty_min, ty_max, ty_interval)
n_transmitters = len(tx_points) * len(ty_points)
tx, ty = numpy.meshgrid(tx_points, ty_points)
tx = tx.flatten()
ty = ty.flatten()
fiducial_id = numpy.arange(len(tx))
line_id = numpy.zeros(len(tx), dtype=int)
line_id[ty==ty_points[0]] = 0
line_id[ty==ty_points[1]] = 1
line_id[ty==ty_points[2]] = 2
survey_setup = {
"tx": tx, # transmitter easting/x-position
"ty": ty, # transmitter northing/y-position
"tz": numpy.array([50]*n_transmitters), # transmitter height/z-position
"tazi": numpy.deg2rad(numpy.array([90]*n_transmitters)), # transmitter azimuth
"tincl": numpy.deg2rad(numpy.array([6]*n_transmitters)), # transmitter inclination
"rx": tx, # receiver easting/x-position
"ry": numpy.array([100]*n_transmitters), # receiver northing/y-position
"rz": numpy.array([50]*n_transmitters), # receiver height/z-position
"trdx": numpy.array([0]*n_transmitters), # transmitter receiver separation inline
"trdy": numpy.array([0]*n_transmitters), # transmitter receiver separation crossline
"trdz": numpy.array([0]*n_transmitters), # transmitter receiver separation vertical
"fiducial_id": fiducial_id, # unique id for each transmitter
"line_id": line_id # id for each line
}
true_model = {
"res": numpy.array([300, 1000]),
"thk": numpy.array([20]),
"peast": numpy.array([175]),
"pnorth": numpy.array([100]),
"ptop": numpy.array([30]),
"pres": numpy.array([0.1]),
"plngth1": numpy.array([100]),
"plngth2": numpy.array([100]),
"pwdth1": numpy.array([0.1]),
"pwdth2": numpy.array([90]),
"pdzm": numpy.array([75]),
"pdip": numpy.array([60])
}
forward = ForwardWrapper(true_model, problem_setup, system_spec, survey_setup,
["pdip","pdzm", "peast", "ptop", "pwdth2"], data_returned=["vertical"])
['pdip', 'pdzm', 'peast', 'ptop', 'pwdth2']
# check the order of parameters in a model vector
forward.params_to_invert
['pdip', 'pdzm', 'peast', 'ptop', 'pwdth2']
true_param_value = numpy.array([60.,65., 175., 30., 90.])
Ensemble method using the surrogate model#
sm = None
filename = "kriging_surrogate_model.pkl"
with open(filename, "rb") as f:
sm = pickle.load(f)
Initialise a model for inversion
init_param_value = numpy.array([45, 90, 160, 35, 80])
m_min = numpy.array([15, 35, 155, 30, 65])
m_max = numpy.array([75, 145, 185, 40, 115])
Define helper functions for CoFI
def my_objective(model):
val=sm.predict_values(numpy.array([model]))[0][0]
if val<1e-3:
return 1e-3
else:
return val
def my_log_likelihood(model):
return -0.5 * my_objective(model)
def my_log_prior(model): # uniform distribution
for i in range(len(model)):
if model[i] < m_min[i] or model[i] > m_max[i]: return -numpy.inf
return 0.0 # model lies within bounds -> return log(1)
Define CoFI problem
my_problem = cofi.BaseProblem()
my_problem.set_objective(my_objective)
Define CoFI options
inv_options = cofi.InversionOptions()
inv_options.set_tool("neighpy")
inv_options.suggest_solver_params()
Current backend tool neighpy has the following solver-specific parameters:
Required parameters:
{'n_cells_to_resample', 'n_initial_samples', 'n_walkers', 'bounds', 'n_resample', 'n_iterations', 'n_samples_per_iteration'}
Optional parameters & default settings:
{'serial': False}
inv_options.set_params(
n_samples_per_iteration=100,
n_initial_samples=10,
n_resample=8000,
n_iterations=100,
bounds=numpy.array([m_min, m_max]).T,
n_cells_to_resample=10,
n_walkers=4
)
######## Run it
inv = cofi.Inversion(my_problem, inv_options)
inv_result = inv.run()
######## Check result
print(f"The inversion result from `neighpy`:")
inv_result.summary()
NAI - Initial Random Search
NAI - Optimisation Loop: 0%| | 0/100 [00:00<?, ?it/s]
NAI - Optimisation Loop: 1%| | 1/100 [00:00<00:36, 2.73it/s]
NAI - Optimisation Loop: 3%|▎ | 3/100 [00:00<00:13, 7.46it/s]
NAI - Optimisation Loop: 5%|▌ | 5/100 [00:00<00:08, 10.62it/s]
NAI - Optimisation Loop: 7%|▋ | 7/100 [00:00<00:07, 12.86it/s]
NAI - Optimisation Loop: 9%|▉ | 9/100 [00:00<00:06, 14.38it/s]
NAI - Optimisation Loop: 11%|█ | 11/100 [00:00<00:05, 15.48it/s]
NAI - Optimisation Loop: 13%|█▎ | 13/100 [00:01<00:05, 16.19it/s]
NAI - Optimisation Loop: 15%|█▌ | 15/100 [00:01<00:05, 16.64it/s]
NAI - Optimisation Loop: 17%|█▋ | 17/100 [00:01<00:04, 16.87it/s]
NAI - Optimisation Loop: 19%|█▉ | 19/100 [00:01<00:04, 17.10it/s]
NAI - Optimisation Loop: 21%|██ | 21/100 [00:01<00:04, 17.65it/s]
NAI - Optimisation Loop: 23%|██▎ | 23/100 [00:01<00:04, 17.72it/s]
NAI - Optimisation Loop: 25%|██▌ | 25/100 [00:01<00:04, 17.70it/s]
NAI - Optimisation Loop: 27%|██▋ | 27/100 [00:01<00:04, 17.67it/s]
NAI - Optimisation Loop: 29%|██▉ | 29/100 [00:01<00:04, 17.61it/s]
NAI - Optimisation Loop: 31%|███ | 31/100 [00:02<00:03, 17.50it/s]
NAI - Optimisation Loop: 33%|███▎ | 33/100 [00:02<00:03, 17.76it/s]
NAI - Optimisation Loop: 35%|███▌ | 35/100 [00:02<00:03, 17.66it/s]
NAI - Optimisation Loop: 37%|███▋ | 37/100 [00:02<00:03, 17.50it/s]
NAI - Optimisation Loop: 39%|███▉ | 39/100 [00:02<00:03, 17.25it/s]
NAI - Optimisation Loop: 41%|████ | 41/100 [00:02<00:03, 17.20it/s]
NAI - Optimisation Loop: 43%|████▎ | 43/100 [00:02<00:03, 17.24it/s]
NAI - Optimisation Loop: 45%|████▌ | 45/100 [00:02<00:03, 17.27it/s]
NAI - Optimisation Loop: 47%|████▋ | 47/100 [00:02<00:03, 17.23it/s]
NAI - Optimisation Loop: 49%|████▉ | 49/100 [00:03<00:02, 17.22it/s]
NAI - Optimisation Loop: 51%|█████ | 51/100 [00:03<00:02, 17.12it/s]
NAI - Optimisation Loop: 53%|█████▎ | 53/100 [00:03<00:02, 17.04it/s]
NAI - Optimisation Loop: 55%|█████▌ | 55/100 [00:03<00:02, 16.99it/s]
NAI - Optimisation Loop: 57%|█████▋ | 57/100 [00:03<00:02, 17.06it/s]
NAI - Optimisation Loop: 59%|█████▉ | 59/100 [00:03<00:02, 17.07it/s]
NAI - Optimisation Loop: 61%|██████ | 61/100 [00:03<00:02, 16.96it/s]
NAI - Optimisation Loop: 63%|██████▎ | 63/100 [00:03<00:02, 16.87it/s]
NAI - Optimisation Loop: 65%|██████▌ | 65/100 [00:04<00:02, 16.75it/s]
NAI - Optimisation Loop: 67%|██████▋ | 67/100 [00:04<00:01, 16.70it/s]
NAI - Optimisation Loop: 69%|██████▉ | 69/100 [00:04<00:01, 16.76it/s]
NAI - Optimisation Loop: 71%|███████ | 71/100 [00:04<00:01, 16.59it/s]
NAI - Optimisation Loop: 73%|███████▎ | 73/100 [00:04<00:01, 16.59it/s]
NAI - Optimisation Loop: 75%|███████▌ | 75/100 [00:04<00:01, 16.58it/s]
NAI - Optimisation Loop: 77%|███████▋ | 77/100 [00:04<00:01, 16.68it/s]
NAI - Optimisation Loop: 79%|███████▉ | 79/100 [00:04<00:01, 16.60it/s]
NAI - Optimisation Loop: 81%|████████ | 81/100 [00:05<00:01, 16.22it/s]
NAI - Optimisation Loop: 83%|████████▎ | 83/100 [00:05<00:01, 16.37it/s]
NAI - Optimisation Loop: 85%|████████▌ | 85/100 [00:05<00:00, 16.60it/s]
NAI - Optimisation Loop: 87%|████████▋ | 87/100 [00:05<00:00, 16.54it/s]
NAI - Optimisation Loop: 89%|████████▉ | 89/100 [00:05<00:00, 16.56it/s]
NAI - Optimisation Loop: 91%|█████████ | 91/100 [00:05<00:00, 16.75it/s]
NAI - Optimisation Loop: 93%|█████████▎| 93/100 [00:05<00:00, 16.75it/s]
NAI - Optimisation Loop: 95%|█████████▌| 95/100 [00:05<00:00, 16.33it/s]
NAI - Optimisation Loop: 97%|█████████▋| 97/100 [00:05<00:00, 15.99it/s]
NAI - Optimisation Loop: 99%|█████████▉| 99/100 [00:06<00:00, 15.51it/s]
NAI - Optimisation Loop: 100%|██████████| 100/100 [00:06<00:00, 16.18it/s]
The inversion result from `neighpy`:
============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [ 60.54890045 68.21943427 179.34874008 38.03067502 102.4513261 ]
direct_search_samples: [[ 56.6912376 61.72943964 165.30224057 39.25708288 65.24360417]
[ 64.07747464 94.20427513 166.36617511 32.13264507 111.92329337]
[ 39.33555312 55.25951694 156.03017994 39.42022733 84.52241171]
...
[ 60.32016867 68.21948407 179.35088011 38.03887313 102.45475841]
[ 60.32017207 68.21948739 179.35087916 38.03887312 102.45475954]
[ 60.32017157 68.21948501 179.35088322 38.03887402 102.45475897]]
direct_search_objectives: [4.12243928e+04 2.48002996e+04 3.03365076e+04 ... 1.00000000e-03
1.00000000e-03 1.00000000e-03]
appraisal_samples: [[ 56.72072038 52.22306732 176.15083527 38.19021261 111.16926771]
[ 55.22087425 49.68086801 176.86884413 38.41193374 110.4078104 ]
[ 56.62957613 55.41558994 173.9893494 35.65224388 108.4636213 ]
...
[ 45.08009491 61.27801175 181.71420978 38.51149541 97.89571622]
[ 52.81663461 50.63049815 184.57462469 34.88042899 98.21648662]
[ 53.89568914 60.7444161 181.90052107 33.99459004 107.94323154]]
Plotting#
import arviz_base
arviz.style.use("arviz-variat")
display_names = [
"Dip (°)",
"Dip azimuth (°)",
"Easting (m)",
"Depth (m)",
"Width (m)"
]
clean_names = ["Dip_deg", "Dip_azimuth_deg", "Easting_m", "Depth_m", "Width_m"]
true_values = [60, 65, 175, 30, 90]
d = {k: v[numpy.newaxis, :] for k, v in zip(clean_names, inv_result.appraisal_samples.T)}
az_idata = arviz_base.from_dict({"posterior": d})
pc = arviz.plot_trace_dist(
az_idata.sel(draw=slice(2000, None)),
visuals={"xlabel_trace": False, "trace": {"color": "C0"}, "dist": {"color": "C0"}},
figure_kwargs={"figsize": (12, 20), "constrained_layout": True},
)
for i, dname in enumerate(display_names):
ax_kde = pc.iget_target(i, 0)
ax_trace = pc.iget_target(i, 1)
ax_kde.set_title(dname)
ax_trace.set_title(dname)
ax_kde.axvline(true_values[i], color="green", linestyle="--", lw=1, alpha=0.5)
ax_trace.axhline(true_values[i], color="green", linestyle="--", lw=1, alpha=0.5)
ax_trace.margins(x=0)
from scipy.stats import gaussian_kde
true_values_dict = {
display_names[i]: true_param_value[i] for i in range(init_param_value.size)
}
# Set to True for KDE contours (old style), False for scatter (arviz 1.0 default)
USE_KDE_CONTOURS = True
arviz_base.rcParams["plot.max_subplots"] = 80
if USE_KDE_CONTOURS:
pm = arviz.plot_pair(
az_idata.sel(draw=slice(4000, None)),
marginal=True,
triangle="lower",
visuals={"scatter": False},
)
# Add KDE contours to off-diagonal panels
posterior = az_idata.posterior.sel(draw=slice(4000, None))
n = len(clean_names)
for i in range(n):
for j in range(n):
if i <= j:
continue
try:
ax = pm.iget_target(i, j)
except (ValueError, IndexError):
continue
x = posterior[clean_names[j]].values.flatten()
y = posterior[clean_names[i]].values.flatten()
kde = gaussian_kde(numpy.vstack([x, y]))
xmin, xmax = x.min(), x.max()
ymin, ymax = y.min(), y.max()
xx, yy = numpy.mgrid[xmin:xmax:100j, ymin:ymax:100j]
zz = kde(numpy.vstack([xx.ravel(), yy.ravel()])).reshape(xx.shape)
ax.contourf(xx, yy, zz, levels=10, cmap="Blues")
ax.contour(xx, yy, zz, levels=10, colors="grey", linewidths=0.5, alpha=0.5)
else:
pm = arviz.plot_pair(
az_idata.sel(draw=slice(4000, None)),
marginal=True,
triangle="lower",
)
# Add reference values
ref_vals = list(true_values_dict.values())
n = len(ref_vals)
for i in range(n):
for j in range(n):
try:
ax = pm.iget_target(i, j)
except (ValueError, IndexError):
continue
if i == j:
ax.axvline(ref_vals[i], color="green", linestyle="--", lw=1, alpha=0.5)
elif i > j:
ax.plot(
ref_vals[j], ref_vals[i], "o",
color="yellow", markeredgecolor="k", ms=10, zorder=5,
)
arviz.style.use("arviz-variat")
_, axes = plt.subplots(2, 2)
axes[1,1].axis("off")
plot_plate_faces(
"plate_true", forward, true_param_value,
axes[0,0], axes[0,1], axes[1,0], color="purple", label="True model"
)
plot_plate_faces(
"plate_init", forward, init_param_value,
axes[0,0], axes[0,1], axes[1,0], color="green", label="Starting model"
)
plt.tight_layout()
posterior = az_idata.posterior
has_chain = "chain" in posterior.dims
n_chains = int(posterior.sizes["chain"]) if has_chain else 1
n_draws = int(posterior.sizes["draw"])
ichain = 0
idraw = min(2500, n_draws - 1)
sample = numpy.zeros(5)
if has_chain:
for idx, cn in enumerate(clean_names):
sample[idx] = float(posterior[cn].isel(chain=ichain, draw=idraw))
else:
for idx, cn in enumerate(clean_names):
sample[idx] = float(posterior[cn].isel(draw=idraw))
plot_plate_faces(
"plate_inverted", forward, sample,
axes[0,0], axes[0,1], axes[1,0], color="red", label="Posterior sample", linestyle="dotted"
)
point = Line2D([0], [0], label='Fiducial', marker='o', markersize=5,
markeredgecolor='orange', markerfacecolor='orange', linestyle='')
handles, labels = axes[1,0].get_legend_handles_labels()
handles.extend([point])
axes[1,0].legend(handles=handles,bbox_to_anchor=(1.04, 0), loc="lower left")
# plot 10 randomly selected samples of the posterior distribution
for i in range(10):
idraw = numpy.random.randint(min(2000, n_draws), n_draws)
if has_chain:
for idx, cn in enumerate(clean_names):
sample[idx] = float(posterior[cn].isel(chain=ichain, draw=idraw))
else:
for idx, cn in enumerate(clean_names):
sample[idx] = float(posterior[cn].isel(draw=idraw))
plot_plate_faces(
"plate_inverted", forward, sample,
axes[0,0], axes[0,1], axes[1,0], color="red", label="Posterior sample", linestyle="dotted"
)
Watermark#
watermark_list = ["cofi", "numpy", "scipy", "matplotlib"]
for pkg in watermark_list:
pkg_var = __import__(pkg)
print(pkg, getattr(pkg_var, "__version__"))
cofi 0.2.11+71.gb28b5b0
numpy 2.2.6
scipy 1.17.1
matplotlib 3.10.9
sphinx_gallery_thumbnail_number = -1
Total running time of the script: (0 minutes 44.679 seconds)