.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "galley_examples/MAF/plot_2d_5_CaO.SiO2.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_MAF_plot_2d_5_CaO.SiO2.py: 2D MAF data of CaO.SiO2 glass ============================= .. GENERATED FROM PYTHON SOURCE LINES 7-16 The following example illustrates an application of the statistical learning method applied in determining the distribution of the nuclear shielding tensor parameters from a 2D magic-angle flipping (MAF) spectrum. In this example, we use the 2D MAF spectrum [#f1]_ of :math:`\text{CaO}\cdot\text{SiO}_2` glass. Before getting started ---------------------- Import all relevant packages. .. GENERATED FROM PYTHON SOURCE LINES 16-25 .. code-block:: Python import csdmpy as cp import matplotlib.pyplot as plt import numpy as np from mrinversion.kernel.nmr import ShieldingPALineshape from mrinversion.linear_model import SmoothLasso, TSVDCompression from mrinversion.utils import plot_3d, to_Haeberlen_grid .. GENERATED FROM PYTHON SOURCE LINES 27-28 Setup for the matplotlib figures. .. GENERATED FROM PYTHON SOURCE LINES 28-41 .. code-block:: Python # function for plotting 2D dataset def plot2D(csdm_object, **kwargs): plt.figure(figsize=(4.5, 3.5)) ax = plt.subplot(projection="csdm") ax.imshow(csdm_object, cmap="gist_ncar_r", aspect="auto", **kwargs) ax.invert_xaxis() ax.invert_yaxis() plt.tight_layout() plt.show() .. GENERATED FROM PYTHON SOURCE LINES 42-49 Dataset setup ------------- Import the dataset '''''''''''''''''' Load the dataset. Here, we import the dataset as the CSDM data-object. .. GENERATED FROM PYTHON SOURCE LINES 49-60 .. code-block:: Python # The 2D MAF dataset in csdm format filename = "https://zenodo.org/record/3964531/files/CaO-SiO2-MAF.csdf" data_object = cp.load(filename) # For inversion, we only interest ourselves with the real part of the complex dataset. data_object = data_object.real # We will also convert the coordinates of both dimensions from Hz to ppm. _ = [item.to("ppm", "nmr_frequency_ratio") for item in data_object.dimensions] .. GENERATED FROM PYTHON SOURCE LINES 61-65 Here, the variable ``data_object`` is a `CSDM `_ object that holds the real part of the 2D MAF dataset. The plot of the 2D MAF dataset is .. GENERATED FROM PYTHON SOURCE LINES 65-67 .. code-block:: Python plot2D(data_object) .. image-sg:: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_001.png :alt: plot 2d 5 CaO.SiO2 :srcset: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 68-80 There are two dimensions in this dataset. The dimension at index 0 is the isotropic chemical shift dimension, whereas the dimension at index 1 is the pure anisotropic dimension. Prepping the data for inversion ''''''''''''''''''''''''''''''' **Step-1: Data Alignment** When using the csdm objects with the ``mrinversion`` package, the dimension at index 0 must be the dimension undergoing the linear inversion. In this example, we plan to invert the pure anisotropic shielding line-shape. In the ``data_object``, the anisotropic dimension is at index 1. Transpose the dataset before proceeding. .. GENERATED FROM PYTHON SOURCE LINES 80-82 .. code-block:: Python data_object = data_object.T .. GENERATED FROM PYTHON SOURCE LINES 83-89 **Step-2: Optimization** Also notice, the signal from the 2D MAF dataset occupies a small fraction of the two-dimensional frequency grid. For optimum performance, truncate the dataset to the relevant region before proceeding. Use the appropriate array indexing/slicing to select the signal region. .. GENERATED FROM PYTHON SOURCE LINES 89-92 .. code-block:: Python data_object_truncated = data_object[30:-30, 110:145] plot2D(data_object_truncated) .. image-sg:: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_002.png :alt: plot 2d 5 CaO.SiO2 :srcset: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 93-103 Linear Inversion setup ---------------------- Dimension setup ''''''''''''''' **Anisotropic-dimension:** The dimension of the dataset that holds the pure anisotropic frequency contributions. In ``mrinversion``, this must always be the dimension at index 0 of the data object. .. GENERATED FROM PYTHON SOURCE LINES 103-105 .. code-block:: Python anisotropic_dimension = data_object_truncated.dimensions[0] .. GENERATED FROM PYTHON SOURCE LINES 106-108 **x-y dimensions:** The two inverse dimensions corresponding to the `x` and `y`-axis of the `x`-`y` grid. .. GENERATED FROM PYTHON SOURCE LINES 108-113 .. code-block:: Python inverse_dimensions = [ cp.LinearDimension(count=25, increment="400 Hz", label="x"), # the `x`-dimension. cp.LinearDimension(count=25, increment="400 Hz", label="y"), # the `y`-dimension. ] .. GENERATED FROM PYTHON SOURCE LINES 114-121 Generating the kernel ''''''''''''''''''''' For MAF datasets, the line-shape kernel corresponds to the pure nuclear shielding anisotropy line-shapes. Use the :class:`~mrinversion.kernel.nmr.ShieldingPALineshape` class to generate a shielding line-shape kernel. .. GENERATED FROM PYTHON SOURCE LINES 121-131 .. code-block:: Python lineshape = ShieldingPALineshape( anisotropic_dimension=anisotropic_dimension, inverse_dimension=inverse_dimensions, channel="29Si", magnetic_flux_density="9.4 T", rotor_angle="90°", rotor_frequency="10.4 kHz", number_of_sidebands=4, ) .. GENERATED FROM PYTHON SOURCE LINES 132-148 Here, ``lineshape`` is an instance of the :class:`~mrinversion.kernel.nmr.ShieldingPALineshape` class. The required arguments of this class are the `anisotropic_dimension`, `inverse_dimension`, and `channel`. We have already defined the first two arguments in the previous sub-section. The value of the `channel` argument is the nucleus observed in the MAF experiment. In this example, this value is '29Si'. The remaining arguments, such as the `magnetic_flux_density`, `rotor_angle`, and `rotor_frequency`, are set to match the conditions under which the 2D MAF spectrum was acquired. The value of the `number_of_sidebands` argument is the number of sidebands calculated for each line-shape within the kernel. Unless, you have a lot of spinning sidebands in your MAF dataset, four sidebands should be enough. Once the ShieldingPALineshape instance is created, use the :meth:`~mrinversion.kernel.nmr.ShieldingPALineshape.kernel` method of the instance to generate the MAF line-shape kernel. .. GENERATED FROM PYTHON SOURCE LINES 148-151 .. code-block:: Python K = lineshape.kernel(supersampling=1) print(K.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none (68, 625) .. GENERATED FROM PYTHON SOURCE LINES 152-155 The kernel ``K`` is a NumPy array of shape (32, 784), where the axes with 32 and 784 points are the anisotropic dimension and the features (x-y coordinates) corresponding to the :math:`28\times 28` `x`-`y` grid, respectively. .. GENERATED FROM PYTHON SOURCE LINES 157-162 Data Compression '''''''''''''''' Data compression is optional but recommended. It may reduce the size of the inverse problem and, thus, further computation time. .. GENERATED FROM PYTHON SOURCE LINES 162-168 .. code-block:: Python new_system = TSVDCompression(K=K, s=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.5454545454545454 truncation_index = 44 .. GENERATED FROM PYTHON SOURCE LINES 169-182 Solving the inverse problem --------------------------- Smooth LASSO cross-validation ''''''''''''''''''''''''''''' Solve the smooth-lasso problem. Ordinarily, one should use the statistical learning method to solve the inverse problem over a range of α and λ values and then determine the best nuclear shielding tensor parameter distribution for the given 2D MAF dataset. Considering the limited build time for the documentation, we skip this step and evaluate the distribution at pre-optimized α and λ values. The optimum values are :math:`\alpha = 2.8\times 10^{-5}` and :math:`\lambda = 8.85\times 10^{-6}`. The following commented code was used in determining the optimum α and λ values. .. GENERATED FROM PYTHON SOURCE LINES 184-216 .. code-block:: Python # from mrinversion.linear_model import SmoothLassoCV # import numpy as np # # setup the pre-defined range of alpha and lambda values # lambdas = 10 ** (-4 - 2 * (np.arange(20) / 19)) # alphas = 10 ** (-3.5 - 2 * (np.arange(20) / 19)) # # setup the smooth lasso cross-validation class # s_lasso = SmoothLassoCV( # alphas=alphas, # A numpy array of alpha values. # lambdas=lambdas, # A numpy array of lambda values. # sigma=0.0012, # The standard deviation of noise from the MAF data. # folds=10, # The number of folds in n-folds cross-validation. # inverse_dimension=inverse_dimensions, # previously defined inverse dimensions. # verbose=1, # If non-zero, prints the progress as the computation proceeds. # max_iterations=20000, # maximum number of allowed iterations. # ) # # run fit using the compressed kernel and compressed data. # s_lasso.fit(compressed_K, compressed_s) # # the optimum hyper-parameters, alpha and lambda, from the cross-validation. # print(s_lasso.hyperparameters) # # {'alpha': 3.359818286283781e-05, 'lambda': 5.324953129837531e-06} # # the solution # f_sol = s_lasso.f # # the cross-validation error curve # CV_metric = s_lasso.cross_validation_curve .. GENERATED FROM PYTHON SOURCE LINES 217-218 If you use the above ``SmoothLassoCV`` method, skip the following code-block. .. GENERATED FROM PYTHON SOURCE LINES 218-225 .. code-block:: Python s_lasso = SmoothLasso( alpha=2.8e-5, lambda1=8.85e-6, inverse_dimension=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 226-231 The optimum solution '''''''''''''''''''' The :attr:`~mrinversion.linear_model.SmoothLasso.f` attribute of the instance holds the solution, .. GENERATED FROM PYTHON SOURCE LINES 231-233 .. code-block:: Python f_sol = s_lasso.f # f_sol is a CSDM object. .. GENERATED FROM PYTHON SOURCE LINES 234-241 where ``f_sol`` is the optimum solution. The fit residuals ''''''''''''''''' To calculate the residuals between the data and predicted data(fit), use the :meth:`~mrinversion.linear_model.SmoothLasso.residuals` method, as follows, .. GENERATED FROM PYTHON SOURCE LINES 241-247 .. code-block:: Python residuals = s_lasso.residuals(K, data_object_truncated) # residuals is a CSDM object. # The plot of the residuals. plot2D(residuals, vmax=data_object_truncated.max(), vmin=data_object_truncated.min()) .. image-sg:: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_003.png :alt: plot 2d 5 CaO.SiO2 :srcset: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 248-249 The standard deviation of the residuals is .. GENERATED FROM PYTHON SOURCE LINES 249-251 .. code-block:: Python residuals.std() .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 252-257 Saving the solution ''''''''''''''''''' To serialize the solution to a file, use the `save()` method of the CSDM object, for example, .. GENERATED FROM PYTHON SOURCE LINES 257-260 .. code-block:: Python f_sol.save("CaO.SiO2_inverse.csdf") # save the solution residuals.save("CaO.SiO2_residue.csdf") # save the residuals .. GENERATED FROM PYTHON SOURCE LINES 261-272 Data Visualization ------------------ At this point, we have solved the inverse problem and obtained an optimum distribution of the nuclear shielding tensor parameters from the 2D MAF dataset. You may use any data visualization and interpretation tool of choice for further analysis. In the following sections, we provide minimal visualization to complete the case study. Visualizing the 3D solution ''''''''''''''''''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 272-286 .. code-block:: Python # Normalize the solution f_sol /= f_sol.max() # Convert the coordinates of the solution, `f_sol`, from Hz to ppm. [item.to("ppm", "nmr_frequency_ratio") for item in f_sol.dimensions] # The 3D plot of the solution plt.figure(figsize=(5, 4.4)) ax = plt.subplot(projection="3d") plot_3d(ax, f_sol, x_lim=[0, 140], y_lim=[0, 140], z_lim=[-50, -120]) plt.tight_layout() plt.show() .. image-sg:: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_004.png :alt: plot 2d 5 CaO.SiO2 :srcset: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 287-292 Convert the 3D tensor distribution in Haeberlen parameters ---------------------------------------------------------- You may re-bin the 3D tensor parameter distribution from a :math:`\rho(\delta_\text{iso}, x, y)` distribution to :math:`\rho(\delta_\text{iso}, \zeta_\sigma, \eta_\sigma)` distribution as follows. .. GENERATED FROM PYTHON SOURCE LINES 292-300 .. code-block:: Python # Create the zeta and eta dimensions,, as shown below. zeta = cp.as_dimension(np.arange(40) * 8 - 150, unit="ppm", label="zeta") eta = cp.as_dimension(np.arange(16) / 15, label="eta") # Use the `to_Haeberlen_grid` function to convert the tensor parameter distribution. fsol_Hae = to_Haeberlen_grid(f_sol, zeta, eta) .. GENERATED FROM PYTHON SOURCE LINES 301-303 The 3D plot ''''''''''' .. GENERATED FROM PYTHON SOURCE LINES 303-309 .. code-block:: Python plt.figure(figsize=(5, 4.4)) ax = plt.subplot(projection="3d") plot_3d(ax, fsol_Hae, x_lim=[0, 1], y_lim=[-150, 150], z_lim=[-50, -120], alpha=0.05) plt.tight_layout() plt.show() .. image-sg:: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_005.png :alt: plot 2d 5 CaO.SiO2 :srcset: /galley_examples/MAF/images/sphx_glr_plot_2d_5_CaO.SiO2_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 310-317 References ---------- .. [#f1] Zhang, P., Grandinetti, P. J., Stebbins, J. F., Anionic Species Determination in CaSiO3 Glass Using Two-Dimensional 29Si NMR, J. Phys. Chem. B, **101**, 4004-4008 (1997). `doi:10.1021/jp9700342. `_ .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 3.857 seconds) .. _sphx_glr_download_galley_examples_MAF_plot_2d_5_CaO.SiO2.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_2d_5_CaO.SiO2.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_2d_5_CaO.SiO2.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_2d_5_CaO.SiO2.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_