.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "fitting/1D_fitting/plot_7_27Al_Extended_Czjzek.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_fitting_1D_fitting_plot_7_27Al_Extended_Czjzek.py: Fitting a Czjzek Model ^^^^^^^^^^^^^^^^^^^^^^ .. GENERATED FROM PYTHON SOURCE LINES 7-19 In this example, we illustrate an application of mrsimulator in fitting a lineshape from a Czjzek distribution of quadrupolar tensors to an experimental :math:`^{27}\text{Al}` MAS spectrum of a phospho-aluminosilicate glass. Setting up the least-squares minimization for a distribution of spin systems is slightly different than that of a crystalline solid. There are 4 steps involved in the processs: - Importing the experimental dataset, - Generating a pre-computed line shape kernel of subspectra, - Creating parameters for the Czjzek distribution model from an initial guess, - Minimizing and visualizing. .. GENERATED FROM PYTHON SOURCE LINES 19-34 .. code-block:: Python import numpy as np import csdmpy as cp import matplotlib.pyplot as plt import mrsimulator.signal_processor as sp from mrsimulator.simulator.config import ConfigSimulator import mrsimulator.utils.spectral_fitting as sf from lmfit import Minimizer from mrsimulator.method.lib import BlochDecayCTSpectrum from mrsimulator.utils import get_spectral_dimensions from mrsimulator.models.czjzek import CzjzekDistribution from mrsimulator.models.utils import LineShapeKernel .. GENERATED FROM PYTHON SOURCE LINES 36-40 Import the experimental dataset ------------------------------- Below we import and visualize the experimental dataset. .. GENERATED FROM PYTHON SOURCE LINES 40-53 .. code-block:: Python host = "http://ssnmr.org/sites/default/files/mrsimulator/" filename = "20K_20Al_10P_50Si_HahnEcho_27Al.csdf" exp_data = cp.load(host + filename).real exp_data.x[0].to("ppm", "nmr_frequency_ratio") exp_data /= exp_data.max() plt.figure(figsize=(4, 3)) ax = plt.subplot(projection="csdm") ax.plot(exp_data) plt.tight_layout() plt.show() .. image-sg:: /fitting/1D_fitting/images/sphx_glr_plot_7_27Al_Extended_Czjzek_001.png :alt: plot 7 27Al Extended Czjzek :srcset: /fitting/1D_fitting/images/sphx_glr_plot_7_27Al_Extended_Czjzek_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 54-81 Generating a line shape kernel ------------------------------ The Czjzek distribution is a statistical model that represents a five-dimensional multivariate normal distribution of second-rank tensor components. These random second-rank tensors are converted into corresponding anisotropy and asymmetry parameters using the Haeberlen convention. The resulting parameters are then binned on a two-dimensional grid, forming the Czjzek probability distribution. The Czjzek spectra is which are the probability weighted sum of lineshapes at each grid point. However, simulating the spectra of spin systems at each minimization step for the least-squares fitting is computationally expensive. A more efficient way is to pre-define a grid system for the tensor parameters, simulate a library of sub-spectra for each grid point, and only update the probability distribution during each minimization step. To simulate the spectra of the given experiment, we create a Method object. Using this method, we generate a kernel of lineshape sub-spectra defined on a polar grid by using the LineShapeKernel class. The argument "method" is the mrsimulator method object used for the simulation, "pos" is a tuple of coordinates defining the two-dimensional grid, "polar" defines a polar parameter coordinate space, and "config" is the Simulator config object used in lineshape simulation. Here, we use the "config" argument to set the number of sidebands as 8. The kernel is generated with the "generate_lineshape()" function, where we pass the "tensor_type='quadrupolar'" argument to specify a quadrupolar parameter grid. A symmetric shielding lineshape kernel can also be generated by specifying .. GENERATED FROM PYTHON SOURCE LINES 81-104 .. code-block:: Python # Create a Method object to simulate the spectrum spectral_dimensions = get_spectral_dimensions(exp_data) method = BlochDecayCTSpectrum( channels=["27Al"], rotor_frequency=14.2e3, spectral_dimensions=spectral_dimensions, experiment=exp_data, ) # Define a polar grid for the lineshape kernel x = np.linspace(0, 2e7, num=36) y = np.linspace(0, 2e7, num=36) pos = (x, y) # Generate the kernel sim_config = ConfigSimulator(number_of_sidebands=8) quad_kernel = LineShapeKernel(method=method, pos=pos, polar=True, config=sim_config) quad_kernel.generate_lineshape(tensor_type="quadrupolar") print("Desired Kernel shape: ", (x.size * y.size, spectral_dimensions[0]["count"])) print("Actual Kernel shape: ", quad_kernel.kernel.shape) .. rst-class:: sphx-glr-script-out .. code-block:: none Desired Kernel shape: (1296, 512) Actual Kernel shape: (1296, 512) .. GENERATED FROM PYTHON SOURCE LINES 105-109 Create a Parameters object -------------------------- Next, we create an instance of the `CzjzekDistribution` class with initial guess values along with a `SignalProcessor` object. .. GENERATED FROM PYTHON SOURCE LINES 109-125 .. code-block:: Python # Create initial guess CzjzekDistribution cz_model = CzjzekDistribution( mean_isotropic_chemical_shift=60.0, sigma=1.4e6, polar=True ) all_models = [cz_model] processor = sp.SignalProcessor( operations=[ sp.IFFT(), sp.apodization.Gaussian(FWHM="600 Hz"), sp.FFT(), sp.Scale(factor=0.3), ] ) .. GENERATED FROM PYTHON SOURCE LINES 126-130 Make the Parameters object and simulate a guess spectrum. Note that the variable `sf_kwargs` holds some additional keyword arguments that many of the spectral fitting function takes in. This dictionary needs to be updated to reflect any changes made in the minimization. .. GENERATED FROM PYTHON SOURCE LINES 130-157 .. code-block:: Python params = sf.make_LMFIT_params(spin_system_models=all_models, processors=[processor]) # Additional keyword arguments passed to best-fit and residual functions. sf_kwargs = dict( kernel=quad_kernel, spin_system_models=all_models, processor=processor, ) # Make a guess and residuals spectrum from the initial guess guess = sf.bestfit_dist(params=params, **sf_kwargs) residuals = sf.residuals_dist(params=params, **sf_kwargs) plt.figure(figsize=(4, 3)) ax = plt.subplot(projection="csdm") ax.plot(exp_data, "k", alpha=0.5, label="Experiment") ax.plot(guess, "r", alpha=0.3, label="Guess") ax.plot(residuals, "b", alpha=0.3, label="Residuals") plt.legend() plt.grid() plt.title("Initial Guess") plt.tight_layout() plt.show() # Print the Parameters object params .. image-sg:: /fitting/1D_fitting/images/sphx_glr_plot_7_27Al_Extended_Czjzek_002.png :alt: Initial Guess :srcset: /fitting/1D_fitting/images/sphx_glr_plot_7_27Al_Extended_Czjzek_002.png :class: sphx-glr-single-img .. raw:: html


.. GENERATED FROM PYTHON SOURCE LINES 158-165 Create and run a minimization ----------------------------- Finally, a `Minimizer` object is created and a minimization run using least-squares. The same arguments defined in the `addtl_sf_kwargs` variable are also passed to the minimizer. Sinice the probabilty distribution is generated from a sparsely sampled from a 5D second rank tensor parameter space, we increase the `diff_step` size from machine precession to avoid approaching local minima from noise. .. GENERATED FROM PYTHON SOURCE LINES 165-182 .. code-block:: Python scipy_minimization_kwargs = dict( diff_step=1e-4, # Increase step size from machine precesion. gtol=1e-10, # Decrease global convergence requirement (default 1e-8) xtol=1e-10, # Decrease variable convergence requirement (default 1e-8) verbose=2, # Print minimization info during each step loss="linear", ) minner = Minimizer( sf.LMFIT_min_function_dist, params, fcn_kws=sf_kwargs, **scipy_minimization_kwargs, ) result = minner.minimize(method="least_squares") result .. rst-class:: sphx-glr-script-out .. code-block:: none Iteration Total nfev Cost Cost reduction Step norm Optimality 0 1 1.9057e+00 4.06e+00 1 2 4.6019e-01 1.45e+00 1.00e+05 5.95e+00 2 3 8.4380e-02 3.76e-01 5.94e+04 6.06e-01 3 4 1.3069e-02 7.13e-02 2.66e+03 4.21e-01 4 5 6.4617e-03 6.61e-03 4.65e+04 2.81e-02 5 6 6.4205e-03 4.12e-05 4.78e+03 1.17e-03 6 7 6.4204e-03 6.98e-08 6.18e+01 1.38e-05 7 8 6.4204e-03 3.96e-10 8.39e+00 6.97e-07 8 10 6.4204e-03 4.41e-10 9.13e-02 3.52e-07 9 12 6.4204e-03 0.00e+00 0.00e+00 3.52e-07 `xtol` termination condition is satisfied. Function evaluations 12, initial cost 1.9057e+00, final cost 6.4204e-03, first-order optimality 3.52e-07. .. raw:: html

Fit Result



.. GENERATED FROM PYTHON SOURCE LINES 183-185 Plot the best-fit spectrum -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 185-199 .. code-block:: Python bestfit = sf.bestfit_dist(params=result.params, **sf_kwargs) residuals = sf.residuals_dist(params=result.params, **sf_kwargs) plt.figure(figsize=(4, 3)) ax = plt.subplot(projection="csdm") ax.plot(exp_data, "k", alpha=0.5, label="Experiment") ax.plot(bestfit, "r", alpha=0.3, label="Fit") ax.plot(residuals, "b", alpha=0.3, label="Residuals") plt.legend() plt.grid() plt.title("Best Fit") plt.tight_layout() plt.show() .. image-sg:: /fitting/1D_fitting/images/sphx_glr_plot_7_27Al_Extended_Czjzek_003.png :alt: Best Fit :srcset: /fitting/1D_fitting/images/sphx_glr_plot_7_27Al_Extended_Czjzek_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 200-203 Plot the best-fit distribution ------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 203-215 .. code-block:: Python for i, model in enumerate(all_models): model.update_lmfit_params(result.params, i) amp = cz_model.pdf(pos=pos, pack_as_csdm=True) plt.figure(figsize=(4, 3)) ax = plt.subplot(projection="csdm") ax.imshow(amp, cmap="gist_ncar_r", interpolation="none", aspect="auto") ax.set_xlabel("x / Hz") ax.set_ylabel("y / Hz") plt.tight_layout() plt.show() .. image-sg:: /fitting/1D_fitting/images/sphx_glr_plot_7_27Al_Extended_Czjzek_004.png :alt: plot 7 27Al Extended Czjzek :srcset: /fitting/1D_fitting/images/sphx_glr_plot_7_27Al_Extended_Czjzek_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 6.383 seconds) .. _sphx_glr_download_fitting_1D_fitting_plot_7_27Al_Extended_Czjzek.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_7_27Al_Extended_Czjzek.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_7_27Al_Extended_Czjzek.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_