.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/plot_seismic_regularization.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_gallery_plot_seismic_regularization.py: 5. Seismic Regularization ========================= This example shows how to use the Curvelet transform to condition a missing-data seismic regularization problem. .. GENERATED FROM PYTHON SOURCE LINES 7-9 .. code-block:: Python # sphinx_gallery_thumbnail_number = 2 .. GENERATED FROM PYTHON SOURCE LINES 10-25 .. code-block:: Python import warnings warnings.filterwarnings("ignore") import matplotlib.pyplot as plt import numpy as np import pylops from pylops.optimization.sparsity import fista from scipy.signal import convolve from curvelops import FDCT2D np.random.seed(0) warnings.filterwarnings("ignore") .. GENERATED FROM PYTHON SOURCE LINES 26-28 Setup ===== .. GENERATED FROM PYTHON SOURCE LINES 28-64 .. code-block:: Python inputfile = "../testdata/seismic.npz" inputdata = np.load(inputfile) x = inputdata["R"][50, :, ::2] x = x / np.abs(x).max() taxis, xaxis = inputdata["t"][::2], inputdata["r"][0] par = {} par["nx"], par["nt"] = x.shape par["dx"] = inputdata["r"][0, 1] - inputdata["r"][0, 0] par["dt"] = inputdata["t"][1] - inputdata["t"][0] # Add wavelet wav = inputdata["wav"][::2] wav_c = np.argmax(wav) x = np.apply_along_axis(convolve, 1, x, wav, mode="full") x = x[:, wav_c:][:, : par["nt"]] # Gain gain = np.tile((taxis**2)[:, np.newaxis], (1, par["nx"])).T x *= gain # Subsampling locations perc_subsampling = 0.5 Nsub = int(np.round(par["nx"] * perc_subsampling)) iava = np.sort(np.random.permutation(np.arange(par["nx"]))[:Nsub]) # Restriction operator Rop = pylops.Restriction((par["nx"], par["nt"]), iava, axis=0, dtype="float64") y = Rop @ x xadj = Rop.H @ y # Apply mask ymask = Rop.mask(x) .. GENERATED FROM PYTHON SOURCE LINES 65-67 Curvelet transform ================== .. GENERATED FROM PYTHON SOURCE LINES 69-74 .. code-block:: Python DCTOp = FDCT2D((par["nx"], par["nt"]), nbscales=4) yc = DCTOp @ x xcadj = DCTOp.H @ yc .. GENERATED FROM PYTHON SOURCE LINES 75-90 .. code-block:: Python opts_plot = dict( cmap="gray", vmin=-0.1, vmax=0.1, extent=(xaxis[0], xaxis[-1], taxis[-1], taxis[0]), ) fig, axs = plt.subplots(1, 2, sharey=True, figsize=(10, 7)) axs[0].imshow(x.T, **opts_plot) axs[0].set_title("Data") axs[0].axis("tight") axs[1].imshow(np.real(xcadj).T, **opts_plot) axs[1].set_title("Adjoint curvelet") axs[1].axis("tight") .. image-sg:: /gallery/images/sphx_glr_plot_seismic_regularization_001.png :alt: Data, Adjoint curvelet :srcset: /gallery/images/sphx_glr_plot_seismic_regularization_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (0.0, 3000.0, 1.995, 0.0) .. GENERATED FROM PYTHON SOURCE LINES 91-93 Reconstruction based on Curvelet transform ########################################## .. GENERATED FROM PYTHON SOURCE LINES 95-96 Combined modelling operator .. GENERATED FROM PYTHON SOURCE LINES 96-104 .. code-block:: Python RCop = Rop @ DCTOp.H RCop.dims = (RCop.shape[1],) # flatten RCop.dimsd = (RCop.shape[0],) # Inverse pl1, _, cost = fista(RCop, y.ravel(), niter=100, eps=1e-3, show=True) xl1 = (DCTOp.H @ pl1).real.reshape(x.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none FISTA (soft thresholding) -------------------------------------------------------------------------------- The Operator Op has 20000 rows and 305683 cols eps = 1.000000e-03 tol = 1.000000e-10 niter = 100 alpha = 1.000000e+00 thresh = 5.000000e-04 -------------------------------------------------------------------------------- Itn x[0] r2norm r12norm xupdate 1 -7.50e-04+7.58e-19j 1.067e-02 2.690e-01 1.383e+00 2 -9.12e-04+5.16e-19j 9.923e-03 2.541e-01 8.676e-02 3 -1.04e-03+1.27e-19j 9.270e-03 2.408e-01 9.289e-02 4 -1.17e-03-1.43e-19j 8.703e-03 2.293e-01 9.611e-02 5 -1.24e-03-1.52e-19j 8.240e-03 2.191e-01 9.767e-02 6 -1.25e-03-1.53e-19j 7.846e-03 2.104e-01 9.769e-02 7 -1.24e-03-1.52e-19j 7.534e-03 2.028e-01 9.678e-02 8 -1.23e-03-1.50e-19j 7.286e-03 1.963e-01 9.516e-02 9 -1.21e-03-1.49e-19j 7.091e-03 1.906e-01 9.344e-02 10 -1.21e-03-1.48e-19j 6.910e-03 1.856e-01 9.125e-02 11 -1.21e-03-1.49e-19j 6.755e-03 1.814e-01 8.857e-02 21 -1.24e-03-1.52e-19j 6.198e-03 1.592e-01 6.407e-02 31 -1.06e-03+1.30e-19j 6.083e-03 1.522e-01 4.776e-02 41 -1.05e-03+1.29e-19j 6.050e-03 1.493e-01 3.752e-02 51 -1.12e-03+1.37e-19j 6.022e-03 1.478e-01 3.129e-02 61 -1.06e-03+1.29e-19j 6.033e-03 1.469e-01 2.793e-02 71 -1.06e-03-1.30e-19j 6.018e-03 1.463e-01 2.484e-02 81 -1.05e-03-1.29e-19j 6.023e-03 1.459e-01 2.252e-02 91 -1.04e-03+1.28e-19j 6.031e-03 1.456e-01 1.938e-02 92 -1.04e-03+1.28e-19j 6.034e-03 1.456e-01 1.921e-02 93 -1.04e-03+1.28e-19j 6.034e-03 1.456e-01 1.903e-02 94 -1.04e-03+1.27e-19j 6.033e-03 1.456e-01 1.884e-02 95 -1.04e-03+1.27e-19j 6.033e-03 1.455e-01 1.867e-02 96 -1.04e-03+1.27e-19j 6.033e-03 1.455e-01 1.851e-02 97 -1.04e-03+1.27e-19j 6.033e-03 1.455e-01 1.830e-02 98 -1.04e-03+1.27e-19j 6.029e-03 1.455e-01 1.802e-02 99 -1.04e-03+1.27e-19j 6.026e-03 1.455e-01 1.775e-02 100 -1.03e-03+1.27e-19j 6.024e-03 1.455e-01 1.753e-02 Iterations = 100 Total time (s) = 18.18 -------------------------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 105-119 .. code-block:: Python fig, axs = plt.subplots(1, 4, sharey=True, figsize=(16, 7)) axs[0].imshow(x.T, **opts_plot) axs[0].set_title("Data") axs[0].axis("tight") axs[1].imshow(ymask.T, **opts_plot) axs[1].set_title("Masked data") axs[1].axis("tight") axs[2].imshow(xl1.T, **opts_plot) axs[2].set_title("Reconstructed data") axs[2].axis("tight") axs[3].imshow((x - xl1).T, **opts_plot) axs[3].set_title("Reconstruction error") axs[3].axis("tight") .. image-sg:: /gallery/images/sphx_glr_plot_seismic_regularization_002.png :alt: Data, Masked data, Reconstructed data, Reconstruction error :srcset: /gallery/images/sphx_glr_plot_seismic_regularization_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (0.0, 3000.0, 1.995, 0.0) .. GENERATED FROM PYTHON SOURCE LINES 120-124 .. code-block:: Python fig, ax = plt.subplots(figsize=(16, 2)) ax.plot(range(1, len(cost) + 1), cost, "k") ax.set(xlim=[1, len(cost)]) fig.suptitle("FISTA convergence") .. image-sg:: /gallery/images/sphx_glr_plot_seismic_regularization_003.png :alt: FISTA convergence :srcset: /gallery/images/sphx_glr_plot_seismic_regularization_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 0.98, 'FISTA convergence') .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 19.310 seconds) .. _sphx_glr_download_gallery_plot_seismic_regularization.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_seismic_regularization.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_seismic_regularization.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_seismic_regularization.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_