.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/generated/scripts_synth_data/caber_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_examples_generated_scripts_synth_data_caber_inversion.py: Inversion of the Caber VTEM max survey for a thin plate target ============================================================== .. raw:: html .. GENERATED FROM PYTHON SOURCE LINES 13-18 |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/examples/airborne_em/airborne_em_caber_dataset_verification.ipynb .. GENERATED FROM PYTHON SOURCE LINES 21-40 .. raw:: html .. raw:: html .. If you are running this notebook locally, make sure you’ve followed `steps here `__ to set up the environment. (This `environment.yml `__ file specifies a list of packages required to run the notebooks) .. GENERATED FROM PYTHON SOURCE LINES 40-52 .. code-block:: Python # -------------------------------------------------------- # # # # Uncomment below to set up environment on "colab" # # # # -------------------------------------------------------- # # !pip install -U cofi # !pip install git+https://github.com/JuergHauser/PyP223.git # !git clone https://github.com/inlab-geo/cofi-examples.git # %cd cofi-examples/examples/vtem_max .. GENERATED FROM PYTHON SOURCE LINES 54-59 .. code-block:: Python # If this notebook is run locally PyP223 needs to be installed separately by uncommenting the following line, # that is by removing the # and the white space between it and the exclamation mark. # !pip install git+https://github.com/JuergHauser/PyP223.git .. GENERATED FROM PYTHON SOURCE LINES 61-84 .. code-block:: Python import numpy import pandas import matplotlib.pyplot as plt from matplotlib.lines import Line2D import cofi import json from vtem_max_forward_lib import ( ForwardWrapper, problem_setup, plot_field_data, plot_transient, plot_predicted_profile, plot_observed_profile, get_subset_of_survey, plot_plate_faces, plot_plate_faces_single, plot_survey_map ) numpy.random.seed(42) .. GENERATED FROM PYTHON SOURCE LINES 89-113 Background ---------- (Prikhodko et. al,2010) The Caber orebody is a volcanogenic massive sulfide (VMS) deposit in western Quebec, Canada, that hosts significant Zn-Cu-Ag mineralisation and VTEM max surveys have been collected over the target and previously been inverted recovering subvertical body dipping at about 80 degrees to the southwest. Here we invert the vertical compnent for a starting model with a dip azimuth of :math:`180 \degree` and expect the inversion to change the dip azimuth to a southwesterly direction. Further Reading ''''''''''''''' 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. McMillan, M. S., Schwarzbach, C., Haber, E., & Oldenburg, D. W. (2015). 3D parametric hybrid inversion of time-domain airborne electromagnetic data. Geophysics, 80(6), K25-K36. .. GENERATED FROM PYTHON SOURCE LINES 116-119 Decimated survey data --------------------- .. GENERATED FROM PYTHON SOURCE LINES 119-126 .. code-block:: Python # Load the data files with open('caber_survey.npy', 'rb') as f: survey_setup = numpy.load(f,allow_pickle=True) system_spec = numpy.load(f,allow_pickle=True) data_obs = numpy.load(f,allow_pickle=True) .. GENERATED FROM PYTHON SOURCE LINES 128-132 .. code-block:: Python survey_setup=survey_setup.tolist() system_spec=system_spec.tolist() .. GENERATED FROM PYTHON SOURCE LINES 134-137 .. code-block:: Python plot_survey_map(survey_setup) .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_001.png :alt: caber inversion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 139-147 .. code-block:: Python _, axes = plt.subplots(5, 1, figsize=(7, 10)) for i in range(5): plot_observed_profile(survey_setup, data_obs, label="observed data", gate_idx=numpy.arange(0, 34, 2), line_id=[i], ax=axes[i], color="pink") plt.tight_layout() .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_002.png :alt: caber inversion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 149-166 .. code-block:: Python # initial guess initial_model = { "res": numpy.array([1000, 1000]), "thk": numpy.array([25]), "peast": numpy.array([710200]), "pnorth": numpy.array([5513450]), "ptop": numpy.array([150]), "pres": numpy.array([0.008]), "plngth1": numpy.array([150]), "plngth2": numpy.array([150]), "pwdth1": numpy.array([0.1]), "pwdth2": numpy.array([500]), "pdzm": numpy.array([225]), "pdip": numpy.array([80]) } .. GENERATED FROM PYTHON SOURCE LINES 168-171 .. code-block:: Python problem_setup['cellw']=100 .. GENERATED FROM PYTHON SOURCE LINES 173-177 .. code-block:: Python forward = ForwardWrapper(initial_model, problem_setup, system_spec, survey_setup, ["pdip","pdzm", "pwdth2","ptop","plngth1","plngth2"], data_returned=["vertical"]) .. rst-class:: sphx-glr-script-out .. code-block:: none ['pdip', 'pdzm', 'pwdth2', 'ptop', 'plngth1', 'plngth2'] .. GENERATED FROM PYTHON SOURCE LINES 179-182 .. code-block:: Python forward.params_to_invert .. rst-class:: sphx-glr-script-out .. code-block:: none ['pdip', 'pdzm', 'plngth1', 'plngth2', 'ptop', 'pwdth2'] .. GENERATED FROM PYTHON SOURCE LINES 184-188 .. code-block:: Python init_param_value=[90,200,150,150,150,500] init_param_value=[90,180,150,150,150,500] .. GENERATED FROM PYTHON SOURCE LINES 193-196 Parameter estimation -------------------- .. GENERATED FROM PYTHON SOURCE LINES 199-202 Define helper functions for CoFI ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 202-229 .. code-block:: Python def my_objective(model): dpred = forward(model) residual = dpred - data_obs return residual.T @ residual def my_gradient(model): dpred = forward(model) jacobian = forward.jacobian(model) residual = dpred - data_obs return 2 * jacobian.T @ residual def my_hessian(model): jacobian = forward.jacobian(model) return 2 * jacobian.T @ jacobian class CallbackFunction: 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 234-237 Define CoFI problem ~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 237-244 .. 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 249-252 Define CoFI options ~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 252-258 .. code-block:: Python my_options = cofi.InversionOptions() my_options.set_tool("scipy.optimize.minimize") #my_options.set_params(method="newton-cg", options={"maxiter": 5}, callback=CallbackFunction()) my_options.set_params(method="trust-ncg", options={"maxiter": 10}, callback=CallbackFunction()) .. GENERATED FROM PYTHON SOURCE LINES 263-266 Run CoFI inversion ~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 266-271 .. code-block:: Python my_inversion = cofi.Inversion(my_problem, my_options) my_result = my_inversion.run() print(my_result.model) .. rst-class:: sphx-glr-script-out .. code-block:: none Iteration #1 objective value: 3682.494724069316 Iteration #2 objective value: 3652.5047048741376 Iteration #3 objective value: 3598.26112692279 Iteration #4 objective value: 3510.5582707397016 Iteration #5 objective value: 3393.85818986716 Iteration #6 objective value: 3337.42629970419 Iteration #7 objective value: 3334.662725587539 Iteration #8 objective value: 3334.662722409527 Iteration #9 objective value: 3334.6627235490664 Iteration #10 objective value: 3334.6627247255874 [ 89.08492569 206.93166544 135.4257979 142.5526461 199.57746436 505.35124188] .. GENERATED FROM PYTHON SOURCE LINES 273-276 .. code-block:: Python my_result.message .. rst-class:: sphx-glr-script-out .. code-block:: none 'Maximum number of iterations has been exceeded.' .. GENERATED FROM PYTHON SOURCE LINES 278-286 .. code-block:: Python _, ax = plt.subplots() plot_transient(init_param_value, forward, "data from init model", ax, color="green", linestyle=":") plot_transient(my_result.model, forward, "data from inverted model", ax, color="red", linestyle="-.") ax.legend(loc="upper center") ax.set_title("vertical") plt.tight_layout() .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_003.png :alt: vertical :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 288-301 .. code-block:: Python idx_to_plot = numpy.arange(8, 30) # (0, 44) _, axes = plt.subplots(5, 1, figsize=(7,10)) for i in range(5): plot_predicted_profile(init_param_value, forward, "data from init model", gate_idx=idx_to_plot, line_id=[i], ax=axes[i], color="green", linestyle=":") plot_predicted_profile(my_result.model, forward, "data from inverted model", gate_idx=idx_to_plot, line_id=[i], ax=axes[i], color="red", linestyle="-.") ax.legend(loc="upper center") plt.tight_layout() .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_004.png :alt: caber inversion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 306-309 Despite only a small improvement in fit to the data the recovered thing plate shows the expected southwesterly dip direction. .. GENERATED FROM PYTHON SOURCE LINES 309-334 .. code-block:: Python _, axes = plt.subplots(2, 2) axes[1,1].axis("off") plot_plate_faces( "plate_init", forward, init_param_value, axes[0,0], axes[0,1], axes[1,0], color="green", label="Starting model", surface_elevation=0 ) 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", surface_elevation=0 ) 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:: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_005.png :alt: caber inversion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_caber_inversion_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 339-352 -------------- Watermark ========= .. raw:: html .. raw:: html .. GENERATED FROM PYTHON SOURCE LINES 352-358 .. code-block:: Python watermark_list = ["cofi", "numpy", "scipy", "matplotlib"] 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 numpy 2.3.5 scipy 1.17.0 matplotlib 3.10.8 .. GENERATED FROM PYTHON SOURCE LINES 359-359 sphinx_gallery_thumbnail_number = -1 .. rst-class:: sphx-glr-timing **Total running time of the script:** (25 minutes 4.009 seconds) .. _sphx_glr_download_examples_generated_scripts_synth_data_caber_inversion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: caber_inversion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: caber_inversion.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: caber_inversion.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_