{
  "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 17O MAS NMR of crystalline Na2SiO3\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this example, we illustrate the use of the mrsimulator objects to\n\n- create a spin system fitting model,\n- use the fitting model to perform a least-squares fit on the experimental, and\n- extract the tensor parameters of the spin system model.\n\nWe will be using the `LMFIT <https://lmfit.github.io/lmfit-py/>`_ methods to\nestablish fitting parameters and fit the spectrum. The following example illustrates\nthe least-squares fitting on a $^{17}\\text{O}$ measurement of\n$\\text{Na}_{2}\\text{SiO}_{3}$ [#f5]_.\n\nWe will begin by importing relevant modules and presetting figure style and layout.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import csdmpy as cp\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport mrsimulator.signal_processing as sp\nimport mrsimulator.signal_processing.apodization as apo\nfrom lmfit import Minimizer, report_fit\nfrom mrsimulator import Simulator, SpinSystem, Site\nfrom mrsimulator.methods import BlochDecayCentralTransitionSpectrum\nfrom mrsimulator.utils import get_spectral_dimensions\nfrom mrsimulator.utils.spectral_fitting import LMFIT_min_function, make_LMFIT_params\n\nfont = {\"size\": 9}\nmpl.rc(\"font\", **font)\nmpl.rcParams[\"figure.figsize\"] = [4.25, 3.0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Import the dataset\n\nImport the experimental data. In this example, we will import the dataset file\nserialized with the CSDM file-format, using the\n`csdmpy <https://csdmpy.readthedocs.io/en/stable/index.html>`_ module.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "filename = \"https://sandbox.zenodo.org/record/687656/files/Na2SiO3_O17.csdf\"\noxygen_experiment = cp.load(filename)\n\n# For spectral fitting, we only focus on the real part of the complex dataset\noxygen_experiment = oxygen_experiment.real\n\n# Convert the dimension coordinates from Hz to ppm.\noxygen_experiment.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n\n# Normalize the spectrum\noxygen_experiment /= oxygen_experiment.max()\n\n# plot of the dataset.\nax = plt.subplot(projection=\"csdm\")\nax.plot(oxygen_experiment, color=\"black\", linewidth=1)\nax.set_xlim(-50, 100)\nax.invert_xaxis()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Create a fitting model\n\nNext, we will create a ``simulator`` object that we use to fit the spectrum. We will\nstart by creating the guess ``SpinSystem`` objects.\n\n**Step 1:** Create initial guess sites and spin systems\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "O17_1 = Site(\n    isotope=\"17O\",\n    isotropic_chemical_shift=60.0,  # in ppm,\n    quadrupolar={\"Cq\": 4.2e6, \"eta\": 0.5},  # Cq in Hz\n)\n\nO17_2 = Site(\n    isotope=\"17O\",\n    isotropic_chemical_shift=40.0,  # in ppm,\n    quadrupolar={\"Cq\": 2.4e6, \"eta\": 0},  # Cq in Hz\n)\n\nsystem_object = [SpinSystem(sites=[s], abundance=50) for s in [O17_1, O17_2]]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 2:** Create the method object. Note, when performing the least-squares fit, you\nmust create an appropriate method object which matches the method used in acquiring\nthe experimental data. The attribute values of this method must match the\nexact conditions under which the experiment was acquired. This including the\nacquisition channels, the magnetic flux density, rotor angle, rotor frequency, and\nthe spectral/spectroscopic dimension. In the following example, we set up a central\ntransition selective Bloch decay spectrum method, where we obtain the\nspectral/spectroscopic information from the metadata of the CSDM dimension. Use the\n:func:`~mrsimulator.utils.get_spectral_dimensions` utility function for quick\nextraction of the spectroscopic information, `i.e.`, count, spectral_width, and\nreference_offset from the CSDM object. The remaining attribute values are set to the\nexperimental conditions.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# get the count, spectral_width, and reference_offset information from the experiment.\nspectral_dims = get_spectral_dimensions(oxygen_experiment)\n\nmethod = BlochDecayCentralTransitionSpectrum(\n    channels=[\"17O\"],\n    magnetic_flux_density=9.4,  # in T\n    rotor_frequency=14000,  # in Hz\n    spectral_dimensions=spectral_dims,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Assign the experimental dataset to the ``experiment`` attribute of the above method.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "method.experiment = oxygen_experiment"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 3:** Create the Simulator object and add the method and spin system objects.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sim = Simulator()\nsim.spin_systems = system_object\nsim.methods = [method]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 4:** Simulate the spectrum.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for iso in sim.spin_systems:\n    # A method object queries every spin system for a list of transition pathways that\n    # are relevant for the given method. Since the method and the number of spin systems\n    # remain the same during the least-squares fit, a one-time query is sufficient. To\n    # avoid querying for the transition pathways at every iteration in a least-squares\n    # fitting, evaluate the transition pathways once and store it as follows\n    iso.transition_pathways = method.get_transition_pathways(iso)\n\n# Now simulate as usual.\nsim.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 5:** Create the SignalProcessor class object and apply the post-simulation\nsignal processing operations.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "processor = sp.SignalProcessor(\n    operations=[\n        sp.IFFT(),\n        apo.Gaussian(FWHM=\"100 Hz\"),\n        sp.FFT(),\n        sp.Scale(factor=20000.0),\n    ]\n)\nprocessed_data = processor.apply_operations(data=sim.methods[0].simulation)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 6:** The plot of initial guess simulation (black) along with the experiment\n(red) is shown below.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = plt.subplot(projection=\"csdm\")\nax.plot(processed_data.real, color=\"black\", linewidth=1, label=\"guess spectrum\")\nax.plot(oxygen_experiment, c=\"r\", linewidth=1.5, alpha=0.5, label=\"experiment\")\nax.set_xlim(-50, 100)\nax.invert_xaxis()\nplt.legend()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Least-squares minimization with LMFIT\n\nOnce you have a fitting model, you need to create the list of parameters to use in the\nleast-squares fitting. For this, you may use the\n`Parameters <https://lmfit.github.io/lmfit-py/parameters.html>`_ class from *LMFIT*,\nas described in the previous example.\nHere, we make use of a utility function,\n:func:`~mrsimulator.utils.spectral_fitting.make_LMFIT_params`, that considerably\nsimplifies the LMFIT parameters generation process.\n\n**Step 7:** Create a list of parameters.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "params = make_LMFIT_params(sim, processor)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The `make_LMFIT_params` parses the instances of the ``Simulator`` and the\n``PostSimulator`` objects for parameters and returns an LMFIT `Parameters` object.\n\n**Customize the Parameters:**\nYou may customize the parameters list, ``params``, as desired. Here, we remove the\nabundance of the two spin systems and constrain it to the initial value of 50% each.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "params.pop(\"sys_0_abundance\")\nparams.pop(\"sys_1_abundance\")\nparams.pretty_print()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 8:** Perform least-squares minimization. For the user's convenience, we also\nprovide a utility function,\n:func:`~mrsimulator.utils.spectral_fitting.LMFIT_min_function`, for evaluating the\ndifference vector between the simulation and experiment, based on\nthe parameters update. You may use this function directly as the argument of the\nLMFIT Minimizer class, as follows,\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": [
        "**Step 9:** The plot of the fit, measurement and the residuals is shown below.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figsize = (4, 3)\nx, y_data = oxygen_experiment.to_list()\nresidual = result.residual\nplt.plot(x, y_data, label=\"Spectrum\")\nplt.plot(x, y_data - residual, \"r\", alpha=0.5, label=\"Fit\")\nplt.plot(x, residual, alpha=0.5, label=\"Residual\")\n\nplt.xlabel(\"$^{17}$O frequency / ppm\")\nplt.xlim(-50, 100)\nplt.gca().invert_xaxis()\nplt.grid(which=\"major\", axis=\"both\", linestyle=\"--\")\nplt.legend()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. [#f5] T. M. Clark, P. Florian, J. F. Stebbins, and P. J. Grandinetti,\n      An $^{17}\\text{O}$ NMR Investigation of Crystalline Sodium Metasilicate:\n      Implications for the Determination of Local Structure in Alkali Silicates,\n      J. Phys. Chem. B. 2001, **105**, 12257-12265.\n      `DOI: 10.1021/jp011289p  <https://doi.org/10.1021/jp011289p>`_\n\n"
      ]
    }
  ],
  "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
}