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
# !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_log_prior(my_log_prior)
my_problem.set_log_likelihood(my_log_likelihood)
my_problem.set_model_shape(len(init_param_value))
Define CoFI options
nwalkers = 12
ndim = len(init_param_value)
nsteps = 5000
walkers_start = init_param_value + 0.5 * numpy.random.randn(nwalkers, ndim)
numpy.array([walkers_start[0]])
array([[ 45.24835708, 89.93086785, 160.32384427, 35.76151493,
79.88292331]])
inv_options = cofi.InversionOptions()
inv_options.set_tool("emcee")
inv_options.set_params(nwalkers=nwalkers, nsteps=nsteps, initial_state=walkers_start, progress=True)
######## Run it
inv = cofi.Inversion(my_problem, inv_options)
inv_result = inv.run()
######## Check result
print(f"The inversion result from `emcee`:")
inv_result.summary()
0%| | 0/5000 [00:00<?, ?it/s]
1%| | 52/5000 [00:00<00:09, 508.24it/s]
2%|▏ | 105/5000 [00:00<00:09, 519.80it/s]
3%|▎ | 157/5000 [00:00<00:09, 516.37it/s]
4%|▍ | 213/5000 [00:00<00:09, 530.95it/s]
5%|▌ | 273/5000 [00:00<00:08, 553.50it/s]
7%|▋ | 329/5000 [00:00<00:08, 542.74it/s]
8%|▊ | 384/5000 [00:00<00:08, 530.67it/s]
9%|▉ | 438/5000 [00:00<00:08, 531.68it/s]
10%|▉ | 492/5000 [00:00<00:08, 521.83it/s]
11%|█ | 545/5000 [00:01<00:08, 510.66it/s]
12%|█▏ | 597/5000 [00:01<00:08, 510.57it/s]
13%|█▎ | 649/5000 [00:01<00:08, 508.50it/s]
14%|█▍ | 702/5000 [00:01<00:08, 512.76it/s]
15%|█▌ | 754/5000 [00:01<00:08, 501.87it/s]
16%|█▌ | 805/5000 [00:01<00:08, 493.29it/s]
17%|█▋ | 857/5000 [00:01<00:08, 498.97it/s]
18%|█▊ | 909/5000 [00:01<00:08, 503.22it/s]
19%|█▉ | 961/5000 [00:01<00:07, 507.24it/s]
20%|██ | 1016/5000 [00:01<00:07, 518.19it/s]
21%|██▏ | 1068/5000 [00:02<00:07, 504.89it/s]
22%|██▏ | 1119/5000 [00:02<00:07, 501.97it/s]
23%|██▎ | 1171/5000 [00:02<00:07, 505.29it/s]
24%|██▍ | 1222/5000 [00:02<00:07, 497.47it/s]
25%|██▌ | 1272/5000 [00:02<00:07, 488.70it/s]
26%|██▋ | 1321/5000 [00:02<00:07, 485.88it/s]
27%|██▋ | 1370/5000 [00:02<00:07, 477.79it/s]
28%|██▊ | 1418/5000 [00:02<00:07, 476.51it/s]
29%|██▉ | 1466/5000 [00:02<00:07, 475.49it/s]
30%|███ | 1520/5000 [00:03<00:07, 493.94it/s]
32%|███▏ | 1577/5000 [00:03<00:06, 513.77it/s]
33%|███▎ | 1632/5000 [00:03<00:06, 522.97it/s]
34%|███▎ | 1687/5000 [00:03<00:06, 529.19it/s]
35%|███▍ | 1740/5000 [00:03<00:06, 519.78it/s]
36%|███▌ | 1793/5000 [00:03<00:06, 506.00it/s]
37%|███▋ | 1844/5000 [00:03<00:06, 495.70it/s]
38%|███▊ | 1894/5000 [00:03<00:06, 488.54it/s]
39%|███▉ | 1943/5000 [00:03<00:06, 485.46it/s]
40%|███▉ | 1992/5000 [00:03<00:06, 482.20it/s]
41%|████ | 2041/5000 [00:04<00:06, 479.21it/s]
42%|████▏ | 2090/5000 [00:04<00:06, 481.74it/s]
43%|████▎ | 2139/5000 [00:04<00:05, 481.75it/s]
44%|████▍ | 2188/5000 [00:04<00:05, 482.06it/s]
45%|████▍ | 2237/5000 [00:04<00:05, 480.63it/s]
46%|████▌ | 2286/5000 [00:04<00:05, 477.28it/s]
47%|████▋ | 2334/5000 [00:04<00:05, 476.52it/s]
48%|████▊ | 2382/5000 [00:04<00:05, 476.73it/s]
49%|████▊ | 2430/5000 [00:04<00:05, 476.65it/s]
50%|████▉ | 2478/5000 [00:04<00:05, 466.54it/s]
50%|█████ | 2525/5000 [00:05<00:05, 461.25it/s]
51%|█████▏ | 2573/5000 [00:05<00:05, 465.88it/s]
52%|█████▏ | 2621/5000 [00:05<00:05, 467.91it/s]
53%|█████▎ | 2669/5000 [00:05<00:04, 471.27it/s]
54%|█████▍ | 2717/5000 [00:05<00:04, 470.54it/s]
55%|█████▌ | 2765/5000 [00:05<00:04, 471.93it/s]
56%|█████▋ | 2815/5000 [00:05<00:04, 477.78it/s]
57%|█████▋ | 2867/5000 [00:05<00:04, 487.33it/s]
58%|█████▊ | 2916/5000 [00:05<00:04, 487.20it/s]
59%|█████▉ | 2967/5000 [00:05<00:04, 492.92it/s]
60%|██████ | 3018/5000 [00:06<00:03, 496.88it/s]
61%|██████▏ | 3070/5000 [00:06<00:03, 503.17it/s]
62%|██████▏ | 3121/5000 [00:06<00:03, 499.10it/s]
63%|██████▎ | 3173/5000 [00:06<00:03, 503.65it/s]
64%|██████▍ | 3225/5000 [00:06<00:03, 505.83it/s]
66%|██████▌ | 3277/5000 [00:06<00:03, 509.31it/s]
67%|██████▋ | 3329/5000 [00:06<00:03, 511.73it/s]
68%|██████▊ | 3381/5000 [00:06<00:03, 508.77it/s]
69%|██████▊ | 3433/5000 [00:06<00:03, 510.26it/s]
70%|██████▉ | 3485/5000 [00:07<00:02, 511.26it/s]
71%|███████ | 3537/5000 [00:07<00:02, 507.20it/s]
72%|███████▏ | 3588/5000 [00:07<00:02, 494.40it/s]
73%|███████▎ | 3638/5000 [00:07<00:02, 488.91it/s]
74%|███████▍ | 3690/5000 [00:07<00:02, 496.39it/s]
75%|███████▍ | 3741/5000 [00:07<00:02, 498.35it/s]
76%|███████▌ | 3794/5000 [00:07<00:02, 504.67it/s]
77%|███████▋ | 3847/5000 [00:07<00:02, 509.77it/s]
78%|███████▊ | 3899/5000 [00:07<00:02, 511.82it/s]
79%|███████▉ | 3951/5000 [00:07<00:02, 513.87it/s]
80%|████████ | 4004/5000 [00:08<00:01, 516.80it/s]
81%|████████ | 4056/5000 [00:08<00:01, 511.00it/s]
82%|████████▏ | 4108/5000 [00:08<00:01, 507.88it/s]
83%|████████▎ | 4160/5000 [00:08<00:01, 510.48it/s]
84%|████████▍ | 4214/5000 [00:08<00:01, 518.00it/s]
85%|████████▌ | 4266/5000 [00:08<00:01, 517.33it/s]
86%|████████▋ | 4318/5000 [00:08<00:01, 516.87it/s]
87%|████████▋ | 4373/5000 [00:08<00:01, 523.93it/s]
89%|████████▊ | 4426/5000 [00:08<00:01, 525.09it/s]
90%|████████▉ | 4479/5000 [00:08<00:00, 524.87it/s]
91%|█████████ | 4532/5000 [00:09<00:00, 520.51it/s]
92%|█████████▏| 4585/5000 [00:09<00:00, 519.60it/s]
93%|█████████▎| 4637/5000 [00:09<00:00, 512.62it/s]
94%|█████████▍| 4689/5000 [00:09<00:00, 513.23it/s]
95%|█████████▍| 4741/5000 [00:09<00:00, 510.75it/s]
96%|█████████▌| 4793/5000 [00:09<00:00, 509.30it/s]
97%|█████████▋| 4845/5000 [00:09<00:00, 509.62it/s]
98%|█████████▊| 4896/5000 [00:09<00:00, 508.16it/s]
99%|█████████▉| 4947/5000 [00:09<00:00, 507.80it/s]
100%|█████████▉| 4998/5000 [00:09<00:00, 507.71it/s]
100%|██████████| 5000/5000 [00:09<00:00, 501.47it/s]
The inversion result from `emcee`:
============================
Summary for inversion result
============================
SUCCESS
----------------------------
sampler: <emcee.ensemble.EnsembleSampler object>
blob_names: ['log_likelihood', 'log_prior']
Plotting#
arviz.style.use("arviz-variat")
var_names = [
"Dip (°)",
"Dip azimuth (°)",
"Easting (m)",
"Depth (m)",
"Width (m)",
]
true_values = [60, 65, 175, 30, 90]
sampler = inv_result.sampler
az_idata = inv_result.to_arviz(var_names=var_names)
pc = arviz.plot_trace_dist(
az_idata.sel(draw=slice(2000, None)),
visuals={"xlabel_trace": False, "trace": {"color": "C0", "lw": 0.5}, "dist": {"color": "C0", "lw": 0.5}},
figure_kwargs={"figsize": (12, 20), "constrained_layout": True},
)
var_list = list(az_idata.posterior.data_vars)
for i, vname in enumerate(var_list):
ax_kde = pc.iget_target(i, 0)
ax_trace = pc.iget_target(i, 1)
ax_kde.set_title(vname)
ax_trace.set_title(vname)
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
import arviz_base
true_values_dict = {
f"{var_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))
var_list = list(posterior.data_vars)
n = len(var_list)
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[var_list[j]].values.flatten()
y = posterior[var_list[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
var_list = list(posterior.data_vars)
n_chains = int(posterior.sizes["chain"])
n_draws = int(posterior.sizes["draw"])
ichain = 0
idraw = min(2500, n_draws - 1)
sample = numpy.zeros(5)
for idx, vn in enumerate(var_list):
sample[idx] = float(posterior[vn].isel(chain=ichain, 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):
ichain = numpy.random.randint(0, n_chains)
idraw = numpy.random.randint(min(2000, n_draws), n_draws)
for idx, vn in enumerate(var_list):
sample[idx] = float(posterior[vn].isel(chain=ichain, 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 26.419 seconds)