.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "galley_examples/relaxation/plot-5Cs-95Si.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_galley_examples_relaxation_plot-5Cs-95Si.py: 0.05 Cs2O • 0.95 SiO2 MAS-ETA ============================= .. GENERATED FROM PYTHON SOURCE LINES 7-11 The following example is an application of the statistical learning method in determining the distribution of the Si-29 echo train decay constants in glasses. Import all relevant packages. .. GENERATED FROM PYTHON SOURCE LINES 11-31 .. code-block:: Python import csdmpy as cp import matplotlib.pyplot as plt import numpy as np from mrinversion.kernel import relaxation from mrinversion.linear_model import LassoFistaCV, TSVDCompression from csdmpy import statistics as stats plt.rcParams["pdf.fonttype"] = 42 # For using plots in Illustrator plt.rc("font", size=9) def plot2D(csdm_object, **kwargs): plt.figure(figsize=(4, 3)) csdm_object.plot(cmap="gist_ncar_r", **kwargs) plt.tight_layout() plt.show() .. GENERATED FROM PYTHON SOURCE LINES 33-38 Dataset setup ------------- Import the dataset '''''''''''''''''' Load the dataset as a CSDM data-object. .. GENERATED FROM PYTHON SOURCE LINES 38-52 .. code-block:: Python # The 2D SE-PIETA-MAS dataset in csdm format domain = "https://www.ssnmr.org/sites/default/files/mrsimulator" filename = f"{domain}/MAS_SE_PIETA_5Cs_95Si_FT.csdf" data_object = cp.load(filename) # Inversion only requires the real part of the complex dataset. data_object = data_object.real sigma = 1407.443 # data standard deviation # Convert the MAS dimension from Hz to ppm. data_object.dimensions[0].to("ppm", "nmr_frequency_ratio") plot2D(data_object) .. image-sg:: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_001.png :alt: plot 5Cs 95Si :srcset: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 53-55 Prepping the data for inversion ''''''''''''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 55-59 .. code-block:: Python data_object = data_object.T data_object_truncated = data_object[:, 1250:-1250] plot2D(data_object_truncated) .. image-sg:: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_002.png :alt: plot 5Cs 95Si :srcset: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 60-64 Linear Inversion setup ---------------------- Dimension setup ''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 64-67 .. code-block:: Python data_object_truncated.dimensions[0].to("s") # set coordinates to 's' kernel_dimension = data_object_truncated.dimensions[0] .. GENERATED FROM PYTHON SOURCE LINES 68-70 Generating the kernel ''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 70-83 .. code-block:: Python relaxT2 = relaxation.T2( kernel_dimension=kernel_dimension, inverse_dimension=dict( count=32, minimum="1e-3 s", maximum="1e4 s", scale="log", label=r"log ($\lambda^{-1}$ / s)", ), ) inverse_dimension = relaxT2.inverse_dimension K = relaxT2.kernel(supersampling=20) .. GENERATED FROM PYTHON SOURCE LINES 84-86 Data Compression '''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 86-92 .. code-block:: Python new_system = TSVDCompression(K, data_object_truncated) compressed_K = new_system.compressed_K compressed_s = new_system.compressed_s print(f"truncation_index = {new_system.truncation_index}") .. rst-class:: sphx-glr-script-out .. code-block:: none compression factor = 1.3333333333333333 truncation_index = 18 .. GENERATED FROM PYTHON SOURCE LINES 93-97 Solving the inverse problem --------------------------- FISTA LASSO cross-validation ''''''''''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 97-112 .. code-block:: Python # setup the pre-defined range of alpha and lambda values lambdas = 10 ** (-4 + 5 * (np.arange(32) / 31)) # setup the smooth lasso cross-validation class s_lasso = LassoFistaCV( lambdas=lambdas, # A numpy array of lambda values. sigma=sigma, # data standard deviation folds=5, # The number of folds in n-folds cross-validation. inverse_dimension=inverse_dimension, # previously defined inverse dimensions. ) # run the fit method on the compressed kernel and compressed data. s_lasso.fit(K=compressed_K, s=compressed_s) .. GENERATED FROM PYTHON SOURCE LINES 113-115 The optimum hyper-parameters '''''''''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 115-117 .. code-block:: Python print(s_lasso.hyperparameters) .. rst-class:: sphx-glr-script-out .. code-block:: none {'lambda': np.float64(0.012496091412919868)} .. GENERATED FROM PYTHON SOURCE LINES 118-120 The cross-validation curve '''''''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 120-125 .. code-block:: Python plt.figure(figsize=(4, 3)) s_lasso.cv_plot() plt.tight_layout() plt.show() .. image-sg:: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_003.png :alt: plot 5Cs 95Si :srcset: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 126-128 The optimum solution '''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 128-144 .. code-block:: Python f_sol = s_lasso.f levels = np.arange(15) / 15 + 0.1 plt.figure(figsize=(3.85, 2.75)) # set the figure size ax = plt.subplot(projection="csdm") cb = ax.contourf(f_sol / f_sol.max(), levels=levels, cmap="jet_r") ax.set_ylim(-70, -130) ax.set_xlim(-3, 2.5) plt.title("5Cs:95Si") ax.set_xlabel(r"$\log(\lambda^{-1}\,/\,$s)") ax.set_ylabel("Frequency / ppm") plt.grid(linestyle="--", alpha=0.75) plt.colorbar(cb, ticks=np.arange(11) / 10) plt.tight_layout() plt.show() .. image-sg:: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_004.png :alt: 5Cs:95Si :srcset: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 145-147 The fit residuals ''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 147-150 .. code-block:: Python residuals = s_lasso.residuals(K=K, s=data_object_truncated) plot2D(residuals) .. image-sg:: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_005.png :alt: plot 5Cs 95Si :srcset: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 151-152 The standard deviation of the residuals is .. GENERATED FROM PYTHON SOURCE LINES 152-154 .. code-block:: Python residuals.std() .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 155-157 Saving the solution ''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 157-160 .. code-block:: Python f_sol.save("5Cs-95Si_inverse.csdf") # save the solution residuals.save("5Cs-95Si-residue.csdf") # save the residuals .. GENERATED FROM PYTHON SOURCE LINES 161-163 Analysis -------- .. GENERATED FROM PYTHON SOURCE LINES 163-176 .. code-block:: Python # Normalize the distribution to 1. f_sol /= f_sol.max() # Get the Q4 and Q3 cross-sections. Q4_coordinate = -110.1e-6 # ppm Q3_coordinate = -99.4e-6 # ppm Q4_index = np.where(f_sol.dimensions[1].coordinates >= Q4_coordinate)[0][0] Q3_index = np.where(f_sol.dimensions[1].coordinates >= Q3_coordinate)[0][0] Q4_region = f_sol[:, Q4_index] Q3_region = f_sol[:, Q3_index] .. GENERATED FROM PYTHON SOURCE LINES 177-178 Plot of the Q4 and Q3 cross-sections .. GENERATED FROM PYTHON SOURCE LINES 178-200 .. code-block:: Python fig, ax = plt.subplots(1, 2, figsize=(7, 2.75), subplot_kw={"projection": "csdm"}) cb = ax[0].contourf(f_sol, levels=levels, cmap="jet_r") ax[0].arrow(1, Q4_coordinate * 1e6, -0.5, 0, color="blue") ax[0].arrow(1, Q3_coordinate * 1e6, -0.5, 0, color="orange") ax[0].set_ylim(-70, -130) ax[0].set_xlim(-3, 2.5) ax[0].set_xlabel(r"$\log(\lambda^{-1}\,/\,$s)") ax[0].set_ylabel("Frequency / ppm") ax[0].grid(linestyle="--", alpha=0.75) ax[1].plot(Q4_region, label="Q4") ax[1].plot(Q3_region, label="Q3") ax[1].set_xlim(-3, 2.5) ax[1].set_xlabel(r"$\log(\lambda^{-1}\,/\,$s)") ax[1].grid(linestyle="--", alpha=0.75) plt.colorbar(cb, ax=ax[0], ticks=np.arange(11) / 10) plt.tight_layout() plt.legend() plt.savefig("5Cs-95Si.pdf") plt.show() .. image-sg:: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_006.png :alt: plot 5Cs 95Si :srcset: /galley_examples/relaxation/images/sphx_glr_plot-5Cs-95Si_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 201-206 Mean and mode analysis '''''''''''''''''''''' The T2 distribution is sampled over a log-linear scale. The statistical mean of the `Q4_region` and `Q3_region` in log10(T2). The mean T2 is 10**(log10(T2)), in units of seconds. .. GENERATED FROM PYTHON SOURCE LINES 206-210 .. code-block:: Python Q4_mean = 10 ** stats.mean(Q4_region)[0] * 1e3 # ms Q3_mean = 10 ** stats.mean(Q3_region)[0] * 1e3 # ms .. GENERATED FROM PYTHON SOURCE LINES 211-212 Mode the argument corresponding to the max distribution. .. GENERATED FROM PYTHON SOURCE LINES 212-225 .. code-block:: Python # index corresponding to the max distribution. arg_index_Q4 = int(np.argmax(Q4_region)) arg_index_Q3 = int(np.argmax(Q3_region)) # log10(T2) coordinates corresponding to the max distribution. arg_coord_Q4 = Q4_region.dimensions[0].coordinates[arg_index_Q4] arg_coord_Q3 = Q3_region.dimensions[0].coordinates[arg_index_Q3] # T2 coordinates corresponding to the max distribution. Q4_mode = 10**arg_coord_Q4 * 1e3 # ms Q3_mode = 10**arg_coord_Q3 * 1e3 # ms .. GENERATED FROM PYTHON SOURCE LINES 226-228 Results ''''''' .. GENERATED FROM PYTHON SOURCE LINES 228-232 .. code-block:: Python print(f"Q4 statistics:\n\tmean = {Q4_mean} ms,\n\tmode = {Q4_mode} ms\n") print(f"Q3 statistics:\n\tmean = {Q3_mean} ms,\n\tmode = {Q3_mode} ms\n") print(f"r_λ (mean) = {Q4_mean/Q3_mean}") print(f"r_λ (mode) = {Q4_mode/Q3_mode}") .. rst-class:: sphx-glr-script-out .. code-block:: none Q4 statistics: mean = 316.1468617571331 ms, mode = 512.4805876960905 ms Q3 statistics: mean = 192.10371720040615 ms, mode = 181.16091942004059 ms r_λ (mean) = 1.6457092364710613 r_λ (mode) = 2.828869434625967 .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.826 seconds) .. _sphx_glr_download_galley_examples_relaxation_plot-5Cs-95Si.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot-5Cs-95Si.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot-5Cs-95Si.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot-5Cs-95Si.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_