{
  "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 Cusipidine\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "After acquiring an NMR spectrum, we often require a least-squares analysis to\ndetermine site populations and nuclear spin interaction parameters. Generally, this\ncomprises of two steps:\n\n- create a fitting model, and\n- determine the model parameters that give the best fit to the spectrum.\n\nHere, we will use the mrsimulator objects to create a fitting model, and use the\n`LMFIT <https://lmfit.github.io/lmfit-py/>`_ library for performing the least-squares\nfitting optimization.\nIn this example, we use a synthetic $^{29}\\text{Si}$ NMR spectrum of cuspidine,\ngenerated from the tensor parameters reported by Hansen `et. al.` [#f1]_, to\ndemonstrate a simple fitting procedure.\n\nWe will begin by importing relevant modules and establishing figure size.\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 mrsimulator import Simulator, SpinSystem\nfrom mrsimulator.methods import BlochDecaySpectrum\nfrom lmfit import Minimizer, Parameters, fit_report\n\nfont = {\"size\": 9}\nmpl.rc(\"font\", **font)\nmpl.rcParams[\"figure.figsize\"] = [4.5, 3.0]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Import the dataset\nUse the `csdmpy <https://csdmpy.readthedocs.io/en/stable/index.html>`_\nmodule to load the synthetic dataset as a CSDM object.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "file_ = \"https://sandbox.zenodo.org/record/687656/files/synthetic_cuspidine_test.csdf\"\nsynthetic_experiment = cp.load(file_)\n\n# convert the dimension coordinates from Hz to ppm\nsynthetic_experiment.dimensions[0].to(\"ppm\", \"nmr_frequency_ratio\")\n\n# Normalize the spectrum\nsynthetic_experiment /= synthetic_experiment.max()\n\n# Plot of the synthetic dataset.\nax = plt.subplot(projection=\"csdm\")\nax.plot(synthetic_experiment, color=\"black\", linewidth=1)\nax.set_xlim(-200, 50)\nax.invert_xaxis()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Create a fitting model\n\nBefore you can fit a simulation to an experiment, in this case, the synthetic dataset,\nyou will first need to create a fitting model. We will use the ``mrsimulator`` objects\nas tools in creating a model for the least-squares fitting.\n\n**Step 1:** Create initial guess sites and spin systems. The initial guess is\noften based on some prior knowledge about the system under investigation. For the\ncurrent example, we know that Cuspidine is a crystalline silica polymorph with one\ncrystallographic Si site. Therefore, our initial guess model is a single\n$^{29}\\text{Si}$ site spin system. For non-linear fitting algorithms, as a\ngeneral recommendation, the initial guess model parameters should be a good starting\npoint for the algorithms to converge.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# the guess model comprising of a single site spin system\nsite = dict(\n    isotope=\"29Si\",\n    isotropic_chemical_shift=-82.0,  # in ppm,\n    shielding_symmetric={\"zeta\": -63, \"eta\": 0.4},  # zeta in ppm\n)\n\nsystem_object = SpinSystem(\n    name=\"Si Site\",\n    description=\"A 29Si site in cuspidine\",\n    sites=[site],  # from the above code\n    abundance=100,\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 2:** Create the method object. The method should be the same as the one used\nin the measurement. In this example, we use the `BlochDecaySpectrum` method. Note,\nwhen creating the method object, the value of the method parameters must match the\nrespective values used in the experiment.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "method = BlochDecaySpectrum(\n    channels=[\"29Si\"],\n    magnetic_flux_density=7.1,  # in T\n    rotor_frequency=780,  # in Hz\n    spectral_dimensions=[\n        {\n            \"count\": 2048,\n            \"spectral_width\": 25000,  # in Hz\n            \"reference_offset\": -5000,  # in Hz\n        }\n    ],\n)"
      ]
    },
    {
      "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]\n\nsim.methods[0].experiment = synthetic_experiment"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 5:** Simulate the spectrum.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "sim.run()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 6:** Create a SignalProcessor class and apply post simulation processing.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "processor = sp.SignalProcessor(\n    operations=[\n        sp.IFFT(),\n        apo.Exponential(FWHM=\"200 Hz\"),\n        sp.FFT(),\n        sp.Scale(factor=1.5),\n    ]\n)\nprocessed_data = processor.apply_operations(data=sim.methods[0].simulation)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "**Step 7:** The plot the spectrum. We also plot the synthetic dataset for comparison.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "ax = plt.subplot(projection=\"csdm\")\nax.plot(processed_data.real, c=\"k\", linewidth=1, label=\"guess spectrum\")\nax.plot(synthetic_experiment.real, c=\"r\", linewidth=1.5, alpha=0.5, label=\"experiment\")\nax.set_xlim(-200, 50)\nax.invert_xaxis()\nplt.legend()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Setup a Least-squares minimization\n\nNow that our model is ready, the next step is to set up a least-squares minimization.\nYou may use any optimization package of choice, here we show an application using\nLMFIT. You may read more on the LMFIT\n`documentation page <https://lmfit.github.io/lmfit-py/index.html>`_.\n\n### Create fitting parameters\n\nNext, you will need a list of parameters that will be used in the fit. The *LMFIT*\nlibrary provides a `Parameters <https://lmfit.github.io/lmfit-py/parameters.html>`_\nclass to create a list of parameters.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "site1 = system_object.sites[0]\nparams = Parameters()\n\nparams.add(name=\"iso\", value=site1.isotropic_chemical_shift)\nparams.add(name=\"eta\", value=site1.shielding_symmetric.eta, min=0, max=1)\nparams.add(name=\"zeta\", value=site1.shielding_symmetric.zeta)\nparams.add(name=\"FWHM\", value=processor.operations[1].FWHM)\nparams.add(name=\"factor\", value=processor.operations[3].factor)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "### Create a minimization function\n\nNote, the above set of parameters does not know about the model. You will need to\nset up a function that will\n\n- update the parameters of the `Simulator` and `SignalProcessor` object based on the\n  LMFIT parameter updates,\n- re-simulate the spectrum based on the updated values, and\n- return the difference between the experiment and simulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def minimization_function(params, sim, processor):\n    values = params.valuesdict()\n\n    # the experiment data as a Numpy array\n    intensity = sim.methods[0].experiment.dependent_variables[0].components[0].real\n\n    # Here, we update simulation parameters iso, eta, and zeta for the site object\n    site = sim.spin_systems[0].sites[0]\n    site.isotropic_chemical_shift = values[\"iso\"]\n    site.shielding_symmetric.eta = values[\"eta\"]\n    site.shielding_symmetric.zeta = values[\"zeta\"]\n\n    # run the simulation\n    sim.run()\n\n    # update the SignalProcessor parameter and apply line broadening.\n    # update the scaling factor parameter at index 3 of operations list.\n    processor.operations[3].factor = values[\"factor\"]\n    # update the exponential apodization FWHM parameter at index 1 of operations list.\n    processor.operations[1].FWHM = values[\"FWHM\"]\n\n    # apply signal processing\n    processed_data = processor.apply_operations(sim.methods[0].simulation)\n\n    # return the difference vector.\n    return intensity - processed_data.dependent_variables[0].components[0].real"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "<div class=\"alert alert-info\"><h4>Note</h4><p>To automate the fitting process, we provide a function to parse the\n      ``Simulator`` and ``SignalProcessor`` objects for parameters and construct an\n      *LMFIT* ``Parameters`` object. Similarly, a minimization function, analogous to\n      the above `minimization_function`, is also included in the *mrsimulator*\n      library. See the next example for usage instructions.</p></div>\n\n### Perform the least-squares minimization\n\nWith the synthetic dataset, simulation, and the initial guess parameters, we are ready\nto perform the fit. To fit, we use the *LMFIT*\n`Minimizer <https://lmfit.github.io/lmfit-py/fitting.html>`_ class.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "minner = Minimizer(minimization_function, params, fcn_args=(sim, processor))\nresult = minner.minimize()\nprint(fit_report(result))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "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 = synthetic_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(\"Frequency / Hz\")\nplt.xlim(-200, 50)\nplt.gca().invert_xaxis()\nplt.grid(which=\"major\", axis=\"both\", linestyle=\"--\")\nplt.legend()\nplt.tight_layout()\nplt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        ".. [#f1] Hansen, M. R., Jakobsen, H. J., Skibsted, J., $^{29}\\text{Si}$\n      Chemical Shift Anisotropies in Calcium Silicates from High-Field\n      $^{29}\\text{Si}$ MAS NMR Spectroscopy, Inorg. Chem. 2003,\n      **42**, *7*, 2368-2377.\n      `DOI: 10.1021/ic020647f <https://doi.org/10.1021/ic020647f>`_\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
}