.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/0003.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_0003.py: .. _gallery:0003: Cosine+Spikes in Dirac-Cosine Basis ====================================== .. contents:: :depth: 2 :local: The Cosine+Spikes signal in this example consists of a mixture of 3 different cosine waves with different amplitudes and 120 different spikes (with normally distributed amplitudes and randomly chosen locations). The Cosine part of the signal has a sparse representation in the Cosine (DCT) basis. The Spikes part of the signal has a sparse representation in the Dirac(Identity/Standard) basis. Thus, the mixture Cosine+Spikes signal has a sparse representation (of 123 nonzero coefficients) in the Dirac-Cosine two ortho basis. Note that the spikes are normally distributed. Some of the spikes have extremely low amplitudes. They may be missed by a recovery algorithm depending on convergence thresholds. See also: * :ref:`api:problems` * :ref:`api:lop` .. GENERATED FROM PYTHON SOURCE LINES 35-46 .. code-block:: default # Configure JAX to work with 64-bit floating point precision. from jax.config import config config.update("jax_enable_x64", True) import jax.numpy as jnp import cr.nimble as crn import cr.sparse.plots as crplot .. GENERATED FROM PYTHON SOURCE LINES 47-51 Setup ------------------------------ We shall construct our test signal and dictionary using our test problems module. .. GENERATED FROM PYTHON SOURCE LINES 51-57 .. code-block:: default from cr.sparse import problems prob = problems.generate('cosine-spikes:dirac-dct', c=3) fig, ax = problems.plot(prob) .. image-sg:: /gallery/images/sphx_glr_0003_001.png :alt: Cosine with spikes, Cosine part of signal, Spikes part of signal, Dirac DCT Model :srcset: /gallery/images/sphx_glr_0003_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 58-59 Let us access the relevant parts of our test problem .. GENERATED FROM PYTHON SOURCE LINES 59-68 .. code-block:: default # The sparsifying basis linear operator A = prob.A # The Cosine+Spikes signal b0 = prob.b # The sparse representation of the signal in the dictionary x0 = prob.x .. GENERATED FROM PYTHON SOURCE LINES 69-71 Check how many coefficients in the sparse representation are sufficient to capture 99.9% of the energy of the signal .. GENERATED FROM PYTHON SOURCE LINES 71-73 .. code-block:: default print(crn.num_largest_coeffs_for_energy_percent(x0, 99.9)) .. rst-class:: sphx-glr-script-out .. code-block:: none 78 .. GENERATED FROM PYTHON SOURCE LINES 74-76 This number gives us an idea about the required sparsity to be configured for greedy pursuit algorithms. .. GENERATED FROM PYTHON SOURCE LINES 78-81 Sparse Recovery using Subspace Pursuit ------------------------------------------- We shall use subspace pursuit to reconstruct the signal. .. GENERATED FROM PYTHON SOURCE LINES 81-85 .. code-block:: default import cr.sparse.pursuit.sp as sp # We will first try to estimate a 100-sparse representation sol = sp.solve(A, b0, 100) .. GENERATED FROM PYTHON SOURCE LINES 86-87 This utility function helps us quickly analyze the quality of reconstruction .. GENERATED FROM PYTHON SOURCE LINES 87-88 .. code-block:: default problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 50.698 reconstruction: 50.695 SNR: 38.38 dB x_norm: original: 50.704 reconstruction: 50.673 SNR: 38.36 dB Sparsity: original: 78, reconstructed: 76, overlap: 76, ratio: 0.974 Iterations: 20 .. GENERATED FROM PYTHON SOURCE LINES 89-91 It takes 20 iterations to converge and 76 of the largest 78 entries have been correctly identified. .. GENERATED FROM PYTHON SOURCE LINES 93-94 We will now try to estimate a 150-sparse representation .. GENERATED FROM PYTHON SOURCE LINES 94-96 .. code-block:: default sol = sp.solve(A, b0, 150) problems.analyze_solution(prob, sol) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 50.698 reconstruction: 50.698 SNR: 294.28 dB x_norm: original: 50.704 reconstruction: 50.704 SNR: 294.33 dB Sparsity: original: 78, reconstructed: 78, overlap: 78, ratio: 1.000 Iterations: 5 .. GENERATED FROM PYTHON SOURCE LINES 97-100 We have correctly detected all the 78 most significant entries Let us check if we correctly decoded all the nonzero entries in the sparse representation x .. GENERATED FROM PYTHON SOURCE LINES 100-102 .. code-block:: default problems.analyze_solution(prob, sol, perc=100) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 50.698 reconstruction: 50.698 SNR: 294.28 dB x_norm: original: 50.704 reconstruction: 50.704 SNR: 294.33 dB Sparsity: original: 123, reconstructed: 123, overlap: 123, ratio: 1.000 Iterations: 5 .. GENERATED FROM PYTHON SOURCE LINES 103-104 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 104-105 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 106-107 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 107-109 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 110-111 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 111-117 .. code-block:: default ax = crplot.h_plots(2) ax[0].stem(x0, markerfmt='.') ax[1].stem(x, markerfmt='.') .. image-sg:: /gallery/images/sphx_glr_0003_002.png :alt: 0003 :srcset: /gallery/images/sphx_glr_0003_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 118-119 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 119-125 .. code-block:: default ax = crplot.h_plots(2) ax[0].plot(b0) ax[1].plot(b) .. image-sg:: /gallery/images/sphx_glr_0003_003.png :alt: 0003 :srcset: /gallery/images/sphx_glr_0003_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 126-129 Sparse Recovery using Compressive Sampling Matching Pursuit --------------------------------------------------------------- We shall now use compressive sampling matching pursuit to reconstruct the signal. .. GENERATED FROM PYTHON SOURCE LINES 129-134 .. code-block:: default import cr.sparse.pursuit.cosamp as cosamp # We will try to estimate a 150-sparse representation sol = cosamp.solve(A, b0, 150) problems.analyze_solution(prob, sol, perc=100) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 50.698 reconstruction: 50.698 SNR: 296.28 dB x_norm: original: 50.704 reconstruction: 50.704 SNR: 295.75 dB Sparsity: original: 123, reconstructed: 123, overlap: 123, ratio: 1.000 Iterations: 6 .. GENERATED FROM PYTHON SOURCE LINES 135-136 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 136-137 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 138-139 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 139-142 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 143-144 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 144-149 .. code-block:: default ax = crplot.h_plots(2) ax[0].stem(x0, markerfmt='.') ax[1].stem(x, markerfmt='.') .. image-sg:: /gallery/images/sphx_glr_0003_004.png :alt: 0003 :srcset: /gallery/images/sphx_glr_0003_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 150-151 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 151-156 .. code-block:: default ax = crplot.h_plots(2) ax[0].plot(b0) ax[1].plot(b) .. image-sg:: /gallery/images/sphx_glr_0003_005.png :alt: 0003 :srcset: /gallery/images/sphx_glr_0003_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 157-159 Sparse Recovery using SPGL1 --------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 159-164 .. code-block:: default import cr.sparse.cvx.spgl1 as crspgl1 options = crspgl1.SPGL1Options() sol = crspgl1.solve_bp_jit(A, b0, options=options) problems.analyze_solution(prob, sol, perc=100) .. rst-class:: sphx-glr-script-out .. code-block:: none m: 1024, n: 2048 b_norm: original: 50.698 reconstruction: 50.698 SNR: 123.14 dB x_norm: original: 50.704 reconstruction: 50.704 SNR: 115.88 dB Sparsity: original: 123, reconstructed: 123, overlap: 123, ratio: 1.000 Iterations: 55 n_times: 58, n_trans: 56 .. GENERATED FROM PYTHON SOURCE LINES 165-166 The estimated sparse representation .. GENERATED FROM PYTHON SOURCE LINES 166-167 .. code-block:: default x = sol.x .. GENERATED FROM PYTHON SOURCE LINES 168-169 Let us reconstruct the signal from this sparse representation .. GENERATED FROM PYTHON SOURCE LINES 169-172 .. code-block:: default b = prob.reconstruct(x) .. GENERATED FROM PYTHON SOURCE LINES 173-174 Let us visualize the original and reconstructed representation .. GENERATED FROM PYTHON SOURCE LINES 174-179 .. code-block:: default ax = crplot.h_plots(2) ax[0].stem(x0, markerfmt='.') ax[1].stem(x, markerfmt='.') .. image-sg:: /gallery/images/sphx_glr_0003_006.png :alt: 0003 :srcset: /gallery/images/sphx_glr_0003_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 180-181 Let us visualize the original and reconstructed signal .. GENERATED FROM PYTHON SOURCE LINES 181-186 .. code-block:: default ax = crplot.h_plots(2) ax[0].plot(b0) ax[1].plot(b) .. image-sg:: /gallery/images/sphx_glr_0003_007.png :alt: 0003 :srcset: /gallery/images/sphx_glr_0003_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 187-196 Comments --------- * Both SP and CoSaMP correctly recover the signal in 5 iterations if the sparsity is specified properly. * SP recovery is slightly inaccurate if the sparsity is incorrectly specified. It also takes more iterations to converge. * SPGL1 converges in 55 iterations but correctly discovers the support. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 10.973 seconds) .. _sphx_glr_download_gallery_0003.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 0003.py <0003.py>` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 0003.ipynb <0003.ipynb>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_