.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/generated/thin_plate_inversion.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_generated_thin_plate_inversion.py: Thin plate inversion ==================== .. GENERATED FROM PYTHON SOURCE LINES 9-14 |Open In Colab| .. |Open In Colab| image:: https://img.shields.io/badge/open%20in-Colab-b5e2fa?logo=googlecolab&style=flat-square&color=ffd670 :target: https://colab.research.google.com/github/inlab-geo/cofi-examples/blob/main/tutorials/thin_plate_inversion/thin_plate_inversion.ipynb .. GENERATED FROM PYTHON SOURCE LINES 17-33 -------------- 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`` `__. -------------- .. GENERATED FROM PYTHON SOURCE LINES 36-45 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 .. GENERATED FROM PYTHON SOURCE LINES 45-53 .. code-block:: Python # 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 .. GENERATED FROM PYTHON SOURCE LINES 55-61 .. code-block:: Python # 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 .. GENERATED FROM PYTHON SOURCE LINES 66-174 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 :math:`2 \mathrm{S}` located in a halfspace with a resistivity of :math:`1000 \mathrm{\Omega m}` with a :math:`20 \mathrm{m}` thick covering layer that has a resistivity of :math:`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 :math:`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 :math:`h` might result in a Jacobian that is not representative of the derivatives at location :math:`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 :math:`q` with :math:`h` defined as :math:`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. Further reading ~~~~~~~~~~~~~~~ Hauser, J., Gunning, J., & Annetts, D. (2016). Probabilistic inversion of airborne electromagnetic data for basement conductors. Geophysics, 81(5), E389-E400. Prikhodko, A., Morrison, E., Bagrianski, A., Kuzmin, P., Tishin, P., & Legault, J. (2010). Evolution of VTEM? technical solutions for effective exploration. ASEG Extended Abstracts, 2010(1), 1-4. Raiche, A., Sugeng, F. and Wilson, G. (2007) Practical 3D EM inversion the P223F software suite, ASEG Extended Abstracts, 2007:1, 1-5 Wheelock, B., Constable, S., & Key, K. (2015). The advantages of logarithmically scaled data for electromagnetic inversion. Geophysical Journal International, 201(3), 1765–1780. https://doi.org/10.1093/GJI/GGV107 .. GENERATED FROM PYTHON SOURCE LINES 174-202 .. code-block:: Python # import libraries needed for this example import pickle import functools import numpy import matplotlib.pyplot as plt from matplotlib.lines import Line2D import arviz import cofi import bayesbay import smt import smt.surrogate_models import tqdm from vtem_max_forward_lib import ( problem_setup, system_spec, survey_setup, true_model, ForwardWrapper, plot_transient, plot_predicted_profile, plot_plate_faces, plot_plate_faces_single ) numpy.random.seed(42) .. GENERATED FROM PYTHON SOURCE LINES 207-222 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. .. GENERATED FROM PYTHON SOURCE LINES 222-237 .. code-block:: Python 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 } .. GENERATED FROM PYTHON SOURCE LINES 242-258 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 :math:`60\degree` .. GENERATED FROM PYTHON SOURCE LINES 258-262 .. code-block:: Python forward = ForwardWrapper(true_model, problem_setup, system_spec, survey_setup, ["pdip"]) true_param_value = numpy.array([60]) .. rst-class:: sphx-glr-script-out .. code-block:: none ['pdip'] .. GENERATED FROM PYTHON SOURCE LINES 267-270 Problem setup ------------- .. GENERATED FROM PYTHON SOURCE LINES 273-279 True model ~~~~~~~~~~ We first the plot the true model and the location of the VTEMmax measurement, which we will be inverting. .. GENERATED FROM PYTHON SOURCE LINES 279-295 .. code-block:: Python #@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") .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_001.png :alt: thin plate inversion :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 300-311 Generate synthetic data ~~~~~~~~~~~~~~~~~~~~~~~ This tutorial uses a simplified noise model that assumes an absolute noise, that is a standard deviation of :math:`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. .. GENERATED FROM PYTHON SOURCE LINES 311-323 .. code-block:: Python # 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)) .. GENERATED FROM PYTHON SOURCE LINES 328-333 Starting model ~~~~~~~~~~~~~~ Set an initial guess for the dip of the thin plate. .. GENERATED FROM PYTHON SOURCE LINES 333-336 .. code-block:: Python init_param_value = numpy.array([45]) .. GENERATED FROM PYTHON SOURCE LINES 341-344 Define helper functions for CoFI -------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 347-369 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* |Upload to Excalidraw_1| .. |Upload to Excalidraw_1| image:: https://img.shields.io/badge/Click%20&%20upload%20your%20results%20to-Excalidraw-lightgrey?logo=jamboard&style=for-the-badge&color=fcbf49&labelColor=edede9 :target: https://excalidraw.com/#room=f4481f8278ad2ddcb96d,i9Ki_GouExK4GylmrsrZ2A .. GENERATED FROM PYTHON SOURCE LINES 372-404 Using the template below se the relative step size :: my_relative_step = 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 .. GENERATED FROM PYTHON SOURCE LINES 404-407 .. code-block:: Python # Copy the template above, Replace with your answer .. GENERATED FROM PYTHON SOURCE LINES 409-440 .. code-block:: Python #@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 .. GENERATED FROM PYTHON SOURCE LINES 445-448 Plot the misfit and gradient as a function of the plate dip ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 448-478 .. code-block:: Python #@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() .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_002.png :alt: thin plate inversion :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none 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] .. GENERATED FROM PYTHON SOURCE LINES 483-529 Parameter estimation -------------------- First we solve the inverse problem using optimisation, that is we seek to find the minimum of the objective function given as .. math:: \chi^2 = (\mathbf{d} - \mathbf{f}(\mathbf{m}))^T\mathbf{C}_d^{-1}(\mathbf{d}-\mathbf{f}(\mathbf{m})), where :math:`\mathbf{d}` is the observed data, :math:`\mathbf{f}(\mathbf{m})` the model prediction and :math:`\mathbf{C}_d` the data covariance matrix. The full Newton step is then given as .. math:: \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 :math:`\mathbf{J}` and Hessian, here approximated by :math:`\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. Further reading ~~~~~~~~~~~~~~~ https://medium.com/intro-to-artificial-intelligence/line-search-and-trust-region-optimisation-strategies-638a4a7490ca .. GENERATED FROM PYTHON SOURCE LINES 532-535 Define CoFI problem ^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 535-542 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 547-550 Define CoFI options ^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 550-555 .. code-block:: Python my_options = cofi.InversionOptions() my_options.set_tool("scipy.optimize.minimize") my_options.set_params(method="Newton-CG",callback=PerIterationCallbackFunction()) .. GENERATED FROM PYTHON SOURCE LINES 560-563 CoFI inversion ^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 563-571 .. code-block:: Python my_inversion = cofi.Inversion(my_problem, my_options) my_result = my_inversion.run() print(f"\nNumber of objective function evaluations: {my_result.nfev}") print(f"Number of gradient function evaluations: {my_result.njev}") print(f"Number of hessian function evaluations: {my_result.nhev}") print(f"Solution vector:\n",my_result.model) .. rst-class:: sphx-glr-script-out .. code-block:: none Iteration #1 objective value: 77.89136366904667 Iteration #2 objective value: 77.67785905000744 Number of objective function evaluations: 30 Number of gradient function evaluations: 11 Number of hessian function evaluations: 3 Solution vector: [59.18480808] .. GENERATED FROM PYTHON SOURCE LINES 576-579 Plotting ~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 582-589 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. .. GENERATED FROM PYTHON SOURCE LINES 589-602 .. code-block:: Python #@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() .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_003.png :alt: vertical, inline :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 607-610 Model ^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 610-636 .. code-block:: Python #@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") .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_004.png :alt: thin plate inversion :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 641-662 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 :math:`-\frac{1}{2}` .. math:: 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 :math:`10 \degree` and an upper boundary of :math:`80 \degree`. .. GENERATED FROM PYTHON SOURCE LINES 662-671 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 676-681 Challenge: Given the objective function define a log likelihood function. ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ In the previous section we defined the following objective function. .. GENERATED FROM PYTHON SOURCE LINES 681-687 .. code-block:: Python def my_objective(model): dpred = forward(model) residual = dpred - data_obs return residual.T @ Cdinv @ residual .. GENERATED FROM PYTHON SOURCE LINES 692-695 This function can be used to now create the log likelihood function, typically needed by the ensemble methods made available in CoFI. .. GENERATED FROM PYTHON SOURCE LINES 695-699 .. code-block:: Python def my_log_likelihood(model): return # DEFINE ME .. GENERATED FROM PYTHON SOURCE LINES 701-706 .. code-block:: Python #@title Solution def my_log_likelihood(model): return -0.5 * my_objective(model) .. GENERATED FROM PYTHON SOURCE LINES 711-720 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. .. GENERATED FROM PYTHON SOURCE LINES 720-725 .. code-block:: Python 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)) .. GENERATED FROM PYTHON SOURCE LINES 730-733 Define CoFI options ~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 733-739 .. code-block:: Python nwalkers = 2 ndim = len(init_param_value) nsteps = 50 walkers_start = init_param_value + 1 * numpy.random.randn(nwalkers, ndim) .. GENERATED FROM PYTHON SOURCE LINES 744-747 CoFI Inversion ~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 747-760 .. code-block:: Python 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() .. rst-class:: sphx-glr-script-out .. code-block:: none 0%| | 0/50 [00:00 blob_names: ['log_likelihood', 'log_prior'] .. GENERATED FROM PYTHON SOURCE LINES 762-787 .. code-block:: Python #@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) .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_005.png :alt: plate dip (°), plate dip (°) :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 792-798 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. .. GENERATED FROM PYTHON SOURCE LINES 801-819 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 `__. .. GENERATED FROM PYTHON SOURCE LINES 819-837 .. code-block:: Python #@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) .. rst-class:: sphx-glr-script-out .. code-block:: none 0%| | 0/25 [00:00 blob_names: ['log_likelihood', 'log_prior'] .. GENERATED FROM PYTHON SOURCE LINES 924-949 .. code-block:: Python #@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) .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_007.png :alt: plate dip (°), plate dip (°) :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 954-962 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. .. GENERATED FROM PYTHON SOURCE LINES 965-970 Problem setup ------------- In the following we define three survey lines covering the thin plate .. GENERATED FROM PYTHON SOURCE LINES 970-984 .. code-block:: Python 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() .. GENERATED FROM PYTHON SOURCE LINES 986-993 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 995-1012 .. code-block:: Python 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 } .. GENERATED FROM PYTHON SOURCE LINES 1017-1020 True model ~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 1020-1036 .. code-block:: Python 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]) } .. GENERATED FROM PYTHON SOURCE LINES 1041-1054 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. .. GENERATED FROM PYTHON SOURCE LINES 1054-1058 .. code-block:: Python forward = ForwardWrapper(true_model, problem_setup, system_spec, survey_setup, ["pdip","pdzm", "peast", "pwdth2","ptop"], data_returned=["vertical"]); .. rst-class:: sphx-glr-script-out .. code-block:: none ['pdip', 'pdzm', 'peast', 'pwdth2', 'ptop'] .. GENERATED FROM PYTHON SOURCE LINES 1060-1064 .. code-block:: Python # check the order of parameters in a model vector forward.params_to_invert .. rst-class:: sphx-glr-script-out .. code-block:: none ['pdip', 'pdzm', 'peast', 'ptop', 'pwdth2'] .. GENERATED FROM PYTHON SOURCE LINES 1066-1069 .. code-block:: Python true_param_value = numpy.array([60,65, 175, 30, 90]) .. GENERATED FROM PYTHON SOURCE LINES 1071-1089 .. code-block:: Python #@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") .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_008.png :alt: thin plate inversion :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_008.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 1094-1100 Generate synthetic data ~~~~~~~~~~~~~~~~~~~~~~~ We again generate a synthetic data set and add a realisation of the noise. .. GENERATED FROM PYTHON SOURCE LINES 1100-1112 .. code-block:: Python # 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)) .. GENERATED FROM PYTHON SOURCE LINES 1117-1149 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 <#Parameter-estimation-applied-to-three-survey-lines>`__ - `Ensemble methods <#Ensemble-methods-applied-to-three-survey-lines>`__ .. container:: alert alert-block alert-info 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* |Upload to Excalidraw_2| .. |Upload to Excalidraw_2| image:: https://img.shields.io/badge/Click%20&%20upload%20your%20results%20to-Excalidraw-lightgrey?logo=jamboard&style=for-the-badge&color=fcbf49&labelColor=edede9 :target: https://excalidraw.com/#room=52f0ac5f10e0111ee085,3nYggMVJOpmqlV1ZbYj0Eg .. GENERATED FROM PYTHON SOURCE LINES 1152-1155 Parameter estimation applied to three survey lines -------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 1158-1160 **Initialise a model for inversion** .. GENERATED FROM PYTHON SOURCE LINES 1160-1163 .. code-block:: Python forward.params_to_invert .. rst-class:: sphx-glr-script-out .. code-block:: none ['pdip', 'pdzm', 'peast', 'ptop', 'pwdth2'] .. GENERATED FROM PYTHON SOURCE LINES 1165-1168 .. code-block:: Python init_param_value = numpy.array([45, 90, 150, 20, 80]) .. GENERATED FROM PYTHON SOURCE LINES 1173-1175 **Define helper functions for CoFI** .. GENERATED FROM PYTHON SOURCE LINES 1175-1202 .. code-block:: Python 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 .. GENERATED FROM PYTHON SOURCE LINES 1207-1209 **Define CoFI problem** .. GENERATED FROM PYTHON SOURCE LINES 1209-1216 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 1221-1255 Challenge: Choose a parameter estimation method ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ CoFI provides access to the parameter estimation methods that are available in ```scipy.optimize.minimize`` `__ 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. .. container:: alert alert-block alert-info 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 |Upload to Excalidraw_2| .. |Upload to Excalidraw_2| image:: https://img.shields.io/badge/Click%20&%20upload%20your%20results%20to-Excalidraw-lightgrey?logo=jamboard&style=for-the-badge&color=fcbf49&labelColor=edede9 :target: https://excalidraw.com/#room=52f0ac5f10e0111ee085,3nYggMVJOpmqlV1ZbYj0Eg .. GENERATED FROM PYTHON SOURCE LINES 1258-1260 **Define CoFI options** .. GENERATED FROM PYTHON SOURCE LINES 1263-1272 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=,callback=PerIterationCallbackFunction(),options={"maxiter": 10}) .. GENERATED FROM PYTHON SOURCE LINES 1272-1275 .. code-block:: Python # Copy the template above, Replace with your answer .. GENERATED FROM PYTHON SOURCE LINES 1277-1283 .. code-block:: Python #@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}) .. GENERATED FROM PYTHON SOURCE LINES 1285-1291 .. code-block:: Python #@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}) .. GENERATED FROM PYTHON SOURCE LINES 1296-1298 **Run CoFI inversion** .. GENERATED FROM PYTHON SOURCE LINES 1298-1302 .. code-block:: Python my_inversion = cofi.Inversion(my_problem, my_options) my_result = my_inversion.run() .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 1307-1310 Plotting ~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 1313-1316 Data - Profiles ^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 1316-1335 .. code-block:: Python #@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() .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_009.png :alt: Crossline 25 (m), Crossline 100 (m), Crossline 175 (m) :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1340-1343 Model ^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 1343-1371 .. code-block:: Python #@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") .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_010.png :alt: thin plate inversion :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_010.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 1376-1393 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 `__ .. GENERATED FROM PYTHON SOURCE LINES 1393-1418 .. code-block:: Python #@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() .. rst-class:: sphx-glr-script-out .. code-block:: none ___________________________________________________________________________ Kriging ___________________________________________________________________________ Problem size # training points. : 100 ___________________________________________________________________________ Training Training ... Training - done. Time (sec): 4.5371859 .. GENERATED FROM PYTHON SOURCE LINES 1423-1425 **Initialise a model for inversion and define prior distribution** .. GENERATED FROM PYTHON SOURCE LINES 1425-1430 .. code-block:: Python 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]) .. GENERATED FROM PYTHON SOURCE LINES 1435-1437 **Define helper functions for CoFI** .. GENERATED FROM PYTHON SOURCE LINES 1437-1454 .. code-block:: Python 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) .. GENERATED FROM PYTHON SOURCE LINES 1459-1489 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>`__ .. container:: alert alert-block alert-info 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 |Upload to Excalidraw_2| .. |Upload to Excalidraw_2| image:: https://img.shields.io/badge/Click%20&%20upload%20your%20results%20to-Excalidraw-lightgrey?logo=jamboard&style=for-the-badge&color=fcbf49&labelColor=edede9 :target: https://excalidraw.com/#room=52f0ac5f10e0111ee085,3nYggMVJOpmqlV1ZbYj0Eg .. GENERATED FROM PYTHON SOURCE LINES 1492-1507 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 `__). |Upload to Excalidraw_2| .. |Upload to Excalidraw_2| image:: https://img.shields.io/badge/Click%20&%20upload%20your%20results%20to-Excalidraw-lightgrey?logo=jamboard&style=for-the-badge&color=fcbf49&labelColor=edede9 :target: https://excalidraw.com/#room=52f0ac5f10e0111ee085,3nYggMVJOpmqlV1ZbYj0Eg .. GENERATED FROM PYTHON SOURCE LINES 1507-1513 .. code-block:: Python 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)) .. GENERATED FROM PYTHON SOURCE LINES 1518-1520 **Define CoFI options** .. GENERATED FROM PYTHON SOURCE LINES 1520-1527 .. code-block:: Python #@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) .. GENERATED FROM PYTHON SOURCE LINES 1532-1548 Use the template below and set the CoFI tool to ``emcee`` :: inv_options = cofi.InversionOptions() inv_options.set_tool() 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() .. GENERATED FROM PYTHON SOURCE LINES 1548-1551 .. code-block:: Python # Copy the template above, Replace with your answer .. GENERATED FROM PYTHON SOURCE LINES 1553-1568 .. code-block:: Python #@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() .. rst-class:: sphx-glr-script-out .. code-block:: none 0%| | 0/5000 [00:00 blob_names: ['log_likelihood', 'log_prior'] .. GENERATED FROM PYTHON SOURCE LINES 1570-1602 .. code-block:: Python #@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) .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_011.png :alt: Dip (°), Dip (°), Dip azimuth (°), Dip azimuth (°), Easting (m), Easting (m), Depth (m), Depth (m), Width (m), Width (m) :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_011.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1604-1635 .. code-block:: Python #@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, ) .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_012.png :alt: thin plate inversion :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_012.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1637-1712 .. code-block:: Python #@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" ) .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_013.png :alt: thin plate inversion :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_013.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1717-1740 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()`` |Upload to Excalidraw_2| .. |Upload to Excalidraw_2| image:: https://img.shields.io/badge/Click%20&%20upload%20your%20results%20to-Excalidraw-lightgrey?logo=jamboard&style=for-the-badge&color=fcbf49&labelColor=edede9 :target: https://excalidraw.com/#room=52f0ac5f10e0111ee085,3nYggMVJOpmqlV1ZbYj0Eg .. GENERATED FROM PYTHON SOURCE LINES 1740-1744 .. code-block:: Python my_problem = cofi.BaseProblem() my_problem.set_objective(my_objective) .. GENERATED FROM PYTHON SOURCE LINES 1749-1757 Using the template below set the CoFI tool to ``neighpy`` :: inv_options = cofi.InversionOptions() inv_options.set_tool() inv_options.suggest_solver_params() .. GENERATED FROM PYTHON SOURCE LINES 1757-1760 .. code-block:: Python # Copy the template above, Replace with your answer .. GENERATED FROM PYTHON SOURCE LINES 1762-1768 .. code-block:: Python #@title Solution inv_options = cofi.InversionOptions() inv_options.set_tool("neighpy") inv_options.suggest_solver_params() .. rst-class:: sphx-glr-script-out .. code-block:: none 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} .. GENERATED FROM PYTHON SOURCE LINES 1770-1789 .. code-block:: Python 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() .. rst-class:: sphx-glr-script-out .. code-block:: none NAI - Initial Random Search NAI - Optimisation Loop: 0%| | 0/100 [00:00`__ 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()`` |Upload to Excalidraw_2| An example for CoFI using `BayesBay `__ to solve a trans-dimensional inverse problem is given `here `__. .. |Upload to Excalidraw_2| image:: https://img.shields.io/badge/Click%20&%20upload%20your%20results%20to-Excalidraw-lightgrey?logo=jamboard&style=for-the-badge&color=fcbf49&labelColor=edede9 :target: https://excalidraw.com/#room=52f0ac5f10e0111ee085,3nYggMVJOpmqlV1ZbYj0Eg .. GENERATED FROM PYTHON SOURCE LINES 1930-1955 .. code-block:: Python #@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) .. GENERATED FROM PYTHON SOURCE LINES 1957-1968 .. code-block:: Python #@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 .. GENERATED FROM PYTHON SOURCE LINES 1970-1979 .. code-block:: Python #@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 .. GENERATED FROM PYTHON SOURCE LINES 1984-2002 Using the template below set the CoFI tool to ``bayesbay`` :: inv_options = cofi.InversionOptions() inv_options.set_tool() 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, ) .. GENERATED FROM PYTHON SOURCE LINES 2002-2005 .. code-block:: Python # Copy the template above, Replace with your answer .. GENERATED FROM PYTHON SOURCE LINES 2007-2023 .. code-block:: Python #@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, ) .. GENERATED FROM PYTHON SOURCE LINES 2025-2029 .. code-block:: Python inv = cofi.Inversion(cofi.BaseProblem(), inv_options) my_result = inv.run() .. GENERATED FROM PYTHON SOURCE LINES 2031-2034 .. code-block:: Python results = my_result.models .. GENERATED FROM PYTHON SOURCE LINES 2036-2084 .. code-block:: Python #@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, ) .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_017.png :alt: thin plate inversion :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_017.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 2086-2162 .. code-block:: Python #@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", ) .. image-sg:: /tutorials/generated/images/sphx_glr_thin_plate_inversion_018.png :alt: thin plate inversion :srcset: /tutorials/generated/images/sphx_glr_thin_plate_inversion_018.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 2167-2182 Where to next? ============== This tutorial is a based on the material available as a a CoFI example under the following link https://github.com/inlab-geo/cofi-examples/tree/main/examples/vtem_max 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 `__ .. GENERATED FROM PYTHON SOURCE LINES 2185-2198 -------------- Watermark ========= .. raw:: html .. raw:: html .. GENERATED FROM PYTHON SOURCE LINES 2198-2204 .. code-block:: Python watermark_list = ["cofi", "numpy", "scipy", "matplotlib","bayesbay","smt","neighpy"] for pkg in watermark_list: pkg_var = __import__(pkg) print(pkg, getattr(pkg_var, "__version__")) .. rst-class:: sphx-glr-script-out .. code-block:: none 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 .. GENERATED FROM PYTHON SOURCE LINES 2205-2205 sphinx_gallery_thumbnail_number = -1 .. rst-class:: sphx-glr-timing **Total running time of the script:** (11 minutes 50.796 seconds) .. _sphx_glr_download_tutorials_generated_thin_plate_inversion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: thin_plate_inversion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: thin_plate_inversion.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: thin_plate_inversion.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_