Thin plate inversion

Thin plate inversion#

What we do in this notebook#

When modelling the airborne electromagnetic response of subvertical bodies such as a volcanogenic massive sulfide deposit they can be approximated using a thin plate in the halfspace of a layered earth, that is the 3D response of a thin plate is being considered. Here we use CoFI to infer such a thin plate target in the basement of a layered earth given airborne electromagnetic data. This tutorial provides a guided tour of a subset of the material in `cofi-examples/examples/vtem_max <inlab-geo/cofi-examples>`__.

Learning outcomes#

An understaning of the pitfalls around using a numerical gradient
An exposé of CoFI’s ability to combine a forward solver and a range of inference methods
An appreciation of the fact that CoFI only requires limited information about the forward problem

# Environment setup (uncomment code below)
# !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/tutorials/thin_plate_inversion

# 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

Problem description#

Given airborne electromagnetic data, a common goal in mineral exploration is to detect and delineate an economic target, and if this target is in the form of a subvertical body in the basement, its electromagnetic response can be approximated using at conductive thin plate located in the halfspace of a layered earth (e.g. Prikhodko et al. 2019).

In this tutorial we look at the inference of such a thin plate with a conductance of $2 \mathrm{S}$ located in a halfspace with a resistivity of $1000 \mathrm{\Omega m}$ with a $20 \mathrm{m}$ thick covering layer that has a resistivity of $300 \mathrm{\Omega m}$. Specifically we develop a thin plate inversion method using CoFI to solve the inverse problem and P223 (Raiche et. al., 2007) to solve the forward problem.

Successful inversion frequently relies on the objective function being smooth and predictable. For the data being inverted here it is advantageous to convert measurements to scale logarithmically to obtain a smoother and more predictable objective function when compared with using the unscaled data. Similarly, plate orientation angles are converted into radians and the remaining model parameters are log-scaled.

Forward solver#

The forward solver is LeroiAir (Raiche et. al, 2007) and the code has been reorganised so that the response measured by an AEM system here VTEMmax is given by a function that can be called from Python. In LeroiAir plates are discretised into cells, with the accuracy of the forward solver being a function of the chosen cell-size. The forward solver is available in a seperate Python package here.

Jacobian via finite differencing#

Parameter estimation methods often rely on the provision of a Jacobian for efficient optimisation. If an analytical Jacobian is not available one can be computed via finite differencing.

$ f’(x_0) = {f(x_0 +h)-f(x_0):raw-latex:`over `h} $

Care must though be taken when choosing the step size $h$, as a too small step size may result in a Jacobian that is affected by a limited accuracy of a forward solver and a too large step size $h$ might result in a Jacobian that is not representative of the derivatives at location $x_0$. Further to this the gradient of the objective function itself is affected by the noise on the data, thus for noisy data choosing a larger step size when computing the Jacobian can be advisable.

In the following we will be using a relative step size $q$ with $h$ defined as $h=x_0*q$.

VTEMmax AEM system#

Airborne electromagnetic systems can be categorised into either helicopter or fixed wing systems. This tutorial is using a VTEMmax system, a helicopter based system developed and operated by Geotech.

https://geotech.ca/services/electromagnetic/vtem-versatile-time-domain-electromagnetic-system/

Plate parametrisation#

The thin plate is parametrised using the parametrisation introduced in (Hauser et. al. 2016). Compared to the commonly employed parametrisation with a plate reference point on the edge of the plate this parametrisation allows for a thin plate to grow and shrink around a plate refrerence point, without the need to move the reference point. This can be advantageous when there is for example a borehole intersecting a thin plate and we seek to determine the extent of the thin plate.

Implementation details#

The problem setup is imported from forward_lib.py but can be adjusted for other applications. The wrapper is created so that we can declare model parameters which are a subset of all the model parameters required by the forward problem. This allows to, for example, invert only for dip of the thin plate with all the other model parameters assumed to be known.

Problem definition#

For convenience we split the problem definition into three objects that are initialised with defaults for our synthetic example introduced in the problem description section.

System specification - system_spec - contains information about the AEM system such as the transmitter waveform as well as start and end times of gates.
Survey setup - survey_setup - contains information about the survey for example the transmitter and receiver locations.
Problem setup - problem_setup - contains the model and exposes the declared model parameters to CoFI.

survey_setup = {
    "tx": numpy.array([205.]),                  # transmitter easting/x-position
    "ty": numpy.array([100.]),                  # transmitter northing/y-position
    "tz": numpy.array([50.]),                   # transmitter height/z-position
    "tazi": numpy.deg2rad(numpy.array([90.])),  # transmitter azimuth
    "tincl": numpy.deg2rad(numpy.array([6.])),  # transmitter inclination
    "rx": numpy.array([205.]),                  # receiver easting/x-position
    "ry": numpy.array([100.]),                  # receiver northing/y-position
    "rz": numpy.array([50.]),                   # receiver height/z-position
    "trdx": numpy.array([0.]),                  # transmitter receiver separation inline
    "trdy": numpy.array([0.]),                  # transmitter receiver separation crossline
    "trdz": numpy.array([0.]),                  # transmitter receiver separation vertical
}

Inverting for the dip of a thin plate#

While a thin plate can not be recovered from a single fiducial, its dip can be recovered from a carefully positioned fiducial. When setting up an inverse problem it is good practice to initially setup the simplest possible problem and experiment with it to ensure that the forward solver and inference methods work as intended. Before attempting an inversion it also is good practice to verify that the misfit function is sensitive to the parameter of interest and that the gradient of the objective function points in the right direction.

The dip of the thin plate thus becomes our declared model parameter and we create the corresponding forward function and set the true value to $60\degree$

forward = ForwardWrapper(true_model, problem_setup, system_spec, survey_setup, ["pdip"])
true_param_value = numpy.array([60])

['pdip']

Problem setup#

True model#

We first the plot the true model and the location of the VTEMmax measurement, which we will be inverting.

#@title plotting function (hidden)
_, 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"
)
plt.tight_layout()
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")

<matplotlib.legend.Legend object at 0x7efda2586210>

Generate synthetic data#

This tutorial uses a simplified noise model that assumes an absolute noise, that is a standard deviation of $0.05$ for the logarithms of the measured and observed data. In the objective and likelihood funtions the noise model is captured in the data covariance matrix, thus a more sophisticated noise model, for example, one that accounts for the known corrleation between time gates in AEM data could easily be implemented.

# The data
absolute_noise= 0.05

# create data and add a realisation of the noise
data_pred_true = forward(true_param_value)
data_obs = data_pred_true + numpy.random.randn(len(data_pred_true))*absolute_noise

# define data covariance matrix
sigma=absolute_noise
Cdinv=numpy.identity(len(data_obs))*(1.0/(sigma*sigma))

Starting model#

Set an initial guess for the dip of the thin plate.

init_param_value = numpy.array([45])

Define helper functions for CoFI#

Challenge: Choose relative step size for the numerical Jacobian#

To compute the numerical Jacobian we need to choose the step size for the perturbation. Here we use a relative step size and if the chosen step size is too large the gradient of the objective function may miss the global minimum, if it is located in a small basins of attraction and thus an inversion may never converge. If the step size is chosen too small the gradient of the objective function will be affected by the noise on the data and the numerical noise of the forward problem. Thus the smallest step size providing a stable gradient is the ideal step size.

Experiment with relative stape size in the range between :math:`0.01` and :math:`0.5` and upload the figure of the gradient and misfit function as a function of plate dip

Using the template below se the relative step size

my_relative_step = <DEFINE_ME>

def my_objective(model):
    dpred = forward(model)
    residual = dpred - data_obs
    return residual.T @ Cdinv @ residual

def my_gradient(model):
    dpred = forward(model)
    jacobian = forward.jacobian(model, relative_step = my_relative_step)
    residual = dpred - data_obs
    return jacobian.T @ Cdinv @ residual

def my_hessian(model):
    jacobian = forward.jacobian(model, relative_step = my_relative_step)
    return jacobian.T @ Cdinv @ jacobian

class PerIterationCallbackFunction:
    def __init__(self):
        self.x = None
        self.i = 0

    def __call__(self, xk):
        print(f"Iteration #{self.i+1}")
        print(f"  objective value: {my_problem.objective(xk)}")
        self.x = xk
        self.i += 1

# Copy the template above, Replace <DEFINE ME> with your answer

#@title Solution

my_relative_step = 0.1

def my_objective(model):
    dpred = forward(model)
    residual = dpred - data_obs
    return residual.T @ Cdinv @ residual

def my_gradient(model):
    dpred = forward(model)
    jacobian = forward.jacobian(model, relative_step = my_relative_step)
    residual = dpred - data_obs
    return jacobian.T @ Cdinv @ residual

def my_hessian(model):
    jacobian = forward.jacobian(model, relative_step = my_relative_step)
    return jacobian.T @ Cdinv @ jacobian

class PerIterationCallbackFunction:
    def __init__(self):
        self.x = None
        self.i = 0

    def __call__(self, xk):
        print(f"Iteration #{self.i+1}")
        print(f"  objective value: {my_problem.objective(xk)}")
        self.x = xk
        self.i += 1

Plot the misfit and gradient as a function of the plate dip#

#@title plotting function (hidden)

all_models = [numpy.array([pdip]) for pdip in range(40, 120, 5)]
all_misfits = []
all_gradients = []
for model in all_models:
    misfit = my_objective(model)
    gradient = my_gradient(model)
    all_misfits.append(misfit)
    all_gradients.append(gradient)
    print(f"pdip: {model}, data misfit: {misfit}, gradient: {gradient}")


fig, ax1 = plt.subplots()
color = 'tab:red'
ax1.plot(all_models, all_misfits,color=color)
ax1.tick_params(axis='y', labelcolor=color)
ax1.set_xlabel("Plate dip")
ax1.set_ylabel("Data misfit",color=color)
plt.grid(axis='x')
ax2 = ax1.twinx()
color = 'tab:blue'
ax2.plot(all_models, all_gradients,color=color)
ax2.set_ylabel('Gradient', color=color)
ax2.tick_params(axis='y', labelcolor=color)
plt.grid(axis='y')
fig.tight_layout()  # otherwise the right y-label is slightly clipped
plt.show()

pdip: [40], data misfit: 231.708835656254, gradient: [-8.38525334]
pdip: [45], data misfit: 161.52059630302878, gradient: [-5.17711959]
pdip: [50], data misfit: 115.89689554155724, gradient: [-3.41393833]
pdip: [55], data misfit: 91.08438559558219, gradient: [-1.79238028]
pdip: [60], data misfit: 78.13208270691166, gradient: [-0.40612889]
pdip: [65], data misfit: 79.47542276986738, gradient: [0.60701868]
pdip: [70], data misfit: 86.57794284916837, gradient: [0.77951888]
pdip: [75], data misfit: 99.82377503892143, gradient: [1.04776418]
pdip: [80], data misfit: 114.82121927725917, gradient: [1.88028555]
pdip: [85], data misfit: 135.80364266137616, gradient: [1.99672184]
pdip: [90], data misfit: 158.91312196239573, gradient: [2.12201565]
pdip: [95], data misfit: 190.36506697630017, gradient: [1.90338506]
pdip: [100], data misfit: 211.34574782298415, gradient: [3.0843773]
pdip: [105], data misfit: 248.95857134110955, gradient: [3.54243063]
pdip: [110], data misfit: 283.9454162116735, gradient: [3.73187192]
pdip: [115], data misfit: 338.73968428973234, gradient: [3.35857177]

Parameter estimation#

First we solve the inverse problem using optimisation, that is we seek to find the minimum of the objective function given as

\[\chi^2 = (\mathbf{d} - \mathbf{f}(\mathbf{m}))^T\mathbf{C}_d^{-1}(\mathbf{d}-\mathbf{f}(\mathbf{m})),\]

where $\mathbf{d}$ is the observed data, $\mathbf{f}(\mathbf{m})$ the model prediction and $\mathbf{C}_d$ the data covariance matrix. The full Newton step is then given as

\[\begin{equation} \Delta \mathbf{m}= (\underbrace{\mathbf{J}^T \mathbf{C}_d^{-1} \mathbf{J}}_{\mathbf{Hessian}})^{-1} (\underbrace{ \mathbf{J}^T\mathbf{C}_d^{-1} (\mathbf{y}-\mathbf{f}(\mathbf{m}))}_\mathbf{Gradient}). \end{equation}\]

The Jacobian $\mathbf{J}$ and Hessian, here approximated by $\mathbf{J}^T \mathbf{C}_d^{-1} \mathbf{J}$, are only local measures of the first and second derivatives of the objective function and given this a non-linear inverse problem and the numerical derivatives can be affected by noise, we can seldom take the full Newton step to compute a model update as we are likely to overshoot and not improve fit to the data.

One strategy is to employ a line search to determine the optimal step length, that means the descent direction is given by the full Newton step with the length adjusted so that it does not overshoot and results in an improvement of the fit to the data. The major alternative to employing a line search is to use a trust region method. Trust regions methods try to estimate the region around the current model within which the assumption of local linearity holds and then limit the model update to stay within that region.

Plotting#

Data#

In the following we plot the vertical and inline component for the true model, the starting model and the MAP model that is the maximum a posterior model, the solution found by the chosen optimisation method.

#@title plotting function (hidden)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))
plot_transient(true_param_value, forward, "Data from true model", ax1, ax2, color="purple")
plot_transient(init_param_value, forward, "Data from starting model", ax1, ax2, color="green", linestyle=":")
plot_transient(my_result.model, forward, "Data from MAP model", ax1, ax2, color="red", linestyle="-.")
ax1.legend(loc="lower center", bbox_to_anchor=(0.5, 1.02), ncol=1, fontsize="small")
ax2.legend(loc="lower center", bbox_to_anchor=(0.5, 1.02), ncol=1, fontsize="small")
ax1.set_title("vertical", pad=55)
ax2.set_title("inline", pad=55)
plt.tight_layout()

Model#

#@title plotting function (hidden)

_, 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"
)
plot_plate_faces(
    "plate_inverted", forward, my_result.model,
    axes[0,0], axes[0,1], axes[1,0], color="red", label="MAP solution", linestyle="dotted"
)
plt.tight_layout()
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")

<matplotlib.legend.Legend object at 0x7efda25860d0>

Ensemble method#

Parameter estimation methods require an objective function while ensemble methods require a likelihood function, typically given in the form of a log likelihood function. The objective function used for the parameter estimation consists of only a data misfit term and thus is closely related to the likelihood function, with the log likelihood being proportional to the value of the objective function multiplied by a factor of $-\frac{1}{2}$

\[p({\mathbf d} | {\mathbf m}) \propto \exp \left\{- \frac{1}{2} ({\mathbf d}-{\mathbf f}({\mathbf m}))^T C_d^{-1} ({\mathbf d}-{\mathbf f}({\mathbf m})) \right\}\]

To use an ensemble method in CoFI we will need to define a likelihood function and prior distribution. The prior distribution we choose here is a uniform distribution with a lower boundary of $10 \degree$ and an upper boundary of $80 \degree$.

m_min=numpy.array([10])
m_max=numpy.array([80])

def my_log_prior(m):    # uniform distribution
    for i in range(len(m)):
        if m[i] < m_min[i] or m[i] > m_max[i]: return -numpy.inf
    return 0.0 # model lies within bounds -> return log(1)

Challenge: Given the objective function define a log likelihood function.#

In the previous section we defined the following objective function.

def my_objective(model):
    dpred = forward(model)
    residual = dpred - data_obs
    return residual.T @ Cdinv @ residual

This function can be used to now create the log likelihood function, typically needed by the ensemble methods made available in CoFI.

def my_log_likelihood(model):
    return # DEFINE ME

#@title Solution
def my_log_likelihood(model):
    return -0.5 * my_objective(model)

Augment the CoFI problem#

To finally be able to use an ensemble method we also need to augment our CoFI problem with the functions we defined for the log of the prior probability and the log of the likelihood function. We will again be using emcee, which we already used in the linear regression tutorial.

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 = 2
ndim = len(init_param_value)
nsteps = 50
walkers_start = init_param_value + 1 * numpy.random.randn(nwalkers, ndim)

CoFI Inversion#

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)
my_result = inv.run()

######## Check result
print(f"The inversion result from `emcee`:")
my_result.summary()

  0%|          | 0/50 [00:00<?, ?it/s]
  2%|▏         | 1/50 [00:00<00:27,  1.79it/s]
  4%|▍         | 2/50 [00:01<00:26,  1.79it/s]
  6%|▌         | 3/50 [00:01<00:26,  1.79it/s]
  8%|▊         | 4/50 [00:02<00:25,  1.79it/s]
 10%|█         | 5/50 [00:02<00:25,  1.79it/s]
 12%|█▏        | 6/50 [00:03<00:24,  1.78it/s]
 14%|█▍        | 7/50 [00:03<00:24,  1.77it/s]
 16%|█▌        | 8/50 [00:04<00:23,  1.77it/s]
 18%|█▊        | 9/50 [00:05<00:23,  1.78it/s]
 20%|██        | 10/50 [00:05<00:22,  1.78it/s]
 22%|██▏       | 11/50 [00:06<00:21,  1.78it/s]
 24%|██▍       | 12/50 [00:06<00:21,  1.78it/s]
 26%|██▌       | 13/50 [00:07<00:20,  1.78it/s]
 28%|██▊       | 14/50 [00:07<00:20,  1.78it/s]
 30%|███       | 15/50 [00:08<00:19,  1.78it/s]
 32%|███▏      | 16/50 [00:09<00:19,  1.77it/s]
 34%|███▍      | 17/50 [00:09<00:18,  1.77it/s]
 36%|███▌      | 18/50 [00:10<00:18,  1.78it/s]
 38%|███▊      | 19/50 [00:10<00:17,  1.78it/s]
 40%|████      | 20/50 [00:11<00:16,  1.78it/s]
 42%|████▏     | 21/50 [00:11<00:16,  1.78it/s]
 44%|████▍     | 22/50 [00:12<00:15,  1.78it/s]
 46%|████▌     | 23/50 [00:12<00:15,  1.78it/s]
 48%|████▊     | 24/50 [00:13<00:14,  1.77it/s]
 50%|█████     | 25/50 [00:14<00:14,  1.77it/s]
 52%|█████▏    | 26/50 [00:14<00:13,  1.77it/s]
 54%|█████▍    | 27/50 [00:15<00:12,  1.78it/s]
 56%|█████▌    | 28/50 [00:15<00:12,  1.78it/s]
 58%|█████▊    | 29/50 [00:16<00:11,  1.78it/s]
 60%|██████    | 30/50 [00:16<00:11,  1.78it/s]
 62%|██████▏   | 31/50 [00:17<00:10,  1.78it/s]
 64%|██████▍   | 32/50 [00:17<00:10,  1.78it/s]
 66%|██████▌   | 33/50 [00:18<00:09,  1.78it/s]
 68%|██████▊   | 34/50 [00:19<00:09,  1.77it/s]
 70%|███████   | 35/50 [00:19<00:08,  1.78it/s]
 72%|███████▏  | 36/50 [00:20<00:07,  1.78it/s]
 74%|███████▍  | 37/50 [00:20<00:07,  1.78it/s]
 76%|███████▌  | 38/50 [00:21<00:06,  1.78it/s]
 78%|███████▊  | 39/50 [00:21<00:06,  1.79it/s]
 80%|████████  | 40/50 [00:22<00:05,  1.78it/s]
 82%|████████▏ | 41/50 [00:23<00:05,  1.78it/s]
 84%|████████▍ | 42/50 [00:23<00:04,  1.78it/s]
 86%|████████▌ | 43/50 [00:24<00:03,  1.78it/s]
 88%|████████▊ | 44/50 [00:24<00:03,  1.78it/s]
 90%|█████████ | 45/50 [00:25<00:02,  1.79it/s]
 92%|█████████▏| 46/50 [00:25<00:02,  1.80it/s]
 94%|█████████▍| 47/50 [00:26<00:01,  1.80it/s]
 96%|█████████▌| 48/50 [00:26<00:01,  1.80it/s]
 98%|█████████▊| 49/50 [00:27<00:00,  1.78it/s]
100%|██████████| 50/50 [00:28<00:00,  1.78it/s]
100%|██████████| 50/50 [00:28<00:00,  1.78it/s]
The inversion result from `emcee`:
============================
Summary for inversion result
============================
SUCCESS
----------------------------
sampler: <emcee.ensemble.EnsembleSampler object>
blob_names: ['log_likelihood', 'log_prior']

#@title plotting function (hidden)
sampler = my_result.sampler
arviz.style.use("arviz-variat")
var_names = [
    "plate dip (°)",
]
az_idata = my_result.to_arviz(var_names=var_names)
true_values = [60]
pc = arviz.plot_trace_dist(
    az_idata.sel(draw=slice(0, None)),
    visuals={"xlabel_trace": False, "trace": {"color": "C0"}, "dist": {"color": "C0"}},
    figure_kwargs={"figsize": (12, 4), "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)

While both chains have moved away from the starting model in the direction of the true dip the 50 samples are not enough to recover the posterior distribution. In practice a longer chain would be needed, but this would increase the number of times we need to call the forward problem which is computationally expensive.

Limited computational resources#

From our tests we know that the objective function appears smooth and well behaved with a large basin of attraction, thus we can create a surrogate model. The basic idea is that given samples of the objective function as a function of the declared model parameters we can interpolate a response surface and use it instead of solving the computationally more expensive forward problem when we need to evaluate the log likelihood function. Evaluating the surrogate model will take a fraction of the time. In the following we thus create a surrogate model for the objective function using the surrogate modelling toolbox.

We first generate random samples of the objective function using latin hypercube sampling.

#@title random samples of the objective function in the range 10 to 90 for the dip.
ndim=len(true_param_value)
xlimits=numpy.array([[10,90]])
ntrain=25
ntest=5
sampling = smt.sampling_methods.LHS(xlimits=xlimits,seed=42)
xtrain=sampling(ntrain)
ytrain=[]
xtest=sampling(ntest)
ytest=[]
for x in tqdm.tqdm(xtrain):
    ytrain.append(my_objective(x))
for x in tqdm.tqdm(xtest):
    ytest.append(my_objective(x))
ytrain=numpy.array(ytrain)
ytest=numpy.array(ytest)

  0%|          | 0/25 [00:00<?, ?it/s]
  4%|▍         | 1/25 [00:00<00:06,  3.46it/s]
  8%|▊         | 2/25 [00:00<00:06,  3.53it/s]
 12%|█▏        | 3/25 [00:00<00:06,  3.57it/s]
 16%|█▌        | 4/25 [00:01<00:05,  3.55it/s]
 20%|██        | 5/25 [00:01<00:05,  3.56it/s]
 24%|██▍       | 6/25 [00:01<00:05,  3.58it/s]
 28%|██▊       | 7/25 [00:01<00:04,  3.61it/s]
 32%|███▏      | 8/25 [00:02<00:04,  3.60it/s]
 36%|███▌      | 9/25 [00:02<00:04,  3.58it/s]
 40%|████      | 10/25 [00:02<00:04,  3.56it/s]
 44%|████▍     | 11/25 [00:03<00:03,  3.58it/s]
 48%|████▊     | 12/25 [00:03<00:03,  3.56it/s]
 52%|█████▏    | 13/25 [00:03<00:03,  3.54it/s]
 56%|█████▌    | 14/25 [00:03<00:03,  3.57it/s]
 60%|██████    | 15/25 [00:04<00:02,  3.57it/s]
 64%|██████▍   | 16/25 [00:04<00:02,  3.59it/s]
 68%|██████▊   | 17/25 [00:04<00:02,  3.57it/s]
 72%|███████▏  | 18/25 [00:05<00:01,  3.60it/s]
 76%|███████▌  | 19/25 [00:05<00:01,  3.60it/s]
 80%|████████  | 20/25 [00:05<00:01,  3.59it/s]
 84%|████████▍ | 21/25 [00:05<00:01,  3.59it/s]
 88%|████████▊ | 22/25 [00:06<00:00,  3.52it/s]
 92%|█████████▏| 23/25 [00:06<00:00,  3.50it/s]
 96%|█████████▌| 24/25 [00:06<00:00,  3.49it/s]
100%|██████████| 25/25 [00:07<00:00,  3.53it/s]
100%|██████████| 25/25 [00:07<00:00,  3.56it/s]

  0%|          | 0/5 [00:00<?, ?it/s]
 20%|██        | 1/5 [00:00<00:01,  3.62it/s]
 40%|████      | 2/5 [00:00<00:00,  3.53it/s]
 60%|██████    | 3/5 [00:00<00:00,  3.55it/s]
 80%|████████  | 4/5 [00:01<00:00,  3.59it/s]
100%|██████████| 5/5 [00:01<00:00,  3.61it/s]
100%|██████████| 5/5 [00:01<00:00,  3.59it/s]

Next we create the surrogate model itself, that is here the construction of a response surface using Kriging.

sm = smt.surrogate_models.KRG(theta0=[1e-2]*ndim,print_prediction = False)
sm.set_training_values(xtrain,ytrain)
sm.train()

___________________________________________________________________________

                                  Kriging
___________________________________________________________________________

 Problem size

      # training points.        : 25

___________________________________________________________________________

 Training

   Training ...
   Training - done. Time (sec):  0.5426669

Finally, we test how well our surrogate model predicts the value of the objective function.

#@title prediction of the validation points
y = sm.predict_values(xtest)
# Estimated variance for the validation points
s2 = sm.predict_variances(xtest)
#plot with the associated interval confidence
yerr= 2*3*numpy.sqrt(s2) #in order to use +/- 3 x standard deviation: 99% confidence interval estimation

# Plot the function, the prediction and the 99% confidence interval based on
# the MSE
fig = plt.figure()
plt.plot(ytest, ytest, '-', label='$y_{true}$')
plt.plot(ytest, y, 'r.', label='$\\hat{y}$')
plt.errorbar(numpy.squeeze(ytest), numpy.squeeze(y), yerr=numpy.squeeze(yerr), fmt = 'none', capsize = 5, ecolor = 'lightgray', elinewidth = 1, capthick = 0.5, label='confidence estimate 99%')
plt.xlabel('$y_{true}$')
plt.ylabel('$\\hat{y}$')

plt.legend(loc='upper left')
plt.title('Validation of the model prediction with confidence estimates')
plt.show()

Now we run emcee again, but we call the surrogate model instead of the P223 kernel and thus we can generate more samples of the posterior distribution in a fraction of the time. It is however important to keep in mind that the surrogate model is only an approximation and likely to have a limited accuracy particularly if we increase the number of declared model parameters.

nwalkers = 5
ndim = len(init_param_value)
nsteps = 500
walkers_start = init_param_value + 1 * numpy.random.randn(nwalkers, ndim)

def my_objective(model):
    val=sm.predict_values(numpy.array([model]))[0][0]
    if val<1e-3:
        return 1e-3
    else:
        return val

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)
my_result = inv.run()

######## Check result
print(f"The inversion result from `emcee`:")
my_result.summary()

  0%|          | 0/500 [00:00<?, ?it/s]
 19%|█▉        | 97/500 [00:00<00:00, 963.57it/s]
 39%|███▉      | 194/500 [00:00<00:00, 955.00it/s]
 58%|█████▊    | 290/500 [00:00<00:00, 956.63it/s]
 77%|███████▋  | 386/500 [00:00<00:00, 954.24it/s]
 96%|█████████▋| 482/500 [00:00<00:00, 946.63it/s]
100%|██████████| 500/500 [00:00<00:00, 949.63it/s]
The inversion result from `emcee`:
============================
Summary for inversion result
============================
SUCCESS
----------------------------
sampler: <emcee.ensemble.EnsembleSampler object>
blob_names: ['log_likelihood', 'log_prior']

#@title plotting function (hidden)
sampler = my_result.sampler
arviz.style.use("arviz-variat")
var_names = [
    "plate dip (°)",
]
az_idata = my_result.to_arviz(var_names=var_names)
true_values = [60]
pc = arviz.plot_trace_dist(
    az_idata.sel(draw=slice(100, None)),
    visuals={"xlabel_trace": False, "trace": {"color": "C0"}, "dist": {"color": "C0"}},
    figure_kwargs={"figsize": (12, 4), "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)

Inverting for a thin plate given three survey lines#

A more realistic synthetic example is the inference of a thin plate target given three survey lines. It now becomes possible to invert for the easting, depth of the plate reference point, the plate dip and plate azimuth and the plate length.

Problem setup#

In the following we define three survey lines covering the thin plate

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#

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])
}

We now increase the number of declared model parameters and they include the plate dip, the plate dip azimuth, the easting of the plate reference point, the depth of the plate reference point and the plate width and will only be using the vertical component.

As a general rule we can only constrain parameters if there are fiducials in the survey that are sensitive to them and also fiducials that are not sensitive to them. In order words the anomaly needs to be closed with respect to the model parameter in question; to for example constrain plate length we would need survey lines to the north and south that are not overflying the thin plate, as this is not the case in our setup we do not invert for plate length.

forward = ForwardWrapper(true_model, problem_setup, system_spec, survey_setup,
                         ["pdip","pdzm", "peast", "pwdth2","ptop"], data_returned=["vertical"]);

['pdip', 'pdzm', 'peast', 'pwdth2', 'ptop']

# 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])

#@title plotting function (hidden)

_, 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"
)
plt.tight_layout()
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")

<matplotlib.legend.Legend object at 0x7efda24a5450>

Generate synthetic data#

We again generate a synthetic data set and add a realisation of the noise.

# The data
absolute_noise= 0.05

# create data and add a realisation of the noise
data_pred_true = forward(true_param_value)
data_obs = data_pred_true + numpy.random.randn(len(data_pred_true))*absolute_noise

# define data covariance matrix
sigma=absolute_noise
Cdinv=numpy.identity(len(data_obs))*(1.0/(sigma*sigma))

Challenge: Implement a parameter estimation or ensemble method in CoFI#

CoFI is about experimentation and given the experiments in the previous section you can now head down a branch of the CoFI tree and infer a thin plate from the synthetic data we just generated. The first choice is between parameter estimation and ensemble methods.

Parameter estimation
Ensemble methods

Colab: If this notebook is used on colab anchor links will currently not work and the table of contents needs to be used to navigate to the relevant section… #3983

The parameter estimation methods call the forward kernel and compute a numerical Jacobian, while the ensemble methods will in the following again use a surrogate model to compute the objective/likelihood function, which takes a fraction of a second. Thus the suggestion is to use an ensemble method.

Once you have performed an inversion using CoFI upload your solution

Parameter estimation applied to three survey lines#

Initialise a model for inversion

forward.params_to_invert

['pdip', 'pdzm', 'peast', 'ptop', 'pwdth2']

init_param_value = numpy.array([45, 90, 150, 20, 80])

Define helper functions for CoFI

def my_objective(model):
    dpred = forward(model)
    residual = dpred - data_obs
    return residual.T @ Cdinv @ residual

def my_gradient(model):
    dpred = forward(model)
    jacobian = forward.jacobian(model, relative_step=0.1)
    residual = dpred - data_obs
    return jacobian.T @ Cdinv @ residual

def my_hessian(model):
    jacobian = forward.jacobian(model)
    return jacobian.T @ Cdinv @ jacobian

class PerIterationCallbackFunction:
    def __init__(self):
        self.x = None
        self.i = 0

    def __call__(self, xk):
        print(f"Iteration #{self.i+1}")
        print(f"  objective value: {my_problem.objective(xk)}")
        self.x = xk
        self.i += 1

Define CoFI problem

my_problem = cofi.BaseProblem()
my_problem.set_objective(my_objective)
my_problem.set_gradient(my_gradient)
my_problem.set_hessian(my_hessian)
my_problem.set_initial_model(init_param_value)

Challenge: Choose a parameter estimation method#

CoFI provides access to the parameter estimation methods that are available in `scipy.optimize.minimize <https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html>`__

For practical application we are interested in a solver that converges with the fewest calls to the forward problem to a model that is acceptably close to the true model and explains the data. The consequence of employing a line search or trust region method or more broadly any method seeking to find the optimal step length is, that typically additional calls to a forward problem need to be made to determine the optimal step length and different approaches require different numbers of calls to the forward problem depending on the shape of the objective function.

Colab: The computational resources offered by colab for free are limited and thus the inverisons performed in the following may take a while if the free resources offered by colab are being used…

Choose one of the following two solvers and perform a parameter estimation using CoFI and upload your solution - newton-cg - https://docs.scipy.org/doc/scipy/reference/optimize.minimize-newtoncg.html - trust-ncg- https://docs.scipy.org/doc/scipy/reference/optimize.minimize-trustncg.html

Define CoFI options

Use the template below and set the CoFI tool to newton-cg or trust-ncg

my_options = cofi.InversionOptions()
my_options.set_tool("scipy.optimize.minimize")
my_options.set_params(method=<DEFINE_ME>,callback=PerIterationCallbackFunction(),options={"maxiter": 10})

# Copy the template above, Replace <DEFINE ME> with your answer

#@title Solution newton-cg
my_options = cofi.InversionOptions()
my_options.set_tool("scipy.optimize.minimize")
my_options.set_params(method="newton-cg",callback=PerIterationCallbackFunction(),options={"maxiter": 10})

#@title Solution trust-ncg
my_options = cofi.InversionOptions()
my_options.set_tool("scipy.optimize.minimize")
my_options.set_params(method="trust-ncg",callback=PerIterationCallbackFunction(),options={"maxiter": 10})

Run CoFI inversion

my_inversion = cofi.Inversion(my_problem, my_options)
my_result = my_inversion.run()

Iteration #1
  objective value: 81782.04925107711
Iteration #2
  objective value: 69010.63281493413
Iteration #3
  objective value: 49533.378739060325
Iteration #4
  objective value: 27942.1106349338
Iteration #5
  objective value: 13518.328094206927
Iteration #6
  objective value: 8809.363813455235
Iteration #7
  objective value: 3944.3620556606566
Iteration #8
  objective value: 2342.7727354235585
Iteration #9
  objective value: 1997.313570697696
Iteration #10
  objective value: 1698.473407253146

Plotting#

Data - Profiles#

#@title plotting function (hidden)

# Select gates to plot
idx_to_plot = numpy.arange(8, 30)

_, axes = plt.subplots(3, 1, figsize=(12,12))

for i in range(3):
    plot_predicted_profile(true_param_value, forward, "Data from true model", gate_idx=idx_to_plot,
                                        line_id=[i], ax=axes[i], color="purple")
    plot_predicted_profile(init_param_value, forward, "Data from starting model", gate_idx=idx_to_plot,
                                        line_id=[i], ax=axes[i], color="green", linestyle=":")
    plot_predicted_profile(my_result.model, forward, "Data from MAP model", gate_idx=idx_to_plot,
                                        line_id=[i], ax=axes[i], color="red", linestyle="-.")
    axes[i].set_title("Crossline {} (m)".format(ty_points[i]))
    axes[i].legend(bbox_to_anchor=(1.04, 0), loc="lower left")
plt.tight_layout()

Model#

#@title plotting function (hidden)

_, 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"
)
plot_plate_faces(
    "plate_inverted", forward, my_result.model,
    axes[0,0], axes[0,1], axes[1,0], color="red", label="MAP model", linestyle="dotted"
)

plt.tight_layout()

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

<matplotlib.legend.Legend object at 0x7efd9290ac10>

Ensemble methods applied to three survey lines#

To speed up the forward computations for this tutorial and to be able to use ensemble methods we again create a surrogate model for the objective function using the surrogate modelling toolbox.

The two steps to create the surrogate model are the sampling of the objective function using latin hypercube sampling and the creation of the surrogate model itself, that is here the construction of a response surface using Kriging. For completness the two notebooks implementing this are given here: - Latin Hypercube Sampling - Surrogate model creation

#@title Create surrogate model given latin hypercube samples

# Different versions of the surrogate modelling toolbox store the model in different
# formats thus it is safer to just create the model on the fly given the latin hypercube
# samples which here are pre-computed.

with open('three_survey_lines_lhs.npy', 'rb') as f:
    ndim=int(numpy.load(f))
    xlimits=numpy.load(f)
    xtrain=numpy.load(f)
    ytrain=numpy.load(f)
    xtest=numpy.load(f)
    ytest=numpy.load(f)

xlimits=xlimits.astype('double')
xtrain=xtrain[0:100].astype('double')
ytrain=ytrain[0:100].astype('double')
xtest=xtest[0:10].astype('double')
ytest=ytest[0:10].astype('double')

sm = smt.surrogate_models.KRG(theta0=[1e-2]*ndim,print_prediction = False)
sm.set_training_values(xtrain,ytrain)
sm.train()

___________________________________________________________________________

                                  Kriging
___________________________________________________________________________

 Problem size

      # training points.        : 100

___________________________________________________________________________

 Training

   Training ...
   Training - done. Time (sec):  4.5371859

Initialise a model for inversion and define prior distribution

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)

Challenge: Select an ensemble method#

The nature of the airborne electromagnetic forward physics is that the deeper a feature of interest the less well it can be recovered. This information is not captured in the solution obtained using a parameter estimation method. CoFI provides access to range of ensemble methods that recover the distribution of models that fit the data and thus allows to estimate model uncertainty.

Using one of the following three ensemble methods available in CofI complete the relevant section and upload your solution.

`emcee Affine Invariant Markov chain Monte Carlo Ensemble sampler <#Affine-Invariant-Markov-chain-Monte-Carlo-Ensemble-sampler>`__
`neighpy Neighbourhood algorithm <#Neighbourhood-algorithm>`__
`bayesbay Metropolis Hastings algorithm <#Metropolis-Hastings-algorithm>`__

Colab: If this notebook is used on colab anchor links will currently not work and the table of contents needs to be used to navigate to the relevant section… #3983

Affine Invariant Markov chain Monte Carlo Ensemble sampler#

Using emcee as the string in inv_options.set_tool() CoFI offers an interface to the Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler by Goodman and Weare 2010 to sample the posterior distribution. (See more details about emcee).

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

#@title defaults for emcee
nwalkers = 12
ndim = len(init_param_value)
nsteps = 5000
walkers_start = init_param_value + 0.5 * numpy.random.randn(nwalkers, ndim)

Use the template below and set the CoFI tool to emcee

inv_options = cofi.InversionOptions()
inv_options.set_tool(<DEFINE_ME>)
inv_options.set_params(nwalkers=nwalkers, nsteps=nsteps, initial_state=walkers_start, progress=True)

######## Run it
inv = cofi.Inversion(my_problem, inv_options)
my_result = inv.run()

######## Check result
print(f"The inversion result from `emcee`:")
my_result.summary()

# Copy the template above, Replace <DEFINE ME> with your answer

#@title Solution

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)
my_result = inv.run()

######## Check result
print(f"The inversion result from `emcee`:")
my_result.summary()

  0%|          | 0/5000 [00:00<?, ?it/s]
  1%|          | 33/5000 [00:00<00:15, 325.08it/s]
  2%|▏         | 85/5000 [00:00<00:11, 437.69it/s]
  3%|▎         | 132/5000 [00:00<00:10, 451.25it/s]
  4%|▎         | 181/5000 [00:00<00:10, 464.47it/s]
  5%|▍         | 229/5000 [00:00<00:10, 467.58it/s]
  6%|▌         | 277/5000 [00:00<00:10, 470.46it/s]
  7%|▋         | 326/5000 [00:00<00:09, 473.97it/s]
  8%|▊         | 375/5000 [00:00<00:09, 477.36it/s]
  8%|▊         | 424/5000 [00:00<00:09, 478.96it/s]
  9%|▉         | 473/5000 [00:01<00:09, 479.88it/s]
 10%|█         | 522/5000 [00:01<00:09, 480.14it/s]
 11%|█▏        | 571/5000 [00:01<00:09, 479.21it/s]
 12%|█▏        | 619/5000 [00:01<00:09, 479.41it/s]
 13%|█▎        | 667/5000 [00:01<00:09, 477.21it/s]
 14%|█▍        | 715/5000 [00:01<00:08, 476.38it/s]
 15%|█▌        | 763/5000 [00:01<00:08, 475.84it/s]
 16%|█▌        | 811/5000 [00:01<00:08, 476.55it/s]
 17%|█▋        | 859/5000 [00:01<00:08, 475.28it/s]
 18%|█▊        | 907/5000 [00:01<00:08, 475.35it/s]
 19%|█▉        | 955/5000 [00:02<00:08, 473.51it/s]
 20%|██        | 1004/5000 [00:02<00:08, 476.53it/s]
 21%|██        | 1056/5000 [00:02<00:08, 486.55it/s]
 22%|██▏       | 1106/5000 [00:02<00:07, 489.39it/s]
 23%|██▎       | 1157/5000 [00:02<00:07, 492.80it/s]
 24%|██▍       | 1209/5000 [00:02<00:07, 498.63it/s]
 25%|██▌       | 1263/5000 [00:02<00:07, 509.26it/s]
 26%|██▋       | 1315/5000 [00:02<00:07, 510.82it/s]
 27%|██▋       | 1367/5000 [00:02<00:07, 509.74it/s]
 28%|██▊       | 1419/5000 [00:02<00:06, 512.75it/s]
 29%|██▉       | 1471/5000 [00:03<00:06, 512.45it/s]
 30%|███       | 1524/5000 [00:03<00:06, 515.85it/s]
 32%|███▏      | 1576/5000 [00:03<00:06, 514.87it/s]
 33%|███▎      | 1628/5000 [00:03<00:06, 510.77it/s]
 34%|███▎      | 1681/5000 [00:03<00:06, 514.25it/s]
 35%|███▍      | 1735/5000 [00:03<00:06, 519.69it/s]
 36%|███▌      | 1787/5000 [00:03<00:06, 518.64it/s]
 37%|███▋      | 1839/5000 [00:03<00:06, 515.81it/s]
 38%|███▊      | 1891/5000 [00:03<00:06, 509.72it/s]
 39%|███▉      | 1942/5000 [00:03<00:06, 508.31it/s]
 40%|███▉      | 1993/5000 [00:04<00:05, 507.50it/s]
 41%|████      | 2046/5000 [00:04<00:05, 511.32it/s]
 42%|████▏     | 2101/5000 [00:04<00:05, 519.26it/s]
 43%|████▎     | 2154/5000 [00:04<00:05, 519.82it/s]
 44%|████▍     | 2206/5000 [00:04<00:05, 514.55it/s]
 45%|████▌     | 2258/5000 [00:04<00:05, 507.43it/s]
 46%|████▌     | 2309/5000 [00:04<00:05, 502.74it/s]
 47%|████▋     | 2360/5000 [00:04<00:05, 500.74it/s]
 48%|████▊     | 2413/5000 [00:04<00:05, 506.98it/s]
 49%|████▉     | 2466/5000 [00:04<00:04, 512.02it/s]
 50%|█████     | 2518/5000 [00:05<00:04, 513.06it/s]
 51%|█████▏    | 2570/5000 [00:05<00:04, 512.89it/s]
 52%|█████▏    | 2622/5000 [00:05<00:04, 511.78it/s]
 54%|█████▎    | 2675/5000 [00:05<00:04, 514.83it/s]
 55%|█████▍    | 2727/5000 [00:05<00:04, 515.39it/s]
 56%|█████▌    | 2779/5000 [00:05<00:04, 515.42it/s]
 57%|█████▋    | 2832/5000 [00:05<00:04, 519.73it/s]
 58%|█████▊    | 2884/5000 [00:05<00:04, 516.19it/s]
 59%|█████▊    | 2936/5000 [00:05<00:03, 516.12it/s]
 60%|█████▉    | 2988/5000 [00:05<00:03, 512.53it/s]
 61%|██████    | 3040/5000 [00:06<00:03, 507.49it/s]
 62%|██████▏   | 3091/5000 [00:06<00:03, 504.41it/s]
 63%|██████▎   | 3144/5000 [00:06<00:03, 510.85it/s]
 64%|██████▍   | 3197/5000 [00:06<00:03, 514.73it/s]
 65%|██████▍   | 3249/5000 [00:06<00:03, 512.01it/s]
 66%|██████▌   | 3302/5000 [00:06<00:03, 515.19it/s]
 67%|██████▋   | 3354/5000 [00:06<00:03, 514.94it/s]
 68%|██████▊   | 3406/5000 [00:06<00:03, 508.19it/s]
 69%|██████▉   | 3457/5000 [00:06<00:03, 506.62it/s]
 70%|███████   | 3510/5000 [00:07<00:02, 511.84it/s]
 71%|███████   | 3562/5000 [00:07<00:02, 508.48it/s]
 72%|███████▏  | 3613/5000 [00:07<00:02, 508.76it/s]
 73%|███████▎  | 3664/5000 [00:07<00:02, 505.83it/s]
 74%|███████▍  | 3716/5000 [00:07<00:02, 509.29it/s]
 75%|███████▌  | 3769/5000 [00:07<00:02, 514.08it/s]
 76%|███████▋  | 3821/5000 [00:07<00:02, 513.72it/s]
 77%|███████▋  | 3874/5000 [00:07<00:02, 516.47it/s]
 79%|███████▊  | 3926/5000 [00:07<00:02, 517.30it/s]
 80%|███████▉  | 3981/5000 [00:07<00:01, 525.55it/s]
 81%|████████  | 4034/5000 [00:08<00:01, 525.35it/s]
 82%|████████▏ | 4087/5000 [00:08<00:01, 517.61it/s]
 83%|████████▎ | 4139/5000 [00:08<00:01, 516.26it/s]
 84%|████████▍ | 4193/5000 [00:08<00:01, 521.22it/s]
 85%|████████▍ | 4246/5000 [00:08<00:01, 522.74it/s]
 86%|████████▌ | 4299/5000 [00:08<00:01, 520.78it/s]
 87%|████████▋ | 4352/5000 [00:08<00:01, 519.42it/s]
 88%|████████▊ | 4405/5000 [00:08<00:01, 521.92it/s]
 89%|████████▉ | 4459/5000 [00:08<00:01, 525.21it/s]
 90%|█████████ | 4514/5000 [00:08<00:00, 530.01it/s]
 91%|█████████▏| 4568/5000 [00:09<00:00, 531.83it/s]
 92%|█████████▏| 4622/5000 [00:09<00:00, 524.44it/s]
 94%|█████████▎| 4675/5000 [00:09<00:00, 523.60it/s]
 95%|█████████▍| 4728/5000 [00:09<00:00, 523.40it/s]
 96%|█████████▌| 4781/5000 [00:09<00:00, 520.48it/s]
 97%|█████████▋| 4834/5000 [00:09<00:00, 520.93it/s]
 98%|█████████▊| 4887/5000 [00:09<00:00, 519.42it/s]
 99%|█████████▉| 4940/5000 [00:09<00:00, 520.45it/s]
100%|█████████▉| 4993/5000 [00:09<00:00, 520.66it/s]
100%|██████████| 5000/5000 [00:09<00:00, 505.90it/s]
The inversion result from `emcee`:
============================
Summary for inversion result
============================
SUCCESS
----------------------------
sampler: <emcee.ensemble.EnsembleSampler object>
blob_names: ['log_likelihood', 'log_prior']

#@title plotting function (hidden)

arviz.style.use("arviz-variat")

var_names = [
    "Dip (°)",
    "Dip azimuth (°)",
    "Easting (m)",
    "Depth (m)",
    "Width (m)",
]

true_values = [60, 65, 175, 30, 90]
sampler = my_result.sampler
az_idata = my_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"}, "dist": {"color": "C0"}},
    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)

#@title plotting function (hidden)

true_values = {
    f"{var_names[i]}": true_param_value[i] for i in range(init_param_value.size)
}
import arviz_base
arviz_base.rcParams["plot.max_subplots"] = 80
pm = arviz.plot_pair(
    az_idata.sel(draw=slice(4000, None)),
    marginal=True,
    triangle="lower",
)
ref_names = list(true_values.keys())
ref_vals = list(true_values.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,
            )

#@title plotting function (hidden)

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:
    sample[0] = float(posterior["Dip (°)"].isel(chain=ichain, draw=idraw))
    sample[1] = float(posterior["Dip azimuth (°)"].isel(chain=ichain, draw=idraw))
    sample[2] = float(posterior["Easting (m)"].isel(chain=ichain, draw=idraw))
    sample[3] = float(posterior["Depth (m)"].isel(chain=ichain, draw=idraw))
    sample[4] = float(posterior["Width (m)"].isel(chain=ichain, draw=idraw))
else:
    sample[0] = float(posterior["Dip (°)"].isel(draw=idraw))
    sample[1] = float(posterior["Dip azimuth (°)"].isel(draw=idraw))
    sample[2] = float(posterior["Easting (m)"].isel(draw=idraw))
    sample[3] = float(posterior["Depth (m)"].isel(draw=idraw))
    sample[4] = float(posterior["Width (m)"].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):
    ichain = numpy.random.randint(0, n_chains)
    idraw = numpy.random.randint(min(2000, n_draws), n_draws)
    if has_chain:
        sample[0] = float(posterior["Dip (°)"].isel(chain=ichain, draw=idraw))
        sample[1] = float(posterior["Dip azimuth (°)"].isel(chain=ichain, draw=idraw))
        sample[2] = float(posterior["Easting (m)"].isel(chain=ichain, draw=idraw))
        sample[3] = float(posterior["Depth (m)"].isel(chain=ichain, draw=idraw))
        sample[4] = float(posterior["Width (m)"].isel(chain=ichain, draw=idraw))
    else:
        sample[0] = float(posterior["Dip (°)"].isel(draw=idraw))
        sample[1] = float(posterior["Dip azimuth (°)"].isel(draw=idraw))
        sample[2] = float(posterior["Easting (m)"].isel(draw=idraw))
        sample[3] = float(posterior["Depth (m)"].isel(draw=idraw))
        sample[4] = float(posterior["Width (m)"].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"
    )

Neighbourhood algorithm#

The Neighbourhood Algorithm is a two-stage ensemble method for non-linear inverse problems with application as a direct search method for global optimisation. The first stage seeks to find points in model space with acceptable values of the objective function. The second stage analysis the ensemble of models generated in the first stage and provides Bayesian measures of properties of the ensemble such as resolution and covariance structure.

Here CoFI is providing an interface the the implementation of the Neighbourhood Algorithm provided by neighpy and it is accessed by using neighpy as the string in inv_options.set_tool()

my_problem = cofi.BaseProblem()
my_problem.set_objective(my_objective)

Using the template below set the CoFI tool to neighpy

inv_options = cofi.InversionOptions()
inv_options.set_tool(<DEFINE_ME>)
inv_options.suggest_solver_params()

# Copy the template above, Replace <DEFINE ME> with your answer

#@title Solution
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)
my_result = inv.run()

######## Check result
print(f"The inversion result from `neighpy`:")
my_result.summary()

NAI - Initial Random Search

NAI - Optimisation Loop:   0%|          | 0/100 [00:00<?, ?it/s]
NAI - Optimisation Loop:   1%|          | 1/100 [00:02<04:44,  2.88s/it]
NAI - Optimisation Loop:   3%|▎         | 3/100 [00:03<01:17,  1.25it/s]
NAI - Optimisation Loop:   5%|▌         | 5/100 [00:03<00:40,  2.35it/s]
NAI - Optimisation Loop:   7%|▋         | 7/100 [00:03<00:25,  3.61it/s]
NAI - Optimisation Loop:   9%|▉         | 9/100 [00:03<00:18,  5.00it/s]
NAI - Optimisation Loop:  11%|█         | 11/100 [00:03<00:13,  6.43it/s]
NAI - Optimisation Loop:  13%|█▎        | 13/100 [00:03<00:11,  7.82it/s]
NAI - Optimisation Loop:  15%|█▌        | 15/100 [00:03<00:09,  9.05it/s]
NAI - Optimisation Loop:  17%|█▋        | 17/100 [00:04<00:08, 10.20it/s]
NAI - Optimisation Loop:  19%|█▉        | 19/100 [00:04<00:07, 11.10it/s]
NAI - Optimisation Loop:  21%|██        | 21/100 [00:04<00:06, 11.69it/s]
NAI - Optimisation Loop:  23%|██▎       | 23/100 [00:04<00:06, 12.24it/s]
NAI - Optimisation Loop:  25%|██▌       | 25/100 [00:04<00:06, 12.30it/s]
NAI - Optimisation Loop:  27%|██▋       | 27/100 [00:04<00:05, 12.60it/s]
NAI - Optimisation Loop:  29%|██▉       | 29/100 [00:04<00:05, 12.87it/s]
NAI - Optimisation Loop:  31%|███       | 31/100 [00:05<00:05, 13.07it/s]
NAI - Optimisation Loop:  33%|███▎      | 33/100 [00:05<00:05, 13.36it/s]
NAI - Optimisation Loop:  35%|███▌      | 35/100 [00:05<00:04, 13.50it/s]
NAI - Optimisation Loop:  37%|███▋      | 37/100 [00:05<00:04, 13.72it/s]
NAI - Optimisation Loop:  39%|███▉      | 39/100 [00:05<00:04, 13.25it/s]
NAI - Optimisation Loop:  41%|████      | 41/100 [00:05<00:04, 13.10it/s]
NAI - Optimisation Loop:  43%|████▎     | 43/100 [00:05<00:04, 13.36it/s]
NAI - Optimisation Loop:  45%|████▌     | 45/100 [00:06<00:04, 13.58it/s]
NAI - Optimisation Loop:  47%|████▋     | 47/100 [00:06<00:03, 13.72it/s]
NAI - Optimisation Loop:  49%|████▉     | 49/100 [00:06<00:03, 13.84it/s]
NAI - Optimisation Loop:  51%|█████     | 51/100 [00:06<00:03, 13.93it/s]
NAI - Optimisation Loop:  53%|█████▎    | 53/100 [00:06<00:03, 13.42it/s]
NAI - Optimisation Loop:  55%|█████▌    | 55/100 [00:06<00:03, 13.20it/s]
NAI - Optimisation Loop:  57%|█████▋    | 57/100 [00:07<00:03, 13.42it/s]
NAI - Optimisation Loop:  59%|█████▉    | 59/100 [00:07<00:03, 13.59it/s]
NAI - Optimisation Loop:  61%|██████    | 61/100 [00:07<00:02, 13.63it/s]
NAI - Optimisation Loop:  63%|██████▎   | 63/100 [00:07<00:02, 13.80it/s]
NAI - Optimisation Loop:  65%|██████▌   | 65/100 [00:07<00:02, 13.89it/s]
NAI - Optimisation Loop:  67%|██████▋   | 67/100 [00:07<00:02, 13.45it/s]
NAI - Optimisation Loop:  69%|██████▉   | 69/100 [00:07<00:02, 12.99it/s]
NAI - Optimisation Loop:  71%|███████   | 71/100 [00:08<00:02, 13.25it/s]
NAI - Optimisation Loop:  73%|███████▎  | 73/100 [00:08<00:02, 13.48it/s]
NAI - Optimisation Loop:  75%|███████▌  | 75/100 [00:08<00:01, 13.53it/s]
NAI - Optimisation Loop:  77%|███████▋  | 77/100 [00:08<00:01, 13.41it/s]
NAI - Optimisation Loop:  79%|███████▉  | 79/100 [00:08<00:01, 13.59it/s]
NAI - Optimisation Loop:  81%|████████  | 81/100 [00:08<00:01, 13.28it/s]
NAI - Optimisation Loop:  83%|████████▎ | 83/100 [00:08<00:01, 12.83it/s]
NAI - Optimisation Loop:  85%|████████▌ | 85/100 [00:09<00:01, 12.66it/s]
NAI - Optimisation Loop:  87%|████████▋ | 87/100 [00:09<00:00, 13.03it/s]
NAI - Optimisation Loop:  89%|████████▉ | 89/100 [00:09<00:00, 12.98it/s]
NAI - Optimisation Loop:  91%|█████████ | 91/100 [00:09<00:00, 12.74it/s]
NAI - Optimisation Loop:  93%|█████████▎| 93/100 [00:09<00:00, 12.83it/s]
NAI - Optimisation Loop:  95%|█████████▌| 95/100 [00:09<00:00, 12.38it/s]
NAI - Optimisation Loop:  97%|█████████▋| 97/100 [00:10<00:00, 12.16it/s]
NAI - Optimisation Loop:  99%|█████████▉| 99/100 [00:10<00:00, 12.23it/s]
NAI - Optimisation Loop: 100%|██████████| 100/100 [00:10<00:00,  9.69it/s]
The inversion result from `neighpy`:
============================
Summary for inversion result
============================
SUCCESS
----------------------------
model: [ 23.72857143  60.1643069  184.59628132  37.00981831  69.35691629]
direct_search_samples: [[ 28.5992124   51.76986253 160.76969216  36.68089787 113.94439725]
 [ 22.32906315  55.99849951 182.97693452  36.60204959  73.38038643]
 [ 26.85229621  73.32979535 160.83702553  32.88803251 106.42735064]
 ...
 [ 23.73103012  60.16619679 184.43034114  37.01024757  69.39008754]
 [ 23.73104021  60.16621931 184.43026621  37.01025554  69.39056737]
 [ 23.7311485   60.16617287 184.43015926  37.0102626   69.3905519 ]]
direct_search_objectives: [167021.88345958  14954.62891599 202454.94361509 ...  11714.95844992
  11714.9105204   11714.9379198 ]
appraisal_samples: [[ 21.23139207  57.1795549  184.72223782  36.83256783  69.83223415]
 [ 24.21589905  68.34014554 184.60095886  36.98444938  69.60130723]
 [ 24.79687323  62.1964014  184.72861247  37.09766029  69.83960096]
 ...
 [ 23.49542463  69.3660225  184.8174731   36.8674073   68.6628274 ]
 [ 28.76544302  58.15223988 184.91913196  37.34169443  67.8928283 ]
 [ 38.40301131  84.28735035 184.98616416  37.11189966  66.66546015]]

#@title plotting function (hidden)
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, my_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)

#@title plotting function (hidden)

arviz.style.use("arviz-variat")

import arviz_base
arviz_base.rcParams["plot.max_subplots"] = 80
pm = arviz.plot_pair(
    az_idata.sel(draw=slice(4000, None)),
    marginal=True,
    triangle="lower",
)

#@title plotting function (hidden)

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

Metropolis Hastings algorithm#

To sample the posterior distribution, BayesBay implements an RJ-MCMC method, which is a generalization of the Metropolis-Hastings algorithm allowing for trans-dimensional inference. Here we use BayesBay to solve a fixed dimensional problem as the number of thin plate targets is one. BayesBay is accessed from CoFI by using bayesbay as the string in inv_options.set_tool()

An example for CoFI using BayesBay to solve a trans-dimensional inverse problem is given here.

#@title defaults for bayesbay

def initialize_param(param, position=None, value=1):
    return numpy.array([value]) + 0.5 * numpy.random.randn()

parameters = []
for iparam, (vmin, vmax) in enumerate(zip(m_min, m_max)):
    parameter = bayesbay.prior.UniformPrior(
        name=f"m{iparam}",
        vmin=m_min[iparam],
        vmax=m_max[iparam],
        perturb_std=(vmax - vmin) / 20,
    )
    custom_init = functools.partial(initialize_param, value=init_param_value[iparam])
    parameter.set_custom_initialize(custom_init)
    parameters.append(parameter)

param_space = bayesbay.parameterization.ParameterSpace(
    name="param_space",
    n_dimensions=1,
    parameters=parameters,
)
parameterization = bayesbay.parameterization.Parameterization(param_space)

#@title wrap log likelihood function for bayesbay
def my_log_likelihood(state, *args, **kwargs):
    model = numpy.array(
        [state["param_space"][f"m{i}"] for i in range(init_param_value.size)]
    )
    return -0.5 * my_objective(model.T[0])

log_likelihood = bayesbay.likelihood.LogLikelihood(log_like_func=my_log_likelihood) # BayesBay version 3.1 and above
#log_likelihood = bayesbay.LogLikelihood(log_like_func=my_log_likelihood)  # BayesBay pre version 3.1

#@title bayesbay initialisation
n_chains = 12
walkers_start = []
for i in range(n_chains):
    walkers_start.append(
        parameterization.initialize()
    )  # A bayesbay.State is appended to walkers_start for each chain

Using the template below set the CoFI tool to bayesbay

inv_options = cofi.InversionOptions()
inv_options.set_tool(<DEFINE ME>)
inv_options.set_params(
    walkers_starting_states=walkers_start,
    perturbation_funcs=parameterization.perturbation_funcs,  # BayesBay version 3.1 and above
    #perturbation_funcs=parameterization.perturbation_functions, # BayesBay pre version 3.1
    log_like_ratio_func=log_likelihood,
    n_chains=n_chains,
    n_iterations=5_000,
    burnin_iterations=500,
    verbose=False,
    save_every=25,
)

# Copy the template above, Replace <DEFINE ME> with your answer

#@title Solution
inv_options = cofi.InversionOptions()
inv_options.set_tool("bayesbay")
inv_options.set_params(
    walkers_starting_states=walkers_start,
    perturbation_funcs=parameterization.perturbation_funcs,  # BayesBay version 3.1 and above
    #perturbation_funcs=parameterization.perturbation_functions, # BayesBay pre version 3.1
    log_like_ratio_func=log_likelihood,
    n_chains=n_chains,
    n_iterations=5_000,
    burnin_iterations=500,
    verbose=False,
    save_every=25,
)

inv = cofi.Inversion(cofi.BaseProblem(), inv_options)
my_result = inv.run()

results = my_result.models

#@title plotting function (hidden)
import arviz_base

arviz.style.use("arviz-variat")
var_names = [
    "Dip (°)",
    "Dip Azimuth (°)",
    "Easting (m)",
    "Depth (m)",
    "Width (m)",
]
clean_names = ["Dip_deg", "Dip_Azimuth_deg", "Easting_m", "Depth_m", "Width_m"]

results = my_result.models
posterior_samples = {
    clean_names[i]: numpy.concatenate(results[f"param_space.m{i}"])[numpy.newaxis, :]
    for i in range(init_param_value.size)
}

true_values = {
    var_names[i]: true_param_value[i] for i in range(init_param_value.size)
}

arviz_base.rcParams["plot.max_subplots"] = 80
az_idata = arviz_base.from_dict({"posterior": posterior_samples})
pm = arviz.plot_pair(
    az_idata,
    marginal=True,
    triangle="lower",
)
ref_vals = list(true_values.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,
            )

#@title plotting function (hidden)

_, 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()
idraw = numpy.random.randint(0, (posterior_samples[clean_names[0]].shape[1]))
sample = numpy.array([posterior_samples[name][0, idraw] for name in clean_names])

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
idraws = numpy.random.choice(
    numpy.arange(0, (posterior_samples[clean_names[0]].shape[1])), 10, replace=False
)
for idraw in idraws:
    sample = numpy.array([posterior_samples[name][0, idraw] for name in clean_names])
    plot_plate_faces(
        "plate_inverted",
        forward,
        sample,
        axes[0, 0],
        axes[0, 1],
        axes[1, 0],
        color="red",
        label="Posterior sample",
        linestyle="dotted",
    )

Where to next?#

This tutorial is a based on the material available as a a CoFI example under the following link

inlab-geo/cofi-examples

The following two notebooks explore the inversion of a field data set collected over the Caber deposit using the methods introduced in this tutorial. - Preprocessing - Inversion

Watermark#

watermark_list = ["cofi", "numpy", "scipy", "matplotlib","bayesbay","smt","neighpy"]
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
bayesbay 0.3.10
smt 2.13.0
neighpy 0.1.9

sphinx_gallery_thumbnail_number = -1

Total running time of the script: (11 minutes 50.796 seconds)

Gallery generated by Sphinx-Gallery

Thin plate inversion

Contents

Thin plate inversion#

What we do in this notebook#

Learning outcomes#

Problem description#

Forward solver#

Jacobian via finite differencing#

VTEMmax AEM system#

Plate parametrisation#

Implementation details#

Further reading#

Problem definition#

Inverting for the dip of a thin plate#

Problem setup#

True model#

Generate synthetic data#

Starting model#

Define helper functions for CoFI#

Challenge: Choose relative step size for the numerical Jacobian#

Plot the misfit and gradient as a function of the plate dip#

Parameter estimation#

Further reading#

Define CoFI problem#

Define CoFI options#

CoFI inversion#

Plotting#

Data#

Model#

Ensemble method#

Challenge: Given the objective function define a log likelihood function.#

Augment the CoFI problem#

Define CoFI options#

CoFI Inversion#

Limited computational resources#

Inverting for a thin plate given three survey lines#

Problem setup#

True model#

Generate synthetic data#

Challenge: Implement a parameter estimation or ensemble method in CoFI#

Parameter estimation applied to three survey lines#

Challenge: Choose a parameter estimation method#

Plotting#

Data - Profiles#

Model#

Ensemble methods applied to three survey lines#

Challenge: Select an ensemble method#

Affine Invariant Markov chain Monte Carlo Ensemble sampler#

Neighbourhood algorithm#

Metropolis Hastings algorithm#

Where to next?#

Watermark#