{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# This cell is added by sphinx-gallery\n\n%matplotlib inline\n\nimport mrsimulator\nprint(f'You are using mrsimulator v{mrsimulator.__version__}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Fitting PASS/MAT cross-sections\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This example illustrates the use of mrsimulator and LMFIT modules in fitting the\nsideband intensity profile across the isotropic chemical shift cross-section from a\nPASS/MAT dataset.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "import numpy as np\nimport csdmpy as cp\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport mrsimulator.signal_processing as sp\nfrom mrsimulator import Simulator, SpinSystem, Site\nfrom mrsimulator.methods import BlochDecaySpectrum\nfrom mrsimulator.utils import get_spectral_dimensions\nfrom mrsimulator.utils.spectral_fitting import LMFIT_min_function, make_LMFIT_params\nfrom lmfit import Minimizer, report_fit\n\n\n# global plot configuration\nmpl.rcParams[\"figure.figsize\"] = [4.5, 3.0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Import the dataset\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "filename = \"https://sandbox.zenodo.org/record/687656/files/1H13C_CPPASS_LHistidine.csdf\"\npass_data = cp.load(filename)\n\n# For the spectral fitting, we only focus on the real part of the complex dataset.\n# The script assumes that the dimension at index 0 is the isotropic dimension.\n# Transpose the dataset as required.\npass_data = pass_data.real.T\n\n# Convert the coordinates along each dimension from Hz to ppm.\n_ = [item.to(\"ppm\", \"nmr_frequency_ratio\") for item in pass_data.dimensions]\n\n# Normalize the spectrum.\npass_data /= pass_data.max()\n\n# The plot of the dataset.\nlevels = (np.arange(10) + 0.3) / 15 # contours are drawn at these levels.\nax = plt.subplot(projection=\"csdm\")\ncb = ax.contour(pass_data, colors=\"k\", levels=levels, alpha=0.5, linewidths=0.5)\nplt.colorbar(cb)\nax.set_xlim(200, 10)\nax.invert_yaxis()\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Extract a 1D sideband intensity cross-section from the 2D dataset using the array\nindexing.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "data1D = pass_data[1100] # sideband dataset\n\n# The plot of the cross-section.\nax = plt.subplot(projection=\"csdm\")\nax.plot(data1D, color=\"k\")\nax.invert_xaxis()\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The isotropic chemical shift coordinate of the cross-section is\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "isotropic_shift = pass_data.x[0].coords[1100]\nprint(isotropic_shift)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a fitting model\n\nThe fitting model includes the Simulator and SignalProcessor objects. First,\ncreate the Simulator object.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Create the guess site and spin system for the 1D cross-section.\nzeta = -70 # in ppm\neta = 0.8\n\nsite = Site(\n isotope=\"13C\",\n isotropic_chemical_shift=0,\n shielding_symmetric={\"zeta\": zeta, \"eta\": eta},\n)\nspin_systems = [SpinSystem(sites=[site])]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For the sideband only cross-section, use the BlochDecaySpectrum method.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Get the dimension information from the experiment. Note, the following function\n# returns an array of two spectral dimensions corresponding to the 2D PASS dimensions.\n# Use the spectral dimension that is along the anisotropic dimensions for the\n# BlochDecaySpectrum method.\nspectral_dims = get_spectral_dimensions(pass_data)\nmethod = BlochDecaySpectrum(\n channels=[\"13C\"],\n magnetic_flux_density=9.4, # in T\n rotor_frequency=1500, # in Hz\n spectral_dimensions=[spectral_dims[0]],\n experiment=data1D, # also add the measurement to the method.\n)\n\n# Optimize the script by pre-setting the transition pathways for each spin system from\n# the method.\nfor sys in spin_systems:\n sys.transition_pathways = method.get_transition_pathways(sys)\n\n# Create the Simulator object and add the method and spin system objects.\nsim = Simulator()\nsim.spin_systems = spin_systems # add the spin systems\nsim.methods = [method] # add the method\nsim.run()\n\n# Add and apply Post simulation processing.\nprocessor = sp.SignalProcessor(operations=[sp.Scale(factor=1)])\nprocessed_data = processor.apply_operations(data=sim.methods[0].simulation).real\n\n# The plot of the simulation from the guess model and experiment cross-section.\nax = plt.subplot(projection=\"csdm\")\nax.plot(processed_data, color=\"r\", label=\"guess\")\nax.plot(data1D, color=\"k\", label=\"experiment\")\nax.invert_xaxis()\nplt.tight_layout()\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Least-squares minimization with LMFIT\nFirst, create the fitting parameters.\nUse the :func:`~mrsimulator.utils.spectral_fitting.make_LMFIT_params` for a quick\nsetup.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "params = make_LMFIT_params(sim, processor)\n\n# Fix the value of the isotropic chemical shift to zero for pure anisotropic sideband\n# amplitude simulation.\nparams[\"sys_0_site_0_isotropic_chemical_shift\"].vary = False\nparams.pretty_print()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Run the minimization using LMFIT\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "minner = Minimizer(LMFIT_min_function, params, fcn_args=(sim, processor))\nresult = minner.minimize()\nreport_fit(result)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Simulate the spectrum corresponding to the optimum parameters\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "sim.run()\nprocessed_data = processor.apply_operations(data=sim.methods[0].simulation).real" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the spectrum\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ax = plt.subplot(projection=\"csdm\")\nax.plot(processed_data, color=\"r\", label=\"fit\")\nax.plot(data1D, color=\"k\", label=\"experiment\")\nax.invert_xaxis()\nplt.tight_layout()\nplt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 0 }