.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/generated/scripts_synth_data/fmm_tomography_regularization_discussion.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_fmm_tomography_regularization_discussion.py: Seismic Travel time Tomography via Fast Marching - Demo on switching regularization and L-curve =============================================================================================== .. 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/examples/fmm_tomography/fmm_tomography.ipynb .. GENERATED FROM PYTHON SOURCE LINES 17-36 .. 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 39-48 .. raw:: html In this notebook, we would like to demonstrate the capability of CoFI to easily switch between different types of regularizations. We will use ``cofi`` to run a seismic tomography example. .. GENERATED FROM PYTHON SOURCE LINES 51-54 Theoretical background ---------------------- .. GENERATED FROM PYTHON SOURCE LINES 54-63 .. code-block:: Python # display theory on travel time tomography from IPython.display import display, Markdown with open("../../theory/geo_travel_time_tomography.md", "r") as f: content = f.read() display(Markdown(content)) .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 68-89 For forward modelling, a fast marching wave front tracker is used, utilizing the Fast Marching Fortran code within the package ```FMTOMO`` `__ by Nick Rawlinson. The Fast Marching code is wrapped in package `pyfm2d `__. Further details can be found in: - Rawlinson, N., de Kool, M. and Sambridge, M., 2006. Seismic wavefront tracking in 3-D heterogeneous media: applications with multiple data classes, Explor. Geophys., 37, 322-330. - Rawlinson, N. and Urvoy, M., 2006. Simultaneous inversion of active and passive source datasets for 3-D seismic structure with application to Tasmania, Geophys. Res. Lett., 33 L24313, 10.1029/2006GL028105. - de Kool, M., Rawlinson, N. and Sambridge, M. 2006. A practical grid based method for tracking multiple refraction and reflection phases in 3D heterogeneous media, Geophys. J. Int., 167, 253-270. - Saygin, E. 2007. Seismic receiver and noise correlation based studies in Australia, PhD thesis, Australian National University, 10.25911/5d7a2d1296f96. .. GENERATED FROM PYTHON SOURCE LINES 92-95 0. Import modules ----------------- .. GENERATED FROM PYTHON SOURCE LINES 95-104 .. code-block:: Python # -------------------------------------------------------- # # # # Uncomment below to set up environment on "colab" # # # # -------------------------------------------------------- # # !pip install -U cofi pyfm2d .. GENERATED FROM PYTHON SOURCE LINES 106-115 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt import pprint import cofi # NB You will need to separately install pyfm2d in your python environment with `pip install pyfm2d' import pyfm2d as wt # import fmm package .. GENERATED FROM PYTHON SOURCE LINES 120-130 Understanding the inference problem ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Before we starting working with ``cofi``, let’s get familiar with the problem itself. Below is a plot of the true model and the paths generated from this model. As you can see, there are two anomalies, one with lower velocity (red, top left) and the other with higher velocity (blue, bottom right). .. GENERATED FROM PYTHON SOURCE LINES 130-136 .. code-block:: Python # read in problem data loaded_dict = np.load('../../data/travel_time_tomography/nonlinear_tomo_example.npz') nonlinear_tomo_example = dict(loaded_dict) loaded_dict.close() .. GENERATED FROM PYTHON SOURCE LINES 138-148 .. code-block:: Python # set up problem good_model = nonlinear_tomo_example["_mtrue"] extent = nonlinear_tomo_example["extent"] sources = nonlinear_tomo_example["sources"] receivers = nonlinear_tomo_example["receivers"] obstimes = nonlinear_tomo_example["_data"] print(' New data set have:\n',len(receivers),' receivers\n',len(sources),' sources\n',len(obstimes),' travel times\n', 'Range of travel times: ',np.min(obstimes),'to',np.max(obstimes),'\n Mean travel time:',np.mean(obstimes)) .. rst-class:: sphx-glr-script-out .. code-block:: none New data set have: 10 receivers 10 sources 100 travel times Range of travel times: 0.009490006 to 0.01558705 Mean travel time: 0.011210016954999999 .. GENERATED FROM PYTHON SOURCE LINES 150-156 .. code-block:: Python # display true model and raypaths options = wt.WaveTrackerOptions(paths=True,cartesian=True) # set wavetracker options result = wt.calc_wavefronts(good_model,receivers,sources,extent=extent, options=options) # track wavefronts wt.display_model(good_model,paths=result.paths,extent=extent,line=0.3,alpha=0.82) .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_001.png :alt: fmm tomography regularization discussion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 161-164 1. Problem setup and utilities ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 164-171 .. code-block:: Python # get problem information model_size = good_model.size # number of model parameters model_shape = good_model.shape # 2D spatial grid shape data_size = data_size = len(obstimes) # number of data ref_start_slowness = nonlinear_tomo_example["_sstart"] # use the starting guess supplied by the nonlinear example .. GENERATED FROM PYTHON SOURCE LINES 173-222 .. code-block:: Python def objective_func(slowness, reg, sigma, reduce_data=None): # reduce_data=(idx_from, idx_to) if reduce_data is None: idx_from, idx_to = (0, data_size) else: idx_from, idx_to = reduce_data if(True): options = wt.WaveTrackerOptions( cartesian=True, ) result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options) # track wavefronts ttimes = result.ttimes residual = obstimes[idx_from:idx_to] - ttimes[idx_from:idx_to] data_misfit = residual.T @ residual / sigma**2 model_reg = reg(slowness) return data_misfit + model_reg def gradient(slowness, reg, sigma, reduce_data=None): # reduce_data=(idx_from, idx_to) if reduce_data is None: idx_from, idx_to = (0, data_size) else: idx_from, idx_to = reduce_data if(True): options = wt.WaveTrackerOptions( paths=True, frechet=True, cartesian=True, ) result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options) # track wavefronts ttimes = result.ttimes A = result.frechet.toarray() ttimes = ttimes[idx_from:idx_to] A = A[idx_from:idx_to] data_misfit_grad = -2 * A.T @ (obstimes[idx_from:idx_to] - ttimes) / sigma**2 model_reg_grad = reg.gradient(slowness) return data_misfit_grad + model_reg_grad def hessian(slowness, reg, sigma, reduce_data=None): # reduce_data=(idx_from, idx_to) if reduce_data is None: idx_from, idx_to = (0, data_size) else: idx_from, idx_to = reduce_data if(True): options = wt.WaveTrackerOptions( paths=True, frechet=True, cartesian=True, ) result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options) A = result.frechet.toarray() A = A[idx_from:idx_to] data_misfit_hess = 2 * A.T @ A / sigma**2 model_reg_hess = reg.hessian(slowness) return data_misfit_hess + model_reg_hess .. GENERATED FROM PYTHON SOURCE LINES 227-233 2. Invert with quadratic smoothing and damping regularization terms ------------------------------------------------------------------- 2.1 Define BaseProblem ~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 233-238 .. code-block:: Python # define CoFI BaseProblem fmm_problem_quadratic_reg = cofi.BaseProblem() fmm_problem_quadratic_reg.set_initial_model(ref_start_slowness.flatten()) .. GENERATED FROM PYTHON SOURCE LINES 240-249 .. code-block:: Python # add regularization: flattening + smoothing smoothing_factor = 5e6 reg_smoothing = smoothing_factor * cofi.utils.QuadraticReg( model_shape=model_shape, weighting_matrix="smoothing" ) reg = reg_smoothing .. GENERATED FROM PYTHON SOURCE LINES 251-257 .. code-block:: Python sigma = 0.000008 # data standard deviation of noise fmm_problem_quadratic_reg.set_objective(objective_func, args=[reg, sigma, None]) fmm_problem_quadratic_reg.set_gradient(gradient, args=[reg, sigma, None]) fmm_problem_quadratic_reg.set_hessian(hessian, args=[reg, sigma, None]) .. GENERATED FROM PYTHON SOURCE LINES 262-265 2.2 Define InversionOptions ~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 265-277 .. code-block:: Python my_options = cofi.InversionOptions() my_options.set_tool("cofi.simple_newton") my_options.set_params( num_iterations=15, step_length=1, obj_tol=1e-16, verbose=True, hessian_is_symmetric=True ) .. GENERATED FROM PYTHON SOURCE LINES 282-285 2.3 Start an inversion ~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 285-290 .. code-block:: Python inv = cofi.Inversion(fmm_problem_quadratic_reg, my_options) inv_result_quadratic_reg = inv.run() inv_result_quadratic_reg.summary() .. rst-class:: sphx-glr-script-out .. code-block:: none Iteration #0, updated objective function value: 37346.04444673342 Iteration #1, updated objective function value: 1496.5260409442917 Iteration #2, updated objective function value: 186.6540128296109 Iteration #3, updated objective function value: 49.72207315667848 Iteration #4, updated objective function value: 44.59552267851312 Iteration #5, updated objective function value: 87.36430757337139 Iteration #6, updated objective function value: 25.36158241794411 Iteration #7, updated objective function value: 36.87090297015299 Iteration #8, updated objective function value: 12.617214770325097 Iteration #9, updated objective function value: 13.857067302487787 Iteration #10, updated objective function value: 7.8961362057078714 Iteration #11, updated objective function value: 10.282043752198193 Iteration #12, updated objective function value: 14.08368705414384 Iteration #13, updated objective function value: 10.513448351701289 Iteration #14, updated objective function value: 5.3676081138584415 ============================ Summary for inversion result ============================ SUCCESS ---------------------------- model: [0.00037221 0.00037399 0.00037577 ... 0.00061149 0.00061299 0.00061453] num_iterations: 14 objective_val: 5.3676081138584415 n_obj_evaluations: 16 n_grad_evaluations: 15 n_hess_evaluations: 15 .. GENERATED FROM PYTHON SOURCE LINES 295-298 2.4 Plotting ~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 298-303 .. code-block:: Python vmodel_inverted = 1./inv_result_quadratic_reg.model.reshape(model_shape) wt.display_model(vmodel_inverted,extent=extent) # inverted model wt.display_model(good_model,extent=extent) # true model .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_002.png :alt: fmm tomography regularization discussion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_002.png :class: sphx-glr-multi-img * .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_003.png :alt: fmm tomography regularization discussion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_003.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 308-326 -------------- 3. Invert with Gaussian prior as regularization term ---------------------------------------------------- Instead of using a smoothing and damping regularization, in this section, we use a model covariance matrix and prior model. :math:`\chi_{P}^{2}=\left(\mathbf{y} -\mathbf{f}(\mathbf{m})\right)^T C_d^{-1} \left(\mathbf{y} -\mathbf{f}(\mathbf{m})\right) + \left( \mathbf{m} - \mathbf{m}_p \right)^T C_p^{-1} \left( \mathbf{m} - \mathbf{m}_p \right)` :math:`\Delta \mathbf{m}= ({J}^T {C}_d^{-1} {J}+{C}_p^{-1})^{-1} ({J}^T{C}_d^{-1} (\mathbf{y}-\mathbf{f}(\mathbf{m}))+{C}_p^{-1}(\mathbf{m}_p-\mathbf{m}))` We can use CoFI’s utility module to help us generate a the Gaussian prior term. 3.1 Define BaseProblem ~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 326-331 .. code-block:: Python # define CoFI BaseProblem fmm_problem_gaussian_prior = cofi.BaseProblem() fmm_problem_gaussian_prior.set_initial_model(ref_start_slowness.flatten()) .. GENERATED FROM PYTHON SOURCE LINES 333-344 .. code-block:: Python # add regularization: Gaussian prior corrx = 3.0 corry = 3.0 sigma_slowness = 0.5**2 sigma_slowness = 2.5E-6 gaussian_prior = 0.01 * cofi.utils.GaussianPrior( model_covariance_inv=((corrx, corry), sigma_slowness), mean_model=ref_start_slowness.reshape(model_shape) ) .. GENERATED FROM PYTHON SOURCE LINES 346-351 .. code-block:: Python fmm_problem_gaussian_prior.set_objective(objective_func, args=[gaussian_prior, sigma]) fmm_problem_gaussian_prior.set_gradient(gradient, args=[gaussian_prior, sigma]) fmm_problem_gaussian_prior.set_hessian(hessian, args=[gaussian_prior, sigma]) .. GENERATED FROM PYTHON SOURCE LINES 356-359 3.2 Start an inversion ~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 359-365 .. code-block:: Python # reuse the previously defined InversionOptions object inv = cofi.Inversion(fmm_problem_gaussian_prior, my_options) inv_result_gaussian_prior = inv.run() inv_result_gaussian_prior.summary() .. rst-class:: sphx-glr-script-out .. code-block:: none Iteration #0, updated objective function value: 2248.1364063962765 Iteration #1, updated objective function value: 92.79030072251801 Iteration #2, updated objective function value: 28.63098462661843 Iteration #3, updated objective function value: 26.08983730748967 Iteration #4, updated objective function value: 26.279799379783867 Iteration #5, updated objective function value: 26.071267115197394 Iteration #6, updated objective function value: 26.298470633668593 Iteration #7, updated objective function value: 26.079080995248482 Iteration #8, updated objective function value: 26.309599744775284 Iteration #9, updated objective function value: 26.07521806350919 Iteration #10, updated objective function value: 26.310385039780268 Iteration #11, updated objective function value: 26.074474718579328 Iteration #12, updated objective function value: 26.31244771548592 Iteration #13, updated objective function value: 26.075721652167612 Iteration #14, updated objective function value: 26.31302076319902 ============================ Summary for inversion result ============================ SUCCESS ---------------------------- model: [0.00049594 0.00049449 0.00049261 ... 0.00050232 0.00050181 0.00050139] num_iterations: 14 objective_val: 26.31302076319902 n_obj_evaluations: 16 n_grad_evaluations: 15 n_hess_evaluations: 15 .. GENERATED FROM PYTHON SOURCE LINES 370-373 3.3 Plotting ~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 373-378 .. code-block:: Python vmodel_inverted = 1./inv_result_gaussian_prior.model.reshape(model_shape) wt.display_model(vmodel_inverted,extent=extent) # inverted model wt.display_model(good_model,extent=extent) # true model .. rst-class:: sphx-glr-horizontal * .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_004.png :alt: fmm tomography regularization discussion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_004.png :class: sphx-glr-multi-img * .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_005.png :alt: fmm tomography regularization discussion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_005.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 383-388 4. L-curve ---------- Now we plot an L-curve for the smoothing regularization case. .. GENERATED FROM PYTHON SOURCE LINES 388-426 .. code-block:: Python lambdas = np.logspace(-4, 4, 10) my_lcurve_problems = [] for lamb in lambdas: my_reg = lamb * reg_smoothing my_problem = cofi.BaseProblem() my_problem.set_objective(objective_func, args=[my_reg, sigma]) my_problem.set_gradient(gradient, args=[my_reg, sigma]) my_problem.set_hessian(hessian, args=[my_reg, sigma]) my_problem.set_initial_model(ref_start_slowness.flatten()) my_lcurve_problems.append(my_problem) my_options.set_params(verbose=False) def my_callback(inv_result, i): m = inv_result.model slowness=m options = wt.WaveTrackerOptions( cartesian=True, ) result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options) # track wavefronts ttimes = result.ttimes res_norm = np.linalg.norm(ttimes - obstimes)/sigma**2 reg_norm = np.sqrt(reg_smoothing(m)) print(f"Finished inversion with lambda={lambdas[i]}: {res_norm}, {reg_norm}") return res_norm, reg_norm my_inversion_pool = cofi.utils.InversionPool( my_lcurve_problems, my_options, my_callback, False ) all_res, all_cb_returns = my_inversion_pool.run() l_curve_points = list(zip(*all_cb_returns)) .. rst-class:: sphx-glr-script-out .. code-block:: none Finished inversion with lambda=0.0001: 211500.44428776175, 0.13483254428280728 Finished inversion with lambda=0.000774263682681127: 260731.1744337502, 0.13496174203337225 Finished inversion with lambda=0.005994842503189409: 271386.6695804191, 0.13502504598278992 Finished inversion with lambda=0.046415888336127774: 451370.3900748823, 0.1327427726565774 Finished inversion with lambda=0.3593813663804626: 453312.5962048154, 0.13102864469172437 Finished inversion with lambda=2.782559402207126: 264044.6349749823, 0.13441418768920918 Finished inversion with lambda=21.54434690031882: 618549.3299736021, 0.13133502515213039 Finished inversion with lambda=166.81005372000558: 723609.8417236926, 0.12674068778075587 Finished inversion with lambda=1291.5496650148827: 316054.01961665374, 0.11347506334192359 Finished inversion with lambda=10000.0: 503616.0297332419, 0.09505699526967784 .. GENERATED FROM PYTHON SOURCE LINES 428-448 .. code-block:: Python # plot the L-curve res_norm, reg_norm = l_curve_points plt.plot(reg_norm, res_norm, '.-') plt.xlabel(r'Norm of regularization term $||Wm||_2$') plt.ylabel(r'Norm of residual $||g(m)-d||_2$') for i in range(len(lambdas)): plt.annotate(f'{lambdas[i]:.1e}', (reg_norm[i], res_norm[i]), fontsize=8) # plot the previously solved model my_inverted_model = inv_result_quadratic_reg.model my_reg_norm = np.sqrt(reg_smoothing(my_inverted_model)) slowness=my_inverted_model options = wt.WaveTrackerOptions(cartesian=True) result = wt.calc_wavefronts(1./slowness.reshape(model_shape),receivers,sources,extent=extent,options=options) # track wavefronts ttimes = result.ttimes my_residual_norm = np.linalg.norm(ttimes - obstimes)/sigma**2 plt.plot(my_reg_norm, my_residual_norm, "x") plt.annotate(f"{smoothing_factor:.1e}", (my_reg_norm, my_residual_norm), fontsize=8); .. image-sg:: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_006.png :alt: fmm tomography regularization discussion :srcset: /examples/generated/scripts_synth_data/images/sphx_glr_fmm_tomography_regularization_discussion_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.1341853290376085, 289115.0941701654, '5.0e+06') .. GENERATED FROM PYTHON SOURCE LINES 453-466 -------------- Watermark --------- .. raw:: html .. raw:: html .. GENERATED FROM PYTHON SOURCE LINES 466-472 .. code-block:: Python watermark_list = ["cofi", "numpy", "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+71.gb28b5b0 numpy 2.2.6 matplotlib 3.10.9 .. GENERATED FROM PYTHON SOURCE LINES 483-483 sphinx_gallery_thumbnail_number = -1 .. rst-class:: sphx-glr-timing **Total running time of the script:** (2 minutes 55.467 seconds) .. _sphx_glr_download_examples_generated_scripts_synth_data_fmm_tomography_regularization_discussion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: fmm_tomography_regularization_discussion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: fmm_tomography_regularization_discussion.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: fmm_tomography_regularization_discussion.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_