diff --git a/.github/workflows/continuous-integration-pip.yml b/.github/workflows/continuous-integration-pip.yml index dbf083260..c1bc9921f 100644 --- a/.github/workflows/continuous-integration-pip.yml +++ b/.github/workflows/continuous-integration-pip.yml @@ -64,6 +64,7 @@ jobs: python -m pip install --upgrade pip pip install pytest cython setuptools pytest-cov sympy pip install -r requirements.txt + pip install numpy -U - name: Build and install package from source run: python setup.py develop @@ -143,6 +144,7 @@ jobs: run: | conda --version which python + pip install numpy -U python setup.py develop - name: Test with pytest shell: pwsh diff --git a/.gitignore b/.gitignore index cc3f58d94..9159f25ed 100644 --- a/.gitignore +++ b/.gitignore @@ -117,6 +117,7 @@ class-diagram* *.csdf *.csdfe *.mrsys +*.json.gz test_os.py note.txt simulate_lineshape.py diff --git a/docs/api_py/py-fitting.rst b/docs/api_py/py-fitting.rst deleted file mode 100644 index 93f85c8e3..000000000 --- a/docs/api_py/py-fitting.rst +++ /dev/null @@ -1,14 +0,0 @@ -.. _fitting_api: - -Fitting Utility API -=================== - -LMFIT supplement functions --------------------------- - -.. currentmodule:: mrsimulator.utils.spectral_fitting - -.. autofunction:: make_LMFIT_params -.. autofunction:: LMFIT_min_function -.. autofunction:: bestfit -.. autofunction:: residuals diff --git a/docs/api_py/py-utilities.rst b/docs/api_py/py-utilities.rst new file mode 100644 index 000000000..d4b34bf13 --- /dev/null +++ b/docs/api_py/py-utilities.rst @@ -0,0 +1,32 @@ +.. _fitting_api: + +Utilities API +============= + +LMFIT Utilities +--------------- + +.. currentmodule:: mrsimulator.utils.spectral_fitting + +.. autofunction:: make_LMFIT_params +.. autofunction:: LMFIT_min_function +.. autofunction:: bestfit +.. autofunction:: residuals + + +mpcontribs Utilities +-------------------- + +.. currentmodule:: mrsimulator.contribs + +.. autofunction:: mpcontribs_export + + +.. currentmodule:: mrsimulator.contribs.base + +.. autoclass:: ChemicalShiftSchema +.. autoclass:: QuadrupolarSchema +.. autoclass:: SiteSchema +.. autoclass:: MethodSchema +.. autoclass:: SimulatorSchema +.. autoclass:: ContribSchema diff --git a/docs/index.rst b/docs/index.rst index 8ba30da26..b431a97ba 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -232,7 +232,7 @@ API and references api_py/py-simulator api_py/py-signal-processing api_py/py-model - api_py/py-fitting + api_py/py-utilities api_c/c_api diff --git a/docs/installation/requirements.rst b/docs/installation/requirements.rst index 3334c649f..67d165db3 100644 --- a/docs/installation/requirements.rst +++ b/docs/installation/requirements.rst @@ -7,7 +7,7 @@ Package dependencies **Required packages** -- `NumPy>=1.17 `_ +- `numpy>=1.17 `_ - openblas - cython>=0.29.14 - typing-extensions>=3.7 diff --git a/docs/notebooks/fitting/1D_fitting/plot_2_xonotlite.ipynb b/docs/notebooks/fitting/1D_fitting/plot_2_xonotlite.ipynb new file mode 100644 index 000000000..d26f38072 --- /dev/null +++ b/docs/notebooks/fitting/1D_fitting/plot_2_xonotlite.ipynb @@ -0,0 +1,230 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# This cell is added by sphinx-gallery\n!pip install mrsimulator --quiet\n\n\n%matplotlib inline\n\nimport mrsimulator\nprint(f'You are using mrsimulator v{mrsimulator.__version__}')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "\n# 29Si 1D MAS spinning sideband (Xonotlite)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The following is an example for submitting the NMR tensor parameters to mpcontribs.\nWe use the $^{29}\\text{Si}$ 1D MAS NMR spectrum of Xonotlite crystal by Hansen\net al. [#f1]_ for demonstration.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "import csdmpy as cp\nimport matplotlib.pyplot as plt\nfrom lmfit import Minimizer\n\nfrom mrsimulator import Simulator, SpinSystem, Site\nfrom mrsimulator.methods import BlochDecaySpectrum\nfrom mrsimulator import signal_processing as sp\nfrom mrsimulator.utils import spectral_fitting as sf\nfrom mrsimulator.utils import get_spectral_dimensions" + ] + }, + { + "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/744498/files/xonotlite.csdf\"\nexperiment = cp.load(filename).real\n\n# standard deviation of noise from the dataset\nsigma = 2.819601\n\n# Convert the coordinates along each dimension from Hz to ppm.\n_ = [item.to(\"ppm\", \"nmr_frequency_ratio\") for item in experiment.dimensions]\n\n# Plot of the synthetic dataset.\nplt.figure(figsize=(4.25, 3.0))\nax = plt.subplot(projection=\"csdm\")\nax.plot(experiment, \"k\", alpha=0.5)\nax.invert_xaxis()\nplt.tight_layout()\nplt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Create a fitting model\n**Guess model**\n\nCreate a guess list of spin systems. There are three crystallographic\n$^{29}\\text{Si}$ sites in Xonotlite.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "s1 = Site(\n isotope=\"29Si\",\n isotropic_chemical_shift=-97.17, # in ppm,\n shielding_symmetric={\"zeta\": 35.0, \"eta\": 0.0}, # zeta in ppm\n)\ns2 = Site(\n isotope=\"29Si\",\n isotropic_chemical_shift=-86.3, # in ppm,\n shielding_symmetric={\"zeta\": 50.0, \"eta\": 0.5}, # zeta in ppm\n)\ns3 = Site(\n isotope=\"29Si\",\n isotropic_chemical_shift=-87.2, # in ppm,\n shielding_symmetric={\"zeta\": 44.0, \"eta\": 0.5}, # zeta in ppm\n)\nspin_systems = [\n SpinSystem(name=\"Q3\", sites=[s1], abundance=25),\n SpinSystem(name=\"Q2 (1)\", sites=[s2], abundance=75 / 2),\n SpinSystem(name=\"Q2 (2)\", sites=[s3], abundance=75 / 2),\n]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Method**\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Get the spectral dimension paramters from the experiment.\nspectral_dims = get_spectral_dimensions(experiment)\n\nmethod = BlochDecaySpectrum(\n channels=[\"29Si\"],\n magnetic_flux_density=14.1, # in T\n rotor_frequency=1800.0, # in Hz\n spectral_dimensions=spectral_dims,\n experiment=experiment, # 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 das method.\nfor sys in spin_systems:\n sys.transition_pathways = method.get_transition_pathways(sys)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Guess Spectrum**\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# Simulation\n# ----------\nsim = Simulator(spin_systems=spin_systems, methods=[method])\nsim.run()\n\n# Post Simulation Processing\n# --------------------------\nprocessor = sp.SignalProcessor(\n operations=[\n sp.IFFT(), # inverse FFT to convert frequency based spectrum to time domain.\n sp.apodization.Exponential(FWHM=\"50 Hz\"), # apodization of time domain signal.\n sp.FFT(), # forward FFT to convert time domain signal to frequency spectrum.\n sp.Scale(factor=500), # scale the frequency spectrum.\n ]\n)\nprocessed_data = processor.apply_operations(data=sim.methods[0].simulation).real\n\n# Plot of the guess Spectrum\n# --------------------------\nplt.figure(figsize=(4.25, 3.0))\nax = plt.subplot(projection=\"csdm\")\nax.plot(experiment, \"k\", linewidth=1, label=\"Experiment\")\nax.plot(processed_data, \"r\", alpha=0.75, linewidth=1, label=\"guess spectrum\")\nax.invert_xaxis()\nplt.grid()\nplt.legend()\nplt.tight_layout()\nplt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Least-squares minimization with LMFIT\nUse the :func:`~mrsimulator.utils.spectral_fitting.make_LMFIT_params` for a quick\nsetup of the fitting parameters.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "params = sf.make_LMFIT_params(sim, processor, include={\"rotor_frequency\"})\n\nparams.pop(\"sys_0_abundance\")\nparams.pop(\"sys_1_abundance\")\nparams.pop(\"sys_2_abundance\")\nparams[\"sys_0_site_0_shielding_symmetric_eta\"].vary = False\nprint(params.pretty_print(columns=[\"value\", \"min\", \"max\", \"vary\", \"expr\"]))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "**Solve the minimizer using LMFIT**\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma))\nresult = minner.minimize()\nresult" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## The best fit solution\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "best_fit = sf.bestfit(sim, processor)[0]\nresiduals = sf.residuals(sim, processor)[0]\n\n# Plot the spectrum\nplt.figure(figsize=(4.25, 3.0))\nax = plt.subplot(projection=\"csdm\")\nax.plot(experiment, \"k\", linewidth=1, label=\"Experiment\")\nax.plot(best_fit, \"r\", alpha=0.75, linewidth=1, label=\"Best Fit\")\nax.plot(residuals, alpha=0.75, linewidth=1, label=\"Residuals\")\nax.invert_xaxis()\nplt.xlabel(\"$^{29}$Si frequency / ppm\")\nplt.grid()\nplt.legend()\nplt.tight_layout()\nplt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Submitting data to MPContribs\n\nTo contribute to MPContribs, we need to export the mrsimulator objects to a list of\nmp-compatible data dictionaries. At present, MPContribs only support data contribution\non a per NMR site basis and, therefore, we only generate mp contributions for\nuncoupled spin systems. Use the ``mrsimulator.contribs`` module to create data\ndictionaries as follows.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "from mrsimulator.contribs import mpcontribs_export\nfrom pprint import pprint\n\nmp_project = \"lsdi_nmr_exp_test\" # this should be your mpcontribs project name\ncards = mpcontribs_export(\n sim,\n [processor],\n project=mp_project,\n identifier=\"Ca6Si6O17(OH)2\",\n exp_dict={\n \"90degreePulseLength\": \"6 \u00b5s\",\n \"relaxationDelay\": \"8 s\",\n \"numberOfScans\": 7224,\n \"referenceCompound\": \"TMS\",\n },\n)\nprint(\"Number of contributions\", len(cards))\npprint(cards[0][\"data\"])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here, ``cards`` hold a list of mp-data dictionaries. In this example, it corresponds\nto three---the number of uncoupled spin systems.\nTo submit contributions, use the mpcontribs client as shown below.\n\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# from mpcontribs.client import Client\n#\n# client = Client() # uses MPCONTRIBS_API_KEY envvar.\n# client.submit_contributions(cards)" + ] + }, + { + "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 `_\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.10" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} \ No newline at end of file diff --git a/environment-dev.yml b/environment-dev.yml index 087f8ea71..6307a2007 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -5,7 +5,7 @@ channels: dependencies: - fftw>=3.3.0 - openblas - - numpy>=1.17 + - numpy>=1.20 - nomkl - setuptools>=27.3 - cython>=0.29 @@ -36,3 +36,4 @@ dependencies: - recommonmark - sphinx-version-warning - plotly + - mpcontribs-client diff --git a/environment.yml b/environment.yml index c86c6ce25..ba0963c75 100644 --- a/environment.yml +++ b/environment.yml @@ -22,3 +22,4 @@ dependencies: - lmfit>=1.0.2 - psutil>=5.4.8 - plotly + - mpcontribs-client diff --git a/fitting_source/1D_fitting/plot_2_xonotlite.py b/fitting_source/1D_fitting/plot_2_xonotlite.py new file mode 100644 index 000000000..b1eae91b3 --- /dev/null +++ b/fitting_source/1D_fitting/plot_2_xonotlite.py @@ -0,0 +1,206 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- +""" +29Si 1D MAS spinning sideband (Xonotlite) +========================================= +""" +# %% +# The following is an example for submitting the NMR tensor parameters to mpcontribs. +# We use the :math:`^{29}\text{Si}` 1D MAS NMR spectrum of Xonotlite crystal by Hansen +# et al. [#f1]_ for demonstration. +import csdmpy as cp +import matplotlib.pyplot as plt +from lmfit import Minimizer + +from mrsimulator import Simulator, SpinSystem, Site +from mrsimulator.methods import BlochDecaySpectrum +from mrsimulator import signal_processing as sp +from mrsimulator.utils import spectral_fitting as sf +from mrsimulator.utils import get_spectral_dimensions + +# sphinx_gallery_thumbnail_number = 3 + +# %% +# Import the dataset +# ------------------ +filename = "https://sandbox.zenodo.org/record/744498/files/xonotlite.csdf" +experiment = cp.load(filename).real + +# standard deviation of noise from the dataset +sigma = 2.819601 + +# Convert the coordinates along each dimension from Hz to ppm. +_ = [item.to("ppm", "nmr_frequency_ratio") for item in experiment.dimensions] + +# Plot of the synthetic dataset. +plt.figure(figsize=(4.25, 3.0)) +ax = plt.subplot(projection="csdm") +ax.plot(experiment, "k", alpha=0.5) +ax.invert_xaxis() +plt.tight_layout() +plt.show() + +# %% +# Create a fitting model +# ---------------------- +# **Guess model** +# +# Create a guess list of spin systems. There are three crystallographic +# :math:`^{29}\text{Si}` sites in Xonotlite. +s1 = Site( + isotope="29Si", + isotropic_chemical_shift=-97.17, # in ppm, + shielding_symmetric={"zeta": 35.0, "eta": 0.0}, # zeta in ppm +) +s2 = Site( + isotope="29Si", + isotropic_chemical_shift=-86.3, # in ppm, + shielding_symmetric={"zeta": 50.0, "eta": 0.5}, # zeta in ppm +) +s3 = Site( + isotope="29Si", + isotropic_chemical_shift=-87.2, # in ppm, + shielding_symmetric={"zeta": 44.0, "eta": 0.5}, # zeta in ppm +) +spin_systems = [ + SpinSystem(name="Q3", sites=[s1], abundance=25), + SpinSystem(name="Q2 (1)", sites=[s2], abundance=75 / 2), + SpinSystem(name="Q2 (2)", sites=[s3], abundance=75 / 2), +] + +# %% +# **Method** + +# Get the spectral dimension paramters from the experiment. +spectral_dims = get_spectral_dimensions(experiment) + +method = BlochDecaySpectrum( + channels=["29Si"], + magnetic_flux_density=14.1, # in T + rotor_frequency=1800.0, # in Hz + spectral_dimensions=spectral_dims, + experiment=experiment, # add the measurement to the method. +) + +# Optimize the script by pre-setting the transition pathways for each spin system from +# the das method. +for sys in spin_systems: + sys.transition_pathways = method.get_transition_pathways(sys) + +# %% +# **Guess Spectrum** + +# Simulation +# ---------- +sim = Simulator(spin_systems=spin_systems, methods=[method]) +sim.run() + +# Post Simulation Processing +# -------------------------- +processor = sp.SignalProcessor( + operations=[ + sp.IFFT(), # inverse FFT to convert frequency based spectrum to time domain. + sp.apodization.Exponential(FWHM="50 Hz"), # apodization of time domain signal. + sp.FFT(), # forward FFT to convert time domain signal to frequency spectrum. + sp.Scale(factor=500), # scale the frequency spectrum. + ] +) +processed_data = processor.apply_operations(data=sim.methods[0].simulation).real + +# Plot of the guess Spectrum +# -------------------------- +plt.figure(figsize=(4.25, 3.0)) +ax = plt.subplot(projection="csdm") +ax.plot(experiment, "k", linewidth=1, label="Experiment") +ax.plot(processed_data, "r", alpha=0.75, linewidth=1, label="guess spectrum") +ax.invert_xaxis() +plt.grid() +plt.legend() +plt.tight_layout() +plt.show() + +# %% +# Least-squares minimization with LMFIT +# ------------------------------------- +# Use the :func:`~mrsimulator.utils.spectral_fitting.make_LMFIT_params` for a quick +# setup of the fitting parameters. +params = sf.make_LMFIT_params(sim, processor, include={"rotor_frequency"}) + +params.pop("sys_0_abundance") +params.pop("sys_1_abundance") +params.pop("sys_2_abundance") +params["sys_0_site_0_shielding_symmetric_eta"].vary = False +print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"])) + +# %% +# **Solve the minimizer using LMFIT** +minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma)) +result = minner.minimize() +result + +# %% +# The best fit solution +# --------------------- +best_fit = sf.bestfit(sim, processor)[0] +residuals = sf.residuals(sim, processor)[0] + +# Plot the spectrum +plt.figure(figsize=(4.25, 3.0)) +ax = plt.subplot(projection="csdm") +ax.plot(experiment, "k", linewidth=1, label="Experiment") +ax.plot(best_fit, "r", alpha=0.75, linewidth=1, label="Best Fit") +ax.plot(residuals, alpha=0.75, linewidth=1, label="Residuals") +ax.invert_xaxis() +plt.xlabel("$^{29}$Si frequency / ppm") +plt.grid() +plt.legend() +plt.tight_layout() +plt.show() + + +# %% +# Submitting data to MPContribs +# ----------------------------- +# +# To contribute to MPContribs, we need to export the mrsimulator objects to a list of +# mp-compatible data dictionaries. At present, MPContribs only support data contribution +# on a per NMR site basis and, therefore, we only generate mp contributions for +# uncoupled spin systems. Use the ``mrsimulator.contribs`` module to create data +# dictionaries as follows. + +# %% +# from mrsimulator.contribs import mpcontribs_export +# from pprint import pprint + +# mp_project = "lsdi_nmr_exp_test" # this should be your mpcontribs project name +# cards = mpcontribs_export( +# sim, +# [processor], +# project=mp_project, +# identifier="Ca6Si6O17(OH)2", +# exp_dict={ +# "90degreePulseLength": "6 µs", +# "relaxationDelay": "8 s", +# "numberOfScans": 7224, +# "referenceCompound": "TMS", +# }, +# ) +# print("Number of contributions", len(cards)) +# pprint(cards[0]["data"]) + +# %% +# Here, ``cards`` hold a list of mp-data dictionaries. In this example, it corresponds +# to three---the number of uncoupled spin systems. +# To submit contributions, use the mpcontribs client as shown below. + +# from mpcontribs.client import Client +# +# client = Client() # uses MPCONTRIBS_API_KEY envvar. +# client.submit_contributions(cards) + +# %% +# .. [#f1] Hansen, M. R., Jakobsen, H. J., Skibsted, J., :math:`^{29}\text{Si}` +# Chemical Shift Anisotropies in Calcium Silicates from High-Field +# :math:`^{29}\text{Si}` MAS NMR Spectroscopy, Inorg. Chem. 2003, +# **42**, *7*, 2368-2377. +# `DOI: 10.1021/ic020647f `_ diff --git a/requirements-dev.txt b/requirements-dev.txt index 7be1343f0..d8d8cc523 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -9,6 +9,7 @@ pandas>=1.1.3 numexpr>=2.7.1 psutil>=5.4.8 joblib>=1.0.0 +mpcontribs-client>=3.10.1 # Formatting requirements black diff --git a/requirements.txt b/requirements.txt index ffcf220ef..7f8ab6c5e 100644 --- a/requirements.txt +++ b/requirements.txt @@ -10,3 +10,4 @@ pandas>=1.1.3 numexpr>=2.7.1 psutil>=5.4.8 joblib>=1.0.0 +mpcontribs-client \ No newline at end of file diff --git a/setup.py b/setup.py index fe6dcb388..d54490bd0 100644 --- a/setup.py +++ b/setup.py @@ -403,7 +403,7 @@ def fftw_info(self): if USE_CYTHON: ext_modules = cythonize(ext_modules, language_level=3) -extras = {"all": ["matplotlib>=3.3.3"]} +extras = {"all": ["matplotlib>=3.3.3", "mpcontribs-client>=3.10.1"]} description = "A python toolbox for simulating fast real-time solid-state NMR spectra." setup( diff --git a/src/mrsimulator/__init__.py b/src/mrsimulator/__init__.py index 538bb3978..8b33dbc09 100644 --- a/src/mrsimulator/__init__.py +++ b/src/mrsimulator/__init__.py @@ -44,6 +44,7 @@ from mrsimulator.signal_processing import SignalProcessor import json from lmfit import Parameters +from monty.io import zopen def save( @@ -67,15 +68,19 @@ def save( """ py_dict = dict(simulator, signal_processors, params, with_units) - with open(filename, "w", encoding="utf8") as outfile: - json.dump( + check = filename.endswith("gz") + operator = zopen if check else open + with operator(filename, "w") as outfile: + target = json.dumps( py_dict, - outfile, ensure_ascii=False, sort_keys=False, allow_nan=False, separators=(",", ":"), ) + target += "\n" + target = target.encode("utf-8") if check else target + outfile.write(target) def dict( diff --git a/src/mrsimulator/contribs/__init__.py b/src/mrsimulator/contribs/__init__.py new file mode 100644 index 000000000..0e680c2c2 --- /dev/null +++ b/src/mrsimulator/contribs/__init__.py @@ -0,0 +1,163 @@ +# -*- coding: utf-8 -*- +import json + +import numpy as np +from monty.io import zopen +from mpcontribs.client import Attachment +from mrsimulator import dict +from mrsimulator import Simulator +from mrsimulator.utils import flatten_dict + +from .base import ContribSchema + +__author__ = "Deepansh Srivastava" +__email__ = "srivastava.89@osu.edu" + + +SITE_MAPPABLE = { + "isotropic_chemical_shift": "isotropic", + "shielding_symmetric.zeta": "zeta", + "shielding_symmetric.eta": "eta", + "shielding_symmetric.alpha": "alpha", + "shielding_symmetric.beta": "beta", + "shielding_symmetric.gamma": "gamma", + "quadrupolar.Cq": "Cq", + "quadrupolar.eta": "eta", + "quadrupolar.alpha": "alpha", + "quadrupolar.beta": "beta", + "quadrupolar.gamma": "gamma", +} + +NEW_MAPPABLE = { + "shielding_symmetric": "ChemicalShift", + "quadrupolar": "Quadrupolar", + "isotropic_chemical_shift": "ChemicalShift", +} + + +def parse_sites(site): + """Parse the site object and generate site contribution for mpcontribs. + + Args: + Site site: Site object. + """ + if site.shielding_symmetric is not None: + site.shielding_symmetric.zeta *= -1 # convert to shift anisotropy + site_ = flatten_dict(site.json()) + + dict_ = {"ChemicalShift": {}, "Quadrupolar": {}} + for k, v in site_.items(): + root = k.split(".")[0] + new_key = SITE_MAPPABLE[k] if k in SITE_MAPPABLE else k + if root in NEW_MAPPABLE: + dict_[NEW_MAPPABLE[root]][new_key] = site_[k] + else: + dict_[new_key] = v + + _ = [dict_.pop(i) for i in ["ChemicalShift", "Quadrupolar"] if dict_[i] == {}] + return dict_ + + +def parse_method(method): + """Parse the method object and generate method contribution for mpcontribs. + + Args: + Method method: Method object. + """ + gamma = method.channels[0].gyromagnetic_ratio + B0 = method.spectral_dimensions[0].events[0].magnetic_flux_density + larmor_frequency = abs(gamma * B0) # MHz + rotor_frequency = method.spectral_dimensions[0].events[0].rotor_frequency # Hz + spectral_width = method.spectral_dimensions[0].spectral_width # Hz + rotor_angle = method.spectral_dimensions[0].events[0].rotor_angle * 180 / np.pi + + return { + "larmorFrequency": f"{larmor_frequency} MHz", + "spinningFrequency": f"{rotor_frequency} Hz", + "spectralWidth": f"{spectral_width} Hz", + "rotorAngle": f"{rotor_angle:.4f} degree", + } + + +def mpcontribs_export( + simulator: Simulator, + signal_processors: list, + project: str, + identifier: str, + exp_dict: dict = {}, + **kwargs, +): + """Generate mpcontribs contribution entries for every site in the Simulator object. + + Arguments + --------- + Simulator simulator + Object from where the site contributions are extracted. + list signal_processors + A list of SignalProcessor objects. + str project + mpcontribs project name (required). + str identifier + mpcontribs identifier (required). + exp_dict + Additional metadata to use in contribution. + **kwargs + Optional keyword arguments from mpcontribs ContributionsSchema. + + Example + ------- + >>> contribution_data = mpcontribs_export( + ... simulator, + ... processors, + ... project='my_project', + ... identifier='mp-test' + ... ) # doctest:+SKIP + """ + + attach = prepare_attachments(simulator, signal_processors, identifier) + + contribs = [ + ContribSchema( + project=project, + identifier=identifier, + data={ + "site": {**parse_sites(site)}, + "method": {**parse_method(simulator.methods[0])}, + }, + ).json() + for sys in simulator.spin_systems + for site in sys.sites + ] + _ = [item.update({"attachments": attach}) for item in contribs] + + for item in contribs: + item["data"]["method"] = {**item["data"]["method"], **exp_dict} + + return contribs + + +def save_obj(filename, data): + with zopen(filename, "w") as f: + json_str = json.dumps(data) + "\n" # 2. string (i.e. JSON) + json_bytes = json_str.encode("utf-8") + f.write(json_bytes) + + +def load_obj(filename): + with zopen(filename) as f: + json_bytes = f.read() + json_str = json_bytes.decode("utf-8") # 2. string (i.e. JSON) + data = json.loads(json_str) + return data + + +def prepare_attachments(simulator, signal_processors, identifier): + """Prepares the dict data as attachments for contribution.""" + filename = f"{identifier}.mrsim" + mrsim = Attachment.from_data(filename, dict(simulator, signal_processors)) + data = [ + Attachment.from_data(f"{identifier}-exp{i}.csdf", mth.experiment.dict()) + for i, mth in enumerate(simulator.methods) + if mth.experiment is not None + ] + return [mrsim, *data] diff --git a/src/mrsimulator/contribs/base.py b/src/mrsimulator/contribs/base.py new file mode 100644 index 000000000..b284e94ba --- /dev/null +++ b/src/mrsimulator/contribs/base.py @@ -0,0 +1,126 @@ +# -*- coding: utf-8 -*- +# import csdmpy as cp +from pathlib import Path +from typing import List +from typing import Union + +from mpcontribs.client import Attachment +from pydantic import BaseModel + +__author__ = "Deepansh Srivastava" +__email__ = "srivastava.89@osu.edu" + + +class Base(BaseModel): + def json(self): + return Base.fullsimplify(super().dict()) + + @staticmethod + def simplify(val): + """Remove value if it is None.""" + return {k: v for k, v in val.items() if v is not None} + + @staticmethod + def fullsimplify(val): + """Iteratively remove None values from a nested dict.""" + initial = { + k: Base.simplify(Base.fullsimplify(v)) if isinstance(v, dict) else v + for k, v in val.items() + } + return Base.simplify(initial) + + +class ChemicalShiftSchema(Base): + """Schema for Chemical shift parameters. + + Args: + str isotropic: Isotropic chemical shift. + str zeta: Chemical shift anisotropy. + float eta: Shift asymmetry. + str alpha: Euler angle alpha. + str beta: Euler angle beta. + str gamma: Euler angle gamma. + """ + + isotropic: str + zeta: str = None + eta: float = None + alpha: str = None + beta: str = None + gamma: str = None + + +class QuadrupolarSchema(Base): + """Schema for Quadrupolar parameters. + + Args: + str Cq: Quadrupolar coupling constant. + float eta: Quadrupolar asymmetry. + str alpha: Euler angle alpha. + str beta: Euler angle beta. + str gamma: Euler angle gamma. + """ + + Cq: str + eta: float + alpha: str = None + beta: str = None + gamma: str = None + + +class SiteSchema(Base): + """Schema for Site parameters. + + Args: + str isotope: The isotope given as atomic number followed by atomic symbol. + ChemicalShiftSchema ChemicalShift: The chemical shift scheme. + QuadrupolarSchema Quadrupolar: The quadrupolar scheme. + """ + + isotope: str + ChemicalShift: ChemicalShiftSchema = None + Quadrupolar: QuadrupolarSchema = None + + +class MethodSchema(Base): + """Schema for Method parameters. + + Args: + str larmorFrequency: The larmor frequency of observed isotope. + str spinningFrequency: The rotor spinning frequency. + str spectralWidth: The spectral width of the spectrum. + str rotorAngle: Angle of sample rotation axis wrt the z-axis. + """ + + larmorFrequency: str + spinningFrequency: str + spectralWidth: str + rotorAngle: str + + +class SimulatorSchema(Base): + """Schema for Simulator parameters. + + Args: + SiteSchema site: The site object. + MethodSchema method: The method object. + """ + + site: SiteSchema = None + method: MethodSchema = None + + class config: + validate_assignment = True + arbitrary_types_allowed = True + + +class ContribSchema(Base): + project: str + identifier: str + formula: str = None + is_public: bool = None + data: SimulatorSchema + structures: list = None + tables: list = None + notebook: dict = None + attachments: List[Union[Attachment, Path]] = None diff --git a/src/mrsimulator/contribs/tests/test_contribs.py b/src/mrsimulator/contribs/tests/test_contribs.py new file mode 100644 index 000000000..d271c9596 --- /dev/null +++ b/src/mrsimulator/contribs/tests/test_contribs.py @@ -0,0 +1,143 @@ +# -*- coding: utf-8 -*- +import os + +import csdmpy as cp +import numpy as np +from mrsimulator import Simulator +from mrsimulator import Site +from mrsimulator import SpinSystem +from mrsimulator.contribs import load_obj +from mrsimulator.contribs import mpcontribs_export +from mrsimulator.contribs import parse_method +from mrsimulator.contribs import parse_sites +from mrsimulator.contribs import save_obj +from mrsimulator.methods import BlochDecaySpectrum + + +def csdm_obj(): + x = np.arange(1024) + y = np.random.rand(1024) + return cp.CSDM( + dimensions=[cp.as_dimension(x)], + dependent_variables=[cp.as_dependent_variable(y)], + ) + + +def test01(): + site = Site(isotope="1H", shielding_symmetric={"zeta": 10, "eta": 0.1}) + method = BlochDecaySpectrum( + channels=["1H"], magnetic_flux_density=9.4, rotor_frequency="15000" + ) + + output_site = parse_sites(site) + output_method = parse_method(method) + + assert output_site == { + "isotope": "1H", + "ChemicalShift": {"isotropic": "0.0 ppm", "zeta": "-10.0 ppm", "eta": 0.1}, + } + + omega_0 = abs(method.channels[0].gyromagnetic_ratio * 9.4) + assert output_method == { + "larmorFrequency": f"{omega_0} MHz", + "spinningFrequency": "15000.0 Hz", + "spectralWidth": "25000.0 Hz", + "rotorAngle": "54.7356 degree", + } + + +def test02(): + site = Site(isotope="27Al", quadrupolar={"Cq": 10e6, "eta": 0.4}) + method = BlochDecaySpectrum(channels=["27Al"], magnetic_flux_density=11.7) + + output_site = parse_sites(site) + output_method = parse_method(method) + + assert output_site == { + "isotope": "27Al", + "ChemicalShift": {"isotropic": "0.0 ppm"}, + "Quadrupolar": {"Cq": "10.0 MHz", "eta": 0.4}, + } + + omega_0 = abs(method.channels[0].gyromagnetic_ratio * 11.7) + assert output_method == { + "larmorFrequency": f"{omega_0} MHz", + "spinningFrequency": "0.0 Hz", + "spectralWidth": "25000.0 Hz", + "rotorAngle": "54.7356 degree", + } + + +def system_setup(): + site = Site(isotope="27Al", quadrupolar={"Cq": 10e6, "eta": 0.4}) + csdm_data = csdm_obj() + method = BlochDecaySpectrum( + channels=["27Al"], magnetic_flux_density=11.7, experiment=csdm_data + ) + + sim = Simulator() + sim.spin_systems = [SpinSystem(sites=[site])] + sim.methods = [method] + return sim + + +def test_contrib_card(): + sim = system_setup() + omega_0 = abs(sim.methods[0].channels[0].gyromagnetic_ratio * 11.7) + + def get_card(project, identifier): + return { + "data": { + "site": { + "isotope": "27Al", + "ChemicalShift": {"isotropic": "0.0 ppm"}, + "Quadrupolar": {"Cq": "10.0 MHz", "eta": 0.4}, + }, + "method": { + "larmorFrequency": f"{omega_0} MHz", + "spinningFrequency": "0.0 Hz", + "spectralWidth": "25000.0 Hz", + "rotorAngle": "54.7356 degree", + "blah": "blah", + }, + }, + "project": project, + "identifier": identifier, + } + + output = mpcontribs_export( + sim, + None, + project="test", + identifier="blah-blah", + exp_dict={"blah": "blah"}, + ) + attachments = output[0].pop("attachments") + assert output == [get_card("test", "blah-blah")] + + assert attachments[0].name == "blah-blah.mrsim.json.gz" + assert attachments[1].name == "blah-blah-exp0.csdf.json.gz" + + site = sim.spin_systems[0].sites[0] + sim.spin_systems = [SpinSystem(sites=[site, site, site])] + output = mpcontribs_export( + sim, + None, + project="test", + identifier="mp-5733", + exp_dict={"blah": "blah"}, + ) + attachments = [item.pop("attachments") for item in output] + card = get_card("test", "mp-5733") + assert output == [card, card, card] + + for att in attachments: + assert att[0].name == "mp-5733.mrsim.json.gz" + assert att[1].name == "mp-5733-exp0.csdf.json.gz" + + +def test_save_load(): + data = csdm_obj() + save_obj("test-data.csdf.json.gz", data.dict()) + assert data.dict() == load_obj("test-data.csdf.json.gz") + os.remove("test-data.csdf.json.gz") diff --git a/src/mrsimulator/simulator/tests/test_simulator.py b/src/mrsimulator/simulator/tests/test_simulator.py index 3a00c54af..ad21b6b0f 100644 --- a/src/mrsimulator/simulator/tests/test_simulator.py +++ b/src/mrsimulator/simulator/tests/test_simulator.py @@ -324,6 +324,7 @@ def test_sites_to_pandas_df(): zeta = [59.8, 52.1, 69.4, 12.4] eta_n = [0.62, 0.68, 0.6, 0.5] Cq = [None, None, None, 5.3e6] + Cq_test = [None, None, None, 5.3] eta_q = [None, None, None, 0.34] spin_systems = single_site_system_generator( @@ -349,7 +350,7 @@ def test_sites_to_pandas_df(): i if i is not None else None for i in eta_n ] assert list(pd_o["quadrupolar.Cq"]) == [ - f"{i} Hz" if i is not None else None for i in Cq + f"{i} MHz" if i is not None else None for i in Cq_test ] # assert list(pd_o["quadrupolar.eta"]) == [ # i if i is not None else None for i in eta_q diff --git a/src/mrsimulator/spin_system/__init__.py b/src/mrsimulator/spin_system/__init__.py index 776d3f030..ad9ebff17 100644 --- a/src/mrsimulator/spin_system/__init__.py +++ b/src/mrsimulator/spin_system/__init__.py @@ -153,7 +153,7 @@ class SpinSystem(Parseable): sites: Union[List[Site], np.ndarray] = [] couplings: Union[List[Coupling], np.ndarray] = [] abundance: float = Field(default=100.0, ge=0.0, le=100.0) - transition_pathways: List = None + transition_pathways: List[Union[TransitionPathway, List]] = None property_unit_types: ClassVar[Dict] = {"abundance": "dimensionless"} property_default_units: ClassVar[Dict] = {"abundance": "pct"} diff --git a/src/mrsimulator/spin_system/site.py b/src/mrsimulator/spin_system/site.py index eb6a30f23..b234fdc83 100644 --- a/src/mrsimulator/spin_system/site.py +++ b/src/mrsimulator/spin_system/site.py @@ -229,6 +229,13 @@ def parse_dict_with_units(cls, py_dict: dict): return super().parse_dict_with_units(py_dict) + def json(self, units=True, **kwargs): + py_dict = super().json(units=units, **kwargs) + if "quadrupolar" in py_dict and units: + value = float(py_dict["quadrupolar"]["Cq"][:-2]) / 1e6 + py_dict["quadrupolar"]["Cq"] = f"{value} MHz" + return py_dict + # Deprecated # def to_freq_dict(self, B0): # """ diff --git a/src/mrsimulator/spin_system/tests/test_site.py b/src/mrsimulator/spin_system/tests/test_site.py index 724743050..be1e6f786 100644 --- a/src/mrsimulator/spin_system/tests/test_site.py +++ b/src/mrsimulator/spin_system/tests/test_site.py @@ -165,7 +165,7 @@ def test_site_object_methods(): "alpha": "0.1 rad", "beta": "2.5 rad", }, - "quadrupolar": {"Cq": "10000000.0 Hz", "eta": 0.6}, + "quadrupolar": {"Cq": "10.0 MHz", "eta": 0.6}, } the_site = Site( isotope="27Al", diff --git a/src/mrsimulator/spin_system/tests/test_spin_systems.py b/src/mrsimulator/spin_system/tests/test_spin_systems.py index f025dcf26..f66067d39 100644 --- a/src/mrsimulator/spin_system/tests/test_spin_systems.py +++ b/src/mrsimulator/spin_system/tests/test_spin_systems.py @@ -178,7 +178,7 @@ def test_direct_init_spin_system(): { "isotope": "17O", "isotropic_chemical_shift": "-10.0 ppm", - "quadrupolar": {"Cq": "5100000.0 Hz", "eta": 0.5}, + "quadrupolar": {"Cq": "5.1 MHz", "eta": 0.5}, }, ], "couplings": [{"site_index": [0, 1], "isotropic_j": "34.0 Hz"}], diff --git a/src/mrsimulator/tests/save_load_test.py b/src/mrsimulator/tests/save_load_test.py index 820f5f5ab..b3877d72c 100644 --- a/src/mrsimulator/tests/save_load_test.py +++ b/src/mrsimulator/tests/save_load_test.py @@ -40,17 +40,20 @@ def test_save(): assert sim.methods[0].simulation is not None assert sim.methods[1].simulation is None - save("test.mrsim", simulator=sim, signal_processors=processors, params=None) + kwargs = dict(simulator=sim, signal_processors=processors, params=None) + save("test.mrsim", **kwargs), "test file" + save("test.mrsim.gz", **kwargs), "test gz" def test_load(): - sim_r, processors_r, report_r = load("test.mrsim") + for file_, tag in zip(["test.mrsim", "test.mrsim.gz"], ["file", "gz"]): + sim_r, processors_r, report_r = load(file_) - sim, processors = setup() - sim_r.methods[0].simulation = None - sim.methods[0].simulation = None - assert sim_r == sim - assert processors_r == processors - assert report_r is None + sim, processors = setup() + sim_r.methods[0].simulation = None + sim.methods[0].simulation = None + assert sim_r == sim, f"test {tag}" + assert processors_r == processors, f"test {tag}" + assert report_r is None, f"test {tag}" - os.remove("test.mrsim") + os.remove(file_) diff --git a/src/mrsimulator/transition/pathway.py b/src/mrsimulator/transition/pathway.py index 555fd4d64..2a4e9f10e 100644 --- a/src/mrsimulator/transition/pathway.py +++ b/src/mrsimulator/transition/pathway.py @@ -114,6 +114,17 @@ def __repr__(self): path = " ⟶ ".join([repr(item) for item in self._list]) return path + f", weight={self.weight}" + def dict(self) -> dict: + """Parse the class object to a python dictionary object. + + Example: + >>> pprint(path.json()) + {'pathway': [{'final': [0.5, -0.5], 'initial': [0.5, 0.5]}, + {'final': [-0.5, 0.5], 'initial': [0.5, 0.5]}], + 'weight': (1+0j)} + """ + return self.json() + def json(self, **kwargs) -> dict: """Parse the class object to a JSON compliant python dictionary object. diff --git a/src/mrsimulator/utils/importer.py b/src/mrsimulator/utils/importer.py index ba3cf744b..f313365db 100644 --- a/src/mrsimulator/utils/importer.py +++ b/src/mrsimulator/utils/importer.py @@ -4,6 +4,7 @@ from urllib.parse import urlparse from csdmpy.dependent_variable.download import download_file_from_url +from monty.io import zopen __author__ = "Deepansh J. Srivastava" @@ -14,6 +15,8 @@ def import_json(filename): res = urlparse(filename) if res[0] not in ["file", ""]: filename = download_file_from_url(filename) - with open(filename, "rb") as f: + + operator = zopen if filename.endswith("gz") else open + with operator(filename, "rb") as f: content = f.read() return json.loads(str(content, encoding="UTF-8")) diff --git a/tests/spectral_integration_tests/utils.py b/tests/spectral_integration_tests/utils.py index 2a89d695a..0b345f5b1 100644 --- a/tests/spectral_integration_tests/utils.py +++ b/tests/spectral_integration_tests/utils.py @@ -1,5 +1,5 @@ # -*- coding: utf-8 -*- -import json +# import json from os import path import numpy as np @@ -7,16 +7,11 @@ from mrsimulator import signal_processing as sp from mrsimulator import Simulator from mrsimulator import SpinSystem +from mrsimulator.utils.importer import import_json from numpy.fft import fft from numpy.fft import fftshift -def _import_json(filename): - with open(filename, "rb") as f: - content = f.read() - return json.loads(str(content, encoding="UTF-8")) - - def _get_header_and_footer(source_file): """ Return the number of rows in the header and footer. @@ -44,7 +39,7 @@ def get_data(filename): """Load a simpson or DMfit output file""" # source data - data_object = _import_json(filename) + data_object = import_json(filename) test_data_object = data_object["test_data"] source_file = test_data_object["filename"]